Multiple stellar-mass black holes in globular clusters: theoretical confirmation
aa r X i v : . [ a s t r o - ph . GA ] N ov Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 19 July 2018 (MN L A TEX style file v2.2)
Multiple stellar-mass black holes in globular clusters:theoretical confirmation
Anna C. Sippel ⋆ , Jarrod R. Hurley Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, VIC 3122, Australia draft 19 July 2018
ABSTRACT
While tens or hundreds of stellar-remnant black holes are expected to form in globularstar clusters, it is still unclear how many of those will be retained upon formation, andhow many will be ejected through subsequent dynamical interactions. No such blackholes have been found in any Milky Way globular cluster until the recent discovery ofstellar-mass black holes in the globular cluster M22 (NGC 6656) with now an estimatedpopulation of 5 −
100 black holes. We present a direct N -body model of a star clusterof the same absolute and dynamical age as M22. Imposing an initial retention fractionof ≈
10% for black holes, 16 stellar-remnant black holes are retained at a cluster ageof 12 Gyr, in agreement with the estimate for M22. Of those 16 BHs, two are in abinary system with a main sequence star each while also one pure black hole binary ispresent. We argue that multiple black holes can be present in any Milky Way clusterwith an extended core radius, such as M22 or the model presented here.
Key words: globular clusters: general - galaxies: star clusters: general - stars: mass-loss - methods: N -body simulations Stellar-mass black holes (BHs) are formed at the endpointsof stellar evolution of very massive stars. Depending on thechemical composition or metallicity of a star, a supernovaresulting in either a neutron star or a black hole is the ul-timate fate for most stars above ≈ − M ⊙ (Pols et al.1998).Tens or hundreds of such stellar-remnant BHs are ex-pected to form in globular clusters (GCs), however at thecurrent stage it is still not entirely clear how many will re-ceive a velocity kick at formation and get ejected imme-diately, and how many BHs might be removed during thesubsequent dyamical evolution of the cluster. While sev-eral extragalactic GCs containing BHs are known at thecurrent stage (e.g. Maccarone et al. 2011), none had beenconfirmed in a Milky Way GC until the recent discovery oftwo such stellar-mass BHs in M22 (NGC 6656, Strader et al.2012). Since only BHs currently undergoing observable ac-cretion can be detected via radio or X -ray emission andStrader et al. (2012) estimate that 2 −
40% of BHs are ex-pected to become members of binary systems with observ-able accretion over 10 Gyr, it is possible that a total popu-lation of ≈ −
100 BHs exists in M22.Theoretical predictions in the past have lead to the as-sumption that well before a cluster age of 12 Gyr all but up ⋆ E-mail: [email protected] to four (possibly all) BHs should be ejected from the clus-ter core (Sigurdsson & Hernquist 1993), or similarly thatnearly all BHs should be ejected, with the possibility of re-maining BHs capturing normal stars to form low-mass X-raybinaries in low-density environments (Kulkarni et al. 1993).Strader et al. (2012) claim their observations to be in con-trast to such prior theoretical predictions. However it wasnot possible to test this in direct N -body models compa-rable to clusters like M22 until very recently as such mod-els are computationally expensive. We note that multipleBHs have already been found to remain in N -body mod-els with smaller numbers of stars (e.g. Mackey et al. 2007,2008) as well as in Monte Carlos simulations of globularcluster evolution (e.g. Downing et al. 2010; Downing 2012),while BHs are ejected from the cluster in N -body modelsby Banerjee et al. (2010) where high numbers of BHs wereadded initially.Based on the recent findings by Repetto et al. (2012)suggesting that BHs receive similar velocity kicks as neu-tron stars (NSs) upon formation, and the observationsthat not all NSs are expected to receive a high kick(Pfahl et al. 2002), we present a direct N -body model for anintermediate-mass globular cluster containing 262 500 starsin total and incorporating a retention fraction of ≈ NBODY6 (Aarseth 1999, 2003) in a Milky Way-like gravitational po-tential field.While the model presented in this paper is by no means c (cid:13) Anna C. Sippel, Jarrod R. Hurley an attempt at a direct model of M22, it is representative ow-ing to comparable dynamical and absolute ages. We evolvethe model up to an age of 20 Gyr, and focus our analysis atthe age of 12 Gyr: the estimated age for many globular clus-ters in the Milky Way including M22 (Salaris & Weiss 2002).At this stage, ≈
180 000 stars are still retained within thecluster (including 16 BHs), implying that this is the largestdirect N -body model of a globular cluster currently available(cf. Baumgardt & Makino 2003; Hurley & Shara 2012). We use the direct N -body code NBODY6 (Aarseth 1999, 2003;Nitadori & Aarseth 2012) to evolve our model and chosea set-up similar to the clusters presented in Sippel et al.(2012). We evolve a star cluster with N i = 250 000 stel-lar systems and an initial binary fraction of b f = 0 . N =262 500 including 25 000 stars in 12 500 binary systems. Thestars are drawn from the initial mass function presented byKroupa et al. (1993) within the limits 0 . ≤ m ≤ M ⊙ and distributed according to a Plummer density profile(Plummer 1911; Aarseth et al. 1974). The initial clustermass adds up to M i = 1 . × M ⊙ . The conversion from N -body units to physical scale sizes leaves the scale radius as afree parameter which is set to give an initial half-mass radiusof r = 6 . d gc = 8 . r t ≈
50 pc for the N -body model, while M22 has a smallertidal radius of r t = 27 pc (Mackey & van den Bergh 2005)owing to the smaller galactocentric distance of d gc = 4 . .
001 (corresponding to [Fe/H] ≈ − . − . M ⊙ , where theresulting remnants are within the mass range of 4 − M ⊙ (Hurley et al. 2000; Belczynski et al. 2010). These massivestars evolve rapidly with most BHs forming during the first ≈
10 Myr of cluster evolution.
It is well established that some NSs will receive velocitykicks upon formation resulting from asymmetries in the starssurface before collapse (Pfahl et al. 2002). To mimic a neu-tron star retention fraction similar to indications from ob-servations of ≈
10% (c.f. Pfahl et al. 2002; Pfahl 2003) weadopt a random flat kick distribution between 0 −
100 km/s.
Figure 1.
Location of BHs (red crosses) at 12 Gyr of cluster age,projected in front of the cluster to increase visibility.
Recent indications suggest that BHs should receive simi-lar kicks to NSs (Repetto et al. 2012), as such we applyan identical kick distribution to BHs. For our model, theescape velocity is v esc = p GM/r = 14 km/s at thevery beginning of cluster evolution ( M is the cluster mass, r the half-mass radius and G the gravitational constant).At an age of 200 Myr, by which time all BHs and NSs willhave safely formed, the escape velocity is 10 . r = 8 . M = 11 . × M ⊙ .We note that in general, for a supernova occurring ina tight binary system, the kick velocity can be significantlyabove the cluster escape velocity and still result in one orboth binary components being retained, with the binaryloosened or broken up. Seven BHs form while part of a bi-nary system, all of them being disrupted either immediatelyor within the next ≈
10 Myr. Ultimately, this means thatout of ≈
450 stellar remnant BHs formed through stellarevolution, 48 are retained beyond the first 200 Myr, i.e. arenot immediately ejected from the cluster via the velocitykick process. Within less than one Gyr, they have settledinto the cluster core owing to mass-segregation and have asteadily decreasing half-mass radius of ≈ The sizes of N -body star cluster models are usually mea-sured by either the half-mass radius or the core radius. Bothquantities are measured in three dimensions, in contrast toobservational size scales. The half-mass radius r is sim-ply the 50% Lagrangian radius, however the N -body coreradius is a density weighted average distance from the clus-ter centre (von Hoerner 1960, 1963; Casertano & Hut 1985;Aarseth 2003) and hence not comparable to the observa-tional core radius resulting from e.g. King model fits to thesurface brightness profile. The procedure to calculate bothobservable half-light as well as core radii is similar to thatdescribed in Sippel et al. (2012). In particular, we cross con-volve the NBODY6 output of mass, luminosity and radius foreach single star with stellar atmosphere model calculations(Kurucz 1979) to obtain V-band magnitudes. Surface bright-ness profiles are produced by taking the average over pro-jection along the x − , y − and z − axes; and we use gridfit(McLaughlin et al. 2008) to fit a King model to the cluster. c (cid:13) , 000–000 tellar-remnant black holes in globular clusters t [Gyr] r [pc] Figure 2.
Size evolution of the cluster over time, from top tobottom these are half-mass r (dashed), N -body core radius r c (dotted), further Lagrangian radii r , r . and r . aswell as the half-mass radius for the BH population (red), alwayssmaller than even the 0 .
1% Lagrangian radius.
From projection effects alone, the half-mass radius is ≈ . The first few Gyr of cluster evolution are dominated byheavy mass loss arising from stellar evolution as well as dy-namical relaxation that results in an initial expansion of thecluster to a half-mass radius of up to r ≈
11 pc (Fig.2). The cluster then slowly starts to contract to r =10 . N -body core radius r tc = 3 . r h = 6 . r oc = 3 . ≈ . light radius to to calculate this number assuggested by Djorgovski (1993) and Harris (1996).Owing to the influence of the external galactic tidal fieldas well as dynamical interactions, stars are continuously lostfrom the cluster, with ≈
180 000 stars remaining at 12 Gyr.This corresponds to ≈
68% with respect to the initial num-ber of stars and a total mass of M = 6 . × M ⊙ ( ≈ . M ⊙ (compared toan overall average of 0 . M ⊙ ). Simultaneously, heating ofthe core from either BHs or NSs can be expected to producesignatures such as an expanded core (Mackey et al. 2007,2008), where both effects can influence the surface bright-ness profile and are measurable by means of the core radius(King 1966). M22 has a large core radius compared to otherglobular clusters in the Milky Way, indicating that a heatingsource may be present. While our model does not enter thephase of core collapse during its evolution, the subsystemof stars contained within the innermost 2 pc (dynamicallydominated by BHs) is entering a phase of core-collapse at ≈ t [Gyr] M [M sun ] mm MS m wd m BH Figure 3.
Contribution of stars to the central parsec of the clus-ter, i.e. d < d is again measured in threedimensions. It can be seen that at 12 Gyr, MS stars as well asWDs and BHs are contributing roughly 1 / ≤ M ⊙ for both together). The ratiois different when looking at the number of stars in each category:at 12 Gyr, MS stars count for roughly 2 / ≈ / ≈ . tion in the core radius (Fig. 2). This results from the factthat the BH population has a much smaller half-mass re-laxation time than the cluster as a whole, resulting the theBHs segregating to the cluster centre. The half-mass radiusof the entire BH population is ≈ . .
1% Lagrangian radius (see Fig. 2). At 12 Gyr, the to-tal mass contained in BHs is ≈ M ⊙ corresponding toonly ≈ .
3% of the overall cluster mass but ≈
32% of themass within the central parsec.This already indicates that BHs are an important com-ponent to the contribution of mass in the cluster centre. In-deed we find that at 12 Gyr, the innermost parsec (measuredin three dimensions) is composed equally of main sequencestars, white dwarfs, and BHs (Fig. 3). For the metallicity ofboth M22 and our model the main sequence turnoff mass is0 . M ⊙ at 12 Gyr, implying that most white dwarfs (witha peak in their mass distribution at 0 . M ⊙ ) will be lessmassive than the main sequence turnoff, while all neutronstars and black holes will be more massive. In total, 245 M ⊙ is contained in the innermost pc in 352 stars (compared to2 322 stars or 1 263 M ⊙ in projection), and this central den-sity increases to 528 stars or 423 M ⊙ at an age of 20 Gyr (or2 582 stars and 1 616 M ⊙ in projection). Binary formation and disruption, as well as hardening of bi-nary systems, are common processes in globular clusters.However it is easier to form new binary systems via ex-change interactions than creating them in two-body encoun- c (cid:13) , 000–000 Anna C. Sippel, Jarrod R. Hurley ters (Heggie 1975; Heggie & Hut 2003). A black hole can bepart of a binary system where the other component can beat any stellar evolution stage. In our simulation, the first dy-namical binary system including a BH forms at ≈
380 Myr(a 25 . M ⊙ BH and a 0 . M ⊙ MS star). At 12 Gyr, the ini-tial binary fraction of b f = 0 .
05 is slightly decreased to b f = 0 .
045 (with 7 850 binary systems present). In M22,one or potentially both observed BHs are in a binary witha low-mass main sequence star, or possibly a white dwarf(Strader et al. 2012). In our model, we find that at around12 Gyr, two BHs are in a long-period binary system with amain sequence star each, where neither of those systems areundergoing mass transfer. Specifically, system a) consists ofa BH with mass 6 . M ⊙ and a main sequence star with mass0 . M ⊙ while system b) consists of a BH with mass 6 . M ⊙ and a main sequence star of 0 . M ⊙ . Just before 12 Gyr athree-body encounter occurs involving system b) and a mainsequence star with mass 0 . M ⊙ , replacing the lighter mainsequence star and forming again a binary system consistingof a black hole and main sequence star.A BH-BH binary system is also present at around12 Gyr, where the components of this binary have massesof 22 M ⊙ each. Just after 12 Gyr, during a three-body en-counter with another BH with 13 M ⊙ , this binary gets dis-rupted, whilst the individual components of the three-bodysystem stay retained within the cluster. A BH-BH binaryin a new combination forms quickly and the total numberof 16 BHs is stable around the 12 Gyr timeframe, as illus-trated in Fig. 4. Dynamical interactions like these are whatcan ultimately cause the depletion of all but one BH or atight BH-BH binary (which might ultimately merge, Aarseth2012), from the cluster. However it is not certain over whichtimescale the evaporation of black holes from globular clus-ters should take place. This depletion timescale is dependenton random encounters, necessary to eject BHs, while thefrequency of such encounters depends on cluster parameterssuch as the core density (Pooley et al. 2003). We concludethat 12 Gyr is not enough time to deplete all or almost allBHs. In fact we have evolved the model up to 20 Gyr, andeven at this late stage, 10 BHs remain in the cluster. Fiveof those are in a binary system with either a main sequencestar, white dwarf or another black hole. We also find a BH-NS system with a total mass of ≈ M ⊙ at ≈
18 Gyr in astable configuration for ≈ . . M ⊙ in a three-bodyencounter. Clausen et al. (2012) have recently presented astudy regarding BH-NS mergers in GCs and we concludethat dynamical interactions don’t efficiently produce BH-NSbinaries in GCs. While a direct model approach for M22 is still out of reachfrom a computational point of view, the N -body model pre-sented in this work is comparable in many aspects. Whileour model has a half-mass relaxation time of t rh = 2 . t rh = 1 . r h = 3 . r oc = 1 . N -bodymodel is less dense than M22. We argue that the higher cen-tral density of M22 does not severly change our predictionsfor stellar remnant BHs in globular clusters with large core a) m [M sun ] maxmean t [Gyr] b) N Figure 4.
Evolution of the mass of the entire black hole pop-ulation a) and the number of BHs within the cluster b). In theupper panel, the solid black line denotes the maximum mass ofthe BHs still retained in the cluster, while the blue-dashed linedenotes the mean mass of the BH population at a given time. Thenumber of BHs is decreasing b) from 48 at 200 Myr to 16 BHs at12 Gyr of cluster age. A yellow background denotes the presenceof one BH-BH binary, while orange means two such systems arepresent, and white none. Note that the disruption of both BH-BH binary systems at ≈ ≈
11 Gyr is accompanied by theejection of the most massive black hole contained in the cluster.The disruption of one BH-BH binary is also accompanied by theejection of one BH at ≈ ≈ . radii since examples of denser model clusters with BHs existin the literature.Most prominent is the N -body model byHurley & Shara (2012) consisting of 200 000 stars ini-tially, evolved at a galactocentric distance of 3 . . ≈
35 000 stars remaining at 12 Gyr. How-ever, the half-mass relaxation time of the Hurley & Shara(2012) model is much below that of M22 resulting in themodel reaching the end of its initial phase of core-collapseat 11 . r oc = 0 .
84 pc with fourBHs still present (two of those in a BH-BH binary). Thisis the case even though the central density is significantlyenhanced at this time (just over 3 000 stars in the innermostparsec, or ≈ N -bodymodel presented in this work – however much closer to thenew model here.Other examples of N -body cluster models with BHscurrently exist in the literature (Mackey et al. 2007, 2008;Hurley & Mackey 2010 just to name a few), as well asMonte Carlo models of more massive globular clustersat various densities (Downing et al. 2010; Downing 2012;Morscher et al. 2012). In addition we note that models in-volving mergers of stellar remnant black holes as well as wellas intermediate mass BHs have been carried out in the past c (cid:13) , 000–000 tellar-remnant black holes in globular clusters (e.g. Portegies Zwart & McMillan 2002; Baumgardt et al.2004; O’Leary et al. 2006; Aarseth 2012). We conclude that all three findings of Strader et al. (2012):a) the observation of two stellar remnant black holes in M22which are b) potentially in binary systems with main se-quence stars and that c) 5 −
100 stellar remnant black holesmight be present in a cluster such as M22 – are in excellentagreement with our direct N -body model. A key point inthis analysis is that M22 has a large core radius compared tothe average for the globular cluster population of the MilkyWay, and is not dynamically old. We conclude that multi-ple stellar remnant black holes in the core of such clusterscan exist even beyond the current time. Furthermore, in the N -body model presented here we find that 1 / We thank Juan Madrid, Marie Martig and Sverre Aarsethfor comments and discussion. The N -body model wasevolved for ≈
250 days on a NVIDIA Tesla C1060 graphicscard at Swinburne University of Technology and the dataanalysis was carried out on the Green Machine at Swin-burne. AS acknowledges a SUPRA scholarship from Swin-burne University.
REFERENCES
Aarseth S., 2012, ArXiv e-printsAarseth S. J., 1999, PASP, 111, 1333Aarseth S. J., 2003, Gravitational N-Body SimulationsAarseth S. J., Henon M., Wielen R., 1974, A&A, 37, 183Banerjee S., Baumgardt H., Kroupa P., 2010, MNRAS, 402,371Baumgardt H., Makino J., 2003, MNRAS, 340, 227Baumgardt H., Makino J., Ebisuzaki T., 2004, ApJ, 613,1143Belczynski K., Bulik T., Fryer C. L., Ruiter A., ValsecchiF., Vink J. S., Hurley J. R., 2010, ApJ, 714, 1217Casertano S., Hut P., 1985, ApJ, 298, 80Clausen D., Sigurdsson S., Chernoff D. F., 2012, ArXiv e-printsDjorgovski S., 1993, in Djorgovski S. G., Meylan G., eds,Structure and Dynamics of Globular Clusters Vol. 50 ofAstronomical Society of the Pacific Conference Series,Physical Parameters of Galactic Globular Clusters. p. 373Downing J. M. B., 2012, arXiv:1204:5363Downing J. M. B., Benacquista M. J., Giersz M., SpurzemR., 2010, MNRAS, 407, 1946Harris W. E., 1996, AJ, 112, 1487 Heggie D., Hut P., 2003, The Gravitational Million-BodyProblem: A Multidisciplinary Approach to Star ClusterDynamicsHeggie D. C., 1975, MNRAS, 173, 729Hurley J. R., Mackey A. D., 2010, MNRAS, 408, 2353Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315,543Hurley J. R., Shara M. M., 2012, MNRAS, 425, 2872Hurley J. R., Tout C. A., Aarseth S. J., Pols O. R., 2004,MNRAS, 355, 1207Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329,897King I. R., 1966, AJ, 71, 276Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262,545Kulkarni S. R., Hut P., McMillan S., 1993, Nat, 364, 421Kurucz R. L., 1979, ApJS, 40, 1Maccarone T. J., Kundu A., Zepf S. E., Rhode K. L., 2011,MNRAS, 410, 1655Mackey A. D., van den Bergh S., 2005, MNRAS, 360, 631Mackey A. D., Wilkinson M. I., Davies M. B., GilmoreG. F., 2007, MNRAS, 379, L40Mackey A. D., Wilkinson M. I., Davies M. B., GilmoreG. F., 2008, MNRAS, 386, 65McLaughlin D. E., Barmby P., Harris W. E., Forbes D. A.,Harris G. L. H., 2008, MNRAS, 384, 563Miyamoto M., Nagai R., 1975, Publ. Astron. Soc. Japan,27, 533Morscher M., Umbreit S., Farr W. M., Rasio F. A., 2012,ArXiv e-printsNitadori K., Aarseth S. J., 2012, arXiv:1205.1222O’Leary R. M., Rasio F. A., Fregeau J. M., Ivanova N.,O’Shaughnessy R., 2006, ApJ, 637, 937Pfahl E., 2003, in KITP Conference: Globular Clusters:Formation, Evolution and the Role of Compact ObjectsTotal Numbers and the NS Retention ProblemPfahl E., Rappaport S., Podsiadlowski P., 2002, ApJ, 573,283Plummer H. C., 1911, MNRAS, 71, 460Pols O. R., Schroder K.-P., Hurley J. R., Tout C. A., Eggle-ton P. P., 1998, MNRAS, 298, 525Pooley D., Lewin W. H. G., Anderson S. F., BaumgardtH., Filippenko A. V., Gaensler B. M., Homer L., Hut P.,Kaspi V. M., Makino J., Margon B., McMillan S., Porte-gies Zwart S., van der Klis M., Verbunt F., 2003, ApJ,591, L131Portegies Zwart S. F., McMillan S. L. W., 2002, ApJ, 576,899Repetto S., Davies M. B., Sigurdsson S., 2012, MNRAS,425, 2799Salaris M., Weiss A., 2002, A&A, 388, 492Sigurdsson S., Hernquist L., 1993, Nat, 364, 423Sippel A. C., Hurley J. R., Madrid J. P., Harris W. E.,2012, ArXiv e-printsStrader J., Chomiuk L., Maccarone T. J., Miller-JonesJ. C. A., Seth A. C., 2012, Nat, 490, 71von Hoerner S., 1960, Z. Astrophys., 50, 184von Hoerner S., 1963, Z. Astrophys., 57, 47 c (cid:13)000