REPEATED MERGERS OF BLACK HOLE BINARIES: IMPLICATIONS FOR GW190521

Single stars and binaries evolve according to the “SSE” and “BSE” algorithms respectively (Hurley et al., 2000, 2013), which are designed to handle complex evolutionary processes, including common envelope evolution and mass transfer, collisions and supernova kicks. The version of SSE available in our simulation code utilizes the Hurley et al. (2000) supernovae scheme, which is based on work done by Eggleton et al. (1989). Under these prescriptions SSE calculates BH birth masses based on the the formulation given by Belczynski et al. (2008), whereby BH mass is soley a function of progenitor star core mass. Note that the code version we use does not include any prescriptions for oversized black holes, such as those presented in Spera et al. (2019). As stars evolve they lose mass and angular momentum due to stellar winds. These winds are calculated using the semi-empirical prescriptions of Belczynski et al. (2010) whereby mass loss rate is a function of metallicity, stellar surface temperature and luminosity. One of the key features of this model is a ”bi-stability jump” that occurs at surface temperatures of $\sim 25000$ K, which leads to a significant increase in mass loss due to an increase in the ratio of Fe III to Fe IV, which is a more effective wind driver. Binary orbital changes due to gravitational waves and subsequent merger times are calculated by integrating the orbit-averaged Peters equations (Peters, 1964), NBODY6 does not include post-Newtonian corrections to the equations of motion.

In this paper we focus on a specific simulation, with the following specifications: initial particle number $\text{N}_{0}=100,000$ particles, primordial binary fraction $f=0$, metallicity $Z$ = 0.001, initial half mass radius $r_{\mathrm{h,0}}=2.5$ pc and dimensionless King concentration parameter $W_{0}=7$. Stellar masses are drawn from a Kroupa (2001) initial mass function, between 0.08 $\text{M}_{\odot}$ and 100 $\text{M}_{\odot}$. The initial core radius is 0.56 pc with a tidal radius of 71.5 pc. The cluster initially under-fills its tidal radius by a factor of three. The tidal radius is calculated by placing the cluster into a circular orbit around a point-mass galaxy with a galactocentric distance of $23.3$ kpc. We also initialize the cluster in dynamical equilibrium and with no mass segregation. The cluster model is representative of an average mid-sized GC (see Heggie & Hut 2003 Table 1.1). The stellar metallicity of Z = 0.001 ([Fe/H] $\approx-1.33$) is approximately the median metallicity for Milky Way GCs (Harris, 1996).The cluster is evolved for a total of 12.5 Gyr.

Natal kicks are applied to neutron stars and black holes (no natal kicks are given to white dwarfs), drawn from a Maxwellian distribution with velocity dispersion $\sigma_{k}\approx 1000$ ms^-1, the initial velocity dispersion of the cluster (see de Vita et al. 2019 for further explanation). Our choice of Maxwellian velocity dispersion is significantly smaller than the default NBODY6 value of $\sigma_{k}\approx 265$ kms^-1, which is based on pulsar proper motion observations (Hobbs et al., 2005). We choose to impart low natal kicks to create conditions such that a large fraction of BHs are retained in the model clusters after formation. This is justified as it somewhat mimics the effects of fallback modulated BH natal kicks on BH retention statistics (Janka, 2013; Repetto & Nelemans, 2015; Banerjee et al., 2020). This natal kick treatment is the same for neutron stars and BHs, as the public version of NBODY6 does not contain prescriptions for mass fallback or momentum conserving BH natal kicks. One caveat of these low natal kicks is that an unusually large number of neutron stars are retained inside the host cluster after formation compared to similar studies (Ye et al., 2020; Kremer et al., 2022, 2020; Chatterjee et al., 2017). However, for most of a clusters lifetime neutron stars are too light to be part of the set of most massive particles that migrate to the cluster core. Hence the presence of neutron stars is unlikely to significantly effect the dynamics of the cluster core until very late in the cluster’s evolution. In addition, in the current study we are only interested in BHB interactions and mergers, and so the effects of the high neutron star retention fraction on our results are expected to be minimal. For information on the complete set of simulations, further detail on the code used, and full discussion of initial conditions see Anagnostou et al. 2020 and de Vita et al. (2019) (see also MacLeod et al. 2016).

where the components perpendicular to the orbital axis are defined as

and

and the component parallel to the orbital axis is

Here $A$, $K_{2}$, $K_{3}$, and $K_{S}$ are constants (units: $\text{kms}^{-1}$), $B$ is a dimensionless constant, and $\xi$ measures the angle between the unequal mass and spin contribution to the recoil velocity in the orbital plane, and is taken to be 145°(Lousto et al., 2010). The mass ratio is defined as $q=m_{2}/m_{1}\leqslant 1$, and the symmetric mass ratio as $\eta=q/(1+q)^{2}$, where $m_{1}$ and $m_{2}$ are the primary and secondary black hole masses respectively. Additionally, $\hat{e}_{1}$ is the unit vector in the direction of the separation between the two BHs at the end of the last orbit prior to the final plunge, $\hat{e}_{3}$ is the unit vector in the direction of the orbital axis, and $\hat{e}_{2}\equiv\hat{e}_{1}\times\hat{e}_{3}$. The symbols $\perp$ and $\parallel$ refer to the components of the dimensionless spin vectors perpendicular and parallel to the orbital angular momentum respectively. For example, $\alpha^{\bot}_{i}$ represents the projection of the i’th black hole’s dimensionless spin onto the orbital plane, where $i=1$ or 2. The symbol $\phi_{i}$ represents the angle between $\alpha^{\bot}_{i}$ and the displacement vector between the black holes, as measured at the end of the last orbit before the final plunge. The symbol ${\Phi_{i}}$ represents the amount by which this angle precesses before the merger. The ${\Phi_{i}}$ parameters are not constants as they depend on the specifics of a BHB inspiral, such as the mass ratio, initial separation and spins. For the purposes of the current work, we set ${\Phi_{1}}$ = ${\Phi_{2}}=0$, as suggested by Lousto et al. (2010) for statistical simulations of BHB mergers. For more details on these analytic formulae, including parameter fitting, see van Meter et al. (2010).

This post-processing implementation of recoil does come with some inherent disadvantages. Since we are not applying recoil velocities directly in NBODY6, our method only approximates the full effects of recoil kicks. Any post-merger BH that would be ejected due to recoil still remains in the cluster (as merger recoils are only calculated in post-processing), we simply ignore any future mergers involving it. Keeping these massive BHs in the cluster when they should be ejected does impact the dynamics of the cluster, particularly in the dense core. It may lower the probability of other BHB mergers, as BHs tend to primarily interact with the most massive BH in the core. However, this method does account for the most important impact of kicks when regarding hierarchical BHB mergers — the prevention of long merger chains due to the ejection of merger products.

However, there are distinct advantages to calculating recoil outside of the N-body code. With this approach, we can enable and disable recoil kicks for our set of simulations - meaning we can directly assess the impact of merger recoils on the distribution of merger chain lengths. Furthermore, our Monte Carlo method for assigning spin magnitudes inherently requires post-processing. In fact, when recoils are calculated within the simulation code, a single spin orientation must be selected, leading to a specific recoil velocity and specific outcome for the post-merger BH. Given the rarity of BH mergers in a single simulation and the computational resources required to run a single realization of a direct N-body simulation, it would be practically impossible to build sufficient information on the expected distribution of BH recoil events from run-time recoil implementation.

Once Snowball migrates to the core by $t\approx 11.7$ Myr, it experiences multiple close interactions with other mass-segregated BHs before forming its first BHB, with an eccentricity of 0.92 and a period of $P=99$ days. The binary remains in the core for 64 Myr, undergoing various three-body interactions, hardening (with the period reducing to 63 days) and altering the eccentricity before being disrupted, exchanging the secondary component for a 26.14$\text{M}_{\odot}$ BH, the most massive BH in the cluster at that time. The lower panel in Figure 3 displays the evolution of Snowball’s binary eccentricity. It is clear that the many scattering events significantly alter the BHB eccentricity, with extremes of $e=0.0008$ (near circular) to $e=0.999$ (maximally eccentric). After five more three-body interactions the binary is again disrupted at $t\approx 230$ Myr, exchanging the 26.14 $\text{M}_{\odot}$ BH for a 25.34 $\text{M}_{\odot}$ BH and ejecting the heavier former companion from the cluster, making Snowball the most massive body in the cluster. The newly formed binary is fairly eccentric, with $e=0.97$.

The process of forming a high eccentricity BHB, undergoing multiple strong interactions and an eventual exchange event occurs twice more, until Snowball forms its first merging BHB with a $m_{2}=23.2$ $\text{M}_{\odot}$ companion, corresponding to a mass ratio of 0.89. The binary forms at $t=258.7$ Myr with $P=1732$ days and $e=0.41$. After undergoing nine hardening events the binary separation is lowered sufficiently to allow for significant gravitational radiation, leading to the inspiral and eventual merger at $t\approx 300$ Myr. This merger nearly doubles Snowball mass to 37.1 $\text{M}_{\odot}$.

Soon after this first merger, at $t=322$ Myr, Snowball forms a new BHB with a 23.34 $\text{M}_{\odot}$ companion at $r=0.2r_{\text{c}}$. The remaining six BHB mergers follow a relatively similar evolution. In the top panel of Figure 3 we display the evolution of Snowball’s binary period. The seven BHB mergers are indicated with vertical black lines. The binaries often form with large periods which rapidly decrease through successive dynamical interactions with other BHs and bodies in the core. We indicate the 26 interactions that cause an exchange by dotted vertical green lines. The rate of BHB hardening tends to decrease as a binary evolves, first hardening rapidly and then slowly reducing the separation until an inspiral occurs, wherein most of a binaries life is spent with small separations. This is partly because the cross-section for three-body interactions of hard binaries is proportional to the semi-major axis (Celoria et al., 2018). Hence the hardening rate decreases as the binary orbit tightens.

Although most scattering events follow Heggie’s law (Heggie, 1975), some soften the BHB involved instead. Most of these softening events occur due to an exchange. The attendant interval of chaotic three-body dynamics can dramatically alter the orbital parameters of the newly formed binary, leading to significant softening. This is most evident in the string of exchange events that occur between $400\lesssim t\lesssim 500$ Myr, leading up to the second BHB merger. The top panel of Figure 3 shows the effects of such interactions, whereby the period varies non-monotonically between 320 and 592 days.

Scattering events harden binaries and can impart significant recoil velocities which send even IMBHs out of the cluster core (de Vita et al., 2018). Anagnostou et al. (2020) found that $\sim 80\%$ of all ejections follow a three-body encounter with another BH, with the rest being ejected via binary-binary scattering. Although the various BHBs which count Snowball as a member experience significant three-body recoil, the resulting velocities are not sufficient to eject the binary. As can be seen in Figure 1, at most, the BHB is flung outside the cluster core, where it migrates back towards the centre to undergo further scattering events. One particular interaction leads to a Snowball binary being flung out to $8.45r_{c}\approx 1.3$ pc when Snowball is $\sim 82$ M_⊙. Although Snowball is never ejected, over the $12.5$ Gyrs approximately $67$ binary-single interactions lead to the ejection of a BH. In particular, Snowball plays a large part in the depletion of BHs in the core, being responsible for ejecting $43$ BHs. As the number of BHs in the core is depleted over time (starting from an initial $N_{BH}=74$ formed from stellar evolution and retained in the cluster after natal kicks) the number of three-body hardening events likewise decreases.

Why is it that no other BHs merge with anything apart from snowball in this cluster?. Throughout the simulation, 10 other BHBs not containing Snowball form in the core, however, none remain bound long enough to either merge or be ejected, with the longest binary only surviving $\sim 31$ Myr. All 10 binaries are disrupted by encounters with Snowball, or with a binary containing Snowball. Seven out of these 10 BHBs undergo an exchange encounter with Snowball, wherein the preexisting binary exchanges a component for Snowball, with the remaining three BHBs simply destroyed with no exchange. As a result of the exchange events, six of the seven newly unbound BHs are ejected from the cluster as a result of the interaction. Snowball’s presence in the core is essentially preventing the survival of any other BHB. Exchange encounters are more frequent within the first few gigayears when there is still a sufficiently high number of relatively massive BHs in the core, with only four out of 28 exchanges occurring after $t=2$ Gyr. This is because exchange events, in general, require BHs of comparable mass (Heggie, 1975). On the other hand, scattering events that don’t result in an exchange remain relatively common, occurring every few Myr at least until $\sim 10$ Gyr, after which the BH population has mostly evaporated (see also Morscher et al. 2013 for a Monte Carlo investigation of the BH population in star clusters).

The specifics of Snowball’s mass gain depend on the mass of the companion BHs involved in the merger. However, our results clearly show that any BHB within a cluster core is likely to undergo multiple exchange events. Previous scattering experiments also show that the lighter of the two binary components is preferentially ejected from the system in favour of the intruder (Sigurdsson & Hernquist, 1993), meaning that BHBs with mass ratios closer to one are dynamically favoured. This process systematically increases the companion mass before merger. Hence, heavier BHs are more likely to merge than lighter ones, which is indeed what we see here. All seven Snowball BHB mergers involve a secondary BH which is in the heaviest 5% of all BHs at the time of the merger. This is why most of the BHs that merge with Snowball have masses in the low to mid 20 M_⊙, as this is approximately the maximum mass of BHs produced due to our cluster models and supernovae prescriptions.

At $t\approx 7.9$ Gyr, this final binary forms a hierarchical triple system. The system undergoes Kozai-Lidov oscillations (Lidov, 1962; Kozai, 1962), driving the inner binary to as high as $e=0.99$ (see the red highlighted section in Figure 3). During this time the inner binary still experiences 10 three-body scattering events. Throughout the next 2 Gyr, the outer component of the triple system is exchanged nine times, ejecting the previously bound outer component from the cluster in all but one of these cases. Eventually, the BHB triple system is disrupted at $t\approx 9.94$ Gyr.

After 10.8 Gyr Snowball and its companion are the only remaining BHs in the cluster. During the last $\sim 200$ Myr, the BHB once again forms a triple system with a white dwarf. This hierarchical triple remains until the end of the simulation.

BH spins are assigned retroactively in post-processing and are used to calculate BH merger recoils (see section 2.2). The two first-generation BHs involved in the initial merger in the Snowball chain are assigned spin magnitudes of $\alpha_{1}=0.08$ and $\alpha_{2}=0.18$, based on the MESA code. From the merger between these two BHs, the resulting second-generation BH is ejected from the cluster in approximately $85.2\%$ of the $10^{4}$ merger recoil calculations. In approximately $14.5\%$ of the $10^{4}$ recoil calculations the second generation BH merges again, and the resulting third generation BH is ejected. In approximately $0.430\%$ of the $10^{4}$ recoil calculations the third generation BH merges and the resulting fourth generation BH is ejected. The longest Snowball chain, comprising four mergers, occurs twice in $10^{4}$ recoil calculations.

Although gravitational recoils eject Snowball after at most only four mergers, longer merger chains can occur rarely. One could foresee a scenario whereby an initial merger forms a BH that is not only retained but is also substantially larger than the remaining BH population. This massive second-generation BH could then merge multiple times with relatively light BHs without being ejected, as the low mass ratios would lead to low recoil velocities.

This is expected, as longer chains are less likely regardless of gravitational recoil, as longer chains are multiplicative in nature. At a most basic level, we expect the probability of a certain chain length to follow

Without recoil, approximately half of the 97 in-cluster mergers only involve BHs produced directly through core collapse. 18 of the resulting second-generation BHs never merge again, increasing to 33 when merger recoils are accounted for (without BH spin). 47 mergers involve a BH with $\text{BH}_{\text{gen}}\geqslant 2$. In 36 of these mergers $m_{2}$ is a first-generation BH, with $m_{1}$ being the multi-generational object. This is expected, as hierarchically formed BHs are on average heavier than those formed through core collapse, and thus are more likely to be the heavier element in a BHB. Mergers involving an object with BH${}_{\text{gen}}\geqslant 2$ involve a first-generation BH merging with the higher generation BH, with only three exceptions. One exception is a merger between two second-generation BHs, and the other two are between a second and third-generation BH. These three mergers represent the only such mergers between two previous merger remnants. When we take into account BH spins, the only merger between two multi-generational BHs is the merger between the two second-generation BHs. Hence we conclude two or more chains are unlikely to exist simultaneously within a mid-sized GC ($\sim 10^{5}$ stars). This is expected, as once a BHB merges, its product is likely to be the most massive body in the cluster, and thus sink to the cluster centre due to dynamical friction, even if merger recoil is imparted (as long as the resulting BH remains bound). The massive BH dominates dynamically within the dense core, likely undergoing many dynamical interactions with other single and binary BHs. Hence, any BHB that forms is likely to interact with the second-generation BH, in which case an exchange event will preferentially unbind $m_{2}$, forming a new BHB with the second-generation BH (Sigurdsson & Hernquist, 1993). It is important to note that recoil velocity is approximately independent of the total binary mass, but does depend on $q$. Hence light and heavy mergers have the same chance of reaching a given escape velocity.

Including recoil significantly increases the number of post-merger BH ejections, resulting in most larger merger chains no longer being viable. The $L=7$ chain no longer occurs, now only reaching a chain length of four at most before the BH is ejected. Even without spin, the recoil is often large enough to eject the majority of the 50 second-generation BHs, with only 23 being retained, 17 of which undergo further mergers. Only five of the 13 third-generation BHs are retained, all of which merge again. This suggests that any retained BHs with BH${}_{\text{gen}}>1$ have a high chance of merging again. Including spins in the recoil calculations tends to increase the size of the kicks, with now only approximately $0.731\%$ of the second-generation BHs retained after formation, $\approx 73.9\%$ of which merge again. Approximately $85\%$ of the resulting $L=3$ chains correspond to a merger between two BH${}_{\text{gen}}=2$ objects, forming a $58\text{M}_{\odot}$ BH. This is likely the main way an $L=3$ chain can occur in GCs of these sizes, as the other pathway involves a BH ${}_{\text{gen}}=1$ object merging with a BH${}_{\text{gen}}=3$ object (as in the Snowball merger chain). This pathway is less likely as it requires the retention of a BH${}_{\text{gen}}=3$ object, which only occurs approximately $0.450\%$ of the time in the case of Snowball’s merger chain. If we take the most likely recoil for each merger based on all spin orientations, we find only 58 BHB mergers, meaning that only $3.3\%$ of all the BHs that remain bound after natal kicks end up merging.

where $R_{*}$ is the comoving star formation rate [units: $M_{\odot}\text{yr}^{-1}\text{Mpc}^{-3}$] at the characteristic time of star cluster formation and within the volume observed by LIGO, $f_{\text{GC}}\approx 10^{-3}$ is an order-of-magnitude estimate of the fraction of star formation that occurs within a dense star cluster, $f_{\text{Snowball}}$ is the fraction of those clusters that are likely to produce Snowball-like merger chains, estimated from our simulations, and $M_{GC}\approx 6.5\times 10^{4}M_{\odot}$ is the characteristic mass of a dense star cluster in our runs.

A detailed estimation of $R_{\text{snowball}}$ would require to take into account both the redshift dependence of $R_{*}(z)$ and the distribution of time delays between the formation of the cluster and a chain merger event. However, uncertainties in $f_{\text{GC}}$, $f_{\text{Snowball}}$, and $M_{\text{GC}}$ dominate Eq. 6, thus to remain in the spirit of an estimate to only an order of magnitude, we can neglect redshift dependency for $R_{*}$, adopting an appropriate average value for it, and assume the merger time delay distribution is a delta function at the cluster’s age. Given that LIGO is sensitive to high mass merger events out to $z\sim 1-1.5$ (Chen et al., 2021), and that most of the volume is thus at cosmological distances, we assume $R_{*}\approx 0.1M_{\odot}\text{yr}^{-1}\text{Mpc}^{-3}$, corresponding to the comoving star formation rate at $z\approx 1.2$ from Madau & Dickinson (2014, Figure 9, left panel). We also neglect the effects of a cluster birth mass function term. For our basic estimate, we simply ignore clusters with stellar masses smaller than our simulations, and use characteristic values for the population of globular clusters in the Milky Way. To an order of magnitude, a more detailed calculation would give the same result as although $f_{\text{Snowball}}$ would be smaller in lower mass clusters, one would have to increase $f_{GC}$ to reflect the fact that a greater proportion of star formation is considered in the calculation if lower mass clusters are included.

Based on the results presented in Table 1, we assume $f_{\text{Snowball}}\approx 0.02$ (one IMBH chain event out of over 50 simulations). A more precise characterisation is challenging from a single chain instance discovered in our set of simulations. Given these assumptions, we estimate the Snowball event rate to be

This is within the 90% confidence interval of the inferred rate of mergers similar to GW190521, estimated by Abbott et al. 2020 as $0.13\pm^{0.3}_{0.11}\text{Gpc}^{-3}\text{yr}^{-1}$. We can take $f_{\text{Snowball}}\approx 0.02$ as an estimate for the upper limit for our cluster models. The fact that this upper limit does not lead to an event rate exceeding the observed LIGO rate is reassuring to qualitatively support the plausibility of IMBH formation through chain mergers in dense stellar systems, but of course, a deeper quantitative investigation is needed to further assess the scenario we propose.

Of course, the above agreement is only for the special baseline case without merger recoil. Assuming the rough estimate for a Snowball chain rate of $P(L=7)\approx 10^{-12}$ with spin recoil, presented in Section 5.1, we obtain $f_{\text{Snowball}}\approx 0.02/(8.8\times 10^{9})=2.3\times 10^{-12}$, and hence

Although this event rate is orders of magnitudes bellow even the lower bound of the inferred LIGO rate, an under-prediction is expected, as there are other potential GW190521 formation channels. Our estimates are also based on our cluster models, and one would expect $f_{\text{Snowball}}$ to be larger for more massive clusters with larger escape velocities. Note that Equation 8 is simply based on Equation 5, extrapolated out to $L=7$.