Binary Black Hole Mergers from Globular Clusters: Implications for Advanced LIGO
Abstract
The predicted rate of binary black hole mergers from galactic fields can vary over several orders of magnitude and is extremely sensitive to the assumptions of stellar evolution. But in dense stellar environments such as globular clusters, binary black holes form by wellunderstood gravitational interactions. In this letter, we study the formation of black hole binaries in an extensive collection of realistic globular cluster models. By comparing these models to observed Milky Way and extragalactic globular clusters, we find that the mergers of dynamicallyformed binaries could be detected at a rate of per year, potentially dominating the binary black hole merger rate. We also find that a majority of clusterformed binaries are more massive than their fieldformed counterparts, suggesting that Advanced LIGO could identify certain binaries as originating from dense stellar environments.
I Introduction
By the end of this decade, the Advanced LIGO and Virgo detectors are expected to observe gravitational waves (GWs), ushering in a new postelectromagnetic era of astrophysics Aasi et al. (2015); Acernese et al. (2015). The most anticipated sources of observable GWs will be the signals generated by mergers of binaries with compact object components, such as binary neutron stars (NSs) or binary black holes (BHs). While coalescence rates of NSNS or BHNS systems can be constrained from observations, it is not currently possible to produce observationallymotivated rate predictions for BHBH mergers Belczynski et al. (2002). Typical detection rates of binary BH (BBH) mergers in galaxies can span several orders of magnitude from to with a fiducial value of Abadie et al. (2010); however, these estimates typically ignore the large numbers of BBHs that are formed through dynamical interactions in dense star clusters Morscher et al. (2015); Sadowski et al. (2008).
The dynamical formation of BBHs is a probabilistic process, requiring a very high stellar density. These conditions are believed to exist within the cores of globular clusters (GCs), very old systems of stars with radii of a few parsecs. Approximately 10 Myr after the formation of a GC, the most massive stars explode as supernovae, forming a population of single and binary BHs with individual masses from to Belczynski et al. (2006). The BHs, being more massive than the average star in the cluster, sink to the center of the GC via dynamical friction, until the majority of the BHs reside in the cluster core Spitzer (1969). After this “mass segregation” is complete, the core becomes sufficiently dense that threebody encounters can frequently occur Ivanova et al. (2005), producing BBHs at high rates. In effect, GCs are dynamical factories for BBHs: producing large numbers of binaries within their cores and ejecting them via energetic dynamical encounters.
In this letter, we use an extensive and diverse collection of GC models to study the population of BBHs that Advanced LIGO can detect from GCs. We explore how the observed parameters of a presentday GC correlate with the distribution of BBH inspirals it has produced over its lifetime. We then compare our models to the observed population of Milky Way GCs (MWGCs) and use recent measurements of the GC luminosity function to determine a mean number of BBH inspirals per GC. Finally, we combine these estimates with an updated estimate of the spatial density of GCs in the local universe (Appendix I) into a double integral over comoving volume and inspiral masses to compute the expected Advanced LIGO detection rate. We assume cosmological parameters of , , and , consistent with the latest combined Planck results Ade et al. (2015).
Ii Computing the Rate
We use a collection of 48 GC models generated by our Cluster Monte Carlo (CMC) code, an orbitaveraged Hénontype Monte Carlo code for collisional stellar dynamics Joshi et al. (2000). The models span a range of initial star numbers ( to ), initial virial radii (0.5 pc to 4 pc), and consider low stellar metallicities () and high stellar metallicities (). In addition, the code implements dynamical binary formation via threebody encounters, strong three and fourbody binary interactions, and realistic single and binary stellar evolution. See Appendix II for a complete description of our code and the models used.
Previous studies have explored the contribution of BBHs from GCs to the Advaned LIGO detection rate Zwart and McMillan (1999); O’Leary et al. (2006); Banerjee et al. (2010); Downing et al. (2011); Bae et al. (2014); Tanikawa (2013); however, the majority of these studies have relied on either approximate analytic arguments or simplified body models with particles and have assumed a single black hole mass of . The one exception is Downing et al. (2011), which used a Monte Carlo approach to model GCs of a realistic size (). However, their study only considered GCs of a single mass, and did not extrapolate that result to the observed distribution of GCs in the local universe. Ours is the first study to compare models with all the relevant physics over a range of masses to observed GCs. This comparison enhances our BBH merger rate by more than an order of magnitude over previous results.
We express the rate of detectable mergers per year, , as the following double integral over binary chirp mass () and redshift:
(1) 
This equation is similar to that found in Belczynski et al. (2013, 2014). The components of Eqn. 1 are as follows:

is the rate of merging BBHs from GCs with chirp mass at redshift .

is the fraction of sources with chirp mass at redshift that are detectable by a single Advanced LIGO detector.

is the comoving volume at a given redshift Hogg (1999).

is the time dilation between a clock measuring the merger rate at the source versus a clock on Earth.
This letter focuses on estimating the rate, , using our collection of GC models. We assume the rate can be expressed as the product of the mean number of inspirals per GC, the distribution of those sources in space, and the density of GCs in the local universe, i.e. . The spatial density of GCs in the local universe is taken to be , based on recent measurements of extragalactic GC systems Harris et al. (2013) and modern nearinfrared Schechter functions Kelvin et al. (2014). Note that this estimate, computed in Appendix I, is substantially lower than the previous estimate of from Zwart and McMillan (1999) that has been used in previous studies. We now estimate the values of and .
Iii mean number of mergers per cluster
To determine the mean number of BBH inspirals produced by a GC, we use the
collection of models to explore how the presentday observable parameters of
GCs relate to the number of BBHs it has produced over its lifetime. To quantify
the realism of a particlar model, we compare the total masses and concentrations of
our models to GCs observed in the Milky Way. The concentrations are measured
by considering the ratio of a cluster’s core radius to its halflight radius,
. This massconcentration space is similar to the “fundamental
plane” of GCs described in McLaughlin (2000), with in place of
the King concentration
Two trends emerge in our models. First, the total number of BBH inspirals over
12 Gyr is nearly linearly proportional to the final cluster mass. Second, the
number of inspirals is higher for more compact clusters (those with smaller
). Since the model coverage of the space is poorer than the
coverage of the mass, and since there are no detailed observations of
extragalactic GC concentrations, we elect to focus on the linear relationship
between a GC’s mass and the number of inspirals it has produced. We perform a
weighted linear regression for both lowmetallicity and highmetallicity
GCs (Fig. 1). The weights are created by generating a kernel
density estimate (KDE)
We compute the mean number of inspirals per GC by multiplying the linear relationships from Fig. 1 by the mass distribution of GCs. Recent work Harris et al. (2014) has suggested that the distribution of GC luminosities is universal and welldescribed by a lognormal distribution:
(2) 
with and . Assuming a masstolight ratio of 2 in solar units Harris et al. (2014); Bell et al. (2003), we convert this luminosity function to a mass function. We then compute the average number of inspirals per GC by integrating the linear relationship over the normalized GC mass function. The results for both metallicities and different highmass cutoffs are shown in Table 1.
Mass Cutoff  

Metallicity  
Low  430  967  1512 
High  830  1954  3103 
Iv Distribution of Inspirals
The numbers quoted in Table 1 provide us with the mean number of BBH inspirals from a GC over 12 Gyr. We could use this average rate to compute a detection rate for Advanced LIGO. However, it is qualitatively obvious that the mass distribution of BBH sources is not constant in time (Fig. 2).
Therefore, we must use the distribution of BBH inspiral events over time from GCs to compute the rate. We select inspirals randomly from each of our models, drawing more inspirals from models with higher weights according to the following scheme:
(3) 
where the weight, , of a model with mass and compactness at 12 Gyr is defined by the ratio of the MWGC KDE at divided by the KDE of the models themselves, evaluated at . The reason for these weights is as follows: we wish to draw more samples from models that are more likely to represent MWGCs, but because our collection of models is drawn from a different distribution (the initial conditions from Morscher et al. (2015)), we cannot simply draw inspirals at random from each model according to how well it represents real GCs. To do so would bias our samples with the distribution that results from our initial conditions. By dividing the probability of a model representing a MWGC by the probability density of our collection of models, our scheme naturally corrects for this. Models unlikely to represent MWGCs have small numerators and low weights. Models with no neighboring models that are likely to represent MWGCs have large numerators and small denominators, yielding high weights. Conversely, models with neighbors that are likely to represent MWGCs will have large numerators and large denominators, yielding smaller weights; however, as we will select some number of inspirals from each of these neighboring models, the cumulative effect is the same.
We show a sample distribution of the chirp masses versus redshift in Fig. 2. We distinguish between two different BBH formation channels: primordial and dynamical. We define primordial BBHs as those that are formed from the supernovae of two main sequence stars in a binary, and whose components were never bound to any other star before merger; conversely, we define dynamical binaries as those that are either formed from two isolated BHs via a threebody encounter, or formed from a higherorder dynamical encounter (a binarysingle or binarybinary interaction forming a new binary pair). Primordial binaries can still have their orbital parameters modified by dynamics (via a strong encounter with another BH or BBH), as long as the encounter leaves the primordial BBH intact. One immediately apparent feature is the bimodality between primordial and dynamical BBHs. Over all of our models, the highest chirp mass that is formed by pure binary stellar evolution is , as systems with larger progenitors are distrupted by the supernova kick. This implies that any source from our models with a detected chirp mass greater than could only have formed dynamically.
To compare this result to BBHs formed in the field, we generated two additional CMC models, each containing binaries and different metallicites ( and ). These models were computed without twobody relaxation, binary formation, or strong encounters, and only considered the physics of binary stellar evolution. In this dynamicsfree environment, the maximum chirp mass of any BBH was . Although this result depends on the metallicity and the physics of stellar evolution, it does suggest that GC dynamics forms BBHs consistently more massive than those in the field
V Detection Rate
We now compute the
expected rate of signals detectable by Advanced LIGO. To compute the fraction
of detectable sources, , we use gravitational waveforms
that cover the inspiral, merger, and ringdown phases of a compact binary merger
(known as IMRPhenomC waveforms Santamaría
et al. (2010)) and compute the
signaltonoise ratio (SNR) using the projected zerodetuning, highpower
configuration of Advanced LIGO
Our distribution of inspirals, , is generated by creating a KDE of the inspirals in Fig. 2. Since each “draw” of inspirals will produce a slightly different distribution, we compute 1000 times, and then take the mean of Eqn. 1 for those 1000 draws. We find that a single Advanced LIGO detector operating at design sensitivity will detect BBHs per year from GCs. Of those about will originate from lowmetallicity GCs and the rest from highmetallicity GCs, assuming 76% of GCs are low metallicity (consistent with the MWGC distribution). Approximately 80% of these sources will have chirp masses greater than , meaning that the majority of BBHs detectable by Advanced LIGO from GCs could only havie formed dynamically.
The majority of these BBH sources will be detected at low redshifts. For lowmetallicity clusters, the distribution of detectible sources in redshift peaks at , while for highmetallicity clusters the distribution peaks at . In both cases, 90% of detectable sources are located at .
To obtain a rough estimate of the uncertainty on this prediction, we perform a simple error analysis that considers the optimistic and conservative rates that would be obtained by varying our assumptions and selecting the estimates of certain quantities. For the conservative estimate, we assume that the GC mass function truncates at the mass of (), and that the spatial density of GCs is (the conservative estimate from Appendix I). We also recompute the rate using the uncertainty from the regression in Fig. 1 and the lower standard deviation of our 1000 draws of . This yields a conservative estimate of BBH inspirals per year. Conversely, if we assume the most optimistic truncation mass for GC mass function (), the most optimistic GC spatial density (, the optimistic estimate from Appendix I, similar to the value used in previous studies), and the uncertainties on the linear regression and , we find an optimistic rate of BBH inspirals per year. This range is primarily influenced by the uncertainty in the GC spatial density and the truncation mass of the GC mass function.
Vi Conclusion
In this letter, we compared new GC models computed with our CMC code to the observed distributions of Milky Way and extragalactic GCs to predict the expected rate of BBH inspirals from realistic GCs. We determined a linear relationship between the presentday mass of a GC and the number of BBH inspirals produced by that cluster. By combining this with the universal GC luminosity function and a new estimate for the spatial density of GCs, we were able to predict the mean density of BBH inspirals from GCs in the local universe. Then by weighting our models according to their similarity to the observed distribution of MWGCs, we created a distribution of inspiral sources in chirp mass and redshift. Finally, by combining this with the anticipated sensitivity of Advanced LIGO, we estimated a detection rate of BBH inspiral events per year from GC sources. With highly conservative assumptions, this rate drops to events per year, while highly optimistic assumptions pushes the rate as high as events per year.
We also found that no BBHs with chirp masses above were formed from a primordial binary. In other words, every inspiral with in our models was formed by dynamical processes alone. This could, in theory, provide an easy way to discriminate between binaries that were formed dynamically versus those formed by binary stellar evolution; however, this result is highly dependent on the physics of supernova kicks and the fraction of ejected supernova material which falls back onto the newly formed BH, both of which remain poorly constrained. In addition, recent work has suggested that the mass distribution of chirp masses for BBHs produced by stellar evolution can reach as high as , depending on the physics of the common envelope Belczynski et al. (2010). As such, this result should be treated as a proofofprinciple, and not a concrete physical claim. Investigations to better understand the relationship between this formation cutoff, the distribution of supernovae kicks, and the fallback fraction, are currently underway.
Finally, the number of BHs formed is entirely dependent on the choice of the initial mass function. Although our choice of IMF is typical for studies of this type, a variation of in the slope of the highmass end of the IMF can produce significant differences in the number of BBHs. Investigations to quantify this effect are also underway.
Acknowledgements.
We thank Ilya Mandel for his assistance in computing the fraction of detectable sources, as well as useful discussions. We also thank Ben Farr, ClaudeAndre FaucherGiguere, William Harris, Kyle Willett, Michael Schmitt, and Tyson Littenberg for useful discussions. CR was supported by an NSF GRFP Fellowship, award DGE0824162. This work was supported by NSF Grant AST1312945 and NASA Grant NNX14AP92G. CJH acknowledges support from an RAS grant.Appendix A Number of Globular Clusters per Mpc
In order to estimate the rate of inspirals from an average GC per , we must compute the average spatial density of GCs in the local universe. This can be accomplished by considering the mean number of GCs per galaxy at a given luminosity, multiplied by the spatial density of galaxies at that luminosity, then summing over all luminosities, as illustrated in the following equation:
(4) 
The number of GCs per galaxy per luminosity can be determined by use of the Harris Globular Cluster System catalog Harris et al. (2013). The catalog provides a list of 422 galaxies, their morphological type, visual and Kband magnitudes (where available), and the estimated total number of GCs. In Figure 3, we plot the 346 galaxies for which Kband photometry is available in the catalog against the estimated number of GCs. For each collection of galaxy morphologies, we perform a Gaussian Process regression with the George package, described in Ambikasaran et al. (2014). The Gaussian processes are generated using a squaredexponential kernel combined with a white noise kernel, and then fit to the log of the number of GCs per galaxy. The kernel hyperparameters are selected by maximizing the marginalized loglikelihood of the Gaussian Process. See Rasmussen and Williams (2005) for a detailed description of regression with Gaussian Processes. The mean and standard deviation of the resulting Gaussian Processes are also shown in Fig. 3. Note that the catalog does not include lowluminosity dwarf and irregular galaxies for which . This suggests that our fitted function can be systematically overestimating the number of GCs from dwarf earlytype galaxies with Peng et al. (2008).
All  

Elliptical  
Lenticular  
Spirals (SaSd)  
Irregular 
The mean in Figure 3 gives us the number of GCs per galaxy per Kband luminosity. To compute a spatial density, we must integrate the mean over the density of galaxies per Kband luminosity per . For this, we can use the well known Schechter function Schechter (1976), which describes the spatial density of galaxies per absolute magnitude interval . We use the recentlycomputed Kband Schechter functions fits from the Galaxy and Mass Assembly (GAMA) Survey Kelvin et al. (2014), which contains individual fits to individual Schechter functions for each galaxy morphology, and a single fit to a double Schechter function for the combined sample of all galaxies. We use both to determine our overall density of GCs, as well as the density contributed by galaxies of each type.
Finally, we multiply each Schechter function from Kelvin et al. (2014) by the estimated number of GCs per galaxy from Figure 3, and integrate over all Kband magnitudes (). This results in a GC density estimate of , which we employ in our rate calculation. For completeness, we also consider different cutoffs for our magnitude integral, and report the contribution to the spatial densities from each galaxy morphology, in Table 2.
Since the Schechter Function diverges at low luminosities, and since our fit systematically overestimates the number of GCs for lowluminosity galaxies, we must pick a reasonable limit at which to truncate our integral. We use a lower limit of , although for completeness we also consider lower () and higher () cutoffs in Table 2.
In addition to comparisons with observations, we can also compute the density of GCs using cosmological simulations. The publication of the GC Systems catalog in Harris et al. (2013) noted a correlation between the dynamical mass of a galaxy and the size of its GC population. This relationship was expanded upon in Harris et al. (2015), which measured a very strong correlation between the mass of the GC population and the galaxy halo mass. This relation takes the following form:
(5) 
where is 7.706(7.405)(7.157), is 1.03(0.96)(1.21), and is 12.3(12.2)(12.2) for all(lowmetallicity)(highmetallicity) GCs. Unlike the relationship, the relationship between halo mass and does not strongly depend on the galaxy morphology.
In order to convert this to a spatial density of GCs, we can multiply this relationship by the dark matter halo mass function, as determined by recent cosmological simulations. We use the functional fit to from Tinker et al. (2010), as calculated at redshift by the HMFcalc website Murray et al. (2013). We then compute the integral
(6) 
where we use , the mean of the GC mass function from Harris et al. (2014) used in the main text, to convert from the mass of a GC system to the number of GCs. This yields a spatial density of , or , assuming the value of used throughout the text. We can also use the similar values for low and highmetallicitiy GCs quoted below Eqn. 5. This yields estimates of and , respectively.
Appendix B Models
This letter considers the BBH inspirals from 48 separate GC models generated with our orbitaveraged Hénontype Monte Carlo code, CMC. The majority of these models were first developed in Morscher et al. (2015). The full details of CMC can be found in previous papers Joshi et al. (2000, 2001); Pattabiraman et al. (2013), but the features most relevant to this letter are as follows:

threebody binary formation, which we implement with a probabilistic analytic prescription Morscher et al. (2015),

strong singlebinary and binarybinary stellar encounters, implemented with the small integrator Fewbody Fregeau et al. (2004), and

single and binary stellar evolution with the SSE and BSE packages Hurley et al. (2000, 2002). Note that our implementation includes several improvements Chatterjee et al. (2010), including the stellar remnant prescription and BH kick physics from Belczynski et al. (2002). For BBHs which merge within the cluster, the GW timescale is calculated by BSE. For ejected binaries, the inspiral time is found by integrating the orbitaveraged Peters equations Peters (1964) using the masses, separation, and eccentricity of the binary at the time of ejection.
The models begin with , , and number of particles, and are evolved to an age of 12 Gyr each. We do not include GC models which dissipate before 12 Gyr, as we have no way to compare these models to observations. However, since these models all begin with low numbers of particles and produce low numbers of BBHs, the effect on the rate estimate will be minimal. For a complete list of GC initial conditions considered, see Table 3.
We also explore the space of initial cluster sizes and concentrations, with initial virial radii of 0.5, 1, 1.5, 2, and 4 pc and initial King concentrations () of 2, 5, 7, and 11. We qualitatively find that the initial King concentration does not have a strong influence on the observable GC properties at 12 Gyr. Each of our models starts with 10% of the objects in primordial binaries, and stellar masses chosen from a universal initial mass function (IMF) Kroupa (2001).
We also explore metallicities of , , and , which are placed at galactocentric distance of 2, 8 and 20 kpc respectively. This is due to the observed correlation between GC metallicity and galactocentric distance Djorgovski and Meylan (1994). Although we explore three distinct metallicities, we separate our models into “low metallicitiy” (those for which [Fe/H] , i.e. and ), and “high metallicity” () GCs. This is chosen to simplify the comparison between our models and the observations of GCs, which show a strong bimodality in metallicity Harris (2010). We also assume that the fraction of lowmetallicity GCs is 0.76, since that is the fraction of MWGCs for which [Fe/H] .
Initial Conditions  Properties (12 Gyr)  

N ()  (pc)  Metallicity (z)  Mass ()  BH Retained  
0.5  0.001  0.55  0.17  0.05  16  
1.0  0.001  0.20  0.34  0.03  11  
1.0  0.0005  0.68  0.41  0.21  8  
1.5  0.001  0.54  0.42  0.11  8  
1.5  0.0005  0.62  0.45  0.18  7  
2.0  0.0005  0.63  0.41  0.20  5  
2.0  0.001  0.55  0.77  0.16  1  
2.0  0.0005  0.64  0.53  0.20  8  
2.0  0.001  0.56  1.12  0.18  6  
2.0  0.0005  0.65  0.45  0.23  5  
2.0  0.001  0.56  0.47  0.20  5  
2.0  0.001  0.50  0.99  0.22  3  
4.0  0.001  0.68  0.58  0.37  2  
0.5  0.0005  2.58  0.27  0.10  145  
1.0  0.001  2.17  0.46  0.18  91  
1.0  0.0005  2.78  0.48  0.24  87  
1.0  0.005  2.00  0.45  0.14  88  
1.5  0.0005  2.73  0.47  0.32  81  
1.5  0.001  2.61  0.54  0.30  67  
1.5  0.005  2.02  0.64  0.23  59  
2.0  0.0005  2.76  0.52  0.43  48  
2.0  0.001  2.62  0.58  0.40  53  
2.0  0.005  2.04  0.63  0.30  45  
2.0  0.0005  2.79  0.59  0.40  49  
2.0  0.005  2.00  0.52  0.32  45  
2.0  0.001  2.66  0.41  0.36  59  
2.0  0.0005  2.81  0.77  0.44  66  
2.0  0.001  2.70  0.5  0.40  57  
2.0  0.005  1.92  0.77  0.32  52  
2.0  0.001  2.60  0.67  0.45  38  
4.0  0.001  2.91  0.61  0.58  27  
1.0  0.0005  5.44  0.5  0.28  324  
1.0  0.001  4.69  0.65  0.27  269  
1.0  0.005  4.49  0.52  0.24  270  
1.5  0.0005  5.59  0.72  0.42  204  
1.5  0.001  5.39  0.71  0.41  175  
1.5  0.005  4.75  0.59  0.33  203  
2.0  0.0005  5.68  0.57  0.56  159  
2.0  0.001  5.50  0.56  0.49  154  
2.0  0.005  4.84  0.75  0.43  162  
2.0  0.0005  5.76  0.6  0.52  172  
2.0  0.001  5.56  0.66  0.52  139  
2.0  0.005  4.97  0.57  0.41  175  
2.0  0.0005  5.76  0.56  0.53  167  
2.0  0.001  5.58  0.67  0.52  140  
2.0  0.005  4.96  0.63  0.41  153  
2.0  0.001  5.47  0.61  0.52  121  
4.0  0.001  5.94  0.79  0.67  99 
Footnotes
 preprint: APS/123QED
 The King concentration is defined as the logarithm of the core radius over the tidal radius. We use the simpler , as the tidal radius can be difficult to determine observationally.
 We select the bandwidths using maximum likelihood crossvalidation, since for multidimensional distributions, the optimal bandwidth is not necessarily equal along each dimension, particularly when using mixed units. This is implemented in the StatsModels Python package.
 Both the noise curve and technical reports describing it can be found under LIGO Document T0900288v3
References
 J. Aasi, B. P. Abbott, R. Abbott, T. Abbott, M. R. Abernathy, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, et al., Class. Quant. Grav. 32, 074001 (2015).
 F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, G. Ballardin, et al., Class. Quant. Grav. 32, 024001 (2015).
 K. Belczynski, V. Kalogera, and T. Bulik, Astrophys. J. 572, 407 (2002).
 J. Abadie, B. P. Abbott, R. Abbott, M. Abernathy, T. Accadia, F. Acernese, C. Adams, R. Adhikari, P. Ajith, B. Allen, et al., Class. Quant. Grav. 27, 173001 (2010), eprint 1003.2480.
 M. Morscher, B. Pattabiraman, C. Rodriguez, F. A. Rasio, and S. Umbreit, Astrophys. J. 800, 9 (2015), eprint 1409.0866.
 A. Sadowski, K. Belczynski, T. Bulik, N. Ivanova, F. A. Rasio, and R. O’Shaughnessy, Astrophys. J. 676, 1162 (2008), eprint 0710.0878.
 K. Belczynski, A. Sadowski, F. A. Rasio, and T. Bulik, Astrophys. J. 650, 303 (2006).
 L. Spitzer, Jr., Astrophys. Lett. 158, L139 (1969).
 N. Ivanova, K. Belczynski, J. M. Fregeau, and F. A. Rasio, Mon. Not. R. Astron. Soc. 358, 572 (2005), eprint 0501131.
 P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, and et al., ArXiv eprints (2015), eprint 1502.01589.
 K. J. Joshi, F. A. Rasio, S. P. Zwart, and S. Portegies Zwart, Astrophys. J. 540, 969 (2000), eprint 9909115.
 S. P. Zwart and S. McMillan, Astrophys. Lett. p. 12 (1999), eprint 9910061.
 R. M. O’Leary, F. A. Rasio, J. M. Fregeau, N. Ivanova, and R. O’Shaughnessy, Astrophys. J. 637, 937 (2006).
 S. Banerjee, H. Baumgardt, and P. Kroupa, Mon. Not. R. Astron. Soc. 402, 371 (2010).
 J. M. B. Downing, M. J. Benacquista, M. Giersz, and R. Spurzem, Mon. Not. R. Astron. Soc. 416, 133 (2011), eprint 1008.5060.
 Y.B. Bae, C. Kim, and H. M. Lee, Mon. Not. R. Astron. Soc. 440, 2714 (2014).
 A. Tanikawa, Mon. Not. R. Astron. Soc. 435, 1358 (2013).
 K. Belczynski, T. Bulik, I. Mandel, B. S. Sathyaprakash, A. A. Zdziarski, and J. Mikołajewska, Astrophys. J. 764, 96 (2013), eprint 1209.2658.
 K. Belczynski, A. Buonanno, M. Cantiello, C. L. Fryer, D. E. Holz, I. Mandel, M. C. Miller, and M. Walczak (2014), eprint 1403.0677.
 D. W. Hogg, ArXiv Astrophysics eprints (1999), eprint astroph/9905116.
 W. E. Harris, G. L. H. Harris, and M. Alessi, Astrophys. J. 772, 82 (2013), eprint 1306.2247.
 L. S. Kelvin, S. P. Driver, A. S. G. Robotham, A. W. Graham, S. Phillipps, N. K. Agius, M. Alpaslan, I. Baldry, S. P. Bamford, J. BlandHawthorn, et al., Mon. Not. R. Astron. Soc. 439, 1245 (2014), eprint 1401.1817.
 D. E. McLaughlin, Astrophys. J. 539, 618 (2000).
 The King concentration is defined as the logarithm of the core radius over the tidal radius. We use the simpler , as the tidal radius can be difficult to determine observationally.
 We select the bandwidths using maximum likelihood crossvalidation, since for multidimensional distributions, the optimal bandwidth is not necessarily equal along each dimension, particularly when using mixed units. This is implemented in the StatsModels Python package.
 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, J. Bailin, et al., Astrophys. J. 797, 128 (2014).
 E. F. Bell, D. H. McIntosh, N. Katz, and M. D. Weinberg, p. 23 (2003), eprint 0302543.
 J. Strader, A. C. Seth, D. A. Forbes, G. Fabbiano, A. J. Romanowsky, J. P. Brodie, C. Conroy, N. Caldwell, V. Pota, C. Usher, et al., Astrophys. J. 775, L6 (2013).
 L. Santamaría, F. Ohme, P. Ajith, B. Brügmann, N. Dorband, M. Hannam, S. Husa, P. Mösta, D. Pollney, C. Reisswig, et al., Phys. Rev. D 82, 064016 (2010).
 Both the noise curve and technical reports describing it can be found under LIGO Document T0900288v3.
 K. Belczynski, M. Dominik, T. Bulik, R. O’Shaughnessy, C. Fryer, and D. E. Holz, Astrophys. J. 715, L138 (2010).
 S. Ambikasaran, D. ForemanMackey, L. Greengard, D. W. Hogg, and M. O’Neil, ArXiv eprints (2014), eprint 1403.6015.
 C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) (The MIT Press, 2005), ISBN 026218253X.
 E. W. Peng, A. Jordan, P. Cote, M. Takamiya, M. J. West, J. P. Blakeslee, C.W. Chen, L. Ferrarese, S. Mei, J. L. Tonry, et al., Astrophys. J. p. 27 (2008), eprint 0803.0330.
 P. Schechter, Astrophys. J. 203, 297 (1976).
 W. Harris, G. Harris, and M. Hudson, ArXiv eprints (2015), eprint 1504.03199.
 J. L. Tinker, B. E. Robertson, A. V. Kravtsov, A. Klypin, M. S. Warren, G. Yepes, and S. Gottlöber, Astrophys. J. 724, 878 (2010).
 S. Murray, C. Power, and A. Robotham, Astronomy and Computing 34, 23 (2013).
 K. J. Joshi, C. P. Nave, and F. A. Rasio, Astrophys. J. 550, 691 (2001).
 B. Pattabiraman, S. Umbreit, W.k. Liao, A. Choudhary, V. Kalogera, G. Memik, and F. A. Rasio, Astrophys. J. Suppl. Ser. 204, 15 (2013), eprint 1206.5878.
 J. M. Fregeau, P. Cheung, S. F. Portegies Zwart, and F. A. Rasio, Mon. Not. R. Astron. Soc. 352, 1 (2004).
 J. R. Hurley, O. R. Pols, and C. A. Tout, Mon. Not. R. Astron. Soc. 315, 543 (2000).
 J. R. Hurley, C. A. Tout, and O. R. Pols, Mon. Not. R. Astron. Soc. 329, 897 (2002).
 S. Chatterjee, J. M. Fregeau, S. Umbreit, and F. A. Rasio, Astrophys. J. 719, 915 (2010).
 P. Peters, Phys. Rev. 136, B1224 (1964).
 P. Kroupa, Mon. Not. R. Astron. Soc. 322, 231 (2001).
 S. Djorgovski and G. Meylan, Astron. J. 108, 1292 (1994).
 W. E. Harris, Royal Society of London Philosophical Transactions Series A 368, 889 (2010), eprint 0911.0798.