Breaching the Limit: Formation of GW190521-like and IMBH Mergers in Young Massive Clusters

The third observational campaign ran by the LIGO-Virgo-Kagra Collaboration (LVC) led to the discovery of four black hole binary (BBH) mergers with a total mass >100 M_⊙, i.e., in the mass range of intermediate-mass black holes (IMBHs) (Abbott et al. 2021). All these mergers are characterized by a primary having a mass >64 M_⊙, in a range of masses where stellar evolution theories predict the absence of such compact objects. Massive stars developing a helium burning core in the range of ∼64–135 M_⊙ end their lives in disruptive pair instability supernovae (PISNe) leaving no BH remnant behind (Fowler & Hoyle 1964; Woosley et al. 2007). If the helium core falls in the 32–64 M_⊙ mass range, the star develops rapid pulses that enhance mass loss before the star goes supernova, a process called pulsational pair instability supernova (PPISN; Woosley 2017). These processes result in the so-called upper-mass gap with no stellar BHs in the mass range of ∼64–120 M_⊙ (see, e.g., Spera & Mapelli 2017; Woosley 2017).

The lower edge of the gap is highly uncertain and can be shifted down to 40 M_⊙ or the gap might not even exist, depending on the metallicity, the rotation of the progenitor star, and the rates of carbon–oxygen nuclear reactions (see, e.g., Farmer et al. 2019; Stevenson et al. 2019; Mapelli et al. 2021; Costa et al. 2021; Woosley & Heger 2021). The physical mechanisms that govern core overshooting and stellar winds can make the picture even more complex, especially for metal-poor stars with masses >90 M_⊙, which could collapse to BHs with masses right in the upper-mass gap (e.g., Tanikawa et al. 2021; Vink et al. 2021)

At least one of the four LVC gravitational wave (GW) sources, GW190521, has a primary mass >65 M_⊙ at 99% credibility and a companion with mass ${66}_{-18}^{+17}$ (LIGO Scientific Collaboration 2020), thus potentially both the GW190521 components fall in the upper-mass gap.

There are several potential formation channels for the components involved in such type of mergers, like formation from Population III stars (Liu & Bromm 2020; Kinugawa et al. 2021; Tanikawa et al. 2021), binary evolution (Belczynski 2020), stellar collisions (Mapelli 2016; Di Carlo et al. 2019; Kremer et al. 2020; Renzo et al. 2020; Rizzuto et al. 2021), and hierarchical mergers (Gerosa & Berti 2017; Antonini et al. 2019; Rodriguez et al. 2019; Arca-Sedda 2020; Arca-Sedda et al. 2020b; Belczynski & Banerjee 2020; Doctor et al. 2020; Fragione et al. 2020; Mapelli et al. 2021) in young massive clusters (YMCs), globular clusters (GCs) and nuclear clusters (NCs), or in active galactic nuclei (McKernan et al. 2018; Yang et al. 2019; Graham et al. 2020; Secunda et al. 2020; Tagawa et al. 2021).

Dense star clusters represent an ideal environment to nurture mass-gap BHs via the dynamical assembly mechanism, either through stellar collisions or hierarchical mergers. However, the modeling of such environments have relied so far on either direct N-body simulations, tailored to low-mass clusters, i.e., lighter than 5 × 10⁴ M_⊙ (e.g., Banerjee 2018; Di Carlo et al. 2019), or adopting outdated stellar evolution prescriptions (Wang et al. 2016; Anagnostou et al. 2020a) on Monte Carlo methods (e.g., Giersz et al. 2015, 2019; Rodriguez et al. 2019; Kremer et al. 2020), or on semi-analytic approaches (e.g., Gerosa & Berti 2017; Arca-Sedda 2020; Fragione et al. 2020). In this work, we present and discuss the properties of BH binary mergers with total masses in the range of 150–350 M_⊙, which we identify in a suite of 80 direct N-body simulation modeling compact clusters composed of 110,000 stars (Rizzuto et al. 2021) that implement relativistic corrections to enable a reliable follow-up of the coalescence of compact objects. One of the mergers discussed here is actually the by-product of a series of four BH mergers that lead to the formation of a final merger with component masses 68 and 70 M_⊙, serendipitously compatible with the measured properties of GW190521. We use a statistical approach to discuss the impact of BH natal spins and recoil kick received upon GW emission in determining the likelihood for this type of merger to take place in young and massive clusters and similar environments.

The paper is organized as follows: In Section 2, we briefly describe the main features of the N-body models; In Section 3, we focus on the evolution of mergers beyond the upper-mass gap, i.e., involving IMBHs; Section 4 is devoted to mergers that contribute to the formation of BHs in the mass gap; In Section 5, we discuss the impact of GW kicks; In Section 6, we compare our results with similar works based on N-body models; Finally, Section 7 summarizes our main conclusions.

In the following, we analyze the database of 80 direct N-body simulations performed with an improved version of the NBODY6++GPU code (Wang et al. 2015), presented in Rizzuto et al. (2021). These models represent young massive star clusters evolved up to ∼300–500 Myr composed of 110,000 stars, including 10% of primordial binaries, with stellar masses in the range of (0.08–100) M_⊙, and typical central densities in the range of 10⁵–3 × 10⁷ M_⊙ pc⁻³. The key features of these models include stellar evolution of single and binary stars (Hurley et al. 2000, 2002; Belczynski et al. 2002), natal kicks for neutron stars (Hobbs et al. 2005), a fallback prescription to compute BH masses (Belczynski et al. 2002), a dedicated treatment for tight binary evolution (Mikkola & Aarseth 1998), and general relativistic corrections (Mikkola & Merritt 2008) to the equation of motions of compact stellar objects. ⁸

Although our models do not include PISNe and PPISNe, we note that the initial mass function (IMF) adopted in our models, in the range of 0.08–100 M_⊙ should be not severely affected by such mechanisms, as they should act on stars heavier than 100–120 M_⊙ for the metallicity value adopted here (see Spera & Mapelli 2017).

The stellar evolution prescriptions in these models is slightly outdated, but they represent among the first fully self-consistent N-body simulations in which all the relevant ingredients listed above are accounted for simultaneously. Similar works are either based on smaller (e.g., Di Carlo et al. 2019, 2020), or sparser clusters (e.g., Banerjee 2017, 2020). An improved version of these simulations that include state-of-the-art stellar evolution prescriptions (A. Kamlah et al. 2021, in preparation) is currently underway (M. Arca-Sedda et al. 2021, in preparation).

Among all the simulations we identify two single mergers involving one IMBH and a stellar BH, one merger in which one IMBH undergoes two subsequent mergers with stellar BHs, and one case in which a sequence of three consecutive BH mergers builds up an IMBH with total mass ∼140 M_⊙. In all the cases presented here, the binary merger is triggered by repeated interactions with passing-by stars and BHs, and in some cases, to the temporary formation of hierarchical triples that trigger an increase of the binary eccentricity and facilitate the merger (see Figures 6–7 in Rizzuto et al. 2021).

Hereafter, we label BH mergers consistently with the naming of runs adopted in Rizzuto et al. (2021). Model names consist of three parts: the first part identifies the initial value of the cluster half-mass–radius in parsecs, the second part identifies the value adopted for the central value of the adimensional potential well according to King (1966) models, whereas the third refers to the number of the simulation. We add a further label to identify the number of previous mergers that the BH underwent, e.g., Gen-1(Gen-2) refers to a first(second) generation merger. For the sake of clarity, we use letter a to identify the heavier component of the binary and b to identify the companion.

In the following, we refer to mergers beyond the mass gap as those involving at least one BH with a mass >150 M_⊙, whereas we define IMBHs all BHs with a mass >100 M_⊙. Similarly, we refer to mergers falling in the mass gap as those involving at least one component in the mass range of 65–120 M_⊙.

In the next section, we study separately mergers with BHs beyond or inside the upper-mass gap.

For all the mergers we record the semimajor axis and eccentricity the first time that the pericenter of the BH binary orbit falls below 100 times the sum of BH Schwarzschild radii. We follow the final inspiral phase integrating the Peters (1964) equations:

$$\begin{eqnarray}&&\displaystyle \frac{{da}}{{dt}}=-\displaystyle \frac{64}{5}\displaystyle \frac{{G}^{3}}{{c}^{5}}\beta ({M}_{1},{M}_{2})\displaystyle \frac{F(e)}{{a}^{3}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=-\displaystyle \frac{304}{15}\displaystyle \frac{{G}^{3}}{{c}^{5}}\beta ({M}_{1},{M}_{2})\displaystyle \frac{{eG}(e)}{{a}^{4}},\end{eqnarray} \tag{ 2 }$$

with

$$\begin{eqnarray}&&F(e)={\left(1-{e}^{2}\right)}^{-7/2}\left(1+\displaystyle \frac{73}{24}{e}^{2}+\displaystyle \frac{37}{96}{e}_{f}^{4}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\beta ({M}_{1},{M}_{2})={M}_{1}{M}_{2}({M}_{1}+{M}_{2}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&G(e)={\left(1-{e}^{2}\right)}^{-5/2}\left(1+\displaystyle \frac{121}{304}{e}^{2}\right).\end{eqnarray} \tag{ 5 }$$

This procedure enables us to determine the orbital properties of the binary during the last stages prior to the merger, a phase that is not captured from the direct N-body simulations. Along with the orbital evolution we calculate the associated GW strain and frequency, which are key ingredients in determining whether these mergers can be seen with GW detectors.

As detailed in Appendix B, we exploit a numerical relativity fitting formula to estimate the remnant BH mass, spin, effective spin parameter (Jiménez-Forteza et al. 2017), and the post-merger GW recoil kick (Campanelli et al. 2007; Lousto & Zlochower 2008; Lousto et al. 2012).

Varying the distribution of stellar BH natal spins and for IMBHs, we place constraints on (i) the spin distribution of remnants, and (ii) the likelihood for a merger remnant to be retained in the parent cluster and possibly undergo a further merger.

In the following, we study separately the mergers involving BHs beyond the mass gap and the one case falling right inside the mass gap, calculating for each merger the time evolution of the binary semimajor axis and eccentricity, the associated GW strain, the mass, spin, and effective spin parameter of the BH remnant, and the GW recoil kick.

Among the 80 simulations, we find a handful of interesting examples of BH binary mergers that led to a final BH with a total mass >100 M_⊙, and thus in the IMBHs' mass range. In two models, R06W9sim3 and R06W6sim6, the primary BH reaches a mass M_a = 285–328 M_⊙ through the accretion of a massive star onto a stellar mass BH, and after 11.3, and 113 Myr, respectively, form a binary with a stellar BH that undergoes coalescence in a few years. In another model (R06W9sim7), the primary BH grows up to a mass M_a = 307 M_⊙ and undergoes two subsequent merger events, first with a BH with mass M_b = 22 M_⊙ (R06W9sim7Gen-1) and later with another BH with mass M_b = 26 M_⊙ (R06W9sim7Gen-2).

In these cases, the heavier BH in the binary forms from the collision between a stellar mass BH and a very massive star (VMS) that assembled via a series of collisions in the cluster core among massive main-sequence stars (mass >50 M_⊙) over a short timescale, <10 Myr.

The amount of mass that is actually accreted onto a BH during a BH-star collision or accretion event is highly uncertain (see, e.g., Shiokawa et al. 2015; Metzger & Stone 2016; Law-Smith et al. 2019). In Rizzuto et al. (2021), this quantity is regulated through an accretion parameter f_acc that represents the fraction of the stellar mass fed to the BH. The whole set of simulations explores the impact of f_acc = 0.1–0.5–1.0 values. ⁹ The actual mass of an IMBH candidate in this simulation is thus given by M_IMBH = M_BH,pro + f_acc M_VMS, where M_BH,pro and M_VMS are the masses of the accretor stellar BH and the donor VMS, respectively. All the mergers discussed in this work take place in simulations with f_acc = 1, thus in which the VMS-BH accretion is maximized.

Whether such high-efficiency accretion can be attained is unclear. In the majority of our models, typically the VMS-BH binary tightens via three-body interactions with fly-by stars up to the point at which the pericentral distance falls below the VMS typical radius. The BH will thus penetrate the outer stellar layer and rapidly sink to its center, accreting mass in its course and after settling in the star center.

Many factors can contribute to determine the actual amount of mass accreted onto a compact object, like the type of the accretor (a neutron star or a BH), the evolutionary stage of the star, the size of the VMS, the structure of the VMS envelope, and the inclusion or not of general relativistic effects (see, e.g., MacLeod & Ramirez-Ruiz 2015; Cruz-Osorio & Rezzolla 2020).

Simulations of the common envelope phase in tight binaries composed of a BH and a more massive star show that the accretion of the core of the companion onto the BH and the expulsion of the envelope can trigger explosive phenomena that resemble SNe-like transients (Schrøder et al. 2020). A VMS formed out of repeated collisions among main-sequence stars is expected to have a compact core containing almost all the VMS mass and a tenuous envelope with densities as low as 10⁻¹⁰ g cm⁻³ (Glebbeek et al. 2009).

Therefore, a value of f_acc = 1 seems reasonable under the assumption that most of the VMS mass is concentrated in its core and that the core is completely accreted onto the BH.

In such a picture, the accretion process would spin up the BH up to nearly extremal values regardless of the initial BH spin (Schrøder et al. 2020), thus suggesting that IMBHs formed out of BH-VMS interactions could be characterized by spins close to unity.

As summarized in Table 1, at formation the semimajor axis (a) and eccentricity (e) of our IMBH-BH binaries span a wide range, with two of the cases exhibiting an extremely large eccentricity of e = 0.995–0.997. Whether the large eccentricity is preserved or not when the binary enters the frequency window accessible to GW detectors depends mainly on the binary orbital properties. The associated GW signal emitted by eccentric sources is characterized by a broad frequency spectrum, with the frequency of the dominant harmonic given by O'Leary et al. (2009),

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{GW}} & = & \displaystyle \frac{1.2\,\mathrm{mHz}}{1+z}{\left(\displaystyle \frac{{M}_{a}+{M}_{b}}{300\,{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{a}{2\,{R}_{\odot }}\right)}^{-3/2}\\ & & \times \,\mathrm{ceil}\left[1.15\displaystyle \frac{{\left(1+e\right)}^{1/2}}{{\left(1-e\right)}^{3/2}}\right]{\left(1-e\right)}^{-3/2}.\end{array}\end{eqnarray} \tag{ 6 }$$

Note that around 90% of the GW emitted power is associated with harmonics at frequencies of 0.2 < f/f_GW < 3 (O'Leary et al. 2009; Kocsis et al. 2012).

At formation, the dominant frequency for all the IMBH-BH mergers described above spans the f_GW ∼ 10⁻³–10⁻² Hz frequency range, thus falling inside the LISA sensitive window.

Figure 1 shows the evolution of the GW strain, the frequency, and eccentricity for all IMBH-BH mergers in our model, assuming that they are located at a redshift z = 0.05. As shown in the figure, at such distance the GW strain of these IMBH-BH mergers will move inside the LISA sensitivity window from 10³ yr to ∼1 day prior to the merger.

Zoom In Zoom Out Reset image size

Within this time span, the IMBH-BH eccentricity decreases from e ≃ 0.5 to e < 0.01 when the time to merger is 4 yr, and reaches values below e < 10⁻⁴ 1 s before the merger.

To characterize the typical eccentricity of the binary as it sweeps through the LISA frequency window, we proceed in two ways. On the one hand, we calculate the average eccentricity weighted with the time to the merger, which provides a measure of the probability of observing the binary in a given orbital configuration. On the other hand, we measure the binary eccentricity 10 and 5 yr prior to the merger, which is a period of time during which LISA can likely follow the binary inspiral.

As summarized in Table 1, we infer an average eccentricity of 〈e〉 ∼ 0.302 for the model with (M₁ + M₂) = (328 + 21) M_⊙, 〈e〉 ∼ 0.288 for the model with (M₁ + M₂) = (329 + 26) M_⊙, and 〈e〉 ∼ 0.296 for the model with (M₁ + M₂) = (285 + 22) M_⊙. The eccentricity 10 and 5 yr prior to the merger, instead, ranges in between e_5,10 = 0.04 and 0.13. Even measuring such a relatively small e values would provide a unique proof of a dynamical origin for this class of mergers, as the typical eccentricities for mergers in field binaries falling in the LISA band are expected to be much smaller, i.e., 10⁻⁶ < e < 10⁻⁴ (Kowalska et al. 2011; Nishizawa et al. 2016).

Detecting BH mergers with such masses using LISA depends on several parameters, like the redshift at which the merger takes place, the S/N required to obtain a clear GW signal, and the observation time. Assuming that the mergers occur closer than ∼44 Mpc, corresponding to a redshift z = 0.01, and adopting an observation time of 4 yr, we find that LISA could detect the last phase of these IMBH-BH mergers with S/N = 100–150, well above the minimum threshold required for detection, i.e., ${{\rm{S}}/{\rm{N}}}_{\min }=15$ (Amaro-Seoane et al. 2017). More in general, we derive the LISA horizon redshift, i.e., the maximum distance that LISA can reach for IMBHs with masses in the range of 100–400 M_⊙ merging with stellar BHs (mass 5–50 M_⊙), as detailed in Appendix A, in the limit of nearly circular sources. This choice is motivated by the fact that, as mentioned above, our IMBH-BH mergers are expected to be circularized in the LISA band. Assuming a minimum ${({\rm{S}}/{\rm{N}})}_{\min }=15$ and a 4 yr observation time, we find that LISA could detect these mergers out to redshift z = 0.01–0.1, depending on the masses of the binary components. For the sake of clarity, the GW strain-frequency evolution shown in Figure 1 assumes that the sources are located at redshift z = 0.05, well within the boundaries imposed by the horizon redshift. We note that taking into account the small residual eccentricity in the calculation of the S/N affects the horizon redshift by less than 10%. Our analysis suggests that mergers occurring in the Milky Way (MW) and its surrounding could be extremely bright low-frequency GW sources, making LISA a unique tool to unravel IMBH formation routes in young dense clusters.

Other important information that can be extracted from the GW signal is the spin of the merging BHs, which in turn provides insights into the BH formation mechanisms. The current picture about the natal spin distribution of stellar BHs is poorly constrained, and the impact of binary evolution, stellar rotation, or supernova explosion on the natal BH spin is unclear (de Mink et al. 2013; Amaro-Seoane & Chen 2016; Belczynski et al. 2020; Qin et al. 2018, 2019). The detection of GWs helped in placing some constraints on the spin distribution of merging BH binaries (Abbott et al. 2019), although there are a number of potential correlations between the spin and the binary formation history that are difficult to disentangle (e.g., Arca-Sedda & Benacquista 2019; Arca-Sedda et al. 2020b; Bavera et al. 2020; Zevin et al. 2021). For IMBHs, the picture is even more complex, as the spin will be inevitably connected with their formation history. For instance, matter accretion onto a stellar BH from a much heavier companion could efficiently spin up the BH and lead to a final IMBH with a spin close to unity (Schrøder et al. 2020). On the other hand, it is unclear whether the accretion of material from a normal stellar companion can efficiently spin up the growing BH (Qin et al. 2018, 2019).

Therefore, we follow a conservative approach to obtain insights insights into the IMBH and BH spin distribution and evolution. For each IMBH-BH merger, we assume that either the IMBH is almost non-spinning (S_1a = 0.01), nearly extremal (S_1a = 0.99), or it follows the same spin distribution of stellar BHs. For stellar BH spins we test either a flat distribution in the range of S_1b = 0–1 ¹⁰ or a Gaussian distribution peaked at S_1b = 0.2–0.5–0.7 with a dispersion of σ_S = 0.2. Upon these hypotheses, we sample 2000 values of the stellar BH spin for each merger and each value of S_1a, calculating the remnant mass and spin via numerical relativity fitting formulae (Jiménez-Forteza et al. 2017). Figure 2 shows the spin distribution of the IMBH remnant for model R06W9sim7 after the first- and the second-generation merger, assuming a Gaussian peaked at s_1b = 0.5 for the stellar BH spin distribution.

Zoom In Zoom Out Reset image size

We find that the spin of the remnant depends critically on the pre-merging IMBH spin, whereas it is weakly dependent on the stellar BH spin distribution. The remnant IMBH spin distribution narrowly peaks around S_2a = 0.2, whereas for a nearly extremal IMBH the remnant spin distribution is widely distributed and exhibits two peaks at S_2a = 0.6 and 0.99. The second merger event changes the IMBH spin further, leading to a double peak distribution with peaks at 0.15 and 0.3 if the IMBH was originally non-spinning, while leading to a wide three peak distribution in the case of an initially nearly extremal IMBH. If the IMBH spin distribution follows the same distribution of stellar BHs, e.g., a Gaussian, we found that the distribution gets skewed toward smaller spin values. The distribution peak shifts toward smaller values at every subsequent merger event, moving from ${S}_{{1}_{a},\mathrm{peak}}=0.5$ for Gen-1a to ${S}_{{2}_{a},\mathrm{peak}}=0.4$ for Gen-2a, and ${S}_{{3}_{a},\mathrm{peak}}=0.3$ for Gen-3a. Even in a single merger event, our analysis outlines that the natal IMBH spin has a crucial impact on the remnant spin. Spin measurements in IMBH-BH mergers via GW detections thus have the potential to unravel IMBH natal spin distribution, and thus, to help constrain IMBH formation routes.

Multiple-generation mergers can of further insights into IMBH formation and evolution, but their development can be prevented by recoil kicks imparted to the merger remnant due to anisotropic GW emission, which can be as large as 10³ km s⁻¹ (Lousto & Zlochower 2008) and thus can kick the IMBH out of the cluster after the first merger (e.g., Holley-Bockelmann et al. 2008).

In clusters with sufficiently large escape velocities, repeated mergers can lead the IMBH spin to slowly decrease and attain final values <0.1–0.3, regardless of the IMBH natal spin (see, e.g., Arca-Sedda et al. 2021). Owing to this, the distribution of IMBH spins formed in high escape velocity clusters should contain little information about the IMBH natal spin, which instead will dominate the spin distribution of IMBHs developed in low escape velocity clusters.

Unfortunately, GW recoil is not implemented in our N-body models, so we resort to post-processing of the data to explore how GW recoil kick would impact the retention of the merger products in dense clusters and how this would impact the spin distribution of merged IMBHs, as detailed in Section 5.

In another model (R06W6sim1), we found a binary BH merger with components in the upper-mass gap and interestingly similar to GW190521. This merger was formed via a peculiar channel, namely, a series of three subsequent merger events, as sketched in Figure 3.

Zoom In Zoom Out Reset image size

The sequence of mergers that built up our GW190521-like merger starts with a first-generation merger of two stellar mass BHs with masses 28 M_⊙ (hereafter Gen-1a) and 17 M_⊙ (Gen-1b), after 39 Myr. The Gen-1a merger product, which we refer to as Gen-2a, captures a companion, Gen-2b, with mass 25 M_⊙ and merges 10 Myr later, leaving behind a product, Gen-3a, with a total mass of 70 M_⊙. Gen-3a forms a binary with Gen-3b, a BH with mass 68 M_⊙ formed out of the accretion of a red giant star with mass ∼50 M_⊙ onto a stellar BH with mass ∼30 M_⊙. The binary Gen-3ab undergoes coalescence within 84 Myr, assembling Gen-4, a fourth-generation BH with total mass 138 M_⊙. ¹¹ Note that Gen-3a(b) represents the primary(secondary) in the final merger of the sequence, corresponding to the naming of GW190521 components.

About 200 yr prior to the merger, the Gen-3ab binary had a semimajor axis of 864 solar radii and eccentricity 0.99982. However, under the assumption that the system is isolated, GW emission induces the circularization of the binary and leads its eccentricity to drop below 10⁻³ by the time the binary enters the LIGO sensitivity band, at frequencies >10 Hz. Measuring the eccentricity for GW sources in the LIGO band is quite difficult, as the merger spends a short time in the band. So far, the analysis of LIGO-Virgo data seems to disfavor a significant eccentricity at the merger (LIGO Scientific Collaboration 2020), although a more robust assessment of the binary eccentricity is complicated by a partial degeneracy between high eccentricity and precession (Romero-Shaw et al. 2020; LIGO Scientific Collaboration 2020).

Assuming a redshift z = 0.8, i.e., the value estimated for GW190521, we found that LISA would not be able to detect any of the three mergers. When the last merger enters the LIGO sensitive band, i.e., the frequency range f = 10–1000 Hz, the binary eccentricity is already e < 10⁻⁴. This is due to the fact that at the last snapshot our binary has a relatively wide orbit, with a semimajor axis of a ∼ 4 au. For instance, reducing the a by a factor of 10(50) would lead the eccentricity measured at 10 Hz to increase up to 0.02(0.06), yet smaller than the maximum eccentricity inferred for GW190521 (Romero-Shaw et al. 2020). An eccentricity above 0.1 could be reached if the a = 0.4 au and the corresponding eccentricity is e = 0.99995. Note that the typical semimajor axis of long-lived (hard) binaries in star clusters scales with the cluster velocity dispersion as a_hard ∝ σ², thus statistically speaking, forming binaries 10 times tighter would require a host cluster with a velocity dispersion $\simeq \sqrt{10}$ times larger.

Following the same scheme adopted for IMBH mergers, we calculate how the remnant mass and spin vary after every new merger event and upon different assumptions of the BH natal spin distribution. For the BH natal spin distribution we adopt the same assumption as in the previous section. We assume that the natal spin of Gen-3b is not affected substantially by the accretion from the previous merger with a red giant companion (Valsecchi et al. 2010; Wong et al. 2012; Qin et al. 2019).

As shown in Figure 4, we find that a chain of three mergers, such as the one described here, leads to final BH spin and effective spin parameter χ_eff distributions that are fully compatible with the 90% credible level measured for GW190521. In Section 4.1, we compare the distribution for S and χ_eff assuming that GW190521 formed through one, two, or three consecutive mergers. We show that a triple merger provides a better match to observations than a double merger, depending on the subtle BH natal spin distribution.

Zoom In Zoom Out Reset image size

However, the successful development of a multiple-generation merger depends critically on the ratio between two quantities, namely, the cluster escape velocity and the GW recoil kick imparted onto the remnant BH after each merger event. In the next section we quantify the probability for such merger chain to occur in our simulated cluster or in denser environments.

4.1. Degeneracy in Multiple-merger Scenarios Behind the GW190521 Origin

Several works in the recent literature have pointed out that GW190521 could be the by-product of a hierarchical merger, i.e., a series of mergers that built up the GW190521 primary first and later gave rise to the detected merger. In our simulations we find one such case in which a series of three mergers led to a final merger similar to GW190521.

In this section, we compare a triple-merger channel for the origin of GW190521 with simpler scenarios, namely, a single merger, although this is challenged by stellar evolution theory of both single and binary stars, and a double merger. Figure 5 compares the remnant spin and effective spin distribution for these three channels (single, double, and triple mergers) assuming that the natal BH spin distribution is either Gaussian with a peak at 0.5 or uniform. Regardless of the BH natal spins, the plot highlights that a double-merger scenario seems less favored over a single or a triple-merger scenarios, as they both produce a remnant spin distribution that nicely fits the 90% credible level measured for GW190521. Quantitatively speaking, the fraction of merger remnants with a remnant spin compatible with GW190521 ranges between ∼0.66 and 0.72 for the single and triple-merger channels, while it is limited to 0.58–0.6 for the single channel. Similarly, a single or triple merger is more likely to produce a remnant with χ_eff compatible with that measured for GW190521. Our analysis suggests that measuring spin and effective spin parameters of massive BBH mergers can help unravel the BH natal spin distribution and constrain further the BH natal mass distribution. For instance, the single merger channel would hint at a scenario in which the GW190521 components grew their mass via either stellar accretion or stellar collision, or would imply that some processes in single and binary stellar evolution are yet not fully understood. As we show in the next section, the triple-merger channel would instead have crucial implications on the properties of the host cluster. Table 2 summarizes the fraction of models compatible with GW190521 remnant spin for the different models explored.

Zoom In Zoom Out Reset image size

Table 2. Occurrence of GW190521-like Merger Remnant Spin and Effective Spin Parameter for Single (First Sub-column in Main Columns 2 and 3), Double (Second Sub-column), and Triple (Third Sub-column) Mergers

BH Spin	Srem			χeff
Distribution	1	2	3	1	2	3
Uniform	0.716	0.579	0.666	0.679	0.587	0.624
Gaussian 0.2	0.929	0.618	0.782	0.936	0.616	0.707
Gaussian 0.5	0.716	0.579	0.666	0.691	0.569	0.614
Gaussian 0.7	0.582	0.582	0.597	0.582	0.558	0.581

Download table as:  ASCIITypeset image

In this section we calculate the GW recoil kick imparted onto our IMBH and BH mergers and the implication for IMBH retention and the development of multiple mergers in YMCs. As detailed in Appendix A, for each merger we calculate the GW recoil via numerical relativity fitting formulae (Campanelli et al. 2007; Lousto & Zlochower 2008; Lousto et al. 2012) and compare it with the cluster velocity dispersion.

5.1. Mergers Beyond the Mass Gap

We focus on IMBH-BH mergers first, particularly on the double merger found in model R06W6sim7. Figure 6 shows the kick distribution imparted on the post-merged IMBH in the case in which the IMBH is initially either nonrotating, nearly extremal, or has a spin drawn from the same distribution adopted for BHs. For stellar BHs, the spin distribution is either a Gaussian peaked on 0.5, 0.7, or a uniform distribution. If the IMBH is initially slowly spinning, the kick received after both the first and second merger attains values around v_kick = 30 km s⁻¹, and in general smaller than 80 km s⁻¹. Larger IMBH spins would instead lead the GW kick to increase considerably, reaching values up to 500 km s⁻¹. Nonetheless, even for spinning IMBHs the GW kick distribution shows a peak at 30–40 km s⁻¹. Our simulated cluster has a total mass of M_c = 7.4 × 10⁴ M_⊙ and half-mass–radius of r_h = 1 pc. If we assume that the central part of such cluster can be described by a power law with slope γ, it is possible to express the central escape velocity as (Dehnen 1993)

$$\begin{eqnarray}&&{v}_{\mathrm{esc}}^{2}=\displaystyle \frac{2{{GM}}_{c}}{(2-\gamma ){r}_{h}}\left({2}^{1/(3-\gamma )}-1\right),\end{eqnarray} \tag{ 7 }$$

which returns v_esc ≃ 45.8–51.6 km s⁻¹ assuming γ = 0–1. We note that the estimate obtained via the equation above is around twice the escape velocity value calculated through the gravitational potential calculated directly from the simulation.

Zoom In Zoom Out Reset image size

In the case of a steeper density profile, the escape velocity for a Dehnen (1993) model diverges in the cluster center. If we assume that the merger takes place in the inner part of the cluster, e.g., r < 0.1r_h , the escape velocity rises to up to 60–112 km s⁻¹ for γ = 1.5–2.

Assuming v_esc = 40 km s⁻¹, i.e., a value compatible with the one measured from the N-body simulation, we find that the remnant of a first-generation IMBH-BH merger has a retention probability of P_Gen−2a = 96.9%–99.1% if the IMBH initial spin is almost null, while the second-generation merger remnant is retained with a P_Gen−3a = 1.8%–2.1% probability. This probability rises up to P_Gen−3a = 50% assuming v_esc = 50 km s⁻¹. If the initial IMBH spin is ∼1, instead, the retention probability decreases down to P_Gen−2a = 3(0.1)% for the first(second)-generation merger product

Therefore, measuring the spin in an IMBH-BH merger could provide crucial insights into the possible host environment in which the merger developed.

5.2. Mergers in the Mass Gap

To place constraints on the GW recoil kick imparted to BH remnants in the triple-merger channel that leads to the GW190521-like merger in one of our simulations, we assign to the merger components a spin extracted from either a Gaussian, centered on either S = 0.2–0.5–0.7, or a uniform distribution. For each merger event, we sample 10,000 values of the spin for each component and calculate the GW recoil imparted on the remnant. We find that in all three mergers the remnant has more than a 90% probability of receiving a kick larger than 80 km s⁻¹, regardless of the BH natal spin distribution.

From Equation (7) we thus can conclude that the retention of remnants receiving such large kicks is maximized in clusters with a steep density distribution. Using the values that characterize the N-body model in which the triple merger develops, i.e., a cluster mass of 7.4 × 10⁴ M_⊙ and a half-mass–radius of 0.6 pc, we compute the Gen-2a retention probability, namely, the percentage of models for which v_kick < v_esc, at varying the cluster density slope γ. Figure 7 shows how the retention probability for the Gen-2a merger varies with increasing the slope of the density profile. Our analysis suggests that a slope larger than γ > 1.8–2 would ensure a retention probability >0.3%–80%, depending on the underlying BH natal spin distribution.

Zoom In Zoom Out Reset image size

This implies that is nearly impossible to develop a long chain of mergers in a YMC similar to the simulated one unless the cluster is characterized by a quite steep density profile. In principle, in order for the R6W6sim1 triple merger to develop, our analysis suggests that the cluster escape velocity should increase by a factor of at least 3–5, a threshold that could be achieved with a heavier or denser cluster. However, the development of the triple-merger chain is intrinsically locked to the evolution of the cluster. For instance, the mass-segregation timescale regulates the time over which the BHs reach the cluster core and start interacting with each other. In our N-body models the segregation of the most massive stars takes place within the first 1–2 Myr. Since both the collision rate and the mass-segregation timescale are tightly linked to the cluster relaxation time (Portegies Zwart & McMillan 2002), the minimum requirement to attempt an extrapolation from our models to heavier and denser clusters is that the mass-segregation time, which scales with the cluster relaxation time ${t}_{\mathrm{rel}}\propto {M}_{c}^{1/2}{r}_{h}^{3/2}$ (Binney & Tremaine 2008), remains constant. Note that this roughly corresponds to the requirement that the mass-segregation time in the cluster remains constant as well (see, e.g., Arca-Sedda 2016). At the same time, we need that the cluster escape velocity, which for Dehnen models scales with the M_c /r_h ratio as shown in Equation (7) above, increases above the threshold imposed by v_kick. Thus, by requiring that the relaxation time is the same as in the N-body model and that v_esc > v_kick enables us to find the (M_c , r_h ) values that describe heavier clusters that should be equivalent to our N-body models from a dynamical point of view. Clearly, this is an oversimplification that neglects the impact of several processes, e.g., a substantial primordial mass segregation or fraction of binaries. Figure 8 shows the locus of points with constant relaxation time compared to the typical mass and half-mass–radius of NCs in the local universe (taken from Georgiev et al. 2016), MW GCs (Harris 2010), YMCs in the MW and MW satellites (Portegies Zwart et al. 2010), and the 11 YMCs detected in the starburst galaxy Henize 2–10 (Nguyen et al. 2014; Arca-Sedda et al. 2015), often referred to as super star clusters (SSCs).

From the plot it is apparent that only clusters with a mass M_c > 3 × 10⁵ M_⊙ and half-mass–radius r_h < 0.6 pc have a central escape velocity v_esc ≳ 90 km s⁻¹ and thus could retain the remnant of the first merger in our triple-merger channel. The comparison with observed clusters suggest that a fraction of NCs (∼5% in the Georgiev et al. 2016 sample of low redshift NCs) and SSCs similar to the one found in Henize 2–10 are sufficiently dense to fulfill these two conditions and potentially harbor a multiple-merger chain. The fact that the triple-merger chain occurs very early in our models (<100 Myr) seems to favor the SSC scenario, as this class of clusters is usually very young—the inferred age of the SSCs in Henize 2–10 is ∼5 Myr—massive and dense (e.g., Nguyen et al. 2014). Alternatively, a young NC formed out of SSC collisions (Arca-Sedda et al. 2015) could represent a potential host for our triple-merger chain.

Note that among all the models explored, we find that the kick imparted to the first merger remnant exceeds 100 km s⁻¹ with a probability of >90%. This would suggest that even a double-merger scenario seems unlikely for GW190521, unless the host cluster was sufficiently dense and massive as shown in Figure 8.

Zoom In Zoom Out Reset image size

5.3. Merger Rates

Although our analysis is based on a handful models, the few IMBH-BH mergers developed in our simulations can be used to infer a rough estimate of the merger rate for this type of GW source. In the following, we infer the merger rate for first-generation IMBH-BH mergers. Note that the development of these mergers arise from dynamics and is not affected at all by the lack of a treatment of GW kicks and spins.

One possible way to calculate the merger rate is to assume that the formation of young clusters proceeds at a fraction α < 1 of the cosmic star formation rate (SFR) ψ(z), which can be expressed as (Madau & Fragos 2017)

$$\begin{eqnarray}&&\psi (z)=0.01\displaystyle \frac{{\left(1.+z\right)}^{2.6}}{\left[1.+{\left(\left(1.+z\right)/3.2\right)}^{6.2}\right]}\,\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}.\end{eqnarray} \tag{ 8 }$$

In the following, we adopt α = 0.01 (e.g., Belczynski et al. 2018).

From our models, we can define a merger efficiency as the ratio between the number of merger N_mer = 3 and the total number of simulations ${N}_{\mathrm{sim}}=80$, divided by the cluster mass M_c = 7.4 × 10⁴ M_⊙.

The fraction of star clusters with a mass similar to the simulated one within a given _M range can be obtained as:

$$\begin{eqnarray}&&{\eta }_{M}=\displaystyle \frac{{\int }_{{M}_{1}}^{{M}_{2}}f(M){dM}}{{\int }_{{M}_{\min }}^{{M}_{\max }}f(M){dM}}=\displaystyle \frac{{M}_{2}^{1-s}-{M}_{1}^{1-s}}{{M}_{\max }^{1-s}-{M}_{\min }^{1-s}},\end{eqnarray} \tag{ 9 }$$

where M_1,2 = M_c (1 ∓ _M ) and ${M}_{\min ,\max }={10}^{4-7}\,{M}_{\odot }$. In the latter equality of the equation above, we assume that the cluster mass function is well described by a power law f(M) ∝ M^−s with slope s = 2 (e.g., Gieles 2009). If we assume an _M = 0.3 value, thus implying that only clusters with a mass deviating less than ∼30% from the simulated one can develop the IMBH-BH mergers discussed here, we get η_M ≃ 0.09.

Combining together the considerations above we can infer the merger rate as a function of redshift as

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Gamma }}}{{dz}}=\displaystyle \frac{{N}_{\mathrm{mer}}}{{N}_{\mathrm{sim}}{M}_{c}}{\eta }_{M}\alpha \psi (z),\end{eqnarray} \tag{ 10 }$$

whose behavior as a function of redshift is shown in the top panel of Figure 9. We note that the peak at redshift z = 2 is intrinsically inherited from the adopted SFR.

Zoom In Zoom Out Reset image size

The overall merger rate within a redshift z could be thus obtained by integrating the equation above over the cosmological volume,

$$\begin{eqnarray}&&{\rm{\Gamma }}(\lt z)={\int }_{0}^{{z}_{\max }}\displaystyle \frac{d{\rm{\Gamma }}}{{dz}}\displaystyle \frac{{dV}}{{dz}}\displaystyle \frac{1}{1+z}{dz},\end{eqnarray} \tag{ 11 }$$

where T_obs is the mission duration time, dV/dz is the comoving cosmological volume element, and (1 + z)⁻¹ is a factor accounting for the dilation time. The bottom panel of Figure 9 shows the cumulative merger rate Γ( < z) as a function of the redshift. Given the horizon redshift for LISA, and assuming an observation time of T_obs = 4 yr and a minimum S/N of 15, we infer a total number of detections ΓT_obs = 10⁻³–10⁻¹.

We note that in order to reach a detection rate of Γ = 0.2–2.5 yr⁻¹ with the same observation time, the LISA horizon should expand beyond redshift z = 0.3–1. Reaching such redshift would require an increase in the LISA sensitivity curve, especially in the 0.01–1 Hz frequency range.

Assuming the most recent sensitivity curve (Amaro-Seoane et al. 2017; Robson et al. 2019), we find that increasing LISA's sensitivity by a factor 2 would permit detecting IMBH-BH mergers out to a redshift z = 0.2, thus implying Γ ∼ 0.2 yr⁻¹.

While detecting these sources with LISA might be tricky, future detectors sensitive to decihertz GWs could clearly pitch their signals. For instance, in the same range of masses of our interest, DECIGO could reach a redshift z > 10², thus enabling the possibility of probing IMBH-BH mergers up to the time of the formation of the first stars. A mid-range detector could fill the gap between LISA and DECIGO. For instance, adopting the sensitivity curve proposed for the conservative decihertz observatory concept design discussed in Arca-Sedda et al. (2020a), we find a horizon redshift z > 4, leading to a detection rate Γ ∼ 10–100 yr⁻¹.

The suite of direct N-body simulations of young dense clusters in which we find the GW190521-like merger is carried out with an improved version of the NBODY6++GPU direct N-body code (Wang et al. 2015). The key features of these models include stellar evolution of single and binary stars (Hurley et al. 2000, 2002), natal kicks for neutron stars (Hobbs et al. 2005), a fallback prescription to compute BH masses (Belczynski et al. 2002), a dedicated treatment for tight binary evolution (Mikkola & Aarseth 1998), and general relativistic corrections to the equations of motion of compact stellar objects (Rizzuto et al. 2021). The current implementation of stellar evolution does not include the latest recipes for supernovae mechanisms (e.g., PISN and PPISN).

We note that the mass function adopted in this work, which is truncated at ${m}_{\max }=100\,{M}_{\odot }$, and the metallicity value assumed (Z = 0.0002), make PISN and PPISN mechanisms inefficient in our models, as they should mostly affect stars with initial masses in the range of 100–200 M_⊙ (e.g., see Figure 2 in Spera & Mapelli 2017).

In this regard, it should also be noted that the actual extension of the upper-mass gap is largely unknown and a number of recent stellar evolution models have shown that it is even possible to populate the mass gap with BHs formed through ordinary stellar evolution (see, e.g., Vink et al. 2021; Woosley & Heger 2021).

Although a handful stars in some of our models can exceed the 100 M_⊙ mass limit imposed by the initial mass function, e.g., through stellar mergers, it must be noted that the remnant BH mass is limited to 30–40 M_⊙, owing to the adopted stellar evolution recipes (see Figure 2 in Rizzuto et al. 2021). The few BHs formed this way will have masses overlapping the overall BH mass spectrum, and given their small number (order of a few in 110,000 stars) are not expected to affect significantly the overall cluster evolution.

Nonetheless, these are among the first N-body models that take into account simultaneously both general relativistic effects, stellar evolution for single and binary stars, and a substantial fraction of primordial binaries. The creation of a new simulation database is currently underway and will feature up-to-date stellar evolution prescriptions, a larger number of stars, and a more abundant population of primordial binaries.

Our models simulate the evolution of young and dense star clusters with typical central densities in the range 10⁵–3 × 10⁷ M_⊙ pc⁻³). ¹²

In comparison to other works in the literature (e.g., Banerjee 2017, 2020; Di Carlo et al. 2019, 2020) our models represent on average heavier clusters, 2–10 times compared to Di Carlo et al. (2019, 2020) and 1–10 times compared to Banerjee (2017), and denser, since our clusters have a half-mass–radius of 0.6–1 pc.

Moreover, Rizzuto et al. (2021) explores the impact of the primordial binary fraction adopting a value f_b = 0.1, thus corresponding to the maximum value explored by Banerjee (2017) and definitely being smaller than the value used in Di Carlo et al. (2019), who assumed f_b = 40%.

The models presented here treat stellar feeding onto BHs through a parameter, f_acc, which regulates the amount of mass that is accreted onto the BH during an accretion or collision event. Rizzuto et al. (2021) explores the range f_acc = 0.1, 0.5, 1.0. In comparison to our models, Di Carlo et al. (2019) adopted the more conservative assumption that the stellar material is completely lost and the BH does not accrete at all, whereas Banerjee (2020) explored the range f_acc = 0.7, 0.9 in their latest models.

Another element of difference is in the maximum mass adopted for the IMF: Rizzuto et al. (2021) adopts a ${m}_{\max }=100\,{M}_{\odot }$ value, whereas Banerjee (2017) uses a cluster mass dependent relation (Kroupa et al. 2013) that leads to ${m}_{\max }\lt 100\,{M}_{\odot }$ in the cluster mass range discussed here, and Di Carlo et al. (2019) assumes ${m}_{\max }=150\,{M}_{\odot }$.

Finally, we note that only Banerjee (2020) took into account GW recoil kick self-consistently in the simulations in the recent literature of direct N-body models. However, the models discussed in Banerjee (2020) cover a different portion of the phase space in terms of cluster half-mass–radius, density profile, metallicity, and initial binary fraction. Out of the 65 models described in Banerjee (2020), only one has structural properties compatible with Rizzuto et al. (2021) models, e.g., a similar cluster mass and metallicity, although crucial differences remain, i.e., a different metallicity, fraction of binaries, and stellar evolution recipes.

We also note that the approach presented in this paper is one of the first adopting a post-processing procedure to quantify the impact of such mechanism (see also Anagnostou et al. 2020b). This represents a fast technique to statistically assess the impact of GW kicks and spins without the need of running the exact same simulation many times, and thus saving considerable computational and human time. In these regards, it should be noted that, despite the limitations of our N-body models, all the first-generation mergers presented in this paper involve a stellar BH that never experienced previous mergers and an IMBH formed out of stellar accretion onto a stellar BH, and thus they can be considered a reliable outcome of stellar dynamics and evolution.

The discussion above highlights how our models occupy a portion of the phase space not fully covered by previous works, providing a substantially new exploration of YMC dynamics and evolution.

We analyzed the outputs from a suite of 80 N-body models of YMCs with N = 110,000 stars, 10% of which initially paired in binaries. Our N-body simulations are among the first in which the number of stars is >10⁵ and the motion of compact remnants takes into account general relativistic corrections, thus ensuring a reliable description of possible merger events. Our main results can be summarized as follows:

• 1.  
  We found three mergers out of 80 models involving a BH beyond the upper-mass gap and one model in which a merger similar to GW190521 develops through a triple-generation merger chain, all taking place over timescales t ∼ 10–300 Myr, at an early stage of the host cluster evolution;
• 2.  
  In mergers labeled as beyond the mass gap, the primary is an IMBH with a mass in the range of M₁ = 205–329 M_⊙, whereas the companion is a stellar BH with M₂ ≃ 20 M_⊙. In one of these cases, the IMBH undergoes two subsequent mergers. These mergers fall into the class of intermediate-mass ratio inspiral GW sources;
• 3.  
  We show that this type of merger can be detected with LISA up to a redshift z = 0.01–0.1 assuming a mission duration of 4 yr and ${{\rm{S}}/{\rm{N}}}_{\min }=15;$
• 4.  
  When the IMBH-BH binary decouples from the cluster dynamics the eccentricity can be very large, e = 0.085–0.997. However, 10(5) yr prior to the merger it reduces to e = 0.041–0.128(0.031–0.098), due to angular momentum loss via GW emission. This is likely due to the large separation at formation of the IMBH-BH mergers, which in turn is due to the relatively low cluster mass. Despite being small, measuring such eccentricity values with LISA would uniquely probe the dynamical channel, since isolated binaries are expected to have 10⁻⁶ < e < 10⁻⁴ in the LISA band;
• 5.  
  We find that the IMBH-BH remnant offers insights into the IMBH natal spin. The spin distribution for an IMBH remnant in an IMBH-BH merger is well defined: if the IMBH is initially nonrotating, the remnant spin distribution would be narrowly peaked around S_rem = 0.2, whereas for nearly extremal IMBHs the potential remnant spin values are broadly distributed between 0.6 and 1. The different distribution is preserved upon a double merger, thus suggesting that detecting IMBH mergers beyond the mass gap is crucial to assessing the IMBH natal spin and thus IMBHs' formation processes;
• 6.  
  In one simulation out of 80, we found a triple-merger chain event that led to a merger between BHs likely in the upper-mass gap, M_1,2 = (70 + 68) M_⊙, quite similar to the LVC source GW190521. Regardless of the BH natal spin distribution adopted, we find that such a triple-merger chain leads to a remnant with a final spin that matches well the value inferred for GW190521;
• 7.  
  At formation, the mass-gap merger has an eccentricity of e = 0.897, reducing to e ≤ 3 × 10⁻⁴ when the binary enters the GW frequency window of 10 Hz (the LVC window);
• 8.  
  Comparing models and the detected final spin and effective spin parameter, we show that a triple-merger scenario has a larger probability (59.7%–78.2%) of producing a remnant with a spin and effective spin parameter matching GW190521 than a double-merger scenario (57.9%–61.8%). The probability is similar to the case in which GW190521 developed in a single merger event (58.2%–92.9%);
• 9.  
  We estimate the recoil kick imparted to the merger remnant in the mergers both beyond and inside the mass-gap. In the case of IMBH-BH mergers, we show that the remnant receives a recoil in the range of v_kick = 30–50 km s⁻¹ for Schwarzschild IMBHs, whereas v_kick is as large as 500 km s⁻¹ for Kerr IMBHs. Therefore, finding signatures of IMBHs in YMCs can tell us more about IMBH natal spins. Statistically speaking, detecting signatures of an IMBH in a YMC could imply three possibilities: (1) The IMBH natal spin is very small, thus the recoil kick is small, (2) the IMBH did not undergo mergers with other BHs, and (3) the birth site is an extremely dense YMC;
• 10.  
  We derive a rough estimate for the merger rate and detection rate of IMBH-BH mergers for LISA and mid-range observatories like DECIGO. In the case of LISA, we find a detection rate of Γ = 0.1 yr⁻¹, assuming a 4 yr long mission and a signal-to-noise threshold value of S/N = 15. We show that detectors sensitive to decihertz have the potential to shed light on IMBH formation mechanisms, reaching redshift of up to z > 10, with inferred detection rates of 10–100 yr⁻¹ for DECIGO-like detectors.
• 11.  
  Regarding the merger-gap merger, we show that the kick imparted to the merger remnant of the first event in the merger chain is generally larger than 70–100 km s⁻¹, thus making it unlikely for such a channel to develop in normal YMCs. Nonetheless, retention can be ensured at an 80% level if the YMC density distribution is particularly steep, i.e., with a slope γ > 1.7. We show that, under optimistic assumptions, sufficiently dense and young NCs could be the nursery of such a triple-merger channel.

The authors acknowledge the anonymous referees for their constructive reports. The authors are grateful to Ataru Tanikawa, Jorick Vink, and Chris Belczynski for useful comments and discussions. M.A.S. acknowledges support from the Alexander von Humboldt Foundation and the Federal Ministry for Education and Research for the research project "Black Holes at All Scales." R.S. acknowledges the Strategic Priority Research Program (Pilot B) "Multi-wavelength gravitational wave universe" of the Chinese Academy of Sciences (grant No. XDB23040100) and National Science Foundation of China grant No. 11673032. M.G. was partially supported by the Polish National Science Center (NCN) through grant No. UMO-2016/23/B/ST9/02732. This work benefited of support from the Volkswagen Foundation Trilateral Partnership project No. I/97778 "Dynamical Mechanisms of Accretion in Galactic Nuclei," the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 138713538—SFB 881 ("The Milky Way System"), and the COST Action CA16104.

To make predictions on the detectability of our IMBH-BH mergers by LISA, we need to infer the GW source horizon, which determines the maximum distance in space, or the redshift z_hor, at which the source signal is detected with a threshold S/N, namely,

$$\begin{eqnarray}&&{\left({\rm{S}}/{\rm{N}}\right)}^{2}={\int }_{{f}_{1}}^{{f}_{2}}\displaystyle \frac{{h}_{c}^{2}\left(f,{z}_{\mathrm{hor}}\right)}{{S}_{n}^{2}\left(f\right)}{df},\end{eqnarray} \tag{ A1 }$$

with f_1,2 the initial and final frequency of the GW signal, h_c (f, z) its characteristic strain, and S_n (f) is the detector sensitivity. To determine z_hor, we calculate f₁ and f₂ by integrating the final stage of the IMBH-BH merger signal assuming an observation time of 4 yr, i.e., the nominal duration time of the LISA mission, and adopting an S/N = 15. We assume the set of cosmological parameters measured by the Planck mission, namely, H₀ = 67.74 km s⁻¹ Mpc⁻³ and Ω_m = 0.3089, Ω_Λ = 0.6911 (Planck Collaboration 2016). Figure 10 shows how the horizon redshift changes for LISA ¹³ at varying the IMBH mass in the range of M_IBH = 100–400 M_⊙ and the stellar BH companion mass in the range of M_BH = 5–30 M_⊙.

Zoom In Zoom Out Reset image size

We find that circular IMBH-BH mergers similar to the one found in our simulations can be observed out to a redshift of z_hor = 0.01–0.15, corresponding to a luminosity distance of D_L ≃ 44–500 Mpc. We note however that this represents a lower limit of z_hor, as the GW power emitted by the higher order harmonics in the case of eccentric orbits sums up to the zero-order harmonic, making the source brighter, and thus, observable from slightly farther away.

The N-body simulations presented here do not implement a treatment for BH spins and GW recoil, and do not account for GW energy loss in the evaluation of the remnant mass, which is calculated as the mass of the merging BHs.

To include these important ingredients we post process the simulated data creating a set of mergers with BH masses given by the N-body model and spins drawn by different BH spin distributions.

In the case of IMBHs (mass > 200 M_⊙) we explore two limiting cases, whether the IMBH is slowly (S_IBH = 0.01) or highly (S_IBH = 0.99) spinning. For stellar BHs, the spin should be somehow connected with the BH formation process, but currently the picture around BH natal spins is poorly constrained. Given the lack of a general consensus around BH natal spin, we explore four different cases: either the BH natal spin distribution is a Gaussian peaked at either 0.2, 0.5, or 0.7, or it is uniform in the 0–1 range of values.

We create 10,000 versions of each merger event as listed in Table 1 and for each spin distribution adopted. The mass and spin of the remnant BH are calculated through numerical fitting formulae derived by Jiménez-Forteza et al. (2017). In the case of multiple-generation mergers we use the remnant mass and spin from the previous generation. For each merger we also calculate the chirp mass ${ \mathcal M }={({m}_{1}{m}_{2})}^{3/5}/{({m}_{1}+{m}_{2})}^{1/5}$ and the effective spin parameters

$$\begin{eqnarray}&&\chi =\displaystyle \frac{{M}_{a}{S}_{a}\cos {\theta }_{1}+{M}_{b}{S}_{b}\cos {\theta }_{2}}{{M}_{a}+{M}_{b}},\end{eqnarray} \tag{ B1 }$$

where $\cos {\theta }_{i}$ represents the angle between the ith component spin vector and the binary angular momentum. Note that these two quantities are particularly important for stellar BH mergers, as they are better constrained in high-frequency GW detectors than total mass and remnant spin.

The GW recoil velocity received by a remnant after a merger can be calculated as in Campanelli et al. (2007), Lousto & Zlochower (2008), and Lousto et al. (2012):

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{\mathrm{GW}}={v}_{m}{\hat{e}}_{\perp ,1}+{v}_{\perp }(\cos \xi {\hat{e}}_{\perp ,1}+\sin \xi {\hat{e}}_{\perp ,2})+{v}_{\parallel }{\hat{e}}_{\parallel },\end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray}&&{v}_{m}=A{\eta }^{2}\sqrt{1-4\eta }(1+B\eta ),\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}&&{v}_{\perp }=\displaystyle \frac{H{\eta }^{2}}{1+{q}_{\mathrm{BBH}}}\left({S}_{2,\parallel }-{q}_{\mathrm{BBH}}{S}_{1,\parallel }\right),\end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray}\begin{array}{rcl}{v}_{\parallel } & = & \displaystyle \frac{16{\eta }^{2}}{1+{q}_{\mathrm{BBH}}}\left[{V}_{11}+{V}_{A}{{\rm{\Xi }}}_{\parallel }+{V}_{B}{{\rm{\Xi }}}_{\parallel }^{2}+{V}_{C}{{\rm{\Xi }}}_{\parallel }^{3}\right]\\ & & \times \,\left|{{\boldsymbol{S}}}_{2,\perp }-{q}_{\mathrm{BBH}}{{\boldsymbol{S}}}_{1,\perp }\right|\cos ({\phi }_{{\rm{\Delta }}}-{\phi }_{1}).\end{array}\end{eqnarray} \tag{ B5 }$$

In the equation above, $\eta \equiv {q}_{\mathrm{BBH}}/{(1+{q}_{\mathrm{BBH}})}^{2}$ is the symmetric mass ratio, ${\boldsymbol{\Xi }}\equiv 2{({{\boldsymbol{S}}}_{2}+{q}_{\mathrm{BBH}}^{2}{{\boldsymbol{S}}}_{1})/(1+{q}_{\mathrm{BBH}})}^{2}$, and subscripts ⊥ and ∥ mark directions (perpendicular and parallel) to the BH spin vector with respect to the direction of the binary angular momentum. The orthonormal basis defined by the unit vectors (${\hat{e}}_{\parallel },{\hat{e}}_{\perp ,1},{\hat{e}}_{\perp ,2}$) has one component directed perpendicular to (${\hat{e}}_{\parallel }$) and two components lying in the binary orbital plane. We set A = 1.2 × 10⁴ km s⁻¹, B = − 0.93, H = 6.9 × 10³ km s⁻¹, and ξ = 145° (González et al. 2007; Lousto & Zlochower 2008), and V_A,B,C = (2.481, 1.793, 1.507) × 10³ km s⁻¹ are scaling constants (Lousto et al. 2012). ϕ_Δ represents the angle between the direction of the infall at merger (which we randomly draw in the binary orbital plane) and the in-plane component of ${\boldsymbol{\Delta }}\,\equiv {({M}_{a}+{M}_{b})}^{2}({{\boldsymbol{S}}}_{b}-{q}_{\mathrm{BBH}}{{\boldsymbol{S}}}_{a})/(1+{q}_{\mathrm{BBH}})$, while ϕ₁ = 0–2π is the phase of the binary, extracted randomly between the two limiting values.

The amplitude of the kick depends critically on the spins of the two components and the binary mass ratio. Figure 11 shows the kick distribution for BH mergers with mass ratio between 0.01 and 0.9. As shown in the plot, the recoil is modest (<100 km s⁻¹) only for highly asymmetric mergers (i.e., with q < 0.1). As shown in Table 1, the mass ratio is of the order of q ∼ 0.07 for mergers involving an IMBH, and q ≃ 0.6 − 0.9 for stellar BH mergers, thus implying kicks of ∼10 km s⁻¹ and >100 km s⁻¹, respectively.

Zoom In Zoom Out Reset image size