# Post-Newtonian Dynamics in Dense Star Clusters:

Highly-Eccentric, Highly-Spinning, and Repeated Binary Black Hole Mergers

###### Abstract

We present models of realistic globular clusters with post-Newtonian dynamics for black holes. By modeling the relativistic accelerations and gravitational-wave emission in isolated binaries and during three- and four-body encounters, we find that nearly half of all binary black hole mergers occur inside the cluster, with about 10% of those mergers entering the LIGO/Virgo band with eccentricities greater than 0.1. In-cluster mergers lead to the birth of a second generation of black holes with larger masses and high spins, which, depending on the black hole natal spins, can sometimes be retained in the cluster and merge again. As a result, globular clusters can produce merging binaries with detectable spins regardless of the birth spins of black holes formed from massive stars. These second-generation black holes would also populate any upper mass gap created by pair-instability supernovae.

^{†}

^{†}preprint: APS/123-QED

## I Introduction

With the recent detections of five binary black hole (BBH) mergers, and one binary neutron star merger, the era of gravitational wave (GW) astrophysics has arrived at last Abbott et al. (2017a, b, c, 2016a, 2016b). Despite significant theoretical work, the origins of these systems, particularly the heavier BBHs, remain an open question. Both stellar evolution in isolated massive binaries (e.g. Belczynski et al., 2010a; Marchant et al., 2016; Podsiadlowski et al., 2003; Mandel and De Mink, 2016; De Mink and Mandel, 2016) and dynamical formation in dense star clusters (e.g. Sigurdsson and Hernquist, 1993; Portegies Zwart and Mcmillan, 2000; OâLeary et al., 2006; Miller and Lauburg, 2009; Downing et al., 2010; Banerjee et al., 2010; Downing et al., 2011; Bae et al., 2014; Ziosi et al., 2014; Rodriguez et al., 2015; Rodriguez et al., 2016a; Banerjee, 2017) have been shown to produce merging BBHs similar to GW150914 Belczynski et al. (2016a); Rodriguez et al. (2016b). Understanding which formation pathways are at play will be critical for the interpretation of GW data. While many signatures of dynamical assembly have been proposed, such as highly-eccentric mergers occurring in strong chaotic encounters Samsing et al. (2014) or anti-alignment of the BH spins with the orbit Rodriguez et al. (2016c), none of the BBH mergers detected so far by LIGO/Virgo have displayed any of those signatures clearly (see Amaro-Seoane and Chen, 2016).

What has been displayed clearly in each BBH merger is the birth of a new rapidly-spinning BH with a mass (almost) equal to the sum of its progenitor masses. Many of these new BHs, particularly the remnants of GW150914, GW170104, and GW170814, are significantly more massive than what is thought to form during the collapse of a single star, where the pair-instability mechanism limits the remnant BH mass to Woosley (2016). Were one of these mergers to occur in a dense star cluster, however, the merger product could easily exchange into another BBH and merge again. Because of the distinct BH masses and spins in such second-generation (2G) mergers, it has been suggested that such a population could be easily identifiable with future LIGO/Virgo detections Fishbach et al. (2017); Gerosa and Berti (2017).

In this letter, we present the first models of realistic globular clusters (GCs) with fully post-Newtonian (pN) stellar dynamics. While relativistic -body dynamics has been studied previously for highly idealized systems (e.g., Lee, 1993, 1992; Shapiro and Teukolsky, 1985a, b; Rasio et al., 1989; Quinlan and Shapiro, 1987, 1989, 1990; Kupi et al., 2006; Brem et al., 2013) or open clusters (e.g., Banerjee, 2017), we show here for the first time using self-consistent dynamical models of massive GCs that pN effects play a key role in assembling dynamically the merging BBHs detectable by LIGO/Virgo. In our new pN models, we observe that roughly half of all BBH mergers occur inside clusters, with a significant fraction of those () merging with eccentricities greater than 0.1 following GW captures. In-cluster mergers produce a second generation of BHs that, if not ejected from the cluster through GW recoil, will dynamically exchange into new binaries only to merge again. These 2G BBH mergers have components with large spins and masses significantly beyond what is possible from the collapse of a single star; they may be quite common, with as many as 20% of BBH mergers from our models having components formed in a previous merger.

Throughout this paper, we assume a flat cosmology with and Planck Collaboration et al. (2015).

## Ii Post-Newtonian Dynamics

We have computed the new GC models presented here using the Cluster Monte Carlo (CMC) code. CMC has been developed over many years Joshi et al. (2000); Pattabiraman et al. (2013), and includes all the necessary physics for the long-term evolution of GCs, including two-body relaxation Hénon (1971); Henon (1971), single and binary stellar evolution Hurley et al. (2000, 2002); Chatterjee et al. (2010), galactic tides, three-body binary formation Morscher et al. (2013), and three- and four-body gravitational encounters via the fewbody package Fregeau et al. (2004); Fregeau and Rasio (2007). We have shown in Rodriguez et al. (2016d) that CMC can reproduce with a high degree of fidelity both the global cluster properties and BBH distributions computed with state-of-the-art direct -body simulations Wang et al. (2016), while at the same time being at least two orders of magnitude faster (essential for the sort of extensive parameter-space study presented here). Furthermore, CMC has been upgraded Rodriguez et al. (2016a) to employ the most recent prescriptions for stellar-wind-driven mass loss Vink et al. (2001); Belczynski et al. (2010b) and compact-object formation Fryer et al. (2012), allowing us to compare our results directly to those of population synthesis studies for isolated binaries (e.g. Belczynski et al., 2016a).

To incorporate pN effects into CMC, we make the following modifications. We account for relativistic effects during three- and four-body encounters by adopting a modified version of the fewbody code with pN accelerations up to and including the 2.5pN order. This code has been described in detail in Antognini et al. (2014); Amaro-Seoane and Chen (2016) and has been shown to conserve energy to 2pN order and to reproduce the inspiral times for compact binaries Peters (1964). For BBHs which merge during an encounter, we perform a standard sticky-sphere merger, using detailed, spin-dependent fitting formulas from analytic and numerical relativity calculations Barausse et al. (2012); Kesden (2008); Lousto and Zlochower (2014); Berti et al. (2007); Campanelli et al. (2007); Tichy and Marronetti (2008); Barausse and Rezzolla (2009); Kesden (2008); Buonanno et al. (2008); Lousto and Zlochower (2014); Tichy and Marronetti (2008); Rezzolla et al. (2008); González et al. (2007); Lousto and Zlochower (2008); Lousto et al. (2012); Lousto and Zlochower (2013). The new masses, spins, and recoil kicks are applied immediately during any merger, allowing us to model the retention of BHs by the cluster self-consistently. See Supplemental Materials A for details, which includes references Wen (2003); Blanchet (2006); Sopuerta et al. (2007); Yunes and Berti (2008); Gerosa and Kesden (2016); Jiménez-Forteza et al. (2017) . We initially assume all BHs from stellar collapse have no spins at birth (, where is the dimensionless Kerr spin parameter), though we relax this assumption in Section V. For BBHs which do not merge during a fewbody encounter, we directly integrate the orbit-averaged Peters equations Peters (1964) for the change in semi-major axis and eccentricity due to GW emission. This represents a departure from our previous work where we relied on the binary stellar evolution module (BSE, Hurley et al., 2002). By default, BSE only applies GW energy loss to binaries with . This assumption leads BSE to significantly underestimate the number of GW-driven mergers for binaries in a typical cluster, which can be highly eccentric and very massive. When accounting for GW energy loss in eccentric binaries, the number of in-cluster mergers becomes comparable to the number of merging binaries that are ejected from the cluster. This is a significant improvement over previous results in the literature (e.g., Rodriguez et al., 2015; Rodriguez et al., 2016a, b; Downing et al., 2010; Askar et al., 2016), where ejected BBHs dominated the merger rate in the local universe. Semi-analytic approaches to cluster dynamics (e.g., Antonini and Rasio, 2016; Giesler et al., 2017) have reported significantly higher fractions of in-cluster mergers, similar to those presented here, and have noted the possibility of multiple mergers in galactic nuclei Antonini and Rasio (2016).

We generate 24 GC models covering a range of masses, metallicities, galactocentric distances, and virial radii, similar to those observed in the Milky Way and beyond. These initial conditions are identical to those from Rodriguez et al. (2016a), allowing us to explicitly compare our pN results to those in the literature. Our physics for single and binary stellar evolution is nearly identical to Rodriguez et al. (2016a). We have added a prescription for stellar mass loss via pulsational-pair-instability supernovae and stellar destruction via pair-instability supernovae. This physics, powered by the rapid production of electron-positron pairs in the stellar core Woosley (2016), places a well-understood upper limit on the masses of BHs that can from from the collapse of a single star. We take the limit from Belczynski et al. (2016b) of which is reduced to via neutrino emission. See Appendix B for details. This limit is in tentative agreement with the BH mass distribution measured by LIGO/Virgo Fishbach and Holz (2017). In our simulations, no BH can be born with a mass above unless the BH or its stellar progenitor has undergone a dynamical merger or mass transfer. Finally, unlike previous studies Rodriguez et al. (2015); Rodriguez et al. (2016a, b), we have not weighted our models according to the distribution of observed GCs. We will explore more realistic sets of models in future work focusing specifically on LIGO/Virgo detection rates. In practice a more realistic weighting should make little difference, as our previously adopted weighting scheme primarily selected BBHs from the most massive clusters, which also contribute the majority of sources in our current grid.

## Iii In-Cluster Mergers

With the addition of the pN physics, we see a significant increase in the number of in-cluster mergers. Whereas before the number of in-cluster mergers was a minor correction to the BBH mergers in the local universe (0.06% of mergers at , see Rodriguez et al., 2016a), we now find that nearly half of mergers now occur inside the cluster. For the 24 models considered here we find a total of 2819 mergers, 55% of which occur in the cluster. At low redshifts (), this number decreases to 45%, as the primordial binaries which merged at early times after a common-envelope phase have merged many Gyr ago. Compared to similar models without pN physics Rodriguez et al. (2016a), the number of ejected BBH mergers at decreases by (496 versus 410). However, the number of in-cluster mergers has jumped significantly, from one to 338. This increases the total number of mergers (in-cluster and ejected) by .

This increase in the number of BBH mergers occurring in the cluster primarily arises from properly accounting for GW emission for binaries regardless of their semi-major axis. For example, a typical BBH is ejected from a GC with (roughly 10 times greater than the cutoff in BSE) after undergoing dynamical encounters Rodriguez et al. (2016b). During a typical encounter, the BBH semi-major axis will characteristically shrink while the orbital eccentricity randomly drawn from the thermal distribution, Heggie (1975). These “hardening” encounters continue, shrinking the binary’s semi-major axis until either the BBH is ejected from the cluster by the third body or until GWs drive the binary to merger. The timescale for each BBH to merge can be roughly approximated by Peters (1964):

(1) |

As (1) makes clear, a large eccentricity can significantly decrease the merger timescale. For (roughly %10 post-encounter binaries) will decrease by more than , leading the BBH to promptly merge in the cluster. On the other hand, for BBHs that never reach a high eccentricity, these encounters will continue to harden the binary until it is ejected from the cluster (where its eccentricity at ejection is set by a single draw from the thermal distribution). Because the dependence in (1) preferentially selects in-cluster mergers from a super-thermal distribution, we expect these mergers to have larger eccentricities than their ejected counterparts by the time they reach the LIGO/Virgo band.

In Figure 1, we show the eccentricity distribution of merging binaries as they enter the LIGO/Virgo band (which we define as a circular GW frequency of 10Hz). We see the expected separation in eccentricity between BBHs which merge in the cluster and those that merge after being ejected from the cluster. For the in-cluster mergers, we also find a clear bimodality, with the lower peak corresponding to isolated binaries that merge after a dynamical encounter and the higher peak () corresponding to sources which merge during the encounter via GW capture. Although previous work Samsing et al. (2014); Amaro-Seoane and Chen (2016); Samsing and Ramirez-Ruiz (2017) has shown through scattering experiments that such mergers are to be expected at the 1% level, this is the first work to show that these mergers occur in realistic GC environments. From our combined 24 models, we find that about 10% of the in-cluster mergers ( of all mergers) at occur during these GW captures, in good agreement with analytic work Samsing (2017).

## Iv Mergers over Cosmic Time

In Figure 2, we show the mergers of BBHs as a function of cosmological redshift. What is immediately striking is that the mass distributions for in-cluster and ejected binaries are significantly different at low redshifts. This arises from the delay times between formation and mergers for ejected BBHs. When a BBH is ejected from the cluster, it may still take several Gyr to merge in the field (see e.g., Benacquista and Downing, 2013, and references therein). Even for the most massive clusters, the median inspiral time for ejected binaries is (e.g., Rodriguez et al., 2016a, Figure 1). In effect, the ejected BBHs which merge today drew their components from the initial distribution of BH masses in the cluster, where the masses varied from to . On the other hand, the in-cluster mergers have effectively no delay time, and their components are drawn from the present-day distribution of BH masses in the cluster. Because old GCs have ejected their most-massive BHs many Gyrs ago Morscher et al. (2015), the BBHs merging in the cluster today are typically lower-mass than those that were ejected many Gyr ago.

Another interesting feature of Figure 2 is the presence of BBH mergers in the upper-mass gap, beyond the mass limit imposed by pair-instability supernovae. The increased number of in-cluster mergers allows the GCs to produce significant numbers of 2G BBH mergers, some of which will have components above the maximum mass for BHs born from a single stellar collapse. As these systems can only be produced through multiple mergers, they will immediately be identifiable as having arisen from a dynamical environment. The rate of such mergers is small, but LIGO/Virgo is more sensitive to mergers with more massive components (the detection horizon scales with the mass of the more massive component as , Fishbach and Holz, 2017). At the expected sensitivity for Advanced LIGO’s third observing run Abbott et al. (2016c), a BBH with component masses of could be detected out to , encompassing a comoving volume of space three times larger than was observed during LIGO’s second science run Chen et al. (2017).

## V Black Hole Spin and Recoil Kicks

As a conservative assumption, we have assumed that all BHs in the cluster are born with no intrinsic spin. This is consistent with all but one (GW151226, Abbott et al., 2016b) of the BBHs detected by LIGO/Virgo so far. However, the presence of high BH spins, suggested by observations of BH X-ray binaries (see Miller and Miller, 2015, for a review), can radically change the results presented here: depending on the spin magnitudes and orientations, merging BBHs can get kicks as high as 5000 km/s (e.g., Campanelli et al., 2007; Lousto and Zlochower, 2011; Lousto et al., 2012), significantly larger than the escape speed of a typical GC. As a result, the 2G mergers shown in the left-hand panel of Figure 2 would not have formed if BHs are born with large spins, since their components would not have been retained in the cluster Miller and Lauburg (2009).

We can estimate how the numbers in Figure 2 would have changed under different assumptions for BH birth spins. For each repeated merger, we calculate the probability that each of the components would have been retained in the cluster given different birth spins. This is done by computing the recoil kicks over 1000 realizations of the spin orientations at merger. The probability of retaining each progenitor is simply the fraction of mergers for which the recoil speed is smaller than the cluster escape speed where the merger occurred. For each 2G BBH merger, we take the product of the retention probabilities for each component as the probability of that 2G merger occurring. We show the retention of these BBHs in the right panel of Figure 2 by weighting each 2G BBH merger by its retention probability. As expected, the number of 2G BBH mergers decreases as the birth spins of the BHs are increased. When , we find that 20% of mergers at are 2G mergers. As the spins are increased, this number decreases, and once , we observe 2G mergers, compared to the 672 first-generation mergers which occur at .

These assumption have significant implications for the measurable spins of BBH mergers. As shown by numerical relativity Berti and Volonteri (2008); Tichy and Marronetti (2008); Lousto et al. (2010) and idealized pN -body simulations with spins Brem et al. (2013), repeated mergers of BBHs in clusters with near-equal masses will tend to produce BHs with , (assuming the initial spins are isotropically distributed Kesden et al. (2010)). But what LIGO/Virgo is most sensitive to is not the spin magnitudes of the BBHs components Ajith et al. (2011); Pürrer et al. (2016), but the effective spin of the BBH, defined as the mass-weighted projection of the two spins onto the orbital angular momentum:

(2) |

where is the direction of the orbital angular momentum and are the dimensionless-spin vectors for the BHs. For dynamically-formed binaries, the isotropic distribution of the orbit and spin vectors means that (2) will be peaked at with symmetric tails whose extent depends on the BH spin magnitudes. We show the distributions of in Figure 3. When the initial BH spins are low, the 2G systems are the only BBHs which merge with observably large spins. The fraction of systems with large spins increases as a function of total mass, since these larger systems (particularly those beyond the pulsational-pair instability limit) are predominantly formed through repeated mergers. As the birth spins are increased, the number of 2G mergers (with their characteristically large spins) decreases as their components are more likely to be ejected from the cluster during their first merger. But the total number of BBH systems with non-zero increases, as the first generation of BHs will now form mergers with observable spins. This result is key: one of the most promising ways for identifying a dynamically-formed BBH merger is by the alignment of the spins, with anti-aligned systems () being a clear indicator of dynamical formation Rodriguez et al. (2016c). These results indicate that dynamical assembly in dense star clusters will inevitability produce a merger with , regardless of the BH birth spins.

## Vi Conclusion

We have shown that the inclusion of pN effects can have significant implications for BBH mergers from dense star clusters detectable by LIGO/Virgo. By accounting for GW emission from isolated binaries and during three- and four-body dynamical encounters, we find that a significant number of mergers occur in the cluster, and that about 3% of all mergers (and of in-cluster mergers) in our models will enter the LIGO/Virgo detection band with high residual eccentricity (). Because of this, GCs can potentially produce a significant number of 2G BBH mergers with detectable spins and with masses larger than those produced through the collapse of single stars. Dynamics in dense star clusters can therefore produce BBH mergers with anti-aligned spins (a clear indicator of a dynamical origin) regardless of the initial spins of first-generation BHs: if natal BH spins are large, then GCs can produce BBH mergers with from first-generation systems. If the spins are initially small (as predicted by e.g., Amaro-Seoane and Chen, 2016), then the BBH merger products can often be retained in the cluster, forming a second generation of BBHs with large spins ().

We thank Carl-Johan Haster, Michael Zevin, Johan Samsing, Davide Gerosa, Salvatore Vitale, Chris Pankow, and Scott Hughes for useful discussions. CR is supported by a Pappalardo Postdoctoral Fellowship at MIT. This work was supported by NASA Grant NNX14AP92G and NSF Grant AST-1716762 at Northwestern University. PAS acknowledges support from the Ramón y Cajal Programme of the Ministry of Economy, Industry and Competitiveness of Spain, the COST Action GWverse CA16104, and the CAS President’s International Fellowship Initiative. CR thanks the Niels Bohr Institute for its hospitality while part of this work was completed, and the Kavli Foundation and the DNRF for supporting the 2017 Kavli Summer Program. CR and FR also acknowledge support from NSF Grant PHY-1607611 to the Aspen Center for Physics, where this work was started.

## References

- Abbott et al. (2017a) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, et al., Physical Review Letters 119, 1 (2017a).
- Abbott et al. (2017b) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, et al., Physical Review Letters 118, 221101 (2017b).
- Abbott et al. (2017c) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, et al., Physical Review Letters 119, 30 (2017c).
- Abbott et al. (2016a) B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, et al., Physical Review Letters 116, 061102 (2016a).
- Abbott et al. (2016b) B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, et al., Physical Review Letters 116, 241103 (2016b).
- Belczynski et al. (2010a) K. Belczynski, M. Dominik, T. Bulik, R. OâShaughnessy, C. Fryer, and D. E. Holz, The Astrophysical Journal 715, L138 (2010a).
- Marchant et al. (2016) P. Marchant, N. Langer, P. Podsiadlowski, T. M. Tauris, and T. J. Moriya, Astronomy & Astrophysics 588, A50 (2016).
- Podsiadlowski et al. (2003) P. Podsiadlowski, S. Rappaport, and Z. Han, Mon. Not. R. Astron. Soc 341, 385 (2003).
- Mandel and De Mink (2016) I. Mandel and S. E. De Mink, MNRAS 458, 2634 (2016).
- De Mink and Mandel (2016) S. E. De Mink and I. Mandel, Mon. Not. R. Astron. Soc 460, 3545 (2016).
- Sigurdsson and Hernquist (1993) S. Sigurdsson and L. Hernquist, Nature 364, 423 (1993).
- Portegies Zwart and Mcmillan (2000) S. F. Portegies Zwart and S. L. W. Mcmillan, The Astrophysical Journal 528, 17 (2000).
- OâLeary et al. (2006) R. M. OâLeary, F. A. Rasio, J. M. Fregeau, N. Ivanova, and R. OâShaughnessy, The Astrophysical Journal 637, 937 (2006).
- Miller and Lauburg (2009) M. C. Miller and V. M. Lauburg, The Astrophysical Journal 692, 917 (2009), eprint 0804.2783.
- Downing et al. (2010) J. M. B. Downing, M. J. Benacquista, M. Giersz, and R. Spurzem, Mon. Not. R. Astron. Soc 407, 1946 (2010).
- Banerjee et al. (2010) S. Banerjee, H. Baumgardt, and P. Kroupa, Mon. Not. R. Astron. Soc 402, 371 (2010).
- Downing et al. (2011) J. M. B. Downing, M. J. Benacquista, M. Giersz, and R. Spurzem, Mon. Not. R. Astron. Soc 407, 133 (2011).
- Bae et al. (2014) Y.-B. Bae, C. Kim, and H. M. Lee, Mon. Not. R. Astron. Soc 440, 2714 (2014).
- Ziosi et al. (2014) B. M. Ziosi, M. Mapelli, M. Branchesi, and G. Tormen, Mon. Not. R. Astron. Soc 441, 3703 (2014).
- Rodriguez et al. (2015) C. L. Rodriguez, M. Morscher, B. Pattabiraman, S. Chatterjee, C.-J. Haster, and F. A. Rasio, Physical Review Letters 115, 051101 (2015).
- Rodriguez et al. (2016a) C. L. Rodriguez, S. Chatterjee, and F. A. Rasio, Physical Review D 93, 084029 (2016a).
- Banerjee (2017) S. Banerjee, Mon. Not. R. Astron. Soc 467, 524 (2017).
- Belczynski et al. (2016a) K. Belczynski, D. E. Holz, T. Bulik, and R. OâShaughnessy, Nature 534, 512 (2016a).
- Rodriguez et al. (2016b) C. L. Rodriguez, C.-J. Haster, S. Chatterjee, V. Kalogera, and F. A. Rasio, The Astrophysical Journal 824, L8 (2016b).
- Samsing et al. (2014) J. Samsing, M. MacLeod, and E. Ramirez-Ruiz, The Astrophysical Journal 784, 71 (2014).
- Rodriguez et al. (2016c) C. L. Rodriguez, M. Zevin, C. Pankow, V. Kalogera, and F. A. Rasio, The Astrophysical Journal 832, L2 (2016c).
- Amaro-Seoane and Chen (2016) P. Amaro-Seoane and X. Chen, Mon. Not. R. Astron. Soc 458, 3075 (2016), eprint 1512.04897.
- Woosley (2016) S. E. Woosley, The Astrophysical Journal 836, 244 (2016).
- Fishbach et al. (2017) M. Fishbach, D. E. Holz, and B. Farr, The Astrophysical Journal 840, L24 (2017).
- Gerosa and Berti (2017) D. Gerosa and E. Berti, Physical Review D 95, 124046 (2017).
- Lee (1993) M. H. Lee, Astrophys J v418 418, 147 (1993).
- Lee (1992) M. H. Lee, Ph.D. Thesis (1992).
- Shapiro and Teukolsky (1985a) S. L. Shapiro and S. A. Teukolsky, Astrophys J 298, 34 (1985a).
- Shapiro and Teukolsky (1985b) S. L. Shapiro and S. A. Teukolsky, Astrophys J 292, L41 (1985b).
- Rasio et al. (1989) F. A. Rasio, S. L. Shapiro, and S. A. Teukolsky, Astrophys J 336, L63 (1989).
- Quinlan and Shapiro (1987) G. D. Quinlan and S. L. Shapiro, Astrophys J 321, 199 (1987).
- Quinlan and Shapiro (1989) G. D. Quinlan and S. L. Shapiro, Astrophys J 343, 725 (1989).
- Quinlan and Shapiro (1990) G. D. Quinlan and S. L. Shapiro, Astrophys J 356, 483 (1990).
- Kupi et al. (2006) G. Kupi, P. Amaro-Seoane, and R. Spurzem, Mon. Not. R. Astron. Soc 371, L77 (2006), eprint astro-ph/0602125.
- Brem et al. (2013) P. Brem, P. Amaro-Seoane, and R. Spurzem, Mon. Not. R. Astron. Soc 434, 2999 (2013), eprint 1302.3135.
- Planck Collaboration et al. (2015) Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, et al. (2015).
- Joshi et al. (2000) K. J. Joshi, F. A. Rasio, S. P. Zwart, and S. Portegies Zwart, The Astrophysical Journal 540, 969 (2000).
- Pattabiraman et al. (2013) B. Pattabiraman, S. Umbreit, W.-k. Liao, A. Choudhary, V. Kalogera, G. Memik, and F. A. Rasio, The Astrophysical Journal Supplement Series 204, 15 (2013).
- Hénon (1971) M. Hénon, Astrophysics and Space Science 14, 151 (1971).
- Henon (1971) M. Henon, Astrophysics and Space Science 13, 284 (1971).
- Hurley et al. (2000) J. R. Hurley, O. R. Pols, and C. A. Tout, Mon. Not. R. Astron. Soc 315, 543 (2000).
- Hurley et al. (2002) J. R. Hurley, C. A. Tout, and O. R. Pols, Mon. Not. R. Astron. Soc 329, 897 (2002).
- Chatterjee et al. (2010) S. Chatterjee, J. M. Fregeau, S. Umbreit, and F. A. Rasio, The Astrophysical Journal 719, 915 (2010).
- Morscher et al. (2013) M. Morscher, S. Umbreit, W. M. Farr, and F. A. Rasio, The Astrophysical Journal 763, L15 (2013).
- Fregeau et al. (2004) J. M. Fregeau, P. Cheung, S. F. Portegies Zwart, and F. A. Rasio, Mon. Not. R. Astron. Soc 352, 1 (2004).
- Fregeau and Rasio (2007) J. M. Fregeau and F. A. Rasio, The Astrophysical Journal 658, 1047 (2007).
- Rodriguez et al. (2016d) C. L. Rodriguez, M. Morscher, L. Wang, S. Chatterjee, F. A. Rasio, and R. Spurzem, Mon. Not. R. Astron. Soc 463, 2109 (2016d).
- Wang et al. (2016) L. Wang, R. Spurzem, S. Aarseth, M. Giersz, A. Askar, P. Berczik, T. Naab, R. Schadow, and M. B. N. Kouwenhoven, Mon. Not. R. Astron. Soc 458, 1450 (2016).
- Vink et al. (2001) J. S. Vink, A. de Koter, and H. J. G. L. M. Lamers, Astronomy and Astrophysics 369, 574 (2001).
- Belczynski et al. (2010b) K. Belczynski, T. Bulik, C. L. Fryer, A. Ruiter, F. Valsecchi, J. S. Vink, and J. R. Hurley, The Astrophysical Journal 714, 1217 (2010b).
- Fryer et al. (2012) C. L. Fryer, K. Belczynski, G. Wiktorowicz, M. Dominik, V. Kalogera, and D. E. Holz, The Astrophysical Journal 749, 91 (2012).
- Antognini et al. (2014) J. M. Antognini, B. J. Shappee, T. A. Thompson, and P. Amaro-Seoane, Mon. Not. R. Astron. Soc 439, 1079 (2014).
- Peters (1964) P. Peters, Physical Review 136, B1224 (1964).
- Barausse et al. (2012) E. Barausse, V. Morozova, and L. Rezzolla, The Astrophysical Journal 758, 63 (2012).
- Kesden (2008) M. Kesden, Physical Review D - Particles, Fields, Gravitation and Cosmology 78, 084030 (2008).
- Lousto and Zlochower (2014) C. O. Lousto and Y. Zlochower, Physical Review D 89, 104052 (2014).
- Berti et al. (2007) E. Berti, V. Cardoso, J. A. Gonzalez, U. Sperhake, M. Hannam, S. Husa, and B. Brügmann, Physical Review D 76, 064034 (2007).
- Campanelli et al. (2007) M. Campanelli, C. Lousto, Y. Zlochower, and D. Merritt, The Astrophysical Journal 659, L5 (2007).
- Tichy and Marronetti (2008) W. Tichy and P. Marronetti, Physical Review D 78, 081501 (2008).
- Barausse and Rezzolla (2009) E. Barausse and L. Rezzolla, The Astrophysical Journal 704, L40 (2009).
- Buonanno et al. (2008) A. Buonanno, L. E. Kidder, and L. Lehner, Physical Review D 77, 026004 (2008).
- Rezzolla et al. (2008) L. Rezzolla, E. Barausse, E. N. Dorband, D. Pollney, C. Reisswig, J. Seiler, and S. Husa, Physical Review D 78, 044002 (2008).
- González et al. (2007) J. A. González, U. Sperhake, B. Brügmann, M. Hannam, and S. Husa, Physical Review Letters 98, 091101 (2007).
- Lousto and Zlochower (2008) C. O. Lousto and Y. Zlochower, Physical Review D 77, 044028 (2008).
- Lousto et al. (2012) C. O. Lousto, Y. Zlochower, M. Dotti, and M. Volonteri, Physical Review D 85, 084015 (2012).
- Lousto and Zlochower (2013) C. O. Lousto and Y. Zlochower, Physical Review D 87, 084027 (2013).
- Wen (2003) L. Wen, The Astrophysical Journal 598, 419 (2003).
- Blanchet (2006) L. Blanchet, Living Reviews in Relativity 9, 4 (2006).
- Sopuerta et al. (2007) C. F. Sopuerta, N. Yunes, and P. Laguna, Astrophys. J. Lett. 656, L9 (2007), eprint astro-ph/0611110.
- Yunes and Berti (2008) N. Yunes and E. Berti, Phys. Rev. D. 77, 124006 (2008), eprint 0803.1853.
- Gerosa and Kesden (2016) D. Gerosa and M. Kesden, Physical Review D 93, 124066 (2016).
- Jiménez-Forteza et al. (2017) X. Jiménez-Forteza, D. Keitel, S. Husa, M. Hannam, S. Khan, and M. Pürrer, Phys. Rev. D 95, 064024 (2017), eprint 1611.00332.
- Askar et al. (2016) A. Askar, M. Szkudlarek, D. Gondek-Rosińska, M. Giersz, and T. Bulik, Mon. Not. R. Astron. Soc. Lett. 464, L36 (2016).
- Antonini and Rasio (2016) F. Antonini and F. A. Rasio, The Astrophysical Journal 831, 187 (2016), eprint 1606.04889.
- Giesler et al. (2017) M. Giesler, D. Clausen, and C. D. Ott (2017), eprint arXiv: astro-ph/1708.05915.
- Belczynski et al. (2016b) K. Belczynski, A. Heger, W. Gladysz, A. J. Ruiter, S. Woosley, G. Wiktorowicz, H.-Y. Chen, T. Bulik, R. OâShaughnessy, D. E. Holz, et al., Astronomy & Astrophysics 594, A97 (2016b).
- Fishbach and Holz (2017) M. Fishbach and D. E. Holz (2017).
- Heggie (1975) D. C. Heggie, Mon. Not. R. Astron. Soc 173, 729 (1975).
- Samsing and Ramirez-Ruiz (2017) J. Samsing and E. Ramirez-Ruiz, The Astrophysical Journal Letters 840, L14 (2017).
- Samsing (2017) J. Samsing (2017), eprint arXiv: astro-ph/1711.07452.
- Benacquista and Downing (2013) M. J. Benacquista and J. M. B. Downing, Living Reviews in Relativity 16, 4 (2013), eprint 1110.4423.
- Morscher et al. (2015) M. Morscher, B. Pattabiraman, C. Rodriguez, F. A. Rasio, and S. Umbreit, The Astrophysical Journal 800, 9 (2015).
- Abbott et al. (2016c) B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, et al., Living Rev Relativ 19, 1 (2016c), eprint 1304.0670.
- Chen et al. (2017) H.-Y. Chen, D. E. Holz, J. Miller, M. Evans, S. Vitale, and J. Creighton (2017), eprint arXiv: astro-ph/1709.08079.
- Miller and Miller (2015) M. C. Miller and J. M. Miller, Physics Reports 548, 1 (2015).
- Lousto and Zlochower (2011) C. O. Lousto and Y. Zlochower, Physical Review Letters 107, 231102 (2011), eprint 1108.2009.
- Berti and Volonteri (2008) E. Berti and M. Volonteri, Astrophys. J. 684, 822-828 (2008), eprint 0802.0025.
- Lousto et al. (2010) C. O. Lousto, H. Nakano, Y. Zlochower, and M. Campanelli, Phys Rev D 81 (2010), eprint 0910.3197.
- Kesden et al. (2010) M. Kesden, U. Sperhake, and E. Berti, Phys. Rev. D. 81, 084054 (2010), eprint 1002.2643.
- Ajith et al. (2011) P. Ajith, M. Hannam, S. Husa, Y. Chen, B. Brügmann, N. Dorband, D. Müller, F. Ohme, D. Pollney, C. Reisswig, et al., Physical Review Letters 106, 241101 (2011), eprint 0909.2867.
- Pürrer et al. (2016) M. Pürrer, M. Hannam, and F. Ohme, Phys. Rev. D 93, 084042 (2016), eprint 1512.04955.
- King (1966) I. R. King, The Astronomical Journal 71, 64 (1966).
- Kroupa (2001) P. Kroupa, Mon. Not. R. Astron. Soc 322, 231 (2001).
- Spera and Mapelli (2017) M. Spera and M. Mapelli, Mon. Not. R. Astron. Soc 4, 4739 (2017).

## Appendix A Post-Newtonian Physics

Our post-Newtonian (pN) scattering code was originally presented in detail in Antognini et al. (2014); Amaro-Seoane and Chen (2016), building upon the well-tested and commonly used fewbody code Fregeau et al. (2004). fewbody computes strong binary-single and binary-binary encounters using an -order RungeâKutta Prince-Dormand integrator. Because the pN corrections to the equations of motion can be expressed in terms of pair-wise position- and velocity-dependent accelerations (see e.g., Blanchet, 2006), these terms are easily incorporated into fewbody. To incorporate both conservative and dissipative relativistic effects (i.e., GW emission), the pN force terms are included in the code up to 3.5pN order, though for the current study we only use the 1pN, 2pN, and 2.5pN terms. This code has been shown to conserve energy when dissipative effects are ignored, and to reproduce accurately the inspiral times predicted by Peters (1964). See Appendices A1 and A2 in Antognini et al. (2014).

We detect a merger of two BHs in fewbody using a standard sticky-sphere collision criterion where the effective radii of the BHs are set to 5 Schwarzschild radii (), since the pN approximation begins to break down at separations of (e.g., Yunes and Berti, 2008). Any two BHs touching within these radii are assumed to merge. For BBH mergers that occur during a fewbody integration, we determine the properties of the merger remnant using the fitting formulas compiled in Gerosa and Kesden (2016), based on extensive work to parameterize the final mass (using the interpolation in Barausse et al. (2012) based on results from Kesden (2008); Lousto and Zlochower (2014); Berti et al. (2007); Tichy and Marronetti (2008)), final spin (using the interpolation in Barausse and Rezzolla (2009) based on results from Kesden (2008); Buonanno et al. (2008); Lousto and Zlochower (2014); Tichy and Marronetti (2008); Rezzolla et al. (2008), though see Jiménez-Forteza et al. (2017) for more recent fits), and the recoil kick (using the kick model from Campanelli et al. (2007) employing coefficients from González et al. (2007); Lousto and Zlochower (2008); Lousto et al. (2012); Lousto and Zlochower (2013)). Due to a typo in the original published version, we reproduce the correct formula for the final mass here, while referring the reader to Gerosa and Kesden (2016) for the prescriptions for the final spin and recoil kick. We note that the precession software package described in Gerosa and Kesden (2016) does not contain this error. For the final mass, we use the interpolated fit between the test-particle limits Kesden (2008) and numerical relativity Lousto and Zlochower (2014); Berti et al. (2007) that was given in Barausse and Rezzolla (2009). This equation is incorrect in Gerosa and Kesden (2016), and should read

(3) |

where is the ratio of the remnant mass to the binary mass, , , (from Barausse and Rezzolla, 2009) and is the energy at the innermost-stable-circular orbit for a Kerr BH:

(4) | ||||

(5) | ||||

(6) | ||||

(7) |

and where is:

(8) |

Note that is not the same as in equation (2) in the main text. We track the spin magnitudes of any BHs and update them accordingly when a merger occurs. The spin directions during mergers (needed for the computation of the remnant mass, spin, and kick) are assumed to be isotropically distributed on the sphere. The recoil kick from Campanelli et al. (2007) is decomposed into a component perpendicular to the orbital angular momentum (to which both an asymmetric mass ratio and the spins will contribute) and a component parallel to to the angular momentum, which arises only from the spins. We draw a random angle in the plane of the binary for the perpendicular component and apply the two kicks in this coordinate frame.

During integrations of three- and four-body encounters, we do not stop the integration when two BHs merge, but instead continue until fewbody’s standard termination criteria (which check for unbound collections of objects and dynamical stability) are satisfied. This allows us to self-consistently track the outcome of mergers during dynamical encounters, which can be important for the retention of BH merger products.

While we are applying the kicks self-consistently during the fewbody integration, we are not taking into account the residual eccentricity of the BBHs which merge during scattering encounters. As eccentric, non-spinning mergers receive higher kicks than circular BBH mergers (e.g., Sopuerta et al., 2007), it is entirely possible that we are overestimating the retention of mergers which occur during three- and four-body encounters. We do not implement these eccentric corrections here, since there exist no corresponding fitting formulae for eccentric BBH mergers with spinning components.

For BH – (non BH) star mergers, we do not consider any accretion of material from the star, and assume that the entire star is tidally disrupted. This is a highly-conservative assumption, but is safest without a better treatment for accretion and spin-up of compact objects (which, while potentially interesting, is beyond the scope of this study).

For isolated BBHs in the cluster, CMC had previously relied on BSE to determine the change in eccentricity () and semi-major axis () at each timestep. However, the evolution of and , when averaged over the binary orbit, depends strongly on the eccentricity Peters (1964):

(9) | ||||

(10) |

As BSE uses a forward-Euler integration for its equations, it can underestimate the merger time for binaries that reach very high eccentricities. As an example, BSE will overestimate the merger time of a binary with a semi-major axis of 0.1 AU and an eccentricity of 0.9 by nearly 20%. As stated in the main text, the orbital eccentricity of a BBH is critical for the proper treatment of in-cluster mergers. For isolated binaries in the cluster, we directly integrate equations (9) and (10) using an -order Runge-Kutta Prince-Dormand integrator. This is done outside of the BSE package, so that any merging binaries use the same final mass, final spin, and recoil fitting-formula that are employed for mergers in the fewbody integrator. As with mergers occurring during encounters, we terminate the integration when the effective radii of the two BHs (5 Schwarzschild radii) touch. To determine the eccentricity for these systems in Figure 1, we integrate , found by dividing equations (9) and (10). We integrate from the last-reported and in the cluster for each binary (which, for mergers in fewbody, is the and at the point of contact) to the semi-major axis corresponding to a circular GW frequency Hz (where ). Although fitting formulae exist for the peak GW frequency of eccentric mergers (e.g., Wen, 2003), we do not use these here, as several of these GW-capture BBHs form inside the LIGO/Virgo band, and therefore did not pass through the eccentric 10Hz threshold as an isolated binary.

## Appendix B GC Models

We evolve a 4x3x2 grid of GC models with different initial masses, metallicities (assumed correlated with galactocentric distance), and virial radii. These clusters have identical initial conditions to those considered and discussed in more detail in Rodriguez et al. (2016a). We consider clusters with , , , and particles (single or binary stars) initially, covering a realistic range of GC masses. The initial positions and velocities are distributed in phase space following a standard King model King (1966) with concentration . The virial radii of the clusters are then scaled to or pc, representing a realistic range of initial sizes. The individual stellar masses are drawn from the range to assuming a Kroupa initial-mass function (IMF) Kroupa (2001):

(11) |

where

(12) |

As in our previous work we consider metallicities and galactocentric distances of , , and , following the correlation of metallicities with distance observed for GCs in the Milky Way. Finally, we select 10% of particles to be binaries with separations consistent with Öpek’s law (flat in from the point of stellar contact to the hard–soft boundary, where the orbital velocity equals the typical velocity of particles in the cluster) and thermal eccentricities (). The primary masses are drawn from the IMF, while the secondary masses are drawn from a uniform distribution, .

Our stellar evolution prescriptions are almost identical to those adopted in Rodriguez et al. (2016a). We have added new physics describing mass loss from pulsational-pair instabilities and pair-instability supernovae. Briefly, we follow Belczynski et al. (2016b), and assume that any star whose pre-explosion helium core mass is between and will undergo pulsations that eject a significant amount of mass, until the final product is at most . This yields a BH with a mass of , assuming that 10% of the baryonic mass is lost during the conversion from baryonic to gravitational mass. Stars with He core masses above are completely destroyed in pair-instability supernovae. There are significant uncertainties on these numbers, and is likely a conservative limit (with many studies suggesting for the lower limit of the mass gap (e.g., Woosley, 2016; Spera and Mapelli, 2017)). However, our results here are insensitive to the specific limit adopted for stellar-mass BH remnants. BHs in GCs will always undergo mass segregation and form binaries from the most massive BHs first, regardless of how large their masses are. Furthermore, the GW recoil kicks applied to merging BBHs are independent of the total mass (depending only on the mass ratio), implying that the retention of merger products also does not depend on the total mass.

One caveat must be addressed here. By selecting as our maximum BH mass, we are create a population of mergers with exactly equal-mass components (). This may artificially increase the retention of BBH merger products, as equal-mass, non-spinning BBHs receive zero kicks. We find that 40 of the 1521 in-cluster mergers are mergers. However, because these BBHs are the most massive in the cluster, they form and merge very early in the cluster lifetime, with the latest merger in any model occurring after cluster formation (at ). GCs were significantly more massive and compact in the early universe, with correspondingly deeper central potentials and escape velocities. The lowest escape speed where a BBH merger occurred in our models was (while of the mergers occurred where ). For non-spinning BHs and our adopted kick prescription, this would ensure the retention of all merger products assuming the true mass ratio was greater than .