Binary Black Hole Mergers from Globular Clusters: Masses, Merger Rates, and the Impact of Stellar Evolution

# Binary Black Hole Mergers from Globular Clusters: Masses, Merger Rates, and the Impact of Stellar Evolution

## Abstract

The recent discovery of GW150914, the binary black hole merger detected by Advanced LIGO, has the potential to revolutionize observational astrophysics. But to fully utilize this new window into the universe, we must compare these new observations to detailed models of binary black hole formation throughout cosmic time. Expanding upon our previous work Rodriguez et al. (2015), we study merging binary black holes formed in globular clusters using our Monte Carlo approach to stellar dynamics. We have created a new set of 52 cluster models with different masses, metallicities, and radii to fully characterize the binary black hole merger rate. These models include all the relevant dynamical processes (such as two-body relaxation, strong encounters, and three-body binary formation) and agree well with detailed direct -body simulations. In addition, we have enhanced our stellar evolution algorithms with updated metallicity-dependent stellar wind and supernova prescriptions, allowing us to compare our results directly to the most recent population synthesis predictions for merger rates from isolated binary evolution. We explore the relationship between a cluster’s global properties and the population of binary black holes that it produces. In particular, we derive a numerically calibrated relationship between the merger times of ejected black hole binaries and a cluster’s mass and radius. With our improved treatment of stellar evolution, we find that globular clusters can produce a significant population of massive black hole binaries that merge in the local universe. We explore the masses and mass ratios of these binaries as a function of redshift, and find a merger rate of Gpcyr in the local universe, with 80% of sources having total masses from to . Under standard assumptions, approximately 1 out of every 7 binary black hole mergers in the local universe will have originated in a globular cluster, but we also explore the sensitivity of this result to different assumptions for binary stellar evolution. If black holes were born with significant natal kicks, comparable to those of neutron stars, then the merger rate of binary black holes from globular clusters would be comparable to that from the field, with approximately of mergers originating in clusters. Finally we point out that population synthesis results for the field may also be modified by dynamical interactions of binaries taking place in dense star clusters which, unlike globular clusters, dissolved before the present day.

1

## I Introduction

With the recent detection of the binary black hole merger GW150914 Abbott et al. (2016a), the era of gravitational-wave (GW) astronomy is here. Advanced LIGO promises to open an entirely new window into the universe, performing detailed measurements of cataclysmic events such as the mergers of compact objects. In particular, the mergers of binary neutron stars (NS) and binary black holes (BHs) will help answer detailed questions about such diverse fields as stellar evolution, gamma-ray bursts, the NS equation-of-state, and even general relativity itself Read et al. (2009); Cornish et al. (2011); Nissanke et al. (2010); Stevenson et al. (2015). But to fully realize the science potential of these instruments, we must be prepared to relate GW observations to detailed models of the compact object merger rate throughout the universe.

To that end, significant work has been done to translate our understanding of star formation and stellar evolution into realistic predictions for compact object merger rates. Population synthesis codes (e.g. Portegies Zwart et al., 2001; Spera et al., 2015; Belczynski et al., 2002, 2008) use simple recipes for stellar evolution to rapidly model the evolution of a large population of stellar binaries. However, this approach suffers from significant uncertainties. Only weak observational constraints for binary NS systems exist from observations of binary pulsars in the Milky Way Kalogera et al. (2004), while the best rate estimates for stellar-mass binary BH (BBH) mergers arise from only one GW detection (GW150914) and one less-significant GW trigger (LVT151012) Abbott et al. (2016b). The prescriptions used to predict such systems from population synthesis methods are tuned to observations of stellar wind-mass loss rates and low-mass X-ray binaries Belczynski et al. (2010a); Dominik et al. (2012) which are few and far between, particularly for the low-metallicity massive stars that will potentially dominate the BBH merger rate Belczynski et al. (2010b). Furthermore, uncertainties in the physics of binary stellar evolution make it difficult to constrain the BBH merger rate from isolated binaries to within several orders of magnitude.

In Rodriguez et al. (2015), we explored the contribution to the BBH merger rate from globular clusters (GCs), a population of dense, old stellar systems observed in most galaxies. BBH formation in GCs is primarily driven by dynamics, with three-body binary formation and strong gravitational encounters creating the majority of BBH systems. Since these clusters have survived to the present day, they do not suffer from the uncertainties in star formation history, and present a population of systems with well-measured masses and metallicities. Furthermore, because most BBHs in GCs are formed dynamically, they are less sensitive to the physics of binary stellar evolution. We analyzed a collection of GC models from Morscher et al. (2015), and found that a single Advanced LIGO observatory could detect 10 to 100 BBH mergers from GCs per year, similar to many estimates of the merger rate from the field; however, these models assumed a distribution of BH masses that made it difficult to directly compare our results to population synthesis studies.

In this paper, we explore the full range of masses, mass ratios, eccentricities, and merger rates of BBH systems created in GCs. To compare our results to recent estimates from the field, we have upgraded the stellar evolution algorithms of our dynamical models with new prescriptions for the metallicity-dependent stellar winds of O/B stars Vink et al. (2001) and a new prescription for the remnant and fallback masses of compact objects after supernova Fryer et al. (2012). These changes bring our code into agreement with the StarTrack population synthesis code Belczynski et al. (2002, 2008), and allow us to create GC models with the same BH mass distribution used recent estimates for the merger rate from the galactic field Dominik et al. (2012, 2013, 2015). With these new stellar-wind prescriptions, we now find that GCs can produce significant numbers of massive BBHs that merge in the local universe, including BBHs, similar to GW150914.

In Section II, we describe the basics of our Monte Carlo approach to modeling dense star clusters, our choice of initial conditions, and the upgrades to our stellar evolution algorithms. In Section III, we explore how the global properties of a GC determine the masses, inspiral times, and merger rates of the BBH population. Then in Section IV, we show the distributions of merging BBH masses, mass ratios, and eccentricities as a function of redshift. We also explore the effects of various assumptions in our stellar evolution algorithm, such as differences in the common-envelope physics, the supernova mechanism, and the natal kicks given to BHs at formation. Finally, in Section V, we calculate the merger rates for BBHs from GCs as a function of redshift, and compare these rates to the most recent estimates of merger rates from the field Dominik et al. (2013). Throughout this paper, we assume cosmological parameters of , , and km s, to compare with the most recent estimates of BBH mergers from galactic fields Dominik et al. (2013).

## Ii Cluster Models

We create 52 GC models with our Cluster Monte Carlo code (CMC), an orbit-sampling Hénon-style Monte Carlo code Hénon (1971a, b). CMC has been developed and described in several previous papers Joshi et al. (2000, 2001); Fregeau et al. (2003); Fregeau and Rasio (2007); Chatterjee et al. (2010); Pattabiraman et al. (2013). By assuming that the dynamics of a GC is primarily driven by two-body relaxation between neighboring particles, CMC can create models of dense star clusters on a significantly faster timescale than a traditional direct -body integrator. This allows us to fully explore the parameter space of GCs, including the massive and compact clusters used in this study. Recent comparisons to state-of-the-art -body simulations Wang et al. (2015) have shown that CMC can correctly model the evolution of realistic clusters with particles over 12 Gyr, including the masses, semi-major axes, and eccentricities of the ejected BBHs Rodriguez et al. (2016).

In addition to two-body relaxation, CMC incorporates all the relevant physics for treating the BBH formation and merger problem. This includes:

• binary-single and binary-binary gravitational encounters computed explicitly with the Fewbody small- integrator Fregeau et al. (2004),

• three-body binary formation, implemented with an analytic prescription Morscher et al. (2013), and

• single and binary stellar evolution, implemented with BSE Hurley et al. (2000, 2002), and improved with the stellar remnant and BH kick prescription from Kiel and Hurley (2009); Chatterjee et al. (2010). For this paper, we also update the prescriptions for stellar winds and supernova fallback, in order to replicate the BH mass distribution of Belczynski et al. (2010b); Dominik et al. (2012, 2013, 2015) (see Section II.2)

The BBH merger time is found by directly integrating the orbit-averaged Peters equations Peters (1964). For mergers in the cluster, this is calculated by BSE. For mergers of BBHs that were ejected, we calculate the inspiral time from the masses, semi-major axis, and eccentricity at ejection. The merger time is then the sum of the inspiral time and the time at which the binary is ejected from the cluster.

### ii.1 Initial Conditions

For our main result, we create 52 GC models spanning a range of different realistic parameters. These are described in Table 1. We consider a 3x2x4 grid of models, with different metallicities, virial radii, and initial particle numbers, for a total of 24 initial conditions, and four additional models in which we varying the assumptions of stellar evolution. For each of the 24 main models, we generate two statistically independent cluster models for each initial condition. We consider three metallicities: , corresponding roughly to the lower-end of the Milky Way GC (MWGC) metallicity distribution Harris (1996), , corresponding to the peak of the MWGC distribution, and , corresponding to the peak of the high-metallicity MWGCs. We will refer to our and models as low-metallicity, and our model as high-metallicity2. We place the , , and clusters at galactocentric distances of 2 kpc, 8 kpc, and 20 kpc respectively, as there exists a strong correlation between a cluster’s galactocentric distance and its metallicity Djorgovski and Meylan (1994). We also consider clusters with initial virial radii of both and . This is motivated by both observational evidence of cluster formation in the local universe, suggesting that is typical Scheepmaker et al. (2007) and theoretical modeling suggesting that smaller produce GCs with properties more similar to observed MWGCs Morscher et al. (2015).

For each metallicity, galactocentric distance, and virial radius, we then consider clusters with , , , and particles. The particles are placed at radii drawn from a King profile King (1966) with (as (Morscher et al., 2015; Rodriguez et al., 2015) found no correlation between the initial GC concentration and either the observational properties of the clusters or the properties of the BBH population). The stellar masses are drawn from a Kroupa initial mass function (IMF) of the form Kroupa (2001):

 P(m)dm∝m−αdm (1)

with

 α={1.30.08M⊙≤M<0.5M⊙2.30.5M⊙≤M (2)

We consider masses from to , in order to directly compare our results to recent results for the field (e.g. Dominik et al. (2012, 2013)). We also assume that all models begin with 10% of the particles as binaries. This is accomplished by randomly selecting 10% of particles and adding a companion mass drawn from a flat distribution in the mass ratio. The binary separations are drawn from a distribution, with a lower-limit near the point of stellar contact and an upper-limit such that the velocity of the binary components is equal to the average velocity of a star in the cluster core. The eccentricities are chosen from a thermal distribution, Heggie (1975). Because we select 10% of our particles to become binaries, our models have 10% more stars than the number of particles (e.g. our model contains stars).

In Rodriguez et al. (2015), we used a series of models from Morscher et al. (2015) to estimate the BBH merger rate. We noted that in our simulations, dynamically-formed BBHs from GCs were characteristically more massive than those formed by stellar evolution of binaries in the field; however, this result was highly dependent on our particular prescription for binary stellar evolution. More recent work Belczynski et al. (2010b); Dominik et al. (2012); Spera et al. (2015) has suggested that massive BBHs (with total masses up to 160 ) can be formed in low-metallicity environments, where the reduced stellar winds from massive stars enable the formation of BHs with masses up to . Furthermore, the reduced wind-driven mass loss increases the chance that the system will remain bound during its evolution to become a BBH.

We have upgraded the binary stellar evolution module in CMC with improved prescriptions for the metallicity-dependent stellar winds and supernova-driven mass loss. For stellar winds, we make two additions to the default BSE implementation: first, for O and B stars, we implement the mass-loss prescriptions developed in Vink et al. (2001), which we will refer to as the “Vink prescription”. This fit determines the mass-loss as a function of the star’s mass, effective temperature, metallicity, and luminosity for any hot () hydrogen-rich star on the main sequence. The Vink prescription covers separately the temperature range from and , while explicitly excluding the range between and , which is complicated by the appearance of Fe ion line-driven winds at approximately . We follow the prescription from Belczynski et al. (2010a), extending the lower temperature prescription (equation (25) from Vink et al. (2001)) up to and the higher temperature prescription (equation (24) from Vink et al. (2001)) down to . In addition to the new prescription for O and B stars, we add a metallicity dependence to the evolution of naked-helium (Wolf-Rayet) stars. This is done by supplementing the original BSE prescription Hamann and Koesterke (1998) with the metallicity dependence from Vink and de Koter (2005).

We also include two new prescriptions for the supernova mechanism. These prescriptions, first developed in Fryer et al. (2012), describe the amount of material that falls back onto the newly-formed compact object for neutrino-driven and convection-enhanced supernovae. They consider two cases, based on the delay between the core bounce and the explosion: the rapid case, which assumes that any explosion occurs within the first 0.25 seconds after the core bounce, and the delayed case, which relaxes this assumption. The rapid prescription replicates the observed “mass-gap” between NSs and BHs Farr et al. (2011); Belczynski et al. (2012), while the delayed prescription allows for BH masses from to . For the default in this study, we use the rapid supernova model, although in Section IV.4 we also explore the effects of the delayed supernova mechanism on BBH production in GCs

In addition to determining the mass of the supernova remnant, the mass of the fallback material determines the velocity kick the newly-formed BH receives from the explosion. Consistent with the original version of BSE, we give all BHs natal kicks drawn from a Maxwellian with a dispersion of km s Hobbs et al. (2005) as is done for NSs formed by core-collapse supernova. However, we assume that the final velocity is lowered by the mass of the fallback material according to

 Vfinal=(1−ffb)Vnatal (3)

where is the fraction of the ejected supernova mass that will fall back onto the newly-formed proto-compact object, determined by equations (16) and (19) in Fryer et al. (2012). For any sufficiently massive BH progenitor (), the fallback completely damps any natal kick, and the BH is retained in the cluster. This is consistent with the direct collapse scenario discussed in Fryer and Kalogera (2001).

The goal of these modifications is to replicate the BH mass spectrum used in the recent studies of field-born BBHs, specifically those employed in the StarTrack population synthesis code. With these new changes to our version of BSE, we find that CMC produces a similar relationship between the zero-age main-sequence mass of a massive star and its final BH mass (Figure 11 of Fryer et al. (2012)). Therefore we can now explicitly compare our results to those from Dominik et al. (2013, 2015). Note that this does not imply that our stellar evolution code would replicate the StarTrack results for binary stars, as we have not modified our binary stellar evolution prescription from the version used in previous papers. However, since the majority of BBHs from GCs are formed dynamically (90% of all BBH mergers, and 99.7% of mergers at ), this does not significantly impact our results. To confirm this, we create an additional 4 GC models varying the physics of binary stellar evolution, and explore the effects of these assumptions in Section V.1.

## Iii BBH Production in GCs

### iii.1 Binary Properties

One of the most interesting results of Rodriguez et al. (2015) was the significant increase in the merger rate of BBHs from GCs as compared to other studies. This was attributable to two key effects: first, the Monte Carlo method allows us to model significantly larger clusters than a direct -body method. Since CMC can integrate clusters with in days to weeks, Rodriguez et al. (2015) was able to explore the contribution of massive GCs that were shown to dominate the BBH merger rate. Previous studies, particularly those using direct -body integrators, were limited to clusters with at most particles Portegies Zwart and McMillan (2000); Banerjee et al. (2010); Tanikawa (2013); Bae et al. (2014). After 12 Gyr of evolution, these clusters would produce at most 10 to 20 BBHs, only of which would merge in a Hubble time. The one exception was the two studies Downing et al. (2010, 2011), which used a Monte Carlo approach similar to CMC. However, they only considered a single, relatively low cluster mass () and did not include direct integration of binary-single and binary-binary encounters, critical for correctly predicting the properties of dynamically-formed BBHs. By covering the full distribution of GC masses, our models can correctly predict the number of binaries ejected from the most massive GCs, significantly increasing the total number of binaries produced by GCs overall.

Not only do massive clusters produce more BBHs, but the BBHs that they produce are more likely to merge within 12 Gyr than BBHs from lower-mass GCs. This result arises from the physics of binary formation and interaction in a cluster environment. A hard binary (with binding energy greater than the typical kinetic energy of particles in the cluster), typically undergoes a series of strong encounters with other single and binary stars. The result of these repeated encounters, known as “Heggie’s law”, is that hard binaries tend to get harder, increasing their binding energy on average with each encounter, and transferring that energy to the center-of-mass velocity of the interacting stars Heggie (1975). The recoil velocity of the binary in the cluster is increased after each encounter, until eventually the binary receives a velocity kick sufficient to eject it from the cluster.

Since conservation of energy demands that the change in be proportional to the post-encounter kinetic energy of the binary, and the potential of the cluster determines the minimum kinetic energy for escaping particles, we can use the cluster potential to determine the minimum binding energy (and maximum semi-major axis) required to eject a binary. First, we assume that increases by roughly the same fraction after every binary-single encounter, and that of that energy gets transferred into the binary recoil. Heggie (1975); Hut et al. (1992) both assumed an average increase of 40% after each encounter, while Portegies Zwart and McMillan (2000) assumed 20%; however, both studies considered the average only for particles with equal masses. For generality, we simply assume that the binding energy is related to the center-of-mass kinetic energy by some factor , such that after a scattering encounter

 ECM=αEbin=αGm1m22a (4)

where is the binary semi-major axis and and are the component masses. If we assume the cluster potential is represented by a Plummer distribution, then the energy required to eject a binary from the center of the cluster is given by (e.g., Heggie and Hut, 2003)

 Eesc=(m1+m2)|Φc|=GMGC(m1+m2)√22/3−1Rh (5)

where is the cluster mass and is the half-mass radius (the radius enclosing half the cluster mass). Combining (4) and (5), we can relate the semi-major axis of ejected binaries to the mass and half-mass radius of the cluster at the time of ejection by

 MGCRh =(α√22/3−12)μbina =1κμbina (6)

where is the reduced mass of the binary and we have defined the coefficient .

Equation (6) reinforces an important point: when averaged over many encounters, it is the global properties of the cluster, not the micro-physics of the dynamics, that primarily determines the semi-major axes of binaries formed dynamically. We show the semi-major axes, inspiral times, and redshift delays for the binaries from each of our models in Figure 1. Although this relationship between the cluster potential and the semi-major axis of ejected binaries has beeen noted before Portegies Zwart and McMillan (2000); Moody and Sigurdsson (2009), the size of the models considered here allows us to extend this analysis to realistically-sized clusters with realistic BBH masses.

Since our GC models span a range of BH masses, cluster masses, and cluster half-mass radii, we can compare (6) directly to the BBHs ejected from our models. In Figure 2, we show the relationship between the semi-major axis and reduced mass of each binary and the half-mass radius and total mass of the cluster from which it was ejected. We show the distribution of from all ejected binaries in the top panel. We find that the distribution of values roughly follow a log-normal distribution. In the bottom panel, we plot the value of for every binary against at the time it was ejected. We also show a line corresponding to the median of the log-normal fit. We note that this relationship only applies to binaries ejected from the cluster during strong encounters with a single object. Binaries ejected during encounters with another binary have significantly more complicated interactions, including a much greater possibility of exchanging components. However, we find that 81% the of binaries ejected from our GC models are ejected immediately following a binary-single encounter (involving either a hardening encounter or an exchange of componenets), while only 13% are ejected following a binary-binary encounter. For completeness, we also show these points in Figure 2, though we exclude them from the fit to equation (6). With the fit to from equation (6), we can express the probability for a BBH to be ejected with a given eccentricity and semi-major axis as:

 P(e)de=2ede (7) P(a|MGC,Rh,μbin)da=1aσ√2π× (8) exp⎡⎢ ⎢ ⎢ ⎢ ⎢⎣−(logμbinRhaMGC−a∗)22σ2⎤⎥ ⎥ ⎥ ⎥ ⎥⎦da

where and are the parameters of the log-normal distribution with mean and . We reiterate that equations (7) and (8) only apply for the 81% of sources which are ejected from the cluster following a binary-single encounter.

With a relationship between the cluster parameters and the semi-major axes of the binaries it ejects, we can show how cluster dynamics determines the inspiral times of BBHs once they are ejected from the cluster. From the Peters equations Peters (1964), we see that the merger time of a binary scales as

 tinsp∝a4m1m2(m1+m2) (9)

If we combine equation (9) with the scaling from equation (6), we can show that the inspiral time of a binary with total mass scales as

 tinsp∝(RhMGC)4M (10)

where we have assumed all binaries to have equal-mass components.

In Figure 3, we show the inspiral times for BBHs ejected from each of our GC models as a function of binary total mass. What is immediately striking is that the median inspiral time appears to decrease with increasing binary mass, in contrast to the scaling derived in equation (10). This is primarily due to the influence of the cluster itself: while the binary mass does determine the inspiral time of the binary post-ejection, it is the cluster mass and half-mass radius that determine the separation at ejection, and the factor in equation (10) increases by several orders-of-magnitude over the 12 Gyr lifetime of a GC.

In Figure 4 we show the median mass of ejected BBHs and the ratio for a single GC model over time. As the cluster ages, it preferentially ejects its most massive BHs first, working its way from most to least massive BHs as time progresses. As the BHs are ejected, the cluster expands and loses mass, increasing and significantly increasing the inspiral time of the lower-mass binaries ejected at late times. This scaling drives the counter-intuitive decrease in inspiral times for high-mass BBHs seen in Figure 3: these BBHs are ejected early in the cluster lifetime, when is low. As the cluster ages, increases significantly, and the average inspiral time for an ejected BBH increases accordingly.

### iii.2 BBH Mergers per Cluster

To determine the mean number of mergers per cluster, we need a functional form for the number of inspirals as a function of cluster mass. In Rodriguez et al. (2015), this was done with a simple linear regression, assuming that the mass of the cluster at 12 Gyr was proportional to the number of inspirals it had produced. However, as we have shown, this scaling is somewhat more complicated: while the mass of a cluster may be proportional to the number of BBHs it produces, the fraction of those sources that will merge in a Hubble time, , is also controlled by the cluster mass. While the number of BBHs a cluster forms should linearly depend on its mass (more stars yield more BHs, yielding more BBHs), the relationship between and the final GC mass is less obvious.

Rather than derive a physically-motivated functional form for , we choose a simpler approach, and require only that whatever we choose increases with mass and assymptotes to 1 as . As this requirement also describes the family of cumulative probability distributions, we elect to use an error function (the cumulative normal distribution). This yields a final relationship for the number of inspirals as a function of GC mass:

 Ninsp(MGC)=NBBH(MGC)×finsp(MGC) (11)

where

 NBBH(MGC) =aMGC+b (12) finsp(MGC) =erf(MGC−M02σ) (13)

We fit and separately, then take the product to be . We fit equation (11) for high-metallicity and low-metallicity clusters. We also consider separate fits to clusters with and to better understand the impact of cluster size on our results (Figure 5).

### iii.3 Sampling our Inspiral distribution

In order to create a representative population of BBHs from GCs, we must compare our models to observations of GC systems. In Rodriguez et al. (2015), this was accomplished by integrating over the GC mass function (GCMF), which we assumed to be a log-normal distribution with mean and width , based on recent observations of the GC luminosity function in Harris et al. (2014) and an assumed mass-to-light ratio of 2 in solar units Bell et al. (2003). This, when combined with the spatial density of GCs per comoving volume (Rodriguez et al., 2015, Supplemental Materials), yielded an effective estimate for the density of BBH mergers from GCs in the universe. This was then converted into a detection rate by multiplying this density by the distribution of inspirals from our models in redshift, reweighed to favor GC models which more closely resembled the MWGC distribution. This statistical machinery was used to account for the fact that the models developed in Morscher et al. (2015) and used in Rodriguez et al. (2015) were designed to probe the space of GC initial conditions, not to span the space of observed GCs.

We proceed in a similar fashion. We integrate equation (11) over the GCMF from 0 to to determine the mean number of inspirals per GC. Note that we do not consider a varying upper-mass cutoff in the integral, as the upper-mass cutoff does not strongly affect the mean (Rodriguez et al., 2015, and its associated erratum). We do this separately for high- and low-metallicity clusters, though we use the same GCMF for both. We present the results in Table 2. To combine the different metallicity results, we use the computed spatial density of high-metallicity () and low-metallicity clusters () from the Supplemental Materials of Rodriguez et al. (2015), which implies that 56% of all clusters are low-metallicity.

In addition to , we need to weigh the models to more closely represent the observed mass distribution of GCs. To accomplish this, we compute the final mass for our , and models. We then use the average mass for each set of simulations to divide the GCMF into discrete bins, such that each set of models sits in the center of its respective mass bin. Each model is then assigned a weight corresponding to the total integrated number of observed GCs in that mass bin. This is done separately for high- and low-metallicity clusters, with 4 bins for the low-metallicity cluster and 3 for high-metallicity clusters3. In addition to the contribution from clusters with different masses, we have also assumed that 56% of clusters in the universe have low metallicities. However, of our models have low metallicity, so we multiply the weights of the low-metallicity clusters by 0.86, to ensure our inspirals are correctly sampled according to the metallicity spread of observed GCs.

We use these weights to create a collection of BBH mergers representative of a full population of GCs. This is accomplished by randomly drawing samples from each model according to its weight. As an example, the high-metallicity model has the largest weight of all our models, so we draw all of its BBH inspirals for our merger population. On the other hand, the low-metallicity model has a weight 0.46 times the weight of the largest model, so we only draw 46% of its binaries.

For the next section, we will use this collection of inspirals for our primary analysis. For any scatter plots showing a representative sample of mergers, we only use a single draw from our models. For percentiles, rate computations, and any quantities involving relative quantities, we perform 10 independent draws from our GC models. This effectively considers each BBH merger from each of our models, weighted by the contribution of that model to the total merger rate.

## Iv BBH Properties

### iv.1 Mass Distributions

With the Vink prescription for stellar winds, our models now produce significantly more massive BHs than was found in previous GC studies. While Rodriguez et al. (2015) found a maximum BH mass of and for high and low-metallicity systems respectively, the reduced stellar winds for low-metallicity massive stars significantly increase these cutoffs. We now find maximum BH masses of , , and for our , , and models respectively (see Figure 6). This is consistent with the latest results from population synthesis codes, such as StarTrack Belczynski et al. (2010a) and SEVN Spera et al. (2015) (though the latter, which explores significantly higher-mass progenitors than we consider here, can produce BHs as massive as for stars with ).

Despite the significant changes to the BH mass spectrum, the behavior of BHs in clusters, developed in Morscher et al. (2015), remains unchanged: after core collapse, the most massive BHs segregate into the center of the cluster, where they immediately form BBHs via three-body encounters. These BBHs then undergo a series of binary-single and binary-binary encounters, increasing their binding energies and shrinking their semi-major axes. Eventually, the recoil from one of these encounters will be sufficient to eject the binary from the cluster, as discussed in Section III.1. Although a significant number of binaries merge in the cluster (), the majority of these in-cluster inspirals occur early in the GC lifetime. At , only 0.06% of binary mergers (one merger from all 48 models) occur in-cluster. Of the ejected sources merging in the local universe, 99.7% were formed dynamically, which we define to be either a BBH formed from three isolated BHs by a three-body interaction, or a BBH formed from a primordial binary which swapped components at least once during a binary-single or binary-binary encounter.

In Figure 7, we show the masses for each of the inspirals from the weighted sample of GC BBH mergers. We break the masses down into two categories: source masses, or the local masses of each BBH, and observed masses, which correspond to the redshifted mass, , measured by an observer on Earth. We also show separate panels for the chirp mass of the source, , the total mass of the source, and the individual components of each binary.

The overall structure of the plots agrees well with our understanding of BH and BBH evolution in GCs: after the formation and core collapse of the cluster (at ), the most massive BHs form binaries and are ejected immediately. The GC processes through its available population of BHs, working its way through the BH population from most to least massive, so that only low-mass BHs () are still present in massive GCs by the present day. In the total-mass panel of Figure 7, this story is obvious. The majority of the most massive inspirals (total mass ) merge soon after GC formation, between and . After these BBHs are ejected, the cluster moves on to less massive BBHs with longer inspiral times. These sources (with total masses from to ) form the predominant population of BBHs detectable in the present day.

The plateaus in the chirp mass, total mass, and component mass distributions are primarily determined by the maximum BH mass at each metallicity, which is in turn determined by the wind-driven mass loss from the Vink prescription (Figure 6). For the highest metallicity models (), this yields a large population of BHs, which in turn forms a large population of equal-mass BBHs with a total mass of . For the clusters, this yields a smaller collection of sources with a total mass of . However, for the lowest-metallicity models () there is no apparent collection of sources at as might be expected.

This behavior can again be explained by the wind-driven mass loss. Each model begins with an identical distribution of stars drawn from (1). For the highest-metallicity models, the mass-loss from these winds brings all stars with birth masses from to down to a final progenitor mass of to before the supernova occurs. Essentially, this truncates a large section of the high-mass end of the IMF to a single BH mass; however, for lower-metallicity models, the decreased efficiency of the stellar winds means a lower number of high-mass stars are being converted into maximum-mass BHs, essentially spreading out the high-mass stars over a wider range of BH masses. The number of maximum-mass BHs between each of our models decreases by roughly a factor of 5 between each of our metallicitiy bins. This also yields a higher number of inspirals with unequal-mass components for lower-metallicity models, which we discuss in the next section.

We also show in Figure 7 the masses and 90% uncertainties associated with the recent detection of GW150914. The masses reported from the parameter estimation The LIGO Scientific Collaboration and the Virgo Collaboration (2016a) of GW150914 are consistent with the masses of BBH mergers from GCs in the local universe, with the specific masses ( and ) being easily formed in lower-metallicity GCs. We also show the less-significant GW trigger LVT151012, with the 90% uncertainties taken from The LIGO Scientific Collaboration and the Virgo Collaboration (2016b). Although this signal was not claimed as a detection, we note that the masses of this event (total mass of ) is also consistent with BBH mergers from GCs in the local universe.

In Figure 8, we convert our scatter plot of BBH mergers into percentiles as a function of redshift. In the local universe (), we find that the median BBH source from a GC has a chirp mass of , with 50% of sources lying in the range , 80% of sources in the range , and 90% of sources in the range . This corresponds to a median total mass of , with a 50% interval of , a 80% interval of , and a 90% interval of . Given the analysis in Section III.1 and Figure 3, we can expect that massive binaries that merge in the local universe are more likely to be ejected from clusters with lower masses and larger half-mass radii. This is supported by our results: 77% of all BBH mergers are formed in our most massive GCs (). However, in the local universe, this fraction changes as a function of mass. 86% of binary mergers with total masses at are preferentially formed in high-mass clusters. On the other hand, for binary mergers with masses above , this fraction drops to 66%. We note that the large contribution from high-mass GCs is a result of our focus on GCs which have survived to the present day. Were we to consider lower mass clusters that disrupted before 12 Gyr, it is likely that the merger rate for massive BBHs in the local universe would increase significantly.

Finally, we note that there exists a small population of BHs significantly more massive than the maximum-BH mass allowed by stellar evolution. These BHs are the result of repeated stellar mergers early in the evolutionary history of the cluster. These collisions are primarily the mergers of massive giant and main-sequence stars during binary-single and binary-binary encounters, with a smaller number resulting from direct collisions of single stars. We note that 3 of these sources were formed from a BBH merger whose remnant remained in the cluster, forming a second BBH; however, as CMC does not incorporate the physics of gravitational-wave recoil, these 3 sources are most likely nonphysical, as typical recoil velocities from a BBH merger greatly exceed the escape velocity of a typical GC Campanelli et al. (2007).

### iv.2 Mass Ratio Distributions

In addition to the individual BH masses, the mass ratios of BBHs may provide an important clue to the formation mechanism of a binary. In Figure 9, we show the mass ratios from our models as a function of chirp mass and total mass. The mass ratios pile up near unity, with a median mass ratio of 0.87 and 68% of sources having mass ratios greater than 0.8. This is to be expected: binary-single and binary-binary scattering experiments show that binaries in dense stellar environments are prone to swapping components, preferentially ejected less-massive components in favor of more massive companions Sigurdsson and Hernquist (1993).

However, we also note that there exists a small population () of sources with mass ratios less than 0.6. These appear at the base of the apparent “V” formations in the top panels of Figure 9, corresponding to binaries with total masses of , , , and . These features are a direct result of the peaks in the distribution of BH masses (Figure 6). The three most-obvious peaks in the BH mass distribution, around , , and , create peaks of equal-mass binaries at total masses near , , and . For instance, the maximum BH mass for the models creates a large population of equal-mass binaries with a total mass of . However, these BHs will also form binaries with less-massive companions from the BH distribution (roughly down to the second peak of the distribution at ). This creates the feature running from a mass ratio of 1 at to a mass ratio of 0.4 at . We reiterate that these features are the result of the peaks in the BH mass distribution, which strongly depends on the stellar metallicity. As we have only considered three metallicities here, it is likely that a collection of models fully spanning the distribution of GC metallicities would find a more continuous distribution of mass ratios.

### iv.3 Eccentricity Distributions

It is a well-known result that dynamically-formed binaries will follow a thermal distribution of eccentricities Heggie (1975). This distribution, , is the result of the thermalization of velocities that occurs through repeated encounters in the dense cluster core. However, once ejected, these binaries evolve in isolation, and are rapidly circularized due to the preferential emission of gravitational radiation at periapsis Peters (1964). Given the difficulties involved in detecting such binaries and the possibility of parameter estimation for such systems Huerta et al. (2014), it is important to quantify the number of sources that will enter the LIGO detection band with non-negligible eccentricities.

In Figure 10, we show the eccentricities of all ejected binaries. The eccentricity at 10Hz is computed by integrating Peters (1964) from ejection to , the Keplerian separation for a binary with an orbital period of Hz (corresponding to a gravitational-wave frequency of 10Hz). We show the relationship between the 10Hz eccentricity and the semi-major axis and eccentricity at ejection in the upper and middle plots, and the distribution of final eccentricities in the lower plot.

This result suggests that the majority of BBHs from GCs will very nearly circular by the time they are detectable, with 99% of sources entering the LIGO band with eccentricities below . Thus, eccentric systems ejected from clusters will not be a significant source for Advanced LIGO. However, it has been shown that a non-negligible number of sources in the cluster will be driven to merger by long-term secular effects, such as the Lidov-Kozai mechanism, in BH triple systems, which we do not consider here. Were that physics properly included in these models, it is likely that 1% of mergers could be detected with significant eccentricity Antonini et al. (2016).

### iv.4 The Impact of Stellar Evolution

The main uncertainty in population synthesis estimates of the BBH merger rate from the field concerns the physics of binary stellar evolution. The most recent studies (e.g., Dominik et al., 2012) have considered a large number of possible outcomes from various stages of binary stellar evolution. These different models are then used to constrain the uncertainties, although the sheer number of parameters makes it difficult to constrain the merger rate even to within several orders-of-magnitude.

More recent work Belczynski et al. (2015) has focused on three specific questions: can BBH progenitors survive the common-envelope (CE) phase when the donor star is in the Hertzsprung gap (HG)? What is the standard mechanism for BH-producing supernova (rapid or delayed)? And what is the impact on the rate if BHs are born with natal kicks drawn from the same distribution as NSs? The standard StarTrack model assumes that stars in the HG cannot be CE donors, that supernovae are well described by the rapid model, and that BH natal kicks are reduced proportionally to the mass of fallback material, as in equation (3).

In contrast, our models for dynamical formation of BBHs are largely insensitive to the uncertainties in binary stellar evolution. Since the majority of BBHs are formed dynamically, their numbers, semi-major axes, eccentricities, and inspiral times are determined by well-understood gravitational physics. Any changes to our binary stellar evolution prescription should not have a significant effect on our BBH population. However, the same cannot be said for single star evolution: any changes that modify the number or masses of single BHs will significantly change both results from both clusters and the field.

To quantify these effects, we rerun our , , and model, varying some of the assumptions in binary and single star evolution. We consider 5 models:

• the standard model, adopted in the previous sections. We use the standard BSE prescription for the CE evolution, employing the formalism4 Webbink (1984), with calculated by a fitting formula similar to that described in Claeys et al. (2014). By default, BSE will allow binaries to survive the CE phase so long as the core of either star does not fill its Roche lobe (including stars in the HG, for which BSE employs a time-dependent density profile Hurley et al. (2002)).

• low CE binding energy model, in which we set for all binaries,

• high CE binding energy model, in which we set ,

• delayed supernova model, in which we use the delayed supernova model for BH masses and fallback, as described in Fryer et al. (2012), and

• full BH natal kicks model, in which we ignore the effects of fallback, and give each BH its full kick velocity from the NS distribution.

We show the impact of these different cases on the distribution of BBH mergers in Figure 11. As expected, the mass spectrum of BHs does not significantly depend on the physics of stellar evolution, with one notable exception for the full BH natal kicks case. Over 12 Gyr of evolution, the standard model produces 84 mergers. The case produces 72 mergers, the case produces 87 mergers, the delayed supernova case produces 79 mergers, but the full kicks case produces only 7 mergers.

This significant decrease in BBH production is to be expected: the distribution of BH velocities in the high-kick case is drawn from a Maxwellian with km s, meaning most BHs will be born with speeds significantly greater than the escape speed from the center of the cluster ( km s). In that case, the majority of BHs are ejected before they have a chance to dynamically form binaries. We explore the implications of this for the BBH merger rate in Section V.1.

## V Cosmological Merger Rates

To compare different models for BBH formation, it is necessary to understand the predicted merger rate for each different formation scenario. We compute the cosmological merger rate of BBH inspirals from our models as a function of redshift. We then assume all GCs to be 12 Gyr old, and we describe the merger rate per unit time and comoving volume as

 Rs(z)=ρGC⟨Ninsp⟩P(z) (15)

where is the spatial density of GCs (Rodriguez et al., 2015, Supplemental Materials), is the mean number of GC inspirals per cluster from Section III.2, and is the normalized merger rate at a given redshift, defined as

 P(z)=Pt(tlookback(z)) (16)

where is the probability distribution of inspiral merger times, computed by generating a kernel density estimate of the merger times from our sample of BBH mergers, and is the cosmological lookback time at a given redshift Hogg (1999). Note that is the distribution of sources in time at redshift , not the distribution of sources in redshift5.

Additionally, we are interested in the total number of sources that merge within a comoving volume out to a given redshift, . We compute the observed rate of mergers per unit time as

 Ro(z)=∫z0dVcdz′Rs(z′)(dtsdto)dz′ (17)

where is the comoving volume at a given redshift, is the source merger rate (15), and is rate of time dilation between a clock measuring the merger rate at redshift and a clock on Earth (e.g., Belczynski et al., 2014).

To gain a handle on the uncertainties in our assumptions, we consider three different cases: our standard case assumes a spatial density Mpc, and that half of those clusters have virial radii of , and half have virial radii of . We also consider an optimistic case, which assumes Mpc (the optimistic estimate from Rodriguez et al. (2015)) and for all clusters, and a conservative case, which assumes Mpc, and for all clusters. The standard, optimistic, and conservative cases are shown in the left panels of Figure 12. The upper-left panel shows the observed cumulative rate as a function of redshift (equation 17), while the lower-left panel shows the source rate (equation 15). In the local universe, our standard assumptions leads to a merger rate of 5 Gpcyr while our optimistic and conservative assumptions predict merger rates of 20 Gpcyr and 2 Gpcyr, respectively. This is consistent with the merger rate for GW150914-like binaries, which was found to be Gpc Abbott et al. (2016b).

In the right panels of Figure 12, we show the merger rates for binary inspirals with specific total masses from our standard case. We compute separately for each mass bin by fitting the inspirals in each bin to equation (11) and integrating over the GCMF; this is the same technique used for the total rate, although for mass bins with only a handful of inspirals, such as the 5 inspirals in the bin, we abandon equation (11) in favor of a simple linear regression. We then compute the distribution of inspirals, , separately for each mass bin, and compute the rates with equations (15) and (17). In the local universe, we can expect a merger rate of Gpcyr for BBHs with a total mass of , and Gpcyr for , with higher masses contributing minimally to the total merger rate. We note that the oscillatory behavior of the high-mass bins at low redshifts is the result of the small number of events. As an example, the bump at for sources is the result of a single inspiral whose BH progenitor underwent repeated mergers. Given the uncertainties in modeling stellar collisions and mergers, particularly with the Monte Carlo method, the merger rates for such rare events should be treated with appropriate skepticism.

### v.1 Comparison to the Field

One of the basic questions for gravitational-wave astrophysics is whether or not GW observations can discriminate between different formation physics for observed sources. Can intrinsic parameters, such as the BH masses and spin, be used to determine whether a particular event was formed by dynamics, or by isolated binary stellar evolution? And can a collection of GW observations be used to answer detailed questions about star formation and stellar evolution Stevenson et al. (2015)? While each of these questions is worthy of study, in this paper we focus on a much simpler, but more critical question: what fraction of detectable sources originated in GCs?

To that end, we make the following assumptions: we assume that all BBH mergers in the universe arise from either isolated binary evolution in the field (as done in Dominik et al., 2012, 2013, 2015) or from GCs. We also assume that the star formation scenarios for the field and for old GCs are distinct, with all GCs being formed in a single burst of star formation 12 Gyr ago, while star formation in the field occurring continuously, with metallicity increasing with the age of the universe. The StarTrack models in Dominik et al. (2013) employed a star-formation rate from Strolger et al. (2004) and a metallicity-redshift relationship found by averaging the results of Pei et al. (1999) and Young and Fryer (2007); however, these assumptions produced unphysically-high metallicities for star formation in the local universe (with metallicities as high as ), so they considered two extreme scenarios: a high-end metallicity scenario, where they divide their metallicity profile by 1.7, to agree with results from Yuan et al. (2013), and a low-end metallicity scenario, where they divide the profile by 3.0, to agree with SDSS observations from Panter et al. (2008).

We compare the cumulative merger rates from our models to the cumulative rates for the StarTrack models examined in Dominik et al. (2013) in Figure 13. The four models considered are:

• the standard model, in which they assume all binaries which enter the CE with a donor star in the HG will immediately merge, use the rapid supernova prescription, and employ a physical description taken from Xu and Li (2010),

• the optimistic model, where they allow BBH progenitors to survive the CE with donor stars in the HG,

• the delayed supernova model, where they use the delayed supernova model, and

• the full BH natal kick model, where they give newly formed BHs natal kicks identical to those of NSs, regardless of the fallback fraction.

Given that our merger rates are insensitive to the assumptions of binary stellar evolution, we assume that our main results will remain unchanged under the assumptions of the standard, optimistic, and delayed supernova cases. For the case of high BH natal kicks, we multiply our standard merger rates by 0.08, proportional to the decreased number of mergers from our high-kick model (Section IV.4). This oversimplifies the true relationship between BH retention and natal kicks, since the fraction of retained BHs will increase for more compact, massive clusters with higher escape speeds, but it is sufficient for this current estimate. We also compare the GC merger rate to the rate from the field for both the high-end and low-end metallicity evolution scenarios. The fraction of mergers from GCs up to a given redshift is defined to be

 FGC(z)=RGCo(z)RGCo(z)+RFieldo(z) (18)

where is from equation (17), and is the results from Dominik et al. (2013) obtained from the Synthetic Universe website 6.

In the local universe, the standard model suggests that 15% (12%) of merging BBHs will have been formed in a GC under the high-end (low-end) metallicity scenarios. For the optimistic model, the field rate increases significantly, with dropping to 4% (2%) in the optimistic case. For the models using the delayed supernova prescription, the fraction increases to 25% (15%) in the local universe.

However, for the high-BH natal kick model, the merger rate from BBHs is comparable to those from the field, with reaching values of 55% (45%) in the local universe. This surprising result is a direct implication of the dynamical formation scenario. It was previously assumed that, were BHs to be formed with large natal kicks, the majority of BHs would be ejected from clusters, with only a small fraction (the low-velocity tail of the kick distribution) being retained and processed into BBHs. Given that, it is surprising that rates from GCs should rival those from isolated binary stellar evolution in the field.

To understand this, assume that there is some small fraction of BHs, , born with sufficiently low kicks that they remain gravitationally bound to whatever system they are member of, be it a GC or a binary star system. For a GC, this means that after 20 Myr, approximately of the total BHs will be retained in the cluster. If we then assume that some fraction, of these BHs will be dynamically processed into BBHs, then the total number of BBHs produced by the cluster is . However, for systems in the field, each BBH must be formed from a binary progenitor, and therefore each binary must survive two natal kicks. From an initial population of binaries, this suggests that the number of BBHs produced by the field in this scenario is . In other words, if the BH natal kicks are inversely proportional to the fraction of retained BHs, then as the kicks are increased, the rates from GCs decrease as , while the rates from the field decrease as .

Additionally, the retention of BHs in our high-kick GC model is aided by the primordial binary fraction. Of the 27 BHs retained after supernova, 18 were in binaries at the time of formation. As such, these systems are bound by both the local gravitational potential of their companions and the full potential of the cluster. We note that our choice of initial binary fraction () is much lower than the fraction assumed by Dominik et al. (2013). As a quick check, we rerun our high-kick GC model with a binary fraction of 50%. We find that the high-kick GC model now retains 48 BHs initially (41 of which were in binaries at formation), and creates 16 BBH mergers over its 12 Gyr lifetime. Assuming these results scale with the total merger rate of BBHs from GCs, then roughly 75% (70%) of all BBH mergers in the local universe would have formed in GCs. This approach is similar to proposed mechanisms for the retention of NSs in GCs Pfahl et al. (2002).

Three caveats must be mentioned at this point. First, it is unlikely that clusters, and particularly the massive GCs studied here, were formed with a binary fraction near 50%. Theoretical studies have shown that the fraction of binaries in a GC remains roughly constant over time Hurley et al. (2007); Ivanova et al. (2005), while observations of GCs suggest a binary fraction between 1% and 5% (e.g. NGC 6397, Davis et al., 2008). Second, increasing the binary fraction significantly changes the evolution of the GC as the primordial binaries heat the cluster core (e.g. Chatterjee et al., 2013). And finally, such a dependence on the binary fraction will also introduce a significant dependence on the initial conditions for binaries in star clusters, such as the distribution of initial mass ratios and semi-major axes. We explore the effects of our choice of initial conditions in Chatterjee et al. (2016).

### v.2 What is the field?

In several papers examining the BBH merger rate of the field, it has been assumed that the field and star clusters exist as separate populations, formed by fundamentally different processes However, observational evidence suggests that the majority of star formation occurs in clusters Lada and Lada (2003). These regions of star formation, consisting of clusters as small as , might dominate the star formation rate, particularly at high redshifts Kruijssen (2015). Furthermore, observational evidence of the peculiar velocities of O stars in the MW suggest that of these massive stars formed in clusters de Wit et al. (2005). If that is the case, then a clean division between “field” and “dynamical” sources could be a significant oversimplification. Studies of young star clusters () have shown that the properties of a BBH population can be significantly altered by dynamics, even for small clusters that disrupt on a short ( Myr) timescale Ziosi et al. (2014).

A proper analysis of the contribution from clusters of all masses to the BBH merger rate is beyond the scope of this paper; however, it is important to ask the question: what fraction of the merger rate of field binaries computed in population synthesis studies is actually affected by dynamics, which these studies neglect entirely? As an order-of-magnitude check, let us assume that, consistent with our IMF, approximately one of every 600 stars will become a BH. If we assume that, at the very least, a cluster must contain two BHs for dynamics to play a role, then combined with an average stellar mass of , the minimum cluster mass where dynamics could be considered is . If the fraction of star formation that occurs in clusters at high redshifts is (c.f. Figure 9, Kruijssen, 2015), and the cluster initial mass function follows a distribution, then, integrating from to suggests that of BBHs are formed in regions where dynamics may play a significant role.

Such an increase could play a significant role in the merger rate of dynamically-formed BBHs. The results presented here have considered only the present-day population of GCs which, at , constitute only of the baryonic matter in the MW Portegies Zwart et al. (2010). If most galaxies formed with of stars in dynamically-relevant clusters, then a significant fraction of the BBH population may have been affected by dynamics at some point.

However, this does not directly imply that the BBH merger rate would be dominated by clusters. The majority of these disrupted clusters would have masses significantly below the population of old, massive GCs studied here. And since we have shown that less-massive clusters not only produce fewer BBHs, but a lower fraction of BBHs that will merge in a Hubble time, it is not immediately obvious whether an early population of disrupted clusters will significantly affect the merger rates. This is consistent with Ziosi et al. (2014), which found that in young star clusters, only 0.3% of BBHs merged within 12 Gyr.

## Vi Conclusion

In this paper, we have explored mergers of BBHs from GCs. Using our cluster Monte Carlo code, CMC, we have created a broad range of GC models with different masses, metallicities, and initial virial radii, designed to approximate the distribution of GCs observed in the present-day universe. To increase the realism of our models, and to facilitate an easier comparison to BBH merger rates from the field (such as those from Dominik et al. (2013)), we have upgraded the stellar evolution prescription used by CMC, with new temperature-dependent stellar winds for O and B stars and new prescriptions for the supernova mechanism. These changes were designed to bring our stellar evolution subroutines into agreement with those currently employed in the StarTrack population synthesis code. With these enhancements, CMC now reproduces the mass distribution for single BHs that was used in the most recent estimates of BBH merger rates from the field Dominik et al. (2012, 2013, 2015).

By considering a wide range of models, we were able to broadly characterize the relationship between the global properties of a GC (its mass and radius), and the inspiral times of the binaries it creates. We showed that, in addition to creating more BBHs, more massive and more compact clusters ejected BBH binaries with higher binding energies and smaller semi-major axes, leading a greater number of BBHs to merge within 12 Gyr. This explained the significant increase in merger rates first reported in Rodriguez et al. (2015): by accurately modeling the median and high-mass end of the GCMF, our realistic cluster models produced significantly more merging BBHs than previous studies of this type. After integrating over the full GCMF, we found that the average GC will have produced 260 BBH mergers throughout its 12 Gyr lifetime.

We then used the GCMF to select a population of inspiral events representative of the population of GCs in the present-day universe. With the new prescriptions for wind-driven mass loss and the rapid supernova mechanism, our GC models can now form BBHs with total masses from to . These BBHs follow the same story that we first observed in Morscher et al. (2015); Rodriguez et al. (2015): the most-massive BHs lead the first period of mass segregation, driving the deep collapse of the cluster core. These BHs then dynamically form BBHs, which are among the first to be ejected from the cluster, and the first to merge. The GC processes through its BH population from most to least massive, continuing to eject BBHs up to the present day. As such, the total mass of the merging BBHs decreases with redshift, with a median BBH total mass in the local universe of and 50% of sources lying between and . However, there exists a significant tail of massive sources, ejected with large semi-major axes (and correspondingly large inspiral times), such that the 90% interval of masses extends from up to . We found that these massive BBHs were more likely to be formed in lower mass GCs, as these clusters ejected binaries with wider separations and longer inspiral times.

We also found that the distribution of mass ratios depended strongly on the BH mass spectrum. Any peaks in the distribution of single BH masses created corresponding peaks at twice that mass in our BBH population, where these BHs formed a large population of equal-mass BBHs. As such, the three most prominent peaks in the BH mass distribution (and 18, 30, and 53) created corresponding peaks in the distribution of BBH total masses at 36, 60, and 106. In between these regions, binaries would tend to form with unequal mass components, often drawing one of their components from the peaks of the BH mass spectrum. This created an overabundance of sources with only one component from the peak of the BH mass spectrum, creating certain regions of the BBH total mass distribution (20, 45, and 80 ) that favored systems with unequal component masses.

We then computed the merger rates from this population of BBHs as a function of redshift. In the local universe, we found that BBHs formed in GCs will merge at a rate of Gpc yr. Under highly optimistic assumptions, this rate becomes Gpc yr, while highly pessimistic assumptions forces the rate down to Gpc yr. For the standard assumptions, we also found a merger rate of Gpc yr for sources with total masses , and Gpc yr for sources with total masses . This is consistent with the merger rate of Gpc Abbott et al. (2016b) for binaries similar to GW150914. When comparing these numbers to estimates of BBH merger rates from the field, we found that roughly 15% of BBH mergers in the local universe will have originated in a GC. However, if BHs are born with large natal kicks, similar to NSs, this fraction increased significantly, with GCs possibly dominating the merger rate for BBHs. In addition to changing BBH merger rates and properties, any such changes effecting BH retention in early GCs (such as the natal kicks or the high-mass slope of the IMF) will significantly change the observational properties of GCs in the present day. A study to fully characterize the relationship between BH retention and GC observational properties is currently underway Chatterjee et al. (2016).

We reiterate that we have studied a conceptually clean, but not physically complete, picture of star formation and BBH formation in the universe. By using the present-day population of observed GCs, we have restricted ourselves to studying a well-posed problem with good observational constraints. However, both the current study and studies of field populations have ignored the fact that most star formation is assumed to occur in clusters, with lower-mass clusters being disrupted into what we today call “the field”. The merger rates and BBH properties of estimates from the field should only be applied to binary systems from star forming regions that were disrupted before they could be dynamically altered; however, it is not at all obvious that the results presented here can simply be scaled down to lower-mass star forming regions. A full exploration of the BBH merger problem will require an exploration of the intermediate regime, where most BBHs are formed from primordial binaries that undergo significant dynamical perturbations (e.g. Ziosi et al. (2014)). Fully understanding these three scenarios–the field, dynamical, and intermediate BBH formation regimes–will be critical to unlocking the potential of gravitational-wave astrophysics.

Finally, we have shown that the masses and redshift associated with the recent detection of GW150914 (and the less-significant GW trigger LVT151012) are consistent with dynamical formation in the low-metallicity environment of a GC as reported in our models. However, we will address the full implications of this detection in the context of GC modeling in a future work Rodriguez et al. (2016). This detection emphasizes the importance of BBH modeling as we transition into the era of GW astrophysics.

###### Acknowledgements.
We would like to thank Jarrod Hurley, Chris Belczynski, Vicky Kalogera, Onno Pols, Francesca Valsecchi, and Cole Miller for useful discussions. CR was supported by an NSF GRFP Fellowship, award DGE-0824162. This work was supported by NSF Grant AST-1312945 and NASA Grant NNX14AP92G.

### Footnotes

1. preprint: APS/123-QED
2. GCs in most galaxies are observed to fall into two distinct metallicity groups: high-metallicity, or “red” GCs, which are typically found in the galactic bulge, and low-metallicity, or “blue” GCs, typically found in the galactic halo Harris (2010).
3. We do not use the high-metallicity models, as these disrupted before 12 Gyr, and we are interested only in the population of clusters which survives to the present day.
4. The term arises from the energy balance between the gravitational potential energy of the envelope and the orbital energy of the binary: (14) where is the change in binding energy of the binary before and after the CE, is the mass of the donor star, is the mass of the donor star’s envelope, and is the donor star’s radius at Roche lobe overflow. In this model parameterizes the binding energy of the envelope, while parameterizes the efficiency of energy transfer from the binary’s shrinking orbit to the envelope.
5. By computing the kernel density estimate with the inspiral times of the bianries, the probability of seeing an inspiral at time is . By converting this to the distribution in redshift with (16), we ensure that is in units of equal time, not equal redshift. This is different than the approach employed in Rodriguez et al. (2015), and we note that properly accounting for the distribution increases the detection rates from that study by .
6. The merger rates from Dominik et al. (2013) and several other studies using StarTrack are available from the Synthetic Universe website (http://www.syntheticuniverse.org/)

### References

1. C. L. Rodriguez, M. Morscher, B. Pattabiraman, S. Chatterjee, C.-J. Haster,  and F. A. Rasio, Phys. Rev. Lett. 115, 051101 (2015).
2. B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso,  and R. X. Adhikari et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016a).
3. J. S. Read, C. Markakis, M. Shibata, K. Uryū, J. D. E. Creighton,  and J. L. Friedman, Phys. Rev. D 79, 124033 (2009).
4. N. Cornish, L. Sampson, N. Yunes,  and F. Pretorius, Phys. Rev. D 84, 062003 (2011).
5. S. Nissanke, D. E. Holz, S. A. Hughes, N. Dalal,  and J. L. Sievers, Astrophys. J. 725, 496 (2010).
6. S. Stevenson, F. Ohme,  and S. Fairhurst, Astrophys. J. 810, 58 (2015).
7. S. F. Portegies Zwart, S. L. W. McMillan, P. Hut,  and J. Makino, Mon. Not. R. Astron. Soc. 321, 199 (2001).
8. M. Spera, M. Mapelli,  and A. Bressan, Mon. Not. R. Astron. Soc. 451, 4086 (2015).
9. K. Belczynski, V. Kalogera,  and T. Bulik, Astrophys. J. 572, 407 (2002).
10. K. Belczynski, V. Kalogera, F. A. Rasio, R. E. Taam, A. Zezas, T. Bulik, T. J. Maccarone,  and N. Ivanova, Astrophys. J. Suppl. Ser. 174, 223 (2008).
11. V. Kalogera, C. Kim, D. R. Lorimer, M. Burgay, N. D’Amico, A. Possenti, R. N. Manchester, A. G. Lyne, B. C. Joshi,  and M. A. McLaughlin et al., Astrophys. J. 601, L179 (2004).
12. B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso,  and R. X. Adhikari et al.,  , 16 (2016b)arXiv:1602.03842 .
13. K. Belczynski, T. Bulik, C. L. Fryer, A. Ruiter, F. Valsecchi, J. S. Vink,  and J. R. Hurley, Astrophys. J. 714, 1217 (2010a).
14. M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel,  and R. O’Shaughnessy, Astrophys. J. 759, 52 (2012).
15. K. Belczynski, M. Dominik, T. Bulik, R. OâShaughnessy, C. Fryer,  and D. E. Holz, Astrophys. J. 715, L138 (2010b).
16. M. Morscher, B. Pattabiraman, C. Rodriguez, F. A. Rasio,  and S. Umbreit, Astrophys. J. 800, 9 (2015).
17. J. S. Vink, A. de Koter,  and H. J. G. L. M. Lamers, Astron. Astrophys. 369, 574 (2001).
18. C. L. Fryer, K. Belczynski, G. Wiktorowicz, M. Dominik, V. Kalogera,  and D. E. Holz, Astrophys. J. 749, 91 (2012).
19. M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel,  and R. O’Shaughnessy, Astrophys. J. 779, 72 (2013).
20. M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. E. Holz, T. Bulik,  and F. Pannarale, Astrophys. J. 806, 263 (2015).
21. M. Hénon, Astrophy. Space Sci. 13, 284 (1971a).
22. M. Hénon, Astrophy. Space Sci. 14, 151 (1971b).
23. K. J. Joshi, F. A. Rasio, S. P. Zwart,  and S. Portegies Zwart, Astrophys. J. 540, 969 (2000).
24. K. J. Joshi, C. P. Nave,  and F. A. Rasio, Astrophys. J. 550, 691 (2001).
25. J. M. Fregeau, M. A. Gurkan, K. J. Joshi,  and F. A. Rasio, Astrophys. J. 593 (2003).
26. J. M. Fregeau and F. A. Rasio, Astrophys. J. 658, 1047 (2007).
27. S. Chatterjee, J. M. Fregeau, S. Umbreit,  and F. A. Rasio, Astrophys. J. 719, 915 (2010).
28. B. Pattabiraman, S. Umbreit, W.-k. Liao, A. Choudhary, V. Kalogera, G. Memik,  and F. A. Rasio, Astrophys. J. Suppl. Ser. 204, 15 (2013).
29. L. Wang, R. Spurzem, S. Aarseth, K. Nitadori, P. Berczik, M. B. N. Kouwenhoven,  and T. Naab, Mon. Not. R. Astron. Soc. 450, 4070 (2015).
30. C. L. Rodriguez, M. Morscher, L. Wang, S. Chatterjee, F. A. Rasio,  and R. Spurzem, eprint arXiv:1601.04227  (2016).
31. J. M. Fregeau, P. Cheung, S. F. Portegies Zwart,  and F. A. Rasio, Mon. Not. R. Astron. Soc. 352, 1 (2004).
32. M. Morscher, S. Umbreit, W. M. Farr,  and F. A. Rasio, Astrophys. J. 763, L15 (2013).
33. J. R. Hurley, O. R. Pols,  and C. A. Tout, Mon. Not. R. Astron. Soc. 315, 543 (2000).
34. J. R. Hurley, C. A. Tout,  and O. R. Pols, Mon. Not. R. Astron. Soc. 329, 897 (2002).
35. P. D. Kiel and J. R. Hurley, Mon. Not. R. Astron. Soc. 395, 2326 (2009).
36. P. Peters, Phys. Rev. 136, B1224 (1964).
37. W. E. Harris, Astron. J. 112, 1487 (1996).
38. W. E. Harris, Philos. Trans. R. Soc. London, Ser. A 368, 889 (2010).
39. S. Djorgovski and G. Meylan, Astron. J 108, 1292 (1994).
40. R. A. Scheepmaker, M. R. Haas, M. Gieles, N. Bastian, S. S. Larsen,  and H. J. G. L. M. Lamers, Astron. Astrophys. 469, 925 (2007).
41. I. R. King, Astron. J 71, 64 (1966).
42. P. Kroupa, Mon. Not. R. Astron. Soc. 322, 231 (2001).
43. D. C. Heggie, Mon. Not. R. Astron. Soc. 173, 729 (1975).
44. W. R. Hamann and L. Koesterke, Astron. Astrophys. 335, 1003 (1998).
45. J. S. Vink and A. de Koter, Astron. Astrophys. 442, 587 (2005).
46. W. M. Farr, N. Sravan, A. Cantrell, L. Kreidberg, C. D. Bailyn, I. Mandel,  and V. Kalogera, Astrophys. J. 741, 103 (2011).
47. K. Belczynski, G. Wiktorowicz, C. L. Fryer, D. E. Holz,  and V. Kalogera, Astrophys. J. 757, 91 (2012).
48. G. Hobbs, D. R. Lorimer, A. G. Lyne,  and M. Kramer, Mon. Not. R. Astron. Soc. 360, 974 (2005).
49. C. L. Fryer and V. Kalogera, Astrophys. J. 554, 548 (2001).
50. S. F. Portegies Zwart and S. L. W. McMillan, Astrophys. J. Lett. 528, L17 (2000)astro-ph/9910061 .
51. S. Banerjee, H. Baumgardt,  and P. Kroupa, Mon. Not. R. Astron. Soc. 402, 371 (2010).
52. A. Tanikawa, Mon. Not. R. Astron. Soc. 435, 1358 (2013).
53. Y.-B. Bae, C. Kim,  and H. M. Lee, Mon. Not. R. Astron. Soc. 440, 2714 (2014).
54. J. M. B. Downing, M. J. Benacquista, M. Giersz,  and R. Spurzem, Mon. Not. R. Astron. Soc. 407, 1946 (2010).
55. J. M. B. Downing, M. J. Benacquista, M. Giersz,  and R. Spurzem, Mon. Not. R. Astron. Soc. 416, 133 (2011).
56. P. Hut, S. McMillan,  and R. W. Romani, Astrophys. J. 389, 527 (1992).
57. D. Heggie and P. Hut, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics (Cambridge University Press, 2003).
58. K. Moody and S. Sigurdsson, Astrophys. J. 690, 1370 (2009).
59. W. E. Harris, W. Morningstar, O. Y. Gnedin, H. O’Halloran, J. P. Blakeslee, B. C. Whitmore, P. Côté, D. Geisler, E. W. Peng,  and J. Bailin et. al, Astrophys. J. 797, 128 (2014).
60. E. F. Bell, D. H. McIntosh, N. Katz,  and M. D. Weinberg, Astrophys. J. Suppl. Ser. 149, 289 (2003).
61. The LIGO Scientific Collaboration and the Virgo Collaboration, ArXiv e-prints  (2016a), arXiv:1602.03840 [gr-qc] .
62. The LIGO Scientific Collaboration and the Virgo Collaboration, ArXiv e-prints  (2016b), arXiv:1602.03839 [gr-qc] .
63. M. Campanelli, C. Lousto, Y. Zlochower,  and D. Merritt, Astrophys. J. 659, L5 (2007).
64. S. Sigurdsson and L. Hernquist, Nature 364, 423 (1993).
65. E. A. Huerta, P. Kumar, S. T. McWilliams, R. OâShaughnessy,  and N. Yunes, Phys. Rev. D 90, 084016 (2014).
66. F. Antonini, S. Chatterjee, C. L. Rodriguez, M. Morscher, B. Pattabiraman, V. Kalogera,  and F. A. Rasio, Astrophys. J. 816, 65 (2016)arXiv:1509.05080 .
67. K. Belczynski, S. Repetto, D. Holz, R. O’Shaughnessy, T. Bulik, E. Berti, C. Fryer,  and M. Dominik, eprint arXiv:1510.04615  (2015).
68. R. F. Webbink, Astrophys. J. 277, 355 (1984).
69. J. S. W. Claeys, O. R. Pols, R. G. Izzard, J. Vink,  and F. W. M. Verbunt, Astron. Astrophys. 563, A83 (2014).
70. D. W. Hogg, ArXiv Astrophysics e-prints  (1999), astro-ph/9905116 .
71. K. Belczynski, A. Buonanno, M. Cantiello, C. L. Fryer, D. E. Holz, I. Mandel, M. C. Miller,  and M. Walczak, Astrophys. J. 789, 120 (2014)arXiv:1403.0677 [astro-ph.HE] .
72. L. Strolger, A. G. Riess, T. Dahlen, M. Livio, N. Panagia, P. Challis, J. L. Tonry, A. V. Filippenko, R. Chornock,  and H. Ferguson et. al, Astrophys. J. 613, 200 (2004).
73. Y. C. Pei, S. M. Fall,  and M. G. Hauser, Astrophys. J. 522, 604 (1999).
74. P. A. Young and C. L. Fryer, Astrophys. J. 670, 584 (2007).
75. T.-T. Yuan, L. J. Kewley,  and J. Richard, Astrophys. J. 763, 9 (2013).
76. B. Panter, R. Jimenez, A. F. Heavens,  and S. Charlot, Mon. Not. R. Astron. Soc. 391, 1117 (2008).
77. X.-J. Xu and X.-D. Li, Astrophys. J. 716, 114 (2010).
78. E. Pfahl, S. Rappaport,  and P. Podsiadlowski, Astrophys. J. 573, 283 (2002).
79. J. R. Hurley, S. J. Aarseth,  and M. M. Shara, Astrophys. J. 665, 707 (2007).
80. N. Ivanova, K. Belczynski, J. M. Fregeau,  and F. A. Rasio, Mon. Not. R. Astron. Soc. 358, 572 (2005).
81. D. S. Davis, H. B. Richer, J. Anderson, J. Brewer, J. Hurley, J. S. Kalirai, R. M. Rich,  and P. B. Stetson, Astron. J 135, 2155 (2008).
82. S. Chatterjee, S. Umbreit, J. M. Fregeau,  and F. A. Rasio, Mon. Not. R. Astron. Soc. 429, 2881 (2013).
83. S. Chatterjee, C. L. Rodriguez,  and F. A. Rasio, ArXiv e-prints  (2016), arXiv:1603.00884 .
84. C. J. Lada and E. A. Lada, Annu. Rev. Astron. Astrophys. 41, 57 (2003).
85. J. M. D. Kruijssen, Mon. Not. R. Astron. Soc. 454, 1658 (2015).
86. W. J. de Wit, L. Testi, F. Palla,  and H. Zinnecker, Astron. Astrophys. 437, 247 (2005).
87. B. M. Ziosi, M. Mapelli, M. Branchesi,  and G. Tormen, Mon. Not. R. Astron. Soc. 441, 3703 (2014).
88. S. F. Portegies Zwart, S. L. McMillan,  and M. Gieles, Annu. Rev. Astron. Astrophys. 48, 431 (2010).
89. C. L. Rodriguez, C.-J. Haster, S. Chatterjee, V. Kalogera,  and F. A. Rasio, In Prep  (2016).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters