Nov 13, 2016 - ABSTRACT. Recent N-body simulations predict that large numbers of stellar black holes (BHs) could remain bound ... synthesis models for the field, where results are very sensitive to many uncertain parameters in the initial binary ....
Apr 17, 2008 - times well under a Hubble time, and escape speeds that are several .... Studies of main sequence binaries in globular clusters, which have.
son & Hernquist 1993) and has been demonstrated recently by direct N-body simulations .... was quite angry when much later in 1969 Lynden-Bell was given credit for the idea after a ... starâ in the Galactic Center Quintuplet cluster (Figer et al.
Feb 22, 2012 - Consequently, favourable configurations are initiated through the Kozai mechanism where ec- centricity growth is possible for long-lived triple ...
Jun 21, 2016 - are the sites of feedback processes, and likely drive the BH-host scaling relations. 3.1. Triggering of ... a few extreme outliers of very massive BHs in low-mass host galaxies have been reported recently .... the blazar OJ 287 with it
Feb 3, 2010 - 1Harvard-Smithsonian Center for Astrophysics, Cambridge, USA, 2University of Oxford, UK,. 3ESO, Garching ... FIGURE 1. MMT 6.5m optical spectrum of the NGC 404 nuclear cluster (black line) with the best fit ... material from the disk of
May 7, 2012 - 2 Interactive Research Center of Science, Tokyo Institute of. Technology .... in the cluster center. We discuss more details in section. 3.1. We additionally performed simulations of a single model with the same parameters as single-w6
Feb 10, 2011 - Auriga association have much higher binary frequencies (see. DuchÃªne 1999, and references therein). This leads to the ques- tion of what ...
Mar 12, 2012 - initially developed to improve the accuracy of black hole bi- nary mergers. The choice of how to modify .... In this paper, we are interested in the behavior of BHNS bi- naries for high mass ratio q = MBH/MNS ..... ment of Ontario; Ont
Dec 24, 2004 - star clusters throughout the local universe. These systems are large and young enough that they contain statistically significant numbers .... detailed evolution and ultimate fate of stars hundreds or thousands of times more massive th
Dec 24, 2004 - relation r = 3m1/2, then for the object to reach m Â» 103 in 3 Myr (to form .... Thrall, A. P., Deneva, J. S., Fleming, S. W., & Grabowski, P. E. 2003,.
May 2, 2012 - ABSTRACT. A promising mechanism to form intermediate-mass black holes (IMBHs) is the runaway merger in dense star clusters, where main-sequence stars collide and form a very massive star (VMS), which then collapses to a black hole. In t
Mar 12, 2012 - eration of ground-based detectors, Advanced LIGO  and. VIRGO . Indeed, whether for ... tion is limited by the minimum spacing of the numerical grid ..... We pre- fer to choose the free data (gij,K) according to the prescrip-.
Jun 29, 2010 - Conversely, the fastest merging model shows that rotation is quickly lost ... 2136/NGC 2137, in the Large Magellanic Cloud. .... server. This research has made use of the WEBDA database operated at the Institute of. Astronomy of the Un
Oct 30, 2012 - Email: [email protected] Devecchi .... For each simulation, we evolved 2000 binaries, one at a time, ...... NS retention fraction Pfahl et al.
The present model is most similar to recent models in McKee & Offner (2010),. OM11, Myers (2010), and ..... present-day protostar mass functions depends on accretion start times as discussed in. Section 2.3. This subsection ... where all stars start
Abstract. A model of core-clump accretion with equally likely stopping describes star formation in the dense parts of clusters, where models of isolated collapsing ...
Oct 12, 2011 - Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA. 2 Department of Astronomy, University of Maryland, College Park, MD. 20742, USA ...... Oka, K., & Manmoto, T. 2003, MNRAS, 340, 543. O'Neill, S. M., Miller, M. C., Bogda
Feb 5, 2008 - In this study we present the results from realistic N-body modelling of massive star clusters in the Magellanic Clouds. We have computed eight simulations with N ~ 105 particles; six of these were evolved for at least a Hubble time. The
Feb 5, 2008 - They measured a very large sample of Magellanic Cloud clusters .... lanic Cloud clusters possess radial surface brightness pro- ...... It is enlight-.
Feb 5, 2008 - GRAPE-6 special purpose computer (Makino et al. 2003) at ..... NGC 330), as well as in the Galaxy (e.g., Orion Nebula. Cluster, Arches ...
Nov 20, 2006 - short summary of recent developments in a key area of research, incorporating the main ..... 2002), and as confirmed by the CN cycled chemistry of the donor envelope, observed by Sabbi et al. (2003b) and ..... i The study of HÎ± lines
May 13, 2002 - ABSTRACT. We study the growth rate of stars via stellar collisions in dense star clusters, calibrating our analytic calculations with direct N-body simulations of up to 65536 stars, performed on the ..... corresponds to King  dim
D RAFT VERSION A PRIL 18, 2018 Preprint typeset using LATEX style emulateapj v. 6/22/04
BINARY MERGERS AND GROWTH OF BLACK HOLES IN DENSE STAR CLUSTERS RYAN M. O’L EARY , F REDERIC A. R ASIO , J OHN M. F REGEAU , NATALIA I VANOVA , AND R ICHARD O’S HAUGHNESSY
ABSTRACT We model the dynamical evolution of primordial black holes (BHs) in dense star clusters using a simplified treatment of stellar dynamics in which the BHs are assumed to remain concentrated in an inner core, completely decoupled from the background stars. Dynamical interactions involving BH binaries are computed exactly and are generated according to a Monte Carlo prescription. Recoil and ejections lead to complete evaporation of the BH core on a timescale ∼ 109 yr for typical globular cluster parameters. Orbital decay driven by gravitational radiation can make binaries merge and, in some cases, successive mergers can lead to significant BH growth. Our highly simplified treatment of the cluster dynamics allows us to study a large number of models and to compute statistical distributions of outcomes, such as the probability of massive BH growth and retention in a cluster. We find that, in most models, there is a significant probability (∼ 20 − 80%) of BH growth with final masses & 100 M⊙ . In one case, a BH formed with mass ≈ 620 M⊙ . However, if the typical merger recoil speed (due to asymmetric emission of gravitational radiation) significantly exceeds the cluster escape speed, no growth ever occurs. Independent of the recoil speed, we find that BH-BH mergers enhanced by dynamical interactions in cluster cores present an important source of gravitational waves for ground–based laser interferometers. Under optimistic conditions, the total rate of detections by Advanced LIGO, could be as high as a few tens of events per year from inspiraling BHs from clusters. Subject headings: galaxies: star clusters—globular clusters: kinematics and dynamics—black hole physics— gravitational waves 1. INTRODUCTION 1.1. Astrophysical Motivation
Many observations of globular clusters suggest the possible existence of intermediate-mass black holes (IMBHs) with masses ∼ 103 − 104 M⊙ in the centers of some cluster cores. The predicted masses of these IMBHs agree well with a simple extrapolation of the M − σ (mass – velocity dispersion) relation for galactic nuclei (Gebhardt et al. 2000). Observations and dynamical modeling of the globular clusters M15 in the Milky Way and G1 in M31 appear to be consistent with such a central massive BH (Gerssen et al. 2002, 2003; Gebhardt, Rich, & Ho 2002). However, N-body simulations by Baumgardt et al. (2003a,b) suggest that the observations of M15 and G1, and, in general, the properties of all corecollapsed clusters, could be explained equally well by the presence of many compact stars near the center without a massive BH (cf. van der Marel 2004; Gebhardt, Rich, & Ho 2005). On the other hand, these N-body simulations also suggest that many, perhaps most, non-core-collapsed clusters (representing about 80% of globular clusters in the Milky Way) could contain a central IMBH (Baumgardt et al. 2004b, 2005). Although the origin of IMBHs is not directly constrained by any observations, one possibility that has received considerable attention is the growth of a very massive object through successive collisions and mergers of ordinary massive stars in the cluster core (Colgate 1967; Ebisuzaki et al. 2001). Recent numerical studies by Portegies Zwart & McMillan (2002), Gürkan, Freitag, & Rasio (2004), and Freitag, Gürkan, & Rasio (2005) demonstrate that mass segregation and core collapse could proceed so quickly that there is rapid, runaway growth through stellar collisions before the most massive stars have evolved to supernovae (in . 3 Myr). Other numerical evidence comes from the
N-body simulations of Portegies Zwart et al. (2004), which demonstrate such growth in some young clusters. An important alternative, which we study in this paper, is the growth of an IMBH through successive mergers of stellar mass (∼ 10 M⊙ ) BHs (Quinlan & Shapiro 1989; Lee 1995, 2001). In this case, the massive stars formed initially in the cluster (with masses & 20 M⊙ ) must have avoided physical collisions and mergers and instead were able to complete their normal stellar evolution and produce BHs. Unlike massive stars, BHs have negligible cross sections for direct collision, so BH mergers can only occur through gravitational wave (GW) emission in binaries, possibly enhanced by dynamical interactions in the cluster (see Miller & Colbert 2004 for a review and §1.2). Even if it does not lead to growth and IMBH formation, the dynamical evolution of BHs in dense star clusters could also produce large numbers of tight BH–BH binaries that will merge through GW emission (possibly outside the cluster following dynamical ejection). These merging BH–BH binaries are likely candidates for direct GW detection and could even dominate the detection rates for ground–based laser interferometers such as LIGO (Portegies Zwart & McMillan 2000). 1.2. BH Formation and Segregation
In a globular cluster, it can be expected that a fraction ∼ 10−6 to 10−4 of the N ∼ 106 initial stars will become stellar-mass BHs via normal stellar evolution (Sigurdsson & Hernquist 1993). Assuming that all stars with initial mass greater than 20 M⊙ become BHs, the most recent Kroupa initial mass function (IMF) (Kroupa & Weidner 2003) gives a slightly higher initial BH fraction of NBH ≈ 1.5 × 10−3N, where we use mmin = .08 M⊙ and mmax = 150 M⊙ as the minimum and maximum stellar masses. When we scale this to the total mass of the cluster, Mcl , we find NBH ≈
O’Leary et al.
3 × 10−3(Mcl /M⊙ ). All of these BHs should have formed before ∼ 10 Myr, with the most massive BHs forming at around 3 Myr (Schaller et al. 1992). The radial distribution of these BHs in the cluster when they form is not known, but it is reasonable to assume that they should be much more centrally concentrated than the remaining main-sequence stars (MSs). This is for three reasons: (1) we expect significant mass segregation of the initial, higher-mass progenitors (Freitag et al. 2005); (2) there may be preferential formation of massive stars near the cluster center (see, e.g., Murray & Lin 1996; Bonnell et al. 2001); (3) BH birth kicks are not expected to be large enough to displace the BHs into the cluster halo (or eject them from the cluster; see White & van Paradijs 1996; Willems et al. 2005). Even if the BHs were initially distributed throughout the cluster, mass segregation would still be able to concentrate them into a central sub-cluster on a relatively short timescale. Indeed, a BH of mass MBH near the half-mass radius will be brought into the cluster core on the timescale tseg ∼
hmi trh , MBH
where trh is the relaxation time at the half-mass radius, and hmi is the average stellar mass (Fregeau et al. 2002). Considering a typical cluster with trh ∼ 1 Gyr, we see that a subcluster of BHs should still assemble near the center after at most ∼ 100 Myr. After a time very short compared to the overall dynamical evolution timescale of the globular cluster, the BHs should then form a self-gravitating subsystem within the core, which is dynamically decoupled from the rest of the cluster. This decoupling is sometimes referred to as Spitzer’s “mass stratification instability” (Spitzer 1969). The physics of this instability is by now very well understood, both for simple twocomponent systems and for clusters with a broad mass function (Watters et al. 2000; Gürkan et al. 2004). The dynamical evolution of the BH subsystem proceeds on a much shorter timescale since its relaxation time is now ∼ NBH /N times smaller than for the parent cluster. Interactions involving hardening of BH binaries will eventually lead to the ejection of BHs from the cluster, until there are so few left that the BH subsystem recouples dynamically (returns to “thermal equilibrium”) with the other cluster stars, and the interaction rate between BHs and cluster stars becomes comparable to the BH–BH interaction rate. For simple twocomponent clusters, Spitzer (1969) derived through analytic methods the condition necessary to reach energy equipartition: 1.5 m2 M2 < 0.16, (2) M1 m1 where M2 < M1 are the total masses of the two components, and m2 > m1 the mass of each individual object. Watters, Joshi, & Rasio (2000) used an empirical approach to find the more accurate condition 2.4 m2 M2 < 0.32. (3) M1 m1 For a cluster with total mass M1 = 106 M⊙ , BH mass m2 = 15 M⊙ , and average star mass m1 = 1 M⊙ , the cluster can be in equipartition, according to equation (3), when there are . 30 BHs in the cluster, i.e., significantly fewer than the number expected to form from the IMF. The minimum number of BHs
required to decouple from the cluster would likely be even less if their mass spectrum was considered (Gürkan et al. 2004). For most of its subsequent dynamical evolution, we expect the BH sub-cluster to remain largely free of other kinds of stars, even MSs with comparable masses. Indeed, for BHs to have formed, the most massive stars in the cluster must have avoided runaway collisions. This requires an initial halfmass relaxation time trh & 30 Myr (Gürkan et al. 2004). On the other hand, by the time the MS turnoff has decreased to 10 M⊙ , driving core collapse to concentrate the remaining 10 M⊙ MSs into the core (where they could then interact with the BHs) before they evolve would require trh . 20 Myr (Gürkan et al. 2004; see especially their Fig. 10). Therefore we see that the parameter space for clusters where both BHs and massive MSs would segregate and decouple simultaneously is probably very small, or nonexistent. This will allow us to concentrate on “pure” BH systems in our numerical simulations (§2). This picture is not altered significantly by the presence of binaries. During mass segregation, as the BHs start concentrating into the denser cluster core, they will likely interact with each other to form BH–BH binaries, if they were not already in binaries. BH–MS binaries will quickly undergo three-body or four-body exchange interactions that replace the lighter MS companion by a heavier BH (Sigurdsson & Phinney 1993). Therefore, as the BH subcluster begins to dynamically decouple from the rest of the cluster, a considerable number of hard BH–BH binaries are expected to remain. This has been demonstrated with direct N-body simulations by Portegies Zwart & McMillan (2000), where they found that almost all BH-binaries ejected were BH–BH. 1.3. Previous Studies Motivated by the absence of BH X-ray binaries, Kulkarni, Hut, & McMillan (1993) and Sigurdsson & Hernquist (1993) discussed the evolution and fate of 10 M⊙ BHs in globular clusters using simple analytic considerations. Their conclusions can be summarized as follows. Through dynamical interactions, hard BH–BH binaries tend to be hardened further (whereas soft binaries tend to be disrupted; Heggie 1975). Eventually, as the BH–BH binaries harden, the recoil produced by interactions becomes so great that the binaries can be ejected from the cluster. The timescale for merger by GW emission usually remains longer than the interaction time in the BH cluster so that hardened binaries are almost always ejected. Eventually the number of BHs is depleted, and no more than a few BHs would remain in the cluster core. This would imply that BH growth through successive mergers in the cluster cannot occur. Even if ∼ 1 BH remained at the center of every globular cluster today, it is unlikely that this would ever become detectable as an X-ray binary (Kalogera, King, & Rasio 2004). In order to better understand the complex interactions of BHs in clusters, Portegies Zwart & McMillan (2000) performed direct N-body simulations of systems with N = 2048 and N = 4096 total stars, including a small fraction (∼ 1%) of equal–mass “BHs” 10 times more massive than the other stars. They found that ∼ 90% of the BHs were ejected from the cluster after 4 − 10 relaxation times of the cluster (less than a few Gyr for most clusters), including ∼ 30% in BH– BH binaries. Similar N-body simulations were most recently performed by Merritt et al. (2004), who studied the formation of a core of BHs in a two–component cluster, with N = 104 . Merritt et al. (2004) found that the BHs completely segregate
Black Holes in Dense Star Clusters to the core within ∼ 100 Myr (ignoring initial mass segregation of the higher mass progenitors). These results are consistent with our qualitative understanding of the interactions of BHs, and strongly supports the assumptions we will make in § 2. Several studies proposed scenarios that could help BH– BH binaries merge inside clusters, opening the possibility of BH growth and IMBH formation through successive mergers. Miller & Hamilton (2002b) suggested that one larger seed BH may help overcome the Newtonian recoil, since the mass of the binary would be too large to have a recoil velocity greater than the escape velocity from the cluster, thus anchoring it in the core. It has also been proposed by Miller & Hamilton (2002a) and Wen (2003) that binary–binary interactions may have a large influence on BH mergers in the cluster core. These authors suggest that binary–binary interactions can produce significant numbers of long-lived hierarchical triple systems in which the outer BH increases the inner binary’s eccentricity via Kozai-type secular perturbations (Ford et al. 2000), thereby increasing the merger rate. Because these hierarchical triples may be driven to merge before their next interaction they should have a higher probability of staying in the cluster, and this can be a mechanism for retaining merging BHs. Gültekin, Miller, & Hamilton (2004) (hereafter GMH04) were the first to look at the possibility of successive mergers of BHs in a cluster core. In their simulations, GMH04 computed successive interactions between a BH–BH binary of varying mass ratio with single 10 M⊙ BHs sampled from an isotropic background. Between interactions the binary was evolved according to general relativity. They concluded that GW emission allowed for more mergers than previously thought possible, but the number of BHs required to form an IMBH would be much greater than believed to exist in a typical globular cluster. In our simulations, we treat not only binary–single interactions, but also four-body (binary–binary) interactions. Most importantly, we compute dynamical interactions for a realistic mixture of single and binary BHs self-consistently within a “pure BH” cluster core. We implement a treatment for the secular evolution of triples, including the Kozai mechanism (Wen 2003; Miller & Hamilton 2002a). We also include, for the first time in any dynamical treatment, a more realistic IMF for the BHs (to be detailed in §2.2). Our paper is organized as follows. In §2 we present our simulation methods and assumptions, as well as all our initial conditions. The main results of our simulations are shown in §3, where we look at evolution of all BH sub-clusters, and in §4, where we look at merging binary BHs as a source of GWs. Finally we conclude our paper in §5 with a summary and discussion of the implications of our results, and suggestions for further studies. 2. METHODS AND ASSUMPTIONS 2.1. Numerical Methods
Ideally, to simulate the evolution of ∼ 103 BHs in a massive cluster, one would like to do a full N-body simulation of the entire cluster, including the BHs and other stars. Or, capitalizing on the fact that the BH subsystem dynamically decouples from the rest of the cluster early on, one could perform an N-body simulation of just the BHs, subject to the cluster potential due to the other stars and the effects of dynamical friction, which tend to bring BHs that have been kicked out of the core—but not out of the cluster—back into the core on a short timescale. Although the second scenario is evidently
computationally much cheaper than the first to treat via direct N-body techniques, we decided to adopt an even faster technique which, although approximate, includes all the vital physics of the problem. This allows us to sample a wider range of initial conditions, and make a more thorough map of the parameter space of the problem. We treat the evolution of a BH subsystem in a background cluster subject to binary interactions, using a Monte Carlo technique to sample interaction rates and treat mass segregation, in conjunction with the small N-body toolkit Fewbody to numerically integrate binary interactions (Ivanova et al. 2005). We assume that the BH sub-cluster has a constant density and velocity dispersion throughout its evolution. The justification for this assumption is that with any reasonable initial binary fraction (& a few percent), the sub-cluster will enjoy a long-lived binary-burning phase in which its core parameters are roughly constant with time (Fregeau et al. 2003). The code we use is a modified version of the one presented in Ivanova et al. (2005), specially adapted to treat a BH system. Each dynamical binary interaction is followed using Fewbody, a numerical toolkit for simulating small-N gravitational dynamics that provides automatic calculation termination and classification of outcomes (for a detailed description see Fregeau et al. 2004). Fewbody numerically integrates the orbits of small-N systems, and automatically classifies and terminates calculations when an unambiguous outcome is reached. Thus it is well-suited for carrying out large numbers of binary interactions, which can be quite complex and long-lived, and thus must be treated as computationally efficiently as possible. Star clusters are characterized by a dense central core, surrounded by a much larger low-density halo. Consistent with this structure of a star cluster, we assume that all strong interactions occur within the core, which has a velocity dispersion σcore (cf. eq. ). If a product of an interaction has a velocity greater than the escape velocity from the core of the cluster, vesc , then it is assumed ejected from the cluster and is removed from the simulation. If it is less than vesc but greater than the escape velocity from the core into the halo vhalo , it is placed in the halo of the cluster, from where it can later reenter the core through dynamical friction with the background stars. Dynamical friction is implemented in our code by sampling from a Poisson distribution with an average timescale given in equation (1) with the average stellar mass, hmi = 1.0 M⊙ (see §3.3 from Ivanova et al. 2005). Between interactions, all BH–BH binaries are evolved according to the standard post-Newtonian equations (Peters 1964), da 64 G3 m21 m22 (m1 + m2 ) 73 2 37 4 (4) =− e + e 1 + dt 5 c5 a3 (1 − e2)7/2 24 96 de 304 G3 m1 m2 (m1 + m2 )e 121 2 (5) =− e , 1+ dt 15 c5 a4 (1 − e2)5/2 304 where m1 and m2 are the masses of the two BHs, a is the binary semimajor axis, and e is the orbital eccentricity. In some simulations, we account for linear momentum kicks imparted to the binary due to the asymmetry of the GW emission (Fitchett 1983). Because of the large theoretical uncertainty (Favata, Hughes, & Holz 2004; Blanchet, Qusailah, & Will 2005) in the recoil velocity of “major mergers” (i.e., when the mass ratio q = m2 /m1 & 0.4; here we assume m1 > m2 ), and the smaller but significant uncertainty from the spins of the BHs, we opted to neglect spin
O’Leary et al.
in determining recoil velocities in our simulations. We determine the overall recoil velocity of the merger remnant, Vrec , by using the form of the equation derived by Fitchett (1983), 4 f (q) 2GM/c2 ˜ Vrec = V0 , (6) fmax risco
where f (q) = q2 (1 − q)/(1 + q)5, fmax ≈ .0179, V0 is the maximum magnitude of recoil, and risco is the radius of innermost stable circular orbit. Fitchett (1983) found for circular orbits V0 ≈ 1480 km s−1 , much greater than the escape velocity from any globular cluster, whereas Favata et al. (2004) found V0 ∼ 10 − 100 km s−1 . In our simulations we set Vrec = V0
f (q) fmax
for ease of comparing V0 with other works. We use V0 near or slightly above the escape velocity of the cluster, up to 80 km s−1 . The form of equation (7) is consistent with the analysis of (Blanchet et al. 2005), who found the recoil velocity for the merger of two non-spinning BHs to be V0 ≈ 250 ± 50 km s−1 at the second post-Newtonian order. We do not account for GW recoil in the merger of the inner binaries that are part of hierarchical triples, because in the simulations where we include GW recoil, mergers in hierarchical triples are insignificant. Lee (1993) shows that for the velocity dispersions ( < 100 km s−1 ) and numbers of BHs ( . 103 ) expected in the star clusters we are investigating, the rate of two-body binary formation from gravitational bremsstrahlung is much less than that of regular (Newtonian) three-body binary formation, whereby a binary is formed with the help of a third BH, which takes away the excess energy needed to form the bound pair. Therefore, for our simulations, we only account for three-body and binary–binary (four-body) Newtonian interactions. In a dense sub-cluster of BHs, three-body binary formation can lead to the formation of a significant number of BH binaries, and eventually can help lead to the disruption of the entire sub-cluster. Ivanova et al. (2005) calculated the threebody binary formation rate for a binary of minimum hardness ηmin =
G m1 m2 , 2 bmax hmi σcore
where bmax is the maximum size of the region in which the three objects interact and σcore is the three-dimensional velocity dispersion of the core. The final rate per star they found for the formation of a binary with hardness η > ηmin is 5
Γ(η > ηmin ) = π
n2c G5 hmi f (m1 , m2 , m3 , η) 9 σcore
where n2 n3 m51 m52 −5 η (1 + 2η) n2c hmi5 hmi5 r v3 −1/2 v12 m1 m2 × 1+ η , 2 σcore (m1 + m2) hmi σcore
f (m1 , m2 , m3 , η) =
nc is the core density of the BHs, n2 is the core density of objects of mass m2 , n3 is the core density of objects of mass m3 , v12 is the relative velocity of the first object with respect to the second, and v3 is the relative velocity of the third object with respect to the center of mass of the first two objects. We follow the exact treatment of dynamical interactions of Ivanova et al. (2005) but include three-body binary formation
with ηmin = 1 in a consistent manner (see, in particular, their §3.4). We note that a more accurate criterion for the minimum binding energy we use should be based on the orbital velocity of the lightest member of the formed binary since we are not looking at equal–mass clusters (Hills 1990). However, since we typically do not have mass ratios above ∼ 10, we find our criterion to be sufficient. Because the code is not yet capable of following the evolution of triples between interactions, we must break them up before the next interaction time-step. In order to determine how to destroy the triple, we check if the binary is likely to merge before its next interaction. As a first step, we integrate numerically equations (4) and (5) (Peters 1964). We also scale the timescale of merger according to Miller & Hamilton (2002a) by calculating the maximum eccentricity from a first order Kozai mechanism approximation without post-Newtonian precession by solving equation (8) of Wen (2003). We then use the smaller merger time of the two methods. It is necessary to consider both methods because the scaling from Miller & Hamilton (2002a) overestimates the merging time in the instances when the Kozai mechanism is insignificant. In this case, the inner binary is merely perturbed and the eccentricity does not fluctuate, therefore the timescale of the merger should be the same as for an unperturbed binary. If the inner binary is likely to merge before the triple would interact with a field BH or BH–BH binary it is immediately merged, otherwise the triple is broken-up keeping the inner binary. Here, we assume that the outer BH is ejected from the triple with a velocity equal to its relative velocity with the center of mass of the triple. We keep the inner binary, but shrink the binary separation to conserve the energy of the system. 2.2. Initial Conditions and Parameters We use the results of Belczynski, Sadowski, & Rasio (2004) (hereafter BSR04), adopting the mass and binary distributions of their standard model at 11.0 Myr for our initial conditions (see, e.g., their Fig. 2, Fig. 4, & Fig. 5). BSR04 used a population synthesis approach to follow the evolution of a large number of massive stars and binaries, as would likely form in a massive star cluster. The model we base our calculations on has an initial binary fraction fbin = 50%, and follows the traditional Salpeter IMF for all initial stars with mass > 4 M⊙ . The binary fraction we use for our simulations is, of course, the final binary fraction found in BSR04, fbin ≈ 14%. Although most BH binaries have MS companions, we assume that these binaries will eventually become BH–BH binaries through exchange interactions. We create BH–BH binaries in their place, and select the companion BH, such that the distribution of the mass ratio, q, is uniform throughout the range qmin < q < 1 (qmin is set by the minimum BH mass in our distribution), consistent with observations for q & 0.2 (Woitas, Leinert, & Köhler 2001). We then increase the separation of the BH–BH binary assuming the binary preserves its binding energy in the exchange interaction. All wide binaries with orbital period P > 104 days are destroyed before our simulations begin. For each individual run, the mass of each BH is randomly selected with a distribution that reflects the results of BSR04, so that no two runs of any model contain the exact same BH population. The parameters used in all our simulations can be found in Table 1. For our simulations, aside from the exceptions noted in the table, we use self-consistent parameters determined by a King model, with W0 = 7, 9, and 11. Given a total cluster mass, Mcl , and core density, nc , we can calculate the one-
N OTE. — Starting from the top, the table is sorted according to the following parameters in their respective order: W0 , nc , σ1,BH , Mcl , NBH . a These models include GW recoil, as described in §3.4. Models v3e5k7ej54 and v3e5k7ej75 have maximum recoil velocities V = 54 and 75 km s−1 respectively. 0
Models v2e55k9e6, v2e55k9e65, v2e55k9e7, and v2e55k9e8 correspond to V0 = 60, 65, 70, and 80 km s−1 b These models include one primordial seed BH of mass 100 M and 200 M , respectively, as detailed in §3.3. ⊙ ⊙ c The GMH model parameters correspond to the initial conditions used in GMH04.
dimensional velocity dispersion, σ1,core , the escape velocity from the center of the potential to the half-mass radius, vhalo , the escape velocity of a BH in the core from the entire cluster, vesc , and the escape velocity of a BH at the half-mass radius from the entire cluster, vhalo,esc . We analyzed relatively massive and dense clusters, precisely the types of clusters where BH growth would be expected. Specifically, we systematically varied Mcl between 5 × 105 and 2 × 106 M⊙ , and nc between 105 and 107 pc−3 . Most of our simulations had 512 BHs (NBH = 512), but in two simulations we looked at clusters with smaller and larger numbers of BHs. All of these parameters are listed in Table 1. The escape velocities are used to determine whether the product of an interaction is to remain in
the cluster as prescribed in §2.1. The velocity dispersion of the BHs, σBH , can be related √to the one-dimensional velocity dispersion simply as σBH = 3σ1,BH . The properties of the decoupled sub-cluster in relation to the properties of the cluster core are not so apparent. In our simulations we need to know the initial velocity dispersion of the BHs, σBH . Numerical simulations suggest that the mean kinetic energy of the dynamically decoupled massive objects is only a few times larger than for other objects in the core (Gürkan et al. 2004). When decoupling, the BHs will have an
O’Leary et al.
F IG . 1.— Sub-cluster evolution. The data for each panel are binned and then averaged over several simulations in order to reduce the level of noise in the graph. The clusters reach equipartition when the total number of BHs in the core (solid) is . 40. Notice how few BHs are in the cluster halo (long-dashed) throughout the evolution. Only in v2e55k11, panel (c), do the numbers of BHs in the core and halo become comparable, but only when the cluster has already approximately reached equipartition. The more massive cluster v22e6e5k9, panel (a), does not reach equipartition in a Hubble time, suggesting massive clusters with low W0 could still contain significant numbers of single BHs in their cores. Each model has the same core density, nc = 5 × 105 pc−3 , and K = 4, except model v3e55k9, panel (d), which has K = 1.77. Model v22e6e5k9 is a massive W0 = 9 King model with Mcl = 2 × 106 M⊙ . Model v2e55k11 is a W0 = 11 King model with Mcl = 5 × 105 M⊙ . Models v2e55k9, panel (b), and v3e55k9 are the same as v2e55k11, except they have W0 = 9.
effective velocity dispersion which we write as 1/2 hmi σBH = K σcore , (11) hMBH i where σcore is the velocity dispersion of the core, hMBH i is the average BH mass, and K is the ratio of the mean kinetic energy of the BHs to that of the rest of the core, with K = 1 corresponding to complete energy equipartition. Because the BHs are dynamically decoupled from the rest of the cluster, they will undergo their own independent evolution as a small cluster of stars, increasing their density and velocity dispersion in the process. Therefore we look at sub-clusters with K = .64, 1.77, 4, and 16. 3. CLUSTER DYNAMICS AND BH GROWTH
In this section, we discuss our results for the evolution of the BH sub-cluster and we analyze the conditions under which successive mergers of BHs lead to the growth and retention of candidate IMBHs. In most of our simulations the BH sub-clusters reach equipartition (our standard for the complete disruption of the sub-cluster) in a few Gyr. In the most extreme cases the subcluster may dissipate in less than 100 Myr, as in the very dense
model v3e5e7k11; or not at all, as in all models with K = 16. We also find there is a significant probability a BH with mass & 100 M⊙ can form. For clusters that reach equipartition, massive BHs are most likely to form in clusters with high core densities and temperatures. This also leads to more growth for a given King profile. Also, King models with lower W0 are more likely to form candidate IMBHs. The growth of massive BHs is highly dependent on the maximum GW recoil velocity, nearly stopping it completely even when it is only a few percent larger than the core escape speed. 3.1. Fate of the BHs
The fate of the BHs in the cluster is determined mainly by their characteristic interaction rate in the core. However, two main questions remain regarding how the space of initial cluster parameters is divided among the different possible outcomes: does the BH sub-cluster reach equipartition with the cluster stars within a Hubble time? Do the BHs experience enough successive mergers to form an IMBH? Because of the simplicity of our model and speed of our code, the trends associated with varying one parameter in the simulation are evident. Figure 1 shows the number of BHs located in the cluster
Black Holes in Dense Star Clusters
F IG . 2.— Mass distribution of largest remaining BH. The largest remaining BH is simply the most massive BH remaining in the cluster at equipartition. Panel (a) is the mass distribution from a few different cluster types, with the model name in the upper right-hand corner. Panel (b) is the mass distribution from cluster model v2e55k9, with varying GW recoil kick velocities, V0 , corresponding to, from the top of the figure down, models v2e55k9, v2e55k9e6, v2e55k9e65, and v2e55k9e7.
core and halo as a function of time for a variety of different cluster types. We see that most BHs remain in the cluster core rather than in the halo of the cluster. A fraction . 3% of all the BH–BH mergers from cluster occur within the halo. The sudden drop in the binary fraction, as seen again in Figure 1, can be attributed to strong binary–binary interactions. Even though this causes a low binary fraction for the majority of the dynamical evolution of the sub-cluster, approximately 10− 15% of the BHs are ejected in binaries when the BH subcluster reaches equipartition in a Hubble time. For sub-clusters with low densities and high velocity dispersions, the three-body binary formation rate can become so small that a newly formed binary is disrupted before the next binary formation. Despite this, our simulations suggest that there always exists at least one binary in the core that may be able to keep the sub-cluster from undergoing core collapse. This binary, among others, is created through three-body binary formation, which is the underlying mechanism for the entire evolution of the sub-cluster. Since BHs are only ejected from the cluster after three- or four-body interactions, BH sub-clusters that do not produce binaries at a high rate do not dissolve within a Hubble time (see, e.g., models e5e5king7, e5king9, and v22e6e5k9). If the sub-clusters’ parameters have not changed significantly over the evolution of the entire cluster, then a significant number of single BHs could exist in clusters similar to those considered in our simulations. 3.2. Formation of an IMBH
Our simulations uniquely allow us to follow the interactions of a realistic mix of single– and binary–BHs and monitor the growth of BHs through successive mergers. We find that in a given King model, clusters with greater core densities and larger masses have a greater probability of forming an IMBH. They usually grow to even larger masses. For a fixed mass cluster that fits a given King model profile larger core densities result in a smaller half-mass relaxation time, trh . If this timescale is small enough, the cluster may actually undergo runaway growth through stellar collisions. This suggests that
one is more likely to find an IMBH, whether formed through successive BH mergers or stellar collisions, in clusters with dense cores. However, clusters with lower degrees of central concentration (smaller W0 , and higher σcore ) seem to have more BH growth. Figure 2 shows the mass distributions of the largest BHs formed in a few different cluster types. We see that when clusters completely disrupt in less than a Hubble time, larger core temperatures directly correlate with more growth of a single massive BH. This is evident in the directly comparing models v2e55k9 and v3e55k9 (where K = 4, and 1.77 respectively). Although the difference in the average number of mergers in the cluster is small—v2e55k9 had about 14 mergers whereas v3e55k9 had about 13—the number and sizes of the large BHs formed are dramatically different. Model v2e55k9 had about twice the number of BHs as v3e55k9 with mass > 100 M⊙ remain in the core at equipartition, and an average mass of the largest BH remaining in the core about 60% larger as well. In model v3e55k9, the velocity dispersion of the BHs, σBH , is lower than v2e55k9. As is evident in equation (9), the interaction rate increases rapidly with lower velocity dispersions, and hence leads to a quicker dissipation of the BH sub-cluster. This is consistent with our understanding that three-body binary formation drives the evolution to the eventual equipartition of the system. In clusters where there are successive BH mergers, the first formation of BHs with mass & 100 M⊙ occurs at ∼ 10 Myr after the sub-cluster formed, roughly independent of the cluster model. The probability that the cluster has a BH with mass & 100 M⊙ is then proportional to the logarithm of time (until the cluster reaches equipartition). On such a short timescale, one may worry that the BHs could also be colliding with massive stars, counter to our assumption of a "pure BH" system. However, as we discussed in § 1, for clusters in which runaway collisions of massive stars were avoided, the massive stars are not expected to enter the BH core before exploding as supernovae.
O’Leary et al.
3.3. Proposed Mechanisms for BH Retention One explanation for why BHs would stay in the cluster and not be ejected by hardening and eventual recoil is the Kozai mechanism in stable hierarchical triples (Miller & Hamilton 2002a; Wen 2003). If the relative inclination between the orbital planes in a given triple is large enough, the inner binary’s eccentricity can be driven to a high value (formally to 1). The high eccentricity necessarily means a decrease in the merging time of the inner binary so that the inner binary will merge before the triple’s next likely interaction. We find that, typically, this has no significant effect on the formation of large BHs. Overall, mergers caused by the Kozai mechanism account for less than ≈ 10% of the total mergers, and a negligible amount in most models. For example, the model with the largest percentage of triple mergers, v3e5k11, exhibits almost no growth at all, with none of the 14 runs having a BH with mass greater than even 80 M⊙ . Another possibility that favors the retention of large BHs in clusters is the introduction of a massive seed BH (Miller & Hamilton 2002b). The mass spectrum of BHs given in BSR04 generally includes a BH of ∼ 50 M⊙ in each cluster. The presence of a BH of this mass does not always mean that this BH will remain in the cluster. In fact, the largest initial BH is almost always ejected from the cluster. Although Miller & Hamilton (2002b) suggest that the introduction of a seed BH with mass 50 M⊙ may be sufficient to cause significant growth, we find that this is still not massive enough, in agreement with GMH04. In model v2e55k9–100 we introduce an initial seed BH of mass 100 M⊙ ; in model v2e55k9–200 a 200 M⊙ seed BH. We find that even BHs with these large masses can easily be ejected from the cluster. For example, in model v2e55k9– 100, the seed BH was retained only 18% of the time. Even in model v2e55k9–200, with a seed BH of 200 M⊙ , the seed BH is again retained only 35% of the time. These probabilities are still lower than those in GMH04, where they found the BHs to be retained in a similar cluster ∼ 40% and ∼ 90% of the time, respectively. This discrepancy can likely be attributed to the mass distribution of BHs in our simulations. The average mass of the BHs in our simulation is 50% higher than in GMH04, and thus there is an increased probability that the large seed BHs will be ejected. 3.4. GW Recoil One key question regarding coalescing binary BHs is the magnitude of linear recoil caused by the asymmetry of the GW emitted by the binary. Our code allows us to prescribe a systemic recoil velocity for every BH–BH merger. Through this we can follow the effects of GW recoil in a cluster environment. We are able to understand the possible effects of gravitational recoil by varying the maximum magnitude of the recoil velocity in the base model v2e55k9 (Favata et al. 2004; Blanchet et al. 2005). As expected, the number of successive mergers has a strong dependence on the maximum recoil velocity; however, it seems that it has only a small effect on the overall dynamics of the rest of the BH sub-cluster. Increasing the recoil has almost no noticeable effect on the total number of mergers in the cluster, but can have significant consequences on the rate of BH–BH inspirals (see §4). Even a maximum recoil velocity slightly larger than the escape velocity of the core (V0 = 1.042 × vesc = 60 km s−1 ) dramatically changes the number of large BHs formed in this model, as seen in Figure 2.
F IG . 3.— Eccentricity distribution of merging BH binaries. For this log − log plot, we show the eccentricity distribution of all BH binary mergers for model v2e55k9. The distribution of eccentricities is almost entirely independent of the model used. The two frequencies of the GWs were chosen to show the expected eccentricity distribution of a binary as it enters the observable bands of both ground–based, ≈ 10 Hz, and space–based interferometers, ∼ 10−3 Hz. The low eccentricity of most binaries entering the ground–based interferometers detection band suggests almost no loss of detectable BH–BH binary signals if only circular templates are used for analysis.
With this recoil velocity, the probability of forming a 100 M⊙ is cut in half. When we look at higher recoil velocities, the possibility of growing a large BH gets only smaller. This is not such a surprising result when one considers our simple prescription for modeling GW recoil. In mergers with a mass ratio, q = m2 /m1 ≈ .38, or when f (q) ≈ fmax , the coalescing binary will generally be ejected when V0 is close to the core escape speed (see eq. ). For BHs which go through successive mergers, it is likely that before reaching masses > 100 M⊙ the BH will have already been ejected. However, since this only affects BHs after they merge, the number of mergers remains relatively unaffected. Our simple model of GW recoil neglects some aspects of the process, which could alter results. In particular, we do not follow the evolution of the BHs’ spins. Because of this, we must neglect the effects of spin on recoil, and therefore do not look at alternative paths to large BH retention. One possibility is that if the BHs with the appropriate spins merge, some clusters could still retain larger BHs, even if the recoil velocity of the BHs with no spin is much larger than the escape speed. On the other hand, spin breaks the symmetry of the binary and can lead to large recoil velocities even for equal– mass binaries (Favata et al. 2004). 3.5. Properties of Merging BHs Inspiraling BHs with high eccentricity have very complex gravitational waveforms, and their additional free parameters make computational searches even more expensive. Therefore we want to know the eccentricity of the merging BHs’ orbits as they enter possible detection bands of current and future gravitational wave detectors. Many ground– based laser interferometers are currently in operation. These detectors, such as LIGO, are sensitive to GWs above ≈ 10 Hz, below which seismic noise dominates the noise curve
Black Holes in Dense Star Clusters (Cutler & Flanagan 1994). Also, in the planning stages is a space–based laser interferometer, LISA, which, with its longer baseline, is expected to be sensitive to GWs with frequencies ∼ 10−3 − 1 Hz (Bender et al. 1998; Bender & Hils 1997). Dynamical interactions, such as those in a cluster core, coupled with the strong dependence of merging time with eccentricity, suggest that many binaries will have highly eccentric orbits after their last strong interaction. Of course, GW emission reduces the eccentricity of the orbit by circularizing the binaries (see eq. ), and therefore it is often assumed that most binaries can be fitted with GW templates for zero eccentricity. Figure 3 shows that there will be almost no loss in possible LIGO BH–BH binary sources with this assumption because, at such a high frequency, almost all binaries are circular. However, the eccentricity distribution will likely matter for space–based detectors since the inspiraling binaries have not been entirely circularized by the GW emission. This is consistent with the previous study by GMH04. Another important factor in detecting inspirals is the chirp mass of a binary, Mchirp =
(m1 m2 )3/5 , (m1 + m2 )1/5
which solely determines the overall magnitude of the GWs emitted by a coalescing circular binary. Because our simulations allow for successive mergers of BHs, as shown in Figure 4, some mergers have chirp masses well above those expected if dynamics were not included. Because of our realistic initial distribution of BH masses, the chirp masses of most mergers is above the expected value for two 10 M⊙ BHs merging, 8.7 M⊙ . Although, as we discussed in § 3.2, successive mergers of BH–BH binaries can lead to more inspirals of massive BHs, in most cases the chirp mass distribution at the end of the evolution of the cluster is not significantly different from early in the evolution. For one, about half of the mergers occur outside of the cluster, which were ejected through dynamic interactions early in the cluster’s history before significant growth had occurred. Also, only in a few cluster models is the probability of BH growth near unity, in which the massive BHs dominate the mergers in the cluster. In Figure 5, we plot the average merger rate of model v2e55k9 as a function of time. The merger rate can clearly be broken into two regimes. The first occurs early in the cluster evolution, when the number of binaries has not been completely depleted. The second regime occurs after the cluster is nearly all single BHs with only a few binaries. During this later time the merger rate of the cluster falls off inversely proportional to the age of the cluster, independent of cluster model. 3.6. Comparison with Previous Studies
Overall our results agree with those of previous direct N-body simulations. In the simulations of Portegies Zwart & McMillan (2000), where N = 2048, 4096 and NBH = 20, 40, approximately 60% of the BHs where ejected as single BHs and 30% ejected as binary BHs. The lower values of binary ejection, which we found, can be explained by in-core mergers and the lower order interactions followed compared to N-body simulations (the calculations of the authors were Newtonian only, and so did not allow for mergers via GW emission). As shown in Table 2, in
most runs 70 − 85% of the BHs were ejected as singles and 10 − 15% were ejected in binaries. The biggest divergence between our simulations and those conducted by Portegies Zwart & McMillan (2000) is the energy distribution of ejected BH binaries. Portegies Zwart & McMillan (2000) found that the ejected BH binaries had a binding energy distribution more or less uniform in the logarithm. In Figure 6 we show a typical binding energy distribution of the ejected BHs. Our simulations produce a distribution much more log–normal. Every model we analyze has a similar distribution of ejected BH binary binding energy, with values between ≈ 103 − 105 kT , consistent with analytic considerations (Kulkarni et al. 1993). This divergence can possibly be explained by the low number of BHs used in Portegies Zwart & McMillan (2000). To see what effect the mass spectrum has on our simulations, and also to compare our data to the results of GMH04, we use three different models. Each simulation uses the same velocity dispersion and escape velocities as in GMH04 (see model GMH in Table 1), but starts with a different mass function. GMHA and GMHC both have 10 M⊙ and 15 M⊙ equal– mass BHs respectively. GMHB has BHs with a mass distribution consistent with our other simulations following BSR04. This distribution has a corresponding average mass of about 15 M⊙ (see §2.1 & §2.2). As can be seen in Table 2, using equal–mass BHs in our simulations increases the number of binaries ejected by almost a factor of two. The models with equal–mass BHs, GMHA and GMHC, have about 20% of their BHs ejected as binaries, whereas GMHB has only about 10%. Figure 7 shows how the mass spectrum we use causes the cluster to reach equipartition at an earlier time and changes the timescales of when mergers occur, compared with a cluster of equal–mass BHs. 4. BINARY BHS AS SOURCES OF GRAVITATIONAL WAVES
BH–BH binaries formed through dynamical interactions may be some of the best sources of GWs detectable by ground–based laser interferometers. Previous studies of detection rates have led to a large range of possible values, with some of the greatest uncertainty coming from the dynamics of interacting BH binaries in massive clusters (Tutukov & Yungelson 1993; Portegies Zwart & McMillan 2000). In this section we determine possible maximum detection rates of BH–BH binary inspirals from some globular cluster models. 4.1. Calculation of the LIGO Detection Rate
To calculate accurately the net detection rate, one should convolve self–consistent densities and birth rates of observed star clusters throughout the universe with the mergers rates in our simulations. Because there is great uncertainty in these distributions, we look at each of our cluster models, and determine how cluster densities and masses may affect the distribution of BH–BH inspirals, leaving a detailed analysis for further study. Although current ground–based interferometers are not sensitive enough to detect mergers of inspiraling BH–BH mergers to any cosmologically significant distance, as these detectors become more sensitive they will be able to make detections to ever farther distances. It then becomes important that a consistent cosmological model is used to have a better understanding of detectable inspirals. In all of our calculations we use the best–fitting cosmological parameters found by Spergel et al. (2003) from combining WMAP data
O’Leary et al.
F IG . 4.— Chirp mass of mergers versus time. This is a comparison of the two models v2e5k11 and v22e6e5k9, in panels (a) and (b) respectively. Plotted is the chirp mass versus time of all mergers of 46 random runs of model v2e5k11 and all 46 runs of v22e6e5k9. Model v2e5k11 is one of the least efficient clusters in producing large BHs and BH–BH binary mergers in general. Therefore, the distribution is most nearly that expected from the initial mass distribution of BSR04. Because of how quickly v2e5k11 evolves (teq ≈ 200 Myr) almost all mergers in later times occur outside the cluster. In comparison, v22e6e5k9 is a massive cluster that does not reach equipartition before a Hubble time. There is still a significant fraction of BHs in the cluster at the end of the simulation, which allows for more growth, and also more massive BH mergers.
F IG . 5.— Merger rate vs. time. The solid curve is the average merger rate of model v2e55k9 as a function of time. The dotted line is a power– law ∝ time−1 . After ∼ 108 yr, the merger rate is inversely proportional to the age of the cluster. The evolution of the merger rates can be split into two phases. The first when the cluster is undergoing many binary interactions, and the second, when the binary fraction is depleted and nearly zero. These two phases of merger rates appear consistently in all cluster models.
F IG . 6.— Energy distribution of ejected BH binaries. Plotted is the probability distribution of the energy of all BH–BH binaries ejected before equipartition in 117 runs of model v2e55k9. The energy is plotted in units of the mean kinetic energy kT , where 3/2 kT is the mean stellar kinetic energy of the MS stars in the core of a cluster of this type. We find that all other models have a distribution very similar to the one shown above.
with other measurements: H0 = 71 km s−1 Mpc−1 , Ωm = 0.27, Ωγ = 5 × 10−5, and ΩΛ = .73. In our calculations, we assume that the globular cluster model was formed uniformly through the universe at a given cosmological time corresponding to redshift zform . We then record each detectable merger into one of 100 bins each with time width δt = t0 /100, where t0 is the current age of the universe, based on when the merger occurred. If di is the number
of detections in bin i, we sum over the rate of each bin giving the final rate: Rzform =
100 X di 4π ρ0 (D3i − D3i−1 )(1 + zi)−1 , δt 3
where ρ0 is the current density of a given cluster model and zi is the redshift to bin i. With te = t0 − iδt, the proper distance to
Black Holes in Dense Star Clusters
F IG . 7.— A comparison of the models GMHA, GMHB, and GMHC. In panel (a), the number of BHs as a function of time is plotted. Each pair of lines is labeled above by its associated model that we use to compare our results with GMH04. Model GMHA, with only 10 M⊙ equal–mass BHs, on average reaches equipartition at teq ≈ 5.2 Gyr. Model GMHC, a cluster with only 15 M⊙ equal–mass BHs, reaches equipartition at teq ≈ 830 Myr. Finally, GMHB, with a varied mass spectrum as in BSR04 and a corresponding average mass ≈ 15 M⊙ , has not only a significantly smaller number of ejected binaries, but also reaches equipartition at the earliest time, teq ≈ 670 Myr. In panel (b), the number of mergers at a given time, denoted by redshift, is plotted. We assume that the clusters would currently be 13 Gyr old. The model with the varied mass, GMHB, not only dissociates more quickly, but has most of its mergers at an earlier time in the cluster’s evolution, i.e., at a higher redshift.
the edge of bin i is Di =
dt , a(t)
where a(t) is the scale factor of the Friedman-RobertsonWalker metric that satisfies the Einstein equation, and a(t0 ) ≡ 1. The factor (1 + zi )−1 in equation (13) comes from the cosmological time dilation of the merger rate. The final detection rate is directly proportional to the density of such clusters in the universe. For our calculations we assume that ρ0 = 1 Mpc−3 , independent of cluster model, for ease of comparing our results with other works. To put this in perspective, Portegies Zwart & McMillan (2000) found the number density of all globular clusters to be ρ0 ≈ 8.4 h3 Mpc−3 based on rough fits to observations. Our analysis doesn’t include the full distribution of cluster parameters, but instead looks at each cluster individually. The Milky Way, for example, contains hundreds of globular clusters, with only a fraction of clusters similar to those we look at in this study (Harris 1996). In order to be as precise as possible, we must determine which mergers could actually be detected by a given version of LIGO. To do this, we must look at the accumulated signal– to–noise ratio (SNR) of a given merger at redshift zm . Since, gravitational waveforms are invariant under redshift, we use the redshifted chirp mass of a given merger Mchirp = (1 + zm )Mchirp ,
DL = (1 + zm)D,
and its luminosity distance
to calculate the SNR of the merger, where D is the proper distance of the merger as calculated in equation (14). We look at all BH–BH binary mergers caused by interactions before the cluster reaches equipartition, and determine
if it would be detected by a given ground–based GW interferometer by comparing its SNR to that of an inspiraling neutron star–neutron star (NS–NS) binary at luminosity distance DL,0 . In practice, we determine the merger to be detectable if 5/6 s s( foff ) (S/N)( foff ) DL,0 Mchirp = > 1, (17) (S/N)(DL,0) DL Mchirp,0 s( foff,0 ) with s( foff ) =
( f ′ )−7/3 ′ df , Sn ( f ′ )
where Sn ( f ′ ) is the detector’s noise spectrum, Mchirp,0 is the effective chirp mass of the inspiraling NS–NS binary, and foff,0 is the cut–off frequency of the NS–NS merger (Cutler & Flanagan 1994). In this study, we assume that the mergers are isotropic and neglect the orientation of the merger relative to the detector. We use an analytic fit for the shape of Advanced LIGO’s noise spectrum found in Cutler & Flanagan (1994), ∞ f ′ < 10 Hz, ′ (19) Sn ( f ) ∝ ( f0 / f ′ )4 + 2[1 + ( f ′/ f0 )2 ] f ′ ≥ 10 Hz, where f0 = 70 Hz. We do not include the coefficient as calculated by Cutler & Flanagan (1994) in this equation since it is completely canceled out in equation (17) and rescaled by DL,0 . The final piece of equation (17) is the cut–off frequency of the merger, foff , after which the waveform of the GWs can not be accurately modeled. This frequency is generally regarded to be the frequency of GWs at the binary’s last circular orbit (LCO), after which the BHs plunge into each other in a time less than the orbital period. In their calculations, Kidder, Will, & Wiseman (1993) found that for two 10 M⊙ Schwarzschild BHs the frequency of the orbit at the LCO is ≈ 100 Hz, using what they called “hybrid” equations that
O’Leary et al.
were valid through (post)5/2 –Newtonian order for arbitrary masses. Because their calculations also showed that this is a lower limit of the orbital frequency of the binary for an arbitrary mass ratio, we choose to use this limit for calculating foff . For circular orbits, the GW frequency is twice the orbital frequency, therefore in our calculations we use 20M⊙ 1 foff ≈ 200 Hz (20) M 1 + zm
as the cut–off frequency of detectable GWs, where M = m1 + m2 and zm is the redshift of the merger. The location of the LCO and its corresponding orbital frequency are far from well established. Nevertheless, numerical simulations and other analytic approximations have shown the orbital frequencies of equal–mass BHs to be only larger than the value we use here (Blanchet 2002; Grandclément, Gourgoulhon, & Bonazzola 2002). 4.2. Results
We see, in Figure 8, the expected detection rates of BH binary mergers for cluster v2e55k9. Here we assume that all clusters of this model formed when zform = 7.84, or 13 Gyr ago, using the assumptions and equations of § 4.1. As can be seen in the figure, for the luminosity distances Advanced LIGO is expected to reach, the detection rate scales well by the power–law D−3 L . This can be explained by the time evolution of the merger rate of BHs. We find, for all simulations, that the merger rate scales as t −1 after the disruption of the primordial binaries. For the distances of mergers Advanced LIGO is expected to be able to detect, the clusters are relatively old and hence, the rate changes very little. The detection rates and theoretical uncertainty of all models can be found in the last two columns of Table 2. We find for model v2e55k9, for a version of LIGO capable of detecting NS-NS mergers at a luminosity distance, DL,0 ≈ 190 Mpc, the net detection rate is ≈ 2.7 yr−1 . For v2e55k9e8, which has the same conditions as v2e55k9, but a recoil velocity V0 = 80 kms−1 , the net detection rate is actually higher than v2e55k9 at ≈ 4.1 yr−1 . One cause for this may seem to be the lower cut–off frequency of high mass binary mergers, but this is wrong. When analyzing the detection rate assuming a universal chirp mass for all mergers, the rate was still higher. It can actually be attributed to the higher merger rate later in the evolution of the cluster. This delayed merger rate may be a result of a few mechanisms. The dependence of merger rate in models v2e55k9-100 and v2e55k9-200 suggests that having more massive objects in the cluster results in a lower rate overall at the end of the evolution. Model v2e55k9-200, which has a 200 M⊙ seed BH, has a merger about half that of model v2e55k9. Another, slightly less significant mechanisms can still possibly be the masses of the merging binaries. The timescale for merger is longer when the masses of the objects are less (see eqs.  and ). Within the level of theoretical uncertainty, the final detection rate is roughly proportional to the initial number of BHs. This, of course, is not exact, and isn’t expected to be. The timescale for cluster disruption is also proportional to the number of BHs, so it should be expected that there would be some variance in the merger rate at later times depending on the initial number of BHs. In order to get an estimated value of the actual Advanced LIGO detection rate, we must consider not only the number density of the massive globular clusters we look at here, but
F IG . 8.— Detection rate of BH–BH inspirals. The solid curve represents the expected detection rate of mergers from cluster v2e55k9 if it has a current density ρ0 = 1 Mpc−3 and formed at zi = 7.84, or 13 Gyr ago. The method and assumptions of our calculation are detailed in § 4. The dotted line is a power–law for R ∝ D−3 L,0 .
also the density of the low mass clusters. Our results suggest that the expected detection rate scales proportionally to the initial number of BHs in the cluster. Considering the number of BHs that Portegies Zwart & McMillan (2000) quoted as being in the globular clusters they analyzed, we would expect the actual rate to lie within the range of rates for the clusters reported in Table 2, rescaled to the number density of globular clusters in the universe (∼ 3 Mpc−3 ). 5. SUMMARY AND DISCUSSION
Our work is based on the reasonable assumption that, in a sufficiently large and dense star cluster, BHs created via stellar evolution concentrate in an inner core and effectively decouple from the rest of the cluster following mass segregation and the development of the Spitzer instability. Taking advantage of this decoupling, we have computed the dynamical evolution of the BHs in a highly simplified treatment of the stellar dynamics but covering a wide range of cluster models and repeating calculations for each model in order to obtain a complete statistical description of outcomes. Our assumed initial distributions of BH masses and binary parameters are based on the most recent population synthesis calculations for young stellar systems. In our simulations, we use a simple Monte Carlo method to follow the evolution of the BH subsystem in a fixed background cluster described by a King model. The BHs are allowed to interact only in the core, described by a constant density and velocity dispersion. All interactions involving binaries are computed exactly by direct (Fewbody) integrations but we implement three-body binary formation using a simple analytic rate formula (eq. ). Binary formation through dissipative two-body encounters is negligible in the systems we consider here (with velocity dispersions σ < 100 kms−1 ). We allow for Newtonian recoil of BHs into the halo, reintroducing them in the core following mass segregation. Between interactions, BH–BH binaries are evolved taking into account gravitational radiation and the possibility of a merger (eqs.  & ). In some simulations, we attempt to model kicks due to gravitational radiation recoil in merging
Black Holes in Dense Star Clusters binaries, parameterizing the large theoretical uncertainty in a simple analytic formula that depends on the binary mass ratio only (eq. ). We present the results of our simulations in § 3 and 4, and, in particular, we derive the probability of massive BH growth and retention in clusters (∼ 20 − 80%). We also show that the Kozai mechanism in triples has almost no significant effect on the merger rate, BH growth, or retention, in contrast to previous suggestions in the literature. In addition, we derive a net rate of BH–BH binary mergers detectable by current and future ground–based GW interferometers. We find, under extremely optimistic assumptions, that this can be up to a few tens of events per year, if we assume that the globular clusters likely formed over a time span zform . 8. When including recoil from the gravitational “rocket” effect (Fitchett 1983) this rate is nearly doubled due to the increased merger rate from a population of less massive BHs. The simulations presented here improve upon previous studies in that we not only account for binary evolution due to GW emission between successive interactions of a fixed group of BHs, but we also use a realistic BH IMF based on the most recent population synthesis models (BSR04). Most analytic studies and numerical simulations so far assumed that all BHs were 10 M⊙ (see, e.g., Miller & Hamilton 2002b and Portegies Zwart & McMillan 2000), or included just one large BH (GMH04). The mass spectrum from BSR04 gives a slightly higher average mass (≈ 15 M⊙ ), and also includes at least one significantly large BH, which is almost always ejected early in our simulations. Even some of the largest BHs formed from mergers, with mass ≈ 120 M⊙ , are ejected from the cluster when they are formed early enough in the simulation to interact often with other large BHs. The distribution of the binary parameters also has a significant effect on the interactions. The low primordial binary fraction results in very few formed triples, and even fewer mergers in triples, whether enhanced by the Kozai mechanism or not. The results of our simulations indicate a greater likelihood of moderate growth in globular clusters than previous N-body simulations have suggested possible, but they show varied results compared with GMH04. Our wider BH mass function makes the ejection of BHs with masses as high as ∼ 100 M⊙ very likely. On the other hand it also encourages the formation of even larger BHs through successive mergers, and the chirp masses of merging binaries are larger than would be expected for a cluster with equal–mass ∼ 10 M⊙ BHs. Despite the demonstrated growth, the probability for an IMBH of ∼ 103 M⊙ to form directly via successive BH–
BH mergers remains extremely small. However, we should not neglect the possibility of significant further growth of the final remaining BH through stellar collisions after complete evaporation of the BH sub-cluster. Results from Baumgardt, Makino, & Ebisuzaki (2004a) suggest that even a moderately massive ∼ 200 M⊙ BH could grow into a 103 M⊙ IMBH after only a few Gyr, well within the current ages of globular clusters. Initial conditions play a crucial role in determining the probability of IMBH growth. In general, massive, dense clusters with high core temperatures have the greatest likelihood of BH growth. Eccentricity growth through the Kozai mechanism, although it makes it possible for BHs to merge with very little (Newtonian) recoil, does not occur often enough to affect IMBH formation. Note, however, that our basic assumptions break down at late times, when the growth of BHs is most significant. Eventually the number of remaining BHs becomes small enough that the constant core conditions assumed in our simulations are no longer justified. In particular, the BH subsystem will start interacting with lower-mass stars at a significant rate and it will then likely return to energy equipartition. This is why we emphasize the values computed by our code at equipartition, as opposed to the values at the end of the evolution (§ 1.2). Beyond the return to equipartition, the evolution of the remaining BHs will be strongly coupled to the overall evolution of the whole star cluster. Our simulations are far from complete, yet the results show great promise for future work. For one, we can expect that the most massive clusters will have the greatest likelihood of BH growth. To enter the next regime of even more massive clusters, such as galactic nuclei, one must account for binary formation from gravitational bremsstrahlung (Lee 1993). Also, because of the simplifications we make, it is also possible to explore an even larger parameter space, especially of the more numerous and less massive clusters. By convolving these data with the cluster formation history of the universe, one could determine a BH–BH merger detection rate to an ever more accurate degree.
We would like to thank Marc Freitag for providing us with the code to determine the King model parameters used in our simulations. This work was supported by NASA Grant NAG5-12044 and NSF Grant NSF PHY-0245028. F.A.R. thanks the Center for Gravitational Wave Physics at Penn State for hospitality and support.
REFERENCES Baumgardt, H., Hut, P., Makino, J., McMillan, S., & Portegies Zwart, S. 2003a, ApJ, 582, L21 Baumgardt, H., Makino, J., & Ebisuzaki, T. 2004a, ApJ, 613, 1133 —. 2004b, ApJ, 613, 1143 Baumgardt, H., Makino, J., & Hut, P. 2005, ApJ, 620, 238 Baumgardt, H., Makino, J., Hut, P., McMillan, S., & Portegies Zwart, S. 2003b, ApJ, 589, L25 Belczynski, K., Sadowski, A., & Rasio, F. A. 2004, ApJ, 611, 1068 Bender, P., Brillet, A., Ciufolini, I., Cruise, A. M., Cutler, C., Danzmann, K., Fidecaro, F., Folkner, W. M., Hough, J., McNamara, P., Peterseim, M., Robertson, D., Rodrigues, M., Rüdiger, A., Sandford, M., Schäfer, G., Schilling, R., Schutz, B., Speake, C., Stebbins, R. T., Sumner, T., Touboul, P., Vinet, J.-Y., Vitale, S., Ward, H., & Winkler, W. 1998, LISA Pre-Phase A Report, 2nd ed. Bender, P. L., & Hils, D. 1997, Classical and Quantum Gravity, 14, 1439 Blanchet, L. 2002, Phys. Rev. D, 65, 124009 Blanchet, L., Qusailah, M. S. S., & Will, C. M. 2005, ApJ, accepted (astro-ph/0507692)
Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 2001, MNRAS, 323, 785 Colgate, S. A. 1967, ApJ, 150, 163 Cutler, C., & Flanagan, É. E. 1994, Phys. Rev. D, 49, 2658 Ebisuzaki, T., Makino, J., Tsuru, T. G., Funato, Y., Portegies Zwart, S., Hut, P., McMillan, S., Matsushita, S., Matsumoto, H., & Kawabe, R. 2001, ApJ, 562, L19 Favata, M., Hughes, S. A., & Holz, D. E. 2004, ApJ, 607, L5 Fitchett, M. J. 1983, MNRAS, 203, 1049 Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385; erratum 605, 966 Fregeau, J. M., Cheung, P., Portegies Zwart, S. F., & Rasio, F. A. 2004, MNRAS, 352, 1 Fregeau, J. M., Gürkan, M. A., Joshi, K. J., & Rasio, F. A. 2003, ApJ, 593, 772 Fregeau, J. M., Joshi, K. J., Portegies Zwart, S. F., & Rasio, F. A. 2002, ApJ, 570, 171
O’Leary et al.
Freitag, M., Gürkan, M. A., & Rasio, F. A. 2005, MNRAS, submitted (astro-ph/0503130) Gültekin, K., Miller, M. C., & Hamilton, D. P. 2004, ApJ, 616, 221 Gürkan, M. A., Freitag, M., & Rasio, F. A. 2004, ApJ, 604, 632 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., Richstone, D., & Tremaine, S. 2000, ApJ, 539, L13 Gebhardt, K., Rich, R. M., & Ho, L. C. 2002, ApJ, 578, L41 —. 2005, ApJ, submitted Gerssen, J., van der Marel, R. P., Gebhardt, K., Guhathakurta, P., Peterson, R. C., & Pryor, C. 2002, AJ, 124, 3270 —. 2003, AJ, 125, 376 Grandclément, P., Gourgoulhon, E., & Bonazzola, S. 2002, Phys. Rev. D, 65, 044021 Harris, W. E. 1996, AJ, 112, 1487 Heggie, D. C. 1975, MNRAS, 173, 729 Hills, J. G. 1990, AJ, 99, 979 Ivanova, N., Belczynski, K., Fregeau, J. M., & Rasio, F. A. 2005, MNRAS, 358, 572 Kalogera, V., King, A. R., & Rasio, F. A. 2004, ApJ, 601, L171 Kidder, L. E., Will, C. M., & Wiseman, A. G. 1993, Phys. Rev. D, 47, 3281 Kroupa, P., & Weidner, C. 2003, ApJ, 598, 1076 Kulkarni, S. R., Hut, P., & McMillan, S. 1993, Nature, 364, 421 Lee, H. M. 1995, MNRAS, 272, 605 —. 2001, Classical and Quantum Gravity, 18, 3977 Lee, M. H. 1993, ApJ, 418, 147 Merritt, D., Piatek, S., Portegies Zwart, S., & Hemsendorf, M. 2004, ApJ, 608, L25
Miller, M. C., & Colbert, E. J. M. 2004, International Journal of Modern Physics D, 13, 1 Miller, M. C., & Hamilton, D. P. 2002a, ApJ, 576, 894 —. 2002b, MNRAS, 330, 232 Murray, S. D., & Lin, D. N. C. 1996, ApJ, 467, 728 Peters, P. C. 1964, Physical Review, 136, 1224 Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724 Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17 —. 2002, ApJ, 576, 899 Quinlan, G. D., & Shapiro, S. L. 1989, ApJ, 343, 725 Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269 Sigurdsson, S., & Hernquist, L. 1993, Nature, 364, 423 Sigurdsson, S., & Phinney, E. S. 1993, ApJ, 415, 631 Spergel, D. N., Verde, L., Peiris, H. V., Komatsu, E., Nolta, M. R., Bennett, C. L., Halpern, M., Hinshaw, G., Jarosik, N., Kogut, A., Limon, M., Meyer, S. S., Page, L., Tucker, G. S., Weiland, J. L., Wollack, E., & Wright, E. L. 2003, ApJS, 148, 175 Spitzer, L. J. 1969, ApJ, 158, L139 Tutukov, A. V., & Yungelson, L. R. 1993, MNRAS, 260, 675 van der Marel, R. P. 2004, in Coevolution of Black Holes and Galaxies, 37 Watters, W. A., Joshi, K. J., & Rasio, F. A. 2000, ApJ, 539, 331 Wen, L. 2003, ApJ, 598, 419 White, N. E., & van Paradijs, J. 1996, ApJ, 473, L25 Willems, B., Henninger, M., Levin, T., Ivanova, N., Kalogera, V., McGhee, K., Timmes, F. X., & Fryer, C. L. 2005, ApJ, 625, 324 Woitas, J., Leinert, C., & Köhler, R. 2001, A&A, 376, 982
a In these simulations the cluster never reaches equipartition before a Hubble time. All numbers are calculated at the end of the simulation at log t = 10.13. 10 eq
N OTE. — Results of all simulations. Starting from the left the columns are the model name, number of total runs made, the average mass of the largest remaining BH at equipartition, the standard deviation of the largest mass, the total number of successive mergers the average largest remaining BH went through before equipartition, the largest BH formed that remained in the cluster, the fraction of simulations with a BH off mass > 100 M⊙ present at equipartition, the fraction of BHs ejected as singles ( fe,sin ), the fraction of BHs ejected in binaries ( fe,bin ), the log of time at which the cluster reached equipartition(teq ), the fraction of mergers in the cluster that happened in triples (Ntrip /Ncluster ), the fraction of total mergers that occurred in the cluster (Ncluster /Nmerge ), the average total number of mergers for a single simulation (Nmerge ), the expected Advanced LIGO rate of detection if the cluster formed when zform = 7.8 (R7.8 ), and the expected Advanced LIGO detection rate if the cluster formed when zform = 1 (R1 ). For our calculation of the expected detectable merger rate, we assume that the given cluster model has a current uniform number density ρ0 = 1 Mpc−3 and that the interferometer is capable of detecting a NS–NS inspiral up to DL,0 = 190 Mpc (see §4 for our detailed calculations and assumptions).