Binary Black Hole Mergers and Intermediate-mass Black Holes in Dense Star Clusters with Collisional Runaways

Intermediate-mass black holes (IMBHs) are believed to be the missing link between the supermassive black holes (BHs) found at the centers of massive galaxies and BHs formed through stellar core collapse. One of the proposed mechanisms for their formation is a collisional runaway process in high-density young star clusters, where an unusually massive object forms through repeated stellar collisions and mergers, eventually collapsing to form an IMBH. This seed IMBH could then grow further through binary mergers with other stellar-mass BHs. Here we investigate the gravitational-wave (GW) signals produced during these later IMBH–BH mergers. We use a state-of-the-art semi-analytic approach to study the stellar dynamics and to characterize the rates and properties of IMBH–BH mergers. We also study the prospects for detection of these mergers by current and future GW observatories, both space-based (LISA) and ground-based (LIGO Voyager, Einstein Telescope, and Cosmic Explorer). We find that most of the merger signals could be detected, with some of them being multiband sources. Therefore, GWs represent a unique tool to test the collisional runaway scenario and to constrain the population of dynamically assembled IMBHs.


INTRODUCTION
Black holes (BHs) are observed commonly in two mass ranges: stellar-mass BHs, formed through massive star evolution to core collapse, with masses M BH ≲ 100 M ⊙ (e.g., Celotti et al. 1999;Remillard & McClintock 2006; The LIGO Scientific Collaboration et al. 2021), and supermassive BHs found at the centers of most massive galaxies, with M BH ≳ 10 5 M ⊙ (e.g., Tremaine et al. 2002;Marconi & Hunt 2003;Kormendy & Ho 2013).BHs with masses in between these two regimes are labeled intermediate-mass BHs (IMBHs); for a review see Greene et al. (2020).While IMBHs could play a fundamental role in the evolution of galaxies and could be a source of tidal disruption events and gravitational waves (GWs), their existence has not been confirmed beyond a reasonable doubt (e.g., Jardel & Gebhardt 2012;Neu-mayer & Walcher 2012;Graham & Scott 2013;Mezcua 2017;Nguyen et al. 2018;Perley et al. 2019;Smith et al. 2023).
IMBHs have masses beyond the most massive BH that can be produced as a result of direct stellar core collapse.Current stellar evolution models predict a dearth of BHs with masses larger than about 50 M ⊙ as a result of pairinstability physics, where pair production removes pressure support in the core (Heger et al. 2003;Woosley 2017).Whenever the pre-explosion stellar core is in the mass range M ⋆ ∼ 45−65 M ⊙ at the onset of their carbon burning, large amounts of mass can be ejected, leaving a BH remnant with a maximum mass around 50 M ⊙ .Larger stellar cores can trigger thermonuclear reactions that can completely destroy the star and leave no remnant behind.
IMBHs are primary sources for present and future GW observatories.Recent detections by the LIGO-Virgo-KAGRA (LVK) Collaboration have found binary BH (BBH) mergers where one or both components of the merging binary have masses above 50 M ⊙ .Among them, GW190521 is the most interesting event since its remnant has a total mass of ∼ 150 M ⊙ , nominally in the IMBH regime (Abbott et al. 2020).One of the main venues for the assembly of these binaries is the core of dense star clusters, where a massive BH could form through collisions and mergers of massive stars, and through hierarchical BH mergers (Antonini et al. 2016;González et al. 2021;Chattopadhyay et al. 2023;Fragione & Rasio 2023;Atallah et al. 2023).With hundreds of new detections expected in the current O4 run by the LVK Collaboration, studying BBH mergers across all relevant mass ranges has become a key priority.
In this paper, we model repeated mergers of BBHs in dense star clusters that we assume to have undergone at early times (typically within just a few Myr after cluster formation) a collisional runaway process leading to the formation of a very massive star (with M ⋆ ∼ 10 2 − 10 3 M ⊙ ) and, ultimately an IMBH.We perform simulations of merging BBHs formed in these dense clusters using the semi-analytic method developed in Fragione & Rasio (2023).The remnants of the merging BBHs grow the IMBH to final masses ≥ 10 3 M ⊙ .We explore for the first time the implications of such a runaway process for GW detections by present and future observatories.
Our paper is organized as follows: in Section 2 we discuss our semi-analytic method to study BH mergers in dense star clusters hosting an IMBH.In Section 3 we present our numerical results, and we summarize our key findings and implications for future work in Section 4.

METHODS
In this section, we include a summary of the semianalytic approach we use to derive our results.For a detailed discussion, see § 2 of Fragione & Rasio (2023).
We consider a dense stellar cluster of initial mass M CL , in the range [10 5 , 10 7 ] M ⊙ , and half-mass-radius r h = 1 pc, described by a King model with initial moderate concentration (King 1962).This is representative of the typical size of a young massive cluster (Portegies Zwart et al. 2010).We sample stellar masses from the canonical Kroupa initial mass function (Kroupa 2001): (1) Given this IMF, we produce a total of BH progenitors, in the mass range [20M ⊙ , 150M ⊙ ].We consider cluster metallicities Z ∈ {0.02, 0.002, 0.0002}, and evolve the BH progenitors at a given metallicity using the stellar evolution code sse (Hurley et al. 2000), with all the necessary updated prescriptions for stellar winds, stellar interactions, and formation of remnants (for details see Banerjee et al. 2020).
We assume that all BHs are born with negligible natal spins, consistent with the findings of Fuller & Ma (2019).In our simulations, each BH is imparted a natal kick at birth sampled from a Maxwellian distribution with velocity dispersion 265 km s −1 (Hobbs et al. 2005), scaled by a factor of 1.4M ⊙ /M BH to account for momentum conservation (Fryer & Kalogera 2001).If the natal kick exceeds the cluster escape velocity (Fragione & Rasio 2023) we consider the BH ejected and remove it from our simulation.The same applies for dynamical kicks in threebody encounters and recoil kicks imparted to the remnant of a BBH merger, caused by anisotropic GW emission (see e.g., Lousto et al. 2010).We model cluster evolution by following Antonini & Gieles (2020a,b).Briefly, the cluster is assumed to reach a state of balanced evolution, so that the heat generated by the BBHs in the core and the cluster global properties are related.Our simulations include all the fundamental elements of cluster evolution (cluster mass loss and expansion) and the fundamental processes of formation and evolution of BBHs (3-body interactions, mergers, recoil kicks, etc.).In Fragione & Rasio (2023), it was shown that this semi-analytic method performs well at reproducing the essential elements of BBH mergers, especially when compared to N -body simulations.For more details, see Antonini & Gieles (2020a,b) and Fragione & Rasio (2023).
In our study, we build on the scheme presented by Fragione & Rasio (2023) by adding a new parameter f run , which represents the fraction of total cluster mass that participates in the initial runaway process.The mass of the star formed as a result of the runaway is then simply (4) Many previous works have tried to estimate f run from Monte Carlo N-body simulations, finding a typical value of ∼ 10 −3 (Freitag 2001;Gürkan et al. 2004).We assume that every cluster goes through a runaway process initially, regardless of its initial mass and density.Note that however runaways typically develop in clusters with high densities, with other parameters such as the slope of the initial mass function, primordial binary fraction in massive stars, etc. that could play a critical role (e.g., Fregeau et al. 2002;Ivanova et al. 2005;González et al. 2021).
In our simulations, we consider four different cases, with f run ∈ {0, 0.0005, 0.001, 0.005, 0.01}.For each value of f run , we run 5000 simulations for each metallicity.The clusters have masses in the range M CL ∈ [10 5 , 10 7 ] M ⊙ assuming a distribution of cluster masses ∝ M −2 CL (Portegies Zwart et al. 2010).With our simulations, our goal is to measure the rates of BBH mergers for BHs formed through direct collapse of large stars, study the population of the merger remnants by tracking the growth of IMBHs, and the probability of the detection of these merger events using GW instruments.

RESULTS
In this Section, we summarize the results of our simulations and discuss the detectability of the BBH mergers by present and upcoming GW observatories.

Growth of an IMBH
In our simulations, we initialized our star clusters with a massive BH remnant of the collapse of a large star resulting from the runaway process (e.g., Portegies Zwart & McMillan 2002;Portegies Zwart et al. 2004;Gürkan et al. 2006;Giersz et al. 2015a;Mapelli 2016;González et al. 2021).This BH may grow to an IMBH through repeated mergers with other BHs over time, whenever not ejected as a consequence of dynamical kicks in fewbody interactions or recoil kicks after a merger via GW emission.
We track the growth of the most massive BH in each cluster as a result of repeated mergers with other stellar BHs.In Figure 1, we show the growth of IMBHs (ratio of its mass to the initial mass) in clusters with initial cluster mass M CL = 10 6 M ⊙ as a function of time for different values of f run .For clusters with higher values of f run , the ratio does not increase significantly; IMBHs  in clusters with smaller values of f run grow significantly over time, reaching 2-3 times their initial mass.
In Figure 2, we show the most massive BH formed in a star cluster as a function of the initial cluster mass Fractional number of clusters for various values of the fraction of cluster mass that undergoes a runaway process.For simulations sampling all three metallicities, we make bins of cluster masses and select the most massive BH in each bin.As expected, a larger value of f run implies a larger most massive BH across all initial cluster masses.For example, we find that a cluster of initial mass of 10 5 M ⊙ , f run = 0 leads to a most massive BH of about 80M ⊙ , f run = 0.001 of about 1.1 × 10 3 M ⊙ , and f run = 0.01 of about 4.5 × 10 3 M ⊙ .The value of f run is particularly crucial for small clusters.In particular, this is because for smaller values of f run the smaller IMBHs are ejected from the cluster due to dynamical or recoil kicks after the mergers.For example, in the case there is no runaway (f run = 0), clusters with masses ≲ 10 6 M ⊙ do not form an IMBH; however, the most massive IMBH for a cluster of mass ∼ 10 6 M ⊙ has a mass of about 5 × 10 3 M ⊙ for f run = 0.001, while a mass of about 3 × 10 4 M ⊙ for f run = 0.01.
To provide additional insight into the formation of massive BHs within clusters, Figure 3 illustrates the fraction of clusters, or likelihood of the most massive BH (referred to as m 1 in Figure 3) surpassing a specified mass threshold in our simulations.When f run = 0 (no runaway), only very massive clusters exhibit a nonnegligible probability of forming massive BHs.Conversely, with higher f run values, the initial mass of the IMBH is larger, increasing the likelihood of their appearance even in less massive clusters, as ejections also become less important.

BBH merger rates
We follow Fragione & Rasio (2023) and calculate the rates (in units of Gpc −3 yr −1 ) of BBH mergers as 1 where t lb is the look-back time at redshift z, N events is the number of events as a function of the initial cluster mass M CL , initial half-mass radius r h , metallicity Z, and formation redshift z f , and Ψ(M CL , r h , Z, z f ) is a probability function that weighs the the previous cluster properties.In our model, we take the distribution of cluster masses to be ∝ M −2 CL , with the maximum possible cluster mass being M CL = 10 7 M ⊙ (e.g., Portegies Zwart & McMillan 2000).As discussed above, we fix the half-mass radius at r h = 1 pc, which follows the typical distribution of observed values for local, young stellar clusters (Portegies Zwart et al. 2010).We take the formation times to be proportional to exp[−(z−z f ) 2 /(2σ 2 f )] where z f = 3.2 and σ f = 1.5, reminiscent of cluster formation times as inferred from cosmological simulations (Mapelli et al. 2021).The cluster masses normalized such that the cluster density is ∼ 1 Mpc −3 in the local Universe Portegies Zwart et al. (2010).Metallicities are sampled from a log-normal distribution with mean given by (Madau & Fragos 2017) log⟨Z/Z ⊙ ⟩ = 0.153 − 0.074 z 1.34  (6) and a standard deviation of 0.5 dex.In Equation 5, K accounts for the cluster density evolution, considering that a fraction of the star clusters that are formed in the Universe evaporate across cosmic time.We fix K = 32.5, consistent with Antonini & Gieles (2020a) and with the value needed to reproduce the merger rates of BBH mergers in the latest LVK catalog (Fishbach & Fragione 2023).
We calculate the BBH merger rate as a function of redshift for each values of the runaway fraction.Figure 4 shows the overall merger rate for our models, which we also break down to show the contribution of IMBH of different masses (m 1 > 100, 200, 500, 1000 M ⊙ ).We find that relative contribution of massive BHs is larger for larger values of f run , since the IMBH starts dominating the merger rate, as it keeps merging repeatedly with the other stellar-mass BHs.Indeed, we find that f run = 0 has a merger rate of 19 Gpc −3 yr −1 , f run = 0.001 of 6.2 Gpc −3 yr −1 , f run = 0.005 of 2.9 Gpc −3 yr −1 , and f run = 0.01 of 2.24 Gpc −3 yr −1 .At the same time, the merger rate for primary masses larger than 1000M ⊙ for f run = 0 is of 0.09 Gpc −3 yr −1 , for f run = 0.001 of 1.66 Gpc −3 yr −1 , for f run = 0.005 of 2.03 Gpc −3 yr −1 , and for f run = 0.01 of 2.24 Gpc −3 yr −1 .

Detectability of merging BBH with GW instruments
Given a population of merging BBHs, we must ask whether current or future ground-based and space-based instruments are likely or not to detect some of them.We calculate the detection fraction of these merger events for planned ground-and space-based GW instruments, such as the Laser Interferometer Space Antenna (LISA; Robson et al. 2019), Voyager (the upgraded detector of the current LIGO facility, The LIGO Scientific collaboration 2019), the Einstein Telescope (ET; Punturo et al. 2010), and the Cosmic Explorer (CE; Reitze et al. 2019).While LVK could detect the mergers of BBHs with ∼ 100M ⊙ out to redshifts of z ∼ 1 (Abbott et al. 2019), the upcoming space-based LISA and ground-based instruments like ET and CE offer the opportunity of detecting the formation of IMBHs up to z ∼ 100 (Amaro-Seoane et al. 2017;Jani et al. 2020;Fragione & Loeb 2023).Thus, understanding the probability of detection of these events will guide searches for massive BBH systems and help plan future observing cycles.
Following Fragione & Loeb (2023), we design the detection fraction as an instrument-dependent function, which encodes the ability of observing the merger of a binary with primary mass m 1 and mass ratio q at a redshift z.The masses and merger times are obtained directly form our simulations.However, the merger times need to be corrected for the cluster formation time.As for the merger rates, we take the formation times to be proportional to exp[−(z − z f ) 2 /(2σ 2 f )] where z f = 3.2 and σ f = 1.5 described in §3.2.Once we correct for the cluster formation, we discard those binaries that have a merger time larger than the age of the Universe.The sources not discarded are kept in our data and their merger times are converted to redshifts.
The detection of a merger event is modeled as an instrument-dependent function F det (z, m 1 , q) as where H is the Heaviside function and ⟨ρ(z, m 1 , q)⟩ is the averaged signal-to-noise (SNR) ratio.The threshold for a detection is set to an SNR of 8. From ) f run = 0.01 Figure 4.The binary black hole (BBH) merger rate as a function of redshift for 4 different fractions of initial cluster mass that undergo the collisional runaway.The total merger rate is shown, along with rates for a primary mass over a given threshold for the primary mass: 100, 200, 500, 1000M⊙.
Loeb (2023), we compute the average SNR as where C is obtained after averaging over various sky relations, with C = 2/ √ 5 and C = 2/5 for space-based and ground-based detectors, respectively (Robson et al. 2019); f min and f max are the minimum and maximum frequency of the binary in the detector band, respectively; S n (f ) is the noise power spectral density; | h(f )| is the frequency-domain waveform amplitude for a face-on binary.We use pyCBC developed by Nitz et al. (2019) with the IMRPhenomD approximant (Husa et al. 2016) to model the waveform of the merging binary BHs.
We study the detection probability of these merger events using four different gravitational wave observatories: space-based LISA and ground-based Voyager, ET, and CE.The power spectral density of LISA is derived as in Robson et al. (2019), of Voyager as in The LIGO Scientific collaboration (2019), of ET as in Punturo et al. (2010), and of CE as in Reitze et al. (2019).LISA has a planned duration of 5 years, and a binary will evolve towards a merger event as it starts from the frequency f in which is the GW frequency at the start of the observation This initial frequency is higher than the minimum detectable frequency for LISA which is conventionally 10 −5 Hz.Thus, we take, f min = f ini for LISA.For the groundbased instruments in the study namely Voyager, ET, and CE, this initial frequency f ini is typically smaller than the frequency at which the instruments start operating, 5 Hz, 1 Hz, and 5 Hz, respectively.
In Figure 5, we show the detection probability of BBH merger events in clusters, with f run = 0.001.We see that LISA is able to detect most of the mergers with primary mass ≳ 10 3 M ⊙ and mass ratio 10 −2 ≲ q ≲ 10 −1 .For the smaller masses, LISA is only able to detect a few of the merging IMBH binaries.Conversely, the population of merging binaries with m 1 ≲ 10 3 M ⊙ can be studied using ground-based instruments.LIGO's Voyager is only able to detect a few of the merging binaries in this regime, in particular when the mass ratio is close to unity.ET and CE are able to detect most binaries with primary mass m 1 ≲ 10 3 M ⊙ and mass ratio ∼ 10 −1 .
The plots showing detection fractions for other values of f run are included in the Appendix.

CONCLUSION AND DISCUSSION
In this paper, we analyzed simulations of merging BBHs in dense clusters, where an IMBH is assumed to have formed as result of the collapse of a very massive star created from repeated stellar mergers in a collisional runaway.We can summarize our key findings as follows: • The initial IMBH, born as the remnant of a very massive star formed through the collisional runaway, grows in time through repeated mergers with stellar-mass BHs, up to about 2-3 times its initial mass for small values of f run (the fraction of initial cluster mass that goes into the collisional runaway).
• Through collisional runaways star clusters may produce IMBHs with initial masses ∼ 10 2 -10 3 M ⊙ , which further grow to ≳ 10 3 M ⊙ through repeated mergers with other stellar-mass BHs.This is not seen in any of our models when we set f run = 0, i.e., when BHs grow solely through BH mergers.
• The merger rate of BBHs dynamically assembled in dense star clusters tends to decrease if the mass of the IMBH formed as a result of a runaway process is larger.For sufficiently large masses, the IMBH dominates the merger rate and the merger rate of stellar-mass BHs becomes negligible.
• The merging binaries that LISA can potentially detect map the underlying astrophysical population for primary masses m 1 > 100 M ⊙ and mass ratios 1 > q > 10 −2 .The ground-based instruments ET and CE will be capable of observing the higher mass ratio binaries, with primary mass m 1 < 10 3 M ⊙ and mass ratio q > 10 −2 , while the detection efficiency of LIGO Voyager is typically very small.
Our work considered a potentially important effect in the formation of massive BHs and their mergers that previous models have neglected.However, there are some caveats that we leave to future work.For example, we have not included in our treatment the interaction of the IMBH with ordinary stars, which could lead to interesting phenomena, such as tidal disruption events (e.g., Liu et al. 2009;Fragione et al. 2021;Angus et al. 2022;Kıroglu et al. 2023).We have also fixed some of the distributions that describe the birth properties of star clusters, consistent with observed properties in the local universe (Portegies Zwart et al. 2010), but these are poorly constrained elsewhere.Finally, clusters with low densities are less likely produce runaways, while we assume that all clusters undergo an initial runaway.
Furthermore, we have used the IMRPhenomD as our approximation for the waveform, which only includes the dominant harmonic (l, m) = (2, 2) of the GW signal.However, higher-order harmonics could contribute to the GW signal significantly, especially for IMBH binaries (Jani et al. 2020).Finally, another important yet currently uncertain question is the fate of the very massive star (with mass ∼ 10 2 − 10 3 M ⊙ ) that is produced as result of the runaway process, as stellar evolutionary models have not been calibrated over that mass range.
As Mapelli (2016) notes, the mass of the final remnant of the massive star is in the IMBH mass range only if the mass loss due to stellar winds or hydrodynamic processes is moderate, and only if the very massive star undergoes direct collapse to a BH.These conditions impose many new restrictions.GW detection offers an unparalleled opportunity to survey the Universe and detect IMBHs in various mass ranges, making it possible for the first time to constrain their formation, growth, and merger history across cosmic time.
From our analysis of merging BH-IMBH binaries in dense star clusters, we find that the merger rates are ∼ 10 Gpc −3 yr −1 and peak around redshift z ≃ 2. Thus, we expect to detect several IMBH mergers per year with upcoming GW observatories like LISA, LIGO's Voyager, the Einstein Telescope, and the Cosmic Explorer.Analysis of the merger events will provide important constraints on the underlying astrophysical populations of the merging binaries and will help constrain the IMBH formation and growth processes in dense star clusters.

Figure 1 .
Figure1.The growth of IMBHs (ratio of its mass to the initial mass) in clusters with initial cluster mass MCL= 10 6 M⊙ as a function of time for different values of frun.For clusters with higher values of frun, the ratio does not increase significantly, while IMBHs in clusters with smaller values of frun grow significantly over time.

Figure 2 .
Figure 2. Solid lines: most massive BH formed in a cluster as a function of the initial cluster mass MCL for different values of the fraction of cluster mass that undergoes initial collisional runaway, frun.The half-mass radius is fixed at r h = 1 pc representative of the typical size of young massive clusters.Dashed lines: initial IMBH mass, mrun = frunMCL, as a function of MCL.The simulations contain three metallicities as Z ∈ {0.02, 0.002, 0.0002}.

Figure 3 .
Figure3.The fraction of clusters that form a massive BH as a function of the initial cluster mass.As expected for larger values of frun, this fraction increases, especially for lower cluster mass.We show 5 different mass thresholds for what we call a massive BH: 100, 500, 10 3 , 5 × 10 3 , 10 4 M⊙.

Figure 5 .
Figure 5. Detection fractions for merging IMBH binaries as a function of primary mass and binary mass ratio.This figure shows our results for clusters with frun= 0.001 (see Appendix for other values).Top left: LISA; top right: Voyager; bottom left: Einstein Telescope; bottom right: Cosmic Explorer.