DYNAMICAL FORMATION SIGNATURES OF BLACK HOLE BINARIES IN THE FIRST DETECTED MERGERS BY LIGO

After over two decades of development, the Advanced Laser Interferometer Gravitational-Wave Observatory³ (aLIGO) has directly detected gravitational waves from an inspiraling black hole–black hole (BH–BH) binary (Abbott et al. 2016b). Through the precise measurement of the gravitational waves, aLIGO is capable of characterizing many properties of inspiraling binaries, including the total mass of the binary, M_tot, the mass ratio, q, and the black holes' spins. The first detected BH–BH binary, GW150914, had unusually high component masses $({m}_{1},{m}_{2})=({36}_{-4}^{+5}{M}_{\odot },{29}_{-4}^{+4}{M}_{\odot })$ in comparison to BH masses inferred for star–BH X-ray binaries (Farr et al. 2011; Fragos & McClintock 2015). A second, less significant event⁴ , LVT151012, also had high inferred masses $({m}_{1},{m}_{2})=({23}_{-5}^{+18}{M}_{\odot },{13}_{-5}^{+4}{M}_{\odot })$ (The LIGO Scientific Collaboration & the Virgo Collaboration 2016). aLIGO has finally opened a new window to our universe. Along with other upcoming instruments VIRGO⁵ and KAGRA⁶ , aLIGO will allow us to probe the demographics of potentially hundreds of BH–BH binaries (Abbott et al. 2016a).

There are three primary pathways that lead to BH–BH binaries that can merge within the age of the universe, through binary evolution, through gas dynamics (Davies et al. 2011; Bartos et al. 2016; Stone et al. 2016), and through stellar dynamics (see Benacquista & Downing 2013 for a review). First, such binaries can form through the evolution of isolated, massive binary stars.⁷ A major bottleneck in our understanding of this channel is the complex tidal (de Mink & Mandel 2016; Mandel & de Mink 2016; Marchant et al. 2016) and common envelope (e.g., Belczynski et al. 2016) evolution such a binary must go through in order to produce two BHs that can merge within a Hubble time from the loss of gravitational waves. This is in addition to uncertainties in the details of massive star evolution, supernova explosions, and the birth kicks of black holes. Nevertheless, sophisticated population synthesis routines have been developed that incorporate many of these uncertainties to make predictions about the properties of the first gravitational-wave sources (e.g., Belczynski et al. 2008, 2010b; Dominik et al. 2012; Mennekens & Vanbeveren 2014). A second possibility to get BH binary mergers is through gas-assisted mergers (Escala et al. 2005; Davies et al. 2011); however, simple rate estimates suggest that gas-assisted stellar BH mergers are probably relatively uncommon (Bartos et al. 2016; Stone et al. 2016).

Dynamical interactions of BHs in dense stellar environments, such as globular clusters, present another promising method to produce tight BH–BH binaries, whether through exchange interactions (Portegies Zwart & McMillan 2000), three-body dynamics (Miller & Hamilton 2002a; Wen 2003; Antognini et al. 2014; Antonini et al. 2014), or direct dynamical capture⁸ (Quinlan & Shapiro 1989; Lee 1993; O'Leary et al. 2009). In these scenarios, the black holes that form at the end stage of stellar evolution collect near the center of the cluster through dynamical friction (Kulkarni et al. 1993; Sigurdsson & Hernquist 1993). Because of gravitational focusing, the most massive BHs are preferentially involved in subsequent gravitational encounters and form BH–BH binaries. These binaries may merge within the cluster (Miller & Hamilton 2002b) or be ejected from the cluster and merge on much longer timescales (O'Leary et al. 2006). Such models have presented their own theoretical obstacles, the initial mass function of BHs perhaps the largest, but ever sophisticated simulations over nearly a decade have generally found similar estimates for the expected merger rate of the binaries as well as their characteristics (Gültekin et al. 2006; Moody & Sigurdsson 2009; Banerjee et al. 2010; Downing et al. 2010, 2011; Morscher et al. 2013; Bae et al. 2014; Ziosi et al. 2014; Morscher et al. 2015; Rodriguez et al. 2015). These results remain even in simulations that have shown a substantial fraction of BHs remain in the cluster (e.g., Mackey et al. 2008; Morscher et al. 2015; Rodriguez et al. 2015)

In this work, we present a number of observational signatures of the dynamical formation of black hole binaries. In particular, we focus on signatures that are independent of the poorly known black hole initial mass function (IMF; Belczynski et al. 2004; Özel et al. 2010; Farr et al. 2011; Belczynski et al. 2016). A number of studies have qualitatively discussed that dynamical interactions preferentially form binaries with the most massive components in the cluster (e.g., O'Leary et al. 2006, 2007; Morscher et al. 2015). Yet few studies focused on the expected mass ratio distribution of the BH–BH binaries that merge. In this work, we use the original Monte Carlo results of O'Leary et al. (2006), a series of new Monte Carlo simulations, as well as a series of new direct N-body simulations to explore the expected mass distribution of the components of BH–BH binaries. We argue that the mass distribution of the BH binaries found by aLIGO will present a unique signature of dynamically formed binaries and their underlying mass function. After we submitted our manuscript, two papers appeared on the expected rates of stellar black hole mergers in globular clusters with independent methodologies, which confirm our findings (Chatterjee et al. 2016; Rodriguez et al. 2016)

After the first supernovae, the more massive BHs collect near the center of the cluster owing to dynamical friction from the low mass stellar background. In contrast to previous expectations (Kulkarni et al. 1993; Sigurdsson & Hernquist 1993), however, these BHs do not interact exclusively among themselves. After the BH subcluster experiences core collapse and produces a binary, that binary acts through the BH population to heat the entire star cluster until it is ejected from the system (Mackey et al. 2008; cf. O'Leary et al. 2014). Nevertheless, the strongest dynamical interactions that eventually lead to binary ejection rarely include stars (Morscher et al. 2015). In this work, we are interested in the properties of the BH–BH binaries, and only follow the evolution of the BHs separately from the cluster. Although this does not follow the proper evolution of the entire BH cluster, it does capture the important dynamical interactions that lead to binary formation and, eventually, binary heating. In this work, we use two methods for following the evolution of the BHs in dense stellar clusters. We use the Monte Carlo method presented in O'Leary et al. (2006) as well as a direct N-body method (Wang et al. 2015) that follows the evolution of only the BHs. Both methods focus on the evolution of the BHs in isolation from the stellar cluster.

Our Monte Carlo simulations of BHs in globular clusters follow directly from O'Leary et al. (2006) and are based on the original method by Ivanova et al. (2005, 2010). Globular clusters are treated as two-component systems, with a constant density core, where all dynamical interactions take place, and a low-density halo. BHs that are kicked into the halo from dynamical interactions remain there until dynamical friction returns the BHs to the cluster center. The three- and four-body interactions in the cluster are sampled from their likelihood distributions and directly integrated with fewbody (Fregeau et al. 2004). Binaries that are present in the simulation are evolved following Peters (1964) equations that describe the inspiral due to gravitational-wave (GW) emission. If two nonspinning BHs merge within the cluster their merger product receives a kick due to the asymmetric emission of GWs with a maximum velocity of 175 km s⁻¹ near q ≡ M₁/M₂ ≈ 1/5 (González et al. 2007). Equal mass mergers receive no kick as they are symmetric. The BH mass function changes in the simulation due to BH mergers and ejections.

We also run a series of direct N-body simulations to follow the evolution of the BHs that form in a much larger stellar cluster. We integrate the orbits of all the BHs directly using nbody6++ (Wang et al. 2015) and keep track of the dynamically formed binary population. Future work will include the full post-Newtonian treatment throughout the cluster evolution, as mergers in the cluster likely constitute about ∼15% of all mergers (Morscher et al. 2015). These simulations, due to their initial conditions, produce too few BH–BH binaries that merge in a Hubble time, so we include all binaries ejected from the system in our analysis to compare with the Monte Carlo results. Two model suites were produced. First, the model that includes a potential that mimics the underlying stellar cluster using a static Plummer sphere with half-mass radius R_h = 13 pc and mass M_cl = 2 × 10⁵ ${M}_{\odot }$ and a BH subcluster made of 492 single particles and 10 "primordial" binary pairs with half-mass radius of 3.3 pc. We also simulate the cluster of BHs without an external potential. This model has only 512 single particles taken from a Plummer distribution with half-mass radius of 0.5 pc. The integration is performed to 10 Gyr in all models.

In both Monte Carlo and N-body simulations, we generate our initial BH mass function from a power-law distribution

$$\begin{eqnarray}&&f(M)\;\equiv \;\displaystyle \frac{{dN}}{{dM}}\;\propto {M}^{-\alpha },\end{eqnarray} \tag{ 1 }$$

between minimum and maximum BH mass of M_min = 5 ${M}_{\odot }$ and M_max = 40 ${M}_{\odot }$ with α = 0, 1, 2, 3, and 4. While this range of masses is larger than the inferred masses of BHs in X-ray binaries (Özel et al. 2010; Farr et al. 2011), it is consistent with models of BHs that form with fallback accretion and direct collapse (Belczynski et al. 2004, 2010a) and is minimally required to explain GW150914. We also run a set of simulations using the BH mass function found in the population synthesis studies of Belczynski et al. (2004) for low-metallicity stars. This mass function has three peaks in its distribution, near 7 ${M}_{\odot }$, 14 ${M}_{\odot }$, and 24 ${M}_{\odot }$, and represents a more complicated and perhaps more realistic distribution of masses. In the Monte Carlo simulations, we assume that 10% of black holes are in primordial binaries (cf. the initial conditions for the N-body above), with a period distribution that is flat in the log of period.

These two methods are complementary to each other, and both capture the strongest dynamical encounters between BHs. Nevertheless they have their limitations and benefits. In both series of simulations, we do not directly follow the entire coupling between the main stellar population in the cluster and the black holes, which is important to determine the present-day population of BHs in globular clusters (e.g., Mackey et al. 2008; Morscher et al. 2015; Rodriguez et al. 2016; Wang et al. 2016). However, the latter simulations confirm that the dynamical interactions that lead to binary ejections are driven exclusively by the BH subcluster and are dominated by BH–BH binaries. With the fast Monte Carlo simulations, we can run a much larger number of simulations with different compositions and initial conditions than studies that follow the evolution of the entire stellar population, and we can turn important processes, such as three-body binary formation, off in order to understand which aspects of the dynamical interactions are most important.

To leading (2.5 post-Newtonian) order, the root mean square (rms) detection signal to noise ratio for a circular inspiraling binary is given as

$$\begin{eqnarray}&&\displaystyle \frac{S}{N}=k\displaystyle \frac{{\eta }^{1/2}{M}_{{\rm{tot}}}^{5/6}}{{d}_{{\rm{L}}}}{I}_{7}^{1/2}({M}_{{\rm{tot}}})\end{eqnarray} \tag{ 2 }$$

where $k={\pi }^{-2/3}\sqrt{2/15}\approx 0.17$ for an isotropic binary orientation relative to the detector, η = q/(1 + q)² is the symmetric mass ratio, d_L is the luminosity distance, and

$$\begin{eqnarray}&&{I}_{7}(M)={\int }_{{f}_{{\rm{min}}}}^{{f}_{{\rm{ISCO}}}(M)}\displaystyle \frac{{f}^{-7/3}}{{S}_{h}(f)}\;{df}\end{eqnarray} \tag{ 3 }$$

where ${f}_{{\rm{ISCO}}}({M}_{{\rm{tot}}})={\pi }^{-1}{6}^{-3/2}{M}^{-1}$ is the GW frequency at the innermost stable circular orbit, f_min is the lowest GW frequency that the instrument can detect, and S_h is the one-sided noise spectral density (Cutler & Flanagan 1994; Dalal et al. 2006).⁹ We set f_min = 10 Hz and S_h(f) using the calibrated sensitivity spectra of aLIGO-Hanford on 2015 October 1 (Kissel 2015). For a fixed detection threshold (e.g., S/N ≥ 7), the maximum distance range d_L,max may be obtained from Equation (2). Assuming that the source population is uniformly distributed in volume, the detection rate is biased by approximately

$$\begin{eqnarray}&&{V}_{{\rm{det}}}\propto {d}_{{\rm{L}},\mathrm{max}}^{3}\propto \displaystyle \frac{{q}^{3/2}}{{(1+q)}^{3}}{M}_{{\rm{tot}}}^{5/2}{I}_{7}^{3/2}({M}_{{\rm{tot}}}).\end{eqnarray} \tag{ 4 }$$

This function grows as ${M}_{{\rm{tot}}}^{2.3}$ between M ∼ 10 and 20 ${M}_{\odot }$, has a maximum at M_tot = 77 ${M}_{\odot }$, and decreases to zero at 439 ${M}_{\odot }$, where f_ISCO = f_min = 10 Hz. Compared to 10 ${M}_{\odot }$ and fixed q, the detectable volume is a factor of 27 larger for M_tot = 77 ${M}_{\odot }$ and a factor of 4.5 smaller at 400 ${M}_{\odot }$.

Dynamical interactions between BHs in stellar clusters are primarily driven by mass segregation, gravitational focusing, and multibody dynamics. We therefore expect the most massive BHs to be preferentially found in merging BH–BH binaries.

In Figure 1, we compare the mass distribution of merged BH–BH binaries from our simulations to that expected from random pairings of BHs from the BH mass function. We find that dynamical interactions enhance the merger rate of BH binaries relative to random pairings in a star cluster by a boost factor that scales as a power law,

$$\begin{eqnarray}&&B={k}_{1}{M}_{{\rm{tot}}}^{\beta }\end{eqnarray} \tag{ 5 }$$

with β ≳ 4, and we assume that k₁ is independent of (q, M_tot). All of our Monte Carlo simulations have β ≈ 4. This is even true for the Belczynski et al. (2004) mass function, which is not monotonically decreasing as a function of mass and has multiple peaks. We have also found similar results for our Monte Carlo simulations with a variety of cluster models, with differing velocity dispersions and escape velocities. Although our N-body simulations with an external potential still had β ≈ 4, we found that the binaries that form in the cluster without the external potential have β ≫ 4. We discuss this in more detail in Section 4.

Zoom In Zoom Out Reset image size

The total event rate per star cluster is then

$$\begin{eqnarray}&&{\rm{\Gamma }}={\int }_{{M}_{\mathrm{min}}}^{{M}_{\mathrm{max}}}{\int }_{{M}_{\mathrm{min}}}^{{M}_{\mathrm{max}}}{Bf}({M}_{1})f({M}_{2})\;{{dM}}_{1}\;{{dM}}_{2}\end{eqnarray} \tag{ 6 }$$

which depends on the total binary mass as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Gamma }}}{\partial {M}_{{\rm{tot}}}}\propto {M}_{{\rm{tot}}}^{\beta }{\int }_{\mathrm{max}({M}_{{\rm{tot}}}-{M}_{\mathrm{max}},{M}_{\mathrm{min}})}^{\mathrm{min}({M}_{{\rm{tot}}}-{M}_{\mathrm{min}},{M}_{\mathrm{max}})}f({M}_{1})f({M}_{{\rm{tot}}}-{M}_{1}){{dM}}_{1}\end{eqnarray} \tag{ 7 }$$

This integral may be evaluated analytically for f(M) ∝ M^−α using Gauss's hypergeometric function ${}_{2}{F}_{1}(a,b;c;z)$ as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Gamma }}}{\partial {M}_{{\rm{tot}}}}\propto {M}_{{\rm{tot}}}^{1+\beta -2\alpha }{\left.z\ {}_{2}{F}_{1}\left(\displaystyle \frac{1}{2},\alpha ;\displaystyle \frac{3}{2};{z}^{2}\right)\right|}_{{z}_{\mathrm{min}}}^{{z}_{\mathrm{max}}}\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&{z}_{\mathrm{min}}=\mathrm{max}\left(1-\displaystyle \frac{2{M}_{\mathrm{max}}}{{M}_{{\rm{tot}}}},\displaystyle \frac{2{M}_{\mathrm{min}}}{{M}_{{\rm{tot}}{\rm{}}}}-1\right)\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{z}_{\mathrm{max}}=\mathrm{min}\left(1-\displaystyle \frac{2{M}_{\mathrm{min}}}{{M}_{{\rm{tot}}}},\displaystyle \frac{2{M}_{\mathrm{max}}}{{M}_{{\rm{tot}}{\rm{}}}}-1\right).\end{eqnarray} \tag{ 10 }$$

This may be expressed with elementary functions for integer and half-integer α. Plots of Equation (8) show that the merger rate peaks at ${M}_{{\rm{tot}}}\sim {M}_{\mathrm{max}}$ for 1 ≲ α ≲ 3.5 and near 2M_min for larger α. Accounting for the correlation between M_tot and q seen in Figure 2 skews ∂Γ/∂M_tot to even higher M_tot, discussed next. We discuss how this relationship may be used to probe the underlying mass distribution in Section 4.

Zoom In Zoom Out Reset image size

Dynamical encounters not only preferably form high-mass BH–BH binaries, but as M_tot approaches 2M_max the binaries necessarily approach q = 1. In Figure 2, we show two-dimensional (M_tot–q) distribution functions for BH–BH mergers from the cluster for two BH mass functions using our Monte Carlo results. These distributions were recreated from the simulation results using a Gaussian kernel density estimator to generate a smooth distribution. The top panel shows the results from our power-law mass function M⁻² and M_max = 40 ${M}_{\odot }$. The bottom panel shows the results from the Belczynski et al. (2004) mass function. Both distributions peak near q = 1.0. The peak is much more pronounced in mass functions with α < 2. For mass functions with α ≳ 2, the q distribution marginalized over M_tot rises rapidly near q = 0.4 and is flat for higher q until q = 1.0.

Since the maximum distance to which aLIGO may detect a binary is different for different M_tot and q, the observed rate distribution of dynamically formed binaries will be biased relative to the merger rate, ${ \mathcal R }={V}_{{\rm{det}}}{\rm{\Gamma }}$ (see Equation (4)). O'Leary et al. (2009) have shown that for the planned aLIGO spectral noise density, the detection distance will be largest for M_tot between 60 and 90 ${M}_{\odot }$ for circular binaries (see their Figure 11). For the most recent observing run, O1, aLIGO was most sensitive to BH mergers with M_tot ≈ 77 ${M}_{\odot }$. Further, the detection rate is biased toward equal-mass binaries as shown by Equation (4).

The overall impact of aLIGO's sensitivity can be seen by comparing the left and right panels of Figure 2. In the right panel, we show the distribution of mergers weighted by the observable volume of the universe, respectively, for different (M_tot,q) for the aLIGO O1 observing run. All results have been normalized to peak at 1. By comparing our results to the panels on the left, we can see that aLIGO should expect to see high-mass black hole binaries, with mass ratios peaking near 1. However, in these models the overall distribution is mostly determined by the underlying dynamics of BHs in star clusters, rather then aLIGO's intrinsic sensitivity. This suggests that Earth-based GW detectors will be a useful probe of the underlying demographics of BHs.

So far, we have only compared the binary merger distribution to the BH IMF without regard for the age of the cluster. While young clusters may be an important contributor to the merger rate (O'Leary et al. 2007), it is not clear if they survive long enough to dominate the rate. Most studies therefore focus on an old population of clusters, such as globulars, which are long lived. In Figure 3, we show the merger and detection rate distributions from binaries during the last 2.5 Gyr, roughly the observable horizon of aLIGO. We find that mergers from old clusters is better fit by a broken power law. At high masses M_tot ≳ M_max, the merger rate is still enhanced by β ≳ 4 as in Figure 1. However, for M_tot ≲ M_max, β ≈ 0.

Zoom In Zoom Out Reset image size

In dense stellar systems, the black holes that form at the end stage of stellar evolution collect near the center of the cluster through dynamical friction. Gravitational focusing in subsequent multibody dynamical interactions naturally increase the interaction rate between the most massive black holes in the cluster, eventually leading them to merge. In this work, we have used a suite of numerical simulations to better understand the properties of BH–BH binaries that are formed in these clusters. We have found that dynamical interactions between the BHs strongly enhance the rate of mergers among higher-mass BHs, such that the merger rate of binaries is boosted by a factor $\propto {M}_{{\rm{tot}}}^{\beta }$, with β ≳ 4, over a random selection of BH pairs from the cluster.

This relationship may be a useful tool to probe the underlying mass distribution of BHs in the aLIGO era with multiple detections. The merger distribution can be estimated directly from the underlying BH mass function, assuming that aLIGO can distinguish the two primary channels of producing BH–BH binaries, or that BH–BH binaries formed during dynamical interactions dominate the detection rate. Under this assumption, the aLIGO measurement statistics will give physical meaning to the underlying mass function and dynamics.

In this work, we have minimally explored the underlying reason for β ≳ 4. In the dynamical interactions considered here, the process most sensitive to M_tot is the rate of binary formation through three-body encounters, which scales as ${M}_{{\rm{tot}}}^{5}$ (Ivanova et al. 2005). We have run our Monte Carlo simulations without any primordial binaries, and indeed found β ≈ 5. In these simulations, the binaries that form are sampled from the entire BH population and should follow this distribution.¹⁰ Subsequent interactions between the new binary do not appear to alter the distribution. In our N-body simulations, we have found that β ≈ 4, when we include a background potential. However, we find β ≫ 4 when this potential was absent. We expect that this much higher boosting rate was caused by mass segregation in the core of the BH cluster, which is not present in our Monte Carlo simulations. Currently, our simulations do not include the impact of relaxation between the BH cluster and background stellar population. Nevertheless, the results presented in Morscher et al. (2015), who used a similar initial mass function, appear consistent with our results.

Overall, we find the most likely BH binaries to be detected by aLIGO have total masses slightly less than ∼2M_max when α ≲ 2. For shallow mass functions, with α < 2, the effect is most prominent. For the case where β ≈ 4 and α ≈ 2, the distribution of BH mergers is fairly broad between M_max and 2M_max. Mass segregation in the core of the BH cluster, however, can significantly boost the fraction of mergers near 2M_max. Note that aLIGO is sensitive to inspirals with masses up to 400 ${M}_{\odot }$, and it may measure the ringdown waveform for even higher M_tot (Ajith et al. 2007; Abadie et al. 2011; Kalaghatgi et al. 2015). Thus, LIGO will be capable of measuring the value of M_max over a wide range. A Bayesian analysis of the distribution of LIGO detections will allow one to put constraints on the underlying initial BH mass function and degree of mass segregation in these clusters.

Looking at the distribution of mergers in Figure 2, there is a small population of mergers with M_tot > 2M_max. These binaries must have formed from at least one BH that was involved in a previous merger. For low-spinning BHs, as we consider here, we expect low merger kick velocities when q ≈ 1, as these binaries emit GWs symmetrically (van Meter et al. 2010). These merger remnants then remain in the cluster as the largest BH, and subsequently merge with another high-mass BH in the cluster. These mergers are therefore distributed around (M_tot,q) = (2.5 M_max, 0.5). Since BH spin is expected to reach a = 0.7 ± 0.1 in circular mergers of nonspinning BHs, the heavy binary component in this population may be expected to have such high spin while the other component is nonspinning, provided that high stellar-mass BHs form with nearly zero spin (Amaro-Seoane & Chen 2016). These subsequent mergers are infrequent, ≈7% of mergers for the O1 science run of aLIGO; however. they represent a smoking gun of dynamical interactions. At the final design specifications, aLIGO will be more sensitive at detecting this population of BHs, 11% of all inspirals in dense stellar clusters will constitute subsequent mergers. For an old population of clusters, the aLIGO detection rate of subsequent mergers is closer to 15% of all detected inspirals.

R.M.O. acknowledges the support provided by NSF grant AST-1313021. This work was supported in part by the European Research Council under the European Union's Horizon 2020 Programme, ERC-2014-STG grant GalNUC 638435, and it was completed in part (by B.K.) in the Aspen Center for Physics, which is supported by NSF grant #PHY-1066293. The N-body simulations were carried out on the NIIF HPC cluster at the University of Debrecen, Hungary.