# The star clusters that make black hole binaries across cosmic time

Nick Choksi, Marta Volonteri, Monica Colpi, Oleg Y. Gnedin, Hui Li

► To cite this version:

Nick Choksi, Marta Volonteri, Monica Colpi, Oleg Y. Gnedin, Hui Li. The star clusters that make black hole binaries across cosmic time. The Astrophysical Journal, 2019, 873 (1), pp.100. 10.3847/1538-4357/aaffde. hal-01880796

HAL Id: hal-01880796

https://hal.science/hal-01880796v1

Submitted on 13 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Star Clusters That Make Black Hole Binaries across Cosmic Time

Nick Choksi<fnref n="1" /> Ⓓ, Marta Volonteri<fnref n="2" /> Ⓓ, Monica Colpi<fnref n="3" /> Ⓓ, Oleg Y. Gnedin<fnref n="4" /> Ⓓ, and Hui Li<fnref n="5" /> Ⓓ

Received 2018 August 16; revised 2018 December 26; accepted 2019 January 16; published 2019 March 7

## Abstract

We explore the properties of dense star clusters that are likely to be nurseries of stellar black holes pairing in close binaries. We combine a cosmological model of globular cluster formation with analytic prescriptions for the dynamical assembly and evolution of black hole binaries (BHBs) to constrain which types of clusters are most likely to form binaries tight enough to coalesce within a Hubble time. We find that BHBs that are ejected and later merge ex situ form in clusters of a characteristic mass $M_{\text{cl}} \sim 10^{5.3} M_{\odot}$, whereas binaries that merge in situ form in more massive clusters, $M_{\text{cl}} \sim 10^{5.7} M_{\odot}$. The clusters that dominate the production of BHBs are similar in age and metallicity to the entire population. Finally, we estimate an approximate cosmic black hole merger rate of dynamically assembled binaries using the mean black hole mass for each cluster, given its metallicity. We find an intrinsic rate of $\sim 6 \text{ Gpc}^{-3} \text{ yr}^{-1}$ at $z=0$, a weakly increasing merger rate out to $z=1.5$, and then a decrease out to $z=4$. Our results can be used to provide a cosmological context and choose initial conditions in numerical studies of BHBs assembled in star clusters.

**Key words:** globular clusters: general – gravitational waves – stars: black holes – stars: kinematics and dynamics

## 1. Introduction

The detections of gravitational wave (GW) emission from ten distinct black hole (BH) binaries (BHBs) have conclusively demonstrated both that stellar BHBs form in nature, and that they can coalesce within a Hubble time (Abbott et al. 2016a, 2016b, 2017a, 2017b, 2017c; The LIGO Scientific Collaboration et al. 2018). However, dissipation of orbital energy and angular momentum via GW emission is only significant at very small separations ($\lesssim 10R_{\odot}$ for circular orbits). While it is possible that binaries form with separations small enough for GW emission to result in coalescence within a Hubble time, such a formation scenario would require that the parent stars not expand significantly during stellar evolution, so as to avoid a stellar merger before both objects collapse into BHs (e.g., de Mink & Mandel 2016; Mandel & de Mink 2016; Marchant et al. 2016). Alternatively, binaries may form with wide separations and shrink by a few orders of magnitude after a common envelope phase of stellar evolution (e.g., Dominik et al. 2012; Ivanova et al. 2013; Postnov & Yungelson 2014; Belczynski et al. 2016).

Yet another possibility, which does not depend on the details of stellar evolution, is that BHBs form in dense star clusters, such as in globular cluster (GC) progenitors (Morscher et al. 2015; Rodriguez et al. 2015, 2016a, 2016b) or in the more massive nuclear star clusters (NSCs; Portegies Zwart & McMillan 2000; Miller & Lauburg 2009). Because stellar densities are much higher in clusters than in the field, and because BHs segregate in the cluster core (Spitzer 1969), the rate of three-body interactions is enhanced significantly. As a result, BHBs can be assembled dynamically, either via direct three-body interactions (Lee 1995) or via exchange interactions in close encounters with binary stars which may ultimately lead to the formation of a BHB (Goodman & Hut 1993). After their assembly, further hardening via binary-single scatterings provides a mechanism by which BHBs can contract to smaller

separations. Depending on the magnitude of the velocity dispersion of the surrounding stars, the BHBs can either be ejected by dynamical recoil or be retained inside the cluster.

This dynamical formation channel has been explored extensively, both analytically (e.g., Lee 1995; Breen & Heggie 2013a, 2013b) and numerically via Monte Carlo (Gürkan et al. 2004; Morscher et al. 2015; Rodriguez et al. 2015, 2016a, 2016b; Antonini & Rasio 2016; Chatterjee et al. 2017a, 2017b; Hong et al. 2018) and N-body simulations (Lee 1995; O'Leary et al. 2007; Miller & Lauburg 2009; Bae et al. 2014; Park et al. 2017). These studies have confirmed that GCs can indeed be the formation sites for a large population of BHBs (Rodriguez et al. 2016a; see also Benacquista & Downing 2013 for a review on the dynamical assembly of binaries).

Thus far, previous works have lacked detailed information about the initial conditions of the clusters being modeled. They typically adopt initial conditions based on observations of the Galactic GC population. However, GCs are expected to evolve significantly over $\sim 13$ Gyr of cosmic history. Not only do the current properties of GCs differ from their state at formation, but also a large population of clusters that hosted BHB mergers may have been either disrupted by the local Galactic tidal field (Gnedin et al. 1999, 2014; Fall & Zhang 2001; Kravtsov & Gnedin 2005; Gieles & Baumgardt 2008) or spiraled into the central NSC before the present day (e.g., Tremaine et al. 1975; Capuzzo-Dolcetta & Miocchi 2008a, 2008b; Antonini et al. 2012; Antonini 2013; Gnedin et al. 2014).

In this work, we provide constraints on the properties of the clusters most likely to be the nurseries of BHBs. To do so, we extract clusters from a cosmological model of GC formation and apply analytic estimates of relevant dynamical timescales. Our results are complementary to the study of BHB dynamics via more computationally intensive methods which are limited

<footnote n="1">¹University of California, Berkeley, Department of Astronomy, Campbell Hall, Berkeley, CA 94704, USA; nchoksi@berkeley.edu</footnote>

<footnote n="2">²Institut d'Astrophysique de Paris, UMR 7095 CNRS, Université Pierre et Marie Curie, 98bis Blvd Arago, F-75014 Paris, France</footnote>

<footnote n="3">³Department of Physics, University of Milano Bicocca and Istituto Nazionale di Fisica Nucleare-Sezione di Milano Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="4">⁴Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="5">⁵Department of Physics, Kavli Institute for Astrophysics and Space Research, MIT, Cambridge, MA 02139, USA</footnote>
to modeling a small subset of the parameter space of initial conditions.

## 2. Methodology

We begin this section with a brief outline of our GC formation and evolution model. Then we describe the cluster properties relevant to the evolution of BHBs and how we obtain each property (Section 2.1). Given these cluster properties, we list our analytic prescriptions for the treatment of BHB formation and hardening (Section 2.2). We finish with a description of the timescales over which the cluster may be disrupted (Section 2.3). We adopt a flat, $\Lambda$CDM cosmology with $h = 0.7$ and $\Omega_{\Lambda} = 0.7$. We use “cosmic time” to refer to the time from the big bang to a redshift $z$, and define the cosmic time at $z = 0$ as $t_H \approx 13.7$ Gyr. Throughout, we refer to the initial mass of the cluster as $M_{cl}$ and the mass at any other time as $M_{cl}(t)$.

### 2.1. Determination of Cluster Properties

We use the properties of GCs from the model of Choksi et al. (2018, hereafter CGL18). The model assumes that GCs form in periods of rapid accretion onto dark matter halos (e.g., major mergers of galaxies). When such events occur, the total mass that forms in GCs scales linearly with the cold gas mass in the galaxy, which is in turn set by empirical galactic scaling relations. Clusters are drawn from a power-law cluster initial mass function (IMF) $dN/dM_{cl} \propto M_{cl}^{-2}$ (see Choksi & Gnedin 2018 for a discussion of the impact of the adopted mass function on the properties of GC systems) and are assigned metallicities based on the metallicity of their host galaxy using an empirical galaxy stellar mass–metallicity relation. GCs are evolved using prescriptions for both dynamical disruption and stellar evolution. The model was applied to $\approx 200$ halo merger trees in the mass range $10^{11} M_{\odot} \lesssim M_h \lesssim 10^{15} M_{\odot}$ from the Illustris dark-matter-only simulation (Vogelsberger et al. 2014; Nelson et al. 2015). The resulting GC populations at $z = 0$ are shown in CGL18 to match a wide range of the observed properties of GC systems. As a reference, we note that in the fiducial model, half of all clusters are predicted to form in the range $5 < z < 2.3$, corresponding to ages of 10.8–12.5 Gyr, in halos of masses between $10^{11}$ and $10^{12.5} M_{\odot}$. These predictions are consistent with other recent GC formation models (Pfeffer et al. 2018; El-Badry et al. 2019).

We extract directly from the CGL18 model the initial cluster masses, [Fe/H] values, formation times, and host galaxy properties for all clusters formed in the model. To relate [Fe/H] to a total metallicity, we apply the conversion $\log_{10} Z/Z_{\odot} \approx [Fe/H] + 0.2$, based on the simulations of Ma et al. (2016). Below, we detail how we draw the remaining relevant cluster properties that the CGL18 model does not set.

The location of the cluster within the host galaxy is important because clusters may either inspiral due to dynamical friction and merge into the central NSC or be gradually disrupted by the local tidal field before BHBs can be dynamically assembled. Observations of H$\alpha$ emission in high-redshift star-forming galaxies show that most of the star formation occurs near the galaxy’s half-mass-radius, $R_h$ (Förster Schreiber et al. 2018). We therefore assume that GCs form with galactocentric radii distributed uniformly in the range $(0.5-2)R_h$. We correlate $R_h$ to the scale radius $R_d$ of an exponential disk that has the same specific angular momentum

as the dark matter halo. This results in a value of $R_d \approx 2^{-1/2}\lambda R_{vir}$, where $R_{vir}$ is the virial radius and $\lambda$ is the dimensionless spin parameter of the host halo at formation (Mo et al. 1998). The scale length relates to the half-mass-radius as $R_d = 0.58 R_h$. We draw $\lambda$ from a log-normal distribution centered on $\lambda = 0.04$ with scatter $\sigma_{\lambda} = 0.25$ dex, typical for cosmological $\Lambda$CDM simulations (e.g., Rodríguez-Puebla et al. 2016).

The initial size of the cluster will strongly affect the BHB merger rate (e.g., Morscher et al. 2015) because it sets the density, and in turn the BH interaction rate. Furthermore, GCs are expected to evolve in size due to relaxation effects and tidal stripping (e.g., Gnedin et al. 1999; Gieles & Baumgardt 2008; Muratov & Gnedin 2010), so GC radii today cannot be used to estimate the initial properties of the cluster. Therefore, we turn to observations of “young massive clusters” in the local universe, which have masses and sizes consistent with objects that could evolve into GC-like systems after a few Gyr of dynamical and stellar evolution (Chatterjee et al. 2010, 2013; Portegies Zwart et al. 2010; Bastian 2016). These clusters show an approximately log-normal distribution in both their half-light radii $r_h$ (which we take as a proxy for the half-mass-radius) and their core radii $r_c$, with peaks at $r_h \approx 2.8$ pc and $r_c \approx 1$ pc (Portegies Zwart et al. 2010; Bastian et al. 2012; Ryon et al. 2017). We draw values of $r_c$ and $r_h$ from log-normal distributions centered on these values with standard deviations of 0.3 dex. Based on the observed correlation between $r_h$ and $r_c$ of young clusters, we also impose the requirement that $r_h \ge 2r_c$; if this condition is not met, we simply redraw $r_h$ until it is.

We adopt the density and velocity profiles of Stone & Ostriker (2015) for all clusters. The profiles were designed to resemble single-mass King models, while also being analytically simple. They are fully described by three parameters: the half-mass-radius $r_h$, the core radius $r_c$, and cluster mass $M_{cl}$. The density of stars near the center of the cluster for this profile is then given by $\rho_{*,c} = M_{cl}(r_c + r_h)/(2\pi^2 r_c^2 r_h^2)$. Figure 1 shows the distribution of stellar densities and velocity dispersions near the cluster centers for the adopted profiles and cluster sizes.

Two other relevant quantities are the average mass of a BH in the cluster, $\langle m_{\bullet} \rangle$, and the fraction of the total cluster mass locked in BHs, $M_{cl}/M_{cl}$. We calculate both quantities for a Kroupa (2001) stellar IMF, which has the form:

$$ \frac{dN}{dm_i} \propto \begin{cases} m_i^{-1.3} & \text{if } m_{\min} \le m_i \le 0.5 M_{\odot}, \\ m_i^{-2.3} & \text{if } 0.5 M_{\odot} < m_i < m_{\max}, \end{cases} \quad (1) $$

with lower and upper cutoff masses of $m_{\min} = 0.08 M_{\odot}$ and $m_{\max} = 150 M_{\odot}$. We have verified that our results are not sensitive to the particular choice of IMF. The IMF must be convolved with the relation between the mass of the progenitor star and compact remnant, $m_{rem}(m_i, Z)$, which depends on metallicity due to line-driven winds in the final stages of stellar evolution (Woosley et al. 2002; Belczynski et al. 2010a, 2010b; Spera et al. 2015; Vink 2017). For this, we adopt the results from the SEVN code, which also accounts for the gap in the BH mass distribution in low-metallicity environments due to pair instability supernovae (Spera & Mapelli 2017). We can then
<figure><img src="image_3.png" /><figcaption>Figure 1. Distribution of central stellar densities, $\rho_{\star,c}$, and central velocity dispersions, $\sigma_{\star}$, calculated for all clusters extracted from the CGL18 model. The distributions are weighted by the halo mass function, so as to be cosmologically representative. The height of each bin gives the number density of clusters in the bin (i.e., the total number density is given by the sum of bin values and not the integral under the histogram).</figcaption></figure>

calculate $M_{\star}$ and $\langle m_{\star} \rangle$ as

$$ M_{\star} = \int_{m_{i,\min,\star}}^{m_{\max,\star}} m_{\text{rem}}(m_i, Z) \frac{dN}{dm_i} dm_i \\ \langle m_{\star} \rangle = \frac{M_{\star}}{\int_{m_{i,\min,\star}}^{m_{\max,\star}} \frac{dN}{dm_i} dm_i} \qquad (2) $$

The bounds of both integrals are taken from the minimum initial stellar mass that will ultimately form a BH (set by the initial–final mass relation) to the upper mass cutoff of the IMF. In the following, for each cluster in the model, $\langle m_{\star} \rangle$ and $M_{\star}$ are

calculated based on its metallicity; as a reference, for our entire cluster population we find median values of $\langle m_{\star} \rangle \approx 20 M_{\odot}$ and $M_{\star}/M_{\text{cl}} \approx 5 \times 10^{-2}$.

## 2.2. Formation and Coalescence of BHBs

Below, we describe the time and spatial scales relevant to the formation of BHBs. Our prescription follows closely the approach of Antonini & Rasio (2016, hereafter AR16), who developed an analytic model for BHB evolution in star clusters. Unlike AR16, we do not perform a Monte Carlo draw of binaries in each cluster, but instead calculate typical quantities for the BHBs in each cluster, with the goal of understanding which types of clusters, on average, efficiently form mer- ging BHBs.

As the cluster evolves, the system seeks equipartition of kinetic energy, resulting in mass-segregation of the heavier BHs toward the cluster center on a timescale $\tau_{MS} \approx \tau_{rlx} \langle m_{\star} \rangle / \langle m_{\star} \rangle$, where $\tau_{rlx}$ is the cluster relaxation time at the half-mass-radius. Furthermore, if the cluster is susceptible to the Spitzer mass-stratification instability (Spitzer 1969; Watters et al. 2000), the BHs become confined to an ever-smaller core and form a decoupled self-gravitating sub-system. We find that all of our clusters are Spitzer–unstable.

In reality, binary formation in the core provides an energy source that halts core collapse. BHB formation can occur through interactions between three BHs. The rate of BHB formation through this process is $\dot{n}_{\text{BHB}} = 126G^5 \langle m_{\star} \rangle^5 n_{\star}^{3.3} \sigma_{\star}^{-9}$, where $n_{\star}$ and $\sigma_{\star}$ are the BH number density and velocity dispersion (Goodman & Hut 1993; Lee 1995). The velocity dispersion and densities of BHs in the dynamically decoupled core can be related to the corresponding quantities for the background stars. To do so, we begin by defining $\xi \gtrsim 1$ as the ratio of the mean kinetic energy of BHs to that of stars. This gives us the expression relating the stellar and BH velocity dispersions:

$$ \sigma_{\star} = \left( \xi \frac{\langle m_{\star} \rangle}{\langle m_{\star} \rangle} \right)^{1/2} \sigma_{\star} \qquad (3) $$

When the BHs dominate the mass in the core, the densities are related by (Lee 1995):

$$ \begin{aligned} \rho_{\star,c} &\approx \rho_h \frac{M_{\star} r_h^3}{M_{\text{cl}} r_{\star}^3} \\ &\approx \frac{1}{2} \left(\frac{r_c}{r_h}\right)^2 \left(\frac{M_{\star}}{M_{\text{cl}}}\right)^{-2} \rho_{\star,c} \end{aligned} \qquad (4) $$

where $\rho_h$ is the stellar density within $r_h$ and $r_{\star}$ is the half-mass–radius of BHs. In computing the above expression, we assumed $r_h \gg r_c$ and that the BHs are confined within a radius ($M_{\star}/M_{\text{cl}}$) smaller than the cluster half-mass–radius. In reality, the density of BHs in the core evolves with time as mass segregation proceeds. Because our model includes no time evolution, we have here taken the limit that mass segregation is complete, which matches well results from N-body cluster simulations.

Numerical simulations find that $\xi$ is typically a factor of a few (Gürkan et al. 2004). We adopt $\xi = 5$, which gives timescales for BHB formation consistent with the numerical study of Morscher et al. (2015). Using these relations, we write
the corresponding timescale for BHB formation as

$$ \tau_{3bb} = 5 \times 10^6 \left( \frac{\rho_{*,c}}{10^4 M_{\odot} \text{pc}^{-3}} \right)^{-2} \left( \frac{\langle m_{\bullet} \rangle}{20 M_{\odot}} \right)^{-3} \left( \frac{r_h/r_c}{4} \right)^4 \\ \times \left( 4\xi \frac{\langle m_{\star} \rangle}{\langle m_{\bullet} \rangle} \right)^{4.5} \left( \frac{\sigma_{\star}}{15 \text{ km s}^{-1}} \right)^{9} \left( \frac{M_{\star}/M_{\text{cl}}}{10^{-2}} \right)^4 \text{yr.} \quad (5) $$

BHBs can also form by exchanging into stellar binaries that are born at the time of cluster formation or formed dynamically during cluster evolution. An estimate of the timescale for the formation of BHBs through single-binary exchanges is given by (Miller & Lauburg 2009)

$$ \tau_{1-2} = 3.75 \times 10^9 \left(\frac{f_{\text{bin}}}{0.1}\right)^{-1} \left(\frac{\rho_{*,c}}{10^4 M_{\odot} \text{pc}^{-3}}\right)^{-1} \left(\frac{\langle m_{\star} \rangle}{1 M_{\odot}}\right) \\ \times \left(\frac{\sigma_{\star}}{15 \text{ km s}^{-1}}\right) \left(4\xi \frac{\langle m_{\star} \rangle}{\langle m_{\bullet} \rangle}\right)^{1/2} \left(\frac{m_{123}}{20 M_{\odot}}\right)^{-1} \left(\frac{a_{\text{hard}}}{1 \text{ au}}\right)^{-1} \text{yr.} \quad (6) $$

Here $f_{\text{bin}}$ is the stellar binary fraction, for which we adopt a value of 10%, and $m_{123}$ is the mass of the stellar binary plus the BH, which we approximate as $2\langle m_{\star} \rangle + \langle m_{\bullet} \rangle$. Finally, $a_{\text{hard}}$ is the typical separation of a hard stellar binary. We use as an estimate for $a_{\text{hard}}$ the maximum separation of a hard stellar binary in the cluster core (Heggie 1975; Quinlan & Shapiro 1989):

$$ a_{\text{hard}} = 1.5 \left( 3 \frac{q_{\star}}{1+q_{\star}} \right) \left( \frac{r_{\text{h}}}{3 \text{ pc}} \right) \left( \frac{\langle m_{\star} \rangle / M_{\text{cl}}}{10^{-5}} \right) \text{au}, \quad (7) $$

where $q_{\star} \le 1$ is the mass ratio of the stellar binary. We adopt the average value of $q_{\star} \approx 0.5$ for a Kroupa (2001) IMF. We define the binary formation timescale as $\tau_{\text{BHB}} = \min(\tau_{3bb}, \tau_{1-2})$.<fnref n="6" />

After BHBs form inside the cluster core, each BH-BHB interaction will cause them to both harden further and recoil. Denoting now the components of the BH binary as $m_1$ and $m_2$, and the mass of the BH interloper as $m_3$, the maximum separation of the binary below which interactions will lead to ejection from the cluster is (AR16)

$$ \begin{aligned} a_{\text{ej}} &= \frac{0.2G}{v_{\text{esc}}^2} \frac{q_{12}}{(1+q_{12})^2} \frac{q_3}{1+q_3} m_3 \\ &= 0.26 \left( \frac{v_{\text{esc}}}{30 \text{ km s}^{-1}} \right)^{-2} \left( \frac{m_3}{20 M_{\odot}} \right) \\ & \quad \times \left( 5 \frac{q_{12}}{(1+q_{12})^2} \right) \left( 3 \frac{q_3}{1+q_3} \right) \text{au}, \end{aligned} \quad (8) $$

where $q_{12} \le 1$ is the mass ratio of the BHB, $q_3 \equiv m_3/(m_1+m_2)$ and $v_{\text{esc}}$ is the escape velocity from the cluster center. For a Stone & Ostriker (2015) potential-density pair, this is

$$ v_{\text{esc}} = \frac{2}{\sqrt{\pi}} \sqrt{\frac{GM_{\text{cl}}}{r_{\text{h}} - r_{\text{c}}}} \sqrt{\ln(r_{\text{h}}/r_{\text{c}})}. \quad (9) $$

The separation $a_{\text{ej}}$ must be compared to the separation at which the BHB enters the regime in which gravitational waves dominate its evolution, $a_{\text{GW}}$. A reasonable estimate for $a_{\text{GW}}$ is

obtained from setting the rate of dynamically induced shrinking of the binary (Quinlan 1996):

$$ \dot{a}_{\text{dyn}} = -H \frac{G \rho}{\sigma^2} a^2, \quad (10) $$

where $H \approx 20$ is the binary hardening rate, equal to the rate at which the binary shrinks due to gravitational wave emission (Peters 1964), $\dot{a}_{\text{GW}}$. Using the core densities derived in Equation (4) yields the following expression for $a_{\text{GW}}$:

$$ \begin{aligned} a_{\text{GW}} &= 0.03 \left( \frac{m_1 + m_2}{40 M_{\odot}} \right)^{3/5} \left( 5 \frac{q_{12}}{(1 + q_{12})^2} \right)^{1/5} \\ & \quad \times \left( \frac{\sigma_{\star}}{15 \text{ km s}^{-1}} \right)^{1/5} \left( \frac{\rho_{*,c}}{10^4 M_{\odot} \text{pc}^{-3}} \right)^{-1/5} \\ & \quad \times \left( 4\xi \frac{\langle m_{\star} \rangle}{\langle m_{\bullet} \rangle} \right)^{1/10} \left( \frac{r_{\text{h}}/r_{\text{c}}}{4} \right)^{2/5} \left( \frac{M_{\star}/M_{\text{cl}}}{10^{-2}} \right)^{2/5} f^{1/5}(e) \text{au}; \\ f(e) &= \left( 1 + \frac{73}{24}e^2 + \frac{37}{96}e^4 \right) (1 - e^2)^{-7/2}, \end{aligned} \quad (11) $$

where $e$ is the eccentricity of the BHB, for which we adopt the average value $e = 2/3$ of a thermal distribution, $dP/de = 2e$ (Jeans 1919). For each cluster we use the average values of $q_{12}$ and $q_3$, which we calculate based on the metallicity of each cluster, for a Kroupa (2001) IMF convolved with the Spera & Mapelli (2017) initial-final mass relations.

If $a_{\text{GW}} < a_{\text{ej}}$, the BHB will be ejected before it reaches the point where GW emission can quickly coalesce the BHBs. In this case, the BHB is ejected from the cluster with a separation $a_{\text{ej}}$ and can continue to harden ex situ only via GW emission. On the other hand, if $a_{\text{GW}} > a_{\text{ej}}$, GW emission will cause the BHBs to coalesce in situ before dynamical ejection occurs.

Assuming that each interaction extracts 20% of the binary's binding energy (e.g., Quinlan 1996), the time to harden to a separation $a_{\text{crit}} = \max(a_{\text{GW}}, a_{\text{ej}})$, from an initial separation $a \gg a_{\text{crit}}$, can be written as (Miller & Hamilton 2002)<fnref n="7" />

$$ \begin{aligned} \tau_{\text{crit}} = 2.5 \times 10^8 & \left( \frac{q_3}{0.5} \right)^{-1} \left( \frac{\sigma_{\star}}{15 \text{ km s}^{-1}} \right) \left( 4\xi \frac{\langle m_{\star} \rangle}{\langle m_{\bullet} \rangle} \right)^{1/2} \\ & \times \left( \frac{a_{\text{crit}}}{0.1 \text{ au}} \right)^{-1} \left( \frac{m_1 + m_2}{40 M_{\odot}} \right)^{-1} \left( \frac{r_{\text{h}}/r_{\text{c}}}{4} \right)^2 \\ & \times \left( \frac{M_{\star}/M_{\text{cl}}}{10^{-2}} \right) \left( \frac{\langle m_{\bullet} \rangle}{20 M_{\odot}} \right) \left( \frac{\rho_{*,c}}{10^4 M_{\odot} \text{pc}^{-3}} \right)^{-1} \text{yr.} \end{aligned} \quad (12) $$

After reaching $a_{\text{crit}}$, the binary's hardening is dominated by the emission of GWs, which drives the system to coalescence on a timescale (Peters 1964)

$$ \begin{aligned} \tau_{\text{GW}} = 2.8 \times 10^8 & \left( \frac{m_1 m_2 (m_1 + m_2)}{10^4 M_{\odot}^3} \right)^{-1} \left( \frac{a_{\text{crit}}}{0.1 \text{ au}} \right)^4 \\ & \times \left( \frac{1 - e^2}{0.5} \right)^{7/2} \text{yr.} \end{aligned} \quad (13) $$

<footnote n="6">⁶ We neglect the contribution of two-body binary formation through gravitational bremsstrahlung. Lee (1995) showed that this channel is only relevant for $\sigma \gtrsim 100 \text{ km s}^{-1}$, but typical GCs have $\sigma \sim 10 \text{ km s}^{-1}$.</footnote>

<footnote n="7">⁷ AR16 label this timescale $t_{\text{merge}}$, but we find this label somewhat misleading and so opt for a different subscript.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Contours show the distribution of cluster disruption times as a function of the cluster mass $M_{cl}$, taking into account their galactocentric radii. The blue line shows the position-averaged disruption time used in CGL18, which reproduces the $z = 0$ cluster mass function.</figcaption></figure>

## 2.3. Timescales for Cluster Disruption

In this subsection, we introduce the timescales over which a cluster may be destroyed. In Section 3, these will be compared to the timescales for the formation and evolution of BHBs.

Two-body relaxation in star clusters leads to gradual “evaporation” of stars from the cluster and can eventually lead to the complete dynamical disruption of the cluster. To model this effect, we adopt the prescription of Gnedin et al. (2014) for mass loss in the presence of a strong external tidal field<fnref n="8" />:

$$ \frac{dM}{dt} = - \frac{M_{cl}(t)}{\tau_{tid}(M_{cl}, R)}, \qquad (14) $$

where R is the galactocentric radius of the GC and $\tau_{tid}$ is given by (Gieles & Baumgardt 2008)

$$ \tau_{tid}(M_{cl}, R) = 10^{10} \left( \frac{M_{cl}(t)}{2 \times 10^5 M_{\odot}} \right)^{2/3} P(R) \text{ yr} \qquad (15) $$

$$ P(R) = 41.4 \left( \frac{R}{\text{kpc}} \right) \left( \frac{V_{\text{circ}}}{\text{km s}^{-1}} \right)^{-1} \qquad (16) $$

Here $V_{circ}$ is the circular velocity in the galactic potential. We assume a flat rotation curve, with $V_{circ} = \sqrt{0.5GM_{bar}/R_h}$, and $M_{bar} = M_{\star} + M_{gas}$ is the total mass in baryons of the galaxy where the GC formed, which is also extracted from the CGL18 model. Integrating Equation (14) from the initial mass of the cluster, $M_{cl}$, to $M = 0$, gives the time to complete evaporation of the cluster:

$$ \tau_{\text{evap}} = 1.5 \times 10^{10} \left( \frac{M_{\text{cl}}}{2 \times 10^5 M_{\odot}} \right)^{2/3} P(R) \text{yr.} \qquad (17) $$

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Timescale for formation of BHBs via three-body and single-binary exchange, normalized by $\Delta t_{z=0} \equiv t(z=0) - t(z_{\text{form}})$ (upper panel) and the cluster disruption time $\tau_{\text{dis}} \equiv \min(\tau_{\text{evap}}, \tau_{\text{df}})$ (lower panel), as a function of the cluster mass, $M_{\text{cl}}$.</figcaption></figure>

CGL18 ignored the variation in spatial position of clusters, and instead adopted an average value of the factor $P(R)$ which allowed them to match the $z = 0$ GC mass functions. We find that our prescription for assigning cluster positions reproduces the average disruption timescale adopted by CGL18 well (Figure 2).

Finally, the cluster may inspiral into the center of the galaxy due to dynamical friction. We do not integrate the orbit of each cluster within the host galaxy’s potential, but instead use an approximate dynamical friction timescale as calculated by Binney & Tremaine (2008), evaluated at the initial position and mass of the cluster:

$$ \tau_{\text{df}} \approx 2.3 \times 10^{10} \left( \frac{M_{\text{cl}}}{2 \times 10^5 M_{\odot}} \right)^{-1} \left( \frac{R}{\text{kpc}} \right)^2 \times \left( \frac{V_{\text{circ}}}{100 \text{ km s}^{-1}} \right) \text{yr.} \qquad (18) $$

<footnote n="8">⁸ For simplicity, we ignore evaporation in the weak tidal field limit, which is almost always subdominant.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 4. Same as Figure 3, but now showing the timescale to harden the binary separation to $a_{crit} = \max(a_{GW}, a_{ej})$.</figcaption></figure>

After a time $\tau_{df}$ elapses we assume the cluster is tidally disrupted by, and subsequently merges with, the NSC. Finally, we define the timescale for the disruption of the cluster, $\tau_{dis} = \min(\tau_{evap}, \tau_{df})$.

The above prescriptions for the evolution of the GC within the host galaxy potential are necessarily simplifications. For example, GCs may migrate out of the disk and into the stellar halo after a galaxy merger (e.g., Kruijssen 2015; Li & Gnedin 2018) and gradual mass loss of the cluster during its inspiral may increase the total inspiral time (Gnedin et al. 2014). Nevertheless, we believe these prescriptions should provide a reasonable approximation on average.

## 3. Results

In this section, we begin by presenting the distributions of the various timescales discussed in Section 2. Then, we use a combination of these timescales to specify criteria for identifying the clusters most likely to produce merging BHBs.

<figure><img src="image_7.png" /><figcaption>Figure 5. Same as Figures 3 and 4, but now for the timescale to merge via GW emission, evaluated at $a_{crit}$. The change in the scaling of $\tau_{GW}$ with $M_{cl}$ at $M_{cl} \sim 10^6 M_{\odot}$ occurs because $a_{GW}$ begins to exceed $a_{crit}$ at this mass (see Figure 7).</figcaption></figure>

We then analyze the properties of these clusters relative to the overall cluster population. We conclude with an estimate of the cosmic BHB merger rate.

### 3.1. Timescales

The typical time for a BHB in a given cluster to merge is

$$ \tau_{tot} = \tau_{MS} + \tau_{BHB} + \tau_{crit} + \tau_{GW}. \qquad (19) $$

Figures 3–6 show the various timescales<fnref n="9" /> as a function of $M_{cl}$, normalized by both $\Delta t_{z=0} \equiv t(z=0) - t(z_{form})$ (i.e., the cosmic time between cluster formation and $z=0$) and the disruption timescale, $\tau_{dis}$.

For clusters with mass $M_{cl} \lesssim 10^{5.5} M_{\odot}$, the binary formation timescale is set by the three-body interaction timescale because

<footnote n="9">⁹ The mass-segregation timescale, $\tau_{MS}$, is always short compared to $t_H \lesssim 10^8$ yr, so we do not show it here.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Total time to BHB mergers for all clusters (see Equation (19)).</figcaption></figure>

of the low velocity dispersions in these systems (Figure 3). At higher masses ($M_{\text{cl}} \sim 10^{5.5} - 10^6 M_{\odot}$), three-body interactions and single-binary exchange both contribute, while at the very highest masses the timescale for single-binary exchange to pair BHs is always much shorter than the three-body timescale. Although single-binary exchange dominates at high cluster masses, the timescale for it to form BHBs is comparable to $\Delta t_z=0$. Thus, binary formation is difficult at high cluster masses.

Figure 4 shows the timescale $\tau_{\text{crit}}$ over which BHBs dynamically shrink in separation to $a_{\text{crit}}$. The median time spent in this phase is roughly constant at $\sim 10\%$ of $\Delta t_z=0$. After a time $\tau_{\text{crit}}$ has elapsed, BHBs spend a time $\tau_{\text{GW}}$ in the GW dominated regime, which is shown in Figure 5. For the lowest-mass clusters with $M_{\text{cl}} \sim 10^5 M_{\odot}$, we find that $\tau_{\text{GW}}$ is generally very large compared to $\Delta t_z=0$ because BHBs are ejected from the cluster with large separations, since these clusters have low escape speeds. More massive clusters, on the other hand, can retain and dynamically harden BHBs until they reach the regime where hardening via GW emission is effective, and thus

<figure><img src="image_9.png" /><figcaption>**Figure 7.** $a_{\text{ej}}$ and $a_{\text{GW}}$ as a function of $M_{\text{cl}}$. The shaded regions show interquartile ranges (25th to 75th percentile). To the left of the intersection, $a_{\text{ej}} > a_{\text{GW}}$ and BHBs merge ex situ, while to the right, $a_{\text{ej}} < a_{\text{GW}}$ and BHBs merge in situ. The transition mass between the two regimes is $M_{\text{cl}} \approx 10^6 M_{\odot}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Fraction of clusters, at fixed cluster mass, in which BHBs merge efficiently. This is calculated by dividing our cluster sample into narrow bins of cluster mass and calculating the fraction of clusters in each bin that satisfy the criteria put forth in Equations (20) and (21). The peak occurs near $M_{\text{cl}} \approx 10^{5.7} M_{\odot}$, around which 30% of clusters in that mass range produce merging BHBs.</figcaption></figure>

$\tau_{\text{GW}}$ is very short for typical BHBs in these clusters. The vastly different scalings of $\tau_{\text{GW}}$ with $M_{\text{cl}}$ are caused by the switch from evaluating $\tau_{\text{GW}} \propto a_{\text{crit}}^4$ at $a_{\text{ej}} \propto M_{\text{cl}}^{-1}$ to $a_{\text{GW}} \propto M_{\text{cl}}^{-1/10}$, which occurs at a typical mass of $M_{\text{cl}} \sim 10^6 M_{\odot}$ (Figure 7).

The scaling of $\tau_{\text{tot}}$ with cluster mass changes as different components of $\tau_{\text{tot}}$ begin to dominate (Figure 6). At low masses, $\tau_{\text{tot}}$ is typically greater than $\Delta t_z=0$ and dominated by $\tau_{\text{GW}}$. At intermediate and high masses ($M_{\text{cl}} \gtrsim 10^{5.5} M_{\odot}$), $\tau_{\text{tot}}$ is dominated by BHB formation.
<figure><img src="image_11.png" /><figcaption>Figure 9. Distributions of cluster masses, cosmic time at formation, and metallicities. Yellow and red show clusters where BHB mergers occur in situ and ex situ, respectively. Purple shows clusters that both disrupt before z = 0 and whose BHBs merge (either in or ex situ). The distributions are weighted by the halo mass function, so as to be cosmologically representative. The height of each bin gives the number density of clusters in the bin (i.e., the total number density is given by the sum of bin values and not the integral under the histogram). The total number density of clusters for “all,” “in-situ,” “ex-situ,” and “disrupted” are 2.6 Mpc⁻³, 0.10 Mpc⁻³, 0.31 Mpc⁻³, and 0.16 Mpc⁻³.</figcaption></figure>

The time $\tau_{tot}$ does not solely determine whether a cluster can form BHBs that will merge before $z = 0$. If $a_{GW} > a_{ej}$, the BHBs will merge in situ, and it is therefore necessary that the cluster not disrupt before BHB mergers can occur. Thus, we require

$$ 
\begin{align} 
\tau_{tot} &< \Delta t_{z=0} \\ 
\tau_{tot} &< \tau_{evap} \\ 
\tau_{tot} &< \tau_{df}. \tag{20}
\end{align} 
$$

On the other hand, if $a_{GW} < a_{ej}$, BHBs will merge ex situ. In this case, the destruction of the cluster will not affect the mergers of its BHBs, because they have already been ejected from the cluster. Thus, the conditions for a cluster to produce BHBs that merge before $z = 0$ are slightly relaxed compared to Equation (20):

$$ 
\begin{align} 
\tau_{tot} &< \Delta t_{z=0} \\ 
\tau_{tot} - \tau_{GW} &< \tau_{evap} \\ 
\tau_{tot} - \tau_{GW} &< \tau_{df}, \tag{21}
\end{align} 
$$

where $\tau_{tot}-\tau_{GW}$ represents the time to form and harden BHBs to a separation $a_{crit}$.

In both of the above sets of conditions, we implicitly ignore the contribution of clusters that disrupt before their BHBs have shrunk in separation to $a_{crit}$, even though hardening via GW emission will continue and may coalesce BHBs after the cluster disrupts and its constituents join the field or an NSC. For less massive clusters, where $a_{crit} = a_{ej}$, BHBs with separation $a \gtrsim a_{ej}$ would have GW coalescence times orders of magnitude greater than $\Delta t_{z=0}$ because $\tau_{GW} \propto a^4$. For more massive clusters, where $a_{crit} = a_{GW}$, in reality some BHBs may shrink to a separation $a \gtrsim a_{GW}$ which is small enough that GW emission would be sufficient to merge the BHBs after cluster disruption. Thus, we may slightly underestimate the number of clusters that merge BHBs in situ. However, based on the fraction of clusters that disrupt, our estimate can be off by no more than a factor of two.

## 3.2. Properties of Clusters That Make Merging BHBs

We find that $\approx 15\%$ of all clusters produce BHBs that can merge before $z = 0$. This fraction differs from some recent estimates derived by extrapolating results from numerical simulations, which find a much higher fraction of clusters producing merging BHBs. The difference is driven by the fact that our model spans a much wider range in initial cluster sizes than can be sampled in numerical simulations. In Section 3.3, we discuss in detail the implication of larger cluster sizes on the cosmic BHB merger rate.

In Figure 8 we divide the model clusters into bins of cluster mass and show the fraction of clusters in each bin that form merging BHBs. The curve peaks at a value of $\approx 30\%$ for $M_{cl} \approx 10^{5.7} M_{\odot}$, with a decreasing fraction at both lower and higher masses. Thus, clusters near a mass of $10^{5.7} M_{\odot}$ are the most effective, relative to the overall population at the same mass, at forming merging BHBs.

The median mass of all the clusters that form merging BHBs is $M_{cl} = 10^{5.5} M_{\odot}$, which is 0.25 dex higher than the median value for the entire cluster population. We can further divide into clusters that eject BHBs before coalescence and those that do not. The first panel in Figure 9 shows the mass distribution for the two cases, as well as the distribution for all clusters for comparison. In the case of ex situ BHB mergers, the number of clusters peaks and is about constant in the range $10^5 M_{\odot}-10^{5.5} M_{\odot}$, with a sharp truncation at higher masses. For the in situ case, the number of clusters peaks at $M_{cl} \approx 10^{5.7} M_{\odot}$. Approximately three times more clusters will have ex situ BHB mergers than in situ. This is largely due to the fact that in situ mergers occur in clusters that are a few times rarer, because of the bottom-heavy cluster IMF, $dN/dM_{cl} \propto M_{cl}^{-2}$. We note that if BHBs are on circular orbits, in contrast to our adopted eccentricity value $e = 2/3$, then the ratio of ex situ to in situ mergers will increase.

These distributions can be understood as follows. As previously discussed, for low-mass clusters ($M_{cl} \sim 10^5 M_{\odot}$), $a_{ej} \gg a_{GW}$ and BHBs are ejected from the cluster with large separations and correspondingly long $\tau_{GW}$. As a result, only a small fraction of low-mass clusters form BHBs that can merge,
leading to the nearly flat distribution of cluster masses for the case of ex situ mergers, even though low-mass clusters are the most numerous. At intermediate masses ($M_{cl} \sim 10^{5.3} M_{\odot}$), $a_{ej} \gtrsim a_{GW}$ and, although clusters typically eject BHBs before they reach $a_{GW}$, the ejection separation is small enough that the BHBs can coalesce before $z=0$ solely via GW emission after ejection. At higher masses ($M_{cl} \gtrsim 10^{5.7} M_{\odot}$), $a_{GW} \gtrsim a_{ej}$, and the cluster escape velocities are high enough that a large fraction of clusters can merge BHBs in situ. Finally, at the highest cluster masses, BHB formation timescales are comparable to $\Delta t_z=0$ (see Figure 3), so BHB formation and mergers are rare, thus leading to the slow decrease in the number of successful clusters out to the highest cluster masses.

The distribution of cosmic times of formation for clusters that form merging BHBs peaks at $t \approx 1.5$ Gyr, with a long tail extending toward later formation times. Unlike the distribution of formation times, the metallicity distribution is roughly flat between $-2 < \log_{10} Z/Z_{\odot} < 0$. The wide spread in metallicity reflects the fact that GC systems assemble from a diverse range of host galaxies which later merge to combine their respective GC populations.

The metallicities and formation times of clusters that form merging BHBs are similar to those of the overall cluster population. The ratio of the median metallicity of clusters whose BHBs merge to that of the overall population, $\log_{10} Z_{merge}/Z_{all}$, decreases with increasing $M_{cl}$ and varies between +0.3 dex at $M_{cl} \approx 10^5 M_{\odot}$ and -0.2 dex at $M_{cl} \approx 10^8 M_{\odot}$. Likewise, at $M_{cl} \approx 10^5 M_{\odot}$ clusters whose BHBs merge typically form 0.4 Gyr later than the overall population at that mass. The offset decreases with increasing $M_{cl}$ and vanishes for $M_{cl} \gtrsim 10^6 M_{\odot}$. The net result across all cluster masses is that the median metallicity for clusters that form merging BHBs is essentially identical to that of the overall population, and the median formation time is 0.1 Gyr later.

Clusters that evaporate or merge with the NSC before $z=0$ also contribute to the BHB merger rate, despite the fact that disruption limits the available time for BHB formation and hardening (see Figure 9). Approximately 40% of the clusters that produce merging BHBs will also disrupt before $z=0$. The majority of these disrupt due to dynamical evaporation, while only a small fraction merge with the NSC. For comparison, 60% of all clusters formed in the model are disrupted before $z=0$.

The distributions of properties of the clusters that both disrupt and form BHBs that eventually merge are represented by the purple curves in Figure 9. The similar heights of the curves for the disrupted clusters (purple) and the clusters whose BHBs merge ex situ (red) in the mass range $10^5 M_{\odot}-10^{5.5} M_{\odot}$ demonstrate that a large fraction of the hosts of ex situ merging BHBs will have been disrupted by $z=0$. In contrast, almost all of the hosts of in situ merging BHBs will still be around at the present day.

The distributions of core and half-mass-radii of clusters whose BHBs merge also follow log-normal distributions, but these distributions are shifted to smaller radii by $\approx 0.3$ dex relative to the adopted distributions for all clusters. The distribution of formation loci in the galaxy is indistinguishable from that of the whole population.

Table 1 summarizes the fraction of clusters satisfying various conditions.

<table><caption>**Table 1**
Summary of the Fraction of Clusters That Satisfy Various Conditions</caption><thead><tr><td>Description</td><td>Percent of All Clusters</td></tr></thead><tbody><tr><td>BHB mergers (in or ex situ)</td><td>15.8%</td></tr><tr><td>BHB mergers ex situ</td><td>11.8%</td></tr><tr><td>BHB mergers in situ</td><td>4.0%</td></tr><tr><td>BHB mergers (in or ex situ) + disrupt before z = 0</td><td>6.3%</td></tr></tbody></table>

**Note.** The middle two rows split all clusters that form merging BHBs into those in which the BHBs merge ex situ ($a_{ej} > a_{GW}$) and in situ ($a_{GW} > a_{ej}$). The final row refers to clusters that disrupt, either via dynamical evaporation or a merger with the NSC, but also produce BHBs that merge (regardless of where the BHB mergers occur).

### 3.3. Cosmic BHB Merger Rate

While the focus of this paper is to derive the properties of the clusters where BHBs are likely to form, using the framework outlined above for determining whether a cluster produces BHBs that merge by redshift $z$, we can provide a rough estimate of the cosmic merger rate of dynamically assembled BHBs as a function of redshift or cosmic time.

To estimate the merger rate, we select all clusters in the CGL18 model that satisfy the conditions in Equations (20) and (21). For each successful cluster we assume that one binary forms and coalesces after a time $\tau_{tot}$ has elapsed from the time of the cluster formation. Having the merger trees and masses of the halos where the clusters form gives us the number density of GCs as a function of cosmic time. We estimate the weight of each halo using its number density, $W_{HMF}$, from the halo mass function of the parent cosmological N-body simulation. With the information on the time when binaries merge we can directly calculate the rest-frame intrinsic merger rate by summing the number of BHB mergers per cosmic time:

$$ \frac{dn_{\text{com}}}{dt} = \sum_{k=1}^{N_{\text{tree}}} \frac{W_{\text{HMF}}}{N_{\text{h}}^{(k)}} \frac{N_{\text{mergers}}^{(k)}}{dt}, \quad (22) $$

where $N_{mergers}^{(k)}$ is the number of mergers within a time $dt$, $N_{h}^{(k)}$ is the number of $z=0$ halos in each bin, and $N_{tree}$ is the total number of merger trees available or, equivalently, the total number of $z=0$ halos (similar to the approach taken to estimate the merger rate for more massive BHs in the LISA band in Arun et al. 2009).

The number of events per unit observation time over the full sky is then given by

$$ \frac{dN}{dt} = \int 4\pi c \frac{dn_{\text{com}}(z)}{dz} \left[ \frac{D_L(z)}{(1+z)} \right]^2 dz, \quad (23) $$

where $n_{com}(z)$ is the comoving density of events at a given redshift (analogous to Equation (22) but in redshift rather than time bins) and $D_L$ is the luminosity distance (Haehnelt 1994).

To account for multiple BHBs forming in a given cluster, we note that most clusters which can form merging BHBs have $M_{cl} \sim 10^{5.5} M_{\odot}$ (see Figure 9). These clusters typically have $[\text{Fe/H}] \approx -0.7$, which corresponds to $M_{\odot}/M_{cl} \approx 3 \times 10^{-2}$ and $\langle m_{\odot} \rangle \approx 20 M_{\odot}$. Thus, a cluster of this mass will have $\approx 240$ BHBs, assuming all BHs are paired. While not all BHs become merging BHBs, we ignore in this calculation any cluster for which the average BHB does not merge before $z=0$ (following the various timescale conditions). In reality, these
clusters will also contribute to the merger rate. We therefore multiply our lower limit by 240 to obtain a more realistic estimate of the merger rate, and the results can be rescaled accordingly. However, we caution that this calculation represents only a rough estimate of the normalization of the merger rate; other trends such as the dependence on halo mass or redshift are more robust. We also do not apply a signal-to-noise threshold since we only use the average mass of BHs in each cluster.

Including this correction factor, we obtain at $z=0$ an intrinsic merger rate of $5.8 \text{ Gpc}^{-3} \text{ yr}^{-1}$, similar to those estimated in previous work. From run O1 of LIGO the rate of BHB mergers is in the range between 12 and $213 \text{ Gpc}^{-3} \text{ yr}^{-1}$ within the 90% credible region. Thus >10% of LIGO-Virgo events might have a dynamical origin. The total event rate, per unit of observation time, from Equation (23) is $\approx 1975 \text{ yr}^{-1}$.

Figure 10 shows our predicted intrinsic BHB merger rate as a function of cosmic time. After splitting into bins of halo mass and weighting by their respective number densities, we find that halos in the range $10^{12} M_{\odot} \lesssim M_{\text{h}}(z=0) \lesssim 10^{13} M_{\odot}$, similar to or slightly more massive than the Milky Way, have the largest contribution to the merger rate when averaging over a large, unbiased cosmic volume. The most massive halos, with $M_{\text{h}} \gtrsim 10^{14} M_{\odot}$, always contribute the least. Their low contribution is due to their rarity, not because their clusters have significantly different properties. Conversely, low-mass halos, $M_{\text{h}} \lesssim 10^{12} M_{\odot}$, contribute less than Milky Way-sized halos because they have fewer GCs, although the halos themselves are more frequent. However, because the number of successful clusters *per halo* scales linearly with the $z=0$ halo mass, in the case of a specific GW event which is localized to a patch of sky, the most massive halos would be the most likely hosts of the merger, if the binary was dynamically assembled in a star cluster.

### 3.4. Comparison to Other Merger Rate Estimates

The redshift evolution of the estimated merger rate can be compared to those predicted by other works (Rodriguez et al. 2016a; Fragione & Kocsis 2018; Rodriguez & Loeb 2018). We find a weakly increasing rate from $z=0$ out to $z=1.5$ ($t(z=1.5)=4.3$ Gyr), and a steep drop at higher redshift. Rodriguez et al. (2016a) calculated the merger rate evolution by assuming all clusters formed at $z=4$ and using the $z=0$ number density and mass function of GCs combined with the distribution of merger delay times from their suite of Monte Carlo cluster simulations. Fragione & Kocsis (2018) also used fits to the same set of numerical simulations and assumed all clusters formed at $z=3$, taking into account also the evolving number of GCs due to disruption using the GC evolution model of Gnedin et al. (2014) of a single Milky Way-type halo. To compute the total BHB merger rate, they extrapolated their results to all galaxies. Both works found a steady increase in the merger rate with increasing redshift out to their starting epoch. Finally, Rodriguez & Loeb (2018) used GC formation rates from the cosmological GC formation model of El-Badry et al. (2019) combined with delay times from the Rodriguez et al. (2016a) simulations. This study is most similar to ours, and finds a qualitatively similar redshift evolution but with a peak at higher redshift, $z \approx 3$. In Figure 11 we compare the various predicted merger rates.

A major difference between all three of these models and ours is the assumed cluster sizes. The suite of simulations run

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Predicted intrinsic BHB merger rate per comoving volume as a function of redshift. The black solid curve shows the total rate, while the colored dashed curves split the merger rate in bins of halo mass. We emphasize that the normalization of the merger rate here is only a rough approximation, assuming 240 BHBs merge per cluster satisfying the conditions in Equations (20) and (21), and that all other clusters do not contribute BHB mergers. The factor 240 comes from assuming that all clusters producing merging BHBs are of mass $M_{\text{cl}} = 10^{5.5} M_{\odot}$ (see Figure 9).</figcaption></figure>

by Rodriguez et al. (2016a), upon which all three models rely on calculating the BHB merger delay times, only covers clusters with sizes of 1 and 2 pc. Consequently, they find that BHB formation and hardening in the cluster occur quickly, and that the bottleneck is the time spent in the GW-dominated regime.

In contrast, because our model accounts for the full range of cluster sizes, we find that dynamical processes introduce in most clusters a non-negligible delay before BHBs reach the GW regime. Consequently, the typical delay time between cluster formation and BHB mergers in our fiducial model is a few Gyr, pushing the peak of the merger rate to lower redshift and lowering the overall normalization (Figure 12). Using fixed values of $r_{\text{h}} = 1$ pc and $r_{\text{c}} = 0.5$ pc for our entire cluster population, we find a nearly identical redshift evolution to that predicted by RL18. There remains a factor of $\approx 2$ difference in the normalization of the merger rate between our $r_{\text{h}} = 1$ pc model and the RL18 model at all redshifts; the discrepancy is likely due to a combination of effects relating to differences in the adopted GC formation model and treatment of BHB dynamics as well as the simplifying assumptions we made in calculating the merger rate. Given that the technique used to estimate the merger rate is approximate (since this is not the focus of this paper) we defer a more complete discussion of possible differences to future work aimed at this purpose (Choksi et al. 2018).

The redshift evolution of our model is also similar to those of Rodriguez et al. (2016a) and Fragione & Kocsis (2018) at $z \lesssim 3$, but differs at higher redshift. Both these models assume all GCs form in a single burst at their starting epoch; using an extended GC formation history removes the initial peak in the
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison of the merger rate with previous estimates. We report as a black solid curve the same curve of Figure 10, i.e., our full result in which cluster sizes are set from log-normal distributions peaked at $r_h = 2.8$ pc and $r_c = 1$ pc. The dashed line shows the model result when we fix $r_h = 1$ pc and $r_c = 0.5$ pc. Overplotted as points are the rates from the $r_{vir} = 1$ pc model of Rodriguez & Loeb (2018) and the fiducial model of Rodriguez et al. (2016a). The blue shaded region bounds the lower and upper limits of the rates calculated by Fragione & Kocsis (2018).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Distributions of the delay times to BHB mergers, from the time of cluster formation (i.e., $\tau_{tot}$ in our model). The distribution is skewed toward long delay times, resulting in the prediction of Figure 10 that the cosmic BHB merger rate of dynamically assembled BHBs peaks at relatively low redshift, $z \approx 1.5$. The integral under each histogram is normalized to unity.</figcaption></figure>

merger rate predicted by their models and causes a turnover, as seen in Figure 10. Finally, we note that at $z=0$ all recent predictions of the merger rate agree to within a factor of a few,

irrespective of the assumed cluster disruption rate or initial conditions. At higher redshift, the various predictions diverge by up to an order of magnitude.

4. Conclusions

In this work we analyzed the properties of the star clusters most likely to produce merging BHBs. Our model is based on realistic cluster populations extracted from a cosmological model of GC formation (CGL18). To each of these clusters, we applied analytic prescriptions for the formation, hardening, and merging timescales of BHBs, following the method outlined by AR16. This method allowed us to determine which clusters were likely to form merging BHBs and whether or not BHB mergers were likely to occur in the cluster. We extend previous work by modeling a diverse and realistic population of GCs, rather than a small subset of the parameter space. Our main conclusions are as follows.

1. We confirm previous calculations that massive star clusters can effectively create and shrink in separation BHBs. Around 15% of all clusters can produce BHBs that coalesce by $z=0$.

2. Of these, the number of clusters in which BHBs are ejected before merger is approximately three times the number of clusters in which BHBs merge in situ.

3. Ex situ BHB mergers originate from clusters of masses $10^5 M_{\odot}-10^{5.5} M_{\odot}$. In situ mergers occur in clusters over a much wider range of mass, with a peak at $10^{5.7} M_{\odot}$.

4. The distribution of formation times for clusters that form merging BHBs peaks at $t = 1-2$ Gyr, with a long tail extending toward later formation times; their metallicities span a wide range with a roughly flat distribution between $-2 < \log_{10} Z/Z_{\odot} < -0.5$. These distributions are similar to those of the overall cluster population.

5. 40% of the clusters that produce merging BHBs will also disrupt by $z=0$. Most clusters are disrupted due to two-body relaxation-driven evaporation, while a much smaller fraction inspiral into the NSC.

6. Using the mean properties for each cluster, and normalizing the number of BHBs per cluster at $M_{cl} = 10^{5.5} M_{\odot}$, the $z=0$ cosmic merger rate of dynamically assembled BHs is $\sim 6 \text{ Gpc}^{-3} \text{yr}^{-1}$.

7. The merger rate of dynamically assembled BHBs is weakly increasing out to $z \sim 1.5$ and drops at higher redshift. This behavior is driven by dynamical processes within the cluster, which introduce a significant delay between cluster formation and BHB mergers.

We thank Carl Rodriguez and the anonymous referee for comments and discussions that improved this work, as well as Mario Spera and Michaela Mapelli for making the SEVN code publicly available. N.C. acknowledges support from the Centre for Cosmological studies and thanks Joe Silk and Peter Behroozi for helping make this project possible as well as Goni Halevi for support throughout the preparation of this work. M.V. acknowledges funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 grant Agreement no. 614199, project "BLACK"). O.G. was supported in part by the NSF through grant 1412144.
ORCID iDs

Nick Choksi https://orcid.org/0000-0003-0690-1056

Marta Volonteri https://orcid.org/0000-0002-3216-1322

Monica Colpi https://orcid.org/0000-0002-3370-6152

Oleg Y. Gnedin https://orcid.org/0000-0001-9852-9954

Hui Li https://orcid.org/0000-0002-1253-2763

References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, PhRvL, 116, 061102

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, PhRvL, 116, 241103

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, PhRvL, 118, 221101

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, PhRvL, 119, 141101

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, ApJL, 851, L35

Antonini, F. 2013, ApJ, 763, 62

Antonini, F., Capuzzo-Dolcetta, R., Mastrobuono-Battisti, A., & Merritt, D. 2012, ApJ, 750, 111

Antonini, F., & Rasio, F. A. 2016, ApJ, 831, 187

Arun, K. G., Babak, S., Berti, E., et al. 2009, CQGra, 26, 094027

Bae, Y.-B., Kim, C., & Lee, H. M. 2014, MNRAS, 440, 2714

Bastian, N. 2016, EAS, 80, 5

Bastian, N., Adamo, A., Gieles, M., et al. 2012, MNRAS, 419, 2606

Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010a, ApJ, 714, 1217

Belczynski, K., Dominik, M., Bulik, T., et al. 2010b, ApJL, 715, L138

Belczynski, K., Holz, D. E., Bulik, T., & O'Shaughnessy, R. 2016, Natur, 534, 512

Benacquista, M. J., & Downing, J. M. B. 2013, LRR, 16, 4

Binney, J., & Tremaine, S. 2008, Galactic Dynamics (2nd ed.; Princeton, NJ: Princeton Univ. Press)

Breen, P. G., & Heggie, D. C. 2013a, MNRAS, 432, 2779

Breen, P. G., & Heggie, D. C. 2013b, MNRAS, 436, 584

Capuzzo-Dolcetta, R., & Miocchi, P. 2008a, MNRAS, 388, L69

Capuzzo-Dolcetta, R., & Miocchi, P. 2008b, ApJ, 681, 1136

Chatterjee, S., Fregeau, J. M., Umbreit, S., & Rasio, F. A. 2010, ApJ, 719, 915

Chatterjee, S., Rodriguez, C. L., Kalogera, V., & Rasio, F. A. 2017a, ApJL, 836, L26

Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2017b, ApJ, 834, 68

Chatterjee, S., Umbreit, S., Fregeau, J. M., & Rasio, F. A. 2013, MNRAS, 429, 2881

Choksi, N., & Gnedin, O. Y. 2018, arXiv:1810.01888

Choksi, N., Gnedin, O. Y., & Li, H. 2018, MNRAS, 480, 2343

de Mink, S. E., & Mandel, I. 2016, MNRAS, 460, 3545

Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52

El-Badry, K., Quataert, E., Weisz, D. R., Choksi, N., & Boylan-Kolchin, M. 2019, MNRAS, 482, 4528

Fall, S. M., & Zhang, Q. 2001, ApJ, 561, 751

Förster Schreiber, N. M., Renzini, A., Mancini, C., et al. 2018, arXiv:1802.07276

Fragione, G., & Kocsis, B. 2018, PhRvL, 121, 161103

Gieles, M., & Baumgardt, H. 2008, MNRAS, 389, L28

Gnedin, O. Y., Lee, H. M., & Ostriker, J. P. 1999, ApJ, 522, 935

Gnedin, O. Y., Ostriker, J. P., & Tremaine, S. 2014, ApJ, 785, 71

Goodman, J., & Hut, P. 1993, ApJ, 403, 271

Gürkan, M. A., Freitag, M., & Rasio, F. A. 2004, ApJ, 604, 632

Haehnelt, M. G. 1994, MNRAS, 269, 199

Heggie, D. C. 1975, MNRAS, 173, 729

Hong, J., Vesperini, E., Askar, A., et al. 2018, MNRAS, 480, 5645

Ivanova, N., Justham, S., Chen, X., et al. 2013, A&ARv, 21, 59

Jeans, J. H. 1919, MNRAS, 79, 408

Kravtsov, A. V., & Gnedin, O. Y. 2005, ApJ, 623, 650

Kroupa, P. 2001, MNRAS, 322, 231

Kruijssen, J. M. D. 2015, MNRAS, 454, 1658

Lee, H. M. 1995, MNRAS, 272, 605

Li, H., & Gnedin, O. Y. 2018, arXiv:1810.11036

Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140

Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634

Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50

Miller, M. C., & Hamilton, D. P. 2002, MNRAS, 330, 232

Miller, M. C., & Lauburg, V. M. 2009, ApJ, 692, 917

Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319

Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 2015, ApJ, 800, 9

Muratov, A. L., & Gnedin, O. Y. 2010, ApJ, 718, 1266

Nelson, D., Pillepich, A., Genel, S., et al. 2015, A&C, 13, 12

O'Leary, R. M., O'Shaughnessy, R., & Rasio, F. A. 2007, PhRvD, 76, 061504

Park, D., Kim, C., Lee, H. M., Bae, Y.-B., & Belczynski, K. 2017, MNRAS, 469, 4665

Peters, P. C. 1964, PhRv, 136, 1224

Pfeffer, J., Kruijssen, J. M. D., Crain, R. A., & Bastian, N. 2018, MNRAS, 475, 4309

Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJL, 528, L17

Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431

Postnov, K. A., & Yungelson, L. R. 2014, LRR, 17, 3

Quinlan, G. D. 1996, NewA, 1, 35

Quinlan, G. D., & Shapiro, S. L. 1989, ApJ, 343, 725

Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016a, PhRvD, 93, 084029

Rodriguez, C. L., Haster, C.-J., Chatterjee, S., Kalogera, V., & Rasio, F. A. 2016b, ApJL, 824, L8

Rodriguez, C. L., & Loeb, A. 2018, ApJL, 866, L5

Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, PhRvL, 115, 051101

Rodríguez-Puebla, A., Behroozi, P., Primack, J., et al. 2016, MNRAS, 462, 893

Ryon, J. E., Gallagher, J. S., Smith, L. J., et al. 2017, ApJ, 841, 92

Spera, M., & Mapelli, M. 2017, MNRAS, 470, 4739

Spera, M., Mapelli, M., & Bressan, A. 2015, MNRAS, 451, 4086

Spitzer, L., Jr. 1969, ApJL, 158, L139

Stone, N. C., & Ostriker, J. P. 2015, ApJL, 806, L28

The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, B. P., et al. 2018, arXiv:1811.12907

Tremaine, S. D., Ostriker, J. P., & Spitzer, L., Jr. 1975, ApJ, 196, 407

Vink, J. S. 2017, RSPTA, 375, 20160269

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Natur, 509, 177

Watters, W. A., Joshi, K. J., & Rasio, F. A. 2000, ApJ, 539, 331

Woosley, S. E., Heger, A., & Weaver, T. A. 2002, RvMP, 74, 1015