The role of massive black holes in merging star clusters: dynamical evolution, stellar and compact object ejections, and gravitational waves

ABSTRACT

Star clusters can interact and merge in galactic discs, haloes, or centres. We present direct N-body simulations of binary mergers of star clusters with |$M_{\star } = 2.7 \times 10^4 \, \mathrm{M_{\odot }}$| each, using the N-body code bifrost with subsystem regularization and post-Newtonian dynamics. We include 500 |$\mathrm{M_{\odot }}$| massive black holes (MBHs) in the progenitors to investigate their impact on remnant evolution. The MBHs form hard binaries interacting with stars and stellar black holes (BHs). A few Myr after the cluster merger, this produces sizable populations of runaway stars (⁠|$\sim$|800 with |$v_{\mathrm{ej}} \gtrsim 50 \, \mathrm{km\, s^{-1}}$|⁠) and stellar BHs (⁠|$\sim$|30) escaping within 100 Myr. The remnants lose |$\sim 30{{\ \rm per\ cent}}$| of their BH population and |$\sim 3{{\ \rm per\ cent}}$| of their stars, with |$\sim$|30 stars accelerated to high velocities |$\gtrsim 300 \, \mathrm{km\, s^{-1}}$|⁠. Comparison simulations of isolated clusters with central hard MBH binaries and cluster mergers without MBHs show that the process is driven by MBH binaries, while those with a single 1000 |$\mathrm{M_{\odot }}$| MBH in isolated or merging clusters produce fewer runaway stars at lower velocities. Low-eccentricity merger orbits yield rotating remnants (⁠|$v_{\mathrm{rot}} \sim 3 \, \mathrm{km\, s^{-1}}$|⁠), but probing the presence of MBHs via kinematics alone remains challenging. We expect the binary MBHs to merge within a Hubble time, producing observable gravitational-wave (GW) events detectable by future GW detectors such as the Einstein Telescope and Laser Interferometer Space Antenna. The results suggest that interactions with low-mass MBH binaries formed in merging star clusters are an important additional channel for producing runaway and high-velocity stars, free-floating stellar BHs, and compact objects.

1 INTRODUCTION

Numerous escaping stars and compact objects (COs) with velocities large enough to escape the Galaxy have been detected or proposed in recent years. At the same time the dynamical evolution of star clusters in the presence of massive black holes (MBHs) has been widely studied, while evidence about interacting or merging star clusters is growing. In this section, we first outline the observational evidence of escaping stars and COs (Section 1.1), and then present evidence for the presence and formation pathways of MBHs in star clusters (Section 1.2) and the interactions and/or mergers of star clusters (Section 1.3).

1.1 High-velocity escaping stars and COs

The leading scenario of stellar ejections from the Galactic Centre (GC) with velocities exceeding the escape velocity of the Galaxy is based on the ‘Hills mechanism’ (Hills 1988), i.e. when a stellar binary comes sufficiently close to the supermassive black hole (SMBH) Sgr A* (GRAVITY Collaboration 2020; Event Horizon Telescope Collaboration 2022) in the GC, it can be tidally disrupted. This refers to the dynamical breakup of the binary in analogy to the disruption of a single star by an MBH. This results in the capture of one of the stars on a wide eccentric orbit and the ejection of the companion star with extreme velocities of |$1000 \, \mathrm{km} \, \mathrm{s^{-1}}$|⁠, by far larger than the escape velocity of the Milky Way (MW). The ejection due to close encounter with the MBH and velocities higher than escape velocity are the two properties that define a hyper-velocity star (HVS). The first HVSs were discovered by Brown et al. (2005) and Hirsch et al. (2005) with heliocentric radial velocities of |$709$| and |$708 \, \mathrm{km} \, \mathrm{s^{-1}}$| and at 50 and 19 kpc heliocentric distances, respectively, both consistent with GC origin. The very high velocities of those stars allow them to travel long distances following radial trajectories, which makes them ideal probes of the mass distribution of the Galaxy (Gnedin et al. 2005; Kenyon et al. 2008). In order to reveal their origin, those stars are backwards integrated in a Galactic MW-like potential (e.g. Price-Whelan 2017). The respective trajectory can then be used to probe the shape of the potential and dark matter halo profile (Gnedin et al. 2005; Yu & Madau 2007; Perets et al. 2009; Contigiani, Rossi & Marchetti 2019). An additional confirmation of the Hills mechanism was the discovery of an HVS with an extreme velocity of |$1800 \, \mathrm{km} \, \mathrm{s^{-1}}$| and a backwards trajectory pointing towards the GC and Sgr A* (Koposov et al. 2019). The discovery of the first HVS has lead to numerous dedicated searches (Li et al. 2011; Zheng et al. 2014; Huang et al. 2017; Du et al. 2019; Kreuzer, Irrgang & Heber 2020), reporting a large number of HVSs candidates in the Galactic halo (review in Brown (2015) and references therein) but the true origin for many of those still remains a mystery. An additional channel for the production of HVSs was proposed by Begelman, Blandford & Rees (1980) and Yu & Tremaine (2003), extending the classical Hills mechanism: the interaction of single or binary stars with a binary SMBH. Alternatively, the inspiral of an intermediate-mass (in the range of |$100 \, \mathrm{M_{\odot }} \lesssim M_{\bullet } < 10^6 \, \mathrm{M_{\odot }}$|⁠) black hole (IMBH) towards the GC (Baumgardt, Gualandris & Portegies Zwart 2006; Sesana, Haardt & Madau 2007; Sesana, Madau & Haardt 2009; Wang et al. 2018; Rasskazov et al. 2019; Evans et al. 2023) could also explain the origin and ejection mechanism of HVSs.

Another class of high-velocity stars is observed to have large peculiar velocities, in the range of |$40 \, \mathrm{km} \, \mathrm{s^{-1}} \le v_{\mathrm{pec}} \le 200 \, \mathrm{km} \, \mathrm{s^{-1}}$| (Blaauw 1961a; Gies & Bolton 1986; Hoogerwerf, de Bruijne & de Zeeuw 2001; Perets & Šubr 2012), commonly called ‘runaway stars’ (RASs). Those are usually early-type (O-type, B-type) stars with a Galactic disc origin. Their ejection mechanisms are usually divided in two types (Blaauw 1993; Hattori et al. 2019; Carretero-Castrillo, Ribó & Paredes 2023): (a) the dynamical ejection mechanism (DEM), where strong three- and four-body encounters between stars and/or BHs in dense stellar environments lead to the ejection of one of the members (Poveda, Ruiz & Allen 1967; Aarseth & Hills 1972; Hut & Bahcall 1983; Ryu, Leigh & Perna 2017a, b; Weatherford et al. 2023) and (b) the binary ejection mechanism, where the primary of the binary members undergoes a supernova (SN) explosion and the secondary (the RAS) moves on with close to its original orbital velocity. The primary remnant receives a kick from the SN, which might unbind it from the secondary (Zwicky 1957; Blaauw 1961b; Portegies Zwart 2000; Justham et al. 2008; Pakmor et al. 2013; Neunteufel 2020; Rajamuthukumar et al. 2024). Heber et al. (2008) have discovered an unbound star travelling with heliocentric radial velocity of |$>400 \, \mathrm{km} \, \mathrm{s^{-1}}$| with a Galactic disc origin and have classified it as a hyper-runwaway star. Finally, alternative mechanisms for the origin of high-velocity stars include the encounters with satellite galaxies and nearby galaxies such as the Large Magellanic Cloud (LMC) (Gualandris & Portegies Zwart 2007; Boubert et al. 2017; Erkal et al. 2018) and the Andromeda galaxy (Sherwin, Loeb & O’Leary 2008) and from star clusters tidally interacting with a single or a binary SMBH (Capuzzo-Dolcetta & Fragione 2015; Fragione & Capuzzo-Dolcetta 2016; Fragione, Capuzzo-Dolcetta & Kroupa 2017).

Our knowledge on stellar populations in the MW has been vastly expanded by the |$Gaia$| satellite mission of the European Space Agency and the unprecedented quality of astrometric and photometric data it provides. Already from the first data release (DR1) (Gaia Collaboration 2016), the number of HVS candidates has significantly increased (Marchetti et al. 2018a), contributing at least 14 more objects with a total velocity in the Galactic rest frame of |$>400 \, \mathrm{km} \, \mathrm{s^{-1}}$| and one classified as an HVS. A crucial step towards the improvement of modern astrometry came with the second data release (DR2) (Gaia Collaboration 2018) providing us positions and proper motions of more than one billion stars and the radial velocities for about 7 million of those. The number of HVS candidates reported in the literature before Gaia’s DR2 was close to 500 (Hirsch et al. 2005; Brown et al. 2006; Brown, Geller & Kenyon 2008, 2012, 2014; Li et al. 2011, 2015; Huang et al. 2017; Vennes et al. 2017), with about 20 of them being faint and blue stars located in the Galactic halo with very high radial velocities and being classified as HVSs. Utilizing Gaia DR2 Hattori et al. (2018) found 30 |$(>480 \, \mathrm{km} \, \mathrm{s^{-1}})$| old metal-poor stars, with two to three having LMC and one to two GC origin. Three white dwarfs (WDs) with extreme velocities of |$1000 \,\mathrm{km} \, \mathrm{s^{-1}} \le v \le 3000 \, \mathrm{km} \, \mathrm{s^{-1}}$| were also discovered in DR2 from Shen et al. (2018). Marchetti, Rossi & Brown (2018b) found 20 stars unbound (⁠|$> 80 {{\ \rm per\ cent}}$| probability) from the Galaxy, seven of which are consistent with a GC origin and 13 with an origin outside of the MW. By combining Gaia DR2 and the seventh data release of the Large Sky Multi-object Fiber Spectroscopic Telescope (LAMOST), Li et al. (2021) reported a total of 591 high-velocity candidates with 43 of them having more than |$50 {{\ \rm per\ cent}}$| probability of being unbound from the Galaxy, increasing the number of known candidates by a factor of 2. Despite the additional proper motions and radial velocities of 34 million stars that the third Gaia data release (DR3) (Gaia Collaboration 2023) has provided, Marchetti, Evans & Rossi (2022) found no additional HVS candidates |$>400 \, \mathrm{km} \, \mathrm{s^{-1}}$| and Liao et al. (2023) reported only two with radial velocities larger than |$500 \, \mathrm{km} \, \mathrm{s^{-1}}$| without strong evidence for a GC origin. Finally, the recent discovery (Huang et al. 2024) of a high-velocity star (J0731+3717) with |$v_{\mathrm{ej}} \approx 548 \, \mathrm{km} \, \mathrm{s^{-1}}$| and backward trajectory |$21 \, \mathrm{Myr}$| ago reveals its origin from MW globular cluster M15.

1.2 Evidence for MBHs in star clusters

The existence and formation pathway of IMBHs remains an open question. Dense stellar environments like dwarf galaxies and massive star and globular clusters are ideal sites for their formation and growth (Askar, Baldassare & Mezcua 2024). Although limited, observational evidence of IMBHs have been growing recently (see e.g. Mezcua 2017; Greene, Strader & Ho 2020, and references therein). One of the most recent and clear pieces of evidence to date is the detection GW190521 (Abbott et al. 2021) of an IMBH with |$M_{\bullet } \sim 150 \, \mathrm{M_{\odot }}$| via gravitational waves (GWs) from two coalescing stellar-mass BHs of |$m_{\bullet } \sim 66 \, \mathrm{M_{\odot }}$| and |$m_{\bullet } \sim 85 \, \mathrm{M_{\odot }}$| each. Apart from GWs numerous candidates of accreting IMBHs have been proposed due to X-ray emission detected in galaxies nearby, in the mass range of |$200\, \mathrm{M_{\odot }} \lesssim M_{\bullet } \lesssim 10^5 \, \mathrm{M_{\odot }}$| (Matsumoto et al. 2001; Strohmayer & Mushotzky 2003; Farrell et al. 2009; Mezcua et al. 2013, 2015; Mezcua 2017). Moreover, an IMBH of |$M_{\bullet } \sim 10^4 \, \mathrm{M_{\odot }}$| was recently proposed (Wen et al. 2021) as an X-ray source of a tidal disruption event (TDE) in a massive object of |$M_\mathrm{\star } \sim 10^7 \mathrm{M_{\odot }}$| (Lin et al. 2018), which corresponds either to a stripped galactic nucleus or a globular cluster.

Evidence for the existence of IMBHs or central populations of stellar BHs in various stellar and globular clusters is supported by dynamical signatures from observations and dynamical modelling. However, making a clear distinction between the two alternatives remains challenging. Among the top candidates of Galactic (MW) globular clusters hosting an IMBH are |$\omega$|Centauri (⁠|$\omega$|Cen) and 47 Tucanae (47 Tuc). For the former dynamical models (Noyola et al. 2010; van der Marel & Anderson 2010; Jalali et al. 2012; Baumgardt 2016) have proposed a central IMBH in the range |$1.2 \times 10^4 \, \mathrm{M_{\odot }} \lesssim M_{\bullet } \lesssim 5 \times 10^4 \, \mathrm{M_{\odot }}$|⁠. Only recently observations of proper-motions of fast-moving stars in the inner |$0.08 \, \mathrm{pc}$| of |$\omega$|Cen from Hubble Space Telescope (HST) have been analysed by Häberle et al. (2024), who found that the excess velocities can be explained by the presence of an IMBH with a corresponding lower mass limit of |$M_{\bullet } \sim 8.8 \times 10^3 \, \mathrm{M_{\odot }}$|⁠. On the other hand, the situation in 47 Tuc is still under debate. X-ray (Grindlay et al. 2001) and radio (de Rijcke, Buyle & Dejonghe 2006) data observations have provided upper mass limits of |$\sim 470 \, \mathrm{M_{\odot }}$| and |$\sim 2060 \, \mathrm{M_{\odot }}$| for a central IMBH in 47 Tuc, while recent dynamical modelling (Kızıltan, Baumgardt & Loeb 2017a, b) supports this evidence predicting a mass of |$M_{\bullet } \sim 2300 \, \mathrm{M_{\odot }}^{+1500}_{-800}$|⁠. Further observations of millisecond pulsars have shown that the presence of an IMBH is not required to explain the data (Freire et al. 2017; Abbate et al. 2018), in line with multimass dynamical modelling results (Mann et al. 2019; Hénault-Brunet et al. 2020). Updated dynamical models by Croce et al. (2024) place an upper limit of |$\sim 578 \, \mathrm{M_{\odot }}$| and only the latest ultradeep Australia Telescope Compact Array have revealed a central compact radio source associated with a faint X-ray emission corresponding to an IMBH of |$M_{\bullet } \sim 54\!-\!6000 \, \mathrm{M_{\odot }}$|⁠.

Numerous studies have been conducted providing additional evidence of a central IMBH or a dark central cluster (Ibata et al. 2009; Kamann et al. 2014; Baumgardt 2016; Nguyen et al. 2017; Gieles et al. 2018; Vitral et al. 2022) in such environments. A detailed analysis of 3D spectroscopic data by Kamann et al. (2016) revealed the presence of a |$M_{\bullet } \sim 600 \, \mathrm{M_{\odot }}$| IMBH in the core-collapsed cluster NGC6397, while Rui et al. (2021) suggested a central sub-system of stellar remnants supported by Vitral et al. (2022) who utilized Gaia and HST proper motion data, consistent with a dark concentration mass of |$\sim 1000 \, \mathrm{M_{\odot }}$|⁠. A less debated case of whether the central mass corresponds to a single IMBH or to a sub-system of dark objects is that of NGC6388. Combining spectroscopic data from the Very Large Telescope and HST, Lützgendorf et al. (2011) proposed an IMBH of |$M_{\bullet } \sim (1.7 \pm 0.9) \times 10^4 \, \mathrm{M_{\odot }}$|⁠, confirmed by Lanzoni et al. (2013) who placed an upper limit of |$\sim 2000 \, \mathrm{M_{\odot }}$| for the mass of the IMBH. Finally, similar findings have been proposed for Andromeda (M31) galaxy. For example, with the use of HST data, Gebhardt, Rich & Ho (2002) reported the detection of an |$M_{\bullet } \sim 2 \times 10^4 \, \mathrm{M_{\odot }}$| IMBH in G1 cluster of M31 galaxy. Dynamical models by Baumgardt & Makino (2003) proposed that the presence of an IMBH in G1 is not needed to explain the observed data, while X-ray emission detected by XMM–Newton could not distinguish between an accreting IMBH or a low-mass X-ray binary (Pooley & Rappaport 2006). Additionally, there has been recent evidence for an |$M_{\bullet } \sim 10^5 \, \mathrm{M_{\odot }}$| IxsMBH in G078 (Pechetti et al. 2022), which is the most massive globular cluster of M31, originating from the tidally stripped nucleus of a dwarf galaxy (Fuentes-Carrera et al. 2008).

One of the leading mechanisms for the formation of an IMBH is the growth of low-mass BHs (⁠|$m_{\bullet } < 50 \, \mathrm{M_{\odot }}$|⁠) in dense stellar environments via successive stellar and CO collisions (Stone, Küpper & Ostriker 2017; Rizzuto et al. 2020, 2022; Arca-Sedda et al. 2021, 2024a). The high densities in the centre of such environments leads to a sequence of mergers capable of producing BHs on the intermediate mass range |$100 \, \mathrm{M_{\odot }} \lesssim M_{\bullet } < 10^4 \, \mathrm{M_{\odot }}$| (Arca-Sedda et al. 2021; Mapelli et al. 2021; Fragione et al. 2022; Rizzuto et al. 2022; Atallah et al. 2023). Such merger products can receive velocity kicks ranging from tens up to |$5000 \, \mathrm{km} \, \mathrm{s^{-1}}$| (Campanelli et al. 2007; Lousto & Healy 2019) due to the GW-induced recoil. Consequently, it is rather unlikely that the final product can be retained in low escape velocity environments such as stellar and globular clusters (Gerosa & Berti 2019). The extreme densities of nuclear star clusters (NSCs) combined with the increased escape velocities |$v_{\mathrm{esc}} \gtrsim 100 \, \mathrm{km} \, \mathrm{s^{-1}}$| (Neumayer, Seth & Böker 2020) makes them ideal sites for the retention of IMBHs. The hierarchical assembly of massive star clusters and NSCs makes retaining the merger products more likely as the BHs end up residing in environments with higher escape velocities than their original birth clusters (Rantala, Naab & Lahén 2024).

The theory of pair-instability and pulsation pair-instability supernovae predicts a mass gap on BH formation in the range |$\sim 50\!-\!130 \, \mathrm{M_{\odot }}$|⁠, due to single stellar evolution only (Fowler & Hoyle 1964; Woosley, Blinnikov & Heger 2007; Woosley 2017). In young and massive star clusters, a star can grow to very high masses [very massive star (VMS) with |$m_{\star } > 150 \, \mathrm{M_{\odot }}$|] via repeated stellar collisions, which collapses to an IMBH within or above the mass gap, and may further grow through accretion during TDEs. This so-called fast runaway collision scenario has been numerically studied (Portegies Zwart et al. 1999; Portegies Zwart & McMillan 2002; Freitag, Gurkan & Rasio 2006b; Freitag, Rasio & Baumgardt 2006a; Rizzuto et al. 2020, 2023; Arca-Sedda et al. 2024a, b; Prieto et al. 2024; Rantala et al. 2024) where VMSs are efficiently formed, resulting in MBHs of |$M_{\bullet } \gtrsim 10^2\!-\!10^4 \, \mathrm{M_{\odot }}$| when they collapse. Finally, it is important to mention that the mass-loss from VMSs is still unclear and further research in stellar evolution models is needed to understand if and how long such stars could survive even if they do form (Sabhahit et al. 2022, 2023).

The formation of IMBHs and the co-evolution of massive star cluster hosts has been studied mostly in isolated numerical set-ups (Wang et al. 2016; Rizzuto et al. 2023) and only recently in full hierarchical by Rantala et al. (2024). Recent observations by JWST suggest the formation of clumped clusters like the Cosmic Grapes at |$z>6$| redshift (Fujimoto et al. 2024). Furthermore, five massive (⁠|$M_{cl} \sim 10^6 \, \mathrm{M_{\odot }}$|⁠) star clusters, the so-called C

osmic Gems, have been detected by JWST at |$z=10.2$|⁠, providing additional evidence to such hierarchical assembly. This is also supported by previous observations of star cluster forming regions (e.g. Zhang, Fall & Whitmore 2001; Bastian et al. 2005; Grasha et al. 2017; Menon et al. 2021) and simulations of star-burst environments (Lahén et al. 2020), where a monolithic collapse seems unlikely.

The conditions required to produce massive clusters, such as globular clusters, occur in interacting galaxies and star-forming regions during the collapse of giant molecular clouds (see e.g. Krumholz, McKee & Bland-Hawthorn 2019, and references therein). The combined effects of the properties of the collapsing cloud and the tidal field from the forming galaxy may lead to efficient angular momentum transport in the star cluster forming region, resulting in rotating clusters (Lahén et al. 2020). Although debated in the past, there is clear observational evidence nowadays of rotating GCs in the MW (Bellazzini et al. 2012; Fabricius et al. 2014) with rotational velocities of the order of a few |$\mathrm{km \, s^{-1}}$|⁠. Rotation can have significant effects on the overall evolution (Fiestas & Spurzem 2012) of a star cluster, for example, spinning up its core and increasing the stellar ejection rate (Einsel & Spurzem 1997; Ernst et al. 2007) or when combined with stellar evolution of the system (Kamlah et al. 2022).

Similar to the coupled formation and evolution of IMBHs and massive star clusters, the properties of NSCs and SMBHs are closely related to those of their host galaxies (Ferrarese et al. 2006; Leigh et al. 2015; Capuzzo-Dolcetta & e Melo 2017), which indicates that the formation and evolution of the two is tightly related (Antonini, Barausse & Silk 2015; Neumayer et al. 2020). Two scenarios have been proposed for the formation of NSCs. Galactic nuclei with sufficiently high gas densities can trigger in situ star formation (Loose, Kruegel & Tutukov 1982; Mihos & Hernquist 1994; Milosavljević 2004; Nayakshin, Cuadra & Springel 2007; Aharon & Perets 2015), leading to the formation of a NSC. Alternatively, star clusters spiral towards the GC and potentially merge, due to dynamical friction (DF; Tremaine, Ostriker & Spitzer 1975; Capuzzo-Dolcetta 1993; Loose et al. 1982; Agarwal & Milosavljević 2011; Tsatsi et al. 2016). The latter suggests that clusters may interact with each other during their infall towards the GC. In conclusion though, the diversity of stellar age and metallicity observed in galactic nuclei implies that both scenarios contribute to the overall evolution of NSCs (Antonini et al. 2015; Guillard, Emsellem & Renaud 2016; Arca-Sedda et al. 2020; Do et al. 2020; Neumayer et al. 2020).

1.3 Binary star clusters: interactions and mergers

The interaction and merging of star clusters is theoretically expected and observed during their formation (see e.g. de La Fuente Marcos & de La Fuente Marcos 2009; Lahén et al. 2020). Especially for LMC, it has been found that the fraction of binary star clusters is roughly |$\sim 12{{\ \rm per\ cent}}$| (Dieball, Müller & Grebel 2002), while a similar number |$(\sim 10{{\ \rm per\ cent}})$| is expected for the MW (de la Fuente Marcos & de la Fuente Marcos 2010). Bhatia et al. (1991) published a catalogue with numerous binary star cluster candidates in LMC, but their formation path for most of them remains unclear. Apart from forming together at birth, star clusters undergo tidal interactions and mass exchange at later evolutionary phases in the discs of galaxies (Khoperskov et al. 2018; Mastrobuono-Battisti et al. 2019; Camargo 2021; Ishchenko et al. 2024). Recent high-resolution observations of LMC provided evidence of pairs of star clusters in collision course due to tidal capture (TC; Mora, Puzia & Chanamé 2019; Giusti et al. 2023). The kinematic properties of almost the entire population of the MW globular clusters is nowadays available due to the Gaia DR3 release (Gaia Collaboration 2018, 2023). A detailed investigation on 16 Galactic binary cluster candidates from Gaia DR3 (Gaia Collaboration 2023) was recently done by Angelo et al. (2022) who found four pairs of bound open clusters. Ishchenko et al. (2023) conducted a statistical analysis for the interaction probability of the Galactic globular clusters utilizing the Gaia DR2 release (Gaia Collaboration 2018) and found a significant rate of close encounters and an average of 10 intersecting with-each-other trajectories per cluster. Additionally, Chemerynska et al. (2022) integrated the orbits of the 150 clusters in a static MW-like potential up to |$5 \, \mathrm{Gyr}$| and found a probability above |$20{{\ \rm per\ cent}}$| for cluster collisions. Finally, recent studies have shown that many globular clusters inhibit disc-like kinematics (Casetti-Dinescu et al. 2010; VandenBerg et al. 2013), providing evidence of close orbital passages and angular momentum gain leading to rotation. This leads to the formation of IMBH binaries potentially leading to GW-driven coalescence.

In this paper we investigate the presence of single and binary MBHs in star clusters as the potential origin for high-velocity stars. We use direct N-body simulations of idealized isolated and merging star clusters with and without central MBHs. We find that, in particular in the presence of MBH binaries, stars and COs can be accelerated to velocities up to 800 |$\mathrm{km} \, \mathrm{s}^{-1}$|⁠. The paper is structured as follows: in Section 2 we describe the various mechanisms leading to the production of high-velocity stars and COs. The simulation methods and initial conditions (ICs) are presented in Section 3, and the results about cluster evolution, the formation and evolution of MBH binaries, and the demographics of ejections are presented in Sections 4, 5, and 6, respectively. Finally, we summarize our conclusions in Section 7.

2 DEM IN STAR AND GLOBULAR CLUSTERS

In this section we review the basic DEM processes operating in star clusters, where a single or binary MBH and/or a subcluster of stellar-mass BHs are present (Kulkarni, Hut & McMillan 1993; Quinlan 1996; O’Leary et al. 2006; Sesana, Haardt & Madau 2006; Trenti et al. 2007; Morscher et al. 2015; Fragione & Capuzzo-Dolcetta 2016; Fragione & Gualandris 2019; Šubr, Fragione & Dabringhausen 2019; Rasskazov et al. 2019; Weatherford et al. 2023). If runaway and HVSs are observed and a star cluster origin is confirmed, it would be smoking gun evidence for the existence of MBH seeds in the heart of such systems. Additionally, understanding their dynamical production channels can be used to disentangle the tension between the presence of a single massive or a collection of stellar-mass ‘dark’ objects (BHs) in the cores of massive clusters (Anderson & van der Marel 2010; Noyola et al. 2010; Lu & Kong 2011; Strader et al. 2012; Chomiuk et al. 2013; Feldmeier et al. 2013; Peuten et al. 2016; Zocchi, Gieles & Hénault-Brunet 2017, 2018; Bellini et al. 2018; Vitral et al. 2022; Häberle et al. 2024; Paduano et al. 2024; Huang et al. 2024).

2.1 Single stellar encounters and two-body relaxation

The gravitational interaction between two stars in the vicinity of a BH, with masses |$m_1$| and |$m_2$|⁠, results in a change in the velocity of the |$m_1$| (Yu & Tremaine 2003),

where |$v_{1 2}$| is the relative velocity of the two stars and b is the impact parameter (Binney & Tremaine 2008). Although rare, an encounter of this type with small enough impact parameter can result in significant velocity kick. For example, if |$m_1=m_2=1 \mathrm{M_{\odot }}$| and |$b=10^{-6} \mathrm{pc}$|⁠, then |$\delta v_1 \ge 30 \mathrm{~km} \mathrm{~s}^{-1}$|⁠.

The collective motion of objects in an N-body system induces time-varying perturbations on the background potential of the system. A non-smooth potential introduces energy and angular momentum exchange between the cluster members, causing them to randomly diffuse in the phase space. The net effect of all the uncorrelated two-body encounters can lead to a change of each body’s velocity by order of itself |$\Delta v / v \sim 1$| (Chandrasekhar 1942; Spitzer 1987; Heggie & Hut 2003; Aarseth, Tout & Mardling 2008), allowing an object to escape within the relaxation time-scale and velocity of the order |$v_{2b} \approx 2\sqrt{\sigma }$|⁠, where |$\sigma$| is the velocity dispersion the cluster. A star located in the outskirts of a cluster with local velocity dispersion |$\sigma \approx 2 \, \mathrm{km} \, \mathrm{s}^{-1}$| could then escape through that process with |$v_{2b} \approx 2.8 \, \mathrm{km} \, \mathrm{s}^{-1}$|⁠.

2.2 Encounters of singles with binary stars

Interactions between single and binary stars are common in star clusters and play a key role to binary evolution (Heggie 1975). Those interactions can lead to the formation of new binaries from single ones or through the exchange of binary members (Valtonen & Karttunen 2006). Lower mass objects are accelerated to high velocities when encounters between unequal mass stars take place (Spitzer 1987), potentially leading to their escape from the host cluster. The critical velocity of a star (member of the cluster) capable of ionizing (the process where an encounter leads the breakup of a bound system to its individual components) the binary in the centre of mass of a triple system is given by (Hills 1975a; Hut & Bahcall 1983; Sigurdsson & Phinney 1993)

Here |$\mu _{12} = m_1 m_2 / (m_1 + m_2)$|⁠, |$m_{123} = m_1 + m_2 + m_3$|⁠, and a is the semimajor axis of the binary with |$m_1 \ge m_2$|⁠, where the least massive body is the most likely to escape. For a binary with |$m_1= 2 \, \mathrm{M}_{\odot }$|⁠, |$m_2= 1 \, \mathrm{M}_{\odot }$|⁠, and |$a= 2 \, \mathrm{AU}$|⁠, a third star of |$m_3 = 1 \, \mathrm{M}_{\odot }$| would need |$v_c = 33.8 \, \mathrm{km} \, \mathrm{s}^{-1}$| to ionize the binary. The stochastic nature of such process makes it very challenging to predict the ejection velocity of the escaping body from generic ICs, but numerous detailed numerical three-body scattering experiments have been conducted (e.g. Hills 1975a; Hut & Bahcall 1983; Sigurdsson & Phinney 1993) for the equal-mass case and recently from Forastier et al. (2024) for unequal masses. Finally, Samsing, MacLeod & Ramirez-Ruiz (2017) have investigated the role of stellar tides and post-Newtonian (PN) dynamics on such encounters and found that including those effects lead to higher stellar coalescence rates.

2.3 Binary encounters with a single MBH

The strong tidal forces a star experiences once exerted by an SMBH are capable of tearing it apart once it gets sufficiently close. This TDE takes place when the orbital pericentre distance |$r_{\mathrm{peri}} \le r_t$|⁠, where |$r_t$| is the tidal radius given by (Hills 1975b)

Consider now a binary system with mass |$m_{\mathrm{b}}$| and semimajor axis |$a_{\mathrm{b}}$| approaching an SMBH. The same process (replacing |$r_{*}$| with |$a_{\mathrm{b}}$| and |$m=m_{\mathrm{b}}$|⁠) will lead to the tidal breakup of the binary, where one member is captured by the SMBH and the other ejected with |$v_{\mathrm{ej}}$| (Hills 1988). This happens when it passes closer to the SMBH than the tidal radius, i.e. when

For a star cluster with an IMBH at its centre and a number of binaries inside its core, the above mechanism may efficiently produce high-velocity stars escaping the cluster. During the tidal breakup, the stars receive a velocity change of the order of the velocity relative to the centre of mass of the binary, which for |$m_1$| in this case is (Hills 1988; Yu & Tremaine 2003)

or in a more practical form by Bromley et al. (2006), the velocity of the ejected star is

Using equation (6), for a binary of |$m_{\mathrm{b}}=m=1\mathrm{M_{\odot }}$| approaching an IMBH of |$M_{\mathrm{IMBH}}=1000 \, \mathrm{M_{\odot }}$|⁠, the ejection velocity can vary from |$v_{\mathrm{ej}} \approx 136 \mathrm{~km} \mathrm{~s}^{-1}$| for |$a_{\mathrm{b}} = 0.1 \mathrm{mpc}$| up to |$v_{\mathrm{ej}} \approx 1360 \mathrm{~km} \mathrm{~s}^{-1}$| for |$a_{\mathrm{b}} = 0.01 \mathrm{mpc}$|⁠. For star clusters with large binary fractions, this channel acts as the major contribution to the production of RAs and HVSs (Šubr et al. 2019; Fragione & Gualandris 2019; Weatherford et al. 2023). Finally, Ryu et al. (2023a) showed by means of hydrodynamics simulations that in the special case where a tight stellar binary interacts with a stellar-mass BH, both stars can be tidally disrupted if the impact parameter is small enough or end up with one of the members becoming unbound (micro-Hills mechanism). If the binary consists of a star and a BH, then the stellar ejection could also be accompanied by the formation of a binary BH via binary member exchange (Ryu et al. 2023b).

2.4 Encounters with (massive) BH binaries

Star clusters are expected to host a population of binary BHs and COs (Kulkarni et al. 1993; Portegies Zwart & McMillan 2000; Downing et al. 2010; Tanikawa 2013; Morscher et al. 2015; Rodriguez et al. 2015; Park et al. 2017; Anagnostou, Trenti & Melatos 2020; Torniamenti et al. 2022). The presence of bodies with a range of masses in a star cluster leads to the exchange of kinetic energy and ultimately more massive objects tend to shrink to the centre of the cluster. During this stage, the density of the star cluster core is increased leading to the efficient dynamical formation of stellar and stellar remnant binaries. Encounters between single stars and COs with BH binaries can lead to sufficiently high-velocity kicks and result in ejections, especially when their semimajor axis |$a_{\mathrm{BHB}}$| becomes small enough, such that |$a_{\mathrm{b}} \le a_{\mathrm{h}} \equiv G M_2 / 4 \sigma _{c}^2$| (Sesana et al. 2006), where |$\sigma _c$| is the local dispersion velocity (Quinlan 1996). Additionally, a single star approaching a black hole binary (BHB) can be tidally disrupted and further harden the binary by up to |$20 {{\ \rm per\ cent}}$| for retrograde encounters (Ryu, Perna & Wang 2022). The velocity kick of a single or multiple encounters of a star or CO with a BHB is then given by (Yu & Tremaine 2003),

During the hierarchical formation of a massive star cluster (Gieles & Portegies Zwart 2011; Rantala et al. 2024), lower mass seed BHs can also form binaries. Those BH binaries are usually more massive than stellar remnants, since they can grow by successive mergers or as by-products of runaway collisions (Arca-Sedda et al. 2024a, b; Fujii et al. 2024). Moreover, similar to how galaxies merge (e.g. Milosavljević & Merritt 2001) leading to the formation of SMBH binaries, Galactic globular clusters can also interact and merge. If such clusters contain a massive MBH, their merger remnant becomes the host of an MBH binary, which can then scatter single and binary stars around, resulting into their ejection from the cluster. Finally, an MBH-host star or globular cluster on an inspiral orbit towards the GC captured by the SMBH (e.g. Baumgardt et al. 2006; Sesana et al. 2007, 2009) may lead to the formation of an SMBH–IMBH binary, which can eject stars at high speed (Rasskazov et al. 2019) once it gets hard (⁠|$a_{\mathrm{b}}< a_{\mathrm{h}}$|⁠).

In the special but common in dense stellar environments like globular clusters and the GC case where all interacting bodies are BHs, single-binary encounters in star clusters can efficiently lead to the production of GW-driven merger events (Samsing & Ilan 2018; Samsing, Askar & Giersz 2018b; Samsing et al. 2018a), where one of the BHs receives a strong velocity kick enough to be ejected from the cluster.

3 METHODS

To investigate the influence of MBHs on the merging of stellar clusters and the resulting population of escaping stars and COs (i.e. WDs, neutron stars, and BHs), we conduct a series of direct N-body simulations. These simulations involve merging pairs of stellar clusters, where either both, only the primary, or none of them hosts an MBH at their centre. The systems are evolved for 100 Myr using the direct N-body code bifrost (Rantala et al. 2023), which we describe briefly in Section 3.1.

3.1 N-body code

Throughout this study we used the GPU-accelerated direct-summation N-body code bifrost (Rantala et al. 2023), which is the updated version of the frost code (Rantala, Naab & Springel 2021). The code employs a fourth-order forward symplectic integrator (FSI) technique (Chin 1997; Chin & Chen 2005; Dehnen & Hernandez 2017), using a hierarchical implementation (Rantala et al. 2021). Beyond traditional Newtonian accelerations, FSI incorporates gradient accelerations to cancel second-order error terms, resulting in a fourth-order time integration cycle with strictly positive sub-steps. This method was found to be more effective (Chin 2007) e.g. for the Kepler problem than traditional Yoshida-type symplectic integrators (Yoshida 1990), which use negative sub-steps for higher than second integration order. Furthermore, within the hierarchical integration approach, rapidly evolving components of the simulated systems decouple from slowly evolving components, a feature which makes the code particularly efficient when applied to N-body systems characterized with a large dynamical range (Pelupessy, Jänes & Portegies Zwart 2012).

The code handles sub-systems such as binary and multiple systems, fly-bys, and small clusters around SMBHs, through the use of secular and regularized few-body integrators (Rantala et al. 2020). The code also includes the option of enabling PN equations of motion up to order PN3.5 utilizing the formulas given in Thorne & Hartle (1985) and Blanchet, Buonanno & Faye (2006). In binary systems the PN corrections account for relativistic effects such as periastron precession and radiation-reaction (circularization and shrinking of the orbit) due to the emission of GWs. Another feature that few codes share (Wang et al. 2020) is the ability to efficiently simulate massive systems with arbitrary fraction of primordial binary systems.

bifrost makes use of four different criteria for stellar and CO mergers, which we describe in detail Section 3.3. Those involve merging events of COs, disruption of a star by a CO, or stellar collisions. For unbound stars at large distances from the centre of mass of the system, we employ a user-defined threshold at which the particles are removed from the simulation. For a typical star cluster, stars and COs found beyond |$r_\mathrm{esc} \sim 100\, r_\mathrm{h}$| are not considered cluster members anymore, since at such distances the gravitational pull from the parent cluster would be too weak against the Galactic tides. We note that the simulations presented in this study do not for simplicity include any tidal potential. The main user-defined code parameters are listed in Table 1. Overall, the parameters bifrost used are similar to used in previous bifrost studies (e.g. Rantala et al. 2023; Rizzuto et al. 2023).

3.2 Orbital evolution of binary systems

Binary systems that are strongly perturbed or PN are treated in bifrost through algorithmic regularization. We make use of secular binary integration when a binary is sufficiently isolated (no other stars or COs within 10 mpc) and the number of binary orbits per time-step |$\epsilon$|⁠, |$N_{\mathrm{orb}}=P_{\mathrm{bin}}/\epsilon$|⁠, is larger than the user-defined threshold |$N_{\mathrm{orb,sec}}$|⁠, where |$P_{\mathrm{bin}}$| is the orbital period of the binary. We use |$N_{\mathrm{orb,sec}}=0$| throughout the study. The ICs of the simulations presented here do not contain a binary population. We only form binary systems dynamically.

The equations of motion responsible for the secular evolution are characterized by the semimajor axis a, orbital eccentricity e, and the argument of periapsis |$\omega$|⁠, taking into account leading-order relativistic effects that originate from the post-Newtonian PN1.0, PN2.0, and PN2.5. The advance of the orbital periapsis is due to PN1.0 and PN2.0 given by

The PN2.5 term refers to GW radiation reaction due to energy and angular momentum losses, leading to the circularization and shrinkage of the orbit described by the rate of change of binary semimajor axis and eccentricity (Peters 1964) as

in which the auxiliary functions |$\beta (m_\mathrm{1},m_\mathrm{2})$|⁠, |$F(e)$|⁠, and |$G(e)$| are defined as

The secular equations of motion are integrated with the use of a second-order leapfrog integrator, as described in detail in Rantala et al. (2023).

3.3 Merger criteria

During a simulation run two objects can merge if any of the four merger criteria is satisfied. Stars can become gravitationally unbound due to strong encounters with the other cluster members. The type of merging pairs can either be COs (WDs, neutron stars, BHs), tidal disruptions (BH with star), and stellar mergers (stars with stars). The first criterion for a compact binary merger depends on the GW-driven inspiral time-scale |$\mathrm{t_{gw}}$|⁠. If |$\mathrm{t_{gw}}$| is shorter that the binary’s current time-step in the time-step hierarchy, then the particles are merged. COs also merge if the radius of the (relativistic) innermost stable circular orbit (⁠|$R_{\mathrm{ISCO}}$|⁠) is larger than their mutual separation. Finally, stars can be tidally disrupted by a CO (TDE) if the pericentre distance becomes smaller than the tidal disruption radius or in case of a direct collision between two stellar particles, i.e. when their radii overlap.

A common approach for detecting CO mergers (e.g. Rizzuto et al. 2020; Arca-Sedda et al. 2021, 2023; Rizzuto et al. 2022) is to compute the time-scale |$\tau _\mathrm{gw}$| related to the time |$\tau _\mathrm{gw}$| a compact binary needs to reach coalescence. This time-scale can be evaluated through the integral expression (Peters 1964)

where |$a_0$| and |$e_0$| are the initial semimajor axis and eccentricity of the binary, |$\beta (m_\mathrm{1},m_\mathrm{2})$| is the function defined in equation (10), and |$g(e)$| is defined (Maggiore 2007) by

The expression for |$\tau _\mathrm{gw}$| for |$e=0$|⁠, i.e. circular orbits, is then

We tabulate the values of this integral by evaluating it separately (before the actual runs) for a large number of |$e_0$|⁠, in order to avoid computing it every time-step. Then by interpolating for those table values, we have an estimate for |$\tau _\mathrm{gw}$| throughout the simulation. However, extra caution should be taken when using the |$\tau _\mathrm{gw}$| criterion due to the strong dependence of eccentricity in the computed integral value. Even tiny fluctuations of eccentricity (e.g. from weak three-body and fly-by encounters) would result in significant underestimation of |$\tau _\mathrm{gw}$| and false merger identification.

The second criterion for merging COs is based on the innermost stable circular orbit (ISCO) with radius

where c is the speed of light and |$R_\mathrm{sch}$| the Schwarzschild radius of the BH. Two COs are then merged if their mutual distance is lower than |$r_\mathrm{isco}$|⁠.

If one of the particles is a CO and the other one is a star, then the latter could be disrupted. This happens when (Rees 1988; Kochanek 1992)

where |$m_\star$| and |$r_\star$| is the mass and radius of the star and |$M_\bullet$| is the mass of the CO. Two stars are assumed to merge when they overlap i.e. when their distance is shorter than the sum of their radii

3.4 ICs and simulations

We generate ICs for pairs of two identical star clusters of |$N=64\,000$| single stars without primordial binaries and a total mass of |$M_{\star } = 2.7 \times 10^4 \mathrm{M}_{\odot }$|⁠, each using the mcluster code (Küpper et al. 2011). The single star cluster model follows a King density profile (King 1966) with a fixed half-mass radius of |$R_{h}=1 \, \mathrm{pc}$| and central potential parameter related to the compactness of the cluster |$W_{0}=5$|⁠, which represents a moderately compact cluster. The age of the clusters is |$t_{\mathrm{age}}= 1 \, \mathrm{Gyr}$| and the stellar masses are sampled from a Kroupa initial mass function (IMF) (Kroupa 2001) as zero-age main-sequence stars in the range of |$0.08\!-\!100 \mathrm{M}_{\odot }$|⁠. The age of the cluster has been selected such that stellar evolution processes do not contribute to the overall evolution of the system, allowing us to focus on the long-term dynamical evolution. The mass range of stellar particles in the simulations is |$0.08\, \mathrm{M_{\odot }} \le m_{\star } \lesssim 1.89 \, \mathrm{M_{\odot }}$|⁠, |$0.88\, \mathrm{M_{\odot }} \le m_{\mathrm{WD}} \lesssim 1.37 \, \mathrm{M_{\odot }}$| for WDs, |$1.1\, \mathrm{M_{\odot }} \le m_{\mathrm{NS}} \lesssim 1.91 \, \mathrm{M_{\odot }}$| for NSs, and |$5.6\, \mathrm{M_{\odot }} \le m_{\mathrm{BH}} \lesssim 40.5 \, \mathrm{M_{\odot }}$| for stellar-mass BHs. From the ICs, we have initially |$N_{\mathrm{BH}}=250$|⁠, |$N_{\mathrm{NS}}=944$|⁠, and |$N_{\mathrm{WDs}}=4100$|⁠. The number of retained NSs in the ICs is related to the way natal kicks for SN remnants are assigned in Küpper et al. (2011), i.e. when the recoil velocities are low so the remnant cannot escape the cluster. With the use and calibration of pulsar proper motions (Hobbs et al. 2005; Kapil et al. 2023) the recoil kicks reach high velocities resulting in small fraction (only a few) of retained NSs in the initial system. Examples of a more sophisticated prescription for the number of retained NSs in the ICs can be found in Banerjee et al. (2020). The NS population in our simulations does not significantly contribute to the overall ejection demographics, but could potentially lead to false rates of escaping NSs and the formation of a larger number of NS binaries.

We use three different initially bound orbits for the merging star clusters. Each orbit has the same fixed semimajor axis |$a_{\mathrm{semi}}=2 \, \mathrm{pc}$| (and thus a fixed orbital energy, since |$E = -G M_1 M_2/2a_{\mathrm{semi}}$|⁠) and three different values of eccentricity |$e=[0.05,0.5,0.99]$|⁠. The second cluster is located at the apocentre of the orbit with |$r_{\mathrm{apo}}=[1.9,1.0,0.02] \, \mathrm{pc}$| for the three orbits. For each of the merger initial set-ups, we explore three scenarios: neither of the clusters host an MBH, the primary cluster host an MBH of |$M_{\bullet }=M_1=1000 \, \mathrm{M}_{\odot }$|⁠, and when both clusters host an MBH of |$M_1=M_2=500 \, \mathrm{M}_{\odot }$|⁠. Furthermore, we run three additional simulations where we consider a star cluster with |$N=128000$|⁠, |$M_{\star } = 5.4 \times 10^4 \, \mathrm{M}_{\odot }$|⁠, and |$R_{h}=1 \, \mathrm{pc}$| without a prior merger. For these simulations we either place a single, a binary with |$a_{\mathrm{b}}= 0.01 \, \mathrm{pc}$| and |$e_{\mathrm{b}}= \, 0.5$|⁠, or no MBH at all. We follow the evolution of the merger remnants or the single clusters up to |$t=100 \, \mathrm{Myr}$|⁠. An overview of the merger progenitors and their orbital elements is presented on Table 2.

In the next sections, we present the main outcomes of our stellar cluster merger simulations. We discuss how the presence of a single or a binary MBH affects the evolution of the merger remnants, their kinematic and structure properties, and the demographics of the escaping stellar and CO populations, which is the main focus of this work. An example outlook (here for the |$e=0.05$| cluster orbits) of the simulations at different times (⁠|$t = 0.5,1,15 \, \mathrm{Myr}$|⁠) is shown in Fig. 1. In the figure we show the surface stellar density and highlight the velocity vector |$\vec{V}$| field of unbound stars and COs.

If an MBH is present in both of the merging clusters, a massive black hole binary (MBHB) forms. The binary orbit keeps shrinking (monotonic decrease of |$a_{\mathrm{b}}$|⁠) until the end of the simulation, due to interactions with the background stars and COs. We describe the evolution of the MBHBs in detail in Section 5.

4 PROPERTIES OF THE STAR CLUSTERS

Adopting the definition by Trenti, Vesperini & Pasquato (2009), the core size or core radius of a star cluster following a King density profile is given by Casertano & Hut (1985)

where |$\rho _i$| is the local density around a star i defined as

and |$r_n$| is the distance to the nth nearest neighbour i, and |$\mathbf {r}_{\mathrm{d}}$| is the position of the density centre:

When the cores of the two progenitor clusters overlap, i.e. when their relative distance |$d_{12}$| becomes shorter than |$r_c$|⁠, we assume that the clusters have merged. The evolution of the core radii and the respective enclosed number density is shown in Fig. 2. From Fig. 2 we see that the presence of a single or binary MBH in the remnants leads to more rapid expansion due to a higher mass-loss rate. Essentially, particles interacting with the MBHs are radially scattered outwards, leading to constant decrease of the number density inside the core (bottom panel in Fig. 2), which indicates mass-loss inside |$r_c$|⁠.

Since the core radius |$r_c$| is not exactly the same (in the range |$0.6\mathrm{pc}< r_c< 0.8\mathrm{pc}$|⁠) for the various progenitor clusters, we define a threshold distance such that when |$d_{\mathrm{12}}< d_{\mathrm{crit}}=2 \mathrm{pc}$|⁠, the two clusters have merged. The time-scale of this process depends on the orbital eccentricity of the merger and on the number of MBHs in the system. Specifically for remnants with two MBHs, this happens at |$t=1 \, \mathrm{Myr}$|⁠, while for the other cases ranges between |$t=1.25 \, \mathrm{Myr}$| and |$t=2 \, \mathrm{Myr}$|⁠.

The density evolution for |$t>2 \, \mathrm{Myr}$| until the end of the run is only moderately affected by the merger process and is mainly driven by the presence of the single or binary MBHs. Fig. 3 shows the final (at |$t = 100\mathrm{Myr}$|⁠) density profiles of the merger remnants and that of the stellar BHs only. Overall, we observe an increase of the central density, which inversely scales to the presence and number of MBHs in the remnant, with an even more clear signature for the sub-system of BHs (right panel in Fig. 3). The binary transfers its orbital energy,

to the surrounding stars either through DF (Begelman et al. 1980; Varisco et al. 2021) or through three-body encounters (Hills 1991) accompanied by stellar ejections that carry out energy and angular momentum, leading to a decrease of stellar density around the MBHB. This process, known as core scouring, is likely responsible for the presence and formation of flat central cores in luminous galaxies (e.g. Rantala et al. 2018), and can likely affect the central properties of merged star cluster remnants as well. The characteristic time-scale over which the energy is extracted from the binary is given by Merritt (2013),

where |$\sigma$| and |$\rho$| are the three-dimensional velocity dispersion and density of the stellar background and C is a constant that depends on the energy transfer mechanism.¹ In our simulation models we have |$T_{\mathrm{E}} = 1 - 3.7 \, \mathrm{Myr}$|⁠, which is consistent with the time needed |$(\sim 1 \, \mathrm{Myr})$| for the respective progenitor clusters to merge.

4.1 Star cluster evolution

A bound self-gravitating system admits negative heat capacity |$C_{\mathrm{V}}$|⁠, allowing the release of energy from the centre to the outer parts, leading to the formation of a contracting core with radius |$r_c$| (equation 17) and an expanding halo (Lynden-Bell, Wood & Royal 1968; Lynden-Bell 1999).

In the absence of an additional energy source in the centre to halt the contraction the core will eventually collapse, a process called the gravothermal catastrophe (Antonov 1960, 1961, 1962). The core undergoes a series of contraction–expansion phase, known as gravothermal oscillations (Sugimoto & Bettwieser 1983; Heggie 1993). For a single and equal-mass component system in isolation, the core collapse would happen at about |$15t_{\mathrm{rh}}$| (Cohn 1980), where |$t_{\mathrm{rh}}$| is the half-mass relaxation time (Spitzer 1987).

Here N is the star/particle number of the system, |$R_h$| and |$\bar{m}$| the half-mass radius and average stellar mass, and |$\Lambda =\gamma N$| is the Coulomb logarithm. In our simulations we sample the stellar masses from a Kroupa IMF,² which corresponds to |$\gamma =0.02$| (Giersz & Heggie 1997) for multimass clusters. In realistic (multimass) star clusters, higher mass stars have lower velocity dispersion than the low-mass ones, driving the system towards energy equipartition, i.e. evenly distributing the kinetic energy of the system (Spitzer 1969, 1987). The massive stars segregate to the centre by transferring their kinetic energy to low-mass stars, thus expanding the orbits of the latter. The segregation time (Spitzer & Hart 1971; Portegies Zwart et al. 2004) is defined as

where |$M_{\mathrm{cl}}$| is the total mass of the cluster and |$m_{\mathrm{max}}$| is the mass of the most massive object. Since we do not assume any degree of primordial mass segregation in our ICs, we compute the segregation time of the merger remnants at |$t= 2 \, \mathrm{Myr}$|⁠, when all the progenitors have merged. Using equation (23) for the remnant clusters without an MBH, we find |$t_{\mathrm{s}} \sim 8\!-\!10 \, \mathrm{Myr}$|⁠. In time-scales of the order of |$\sim t_{\mathrm{s}}$| COs, the most massive particles sink towards the centre of the remnant (Baumgardt & Makino 2003).

The evolution of the radii enclosing |$3\!-\!80{{\ \rm per\ cent}}$| of the cumulative total mass (Lagrangian radii) of all stellar and COs for all simulations is shown in the top panel of Fig. 4. The left, middle, and right columns show the systems with no, single, or binary MBHs. For the isolated cluster the centre of the system without MBH (left) undergoes core contraction on a time-scale of |$\sim$| 20–30 Myr while the half-mass radius expands. The systems with single and binary MBH (middle and right) show no core collapse as rapidly forming binaries with the central MBH prevent this (see e.g. Rizzuto et al. 2023, and references therein). Stars and COs bound to MBHs interact, heating up the core and lead to its mild but constant expansion. The evolutionary behaviour is very similar for cluster merger remnants (second to bottom row). For mergers without an MBH the one with the highest eccentricity (bottom left panel) shows most core contraction. For systems with one or two MBHs the eccentricities of the merger orbits do not change the structural evolution.

In the bottom panel of Fig. 4 we show the Lagrangian radii evolution of the respective BH subsystems only. Initially, the general behaviour is the same as for the entire cluster (left panels in Fig. 4). However, the BH subsystem contracts faster on a |$\sim$|10 Myr time-scale, which is the expected mass segregation time-scale (see equation 23) in all cases without MBHs (left column). However, the presence of just one MBH does not prevent the collapse in for the BH subsystem. In addition the BHs form a central, stable subcluster, which is not expanding like the entire cluster population. The dynamical interactions in this central BH cluster provide a continuous source of energy for ejections of stars and COs.

The sphere around an MBH of mass |$M_{\bullet }$| where its presence dominates the gravitational potential of the system is called sphere of influence (SOI). The radius of this sphere |$r_{\mathrm{SOI}}$| is then given by (Merritt 2013)

where |$M_{\star }$| refers to the stellar mass.

The interactions of stars within |$r_{\mathrm{SOI}}$| with the MBH/MBHB can significantly affect their growth through TDEs (Rizzuto et al. 2023) or GW-merger events, but also play an important role in the evolution of the remnant cores (e.g. Aros et al. 2020, and references therein). Fig. 5 shows the fraction of bound stars to the single or binary MBH, i.e. stars with negative orbital energy |$E=\frac{m_* v^2}{2}-\frac{G m_* M_{\mathrm{BH}}}{r}<0$|⁠, where |$m_{*}$| is the mass of the individual star and |$r,v$| are the relative distance and velocity w.r.t. the MBH or the centre of mass of the MBH binary. After the cluster mergers the bound fractions are |$f_\mathrm{bound}\sim 0.14$|–0.22 in our simulations. The MBH case behaves in a qualitatively similar manner as the binary MBH set-up: the bound fractions reach values of |$f_\mathrm{bound}\sim 0.1$| in both cases. This reflects the similar overall evolution of central cluster properties such as stellar densities and velocity dispersions. We however note that the bound fraction decreases somewhat more rapidly in the MBHB set-ups. For clusters with a single MBH we find similar behaviour to that of the |$M_{\bullet }=2000 \, \mathrm{M}_{\odot }$| case in Rizzuto et al. (2023). We note though that their clusters are more dense and that stellar-mass BHs were not present in their ICs. The remnants with an MBH binary maintain approximately the same fraction of bound stars until the end of the simulation. The bound stars in the vicinity of the MBH or the MBHB provide enough energy supply to support the expansion of the core (Fig. 2) and the remnant itself. Moreover, encounters with the MBHB lead to three- and four-body encounters (Section 2.4), capable of kicking stars and COs in the outer parts of the cluster. The lower panel of Fig. 5 depicts the total energy loss from the bound cloud of stars inside |$r_\mathrm{SOI}$| leading to its growth. Those encounters drain orbital energy from the MBH binaries, leading to the shrinking of their orbits (decrease of semimajor axis). This hardening process is further discussed in Section 5. The energy inside the cloud in the form of dynamical heating combined with that from the BH sub-system (Mackey et al. 2008) leads to the expansion of the remnant cores and monotonic decrease of the number density (bottom panel in Fig. 2). Subsequently, the drop of the number density is associated with a decrease in the mass of the bound cloud, i.e. the gravitational potential of the core gradually becomes more shallow (i.e. a drop in the local escape velocity |$V_{\mathrm{esc}}$|⁠), allowing more particles to escape to the outer parts of the cluster.

Finally, the merger process itself has an effect on the differential energy distribution of the cluster stars. For a bound gravitational system in virial equilibrium |$2T+V=0$| (kinetic energy T and potential energy V) with |$T=-E$| and |$V=2E$|⁠. The total energy is simply |$E=T+V$|⁠. During the merger, the system is not in equilibrium, which means that although the total energy E remains constant, T and V will oscillate around their values, leading to the widening of the differential energy distribution (Hilz et al. 2012). This results in the most bound particles becoming even more bound, while the weakly bound particles might gain enough energy to escape from the local gravitational potential.

This process, known as violent relaxation, further contributes to the expansion of the core and the cluster itself, comparing to clusters without any prior merger (dashed lines in Fig. 2).

Fig. 6 shows the particle energy distribution at |$t=10$| and |$t=100 \, \mathrm{Myr}$|⁠. We see that clusters with MBHs hold a large fraction of unbound particles at early times corresponding to potential ejections. It is clear that the presence of MBHs broadens the energy distribution, especially for stars with |$E>0$|⁠, we observe a persistent tail in the distribution even after |$t=100\,\,\mathrm{Myr}$|⁠. From the radial velocity |$V_{\mathrm{r}}$| distribution (Fig. 7), we see that single MBH clusters develop a positive |$V_{\mathrm{r}}$| tail that is larger for MBH binaries. This tailed distribution is closely connected to the number and velocity of ejected bodies, as we describe in Section 2.

4.2 Star cluster shapes

In order to determine the inherent three-dimensional structure of the star clusters, we compute the reduced inertia tensor given by Gerhard (1983) and Bailin & Steinmetz (2005)

where k is the total particle number. The eigenvectors and eigenvalues of this tensor correspond to the directions of the principal axes |$a,b, \mathrm{ and }\,\, c$| (major, intermediate, and minor axis, respectively, i.e. |$a>b>c$|⁠) of the remnant. Their ratios can then be used to define the triaxiality parameter T (Jesseit, Naab & Burkert 2005; Binney & Tremaine 2008) given by

Fig. 8 shows the evolution of the principal axes ratios and the triaxiality parameter within the half-mass radius |$r_{\mathrm{h}}$|⁠, where |$T=1$| correspond to prolate, |$T=0.5$| to triaxial, and |$T=0$| to oblate shapes.

The shape of merger remnants is an indicator of the orbital families that are generated during the merger process (for a detailed analysis on a galaxy merger context, see for example Frigo et al. 2021). In Fig. 8 we present the time evolution of the remnant shapes. The size of the markers indicates the time-evolution of cluster shapes, where the largest ones correspond to |$t=100 \, \mathrm{Myr}$| (end of the simulation). Constant triaxiality T is marked with the three dashed lines. This provides a qualitative way to classify the resulting shapes for the various models based on their evolution and initial merger orbit. First of all, we notice that the clusters without prior merger are initially spherical and isotropic and maintain their shape (black markers in the top-right corner of Fig. 8) until the end of the simulation.

The eccentricity of the merging star cluster orbits plays an important role in determining the shape and kinematics of the remnant. Low e orbits carry more angular momentum, since |$\mathcal {L}=\mu \sqrt{G M a\left(1-e^2\right)}$| (where |$M=M_1+M_2$| and |$\mu =M_{1}M_{2}/M$|⁠), which is then transferred to the stars. Our merger simulations result in a narrow range of |$c/a$| with a clear trend of low (⁠|$e=0.05$|⁠) and medium eccentricities (⁠|$e=0.5$|⁠) towards less oblate systems (brown and orange markers in Fig. 8). Very eccentric (⁠|$e>0.9$|⁠) merger orbits on the other hand tend to scatter stars in radial orbits, resulting in more oblate and triaxial shapes as can be seen in Fig. 8 (yellow markers). We conclude that merger remnants have less spherical shapes than isolated clusters, a process driven by the eccentricity of the progenitor cluster orbits.

4.3 Star cluster kinematics

In a spherical coordinate system, the two angular velocity dispersion components |$\sigma _\mathrm{\theta }, \sigma _\mathrm{\phi }$| can be combined into a tangential velocity dispersion |$\sigma _\mathrm{t}$| as

The non-isotropic nature of the kinematic structure of the stellar system can be characterized by the velocity anisotropy parameter |$\beta$| (Binney & Tremaine 2008), defined as

The velocity anisotropy |$\beta$| is closely related to the orbital structure of a stellar population. When the stellar orbits are purely radial, then |$\sigma _\mathrm{t}=0$| and |$\beta =1$|⁠. Conversely, for a stellar population with exclusively circular orbits, |$\sigma _\mathrm{r}=0$| and |$\beta = -\infty$|⁠. The system is isotropic (⁠|$\beta =0$|⁠) if the two components |$\sigma _\mathrm{r} , \sigma _\mathrm{t},$| are equal.

To compute velocity dispersion profiles for our star cluster set-ups and to reduce the noise we construct average dispersion profiles. To do so we use 10 simulation snapshots from the last |$5\times 10^4 \, \mathrm{yr}$| in our runs. Each snapshot is centred in space and velocity to the region of highest stellar density using the shrinking sphere approach (Power et al. 2003). The velocity dispersion and anisotropy profiles (Fig. 9) of stars (thick solid lines) and of COs (thin dotted lines) of the remnants at and |$t=100 \, \mathrm{Myr}$| are presented in Fig. 9. The left panels correspond to simulations without MBHs. Both |$\sigma _\mathrm{r}$| and |$\sigma _\mathrm{t}$| have values of about |$\lesssim 10 \mathrm{~km} \mathrm{~s}^{-1}$| and decrease with increasing radius. The |$\beta \sim 0$| profiles indicate isotropic velocity dispersion in the inner parts of all clusters. The merger orbits do not change the dispersion profiles significantly. The presence of one (middle column) or two (right column) MBHs also does not change the kinematic structure significantly.

4.3.1 Rotation of cluster merger remnants

For the simulated cluster mergers, orbital angular momentum is transferred to internal angular momentum. The effect is strongest for low-eccentricity orbits (i.e. |$e=0.05,0.5$|⁠), which have the highest orbital angular momentum. This corresponds to higher merger remnant rotation velocities.

We present the two-dimensional line-of-sight (LOS) velocity (⁠|$V_{\mathrm{LOS}}$|⁠) maps of the cluster merger simulations after 100 Myr of evolution in Fig. 10. The maps are constructed utilizing the analysis package pygad (Röttgers et al. 2020). Following Naab et al. (2014), these maps are produced similar to the analysis of observational integral field unit data (see e.g. Naab et al. 2014, for more details). We use the clusters’ reduced moment of inertia tensor given by equation (25) to align the major axis of the cluster with the |$x-$|axis and perform the analysis within the central region of about |$\sim 2$| half-mass radii, corresponding to a spatial extent of |$4 \, \mathrm{pc}$|⁠. The velocity data are represented by spaxels of constant signal-to-noise ratio (a fixed number of stars per spaxel) using a Voronoi tessellation algorithm (Cappellari & Copin 2003). We also show contours of constant surface density.

The highly eccentric (e  = 0.99) merger orbits carry the least angular momentum and the respective remnants show no sign of ordered rotation (left panels in Fig. 10) independent of whether the clusters host BHs or not. The mergers with higher angular momentum orbits produce remnants that clearly show ordered rotation with velocities up to |$\sim$|3 |$\mathrm{km} \, \mathrm{s}^{-1}$| (middle and right panels of Fig. 10). In contrast to the velocity dispersion (see Section 4.3), the presence of MBHs in the clusters does not change the rotation features.

5 FORMATION AND EVOLUTION OF MBHBS

In this section we focus on the formation and evolution of the MBH binaries in the merger runs 2MBH_e005, 2MBH_e05, and 2MBH_e099 with two MBHs and the isolated cluster, 2MBH_iso, with a central MBH binary. The different cluster orbits affect the time when the MBH binaries become bound in the cluster mergers, while the 2MBH_iso MBH binary is bound by definition. For the 2MBH_e05 and 2MBH_e005 runs, the MBH binaries become bound at 1 and |$1.3 \, \mathrm{Myr}$|⁠, respectively. For the merger with the highest eccentricity of e  = 0.99 (2MBH_e099) the MBH binary becomes first bound at |$1.1 \, \mathrm{Myr}$| then unbound again at |$1.5 \, \mathrm{Myr}$|⁠. During the first bound phase the semimajor axis of the binary shrinks from |$a_{\mathrm{b}} \approx 1 \, \mathrm{pc}$| down to |$0.1 \, \mathrm{pc}$|⁠. The binary reaches a steady bound state at |$\sim 3 \, \mathrm{Myr}$| with |$a_{\mathrm{b}} \le 0.1 \, \mathrm{pc}$|⁠. Interactions with the stellar background lead to a continuous energy exchange with the MBH binary. This leads to orbital energy being removed from the binary, resulting in decrease of their semimajor axes |$a_{\mathrm{b}}$|⁠. When the MBH binary separation becomes lower than a critical value |$a_\mathrm{h}$| of the semimajor axis, the binary has reached the so-called ‘hard’ phase at Quinlan (1996)

where |$M_{\mathrm{2}}$| is the mass of the secondary MBH (in our case |$M_{1}=M_{2}=500 \, \mathrm{M_{\odot }}$|⁠) and |$\sigma$| is the local velocity dispersion inside the |$r_{\mathrm{SOI}}$| of the MBH binary. A typical value for the hard separation in our simulations is |$a_\mathrm{h} \approx 2\times 10^{-3} \, \mathrm{pc}$|⁠. The first time the MBH binaries in our simulations pass this threshold is at |$t_{\mathrm{hard}}=1.7, \, 4.0,\, 4.8$| and |$7.5 \, \mathrm{Myr}$| for the noMBH_iso, 2MBH_e005, 2MBH_e05, and 2MBH_e099 models, respectively. We note that on galactic scales, the eccentricity |$e_{\mathrm{b}}$| of formed bound SMBH binaries is strongly affected by stochastic effects (Rawlings et al. 2023, ).

5.1 Hardening and eccentricity growth of MBH binaries in a stellar background

The hardening of an MBH binary depends on the properties of its stellar environment, i.e. the velocity dispersion and mass density of the stars. For a fixed stellar background, the hardening rate H and eccentricity growth K are given by (Quinlan 1996)

where |$\rho$| and |$\sigma$| are the stellar density and velocity dispersion inside the SOI |$r_{\mathrm{SOI}}$| (equation 24) of the binary. Equation (30) assumes fixed values for |$\rho$| and |$\sigma$| and has been used in various studies on the efficiency of MBH hardening and its effect on the production of GW merger events and the production of RAs and HVs (Sesana et al. 2006, 2008, 2009; Leigh et al. 2017; Rasskazov et al. 2019; Gualandris et al. 2022; Evans et al. 2023). In our simulations, ongoing dynamical interactions of the MBH binaries with the stellar background change the average central stellar density and velocity dispersion continuously. In Fig. 11 we show that the density inside the SOI (top panel) of the MBH binaries continuously decreases with time, by a factor of 4–5 during the 100 Myr of evolution. Over the same time interval, the stellar velocity dispersion (bottom panel in Fig. 11) is decreased by |$4 \, \mathrm{km}\, \mathrm{s}^{-1}$| for the merger remnants. The decrease in central density is caused by the ejection of stars that have interacted with the central MBH binary whose semimajor axis is continuously shrinking for all simulations with MBH binaries (the top panel of Fig. 12). Individual strong encounters with the central MBH binary result in visible jumps (see e.g. the e = 0.5 merger, orange line, in the top panel of Fig. 12) in the semimajor axis evolution and result in the ejection of high-velocity objects (see discussion in Section 6).

The evolution of the eccentricity of the central MBH binary is not smooth (bottom panel of Fig. 12) and fluctuates significantly. This adds uncertainties on the further PN evolution of the MBH binary as the eccentricity has a strong impact on the process of coalescence through GW emission (e.g. Bonetti et al. 2020; Gualandris et al. 2022) and may also imprint distinct signatures on the number and velocity distribution of captured or ejected stars (Sesana et al. 2006, 2008; Rasskazov et al. 2019). We estimate the hardening rate coefficients H and K from the slope of linear fits on the |$1/a_{\mathrm{b}}$| and |$e_{\mathrm{b}}/\mathrm{ln}(1/a_{\mathrm{b}})$| time evolution. The values are given in Table 3. The coefficients H and K have been typically evaluated (e.g. Sesana et al. 2006) using scattering experiments assuming a fixed stellar background.

Eccentricity growth has been observed again in previous studies (Amaro-Seoane & Freitag 2006; Amaro-Seoane, Miller & Freitag 2009; Arca-Sedda & Mastrobuono-Battisti 2019) in the context of merging star and globular clusters hosting IMBHs. Its origin can be due to the collective effect of many stars interacting with the MBHB or through strong encounters with individual and/or binary stars (Askar, Davies & Church 2021).

5.2 GW-driven evolution and coalescence of MBH binaries

In this section we explore the evolution of our MBH binaries after |$t=100 \, \mathrm{Myr}$| where our simulations are stopped. To do so, we assume that the orbital elements evolve due to interactions of the MBH binary with the stellar background and GW emission that dominates at late times. The differential equations describing the evolution of semimajor axis and eccentricity are therefore written as

where the terms with |$\star$| refer to stellar interactions and |$GW$| to GW-emission, respectively. The former are described by Quinlan (1996),

where G is the gravitational constant, |$\rho$| and |$\sigma$| are, respectively, the stellar density and velocity dispersion inside the SOI |$r_{\mathrm{SOI}}$| (equation 24) of the binary as before, and H and K correspond to the fitted values for the hardening rate and eccentricity growth, respectively, given in Table 3. The GW-driven part is modelled via equation (9), where |$a,e$| correspond to the binary elements |$a_{\mathrm{b}},e_{\mathrm{b}}$|⁠.

The initial values used for solving equation (31) correspond to the values at the end of the simulation, i.e. |$(a_{\mathrm{100}},e_{\mathrm{100}}, \rho _{\mathrm{100}},\sigma _{\mathrm{100}})$|⁠. The results are presented in Fig. 13 (solid lines), where all MBH binaries merge before |$t< 4 \, \mathrm{Gyr}$|⁠. Although |$e_{\mathrm{100}}$| has some fixed value for each simulation but due to the stochastic nature of strong encounters, this value would vary if the simulation runs were slightly longer. For this reason and in order to capture the effect of |$e_{\mathrm{b}}$| in the full range of possible values, we integrate system equation (31) for an almost circular and a highly eccentric binary, whereas |$e_{\mathrm{100}}$| from the actual runs serves as intermediate ones. Assuming a fixed high initial eccentricity |$(e_{\mathrm{b}}=0.9)$| the binaries merge earlier, in less than 2 Gyr. For almost circular orbits |$(e_{\mathrm{b}}=0.01)$| the MBH binaries can take up to |$\sim$|5 Gyr to merge. The above estimates indicate that the MBHs formed in mergers of dense stars clusters will merge within a Hubble time. The main takeaway from this experiment is that given a sufficiently dense environment, hard MBH binaries could potentially merge regardless of how eccentric they are and can be observed with future space-born and ground-based GW detectors like Laser Interferometer Space Antenna (LISA; Amaro-Seoane & Freitag 2006; Amaro-Seoane et al. 2017; Arca-Sedda & Mastrobuono-Battist 2019) and Einstein Telescope (Maggiore et al. 2020).

6 THE EJECTION OF STARS AND COS

6.1 Classification and overview

In each simulation, we monitor the stars and COs that become unbound and eventually escape the cluster at a fixed distance of |$r_\mathrm{esc}=100 \, \mathrm{pc}$|⁠. As outlined in Section 2, the various dynamical ejection channels have characteristic ejection velocities. In this section, we give an overview of the properties of ejected stars and COs in our simulations and connect them to the most likely ejection mechanisms. We classify the ejected objects as follows:

• Low-velocity (LV) escapers (⁠|$v_{\mathrm{ej}} < 50 \mathrm{~km} \mathrm{~s}^{-1}$|⁠), originating mostly from weak two-body encounters and two-body relaxation in the vicinity of an MBH.

• RASs and COs (⁠|$50 \mathrm{~km} \mathrm{~s}^{-1} \le v_{\mathrm{ej}} < 300 \mathrm{~km} \mathrm{~s}^{-1}$|⁠), most likely from three- or four-body encounters of single or binary stars with BH+star and BH+BH or MBH+star and MBH+MBH binaries and/or binary encounters with a single MBH (a Hills-like mechanism) or a binary MBH.

• Hyper-runaway (HR) stars and COs (⁠|$300 \mathrm{~km} \mathrm{~s}^{-1} \le v_{\mathrm{ej}} < 700 \mathrm{~km} \mathrm{~s}^{-1}$|⁠), which are the extreme version of the previous channel from interactions with hard binaries.

• HVSs and COs (⁠|$v_{\mathrm{ej}} \ge 700 \mathrm{~km} \mathrm{~s}^{-1}$|⁠), most likely due to strong three- and four-body interactions with a single and/or binary MBHs.

• Extreme-velocity stars and COs (⁠|$v_{\mathrm{ej}} \ge 1000 \mathrm{~km} \mathrm{~s}^{-1}$|⁠), which are a handful of extreme and rare cases for which we suggest the most likely scenario of ejection.

In Fig. 14 we show the evolution of the cumulative number of escapers (top panel) and the escape rate (bottom panel) for all simulations. The characteristic time-scale for escaping our clusters can be estimated from the escaper removal distance and the cluster centre escape velocity yielding 100 pc/|$22 \mathrm{~km} \mathrm{~s}^{-1}$|  |$\sim 4.4$| Myr. The total number of escapers rises rapidly after 10–15 Myr (the post merger phase) when the first stars cross the 100 pc boundary. After 100 Myr about 1000 stars have escaped for all simulations without MBHs. The escape rates (bottom left panel) are almost constant (though a moderate peak develops for the most circular merger) at |$\sim 0.2/\mathrm{Myr}$|⁠. For simulations with one MBH (middle panels in Fig. 14) the total number of escapers (⁠|$\sim$|2000) is about a factor of 2 higher compared to the run without MBHs, and the escape rates peak at |$\sim 1.0/\mathrm{Myr}$| after |$\sim 20 \, \mathrm{Myr}$|⁠. At the end of the simulation the escape rate drops to similar values as the no MBH simulations. In the presence of two MBHs (right panels) the total number of escapers rises to |$\gtrsim$|3000 objects and the escape rates peak earlier (due to the higher speeds of the escapers) at values of |$\sim 3.0/\mathrm{Myr}$|⁠. After |$100 \, \mathrm{Myr}$| the rates have dropped to similar values as the simulations with one and no MBH. Even though the escape rates are similar towards the end of the simulations, the velocities of the objects are significantly higher in the presence of MBHs. We will discuss more about this fact in the following text.

The number |$N{\mathrm{esc}}$| and rate |$\Gamma _{\mathrm{esc}}$| of ejections do not strongly correlate with the orbital eccentricity of the mergers. The results indicate a weak dependence to our ICs, thus suggesting that our escaper numbers and rate estimates are robust and a generic feature of star cluster mergers with and without MBHs. For clusters without MBHs, only the noMBH_e005 remnant shows an increased number and rate of |$1/\mathrm{Myr}$| of escapers in the first |$15 \mathrm{Myr}$| compared to the (⁠|$\Gamma _{\mathrm{esc}}<0.5$| ejection per |$\mathrm{Myr}$|⁠) other runs, which could be due to the higher number of dynamically formed binaries during that time (Section 6.2.1). When MBHs are present, the merger orbit itself does not significantly affect the overall number and rates of ejection, as we see from the bottom panels (middle, right) in Fig. 14, where |$\Gamma _{\mathrm{esc}}$| from the isolated runs is similar to those of the merger remnants.

The velocity distribution of the escapers in all simulations is shown in Fig. 15, separated into stars (top panels), stellar BHs (middle panel), and WD bottom panels. The assumed velocity limits for RASs, hyper-RASs, and HVSs are indicated. The total numbers of escaping stars, BHs, and WDs for these three classes are given in the first four columns of Table 4 and the total mass and fraction of the respective population lost by escapers is given in the five columns on the right. For the isolated cluster without MBHs all escapers are LV escapers and we do not find stars with velocities in excess of the assumed RAS limit of 50 km s|$^{-1}$|⁠. Merger of clusters without MBHs can produce a small population of RASs (⁠|$\lesssim 10$|⁠) and, in rare cases, an hyper-runaway BH. Already in the presence of one MBH the isolated and merger simulations produce a sizable population of |$\sim 250$| runaway objects, up to five hyper-RASs. The presence of an MBH also results in the ejection of a significant fraction of the initial BH population (from |$\sim$|22 to |$\sim$|37 per cent), some with runaway velocities, and a few per cent of the stellar, WD, and neutron star population in the runaway velocity regime. If both clusters have an MBH that is forming a binary in the remnants, the runaway population increases to |$\sim 800$| and up to 45 HR and three HVSs can be formed. The number of the escaping BHs increases slightly. These results clearly indicate the existence of a single or binary MBH in a star cluster will create a new population of runway stars and COs.

In Fig. 16 we show the distribution of ejection velocities for objects crossing a distance of 100 |$\mathrm{pc}$| as a function of time for all simulations with binary MBHs. The results are not sensitive to the choice of 100 |$\mathrm{pc}$| as the particle removal radius, since this would only slightly affect (increased/decreased for lower/higher values) the low-end of the the |$v_{\mathrm{ej}}$| distribution. After the central MBHs become a hard binary (horizontal red dashed lines) objects are ejected with velocities exceeding 50 km s|$^{-1}$|⁠, which is the velocity limit for runway stars. The ejection velocity of a star ejected via a strong encounter with an MBH binary typically scales with the orbital velocity of the binary |$v_{\mathrm{b}}$| (Valtonen & Karttunen 2006; Merritt 2013)

where |$M_{\bullet }=M_{1}+M_{2}$|⁠. At later times, when the central MBH binary becomes harder with |$a_{\mathrm{b}} \lesssim 10^{-3} \, \mathrm{pc}$| (see Fig. 12) the maximum escape velocities increase into the HR velocity regime. The lower velocity boundary for data points in this figure is determined by the flight times to the 100 pc particle removal distance, with the particle ejection velocities from the central parts of the cluster.

6.2 Ejection mechanisms

6.2.1 Dynamically formed binary systems

Single–binary and binary–binary interactions play a central role in the ejections of stars from their host clusters. Besides from being formed in multiples, stars may form binary systems via dynamical processes, especially in dense cluster environments. These dynamical processes are present in our simulations and they include TCs (Press & Teukolsky 1977; Giersz 1986; Generozov et al. 2018; Rizzuto et al. 2023) or three-body encounters known as three-body binary formation (Mansbach 1970; Aarseth & Heggie 1976; Goodman & Hut 1993; Portegies Zwart & McMillan 2000; Heggie & Hut 2003; Atallah et al. 2024; Ginat & Perets 2024). We briefly discuss the dynamical binary formation in Appendix  A.

We show the number of dynamically formed binaries for snapshots with 4 Myr separation in Fig. 17. The ICs do not contain binary stars or COs. All simulations form stellar (low-mass) binaries. Initially |$(t< 4\, \mathrm{Myr})$| up to 25 such systems form and the numbers drop to up to 10 later on throughout the runs (top panel) without a clear trend with ICs. For stellar BH+star binaries (second panel) the trend is similar with slightly lower (about 10 less on average) numbers. For stellar BH+BH binaries (third panel, Fig. 17) there is a clear trend. Simulations with central MBHs have at most two such systems. Meanwhile, the simulations without central MBHs can have up to 10 BH–BH binaries. In these simulations, more BHs escape and binaries can be disrupted by the interaction with the central MBHs, which leaves one BH directly bound to the central MBHs (see bottom panel). After |$\sim$|Myr up to 10 stars are bound to the single MBHs and up to three to the binary MBHs (fourth panel). This trend also holds for MBH–BH binaries (bottom panel, Fig. 17). Binaries with white dwarfs WD+WD or WD+star do form occasionally with an average number of |$N_{\mathrm{WD+WD}} \approx N_{\mathrm{WD+star}} \approx 1-2$|⁠, but most of them do not survive due to multiple dynamical encounters and the fact that they form in relatively wide configurations.

6.2.2 Ejections from clusters without MBHs

First we examine the ejections from clusters without an MBH. In Fig. 14 we see these clusters mostly produce LV escapers (⁠|$v_{\mathrm{ej}} \lesssim 10 \mathrm{~km} \mathrm{~s}^{-1}$|⁠) corresponding to two-body relaxation and weak encounters between single stars. Still, there is a noteworthy population of ejections above |$10 \mathrm{~km} \mathrm{~s}^{-1}$|⁠, which could be the result of encounters of single stars/COs with star+star, star+BH, or BH + BH binaries (see Fig. 17). An encounter with a star+star binary leads to velocity kicks of the same order as two-body relaxation-driven escapers, whereas encounters with star+BH or BH + BH binaries can explain the higher velocities. The range of velocities for these ejection processes strongly depends on the binary semimajor axis and mass ratio q. For example, for a hard binary with |$a=2 \, \mathrm{au}$| and |$q=1/25$|⁠, equation (7) implies a velocity kick of |$v_{\mathrm{ej}} > 100 \mathrm{~km} \mathrm{~s}^{-1}$|⁠. From Fig. 14 we notice that the noMBH_e005 remnant is able to produce a few RAs even in the absence of an MBH as there is a sufficient number of BH+star and BH + BH binaries present in the simulations (Fig. 17) to fuel the few-body ejection channels. Especially the HR (⁠|$v_{\mathrm{ej}} \approx 473 \mathrm{~km} \mathrm{~s}^{-1}$|⁠) object, which is a BH (see middle-left panel in Fig. 14), is ejected at |$20\, \mathrm{Myr} < t < 25\, \mathrm{Myr}$|⁠. At this time the star cluster remnant contains |$N_\mathrm{BH+star}=15$| and |$N_\mathrm{BH-BH}=10$|⁠, which is the largest number of those types of binaries throughout the run. As such, the ejection time coincides with the most probable moment for a strong few-body ejection.

6.2.3 Ejections from clusters with a single MBH

Next we discuss ejections from clusters and remnants with a single MBH. In addition to the ejection channels in the case without any MBHs, now there is a significant fraction of MBH+star and MBH+BH binaries. The numbers of those binaries vary from |$N_{\mathrm{MBH+star}}\approx 50$| at early times to |$N_{\mathrm{MBH+star}} \gtrsim 10$| at later times and |$N_{\mathrm{MBH+BH}}\approx 3\!-\!6$| (Fig. 17) and essentially generate ejections covering the full range of |$v_{\mathrm{ej}}$| (given a broad distribution of semimajor axes and mass spectrum of stars and BHs). The Hills mechanism also contributes and can explain some rare but strong encounters of binaries with the MBH. The number of stellar and BH binaries is significantly lower compared to the simulations without MBHs. The MBH+star and MBH+BH channel is likely dominant in this case. The number of Hills-like encounters is low due to the lack of available dynamically formed stellar-mass binaries. The number of RAs lies in the range of |$192 < N_\mathrm{RAs} < 280$| for all cases, with a handful |$(0\!-\!5)$| of HRs, and |$0\!-\!2$| HVs. We see from Fig. 14 that none of them is a BH. Most HRs and stars with higher |$v_{\mathrm{ej}}$| are produced in the 1MBH_e005 run, which has a larger number of MBH+star and MBH+BH binaries (Fig. 17), implying stellar encounters with these types of binaries.

A notable (classified as ‘extreme’ and not shown in Fig. 14) ejection is one at |$t \approx 87 \, \mathrm{Myr}$| of the 1MBH_e005 run with |$v_{\mathrm{ej}} > 1000 \mathrm{~km} \mathrm{~s}^{-1}$|⁠. Specifically, a star of |$m_\mathrm{\star }=0.7 \, \mathrm{M_{\odot }}$| accompanied by a |$m_{\bullet }=14 \, \mathrm{M_{\odot }}$| BH, where the two used to be members of a bound binary system before the ejection. An MBH and star+BH interaction, i.e. a Hills ejection, is the most probable reason for this escape event. The binary semimajor axis in the last snapshot about |$1000 \, \mathrm{yr}$| before the ejection was |$a \approx 27 \, \mathrm{au}$|⁠, which gives an estimated ejection velocity using equation (6) of |$v_{\mathrm{ej}} \approx 900 \mathrm{~km} \mathrm{~s}^{-1}$|⁠. This ejection velocity, close to |$1000 \mathrm{~km} \mathrm{~s}^{-1}$|⁠, is very well consistent with the Hills mechanism. We highlight the potential importance of the Hills mechanism for the production of high-velocity stars in the context of star clusters and BHs in the IMBH mass regime as the Hills mechanism is commonly associated with galactic nuclei and their SMBHs.

6.2.4 Ejections from clusters with a binary MBH

Finally, we focus on ejections from clusters with a binary MBH at their centre. We note that although all the previous ejection channels are possible in those clusters, the vast majority does not involve interactions with stellar or BH binaries, since their numbers are very low (see Fig. 17). The dominant additional ejection mechanism here is the interaction of an MBH+MBH binary with a star or a stellar BH. The MBH binaries become hard, i.e. |$a_{\mathrm{b}} \le a_{\mathrm{h}}$|⁠, already in the early |$t < 5 \, \mathrm{Myr}$| post-merger phase as we discussed in Section 5. The MBH binaries efficiently produce additional RAs, HRs, and even HVs. Fig. 16 marks (red dashed lines) the time at which the MBH binaries become hard, while the different horizontal (black) lines correspond to the various ejection velocities (⁠|$v_{\mathrm{ej}}$| is the velocity of the ejected object at |$r_\mathrm{esc}=100 \, \mathrm{pc}$|⁠) classes. High-velocity ejections are initiated only after |$a_{\mathrm{b}} = a_{\mathrm{h}}$|⁠. With enough available stars and COs in the vicinity of the MBH binary (top-right panel in Fig. 5), the number of ejections rises from |$1830 < N_{\mathrm{esc}} < 2450$| for the single MBH clusters to |$3600 < N_{\mathrm{esc}} < 4800$| for the ones with an MBHB (top panels in Fig. 14).

7 SUMMARY AND CONCLUSIONS

Star clusters in galaxies are expected to frequently interact and occasionally merge. This can happen during their hierarchical assembly or for evolved star clusters in the spiral arms of galactic discs, in galactic haloes, or in GCs. In this study we explore the consequences of the merging star clusters hosting MBHs for the cluster evolution and the production of escaping stars and COs. With |$M_\bullet =500 \, \mathrm{M}_\odot$| the MBHs have masses about one order of magnitude higher than the most massive observed massive stellar BH (see e.g. Gaia Collaboration 2024). Such low-mass MBHs might naturally form by rapid stellar collisions in hierarchically forming young and dense star clusters (see e.g. Gieles & Portegies Zwart 2011; Rantala et al. 2024) and be retained in the clusters at later times.

We have performed a suite of simulations of mergers of star clusters with individual cluster masses of |$M_{\star } = 2.7 \times 10^4 \, \mathrm{M}_{\odot }$| and |$N=64\,000$| individual stars and COs. In one set of simulations, each of the merging star clusters has a central MBH of |$M_{\bullet } = 500 \, \mathrm{M}_\odot$|⁠. For comparison we perform simulations without initial MBHs, and with only one MBH of |$M_{\bullet } =1000 \, \mathrm{M}_{\odot }$|⁠. We merge the clusters on a circular, a very radial, and an intermediate eccentricity orbit. To estimate the effect of the merger process itself we also simulate isolated clusters with |$M_{\star } = 5.4 \times 10^4 \, \mathrm{M}_{\odot }$| without any MBHs, with one MBH (1000 |$\mathrm{M}_{\odot }$|⁠), and including a binary (2|$\times 500 \, \mathrm{M}_\odot$|⁠) MBH. The simulations are performed with the hierarchical forward fourth-order GPU-accelerated N-body code bifrost (Rantala et al. 2023), which includes a regularization scheme for binaries and PN dynamics up to order PN3.5. We study the impact of the presence of MBHs on the star cluster merger remnants. Specifically we study how they affect the (i) kinematic and structure of the remnant clusters, (ii) the formation and evolution of MBH binaries and their coalescence time, and (iii) the production and populations of ejected stars and COs.

The clusters merge rapidly and the merger remnants are almost spherical with isotropic velocity dispersions. The two more circular merger orbits with |$e=0.05$| and |$e=0.5$| result in remnant star clusters with rotational velocities of |$v_{\mathrm{LOS}}=v_{\mathrm{rot}}\sim 3 \mathrm{km}\, \mathrm{s^{-1}}$|⁠, similar to observations (e.g. Bellazzini et al. 2012; Fabricius et al. 2014) who find a few |$\mathrm{km}\, \mathrm{s^{-1}}$|⁠. Apart from a slightly reduced central density of the order of |$\sim 300\!-\!700 \, \mathrm{M_{\odot }\, \mathrm{pc^{-3}}}$|⁠, we find no strong evidence for a measurable/observable impact of the MBHs on the structural and kinematic properties of the merger remnants themselves. In conclusion, it is difficult to probe the presence of MBHs in merger remnants based on their kinematic properties alone.

In the cluster merger remnants the sinking MBHs rapidly form binaries and harden by interactions with stars and COs. This process produces a new population of escaping stars and COs with velocities |$\gtrsim$|50 km s|$^{-1}$|⁠, which is absent in the star cluster mergers without MBHs. Within 100 Myr, |$\sim$|800 stars with |$v_{\mathrm{ej}}$|  |$\gtrsim$| 50 km s|$^{-1}$| are ejected and would be classified as RASs. In addition about 30 stellar BHs escape with similar velocities within 100 Myr. About 30 stars can be accelerated to high velocities |$\sim$|300 km s|$^{-1}$|⁠. On average |$\sim$|3000 stars escape at velocities lower than 50 km s|$^{-1}$|⁠. Overall, the remnants lose |$\sim$|30 per cent of their BH population and |$\sim$| 3 to 4 per cent of their WD and star population if MBHs are present.

In the absence of MBHs the fraction of escaping BHs drops to |$\sim$|6 per cent and to below 1 per cent for WDs and stars. Comparison simulations of isolated clusters of the same mass as the merger remnants and initialized with central MBH binaries as well as star cluster mergers without MBHs show that the high-velocity ejection process is driven by the MBH binaries and not the cluster merger process itself. For merger simulations without MBHs about |$\sim$|1000 stars and |$\sim$|10 BHs escape the system at velocities below 50 km s|$^{-1}$|⁠. Comparisons with a single MBH of 1000 |$\mathrm{M}_\odot$| in isolated or merging cluster show a smaller population of runway stars at lower velocities. Recently, Polak et al. (2024) have highlighted the role of star cluster mergers for the formation of RASs. Their simulations feature a more realistic hydrodynamical hierarchical set-up and do not include MBHs. When compared to our idealized star cluster mergers without MBHs we find that one of their ‘runaway’ star group with a velocity of |$\sim$| 39 km s|$^{-1}$| is consistent with our simulations. As we assume a runaway velocity limit of 50 km s|$^{-1}$| they would fall in our LV escaper regime. The second runway group in Polak et al. (2024) has a velocity of |$\sim$|87 km s|$^{-1}$|⁠. We do not reach such high velocities by cluster mergers alone. As discussed in Polak et al. (2024), the more realistic and more dynamical set-up in their simulation might allow for higher escape velocities during cluster mergers.

Based on the evolution of the structural and dynamical cluster properties and the MBH binary hardening rates we expect the binary MBHs in our simulations to merge in less than a Hubble time and produce observable GW emission events, detectable by future GW detectors like Advanced LIGO, the Einstein Telescope, or the LISA. The gravitational recoil kicks exerted on the MBH merger remnants might easily eject the MBHs from the star clusters leaving no detectable evidence for their prior existence apart from a strongly reduced stellar BH sub-population.

Low-mass MBHs might form naturally by stellar collisions in dense stars clusters (see e.g. Portegies Zwart et al. 2004; Rizzuto et al. 2021; Rantala et al. 2024, and references therein), in particular in the dense star-forming low-metallicity interstellar medium in the early Universe. Mergers of these clusters are very likely during their hierarchical formation, during NSC assembly, or in the discs of early compact galaxies. In this case, the results of our study imply that interactions with low-mass MBH binaries formed in merging stars clusters are an important additional channel for producing runaway and high-velocity stars as well as free floating stellar BHs and COs, which are produced at high redshift and still populate the discs and haloes of present-day galaxies.

ACKNOWLEDGEMENTS

TN acknowledges the support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311 of the DFG Cluster of Excellence ‘ORIGINS’. Software: bifrost (Rantala et al. 2023), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), pygad (Röttgers et al. 2020), seaborn (Waskom 2021), pandas (The pandas development Team 2023).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: DYNAMICALLY FORMED BINARIES IN STAR AND GLOBULAR CLUSTERS

The increased densities during the core-collapse phase of star clusters elevate the rate of dynamical encounters, leading to the formation of binaries through three-body interactions. The encounter rate |$\bar{C}$| of three initially unbound bodies leading to binary formation in a system of equal mass objects |$m_*$| is (Goodman & Hut 1993; Binney & Tremaine 2008)

where n is the number density and |$\sigma$| the velocity dispersion of the system. Since the above rate refers to equal-mass bodies only, it is however still uncertain how well it predicts the binary formation rates in multicomponent environments such as star and globular clusters and/or galactic nuclei. Atallah et al. (2024) found that in such encounters, it is highly unlikely that the two most massive bodies become bound, in contrast to what has mostly been suggested in the literature so far. The large number of single bodies in the system leads to frequent fly-by and exchange interactions with the formed binaries, where one of the members is replaced by a field star. The higher the mass of the incoming body, the higher the probability of it replacing a binary member (Valtonen & Karttunen 2006). This inevitably results to massive stars becoming members of binaries.

For a system containing stellar-mass BHs the above mechanisms allow the efficient formation of binary BHs (Portegies Zwart & McMillan 2000; Park et al. 2017; Kremer et al. 2019; Torniamenti et al. 2022; Kritos et al. 2024,; Rantala et al. 2024), potentially leading to their GW-driven coalescence, covering a broad mass-ratio spectrum, detectable with current (Abbott et al. 2016, 2019) and future (Amaro-Seoane et al. 2017) GW detectors. Additionally, BH-star binaries can be formed through exchange in three- or four-body encounters (e.g. Ryu et al. 2023b) in massive star clusters (e.g. Rastello et al. 2023,; Marín Pina et al. 2024; Fantoccoli et al. 2024), potentially explaining the three dormant BHs detected by Gaia (Gaia Collaboration 2016): Gaia BH1 (El-Badry et al. 2022), BH2 (El-Badry et al. 2023), and BH3 (Gaia Collaboration 2024).