# FROST-CLUSTERS – I. Hierarchical star cluster assembly boosts intermediate-
mass black hole formation

Antti Rantala1, Thorsten Naab1, Natalia Lahén1
1Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748,
Garching, Germany
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Observations and high-resolution hydrodynamical simulations indicate that
massive star clusters assemble hierarchically from sub-clusters with a
universal power-law cluster mass function. We study the consequences of such
assembly for the formation of intermediate-mass black holes (IMBHs) at low
metallicities ($Z=0.01\;Z_{\mathrm{\odot}}$) with our updated N-body code
BIFROST based on the hierarchical fourth-order forward integrator. BIFROST
integrates few-body systems using secular and regularized techniques including
post-Newtonian equations of motion up to order PN3.5 and gravitational-wave
recoil kicks for BHs. Single stellar evolution is treated using the fast
population synthesis code SEVN. We evolve three cluster assembly regions with
$N_{\mathrm{tot}}=1.70$–$2.35\times 10^{6}$ stars following a realistic IMF in
$\sim$1000 sub-clusters for $t=50$ Myr. IMBHs with masses up to
$m_{\bullet}\sim{2200}\>\mathrm{M_{\odot}}$ form rapidly mainly via the
collapse of very massive stars (VMSs) assembled through repeated collisions of
massive stars followed by growth through tidal disruption events and BH
mergers. No IMBHs originate from the stars in the initially most massive
clusters. We explain this by suppression of hard massive star binary formation
at high velocity dispersions and the competition between core collapse and
massive star life-times. Later the IMBHs form subsystems resulting in
gravitational-wave BH-BH, IMBH-BH and IMBH-IMBH mergers with a
$m_{\bullet}\sim{1000}\>\mathrm{M_{\odot}}$ gravitational-wave detection being
the observable prediction. Our simulations indicate that the hierarchical
formation of massive star clusters in metal poor environments naturally
results in formation of potential seeds for supermassive black holes.

###### keywords:

gravitation – celestial mechanics – methods: numerical – galaxies: star
clusters: general – stars: black holes

††pubyear: 2023††pagerange: FROST-CLUSTERS – I. Hierarchical star cluster
assembly boosts intermediate-mass black hole formation–B.3.4

##  1 Introduction

The runaway stellar collision scenario in dense, massive star clusters has
been proposed to explain the emergence of massive black hole (MBH) seeds at
high redshifts (Gold et al., 1965; Spitzer & Stone, 1967; Spitzer & Hart,
1971; Peebles, 1972; Begelman & Rees, 1978). If the MBH seeds form but later
fail to further grow in mass, they should be in principle observable in the
present-day Universe as intermediate-mass black holes (IMBHs111In this work we
divide the black hole (BH) population at $Z=0.01\>Z_{\mathrm{\odot}}$ as
follows: stellar-mass BHs with $m_{\bullet}\lesssim{63}\>\mathrm{M_{\odot}}$,
stellar-mass pulsational pair-instability mass-gap BHs in the range of
${63}\>\mathrm{M_{\odot}}\lesssim
m_{\bullet}\lesssim{100}\>\mathrm{M_{\odot}}$, IMBHs with
${100}\>\mathrm{M_{\odot}}\lesssim
m_{\bullet}\lesssim{10^{5}}\>\mathrm{M_{\odot}}$ including the mass gap IMBHs
at ${100}\>\mathrm{M_{\odot}}\lesssim
m_{\bullet}\lesssim{290}\>\mathrm{M_{\odot}}$, and finally MBHs above
$m_{\bullet}\gtrsim{10^{5}}\>\mathrm{M_{\odot}}$. We use the terms IMBH and
MBH seed interchangeably.).

The runaway collision scenario for MBH seed and IMBH formation has been widely
studied using numerical simulations (see e.g. the reviews of Mezcua 2017 and
Greene et al. 2020 and references therein). The so-called slow formation
channel (Miller & Hamilton, 2002) has a stellar mass BH with
$m_{\bullet}>{50}\>\mathrm{M_{\odot}}$ growing in a globular cluster up to
masses of $m_{\bullet}\sim{10^{3}}\>\mathrm{M_{\odot}}$ in 10 Gyr, mainly by
merging with lower-mass BHs. The bottleneck for the slow collision channel is
the initial formation of a stellar mass BH with
$m_{\bullet}>{50}\>\mathrm{M_{\odot}}$. Due to the physics of pair instability
and pulsational pair-instability supernovae, hereafter (P)PISN, it is expected
that no BHs with masses in the range ${50}\>\mathrm{M_{\odot}}\lesssim
m_{\bullet}\lesssim{130}\>\mathrm{M_{\odot}}$ should from through single
stellar evolution at solar metallicity (Fowler & Hoyle, 1964; Woosley et al.,
2007; Woosley, 2017). At lower metallicities this (P)PISN mass gap widens and
moves into somewhat higher BH masses.

The so-called fast runaway collision scenario relies on a very massive star
(VMS222We refer any stars more massive than
$m_{\star}>{150}\>\mathrm{M_{\odot}}$ as VMSs in this work.) growing via
repeated mergers to sufficiently high masses (Lee, 1987; Quinlan & Shapiro,
1990; Bonnell et al., 1998) to produce an IMBH remnant above the (P)PISN mass
gap range. The IMBH may further grow by accreting stellar material in tidal
disruption events (TDEs) and by merging with other BHs. Numerical simulations
of the fast scenario (e.g. Portegies Zwart et al. 1999; Portegies Zwart &
McMillan 2002; Gürkan et al. 2004; Portegies Zwart et al. 2004; Freitag et al.
2006a, b; Rizzuto et al. 2021; Rizzuto et al. 2022; Arca Sedda et al. 2023a;
Arca Sedda et al. 2023c, b) show that the stellar collisions can result in a
formation of a VMS or a supermassive star
($m_{\star}>{10^{3}}\>\mathrm{M_{\odot}}$, e.g. Denissenkov & Hartwick 2014;
Gieles et al. 2018) with a mass from hundreds to potentially tens of thousands
of solar masses. In a true runaway collision sequence, each collision further
increases the collision rate. Rapid growth is however possible even without
fulfilling this criterion. Typically only $1$–$2$ stars per cluster can grow
via collisions (Baumgardt & Klessen, 2011). In the global instability variant
of the fast runaway collision scenario, stellar collisions become relevant if
the collision time-scale is shorter than the characteristic age of the system
(Escala, 2021; Vergara et al., 2023).

There are other ingredients in the runaway VMS scenario in addition to the
star cluster masses, densities and metallicities. Primordial mass segregation
may play a role in the collision dynamics but does not necessarily increase
collision rate (Ardi et al., 2008). An external potential may delay the
cluster evolution and the onset of collisions by increasing the velocities of
the stars (Reinoso et al., 2020). Gas may be very important in the growth of
VMSs (Clarke & Bonnell, 2008; Baumgardt & Klessen, 2011; Moeckel & Clarke,
2011; Dale et al., 2015a), while the initial mass function (IMF) of the
cluster (Goswami et al., 2012) and the primordial binary fraction also likely
play a role (Gürkan et al., 2006) especially due to single-binary interactions
(Fregeau et al., 2004; Gaburov et al., 2008). And early variant of the VMS
growth channel may have occurred in the first clusters (Boekholt et al., 2018;
Reinoso et al., 2018) of Pop-III stars (Schaerer, 2002).

The fast and the slow collisional IMBH formation scenarios are not mutually
exclusive. Even if the stellar collisions do not result in the formation of an
IMBH via the fast VMS formation and collapse channel, repeated tidal
disruption events or BH-BH collisions may still result in a formation of mass
gap BHs or IMBHs (Morscher et al., 2015; Rodriguez et al., 2016; Askar et al.,
2017; Banerjee, 2018; Fragione & Kocsis, 2018; Antonini et al., 2019;
Rodriguez et al., 2019; Belczynski & Banerjee, 2020; Arca-Sedda et al., 2021;
Banerjee, 2021; Rizzuto et al., 2021; Mapelli et al., 2021; Rizzuto et al.,
2022; Rizzuto et al., 2023; Arca Sedda et al., 2023c). IMBH (or MBH seed)
formation via the fast collisional runaway channel most likely occurs in
young, massive, star clusters (Portegies Zwart et al., 2010) with low
metallicity. First, low metallicity environments are favored due to weaker
stellar winds of massive stars resulting in more massive IMBH remnants (e.g.
Mapelli 2016; Di Carlo et al. 2021). There is still uncertainty about the wind
mass-loss rates of massive stars. High wind mass-loss rates (especially at
high metallicities) might suppress the VMS IMBH formation channel or inhibit
it altogether (e.g. Belkus et al. 2007; Yungelson et al. 2008; Glebbeek et al.
2009; Pauldrach et al. 2012).

Next, the massive star clusters need to be dense to have high enough collision
rates. The collision rates of massive stars are strongly enhanced by the mass
segregation and core collapse of the star clusters. Due to the negative heat
capacity of gravitating systems, the centers of star clusters may undergo a
catastrophic core collapse (Lynden-Bell & Wood, 1968; Lynden-Bell & Eggleton,
1980; Sugimoto & Bettwieser, 1983; Bettwieser & Sugimoto, 1984; Heggie, 1993;
Spurzem & Aarseth, 1996; Heggie & Hut, 2003; Kamlah et al., 2022) which is
eventually halted by the dynamical formation of binary stars that act as an
additional energy source at the cluster core. For clusters with a stellar mass
function, mass segregation proceeds on a dynamical friction time-scale
$t_{\mathrm{df}}$ of the most massive stars, and the core collapse time-scale
is only a small fraction of the two-body relaxation time-scale of the cluster.
The core bounce becomes ambiguous in low-mass clusters when the mass of the
most massive star exceeds $0.1\%$ of the cluster mass (Fujii & Portegies
Zwart, 2014). The mass segregation and core collapse of star clusters increase
the stellar collision rate by increasing the number density of massive stars
at their centers. Only $0.1\%$ of the total cluster mass may collapse into the
very center (Gürkan et al., 2004), which is far below the typical definitions
for the core size, density and mass (Casertano & Hut, 1985).

Dense, low-metallicity star clusters have been studied observationally both in
the local Universe and at high redshift. Massive star clusters are
characterized by their flat constant-density cores with core radii
$r_{\mathrm{c}}$ and core densities $\rho_{\mathrm{c}}$ while the outer parts
are described by power-laws optionally including a tidal cut-off (Plummer,
1911; King, 1962; Casertano & Hut, 1985; Elson et al., 1987). In the local
Universe, the Magellanic clouds offer the best-resolved view of the internal
structure of young low-metallicity star clusters (Hunter et al., 2003). The
core-radii and half-mass radii are often less than a parsec (Hill & Zaritsky,
2006; Werchan & Zaritsky, 2011; Sun et al., 2017; Gatto et al., 2021). Extra-
galactic surveys, that offer a larger sample size but limited resolution, have
also reported extremely compact young star clusters (Brown & Gnedin, 2021),
although with a significant scatter in the mass-size relation. The Tarantula
nebula is an extreme local example of the complex hierarchical structure of
young star cluster-forming regions: sub-clusters range from 30 Myr old
clusters in the outer regions (Cignoni et al., 2016) to the 1-2 Myr old core
(Crowther et al., 2010) to still ongoing star formation (Kalari et al., 2018).
The dense core region R136 (Selman & Melnick, 2013) contains several young
stars more massive than $m_{\star}\gtrsim{120}\>\mathrm{M_{\odot}}$ (Massey &
Hunter, 1998) while even higher stellar masses of
$m_{\star}\sim{160}\>\mathrm{M_{\odot}}$–${320}\>\mathrm{M_{\odot}}$ have been
suggested (Crowther et al., 2010). Outside the cluster, massive runaway star
candidates can be observed (Bestenlehner et al., 2011; Gvaramadze &
Gualandris, 2011). For further details about the massive stars in the region
see Vink et al. (2015) and references therein.

Low-mass galaxies that host proto-globular clusters, observed with the aid of
gravitational lensing (Livermore et al., 2017; Vanzella et al., 2019), present
the densest star forming environments in the high-redshift Universe (Vanzella
et al., 2023). The very recent James Webb Space Telescope (JWST) observations
of the Cosmic Gems arc (Adamo et al., 2024) have revealed five massive star
clusters at $z=10.2$ with masses of
$M_{\star}\sim{10^{6}}\>\mathrm{M_{\odot}}$ and lensing-corrected sizes of
$\sim 1$ pc, corresponding to high stellar surface densities of
$\Sigma_{\star}\sim{10^{5}}\>\mathrm{M_{\odot}}\>\mathrm{pc}^{-2}$. The galaxy
GN-z11 observed by the JWST at $z=10.6$ with super-solar N/O abundance at low
metallicity (Bunker et al., 2023) might be an indication of early stages of
globular cluster formation, or even a chemical fingerprint of supermassive
stars (Charbonnel et al., 2023; Marques-Chaves et al., 2024).

The evolution of massive star clusters (Wang et al., 2016) and IMBH formation
and growth within them (Rizzuto et al., 2023; Arca Sedda et al., 2023a) is
typically studied in isolated simulation setups. There are a number of
cosmological setups focusing on the formation of proto-globular clusters and
their IMBH growth (Sakurai et al., 2017; Ma et al., 2020). However, the
cosmological setups typically lack the necessary mass or spatial resolution to
follow the collisional stellar dynamics at the level of the isolated setups.

Observations of star cluster formation regions (Zhang et al., 2001; Lada &
Lada, 2003; Larsen, 2004; Bastian et al., 2005; Grasha et al., 2017; Menon et
al., 2021) suggest that massive star cluster formation resembles a complex
hierarchical assembly rather than a simple monolithic collapse, and gas and
gas expulsion likely has a role in shaping the young star clusters (Dale et
al., 2015b; Krause et al., 2016). Hydrodynamical simulations of star cluster
formation in low-metallicity star-burst environments support this view (Lahén
et al., 2020). Star cluster formation from fractal or otherwise structured
initial conditions has been studied in the literature (e.g. Aarseth & Hills
1972; Smith et al. 2011; Fujii et al. 2012; Howard et al. 2018; Grudić et al.
2018; Sills et al. 2018; Torniamenti et al. 2022; Farias et al. 2024), though
definitely not as widely as the evolution of simple, isolated spherical
systems, partially due to the fact that complex, structured initial conditions
evolve relatively rapidly into spherical configurations. Mass segregation and
core collapse may occur earlier in these setups, especially if the initial
sub-clusters are sub-virial (Allison et al., 2009; Yu et al., 2011) or already
mass segregated (McMillan et al., 2007; Moeckel & Bonnell, 2009). It has been
pointed out by Fujii & Portegies Zwart (2013) that stellar collisions may
occur earlier in ensemble star cluster formation in a filamentary setup as the
sub-clusters have shorter relaxation and mass segregation time-scales compared
to the final assembled star cluster.

Observationally, IMBHs still remain elusive (Greene et al., 2020).
Extrapolating the well-established the $M_{\bullet}$–$\sigma$ relation (e.g.
Ferrarese & Merritt 2000; Gebhardt et al. 2000) for massive galaxies into
systems with lower masses and velocity dispersions, it has been proposed that
dwarf galaxies, especially those with an active low-luminosity nucleus, might
host IMBHs (Greene & Ho, 2007; Dong et al., 2012; Chilingarian et al., 2018;
Greene et al., 2020; Mezcua et al., 2018; Reines, 2022). Several ultra- or
hyper-luminous X-ray emitters (ULXs, HLXs) in nearby galaxies have been
interpreted as candidates for accreting IMBHs with masses ranging from
$m_{\bullet}\sim{200}\>\mathrm{M_{\odot}}$ to
$m_{\bullet}\sim{10^{5}}\>\mathrm{M_{\odot}}$ (e.g. Colbert & Mushotzky 1999;
Kaaret et al. 2001; Matsumoto et al. 2001; Strohmayer & Mushotzky 2003;
Patruno et al. 2006; Farrell et al. 2009; Mezcua et al. 2013; Pasham et al.
2014; Mezcua et al. 2015; Mezcua 2017; Kim et al. 2020). ALMA observations of
gas dynamics in dwarf elliptical NGC 404 (Davis et al., 2020) have revealed a
BH with a mass of $m_{\bullet}\sim{5.5\times 10^{5}}\>\mathrm{M_{\odot}}$
consistent with both a high-mass IMBH or a low-mass MBH. In addition, there is
stellar-dynamical evidence of dark objects in the centers of some globular
clusters (GCs) and stripped dwarf galaxies in the IMBH mass range, consistent
with concentrations of stellar-mass black holes, or possibly IMBHs (e.g.
Gebhardt et al. 2002; Ibata et al. 2009; Noyola et al. 2010; van der Marel &
Anderson 2010; Lützgendorf et al. 2011; Jalali et al. 2012; Lanzoni et al.
2013; Kamann et al. 2016; Baumgardt 2017; Pechetti et al. 2022; Della Croce et
al. 2024). However, a number of studies do not find such stellar-dynamical
evidence of IMBHs (Baumgardt et al., 2003; Baumgardt et al., 2019; Mann et
al., 2019; Hénault-Brunet et al., 2020), and characteristic IMBH accretion
signatures have not been observed in the clusters in question (Tremou et al.,
2018). Still, the most definite evidence of IMBHs to date remains the
gravitational-wave (GW) transient GW190521 (Abbott et al., 2020), a merger
between two stellar-mass black holes with masses of
$m_{\bullet}\sim{66}\>\mathrm{M_{\odot}}$ and
$m_{\bullet}\sim{85}\>\mathrm{M_{\odot}}$ resulting in a formation of an IMBH
with $m_{\bullet}\sim{150}\>\mathrm{M_{\odot}}$. In addition to gravitational
waves, electromagnetic transients such as tidal disruption events may help to
uncover the elusive IMBH population (e.g. Lin et al. 2018; Fragione & Kocsis
2018; Fragione et al. 2018).

The fact that there are no universally accepted solid detections of IMBHs in
globular clusters is not necessarily a problem for the IMBH formation models.
After their formation, IMBHs are initially located at the centers of their
host clusters. In numerical simulations, BHs sink into the center of the star
cluster to form a subsystem of BHs which is an ideal place for both Newtonian
and relativistic BH-BH and IMBH-BH interactions (Breen & Heggie, 2013; Kremer
et al., 2020; González Prieto et al., 2022; Maliszewski et al., 2022). Even
Newtonian three-body interactions can be sufficiently strong to eject (IM)BHs
from their star clusters, and gravitational-wave recoil kicks after BH-BH
mergers with recoil velocities up to several times $1000$ km/s can easily
unbind the merger remnant from its host cluster. Merging binaries with
misaligned spins or highly spinning BH components are more likely to receive
large kicks (Poisson & Will, 2014). As the effective spin parameter of a
multi-generation BH merger product approaches $\sim 0.7$ of the maximal BH
spin, it is likely that the BH growing by mergers will other BHs will
eventually be ejected from its host star cluster, especially if the cluster
mass and thus escape velocity is low. It is thus expected that the BH-BH
channel can efficiently grow IMBHs in the environments with the highest escape
velocities, such as nuclear star clusters (e.g. Fragione & Silk 2020; Fragione
et al. 2022a, b). Even though IMBHs would be ejected from their host clusters
through relativistic recoil kicks, their gravitational-wave signal in the IMBH
mass range should be observable in the near future with either ground-based or
space-borne gravitational-wave detectors (e.g. Jani et al. 2020).

In this first study of our FROST-CLUSTERS project, we examine the formation of
IMBHs during the hierarchical assembly of massive star clusters. IMBH
formation in massive cluster star assembly setups (final cluster masses
$M_{\mathrm{\star}}>{5\times 10^{5}}\>\mathrm{M_{\odot}}$ and total particle
numbers $N_{\mathrm{tot}}=2.35\times 10^{6}$) with a full IMF from
${0.08}\>\mathrm{M_{\odot}}$ to ${150}\>\mathrm{M_{\odot}}$ including post-
Newtonian equations of motion for the BHs and gravitational-wave recoil kicks
for the BH-BH mergers. For simplicity, the first study focuses on single-star
setups. We will examine the effect of binary fraction on our results in a
forthcoming study (FROST-CLUSTERS II).

Our setup relies on two main building blocks. First, the individual sub-
clusters follow an universal young star cluster mass function with a power-law
slope of $\alpha=-2$. This is supported both by observations (Elmegreen &
Efremov, 1996; Zhang & Fall, 1999; Adamo et al., 2020) and numerical
simulations (Lahén et al., 2020). Second, we assume a shallow power-law mass-
size relation slope $\alpha\sim 0.18$ (e.g. Brown & Gnedin 2021) for our sub-
cluster population. The chosen normalization results in small sub-pc birth
radii for the cluster population, motivated by the recent JWST observations of
young, massive and dense star clusters at $z>10$ (Adamo et al., 2024) as well
as massive star cluster formation in state-of-the art star-by-star low-
metallicity hydrodynamical simulations of star-bursting dwarf galaxies (Lahén
et al., 2020).

IMBHs from $m_{\bullet}\sim{100}\>\mathrm{M_{\odot}}$ to
$m_{\bullet}\sim{2200}\>\mathrm{M_{\odot}}$ form during the first $t=50$ Myr
in our simulations. Most interestingly, we find that the initially most
massive star clusters do not form IMBHs from their stars, but rather inherit
IMBHs formed in lower-mass star clusters through star cluster mergers. We
study the reason for the suppressed IMBH formation in isolated, massive star
clusters, in a series of $100$ isolated setups containing up to million stars
per cluster.

The article is structured as follows. After the introduction we review basic
pathways to stellar collisions in dense star clusters in section 2. Next, we
introduce our updated simulation code BIFROST in section 3 as well as in the
Appendix, and the initial conditions for our simulations. We present the
results concerning the IMBH formation in the hierarchical star cluster
assembly simulations, detailing the complex growth histories of the IMBHs, in
section 4, and in isolated star clusters in section 5. After discussing the
implications and caveats of our results in section 6, we summarize our results
and conclude in section 7.

##  2 Stellar collisions in dense star clusters

###  2.1 Stellar collision rate: the cross section argument

Estimates for collision rates of stars can be made based on cross-section
arguments (Spitzer, 1987). Assuming an equilibrium stellar system with a
number density $n$ and velocity dispersion $\sigma$, consisting of stars with
masses $m_{\star}$ and radii $r_{\star}$, the collision rate
$t_{\mathrm{coll}}^{\mathrm{-1}}$ of stars can be estimated as

$\frac{1}{t_{\mathrm{coll}}}=4\sqrt{\pi}n\sigma
r_{\mathrm{\star}}^{2}\left(1+\frac{Gm_{\star}}{\sigma^{2}r_{\mathrm{\star}}}\right),$
(1)

following Binney & Tremaine (2008). If gravitational focusing is unimportant,
i.e. the second term in the parenthesis is small, the collision rate scales
with the velocity dispersion $\sigma$. In a typical star cluster, this is not
the case as

$\frac{Gm_{\star}}{\sigma^{2}r_{\mathrm{\star}}}\sim 1.9\times
10^{3}\left(\frac{m_{\star}}{1\>\mathrm{M_{\odot}}}\right)\left(\frac{r_{\star}}{1\>\mathrm{R_{\odot}}}\right)^{-1}\left(\frac{\sigma}{10\>\mathrm{km/s}}\right)^{-2}\gg
1.$ (2)

Instead, gravitational focusing is important and the collision rate is
inversely proportional to $\sigma$, and the collision rate decreases in high
velocity dispersion environments. The collision rates predicted by the simple
cross section argument of Eq. (1) are typically low. For example, a dense star
cluster following a Plummer (1911) density profile with a total mass of
$M_{\star}={2\times 10^{5}}\>\mathrm{M_{\odot}}$ and half-mass radius of
$r_{\mathrm{h}}=0.25$ pc consisting of equal-mass main-sequence stars of
$m={10}\>\mathrm{M_{\odot}}$ has a collision rate of only
$t_{\mathrm{coll}}^{\mathrm{-1}}\sim 10^{-3}$ Myr-1.

Numerical simulations of evolving star clusters have shown that the collision
rate can be orders of magnitude higher compared to the simple cross section
arguments. This is due to mass segregation and cluster core collapse which
increase the central stellar density $\rho_{\mathrm{c}}$. In addition, after
mass segregation and core collapse, massive stars can form binaries and more
complex hierarchies in the cluster core, enhancing the collision rate (e.g.
Portegies Zwart et al. 1999).

###  2.2 Mass segregation and core collapse

In star clusters consisting of equal-mass stars, core collapse cannot be
driven by mass segregation. Instead, the core collapse proceeds on a time-
scale $t_{\mathrm{cc}}$ by a factor of $15$–$20$ (Spitzer & Shull, 1975; Cohn,
1980; Takahashi, 1995; Makino, 1996) longer than the common two-body
relaxation time-scale $t_{\mathrm{rlx}}$ (Spitzer, 1987; Giersz & Heggie,
1994; Aarseth, 2003; Heggie & Hut, 2003). The expression for the two-body
relaxation time-scale can be written as

$t_{\mathrm{rlx}}=\frac{0.138N}{\ln{(\gamma
N)}}\left(\frac{r_{\mathrm{h}}^{3}}{GM_{\star}}\right)^{1/2}.$ (3)

Here $M_{\star}$ is the total stellar mass, $r_{\mathrm{h}}$ is the half-mass
radius of the cluster and $N=M_{\star}/\tilde{m}$ the total number of stars
with $\tilde{m}$ being the mean stellar mass. The parameter $\gamma$ has a
value of $\gamma\sim 0.11$ for single-mass systems (Giersz & Heggie, 1994).

In multi-mass systems, massive stars will reach the center of the cluster due
to energy equipartition on a time-scale of the order of the segregation time-
scale $t_{\mathrm{seg}}$ (Spitzer & Hart, 1971; Portegies Zwart et al., 2004),
defined as

$t_{\mathrm{seg}}=\frac{M_{\mathrm{\star}}}{m_{\mathrm{max}}}\frac{0.138}{\ln{(\gamma
M_{\mathrm{\star}}/m_{\mathrm{max}})}}\left(\frac{r_{\mathrm{h}}^{3}}{GM_{\mathrm{\star}}}\right)^{1/2},$
(4)

in which $m_{\mathrm{max}}$ is the mass of the most massive individual star in
the cluster. Comparing the expression for $t_{\mathrm{seg}}$ to the definition
of the two-body relaxation time-scale in Eq. 3, we see that the total number
of particles $N=M_{\star}/\tilde{m}$ is replaced by the effective number of
particles $M_{\star}/m_{\mathrm{max}}$ in the expression for
$t_{\mathrm{seg}}$.

For our dense example star cluster fro section 2.1 with $M_{\star}={2\times
10^{5}}\>\mathrm{M_{\odot}}$ and $r_{\mathrm{h}}=0.25$, this time populated
from a Kroupa (2001) IMF with $m_{\mathrm{max}}={150}\>\mathrm{M_{\odot}}$,
$t_{\mathrm{rlx}}\sim 18.5$ Myr and $t_{\mathrm{seg}}\sim 0.15$ Myr. For
realistic IMFs the core collapse time-scale accelerated by mass segregation is
$t_{\mathrm{cc}}\sim 0.2\>t_{\mathrm{rlx}}$ (Portegies Zwart et al., 2004), a
considerably shorter time-scale than for single-mass systems for which
$t_{\mathrm{cc}}\sim 15$–$20\>t_{\mathrm{rlx}}$. For our example cluster, this
yields $t_{\mathrm{cc}}\sim 3.7$ Myr.

Assuming a power-law mass-size relation $r_{\mathrm{h}}\propto
M_{\star}^{\alpha}$ for the star clusters, the central density
$\rho_{\mathrm{c}}$, the two-body relaxation $t_{\mathrm{rlx}}$ and mass
segregation times scale as

$\begin{split}\rho_{\mathrm{c}}&\propto M_{\star}^{\mathrm{1-3\alpha}}\\\
t_{\mathrm{rlx}}&\propto\frac{M_{\mathrm{\star}}^{\frac{1}{2}({3\alpha+1})}}{\ln{(\gamma
M_{\mathrm{\star}}/\tilde{m})}}\\\
t_{\mathrm{seg}}&\propto\frac{M_{\mathrm{\star}}^{\frac{1}{2}({3\alpha+1})}}{\ln{(\gamma
M_{\mathrm{\star}}/m_{\mathrm{max}})}}.\end{split}$ (5)

The central density formula assumes the common Plummer (1911) profile. Cluster
populations with $\alpha<1/3$ have their central density $\rho_{\mathrm{c}}$
increasing when the mass of the clusters $M_{\star}$ increases. For any mass-
size relation with a power-law slope larger than $\alpha>-1/3$, more massive
star clusters have longer two-body relaxation time-scales $t_{\mathrm{rlx}}$
as well as longer mass segregation time-scales $t_{\mathrm{seg}}$. The
observed mass-size relation slopes $\alpha$ for young star clusters in the
local Universe are in general positive but fairly shallow, with $\alpha\sim
0.13$ (Marks & Kroupa, 2012) for the three-dimensional $r_{\mathrm{h}}$ and
$\sim 0.1$–$0.3$ for the projected effective radii (Larsen, 2004; Brown &
Gnedin, 2021). Thus, more massive star clusters have higher central densities
$\rho_{\mathrm{c}}$ but longer mass segregation time-scales $t_{\mathrm{seg}}$
and thus longer core collapse time-scales $t_{\mathrm{cc}}$ than their less
massive counterparts in the local, young star cluster population.

More detailed core collapse time-scale estimates have been proposed in the
literature based on the sinking time-scale for massive stars due to dynamical
friction, the core relaxation time-scale, or the formation time-scale for hard
binaries which eventually halt the core collapse (e.g. Portegies Zwart &
McMillan 2002; Gürkan et al. 2004; Goswami et al. 2012; Fujii & Portegies
Zwart 2014). The mass segregation time-scale $t_{\mathrm{seg}}$ is typically
short, less than $1$-$2$ Myr for dense ($r_{\mathrm{h}}\lesssim 1$ pc) star
clusters. However, mass segregation and core collapse can only enhance the
collision rate of massive stars if the stars are still alive when reaching the
core (Portegies Zwart et al., 2004). The short life-times of massive stars, of
the order of a few Myr, limit the collision rate: for example, a star with an
initial mass of $m_{\star}={150}\>\mathrm{M_{\odot}}$ at a low metallicity of
$Z=0.01Z_{\odot}$ has a life-time of only $t=2.82$ Myr. 333The estimates for
stellar life-times in this work originate from the PARSEC stellar tracks (e.g.
Bressan et al. 2012) used by the fast stellar population synthesis code SEVN
(Iorio et al., 2023).

###  2.3 Binary formation, binary-single interactions and collisions

If there are no primordial binaries in the star cluster, first binaries form
during the core collapse from encounters of three initially unbound stars via
the three-body binary formation channel (hereafter 3BBF). Assuming a system
consisting of identical stars of mass $m_{\star}$ with a number density of $n$
and velocity dispersion of $\sigma$, it can be shown (Goodman & Hut, 1993;
Binney & Tremaine, 2008) that the rate $\bar{C}$444Note that textbook of
Binney & Tremaine (2008) erroneously has $n^{2}$, not $n^{3}$, in their
formula 7.11 discussing the Goodman & Hut (1993) calculation. in which
binaries form in triple encounters is

$\bar{C}\approx\frac{G^{5}m_{\star}^{5}n^{3}}{\sigma^{9}}.$ (6)

Thus, for a system consisting of equal-mass bodies, the 3BBF binary formation
rate $\bar{C}$ steeply increases with increasing stellar mass $m_{\star}$ and
number density $n$, and strongly decreases with increasing velocity dispersion
$\sigma$.

However, the Goodman & Hut (1993) binary formation rate estimate of Eq. (6)
assumes equal-mass stars, and it is unclear to what extent their 3BBF rates
are applicable to realistic multi-mass encounters. Nevertheless, it has been
commonly assumed that two most massive bodies are the most likely to pair.
Very recently, Atallah et al. (2024) showed that this assumption does not
actually hold. Instead, from the three bodies ($m_{\mathrm{1}}$,
$m_{\mathrm{2}}$, $m_{\mathrm{3}}$), the two most massive ($m_{\mathrm{1}}$,
$m_{\mathrm{2}}$) are the least likely to become bound. Wide binaries555Wide
and hard binaries are characterized by their long-term survival in repeated
encounters with the cluster stars. Wide binaries will eventually dissolve
which hard binaries can further shrink (Heggie, 1975). The threshold is
approximately the binary orbital velocity $v_{\mathrm{orb}}\sim\sigma$. For
example, a binary with $m_{\mathrm{b}}=20\>\mathrm{M_{\odot}}$ and $a=10^{-3}$
pc is hard in star clusters with $\sigma\lesssim 9.3$ km/s., which are a
considerably more likely outcome than hard binaries, and stellar collisions
are more common than hard binary formation. The binary eccentricity
distribution is always thermal or steeper, i.e. $f(e)\propto e^{\beta}$ with
$\beta\gtrsim 1$.

Binaries formed through the 3BBF channel subsequently interact with single
stars and other binaries in the star cluster via fly-bys and exchanges, in
which the incoming third body replaces a binary member. The probability for
exchanges strongly increases with the increasing mass of the third body
(Valtonen & Karttunen, 2006), so massive stars prominently end up in binary
systems through the exchange channel. In fly-bys, binaries just exchange
energy and angular momentum with the single cluster members. On average, hard
binaries will further harden, and wide binaries will widen in dynamical
interactions. This fundamental result is know as Heggie’s law or the Heggie-
Hills law (Heggie, 1975; Hills, 1975).

The increased interaction cross-section of the binary stars compared to single
stars enhances the stellar collision rate in the star clusters. It has been
demonstrated that the runaway collision sequence in dense star clusters
typically begins with such a single-binary interaction including massive stars
(Gaburov et al., 2008).

###  2.4 Can the clusters retain their (IM)BHs?

Figure 1: The gravitational wave recoil kick velocities $v_{\mathrm{recoil}}$
(Zlochower & Lousto, 2015) for merging BH binaries as a function of the
primary mass $m_{\mathrm{p}}$ for three secondary masses
$m_{\mathrm{s}}={10}\>\mathrm{M_{\odot}}$ (in blue),
$m_{\mathrm{s}}={60}\>\mathrm{M_{\odot}}$ (in orange) and
$m_{\mathrm{s}}={290}\>\mathrm{M_{\odot}}$ (in yellow). The primary mass range
is $m_{\mathrm{s}}\leq m_{\mathrm{p}}\leq{10^{4}}\>\mathrm{M_{\odot}}$. We
show the two extreme cases of aligned (solid lines) and anti-aligned spins
(dashed lines). The spin magnitudes originate from the Geneva model
(Belczynski et al., 2020). Escape velocities for star clusters with masses up
to $M_{\star}={10^{6}}\>\mathrm{M_{\odot}}$ are shown in horizontal gray
lines. The cluster densities are similar as in numerical simulations in this
work. In general, only equal-mass mergers or mergers with considerable mass
ratios ($m_{\mathrm{s}}/m_{\mathrm{p}}\lesssim 0.04$–$0.1$) lead to weak kicks
after which the remnant BH can remain in the cluster. Increasingly misaligned
spins cause stronger recoil kicks.

Even if IMBHs would form in young, massive star clusters, they might not be
indefinitely retained in their host clusters. Assuming the common Plummer
(1911) model, the escape velocity $v_{\mathrm{esc}}$ from the center of the
cluster is

$v_{\mathrm{esc}}\approx
33.4\>\mathrm{km/s}\left(\frac{M_{\mathrm{\star}}}{10^{5}\>\mathrm{M_{\odot}}}\right)^{1/2}\left(\frac{r_{\mathrm{h}}}{1.0\>\mathrm{pc}}\right)^{-1/2},$
(7)

only of the order of a few tens of km/s to $\sim 100$ km/s for star clusters
in the mass range of
$M_{\star}={10^{5}}\>\mathrm{M_{\odot}}$–${10^{6}}\>\mathrm{M_{\odot}}$. Even
Newtonian three-body and multiple interactions are occasionally strong enough
to unbind BHs and stars from their host clusters, especially in low escape
velocity environments of lower-mass star clusters and in clusters containing
multiple IMBHs (Maliszewski et al. 2022, Souvaitzis et al., in prep.).

Relativistic gravitational-wave recoil kicks (Bekenstein, 1973; Campanelli et
al., 2007b, a; Zlochower et al., 2011) caused by anisotropic emission of
gravitational radiation during BH-BH mergers can easily eject the BH merger
products from their stellar environments. The recoil kick velocity is low for
merging equal-mass BH binaries ($m_{\mathrm{p}}=m_{\mathrm{s}}$) and binaries
with a considerable mass ratio, i.e. $m_{\mathrm{s}}/m_{\mathrm{p}}\ll 1$.
Spinning BHs receive in general stronger recoil kicks compared to non-spinning
BHs, especially if the BH spin directions are anti-aligned. The maximum recoil
kick velocity is of the order of $v_{\mathrm{recoil}}\sim 4000$ km/s for
spinning BHs (Campanelli et al., 2007a), and $v_{\mathrm{recoil}}\sim 175$
km/s for non-spinning ones (González et al., 2007).

We illustrate typical gravitational-wave recoil kick velocities compared to
escape velocities from dense star clusters in Fig. 1. The recoil velocities
are calculated using the fitting formulas by Zlochower & Lousto (2015). We use
two spin direction options (aligned and anti-aligned spins), and three
different secondary BH masses of $m_{\bullet}={10}\>\mathrm{M_{\odot}}$
(common low-mass BHs), $m_{\bullet}={60}\>\mathrm{M_{\odot}}$ (BHs just below
the (P)PISN mass gap) and $m_{\bullet}={290}\>\mathrm{M_{\odot}}$ (IMBHs just
above the (P)PISN mass gap). The primary BH mass $m_{\mathrm{p}}$ ranges in
each model from $m_{\mathrm{s}}$ up to ${10^{4}}\>\mathrm{M_{\odot}}$ in the
IMBH mass range. We assume the Geneva model for BH spins (Belczynski et al.,
2020) as in their Eq. (3) with $Z=0.0004=0.02\>Z_{\odot}$.

From Fig. 1 it is evident that it is difficult to retain a merger remnant
originating from a BH-BH binary that includes a BH below the mass gap
($m_{\mathrm{p}}\lesssim{63}\>\mathrm{M_{\odot}}$), unless the merging binary
is equal-mass or nearly so, and has almost aligned spins. Primary IMBHs above
the mass gap ($m_{\mathrm{p}}\gtrsim{290}\>\mathrm{M_{\odot}}$) can merge with
low-mass BHs ($m_{\mathrm{s}}={10}\>\mathrm{M_{\odot}}$) with any relative
spin orientation and be retained already in star clusters with masses
$M_{\star}\sim{10^{3}}\>\mathrm{M_{\odot}}$–${10^{4}}\>\mathrm{M_{\odot}}$.
However, the primary IMBH masses above
$m_{\bullet}\gtrsim{500}\>\mathrm{M_{\odot}}$ and
$m_{\bullet}\gtrsim{1000}\>\mathrm{M_{\odot}}$ are required to retain the
merger product in the star clusters with masses
$M_{\star}={10^{6}}\>\mathrm{M_{\odot}}$ and
$M_{\star}={10^{5}}\>\mathrm{M_{\odot}}$, respectively, after a merger with
any stellar-mass BH below the mass gap. IMBH primaries with
$m_{\mathrm{p}}\gtrsim{2000}\>\mathrm{M_{\odot}}$ are difficult to remove from
any star cluster in mergers with stellar-mass BHs. Finally, IMBH-IMBH mergers
may remove even massive primary IMBHs from their host clusters, provided that
the IMBH formation channel results in multiple IMBHs in a single star cluster.
A primary mass of $m_{\mathrm{p}}\sim{4500}\>\mathrm{M_{\odot}}$ is required
to retain the remnant after a merger with an above-mass-gap secondary IMBH
with $m_{\bullet}{290}\>\mathrm{M_{\odot}}$ in a star cluster with
$M_{\star}={10^{5}}\>\mathrm{M_{\odot}}$.

##  3 Numerical methods and initial conditions

###  3.1 Numerical methods: the updated BIFROST code

We use an updated version of the BIFROST code (Rantala et al., 2021; Rantala
et al., 2023) for the simulations of this study. In summary, BIFROST is a
novel GPU-accelerated collisional direct-summation N-body code based on the
hierarchical fourth-order forward symplectic integrator. Few-body subsystems
such as binary stars, close fly-bys, hierarchical multiples and small clusters
around black holes are treated using both secular and regularized integration
techniques, including post-Newtonian equations of motion for BHs up to order
PN1.0 in secular techniques and PN3.5 in regularized integration. BIFROST has
already been used to study a number of stellar-dynamical topics involving BHs,
IMBHs and MBHs, such as IMBH growth through repeated TDEs (Rizzuto et al.,
2023), the dynamics of parsec-scale stellar disks in galactic nuclei and
milliparsec-scale S-cluster analogues around MBHs (Rantala & Naab, 2024), and
escapers from merging older star clusters with and without pre-existing IMBHs
(Souvaitzis et al., in prep.).

The main differences of the updated BIFROST code version used in this study
compared to the original code of Rantala et al. (2023) are the following.
First, a key element of the fourth-order forward integration, the calculation
of the term involving the acceleration gradient, is simplified (appendix A.1).
Next, the secular triple integration and the standard algorithmically
regularized few-body integration (LogH) are replaced by an implementation of
the slow-down algorithmic regularization (SDAR) integrator (appendix A.2).
Finally, we present the coupling of the fast stellar population synthesis code
SEVN (Iorio et al., 2023) into BIFROST allowing stellar evolution in the
N-body simulations (appendix B.2). In this study, we focus on single stellar
evolution in the initial stellar mass range of
${2.2}\>\mathrm{M_{\odot}}\lesssim
m_{\star}\lesssim{600}\>\mathrm{M_{\odot}}$. Each of the three code updates is
described in detail in the appendix.

###  3.2 Initial conditions: individual star cluster models

The individual star clusters in our initial models follow the Plummer (1911)
density profile

$\rho(r)=\frac{3M_{\mathrm{\star}}}{4\pi
a^{3}}\left(1+\frac{r^{2}}{a^{2}}\right)^{-5/2}$ (8)

and the corresponding potential. Here $M_{\mathrm{\star}}$ is the total
stellar mass of the model and $a$ is the scale radius related to the three-
dimensional half-mass radius $r_{\mathrm{h}}$ of the model as
$r_{\mathrm{h}}=(2^{2/3}-1)^{-1/2}a\approx 1.3a$. The projected effective
radius $R_{\mathrm{e}}$ is simply $R_{\mathrm{e}}=a$. The spherically
symmetric positions and isotropic velocities of the individual stars are
sampled from the distribution function of the Plummer density-potential pair.

Simulation studies of massive star cluster evolution or collisional IMBH
formation do not typically assume an initial mass-size relation for their sets
of star cluster models, but rather select $M_{\star}$ (or $N$), and the
central $\rho_{\mathrm{c}}$ or the half-mass density $\rho_{\mathrm{h}}$ to
examine the cluster evolution or IMBH formation at different chosen cluster
masses and densities. Our approach in this study is to choose an
observationally motivated mass-size relation which all of our star cluster
models will initially follow.

The observed effective radii $R_{\mathrm{e}}$ of star clusters follow a power-
law star cluster mass-size relation

$\frac{R_{\mathrm{e}}}{\mathrm{pc}}=R_{\mathrm{4}}\left(\frac{M_{\mathrm{\star}}}{10^{4}M_{\odot}}\right)^{\beta}$
(9)

for which Brown & Gnedin (2021) provide values of $\beta=0.180\pm 0.028$ and
$R_{\mathrm{4}}=2.365\pm 0.106$ for star cluster of ages between $1$ Myr and
$10$ Myr in the LEGUS sample (Adamo et al., 2017; Cook et al., 2019). However,
it has been suggested that the birth radii of star clusters may be smaller,
and clusters increase in size during later evolution. Studying young clusters
in their embedded phase, Marks & Kroupa (2012) find that the half-mass radius
scales as the function of the embedded cluster mass $M_{\mathrm{ecl}}$ as

$\frac{r_{\mathrm{h}}}{\mathrm{pc}}=0.10_{\mathrm{-0.04}}^{\mathrm{+0.07}}\times\left(\frac{M_{\mathrm{ecl}}}{M_{\odot}}\right)^{\mathrm{0.13\pm
0.04}}.$ (10)

Thus, a size increase by a factor of up to $\sim 10$ is required from the
embedded phase to the young massive cluster phase. The sizes of star clusters
increase due to internal collisional dynamics, external tides and stellar
evolution (e.g. Arca Sedda et al. 2023a), with gas expulsion suggested as an
additional physical process contributing to the size growth (Banerjee &
Kroupa, 2017).

We adopt small initial birth radii for our star cluster model population. The
relation between the stellar mass $M_{\star}$ to the half-mass radii
$r_{\mathrm{h}}$ of our clusters is

$\frac{r_{\mathrm{h}}}{\mathrm{pc}}=\frac{f_{\mathrm{h}}R_{\mathrm{4}}}{1.3}\left(\frac{M_{\mathrm{\star}}}{10^{4}M_{\odot}}\right)^{\beta}$
(11)

in which $f_{\mathrm{h}}$ is a free parameter and the constant
$r_{\mathrm{h}}/R_{\mathrm{e}}\sim 1.3$ arises from the properties of the
Plummer model. While using the values of $\beta$ and $R_{\mathrm{4}}$ and
their uncertainties from Brown & Gnedin (2021) in Eq. (9), we set
$f_{\mathrm{c}}=1/8$ so that at $M_{\star}={10^{5}}\>\mathrm{M_{\odot}}$ the
half-mass radius of the star cluster ($r_{\mathrm{h}}=0.45$ pc) corresponds to
the size from the relation of Marks & Kroupa (2012) in Eq. (10).

We illustrate the particle numbers and half-mass densities of a number of
direct N-body (Wang et al., 2015; Di Carlo et al., 2021; Banerjee, 2021;
Rastello et al., 2021; Rizzuto et al., 2021; Kamlah et al., 2022; Arca Sedda
et al., 2023a) and cluster Monte Carlo (Askar et al., 2017; Rodriguez et al.,
2016; Rodriguez et al., 2019; Kremer et al., 2020; Maliszewski et al., 2022)
simulations following Fig. 1 of Arca Sedda et al. (2023a) in Fig. 2. We show
both our isolated star cluster models, and the population of star clusters for
the hierarchical cluster assembly simulations. The cluster assembly region
simulations consist of individual star clusters with particle numbers in the
range $1.7\times 10^{2}\lesssim N\lesssim 4.0\times 10^{5}$ particles, up to
$N_{\mathrm{tot}}\sim 2.35\times 10^{6}$ particles per simulation. The
isolated simulations we run have up to $N=1.0\times 10^{6}$ particles. Our
adopted mass-size relation in Eq. (11) results in an initial star cluster
population consistent with other recent simulation studies of massive star
cluster evolution and IMBH formation.

Figure 2: A comparison of the initial conditions of recent direct N-body (blue
filled circles) and Monte Carlo (red filled squares) studies of star clusters
following Fig. 1 of Arca Sedda et al. (2023a). The figure shows the half-mass
density $\rho_{\mathrm{h}}$ of the star clusters as a function of their
particle number $N$. The references to the individual studies are listed in
the main text. Our isolated star clusters (up to $N=10^{6}$) are indicated in
purple cross symbols while the initial star clusters in the hierarchical
cluster assembly simulations (in total $N_{\mathrm{tot}}=2.35\times 10^{6}$
particles) are shown in green. The densities of our star cluster models are
consistent with the other recent star cluster studies in the literature.

We sample the masses of individual stars from the standard piece-wise power-
law stellar initial mass function of Kroupa (2001) as

$\frac{\mathrm{d}N_{\mathrm{\star}}}{\mathrm{d}m_{\mathrm{\star}}}\propto
m_{\star}^{-\alpha}$ (12)

with $\alpha=1.3$ when $0.08M_{\odot}\leq m_{\star}\leq 0.50M_{\odot}$ and
$\alpha=2.3$ for $0.50M_{\odot}<m_{\star}\leq m_{\mathrm{max}}$. The maximum
mass of an individual star $m_{\mathrm{max}}$ in a single star cluster depends
on the mass $M_{\mathrm{\star}}$ of the cluster in question (Weidner & Kroupa,
2006). We adopt a simple piece-wise power-law model

$\log_{\mathrm{10}}\left(\frac{m_{\mathrm{max}}}{M_{\odot}}\right)=a\times\log_{\mathrm{10}}\left(\frac{M_{\mathrm{\star}}}{M_{\odot}}\right)+b.$
(13)

We choose the constant coefficients $a$ and $b$ presented in Table 1 in such a
manner that the simple model reproduces the main features of the cluster mass
– maximum stellar mass relation of Yan et al. (2023). The maximum individual
stellar masses for star clusters with masses
$M_{\mathrm{\star}}=10^{2}M_{\odot}$, $M_{\mathrm{\star}}=10^{3}M_{\odot}$ and
$M_{\mathrm{\star}}=10^{4}M_{\odot}$ are approximately $m_{\mathrm{max}}\sim
11M_{\odot}$, $m_{\mathrm{max}}\sim 25M_{\odot}$ and $m_{\mathrm{max}}\sim
85M_{\odot}$, respectively. Furthermore, we adopt an overall maximum initial
stellar mass of $m_{\mathrm{max}}=150M_{\odot}$ for the most massive star
clusters. The maximum mass of a single star in a star cluster is not typically
strictly enforced when generating star cluster initial conditions for
simulations. Instead, massive stars are simply less probable to occur in low-
mass clusters. In our case, our hierarchical star cluster assembly initial
setup contains a large number of low-mass clusters, and we wish to avoid small
star clusters with very massive individual stars.

Cluster mass range | coeff. a | coeff. b
---|---|---
$\log_{\mathrm{10}}(M_{\mathrm{\star}}/M_{\odot})\leq 1.5$ | $+0.51$ | $-0.17$
$1.5<\log_{\mathrm{10}}(M_{\mathrm{\star}}/M_{\odot})\leq 2.0$ | $+0.95$ | $-0.84$
$2.0<\log_{\mathrm{10}}(M_{\mathrm{\star}}/M_{\odot})\leq 2.6$ | $+0.26$ | $+0.52$
$2.6<\log_{\mathrm{10}}(M_{\mathrm{\star}}/M_{\odot})$ | $+0.53$ | $-0.23$
Table 1: The constant coefficients $a$ and $b$ for the simple model for the
maximum mass of an individual star in a star cluster with a mass of
$M_{\mathrm{\star}}$ for Eq. (13).

We set the stellar metallicity to a constant value of
$Z=0.0002=0.01\;Z_{\odot}$ to emulate the conditions of the low-metallicity
star cluster environment at high redshifts. The initial stellar radii are
obtained by using the stellar evolutionary tracks of SEVN. Each star begins
the simulation at the beginning of its SEVN main sequence stellar evolutionary
phase, i.e. initially phase_star=1. All stars in the initial conditions are
single stars. We will explore the effect of primordial binaries on our setup
in a forthcoming study.

###  3.3 Initial conditions: star cluster assembly regions

The masses $M_{\star}$ of the star clusters are obtained from the power-law
cluster mass function

$\frac{\mathrm{d}N_{\mathrm{cl}}}{\mathrm{d}M_{\mathrm{\star}}}\propto
M_{\mathrm{\star}}^{\alpha}$ (14)

for which both observations (e.g. Elmegreen & Efremov 1996; Zhang & Fall 1999;
Adamo et al. 2020) and solar-mass resolution hydrodynamical simulations (Lahén
et al., 2020) suggest a slope of $\alpha\sim-2$. We sample star cluster masses
from the cluster mass function until the desired total stellar particle number
$N_{\mathrm{tot}}$ is reached. For the hierarchical cluster assembly
simulations of this study $N_{\mathrm{tot}}$ is in the range $1.7\times
10^{6}\leq N_{\mathrm{tot}}\leq 2.35\times 10^{6}$. We show the particle
numbers $N$ and the half-mass densities $\rho_{\mathrm{h}}$ of the individual
star clusters in Fig. 2.

After the star cluster masses $M_{\mathrm{\star}}$ have been obtained, we set
the initial center-of-mass positions and velocities of the clusters. Our setup
is motivated by the structure of star cluster formation regions in the solar-
mass resolution hydrodynamical simulations of star-bursting dwarf galaxies of
Lahén et al. (2020). We choose a representative star formation region with
dimensions of $47$ pc, $77$ pc and $84$ pc from the simulation volume which
contains a large number of young ($t\leq 2$ Myr) star clusters (for an
illustration, see the top panel in Fig. 2 of Lahén et al. 2019). The total
stellar mass in young stars within the region is $M_{\star}\sim 4.8\times
10^{5}M_{\odot}$ while most massive individual gravitationally bound cluster
has a mass of $M_{\mathrm{\star}}\sim 1.5\times 10^{5}M_{\odot}$. Centering on
the most massive cluster, the radial cumulative mass function is close to a
power-law distribution $M_{\mathrm{tot}}(<r)\propto r^{\gamma}$ with
$\gamma\sim 3$, indicating a nearly constant density of the star cluster
region. The region is in a state of collapse with a mean radial infall
velocity of $\langle v_{\mathrm{r}}\rangle\approx-7$ km/s. In addition to the
collapse, the star clusters have random motions up to
$\|\boldsymbol{v}-\langle\boldsymbol{v}_{\mathrm{r}}\rangle\|\sim 14$ km/s.
Not all clusters are gravitationally bound to the bulk of the system. Half of
the stellar mass of the final proto-globular cluster is still yet to form and
the sub-clusters are embedded in a filamentary gaseous structure.

We begin the initial setup construction by placing our most massive sampled
Plummer model at the origin at zero velocity. For simplicity, we choose a
spherically symmetric distribution of cluster positions even though the
spatial distribution of the representative star cluster region in the
hydrodynamical simulation of Lahén et al. (2020) is a somewhat flattened
configuration. The assumption allows us to sample the star cluster positions
within a homogeneous sphere. For the cluster maximum separation from the
origin we set $r_{\mathrm{max}}=50$ pc. For the star cluster velocities we use
an isotropic random component up to $3.5$ km/s combined with a radial
component of $v_{\mathrm{r}}=-3.5$ km/s for each cluster. The radial cluster
infall velocities are by a factor of $\sim 2$ smaller than in the
representative region of Lahén et al. (2020) to compensate for the fact that
half of the stellar mass of the region in the hydrodynamical simulation is
still to be formed. The chosen random velocity component ensures that most of
the stellar mass in the region will end up in a single cluster or its diffuse
envelope.

Simulation | $N_{\mathrm{cl}}$ | $N_{\mathrm{tot}}$ | $N_{\mathrm{min}}$ | $N_{\mathrm{max}}$ | $M_{\mathrm{tot}}$ | $M_{\mathrm{min}}$ | $M_{\mathrm{max}}$
---|---|---|---|---|---|---|---
|  |  |  |  | $[M_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$
H1 | $1.10\times 10^{3}$ | $1.70\times 10^{6}$ | $170$ | $3.45\times 10^{5}$ | $9.42\times 10^{5}$ | $63.1$ | $2.02\times 10^{5}$
H2, H3 | $1.43\times 10^{3}$ | $2.35\times 10^{6}$ | $204$ | $4.07\times 10^{5}$ | $1.29\times 10^{6}$ | $73.2$ | $2.30\times 10^{5}$
Table 2: The three hierarchical star cluster assembly simulation setups H1, H2 and H3. The table lists the total number of star clusters ($N_{\mathrm{cl}}$), the total number of stars ($N_{\mathrm{tot}}$), the minimum and the maximum numbers of particles in a single cluster ($N_{\mathrm{min}}$ and $N_{\mathrm{max}}$), the total stellar mass ($M_{\mathrm{tot}}$), the minimum and the maximum masses of a single cluster ($M_{\mathrm{min}}$ and $M_{\mathrm{max}}$). The setups H2 and H3 are different random realizations from the same list of cluster masses with different particle masses, positions and velocities within the clusters as well as different initial cluster positions and velocities. BIFROST parameter | symbol | value
---|---|---
integration interval / | $\epsilon_{\mathrm{max}}$ |  $5\times 10^{-3}$ Myr
maximum time-step |  |
forward integrator time-step factor |  $\eta_{\mathrm{ff}}$, $\eta_{\mathrm{fb}}$, $\eta_{\mathrm{\nabla}}$ | $0.2$
subsystem neighbor radius | $r_{\mathrm{rgb,\star}}$ |  $0.25$ mpc
SDAR GBS tolerance | $\eta_{\mathrm{GBS}}$ | $5\times 10^{-8}$
SDAR GBS end-time tolerance | $\eta_{\mathrm{endtime}}$ | $10^{-2}$
regularization highest PN order |  | PN3.5
Table 3: The user-given BIFROST parameters used in the simulations of this
study.

We construct in total three hierarchical cluster assembly region setups which
we label H1, H2 and H3. The properties of the cluster assembly setups are
listed in Table 2. The setups H2 and H3 have the same cluster particle numbers
and cluster masses, but are otherwise different random realizations with
differing initial cluster positions and velocities as well as different
randomly drawn particle masses, positions and velocities for stars within the
clusters. Each cluster formation region is evolved using our BIFROST code
until $t=50$ Myr with user-given accuracy parameters listed in Table 3. The
simulations were run between November 2023 and March 2024 using the FREYA
cluster and the supercomputer RAVEN hosted by the Max Planck Computing and
Data Facility (MPCDF), in Garching, Germany using $2$–$4$ computing nodes with
either two Nvidia Tesla P100-PCIE-16GB or four Nvidia Tesla A100-PCIE-40GB
GPUs and $40$-$72$ CPUs per node.

##  4 IMBHs from star cluster assembly

Cluster | IMBH | $M_{\mathrm{\star}}$ | $t_{\mathrm{acc}}$ | $M_{\mathrm{\bullet}}(t_{\mathrm{acc}})$
---|---|---|---|---
|  | $\mathrm{[M_{\odot}]}$ | $\mathrm{[Myr]}$ | $\mathrm{[M_{\odot}]}$
H1-A |  | $2.0\times 10^{5}$ |  |
H1-B | BH1-B | $4.4\times 10^{4}$ | 6.5 | $1194$
H1-C | BH1-C | $1.1\times 10^{5}$ | 3.1 | $536$
H1-D | BH1-D | $5.9\times 10^{4}$ | 2.3 | $348$
H2-A |  | $2.3\times 10^{5}$ |  |
H2-B | BH2-B | $3.2\times 10^{4}$ | 4.8 | $1256$
H2-C | BH2-C | $2.6\times 10^{4}$ | - | $350$
H3-A |  | $2.3\times 10^{5}$ |  |
H3-B | BH3-B | $1.1\times 10^{5}$ | 2.2 | $1748$
H3-C | BH3-C | $1.5\times 10^{4}$ | 8.7 | $347$
H3-D | BH3-D | $2.3\times 10^{4}$ | 3.7 | $338$
H3-E | BH3-E | $1.4\times 10^{4}$ | 6.9 | $289$
H3-F | BH3-F | $2.6\times 10^{4}$ | 7.6 | $105$
Table 4: The initial stellar masses $M_{\mathrm{\star}}$ of the star clusters
which produce an IMBH during the cluster assembly simulations. The most
massive, central initial cluster is indicated with a label A. The accretion
time $t_{\mathrm{acc}}$ indicates when the star cluster merges with the
central cluster, and $M_{\mathrm{\bullet}}(t_{\mathrm{acc}})$ describes the
mass of the most massive IMBH at this time. If $t_{\mathrm{acc}}$ is small
enough, IMBHs have had no time to form yet. In this case the mass of the IMBH
descending from the most massive star at the time of the accretion is given.
Figure 3: The evolution of the star cluster assembly regions into a single,
massive star cluster in the simulation H1. The gray and white colors indicate
the projected stellar density of the clusters. Initially, the star clusters
populate a region with a radius of $50$ pc (top left panel). Most of the in-
falling star clusters interact and merge with the central star cluster between
$\sim 2$ Myr and $\sim 10$ Myr (top right and bottom left panels).
Intermediate-mass black holes (indicated with red arrows) from between $\sim
2$ Myr and $\sim 4$ Myr as collapsed remnants of VMSs formed in sequences of
repeated collisions of massive stars. VMSs in the (P)PISN mass gap leaving
either a small BH or no remnant are indicated in the top-right panel as small,
blue dots in dense star clusters. After $\gtrsim 10$ Myr the region settles
into a single, smooth, spherical star cluster with the remaining sub-clusters
either escaping the region or being disrupted by the main cluster. The bottom
right panel shows the cluster at $18.5$ Myr when an IMBH with
$m_{\bullet}=1498\;\mathrm{M_{\odot}}$ escapes the cluster due to a
gravitational-wave recoil kick of $v_{\mathrm{recoil}}=206$ km/s. Figure 4: A
simulated BVR filter view of the massive stars
($m_{\star}>2\;\mathrm{M_{\odot}}$) in the star cluster H1-A during its
evolution generated using the MYOSOTIS visualization code (Khorrami et al.,
2019). The top left panel displays the initial cluster assembly region as in
Fig. 3 while the remaining three panels focus on the central cluster at $t=5$
Myr (top right), $t=18.5$ Myr (bottom left) and $t=50.0$ Myr (bottom right).
At $t=5.0$ Myr, the cluster assembly is ongoing and substructure such as the
disrupting remains of the star cluster H1-B is clearly visible in the image
near and around the central cluster. At later times the stellar density of the
main cluster decreases due to the combined effects of collisional
gravitational dynamics and stellar evolution. In the simulation H1, no IMBHs
inhabit the star cluster at the end of the simulation due to gravitational-
wave recoil kicks exceeding the cluster escape velocity in the aftermath of
BH-IMBH and IMBH-IMBH mergers.

###  4.1 Star cluster assembly

The collapse of a star cluster assembly region H1 into a single, massive star
cluster is illustrated in Fig. 3. Initially, the in-falling star cluster
population resides in a volume with $r_{\mathrm{max}}=50$ pc in radius with
the most massive star cluster at the center. The massive star clusters begin
to strongly interact and merge with the central cluster after $t\gtrsim 2$
Myr. Before merging, clusters are tidally stripped in close encounters with
the main cluster, in some cases multiple times. Most, though not all, in-
falling star clusters eventually merge into the main cluster with most star
cluster mergers occurring before $t\lesssim 10$ Myr. As massive stars have
short life-times ($t=2.82$ Myr for $m_{\star}={150}\>\mathrm{M_{\odot}}$ at
$Z=0.01Z_{\odot}$), massive stars may already have reached the end of their
lifetimes when their host cluster merges into the main star cluster.

The initial substructure around the main cluster is gradually erased as
unbound material escapes from the vicinity of the cluster while the remaining
bound remnants of the in-falling clusters are tidally disrupted. After
$t\gtrsim 10$ Myr, the cluster assembly region has settled into a single
massive, spherical and smooth star cluster. We use a modified version of the
publicly available MYOSOTIS code (Khorrami et al., 2019) to generate mock BVR
images of the simulated star clusters. The evolution of the star cluster
region H1 and its initially most massive cluster H1-A is illustrated in mock
observations in Fig. 4. The visualizations only involve stars with masses
above $m_{\star}\gtrsim 2.2\>\mathrm{M_{\odot}}$ as this is the minimum
stellar mass in stellar tracks used in this study. IMBH positions and masses
are also indicated in the illustration.

Figure 5: The density profiles of the initial individual star cluster models
(gray lines) as well as the stars (blue lines) and the compact remnants
(orange lines) in the central cluster from $t=10$ Myr to $t=50$ Myr in
intervals of $10$ Myr. At $t=10$ Myr, the density of the central $\lesssim
0.10$ pc is dominated by compact remnants, mostly BHs, while at $t=50$ Myr the
BHs dominate the region within the central $0.22$ pc. Both the central stellar
and compact remnant densities decrease as a function of time while the outer
parts beyond $1.0$ pc (stars) and $0.3$ pc (remnants) remain relatively
unchanged. While the stellar profiles are almost flat at the center, the
compact remnants have a more cuspy central profile.

We present the density profiles of stars and compact remnants (mostly stellar-
mass BHs) in the central cluster in Fig. 5 between $t=10$ Myr and $t=50$ Myr
in intervals of $10$ Myr. IMBHs are excluded from the density profile
calculation. Already at $t=10$ Myr the profiles are smooth with little
indication of remaining substructure. The central stellar profile shows an
almost flat core as in the initial Plummer profile while the BHs have a cuspy
central structure with a power-law slope close to $-1$. The central stellar
density at $t=10$ Myr is $\rho_{\mathrm{c,\star}}\sim{7\times
10^{5}}\>\mathrm{M_{\odot}\>pc^{-3}}$ while the BH density is somewhat higher,
$\rho_{\mathrm{c,\bullet}}\sim{3\times
10^{5}}\>\mathrm{M_{\odot}\>pc^{-3}}$–${5\times
10^{6}}\>\mathrm{M_{\odot}\>pc^{-3}}$. BHs are the dominant density component
within the central regions of the cluster, $r\lesssim 0.10$ pc. At later times
the central densities decrease due to stellar evolution (stellar mass loss, BH
formation) and dynamical interactions for both stars and BHs. At $t=50$ Myr,
the central stellar density is $\rho_{\mathrm{c,\star}}\sim{5\times
10^{4}}\>\mathrm{M_{\odot}\>pc^{-3}}$ while the BH density is
$\rho_{\mathrm{c,\bullet}}\sim{4\times 10^{5}}\>\mathrm{M_{\odot}\>pc^{-3}}$,
and BHs are the dominating density component within the central part
($r\lesssim 0.22$ pc) of the cluster. Both the stellar and BH density profiles
remain unchanged in the outer parts within the first $t=50$ Myr of the
simulations. The outer density profiles slopes are close to a power-law slope
of $-3$ indicated in Fig. 5, considerably less steep than the initial Plummer
outer parts, which have a power-law slope of $-5$.

###  4.2 Merger trees of star clusters, their VMSs and IMBHs

We have modified the FOF (friends-of-friends) and subfind structure finder
routines (Springel et al., 2001; Dolag et al., 2009) of the GADGET-4
simulation code (Springel et al., 2021) to identify the individual star
clusters and their merger histories from BIFROST output during the during the
hierarchical assembly simulations. As the snapshot format of BIFROST closely
resembles the GADGET-4 HDF5 format, this is very straightforward. Given the
subfind catalogs constructed from the BIFROST snapshots (interval $\Delta
t_{\mathrm{output}}=10^{-2}$ Myr), we construct the star cluster merger trees
and identify massive stars and IMBHs associated with each cluster. The IMBH
progenitor merger trees can be constructed from the BIFROST snapshots and
merger data output alone, supplemented by the information about cluster
membership of the IMBH progenitors from the subfind data.

We present a detailed description of the formation and assembly history of the
IMBHs and their host star clusters in the simulations H1, H2 and H3 in the
following sections 4.2.1, 4.2.2 and 4.2.3. In each simulation, at least two
IMBHs form in the mass range from $m_{\bullet}\sim 105\;M_{\odot}$ up to
$m_{\bullet}\sim 2200\>\mathrm{M_{\odot}}$.

Although the formation histories of the IMBHs in our simulations differ in
details, a number of common characteristics can be readily observed. First,
IMBHs form almost always from a collapsing VMS built in a sequence of mergers
of massive stars. The sequence of stellar mergers eventually leading to a
formation of an IMBH through direct collapse always begins from a star more
massive than $m_{\mathrm{\star}}\sim 80\>\mathrm{M_{\odot}}$. The total number
of mergers in a merger tree of a single IMBH can be up to
$N_{\mathrm{merg}}\sim 100$. Most of the stellar mass in the growing VMSs is
acquired in mergers with stars with masses $10\>\mathrm{M_{\odot}}\lesssim
m_{\star}\lesssim 150\>\mathrm{M_{\odot}}$. The branches other than the main
branch of the merger tree before the IMBH formation are almost always short,
i.e. the stars merging with the growing VMS have typically never merged before
with other stars. There are only two mergers occurring in our simulation
sample in which the colliding stars have both masses over
$m_{\star}>{150}\>\mathrm{M_{\odot}}$. Next, most of the IMBH mass originates
from the VMS merger phase, and IMBHs grow in mass by less than $25\%$–$35\%$
by disrupting main sequence stars and evolved stars, and by merging with
stellar-mass BHs and other IMBHs. More massive IMBHs acquire more mass by
disrupting stars and merging with other BHs compared to less massive IMBHs.
Finally, and most interestingly, the main progenitor stars of the IMBHs do not
originate from the initially most massive star cluster. Instead, the main
cluster accretes IMBHs and evolved VMSs formed in the lower-mass star clusters
as they merge into the central cluster. This somewhat unexpected phenomenon is
examined in section 5.

In simulations H1 and H3 more than a single IMBH ends up in the main star
cluster. The interaction histories of the IMBHs in the main cluster are
complex and rich in dynamical phenomena. These involve IMBH binaries, IMBH
binary member exchanges, hierarchical triple IMBH configurations, strong
triple interactions, gravitational-wave driven mergers of BH-IMBH and IMBH-
IMBH binaries as well as IMBHs escaping the star clusters after strong
gravitational-wave recoil kicks.

The initial masses of the star clusters hosting IMBH progenitors are presented
in Table 4. The table also presents the times $t_{\mathrm{acc}}$ when the
clusters merge with the main central star cluster, and the IMBH mass at that
time. If the IMBH has not formed yet at $t=t_{\mathrm{acc}}$, we provide the
mass of the IMBH descending from the most massive VMS of the cluster at the
time of its collapse.

####  4.2.1 Simulation H1

Figure 6: The merger trees of stars leading to IMBHs (black symbols) in the
simulation H1. IMBHs form through collapse of VMSs which have been built up in
a series of collisions of massive main sequence (blue symbols) and giant
(orange symbols) stars. No IMBHs originate from stars initially at the most
massive, central star cluster H1-A. Three IMBHs with masses ranging from
$m_{\bullet}\sim{350}\>\mathrm{M_{\odot}}$ to
$m_{\bullet}\sim{1500}\>\mathrm{M_{\odot}}$ form in star clusters H1-B–H1-D
less massive than the main cluster and are accreted into the main cluster
between $t=3.0$ Myr and $t=6.5$ Myr. In the main cluster, the IMBHs undergo a
series of IMBH binary and triple phases between until $t=18.5$, when IMBHs
with masses of $m_{\bullet}={348}\>\mathrm{M_{\odot}}$ and
$m_{\bullet}={1194}\>\mathrm{M_{\odot}}$ merge. The resulting IMBH with a mass
of $m_{\bullet}={1498}\>\mathrm{M_{\odot}}$ receives a gravitational-wave
recoil kick of $v_{\mathrm{recoil}}=206$ km/s, unbinding it from the star
cluster. The final IMBH of the cluster is also ejected at later times after a
merger with a stellar mass BH, as well as the remnants of two low-mass IMBHs
formed though mergers of two stellar-mass BHs. Thus, no IMBHs remain in the
cluster H1-A at the end of the simulation at $t=50$ Myr. The complex merger
histories of the IMBH progenitors are described in more detail in the text.

The merger trees of stars leading to IMBHs and their host star clusters of the
simulation H1 are presented in Fig. 6. The mass of the central star cluster in
the run is initially $M_{\star}={2.0\times 10^{5}}\>\mathrm{M_{\odot}}$. Only
four stars originating from the main cluster H1-A merge with stars more
massive than $m_{\star}\gtrsim{10}\>\mathrm{M_{\odot}}$, and the mass of these
merger remnant stars never exceeds $m_{\star}={150}\>\mathrm{M_{\odot}}$.
Thus, the most massive stars from the cluster H1-A end up in the (P)PISN mass
gap, leaving behind no massive remnant at the end of their lives.

The IMBHs in the simulation H1 are formed in three smaller star clusters
labeled H1-B (IMBH label BH1-B), H1-C (BH1-C) and H1-D (BH1-D) with decreasing
IMBH mass. The initial stellar masses of the clusters are
$M_{\star}={4.4\times 10^{4}}\>\mathrm{M_{\odot}}$, $M_{\star}={1.1\times
10^{5}}\>\mathrm{M_{\odot}}$ and $M_{\star}={5.9\times
10^{4}}\>\mathrm{M_{\odot}}$, respectively. The progenitor of the IMBH BH1-B
in the star cluster H1-B is a massive star with an initial mass of
$m_{\star}={128}\>\mathrm{M_{\odot}}$. After the mass segregation phase,
stellar collisions between the main progenitor and other massive stars begin
at the center of the cluster at $t\gtrsim 1.8$ Myr. Before $t=2.5$ Myr, the
star has already merged with five massive main sequence (MS) stars with
${63}\>\mathrm{M_{\odot}}\lesssim m_{\star}\lesssim{134}\>\mathrm{M_{\odot}}$.
At this point the star reaches the giant phase with
$m_{\mathrm{VMS}}={620}\>\mathrm{M_{\odot}}$. In the giant phase, the star
further merges with five massive MS stars. When the star reaches the end of
its life at $t=2.73$ Myr, the mass of the resulting IMBH BH1-B is
$m_{\bullet}={1049}\>\mathrm{M_{\odot}}$. Later, the IMBH tidally disrupts a
massive giant and two massive MS stars. The mass of the IMBH BH1-B is
$m_{\bullet}={1194}\>\mathrm{M_{\odot}}$ at $t=6.5$ Myr when the star cluster
H1-B merges with the central cluster.

The second-most-massive IMBH of the simulation (BH1-C) originates from a
progenitor of $m_{\star}={134}\>\mathrm{M_{\odot}}$ in a somewhat more massive
star cluster ($M_{\mathrm{\star}}={1.1\times 10^{5}}\>\mathrm{M_{\odot}}$).
When the star begins its giant phase at $t=2.58$ Myr it has a mass of
$m_{\star}={397}\>\mathrm{M_{\odot}}$, having merged with one MS star with
$m_{\star}={35}\>\mathrm{M_{\odot}}$ and two MS stars with
$m_{\star}>{100}\>\mathrm{M_{\odot}}$. In the giant phase there are no mergers
with massive stars. The star collapses into the IMBH BH1-C
($m_{\bullet}={388}\>\mathrm{M_{\odot}}$) at $t=2.81$ Myr. Before $t=3.0$ Myr,
the IMBH tidally disrupts two massive giants
($m_{\star}={119}\>\mathrm{M_{\odot}}$ and
$m_{\star}={179}\>\mathrm{M_{\odot}}$). The latter giant is one of the two
cases in which the star merging with the main progenitor is itself a stellar
merger remnant with $m_{\star}>{150}\>\mathrm{M_{\odot}}$. This star
originates from an earlier merger of two massive MS stars with
$m_{\star}={101}\>\mathrm{M_{\odot}}$ and
$m_{\star}={81}\>\mathrm{M_{\odot}}$. After these tidal disruption and
accretion events the mass of the BH1-C is
$m_{\bullet}={536}\>\mathrm{M_{\odot}}$. At this point the host star cluster
merges with the main central cluster.

The original progenitor of the lowest-mass IMBH of the simulation (BH1-D) has
a mass of $m_{\star}={130}\>\mathrm{M_{\odot}}$ in a star cluster H1-D with a
mass of $M_{\star}={5.9\times 10^{4}}\>\mathrm{M_{\odot}}$. The star has its
first mergers earlier than the two more massive IMBH progenitors, reaching a
mass of $m_{\star}={313}\>\mathrm{M_{\odot}}$ already by $t=1.6$ Myr after
merging with two massive MS stars of $m_{\star}\sim{90}\>\mathrm{M_{\odot}}$.
The star evolves into the giant phase and later collapses into an IMBH with
$m_{\bullet}={292}\>\mathrm{M_{\odot}}$. The host cluster now merges with the
main star cluster of the simulation, after which BH1-D disrupts a massive
giant reaching the mass of $m_{\star}={348}\>\mathrm{M_{\odot}}$ at $t=2.90$
Myr.

After the star cluster H1-C merges with the main cluster, the IMBHs BH1-C and
BH1-D form a bound IMBH binary at the center of the main star cluster. The
binary remains at the center of the star cluster for about $3$ Myr until the
most massive IMBH of the simulation, BH1-B, is accreted into the main star
cluster at $t=6.5$ Myr. At this point, the IMBHs form a hierarchical triple
system with the inner binary consisting of IMBHs BH1-C
($m_{\bullet}={542}\>\mathrm{M_{\odot}}$) and BH1-D
($m_{\bullet}={348}\>\mathrm{M_{\odot}}$) with the third body being BH1-B
($m_{\bullet}={1194}\>\mathrm{M_{\odot}}$). The hierarchical IMBH triple is
not stable (Mardling & Aarseth, 2001), and an exchange event occurs soon at
$t=7.0$ Myr when BH1-B replaces BH1-C in the inner binary. Next, an
evolutionary phase lasting from $t=7.0$ Myr until $t=18.5$ Myr begins,
illustrated in the bottom-right box in Fig. 6. In this phase, the outer IMBH
($m_{\bullet}={542}\>\mathrm{M_{\odot}}$) is repeatedly ejected into a wide
orbit (typical apocenters at $r_{\mathrm{a}}\sim 0.1$ pc – $1.0$ pc) in the
star cluster by the inner binary in three-body interactions. In the
interaction, the inner binary shrinks towards smaller separations. The outer
third body sinks back into the center of the star cluster on a time-scale of
$t_{\mathrm{sink}}\lesssim 1$ – $2$ Myr due to dynamical friction. The cycle
consisting of interaction, ejection and sinking phases repeats multiple times.
This cluster-aided pump-like process transfers the orbital energy of the inner
binary to the cluster stars until the inner binary finally reaches the regime
in which gravitational-wave emission becomes dynamically important to the
system. The IMBHs BH1-B ($m_{\bullet}={1194}\>\mathrm{M_{\odot}}$) and BH1-D
($m_{\bullet}={348}\>\mathrm{M_{\odot}}$) merge into an IMBH with
$m_{\bullet}={1498}\>\mathrm{M_{\odot}}$ at $t=18.5$ Myr. At the time of the
merger, the remnant receives a gravitational-wave recoil kick of
$v_{\mathrm{recoil}}=206$ km/s, which is enough to unbind the IMBH from the
central star cluster. The third IMBH, BH1-C sinks into the center of the
massive star cluster and remains there until $t=31.1$ Myr at which point it
merges with a stellar-mass BH ($m_{\bullet}={52}\>\mathrm{M_{\odot}}$). Even
though the gravitational-wave recoil kick is relatively weak,
$v_{\mathrm{recoil}}=60$ km/s, it is enough to marginally unbind the merger
remnant from the star cluster. At later times, two low-mass IMBHs
($m_{\bullet}={116}\>\mathrm{M_{\odot}}$ and
$m_{\bullet}={135}\>\mathrm{M_{\odot}}$) form in mergers of stellar-mass BHs
in the cluster, but both are instantly ejected from the cluster with recoil
kick velocities of $v_{\mathrm{recoil}}=1102$ km/s and
$v_{\mathrm{recoil}}=140$ km/s, respectively. At the end of the simulation at
$t=50$ Myr, the star cluster H1-A contains no IMBHs and only evolving stars
and interacting stellar-mass BHs remain.

####  4.2.2 Simulation H2

Figure 7: The merger trees of IMBHs and their host star clusters in the
simulation H2. As in simulation H1, the initially most massive star cluster
does not produce an IMBH, but two IMBHs with masses of
$m_{\bullet}={350}\>\mathrm{M_{\odot}}$ and
$m_{\bullet}={1256}\>\mathrm{M_{\odot}}$ form in lower-mass star clusters. The
host star cluster H2-C of the lower-mass IMBH has a close encounter with the
main cluster at $t\sim 4$ Myr and is tidally stripped, but does not merge with
the main cluster. The more massive IMBH with its host cluster H2-B is accreted
by the main cluster, where it repeatedly merges with stellar-mass black holes.
The resulting gravitational-wave recoil kicks are weak,
$v_{\mathrm{recoil}}<25$ km/s, and the IMBH of mass
$m_{\bullet}={1437}\>\mathrm{M_{\odot}}$ is retained in the main cluster at
the end of the simulation at $t=50$ Myr.

The IMBH and star cluster merger trees of the simulation H2 are presented in
Fig. 7. As in the case of the simulation H1, the initially most massive
central star cluster H2-A ($M_{\mathrm{\star}}={2.3\times
10^{5}}\>\mathrm{M_{\odot}}$) does not itself form an IMBH, but rather
accretes an IMBH formed in the less massive cluster H2-B
($M_{\mathrm{\star}}={3.2\times 10^{4}}\>\mathrm{M_{\odot}}$) with its IMBH
BH2-B at $t=4.8$ Myr when $m_{\bullet}={1256}\>\mathrm{M_{\odot}}$. A third
star cluster H2-C forms an IMBH BH2-C with a mass of
$m_{\bullet}={350}\>\mathrm{M_{\odot}}$, but its host star cluster is never
accreted into the main cluster, although the two clusters have a close
encounter around $t\sim 3.6$ Myr in which the smaller star cluster is tidally
stripped and loses $\sim 30\%$ of its stellar mass.

The progenitor of the BH2-B is a star with
$m_{\star}={149}\>\mathrm{M_{\odot}}$ in the initial conditions. In the MS
phase, this star merges with $10$ other MS stars with masses ranging from
$m_{\star}={28}\>\mathrm{M_{\odot}}$ to $m_{\star}={124}\>\mathrm{M_{\odot}}$
between $t=1.42$ Myr and $t=2.09$ Myr. When the star reaches the giant phase
at $t=2.29$ Myr, it has a mass of
$m_{\mathrm{VMS}}={934}\>\mathrm{M_{\odot}}$. In the giant phase, the star
only merges with a single massive MS star, but accretes $29$ low-mass stars
before collapsing into an IMBH with $m_{\bullet}={1132}\>\mathrm{M_{\odot}}$
at $t=2.52$ Myr. Before the IMBH reaches the main star cluster at $t=4.8$ Myr,
it disrupts three massive MS stars. In the main cluster, the IMBH BH2-B
repeatedly merges with stellar-mass black holes, but the resulting
gravitational-wave recoil kicks are weak ($v_{\mathrm{recoil}}\lesssim 24$
km/s) and unable to unbind the IMBH from its host cluster. The final mass of
the IMBH BH2-B is $m_{\bullet}={1437}\>\mathrm{M_{\odot}}$ at $t=50$ Myr.

The progenitor to the BH2-C is a star with
$m_{\star}={111}\>\mathrm{M_{\odot}}$ which has a relatively simple merger
history. Before the beginning of its giant phase at $t=2.76$ Myr, the star
merges with four MS stars with masses between
$m_{\star}={33}\>\mathrm{M_{\odot}}$ and $m_{\star}={73}\>\mathrm{M_{\odot}}$,
reaching a mass of $m_{\star}={334}\>\mathrm{M_{\odot}}$. The resulting IMBH
forms at $t=2.97$ Myr with a mass of $m_{\bullet}={330}\>\mathrm{M_{\odot}}$.
The IMBH later accretes a single MS star, and its final mass is
$m_{\bullet}={350}\>\mathrm{M_{\odot}}$ at the end of the simulation. At this
point, the IMBH is located at the center of the star cluster
($M_{\mathrm{\star}}={1.4\times 10^{4}}\>\mathrm{M_{\odot}}$) in which it
originally formed.

####  4.2.3 Simulation H3

Figure 8: The merger trees of IMBH progenitors and their host star clusters in
the simulation H3. As in the other two hierarchical assembly simulations, the
stars from the initially most massive star clusters are not the progenitors of
the IMBHs. The IMBH BH3-B with $m_{\bullet}={2201}\>\mathrm{M_{\odot}}$ formed
in a large number of collision with MS and giant stars mainly in only $0.9$
Myr between $t=2.0$ Myr and $t=2.9$ Myr is the most massive IMBH formed in the
simulations of this study. By $t=8.7$ Myr, five IMBHs with eventual final
masses $m_{\bullet}={2201}\>\mathrm{M_{\odot}}$,
$m_{\bullet}={347}\>\mathrm{M_{\odot}}$,
$m_{\bullet}={338}\>\mathrm{M_{\odot}}$,
$m_{\bullet}={289}\>\mathrm{M_{\odot}}$ and
$m_{\bullet}={105}\>\mathrm{M_{\odot}}$ of have ended up in the main cluster.
At $t=20.5$ Myr, a strong three-body interaction among the IMBHs ejects the
IMBH with $m_{\bullet}={338}\>\mathrm{M_{\odot}}$ from the cluster. While a
gravitational-wave recoil kick ($v_{\mathrm{recoil}}=1212$ km/s) after a
merger of two stellar-mass BHs ejects a new low-mass remnant IMBH
($m_{\bullet}={105}\>\mathrm{M_{\odot}}$) from the cluster, mergers between
the most massive IMBH and stellar-mass BHs only result in weak recoil kicks.
At $t=30$ Myr, four IMBHs populate the main star cluster, the most massive
cluster at the center and the other three in the outskirts ($r\lesssim 10$ pc)
of the cluster. At later times, the two most massive IMBHs form a binary, and
by the end of the simulation $t=50$ Myr the two other remaining IMBHs
($m_{\bullet}={105}\>\mathrm{M_{\odot}}$ and
$m_{\bullet}={289}\>\mathrm{M_{\odot}}$) have sunk into separations less than
$r\lesssim 2$ pc from the central binary.

The IMBH and star cluster merger trees of the simulation H3 are presented in
Fig. 8. The simulation H3 has the most complex IMBH formation histories of the
three simulations of this study, as in total five star clusters produce an
IMBH. As in the case of the simulations H1 and H2, the most massive star
cluster H3-A ($M_{\mathrm{\star}}={2.3\times 10^{5}}\>\mathrm{M_{\odot}}$) in
the simulation does not itself produce an IMBH, but rather accretes IMBHs or
their VMS progenitors of IMBHs from lower-mass star clusters (${1.4\times
10^{4}}\>\mathrm{M_{\odot}}\lesssim M_{\mathrm{\star}}\lesssim{1.1\times
10^{5}}\>\mathrm{M_{\odot}}$).

The progenitor of the most massive IMBH of this study (BH3-B) is a star with a
mass of $m_{\star}={117}\>\mathrm{M_{\odot}}$. Its host star cluster H3-B
($M_{\mathrm{\star}}={1.1\times 10^{5}}\>\mathrm{M_{\odot}}$) merges with the
central star cluster already at $t=2.20$ Myr, at which point the star has
reached a mass of $m_{\mathrm{VMS}}\sim{600}\>\mathrm{M_{\odot}}$ after
merging with four massive MS stars with masses exceeding
$m_{\star}={96}\>\mathrm{M_{\odot}}$, but is still in the MS phase. Before
reaching the giant phase at $t=2.62$ Myr, the VMS further merges with seven MS
stars, from which one star is a VMS ($m_{\star}={245}\>\mathrm{M_{\odot}}$)
formed earlier in a merger of two massive MS stars with masses of
$m_{\star}={129}\>\mathrm{M_{\odot}}$ and
$m_{\star}={119}\>\mathrm{M_{\odot}}$. At the beginning of the giant phase,
the BH3-B progenitor VMS has a mass of
$m_{\mathrm{VMS}}={1141}\>\mathrm{M_{\odot}}$. In the giant phase, the star
merges with three more MS stars and three massive giants before collapsing
into an IMBH with a mass of $m_{\bullet}={1748}\>\mathrm{M_{\odot}}$ at
$t=2.84$ Myr. The IMBH BH3-B grows beyond
$m_{\bullet}={2000}\>\mathrm{M_{\odot}}$ by disrupting mainly massive giant
stars.

The progenitors of IMBHs from star clusters H3-C, H3-D and H3-E have similar
merger and evolutionary histories. The initial masses of their host star
clusters are $M_{\star}={1.5\times 10^{4}}\>\mathrm{M_{\odot}}$,
$M_{\star}={2.3\times 10^{4}}\>\mathrm{M_{\odot}}$ and $M_{\star}={1.4\times
10^{4}}\>\mathrm{M_{\odot}}$, respectively. Most mergers of the IMBH
progenitors occur already in the MS phase with relatively few mergers in the
giant and IMBH phases. The masses of the IMBHs are
$m_{\bullet}={347}\>\mathrm{M_{\odot}}$ (BH3-C),
$m_{\bullet}={338}\>\mathrm{M_{\odot}}$ (BH3-D), and
$m_{\bullet}={289}\>\mathrm{M_{\odot}}$ (BH3-E) when their host clusters merge
with the central star cluster at $t=8.7$ Myr, $t=3.2$ Myr and $t=6.9$ Myr,
respectively. The lowest-mass IMBH formed through collisions is BH3-F which
originates from a star cluster with
$M_{\mathrm{\star}}={2.6\times{10^{4}}}\>\mathrm{M_{\odot}}$. The star evolves
in isolation until losing its envelope in the giant phase in a collision and
leaves a stellar-mass remnant of $m_{\bullet}={58}\>\mathrm{M_{\odot}}$.
However, the stellar-mass BH disrupts two massive MS stars, and reaches
$m_{\bullet}={105}\>\mathrm{M_{\odot}}$ before being accreted into the main
cluster at $t=7.6$ Myr.

Between $t=8.7$ Myr and $t=20.5$ Myr, the main star cluster H3-A hosts in
total five IMBHs. In addition, a pair of stellar-mass BHs with masses of
$m_{\bullet}={65}\>\mathrm{M_{\odot}}$ and
$m_{\bullet}={61}\>\mathrm{M_{\odot}}$ merges in the cluster, leaving behind a
low-mass IMBH with $m_{\bullet}={106}\>\mathrm{M_{\odot}}$. However, the
remnant receives a strong gravitational-wave recoil kick of
$v_{\mathrm{recoil}}=1212$ km/s, immediately unbinding the remnant from the
star cluster. At $t=20.5$ Myr, BH3-D ($m_{\bullet}={338}\>\mathrm{M_{\odot}})$
is ejected from the star cluster after a strong three-body encounter with
IMBHs BH3-B ($m_{\bullet}={2094}\>\mathrm{M_{\odot}}$) and BH3-E
($m_{\bullet}={289}\>\mathrm{M_{\odot}}$). The most massive IMBH of the
cluster also merges with two stellar-mass BHs, but the resulting
gravitational-wave recoil kicks are weak ($v_{\mathrm{recoil}}\lesssim 10$
km/s) and thus unable to move the IMBH away from the cluster. At $t=30$ Myr
the star cluster hosts four IMBHs: BH3-B
($m_{\bullet}={2150}\>\mathrm{M_{\odot}}$) at the center and three IMBHs
(BH3-C, BH3-E and BH3-F) orbiting at the outskirts of the cluster ($r\lesssim
10$ pc). The masses of the three IMBHs are at this point
$m_{\bullet}={347}\>\mathrm{M_{\odot}}$,
$m_{\bullet}={289}\>\mathrm{M_{\odot}}$ and
$m_{\bullet}={105}\>\mathrm{M_{\odot}}$, respectively. At later times, the
most massive IMBH BH3-B merges with another stellar-mass BH, reaching its
final mass of $m_{\bullet}={2201}\>\mathrm{M_{\odot}}$. The two most massive
IMBHs (BH3-B and BH3-C) form a binary, and by the end of the simulation the
two other remaining IMBHs in the star cluster have sunk into separations less
than $r\lesssim 2$ pc from the central binary.

###  4.3 Three categories of early IMBH growth histories

Figure 9: The mass growth histories of the ten IMBHs that reached $m_{\bullet}\gtrsim{100}\>\mathrm{M_{\odot}}$ by $t=5$ Myr in the simulations H1, H2 and H3. There are three approximate categories as elaborated in Table 5: the low-mass IMBHs with $m_{\bullet}\sim{100}\>\mathrm{M_{\odot}}$ originating below the mass gap, the numerous just-above-the-mass-gap IMBHs with masses $m_{\bullet}\sim{300}\>\mathrm{M_{\odot}}$–${550}\>\mathrm{M_{\odot}}$, and finally the more massive ($m_{\bullet}\gtrsim{1200}\>\mathrm{M_{\odot}}$) IMBHs experiencing growth also after the MS phase. category | imbh-lowmass | imbh-just-above-massgap | imbh-massive
---|---|---|---
| (IMBH-LM) | (IMBH-JAM) | (IMBH-M)
summary | initial BH mass below the mass gap, | IMBH mass just above the mass gap, | MS end mass well above the mass gap
| growth by TDEs or BH mergers | no substantial later growth | substantial growth in the giant and IMBH phases
MS end mass | $m_{\star}\lesssim{150}\>\mathrm{M_{\odot}}$ |  $m_{\star}\sim{300}\>\mathrm{M_{\odot}}$–${400}\>\mathrm{M_{\odot}}$ | $m_{\star}\gtrsim{600}\>\mathrm{M_{\odot}}$
growth after MS | TDEs and/or BH mergers required | typically no substantial growth | substantial growth by up to a factor of $\sim 2$
until $t=50$ Myr |  |  |
IMBH mass at |  $m_{\bullet}\sim{100}\>\mathrm{M_{\odot}}$–${200}\>\mathrm{M_{\odot}}$ |  $m_{\bullet}\sim{300}\>\mathrm{M_{\odot}}$–${550}\>\mathrm{M_{\odot}}$ | $m_{\bullet}\gtrsim{1200}\>\mathrm{M_{\odot}}$
$t=5$ Myr |  |  |
IMBHs of | BH3-F | BH3-E, BH1-D, BH2-C, | BH1-B, BH2-B, BH3-B, BH1-C (?)
this study |  | BH3-C, BH3-D, BH1-C (?) |
Table 5: The three categories of IMBHs formed in the hierarchical star cluster
assembly simulations and their basic formation pathways. The IMBHs formed in
the recent direct-summation N-body simulation studies of massive, isolated
star clusters of Rizzuto et al. (2021); Rizzuto et al. (2022) and Arca Sedda
et al. (2023c) would correspond the categories IMBH-LM and IMBH-JAM of this
classification.

We present the mass growth histories of stars which became IMBHs
($m_{\bullet}>{100}\>\mathrm{M_{\odot}}$) during the first $t=5$ Myr in the
hierarchical clustering simulations H1–H3 in Fig. 9. There are in total ten
such IMBHs in the three simulations, as detailed in sections 4.2.1–4.2.3.
Although the stars growing by repeated stellar mergers rapidly enter the VMS
phase, their merger histories do not classify as a true collisional runaways
as the merger rate does not monotonically increase after each stellar merger.
The ten growth histories can be qualitatively divided into three categories:
low-mass IMBHs in the mass gap (IMBH-LM), IMBHs just above the mass gap (IMBH-
JAM), and massive IMBHs well above the mass gap (IMBH-M). We summarize these
three categories in Table 5.

In the lowest-mass category IMBH-LM, the IMBH mass threshold of
$m_{\bullet}={100}\>\mathrm{M_{\odot}}$ is reached from below. The only
particular example of the category during the first $t=5$ Myr is the IMBH
BH3-F. The progenitor star of BH3-F experiences sudden mass loss in a common
envelope event before its core collapses into a stellar-mass BH, but later
disrupts two massive stars, reaching the IMBH mass range from below.

In the second category IMBH-JAM, the growing IMBH progenitor stars have grown
by mergers into masses of
$m_{\mathrm{\star}}\sim{290}\>\mathrm{M_{\odot}}$–${400}\>\mathrm{M_{\odot}}$
already in the main sequence. The stars do not grow in mass in their giant
phases, and collapse into IMBHs above the mass gap. IMBHs BH3-E, BH1-D, BH2-C,
BH3-C and BH-D with masses
$m_{\bullet}={289}\>\mathrm{M_{\odot}}$–${350}\>\mathrm{M_{\odot}}$ fall into
this category. The final (the most massive) IMBH of the category IMBH-M,
BH1-C, already shows some characteristics of the third and the most massive
category IMBH-Mas its mass increases by almost $\sim{150}\>\mathrm{M_{\odot}}$
via disrupting stars, reaching a mass of
$m_{\bullet}={536}\>\mathrm{M_{\odot}}$ by $t=5$ Myr.

In the high-mass category IMBH-M, the VMSs have reached a mass of
$m_{\mathrm{VMS}}\gtrsim{620}\>\mathrm{M_{\odot}}$ by the end of their MS
life-time. The category is characterized by high initial mass and substantial
mass growth in the giant and early IMBH phases, up to a factor of $\sim 2$
from their mass at the end of the MS. The category contains three IMBHs:
BH1-B, BH2-B and BH3-B, each being the most massive object in its simulation
with IMBH masses ranging from $m_{\bullet}\sim{1200}\>\mathrm{M_{\odot}}$ to
$m_{\bullet}\sim{2000}\>\mathrm{M_{\odot}}$ at $t=5$ Myr. This highlights the
fact that even though the growth histories of the most massive objects do not
technically classify as runaway growth histories, there is a large difference
between the most massive and the second-most-massive objects formed in the
clusters.

###  4.4 No IMBHs form in isolated comparison simulations

Cluster | $N$ | $M_{\star}$ | $a$ | $\gamma_{\mathrm{EFF}}$ | $\rho_{\mathrm{c,\star}}$
---|---|---|---|---|---
|  | $\mathrm{[M_{\odot}]}$ | [pc] |  | [$\mathrm{M_{\odot}}$ pc-3]
C3 | $10^{6}$ | $5.9\times 10^{5}$ | $0.3$ | $2.25$ | $9.5\times 10^{5}$
C50 | $10^{6}$ | $5.9\times 10^{5}$ | $1.0$ | $2.25$ | $4.1\times 10^{4}$
Table 6: The EFF model parameters $M_{\star}$ (total mass), $a$ (scale radius)
and $\gamma_{\mathrm{EFF}}$ (projected outer power-law slope) with the
resulting central stellar density $\rho_{\mathrm{c,\star}}$ for the isolated
comparison simulations C3 and C50. We set the outer cut-off radius to
$r_{\mathrm{max}}=100$ pc.

We now examine whether IMBHs form in isolated massive star clusters with
masses and central densities comparable to the final, hierarchically assembled
clusters with two comparison simulations C3 and C50. We focus on the final,
central star cluster H1-A from the simulation H1. We select two times, $t=3$
Myr at which point the central stellar density is the highest, and the end of
the simulation $t=50$ Myr. The stellar surface density profiles of the central
clusters at these times are well described by the Elson-Freeman-Fall (EFF)
profile (Elson et al., 1987) with the outer slope of
$\gamma_{\mathrm{EFF}}=2.25$. This corresponds to a three-dimensional density
profile slope of $3.25$. We generate the two EFF initial models C3 and C50
using the McLuster code (Küpper et al., 2011) which has the EFF profile
readily available. We use the Nuker (Lauer et al., 1995) template of McLuster
with $\gamma_{\mathrm{Nuker}}=0.0$, $\beta_{\mathrm{Nuker}}=2.25$ and
$\alpha_{\mathrm{Nuker}}=2.0$. The EFF model parameters for the two comparison
simulation setups are listed in Table 6. The central stellar densities
$\rho_{\mathrm{c,\star}}={9.5\times 10^{5}}\>\mathrm{M_{\odot}\>pc^{-3}}$ (C3)
and $\rho_{\mathrm{c,\star}}={4.1\times 10^{4}}\>\mathrm{M_{\odot}\>pc^{-3}}$
(C50) in the two comparison setups are similar than in Fig. 5 for the cluster
H3-A. The stellar IMF is the same as in the hierarchical assembly simulations.
Even though the initial density profiles of the comparison models correspond
to the cluster H1-A at $t=3$ Myr and $t=50$ Myr, respectively, each star in
the comparison simulation initial conditions is initialized to the beginning
of its MS phase.

We run the two comparison simulations with BIFROST for either $t=20$ Myr (C3)
or $t=30$ Myr (C50) at which point the massive stars required for VMS growth
have reached the end of their lives. We find that neither run forms IMBHs. The
maximum stellar mass reached through stellar collisions in the simulation C3
is $m_{\star}={231}\>\mathrm{M_{\odot}}$, which falls into the (P)PISN mass
gap and the massive star leaves no massive remnant. The maximum BH mass
reached in C3 is $m_{\bullet}={89}\>\mathrm{M_{\odot}}$. In the simulation C50
due to relatively low central density, only five stellar mergers during the
simulation, and only one of which involves a massive star with
$m_{\star}={116}\>\mathrm{M_{\odot}}$ accreting a low-mass star. Thus, the BHs
formed in the simulation all have masses below the mass gap with maximum BH
mass being $m_{\bullet}={63}\>\mathrm{M_{\odot}}$.

Based on the hierarchical star cluster assembly simulations H1, H2, H3 and the
isolated comparison simulations C3 and C50, it seems that key to IMBH
formation in the simulated massive star clusters is the hierarchical formation
process itself. Massive, isolated comparison clusters alone with monolithic
origins could not produce IMBHs, even at corresponding initial central
densities.

##  5 IMBHs from isolated star clusters

Cluster ($N_{\mathrm{seed}}$) | $N$ | $M_{\mathrm{\star}}$ | $\max(m_{\star})$ | $\max(m_{\bullet})$
---|---|---|---|---
|  | $\mathrm{[M_{\odot}]}$ | $\mathrm{[M_{\odot}]}$ | $\mathrm{[M_{\odot}]}$
I1 ($10$) | $1.7\times 10^{4}$ | $9.4\times 10^{3}$ | ${197}$ | ${72}$
I2 ($10$) | $2.6\times 10^{4}$ | $1.5\times 10^{4}$ | ${290}$ | ${87}$
I3 ($10$) | $3.9\times 10^{4}$ | $2.3\times 10^{4}$ | ${419}$ | ${378}$
I4 ($10$) | $6.3\times 10^{4}$ | $3.7\times 10^{4}$ | ${932}$ | ${964}$
I5 ($10$) | $1.0\times 10^{5}$ | $5.9\times 10^{4}$ | ${982}$ | ${1065}$
I6 ($10$) | $1.6\times 10^{5}$ | $9.4\times 10^{4}$ | ${1429}$ | ${1510}$
I7 ($10$) | $2.5\times 10^{5}$ | $1.5\times 10^{5}$ | ${1494}$ | ${1672}$
I8 ($20$) | $4.0\times 10^{5}$ | $2.3\times 10^{5}$ | ${362}$ | ${91}$
I9 ($5$) | $6.2\times 10^{5}$ | $3.7\times 10^{5}$ | ${199}$ | ${81}$
I10 ($5$) | $1.0\times 10^{6}$ | $5.9\times 10^{5}$ | ${150}$ | ${63}$
Table 7: The ten isolated star cluster merger sets I1–I10, in total $100$
isolated simulations with their initial particle numbers $N$, initial masses
$M_{\star}$ and the masses for the most massive star ($\max{m_{\star}}$) and
(IM)BH $(\max{m_{\bullet}}$) formed through collisions in the simulation sets.
Figure 10: The maximum stellar (left panel) and the BH (right panel) masses
$\max{(m_{\star})}$ and $\max{(m_{\bullet})}$ reached in the isolated star
clusters during the first $t=5$ Myr of the simulations. The solid black line
indicates the expectation without mergers from the Kroupa (2001) initial mass
function a maximum mass of $m_{\star}={150}\>\mathrm{M_{\odot}}$ including the
suppression of high-mass stars in low-mass clusters (Yan et al., 2023). The
maximum final individual stellar mass grows linearly as
$\max{(m_{\star})}\approx 0.02\times M_{\star}$ up to star cluster masses of
$M_{\star}\sim{10^{5}}\>\mathrm{M_{\odot}}$, after which the VMS growth in
more massive clusters is strongly suppressed. The maximum BH mass follows the
trends of the maximum stellar mass, processed by stellar evolution. Below
cluster masses $M_{\mathrm{\star}}\lesssim{2\times
10^{4}}\>\mathrm{M_{\odot}}$, no BHs with masses above the (P)PISN mass gap
(in gray shading) form. The maximum reached IMBH masses are
$m_{\bullet}\sim{1630}\>\mathrm{M_{\odot}}$. At cluster masses higher than
$M_{\star}\gtrsim{2\times 10^{5}}\>\mathrm{M_{\odot}}$ there are no IMBHs
above the (P)PISN mass, but BHs in the mass gap below the IMBH limit
$m_{\bullet}={100}\>\mathrm{M_{\odot}}$ still form. Note that the maximum IMBH
mass can exceed the maximum VMS mass as the IMBHs still grow via disrupting
stars in the remnant phase.

###  5.1 Maximum VMS and BH masses in isolated simulations

In the second part of the study we examine, given a density profile shape and
a mass-size relation, which star clusters can form IMBHs. The motivation for
the study is twofold. First, the IMBH formation efficiency
$m_{\bullet}/M_{\star}$ has not been extensively studied in the literature
along a star cluster mass-size relation. Second, the initially most massive
central star clusters in the hierarchical star cluster assembly simulations
did not form IMBHs from their stars but rather accreted them from lower-mass
star clusters, and we wish to understand the physical background of this
phenomenon.

Several studies (e.g. Fujii & Portegies Zwart 2013) have pointed out that the
merger rate of massive stars can be enhanced only if their mass segregation
time-scale is shorter than the life-times of massive stars ($t\sim 2$–$3$
Myr). Even for our most massive star cluster models the mass segregation time-
scale from Eq. (4) is $t_{\mathrm{seg}}\sim 1$ Myr. Thus, massive stars in
massive star clusters should in principle be able to collide with each other.

One can argue that based on virial arguments (e.g. Naab et al. 2009) that the
mean density of the assembling star cluster keeps decreasing, especially if
the accreted systems are small with a fixed total accreted mass. This
consequently leads to fewer stellar collisions. While the decreased mean
cluster density may contribute to the lack of collisions in the central
clusters, we show in this study that VMS growth is suppressed even in
isolated, massive star cluster models.

We setup isolated star cluster initial conditions using the same structural
parameters as for the individual clusters in the hierarchical assembly runs
with ten different masses logarithmically spaced in the particle number range
of $1.65\times 10^{4}\lesssim N\lesssim 1.03\times 10^{6}$. We sample ten
realizations of each of the models I1-I7 with $N\leq 2.53\times 10^{5}$, 20
realizations of the model I8 with $N=4.0\times 10^{5}$ and five realizations
of the two more massive setups I9 and I10. The model I8 is examined the most
as its mass lies close to the threshold after which the central clusters did
not form their own IMBHs in the hierarchical assembly runs. In total, the
sample consists of $100$ isolated star cluster models. The basic properties of
the star cluster models I1-I10 are listed in Table 7. We run each model for
$t=5.0$ Myr using the BIFROST code at which point stars with initial masses
above $m_{\mathrm{\star}}\gtrsim{42}\>\mathrm{M_{\odot}}$ have reached the end
of their lives. In the structured star cluster assembly runs no VMS grew from
a star with an initial mass below $m_{\star}\lesssim{80}\>\mathrm{M_{\odot}}$,
so the most efficient period of VMS growth through stellar collisions is
expected to be over at this point.

Stellar collisions occur in almost all isolated simulations. The maximum
reached stellar ($\max{m_{\star}}$) and black hole masses
($\max{m_{\bullet}}$) in the isolated simulations are presented in Table 7 and
in Fig. 10. The maximum VMS mass is $m_{\star}={197}\>\mathrm{M_{\odot}}$ in
the lowest-mass cluster I1, and increases following a linear trend
$\max{(m_{\star})}\propto M_{\star}$ until the trend flattens around models
I5-I6 with cluster masses $M_{\star}\sim{6\times 10^{4}}\>\mathrm{M_{\odot}}$
– ${9\times 10^{4}}\>\mathrm{M_{\odot}}$. We find that the mass of the most
massive star rarely exceeds $\max{(m_{\star})}\approx 0.02\times M_{\star}$ in
the setups I1-I6, corresponding to a maximum IMBH formation efficiency of
$2\%$.

Fig. 10 also shows the maximum stellar mass in the initial models (solid black
line). At high cluster masses the initial maximum stellar mass is constant,
$m_{\star}={150}\>\mathrm{M_{\odot}}$. At lower star cluster masses high-mass
stars are not initially allowed, following Yan et al. (2023). The
$\max{(m_{\star})}\approx 0.02\times M_{\star}$ trend continues from model I1
to model I6, even though the suppression of massive stars in the initial
conditions only occurs in the models I1–I3. In addition, the relation for the
initial cluster mass – maximum stellar mass is shallower than the linear
trend, so the maximum stellar masses in the initial conditions do not alone
explain the relation between the cluster masses and the mass of the final
collision product. The simple explanation is that the stellar collision rates
also depend on the cluster mass, as will be later shown.

The maximum VMS mass of $\max{(m_{\star})}\sim{1500}\>\mathrm{M_{\odot}}$ is
reached in the setup I7 with $M_{\star}={1.5\times
10^{5}}\>\mathrm{M_{\odot}}$. At higher cluster masses, the suppression of the
VMS growth is apparent, and in the highest-mass cluster model I10
($M_{\star}={5.9\times 10^{5}}\>\mathrm{M_{\odot}}$) the stars do not grow
beyond their initial mass of $m_{\star}={150}\>\mathrm{M_{\odot}}$.

The maximum BH mass $m_{\bullet}$ as a function of their host star cluster
mass $M_{\star}$ follows a similar trend as the maximum stellar mass
$m_{\star}$, additionally processed by stellar evolution, especially the
(P)PISN mass gap. IMBHs above the (P)PISN mass gap
($m_{\bullet}\gtrsim{290}\>\mathrm{M_{\odot}}$) form in cluster setups from I3
($M_{\star}={2.3\times 10^{4}}\>\mathrm{M_{\odot}}$) to I7
($M_{\star}={1.5\times 10^{5}}\>\mathrm{M_{\odot}}$), in total $26$ IMBHs. The
maximum reached IMBH mass is $m_{\bullet}={1672}\>\mathrm{M_{\odot}}$. Note
that most massive IMBHs are more massive than the most massive VMSs as the
IMBHs further grow by disrupting stars after the collapse of the VMS. BHs in
the (P)PISN mass gap form in star clusters across a wide range in cluster
mass. Out of the $100$ simulated isolated clusters, $25$ form a BH in the mass
gap below the IMBH mass limit of $m_{\bullet}={100}\>\mathrm{M_{\odot}}$. Only
two IMBHs are located in the mass gap itself in the IMBH mass range between
${100}\>\mathrm{M_{\odot}}\lesssim
m_{\bullet}\lesssim{290}\>\mathrm{M_{\odot}}$. The VMS growth and collapse
IMBH formation channel is approximately an order of magnitude more efficient
in producing IMBH above the mass gap than in the mass gap above
$m_{\bullet}\gtrsim{100}\>\mathrm{M_{\odot}}$.

###  5.2 Which mergers grow the VMSs?

Figure 11: The normalized cumulative mass function of the stars merged to the
growing VMS in the simulation set I7. The individual mergers are shown in gray
and black symbols, the color hue being different in each ten simulations. The
mass function (yellow solid line) is comparable to the Kroupa IMF (orange
solid line) at low masses below $m_{\star}\lesssim{2}\>\mathrm{M_{\odot}}$.
There are relatively few accreted stars in the range of a few tens of Solar
masses. At high masses beyond $m_{\star}\gtrsim{50}\>\mathrm{M_{\odot}}$, the
mass function has a large overabundance of massive stars (yellow and blue
solid lines) compared to the Kroupa IMF. Only $10\%$ of the total accreted
mass is provided by stars less massive that
$m_{\star}\lesssim{53}\>\mathrm{M_{\odot}}$, and less than $1\%$ by stars with
mass less than $m_{\star}\lesssim{5}\>\mathrm{M_{\odot}}$.

Next, we study the merger histories and the mass budget of the VMSs, focusing
on the simulation set I7 which produced the most massive VMSs (and IMBHs) in
the isolated simulations of this study. As in the structured star cluster
assembly runs, typically only one star grows beyond
$m_{\star}={150}\>\mathrm{M_{\odot}}$ through collisions. The progenitor of
the star growing via stellar collisions is always a massive star with a mass
of at least $m_{\star}\sim{119}\>\mathrm{M_{\odot}}$. The mean mass of initial
VMS progenitors in the set I7 is $m_{\star}\sim{131}\>\mathrm{M_{\odot}}$.

We show the cumulative accreted stellar mass from the accreted stars in Fig.
11. In each ten simulations of the set I7, we extract the masses
($m_{\mathrm{s}}$) of the stars merged into the growing VMS, sort them, and
calculate the normalized cumulative mass function. We also show for comparison
the cumulative mass function of the standard Kroupa (2001) piece-wise power-
law IMF. The masses $m_{\mathrm{s}}$ show a somewhat bimodal distribution,
with the lower-mass end being comparable to the Kroupa IMF up to
$m_{\star}\sim{1}\>\mathrm{M_{\odot}}$–${2}\>\mathrm{M_{\odot}}$. Beyond this,
there are relatively few mergers with stars with masses from
$m_{\star}\sim{1}\>\mathrm{M_{\odot}}$–${2}\>\mathrm{M_{\odot}}$ to
$m_{\star}\sim{40}\>\mathrm{M_{\odot}}$–${50}\>\mathrm{M_{\odot}}$. At large
stellar masses ($m_{\star}\gtrsim{50}\>\mathrm{M_{\odot}}$), the mass
distribution of the accreted stars strongly deviates from the Kroupa IMF as a
consequence of mass segregation in the simulated star clusters.

While the mergers between the growing VMS and low-mass stars are common, they
do not significantly contribute the total mass of the VMS. We find that
mergers with stars with masses less than
$m_{\star}\lesssim{5}\>\mathrm{M_{\odot}}$ only contribute $\sim 1\%$ to the
total mass accreted by the growing VMS, and the contribution of stars with
masses $m_{\star}\lesssim{53}\>\mathrm{M_{\odot}}$ is $\sim 10\%$. In total
$50\%$ of the mass budget of the VMS is provided by massive stars with masses
above $m_{\star}\gtrsim{114}\>\mathrm{M_{\odot}}$. Thus, VMS growth and IMBH
formation is mainly fueled by the most massive stars of the IMF in our
simulations.

###  5.3 Stellar mergers: mass and orbit type classification

Figure 12: Classifying the number of stellar mergers $N_{\mathrm{merg}}$ (per
simulation) in the isolated simulations according to the orbit type and masses
of merging stars. Mergers from bound orbits are shown using solid lines and
mergers from unbound orbits with dashed lines. The mass threshold between the
low-mass (L) and the high-mass stars (H) in this analysis is
$m_{\star}={10}\>\mathrm{M_{\odot}}$. With increasing star cluster mass
$M_{\star}$, mergers from unbound orbits become more common while mergers from
any bound orbits (H-L bound and H-H bound, solid orange and blue lines) are
suppressed above $M_{\star}\gtrsim{1.5\times 10^{5}}\>\mathrm{M_{\odot}}$. As
the numbers of massive stars merging from unbound orbits (H-H unbound, dashed
blue line) is in general low, this leads to a suppression of mergers between
massive stars (H-H). Most common mergers at high cluster masses are unbound
H-L mergers (dashed orange line) while low-mass stars merge with each other
exclusively from unbound orbits (L-L unbound, in yellow).

We classify all the mergers occurred in the isolated star cluster simulations
based on the type of their merger orbit (bound or unbound) and the masses of
the merged stars. The mass threshold between the low-mass (L) and the high-
mass stars (H) in this analysis is $m_{\star}={10}\>\mathrm{M_{\odot}}$. Most
mergers occur from close to parabolic orbits. Approximately $93\%$ of the
mergers occur from orbits with eccentricities between $0.9\lesssim e\lesssim
1.1$, and $\sim 64\%$ from eccentricities between $0.99\lesssim e\lesssim
1.01$. The semi-major axis distribution of the bound merger orbits is well
described by a log-normal distribution with a mean of
$\mu_{\mathrm{log10(a/pc)}}=-2.78$ (corresponding to $a\sim 1.7\times 10^{-3}$
pc) and $\sigma_{\mathrm{log10(a/pc)}}=0.82$, i.e. the typical semi-major axis
of a bound mergers orbit is in the milliparsec range.

The merger mass classification simply divides the stars in to higher-mass (H)
and lower-mass (L) stars (H-L threshold mass
$m_{\star}={10}\>\mathrm{M_{\odot}}$), yielding three merger mass categories:
H-H between massive stars, mixed H-L mergers, and L-L between low-mass stars.
The merger mass and orbit type classification is presented in Fig. 12. First,
we note that the total number of mergers per simulation $N_{\mathrm{merg}}$
increases as a function of cluster mass $M_{\star}$. The trend for
$N_{\mathrm{merg}}$ is not a simple power-law, so the increased number of
mergers in massive clusters cannot be purely explained by the fact that more
massive clusters simply have more stars. Second, mergers from unbound orbits
become increasingly more common compared to mergers from bounds orbits as the
cluster mass $M_{\star}$ increases. At low cluster masses
($M_{\star}\lesssim{4\times 10^{4}}\>\mathrm{M_{\odot}}$) the most common
mergers are bound H-H and H-L mergers, and unbound H-L mergers, up to on
average $N_{\mathrm{merg}}\sim 6.2$ of each of such mergers per simulation.
High-mass stars can also merge with each other from unbound orbits, but this
is more rare, only up to $N_{\mathrm{merg}}\sim 1.7$ unbound H-H mergers per
simulation. Bound L-L mergers never occur in the isolated simulation sample.

Between star cluster masses $M_{\star}={4\times 10^{4}}\>\mathrm{M_{\odot}}$
and $M_{\star}={1.5\times 10^{5}}\>\mathrm{M_{\odot}}$, the number of bound
H-H and H-L mergers only mildly increases, never exceeding
$N_{\mathrm{merg}}\sim 11$ mergers per simulation. Meanwhile, the number of
unbound H-L mergers steadily grows, reaching $N_{\mathrm{merg}}=27$ mergers
per simulation, becoming the most common type of mergers at this cluster mass
range. The total number of unbound mergers exceed the total number of bound
mergers at $M_{\star}\sim{6\times{10^{4}}}\>\mathrm{M_{\odot}}$. Unbound L-L
mergers occur when the mass of the star cluster is $M_{\star}\gtrsim{4\times
10^{4}}\>\mathrm{M_{\odot}}$, and their number monotonically increases as a
function of the cluster mass.

In the massive star clusters
$M_{\star}\gtrsim{1.5\times{10^{5}}}\>\mathrm{M_{\odot}}$ the unbound mergers
follow their trends from the lower-mass clusters, but the number of bound
mergers steeply decreases. The number of bound H-H mergers decreases by over
an order of magnitude, while the H-L merger channel is somewhat less affected.
In the highest-mass ($M_{\star}={5.9\times 10^{5}}\>\mathrm{M_{\odot}}$)
isolated cluster simulations, $\sim 77\%$ of the total mergers are unbound H-L
mergers, $17\%$ unbound L-L mergers, and the remaining $6\%$ bound H-L and
bound H-H mergers. Together with the results in Fig. 10 and Fig. 11, it
appears that the reason for the suppressed VMS growth at high cluster masses
is the suppression of the bound massive-massive merger channel at cluster
masses above $M_{\star}\gtrsim{1.5\times{10^{5}}}\>\mathrm{M_{\odot}}$.

###  5.4 The time of onset of stellar mergers

Figure 13: The tuning-fork diagram for the stellar merger times, i.e. the time
of the first collision $t_{\mathrm{coll}}$ as a function of the isolated star
cluster mass. The first collision times between massive stars
($m_{\mathrm{p}}>{80}\>\mathrm{M_{\odot}}$,
$m_{\mathrm{s}}>{10}\>\mathrm{M_{\odot}}$) are shown in orange while we show
the first collision times of any stars in blue. Collisions occur on average
earlier in more massive clusters due to higher stellar densities. However,
$t_{\mathrm{coll}}$ for the massive stars decouples from this trend at star
cluster masses above $M_{\star}\gtrsim{4\times 10^{4}}\>\mathrm{M_{\odot}}$.

Next, we study the the merger times of stars, specifically the time of the
first stellar collision $t_{\mathrm{coll}}$ of the isolated simulations as a
function of the star cluster mass $M_{\star}$. The results of the first
collision time analysis are presented in Fig. 13.

First, focusing on mergers with any primary ($m_{\mathrm{p}}$) and secondary
($m_{\mathrm{s}}$) stellar masses, we find that the mean time for the first
merger $t_{\mathrm{coll}}$ decreases with increasing cluster mass, i.e.
massive clusters begin their stellar mergers earlier. The mean
$t_{\mathrm{coll}}$ decreases from $t_{\mathrm{coll}}\sim 2.49\pm 1.02$ Myr in
the simulation sample I1 ($M_{\star}=9.4\times 10^{3}$) to
$t_{\mathrm{coll}}\sim 0.19\pm 0.22$ Myr for the most massive cluster I10
($M_{\star}=5.9\times 10^{5}$) in the isolated sample. This can be attributed
to the fact that for the chosen initial star cluster mass-size relation the
central density $\rho_{\mathrm{c}}$ monotonically increases with increasing
cluster mass $M_{\star}$.

However, studying the collision times of massive primary stars
$(m_{\mathrm{p}}>{80}\>\mathrm{M_{\odot}})$ with massive secondaries
$(m_{\mathrm{s}}>{10}\>\mathrm{M_{\odot}})$ we find a qualitatively different
trend. Overall, the collision times of massive stars and all stars in Fig. 13
form a tuning fork like diagram. At low cluster masses
($M_{\star}\lesssim{4\times 10^{4}}\>\mathrm{M_{\odot}}$) the average time for
the first merger of a massive star with a massive secondary corresponds to the
average $t_{\mathrm{coll}}$ for all stars, although with a somewhat smaller
scatter. After $M_{\star}\gtrsim{4\times 10^{4}}\>\mathrm{M_{\odot}}$, the
first massive star collision times increase, decoupling from the decreasing
trend of all stars. The average $t_{\mathrm{coll}}$ for the massive stars
increases from $1.54\pm 0.37$ Myr in the sample I4 ($M_{\star}\lesssim{4\times
10^{4}}\>\mathrm{M_{\odot}}$) to the single merger at $t=2.84$ Myr in the most
massive isolated cluster sample I10. This is very close to the life-time of
the most massive stars in the simulations
($m_{\star}={150}\>\mathrm{M_{\odot}}$ with $t=2.82$ Myr) at the chosen
metallicity $Z=0.01\>Z_{\odot}$.

###  5.5 The stellar-dynamical environments of massive stars

Figure 14: The stellar number densities $n$ and velocity dispersions $\sigma$
around massive ($m_{\star}>{80}\>\mathrm{M_{\odot}}$) stars at $t=2.5$ Myr in
different isolated star cluster simulations I1-I10 indicated by the color
scheme. Massive stars in a binaries ($a<0.01$ pc) with massive secondaries
($m_{\mathrm{s}}>{10}\>\mathrm{M_{\odot}}$) are indicated with open circles
with more binaries occurring in the lower-mass setups. The highest number
densities and velocity dispersions occur at the cluster centers. Increasing
cluster mass from set I1 to I10 results in the increase of the central number
density from $n\sim 1$–${2\times 10^{5}}\>\mathrm{\>pc^{-3}}$ to $n\sim
2$–${3\times 10^{6}}\>\mathrm{\>pc^{-3}}$ and central velocity dispersion from
$\sigma\sim 12$ km/s up to $\sigma\sim 61$ km/s. The four most massive setups
I7-I10 have comparable central number densities $n$, but the central velocity
dispersion in the set I7 is only $\sim 41$ km/s compared to $\sigma\sim 61$
km/s in I10.

Finally, we examine the stellar-dynamical environments of massive stars
characterized by the stellar number density $n$ and the velocity dispersion
$\sigma$ in their immediate vicinity at $t=2.5$ Myr when most massive stars
evolve away from the MS. The stellar number density $n$ and the velocity
dispersion $\sigma$ are calculated from $N_{\mathrm{ngb}}=100$ closest
neighbors of each massive star. The results are displayed in Fig. 14.

The fraction of massive stars in bound binary systems ($a<0.01$ pc) is larger
in lower-mass star clusters at $t=2.5$ Myr. The stellar number densities and
velocity dispersions around the massive stars are the highest at the cluster
centers, as expected. From the lowest mass isolated setup I1
($M_{\star}=9.4\times 10^{3}$) to the most massive isolated clusters I10
($M_{\star}=5.9\times 10^{5}$) the central stellar number density increases
from $n\sim 1$–${2\times 10^{5}}\>\mathrm{\>pc^{-3}}$ to $n\sim 2$–${3\times
10^{6}}\>\mathrm{\>pc^{-3}}$ and central velocity dispersion increases from
$\sigma\sim 12$ km/s up to $\sigma\sim 61$ km/s. The four most massive setups
(I7–I10) have comparable central stellar number densities. However, their
maximum velocity dispersions differ considerably, from $\sim 41$ km/s in the
set I7 to $\sim 61$ km/s in the most massive isolated models I10.

The similar central number densities but differing velocity dispersions
indicate that the centers of simulated star clusters I7 ($M_{\star}={1.5\times
10^{5}}\>\mathrm{M_{\odot}}$) and I10 ($M_{\star}={5.9\times
10^{5}}\>\mathrm{M_{\odot}}$) are very different dynamical environments. This
is due to the steep scaling of few-body process outcomes and collision rates
as a function of the velocity dispersion. For example, for equal stellar
masses and number densities (but differing velocity dispersions) the three-
body binary formation rate estimate of Goodman & Hut (1993) in Eq. (6) yields
a difference of a factor of $\sim 36$ between the sets I7 and I10. In
addition, if binaries form dynamically, they contribute to the stellar
collisions especially through single-binary interactions (Gaburov et al.,
2008). However, the Heggie-Hills (Heggie, 1975; Hills, 1975) law states that
wide binaries will on average eventually dissolve, and it is more difficult
for a binary with a given semi-major axis to be hard in a stellar-dynamical
environment with a high velocity dispersion.

###  5.6 Suppressing collisions between massive stars in massive star clusters

Figure 15: The wide-hard binary threshold ($v_{\mathrm{orb}}=\sigma$) for
binaries of different total masses (background color) as a function of their
semi-major axis $a$ and the velocity dispersion $\sigma$ of their environment.
Dynamically formed binaries typically occur in the range $5\times 10^{-4}$ pc
$\lesssim a\lesssim 10^{-2}$ pc (yellow dashed horizontal lines). The velocity
dispersions $\sigma=10$ km/s and $\sigma=65$ km/s (dash-dotted gray vertical
lines) bracket the central velocity dispersions of isolated models I1-I10. The
wide-hard threshold is shown for two binaries with
$m_{\mathrm{b}}=m_{\mathrm{p}}+m_{\mathrm{s}}={10}\>\mathrm{M_{\odot}}$ (blue
line) and $m_{\mathrm{b}}={300}\>\mathrm{M_{\odot}}$ (orange line). Very
massive ($m_{\mathrm{b}}={300}\>\mathrm{M_{\odot}}$) dynamically formed
binaries in the simulations cannot harden in the region A of the parameter
space, but will harden in region B. With increasing velocity dispersion
$\sigma$ it is increasingly difficult for a binary with a given mass to be
located in the region B in which they can harden. Finally, low-mass binaries
can never harden, explaining the complete lack of mergers between low-mass
star from bound orbits in the isolated simulations.

We have investigated the origin of the suppressed VMS growth and inhibited
IMBH formation initially found in the hierarchical star cluster assembly
simulations through analyzing a set of $100$ isolated star cluster
simulations. We summarize our findings from the isolated simulations as
follows.

  * •

The IMBH formation proceeds through a collapse of a VMS formed in a sequence
of stellar collisions. IMBHs predominantly form in our simulations in the star
cluster mass window between $M_{\star}\sim{2\times
10^{4}}\>\mathrm{M_{\odot}}$ and $M_{\star}\sim{2\times
10^{5}}\>\mathrm{M_{\odot}}$. No more than one star per cluster grows
considerably in mass by mergers.

  * •

The vast majority of the VMS mass is acquired through mergers with massive
stars. Half of the VMS mass originates from stars more massive than
$m_{\star}\gtrsim{114}\>\mathrm{M_{\odot}}$, and only $10\%$ from stars less
massive than $m_{\star}\lesssim{50}\>\mathrm{M_{\odot}}$.

  * •

$93\%$ of all collisions occur from almost parabolic orbits, $0.9\lesssim
e\lesssim 1.1$. The number of collisions between stars of any mass steadily
increases as a function of $M_{\star}$. Collisions between massive stars
($m_{\mathrm{p}}>{80}\>\mathrm{M_{\odot}}$,
$m_{\mathrm{s}}>{10}\>\mathrm{M_{\odot}}$) occur predominantly from bound
orbits. These bound massive-massive stellar collisions are suppressed in
isolated star cluster models with $M_{\star}\gtrsim{1.5\times
10^{5}}\>\mathrm{M_{\odot}}$.

  * •

In more massive star cluster models, the first collision occurs on average
earlier with an almost monotonic trend. At low cluster masses, the average
first collision time $t_{\mathrm{coll}}$ between a massive star
($m_{\mathrm{p}}\gtrsim{80}\>\mathrm{M_{\odot}}$) and a massive secondary
($m_{\mathrm{s}}\gtrsim{10}\>\mathrm{M_{\odot}}$) closely follows the
$t_{\mathrm{coll}}$ between any stars. At $M_{\star}\gtrsim{4\times
10^{4}}\>\mathrm{M_{\odot}}$, the first collision time between massive stars
decouples from the general trend, and actually increases as a function of
$M_{\star}$. At the highest cluster masses, the average time for the first
massive-massive collision is comparable to the life-time of the most massive
stars in our isolated simulations.

  * •

Finally, the central stellar number densities $n$ are approximately similar in
the four most massive isolated cluster setups I7–I10 with
$M_{\star}\gtrsim{1.5\times 10^{5}}\>\mathrm{M_{\odot}}$. Only the set I7 of
the four setups has sustained VMS growth and frequent IMBH formation. However,
given the very different central velocity dispersions from $\sigma\sim 41$
km/s (I7) to $\sigma\sim 61$ km/s (I10), the centers of the isolated cluster
models with different masses are very different stellar-dynamical environments
due to typically steep scaling of the outcomes of the dynamical few-body
processes as a function of $\sigma$.

The snapshot output interval of the simulations of this study is $\Delta
t_{\mathrm{output}}=10^{-2}$ Myr. Three-body processes responsible for
dynamical binary formation and mergers from binary-single interactions occur
on a considerably faster time-scale, and advanced algorithmic techniques such
as continuous simulation data output stream (Hausammann et al., 2022) would be
needed to capture every relevant particle interaction into the output.
Considering typical dynamically formed binary semi-major axis $a=2\times
10^{-3}$ pc and a high velocity dispersion of $\sigma=50$ km/s, the
interaction time-scale for the binary and a third body is of the order of
$a/\sigma\sim 4\times 10^{-5}$ Myr $\ll\Delta t_{\mathrm{output}}$. Although
our snapshot interval is not frequent enough to resolve interactions
responsible for binary formation, we can track the formation, evolution and
fate of binaries of massive stars using the output data as the first strong
interaction rarely leads to a merger (e.g. Atallah et al. 2024).

We track the most bound binaries ($B=-E\propto
m_{\mathrm{p}}m_{\mathrm{s}}a^{-1}$) throughout the simulation sets I7 (VMS
growth) and I8–I10 (no VMS growth). The primary of the most bound binary is
almost immediately a star close to the maximum initial stellar mass
$m_{\star}={150}\>\mathrm{M_{\odot}}$. Binaries with more massive secondaries
form through either the 3BBF or exchanges. By $t=1.5$ Myr, binaries with
massive secondaries of the order of $m_{\mathrm{s}}$ have formed. In the
lower-mass clusters I7 the most bound binary gradually hardens during the next
few $0.1$ Myr–$0.5$ Myr which together with high binary eccentricities enable
collisions. However, in the higher-mass (and higher velocity dispersion)
clusters I8–I10 some massive binaries still form, but they do not typically
harden considerably in the subsequent interactions.

We argue based on the steep 3BBF velocity dispersion scaling and the Heggie-
Hills law that the high velocity dispersion in more massive star clusters both
obstructs the formation of massive binaries, and hinders their later
hardening, inhibiting the VMS growth and thus IMBH formation. We further
elaborate the argument using Fig. 15. Briefly, the binaries that form
dynamically in our isolated simulations have typical semi-major axis in the
range of $5\times 10^{-4}$ pc $\lesssim a\lesssim 10^{-2}$ pc. With the
minimum and maximum central velocity dispersions of our isolated models ($10$
km/s $\lesssim a\lesssim 65$ km/s), the semi-major axis range spans a
rectangular region in the $(\sigma,a)$ space. We also show the hard-wide
binary threshold for a massive equal-mass binary with
$m_{\mathrm{p}}=m_{\mathrm{s}}={150}\>\mathrm{M_{\odot}}$ in Fig. 15. This
line divides the rectangle which the dynamical binaries inhabit into two parts
we label A (above the hard-wide threshold) and B (below the threshold). As a
consequence of the Heggie-Hills law, binaries in region A will on average
further widen and eventually dissolve while binaries in region B can on
average further harden in subsequent encounters. Most importantly, for a fixed
binary mass, there is a high enough threshold $\sigma$ beyond which a typical
dynamically formed binary can only reside in region A, and thus cannot harden.

The velocity dispersion argument is consistent with the earlier findings in
this section, especially the facts that mergers between massive stars from
bound orbits are suppressed at sufficiently high cluster masses $M_{\star}$,
and if collisions involving massive stars take place, they occur on average at
later times. Finally, the life-times of the most massive stars ($t=2.82$ Myr
for $m_{\star}={150}\>\mathrm{M_{\odot}}$) ensure that the collisions are not
just postponed, but completely inhibited as the massive stars reach the ends
the their lives before they manage to collide with each other.

##  6 Discussion

###  6.1 Hierarchical assembly boosts IMBH formation and reduces its
stochasticity

We have shown that in our isolated simulations, massive star clusters with
high central velocity dispersions $\sigma$ do not form IMBHs via the VMS
collapse channel. In the hierarchical star cluster assembly simulations, this
manifests as a fact that the initially most massive star clusters do not form
IMBHs from their stars, but rather by inheriting them from less massive star
clusters during the star cluster mergers. Isolated comparison simulations with
overall density profiles and central densities comparable to the highest and
the final assembling cluster densities do not form IMBHs either. Thus, IMBH
formation is boosted by the hierarchical star cluster assembly. Intermediate-
mass star clusters can better form IMBHs due to their smaller mass segregation
$t_{\mathrm{seg}}$ and core collapse time-scales $t_{\mathrm{cc}}$ compared to
the life-times of massive stars, as well as their lower central velocity
dispersions.

Dense, isolated star clusters may still form IMBHs if their central velocity
dispersion is not too high. This requires lower star cluster masses
$M_{\star}$. In our simulations, the threshold for IMBH formation is around
$M_{\star}\sim{2\times 10^{5}}\>\mathrm{M_{\odot}}$. However, hierarchical
star cluster assembly with this final cluster mass would still most probably
produce more than one IMBH in the star cluster. Consequently, the IMBH content
of a star cluster with a fixed mass can be drastically different depending on
its formation history (monolithic vs hierarchical). In our simulations, young
star clusters finishing their assembly had $1$–$5$ IMBHs in total.

Finally, the hierarchical star cluster assembly reduces the stochasticity of
IMBH formation. In the isolated simulation sets the cluster models which
formed IMBHs did not do so in every tested random realization. In the
hierarchical assembly, tens of intermediate-mass star clusters contribute to
the final mass of the assembling central cluster. If any of the intermediate-
mass clusters end up producing an IMBH, the assembled star cluster will
contain an IMBH, at least initially.

###  6.2 Surviving the gravitational-wave recoil

Due to the hierarchical cluster assembly, IMBHs will end up in environments
with higher escape velocities than the environment they formed in. This
enhances their prospects of being retained in the star clusters in the likely
mergers with stellar-mass black holes and the resulting gravitational-wave
recoil kicks. As discussed in section 2.4, an IMBH needs to have a mass of at
least $m_{\bullet}\gtrsim{500}\>\mathrm{M_{\odot}}$ to be retained after a
merger with most massive stellar BHs below the mass gap
($m_{\bullet}\sim{63}\>\mathrm{M_{\odot}}$) in a star cluster with
$M_{\star}\sim{10^{6}}\>\mathrm{M_{\odot}}$, and
$m_{\bullet}\gtrsim{1000}\>\mathrm{M_{\odot}}$ in a star cluster of