Post-Newtonian Dynamics in Dense Star Clusters: Highly-Eccentric, Highly-Spinning, and Repeated Binary Black Hole Mergers

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; Abbott et al. 2017b; Abbott et al. 2017c; Abbott et al. 2016a; Abbott et al. 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 (Belczynski et al. 2010a; Marchant et al. 2016; Podsiadlowski et al. 2003; Mandel and De Mink 2016; De Mink and Mandel 2016, e.g.) and dynamical formation in dense star clusters (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, e.g.) 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 (Amaro-Seoane and Chen 2016, see).

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 $\lesssim 50M_{\odot}$ 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 $N$-body dynamics has been studied previously for highly idealized systems (Lee 1993; Lee 1992; Shapiro and Teukolsky 1985a; Shapiro and Teukolsky 1985b; Rasio et al. 1989; Quinlan and Shapiro 1987; Quinlan and Shapiro 1989; Quinlan and Shapiro 1990; Kupi et al. 2006; Brem et al. 2013, e.g.,) or open clusters (Banerjee 2017, e.g.,), 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 ($\sim 10\%$) 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 $\sim$20% of BBH mergers from our models having components formed in a previous merger.

Throughout this paper, we assume a flat $\Lambda\rm{CDM}$ cosmology with $h=0.679$ and $\Omega_{M}=0.3065$ 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; Hurley et al. 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 $N$-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 (Belczynski et al. 2016a, e.g.).

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 ($\chi_{b}=0$, where $\chi$ 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 (Hurley et al. 2002, BSE, ). By default, BSE only applies GW energy loss to binaries with $a<10R_{\odot}$. 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 (Rodriguez et al. 2015; Rodriguez et al. 2016a; Rodriguez et al. 2016b; Downing et al. 2010; Askar et al. 2016, e.g.,), where ejected BBHs dominated the merger rate in the local universe. Semi-analytic approaches to cluster dynamics (Antonini and Rasio 2016; Giesler et al. 2017, e.g.,) 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 $\sim 45M_{\odot}$ which is reduced to $\sim 40M_{\odot}$ 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 $40M_{\odot}$ 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; Rodriguez et al. 2016b, 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 (Rodriguez et al. 2016a, 0.06% of mergers at $z<1$, see), 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 ($z<1$), 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 $z<1$ decreases by $\sim 20\%$ (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 $\sim 50\%$.

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 $30M_{\odot}+30M_{\odot}$ BBH is ejected from a GC with $a\sim 0.4~\rm{AU}$ (roughly 10 times greater than the $a<10R_{\odot}$ cutoff in BSE) after undergoing $\mathcal{O}(10)$ 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, $p(e)de=2e~de$ 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:

As (1) makes clear, a large eccentricity can significantly decrease the merger timescale. For $e\gtrsim 0.95$ (roughly %10 post-encounter binaries) $t_{\rm{GW}}$ will decrease by more than $10^{3}$, 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 $(1-e^{2})^{7/2}$ 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 ($e>0.1$) 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 ($\sim 3\%$ of all mergers) at $z<1$ 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 $\sim 10\rm{Gyr}$ (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 $5M_{\odot}$ to $40M_{\odot}$. 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 (Fishbach and Holz 2017, the detection horizon scales with the mass of the more massive component as $m^{2.2}$,). At the expected sensitivity for Advanced LIGO’s third observing run Abbott et al. 2016c, a BBH with component masses of $40M_{\odot}+80M_{\odot}$ could be detected out to $z\sim 1$, 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 (Abbott et al. 2016b, GW151226, ) 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 (Campanelli et al. 2007; Lousto and Zlochower 2011; Lousto et al. 2012, e.g.,), 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 $\chi_{b}=0$, we find that $\sim$20% of mergers at $z<1$ are 2G mergers. As the spins are increased, this number decreases, and once $\chi_{b}=0.4$, we observe $\mathcal{O}(1)$ 2G mergers, compared to the 672 first-generation mergers which occur at $z<1$.

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 $N$-body simulations with spins Brem et al. 2013, repeated mergers of BBHs in clusters with near-equal masses will tend to produce BHs with $\chi\sim 0.7$, (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:

where $\hat{L}$ is the direction of the orbital angular momentum and $\vec{\chi}_{1,2}$ 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 $\chi_{\rm{eff}}=0$ with symmetric tails whose extent depends on the BH spin magnitudes. We show the distributions of $\chi_{\rm{eff}}$ 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 $\chi_{\rm{eff}}$ 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 ($\chi_{\rm{eff}}<0$) 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 $\chi_{\rm{eff}}<0$, 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 $\sim 10\%$ of in-cluster mergers) in our models will enter the LIGO/Virgo detection band with high residual eccentricity ($e>0.1$). 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 $\chi_{\rm{eff}}<0$ from first-generation systems. If the spins are initially small (Amaro-Seoane and Chen 2016, as predicted by e.g.,), then the BBH merger products can often be retained in the cluster, forming a second generation of BBHs with large spins ($\chi\sim 0.7$).

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

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 $8^{th}$-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 (Blanchet 2006, see e.g.,), 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 ($10Gm/c^{2}$), since the pN approximation begins to break down at separations of $r\sim 10M$ (Yunes and Berti 2008, e.g.,). 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

where $M_{f}/M$ is the ratio of the remnant mass to the binary mass, $\eta\equiv m_{1}m_{2}/(m_{1}+m_{2})^{2}$, $p_{0}=0.04827$, $p_{1}=0.01707$ (Barausse and Rezzolla 2009, from) and $E_{\rm{ISCO}}$ is the energy at the innermost-stable-circular orbit for a Kerr BH:

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 (Sopuerta et al. 2007, e.g.,), 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 ($e$) and semi-major axis ($a$) at each timestep. However, the evolution of $a$ and $e$, when averaged over the binary orbit, depends strongly on the eccentricity Peters 1964:

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 $30M_{\odot}+30M_{\odot}$ 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 $8^{th}$-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 $de/da$, found by dividing equations (9) and (10). We integrate from the last-reported $a$ and $e$ in the cluster for each binary (which, for mergers in fewbody, is the $a$ and $e$ at the point of contact) to the semi-major axis corresponding to a circular GW frequency $f_{\rm{GW}}=10\,$Hz (where $f_{\rm{GW}}=2f_{\rm{orbital}}$). Although fitting formulae exist for the peak GW frequency of eccentric mergers (Wen 2003, e.g.,), 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 $2\times 10^{5}$, $5\times 10^{5}$, $10^{6}$, and $2\times 10^{6}$ 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 $W_{0}=5$. The virial radii of the clusters are then scaled to $1$ or $2$ pc, representing a realistic range of initial sizes. The individual stellar masses are drawn from the range $0.08\,M_{\odot}$ to $150\,M_{\odot}$ assuming a Kroupa initial-mass function (IMF) Kroupa 2001:

As in our previous work we consider metallicities $Z$ and galactocentric distances $R_{g}$ of $Z=0.005/R_{g}=2\,\rm{kpc}$, $Z=0.001/R_{g}=8\,\rm{kpc}$, and $Z=0.0002/R_{g}=20\,\rm{kpc}$, 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 $\log(a)$ 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 ($p(e)de=2e\,de$). The primary masses are drawn from the IMF, while the secondary masses are drawn from a uniform distribution, $m_{2}\in U[0,\,m_{1}]$.

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 $45\,M_{\odot}$ and $65\,M_{\odot}$ will undergo pulsations that eject a significant amount of mass, until the final product is at most $45M_{\odot}$. This yields a BH with a mass of $40.5\,M_{\odot}$, assuming that 10% of the baryonic mass is lost during the conversion from baryonic to gravitational mass. Stars with He core masses above $65M_{\odot}$ are completely destroyed in pair-instability supernovae. There are significant uncertainties on these numbers, and $40\,M_{\odot}$ is likely a conservative limit (with many studies suggesting $\sim 50\,M_{\odot}$ for the lower limit of the mass gap (Woosley 2016; Spera and Mapelli 2017, e.g.,)). 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 $40.5M_{\odot}$ as our maximum BH mass, we are create a population of mergers with exactly equal-mass components ($40.5M_{\odot}+40.5M_{\odot}$). 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 $40.5M_{\odot}+40.5M_{\odot}$ 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 $\sim 1\rm{Gyr}$ after cluster formation (at $z\sim 2.3$). 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 $40.5M_{\odot}+40.5M_{\odot}$ BBH merger occurred in our models was $31~\rm{km}/\rm{s}$ (while $\sim 90\%$ of the mergers occurred where $v_{\rm{esc}}>50~\rm{km}/\rm{s}$). For non-spinning BHs and our adopted kick prescription, this would ensure the retention of all $40.5M_{\odot}+40.5M_{\odot}$ merger products assuming the true mass ratio was greater than $\sim 0.9$.