$M_{\star}\sim{10^{5}}\>\mathrm{M_{\odot}}$. The most massive IMBHs in each of
our cluster assembly simulations have masses above
$m_{\bullet}\gtrsim{1400}\>\mathrm{M_{\odot}}$, so they should in principle be
resilient against being ejected from their host clusters when merging with
stellar-mass black holes.

However, the hierarchical star cluster assembly also poses new challenges for
retaining massive IMBHs in their star clusters. In isolated star clusters,
typically no more than one IMBH is formed. The hierarchical assembly commonly
leads into situations with more than a single IMBH in a single star cluster,
enabling IMBH-IMBH mergers in addition to only IMBH-BH mergers. Even though
IMBHs formed through the VMS collapse are assumed to be slowly spinning, the
gravitational-wave recoil kicks can be strong unless the merging IMBHs have a
mass ratio $q=m_{\mathrm{s}}/m_{\mathrm{p}}$ close to unity, or $q$ needs to
be very small. If the secondary is also in the IMBH mass range, this would
require very massive primary IMBHs.

###  6.3 Stellar mergers, merger mass loss and primordial binary stars

In this work we have assumed for simplicity that colliding stars merge without
mass loss, unless a common envelope event is triggered instead of full mixing.
In reality, stellar collision mass loss may be substantial, especially for
stars with low surface gravities, and environments with high collision
velocities. Thus, our VMS masses and consequently IMBH masses are likely
somewhat overestimated.

However, we do not assume any stellar rejuvenation either. In SEVN, the
stellar merger remnant inherits the phase and percentage of life-time of the
most evolved progenitor star, and the interpolation algorithm of SEVN finds
the new post-merger track self-consistently. As the durations of the stellar
evolutionary phases shorten with increasing stellar mass, the conservative
approach does not lead to rejuvenation. It is common in the literature to
numerically rejuvenate a star after a stellar merger. The age of the merger
product $t(m_{\mathrm{p}}+m_{\mathrm{s}})$ can be calculated e.g. as

$t(m_{\mathrm{p}}+m_{\mathrm{s}})=\frac{m_{\mathrm{p}}}{m_{\mathrm{p}}+m_{\mathrm{s}}}\frac{t_{\mathrm{MS}}(m_{\mathrm{p}}+m_{\mathrm{s}})}{t_{\mathrm{MS}}(m_{\mathrm{p}})}t(m_{\mathrm{p}})$
(15)

in which $t_{\mathrm{MS}}$ is the MS life-time of the star (Meurs & van den
Heuvel, 1989; Portegies Zwart et al., 1999). In later evolutionary phases the
MS life-time $t_{\mathrm{MS}}$ can be replaced by the duration of the
particular phase, and $t(m_{\mathrm{p}})$ by the time already spent in the
particular phase (Mapelli, 2016). As extensive rejuvenation will lead to
increasingly long life-times for the collision products enabling more future
collisions, our choice of no explicit rejuvenation most probably
underestimates the final VMS and IMBH masses. Assessing the overall combined
effect of the rejuvenation and the lack of stellar merger mass loss is
complicated and will require careful further model building and additional
simulations beyond this work.

We note that the hierarchical star cluster assembly has implications for the
stellar collision mass loss as well. The common reasoning is that more massive
stellar systems have higher velocity dispersions $\sigma$, higher average
collision velocities thus on average more collision mass loss. At very high
collision velocities the stars can even be destroyed. There are a number of
potential caveats in the picture. First, a large fraction of collisions may
originate from primordial binary stars, in which case the binary population
properties and not the host star cluster determines the collision velocities.
In the absence of primordial binaries, collisions from interactions between
dynamically formed binaries with single stars are important, producing a low-
velocity tail in the collision velocity distribution with relatively little
expected mass loss (Gaburov et al., 2008). The straightforward effect of the
hierarchical star cluster assembly scenario to the stellar collision mass loss
is the lowered mean collision velocity. The stellar collisions leading to VMS
growth and IMBH formation occur in the lower-mass star clusters with lower
velocity dispersions and thus less collision mass loss is expected in the
hierarchical assembly process compared to monolithic clusters.

Finally, we emphasize the probable importance of primordial binary stars for
the IMBH formation in the hierarchical star cluster assembly scenario.
Primordial binaries were not included in the models of this work, but are
instead the focus of the next study (FROST-CLUSTERS II) of this project.
Primordial binaries enhance the interaction cross-sections ($a$ versus
$r_{\mathrm{\star}}$), potentially increasing the number of stellar
collisions. Primordial binaries can ameliorate the suppressed VMS growth at
high cluster masses in several ways. First, massive stars tend to reside in
binary systems, and a sizable fraction of them are already hard. Second, the
3BBF bottleneck of binary formation with its extreme dependence on the
velocity dispersion may vanish if the primordial binary fraction is high
enough. Further simulations in the FROST-CLUSTERS project will explore the
question in the near future.

###  6.4 Gravitational waves

Figure 16: The primary and secondary masses in BH-BH mergers of the
hierarchical star cluster assembly runs of this study (black, open circles)
compared to the primary and secondary masses from the GWTC-3 catalog (in
yellow) detected by the Advanced LIGO and Advanced Virgo gravitational-wave
observatories up to the end of their third observing run (Abbott et al.,
2023). Curves for constant total binary masses
$m_{\mathrm{b}}=m_{\mathrm{p}}+m_{\mathrm{s}}$ are displayed from
$m_{\mathrm{b}}={10}\>\mathrm{M_{\odot}}$ to
$M_{\mathrm{b}}={10^{3}}\>\mathrm{M_{\odot}}$ in gray, dashed lines. Current
ground-based GW detectors cannot observe mergers with
$m_{\mathrm{b}}\gtrsim{500}\>\mathrm{M_{\odot}}$ very well above $z\gtrsim 1$.

We present the primary $m_{\mathrm{p}}$ and secondary masses $m_{\mathrm{s}}$
of merging BH-BH, IMBH-BH and IMBH-IMBH binaries from our hierarchical star
cluster assembly simulations H1, H2 and H3 in Fig 16. In addition, Fig. 16
shows the primary and secondary masses from the GWTC-3 catalog (in yellow)
detected by the Advanced LIGO and Advanced Virgo gravitational-wave
observatories up to the end of their third observing run (Abbott et al.,
2023). The BH-BH mergers (${61}\>\mathrm{M_{\odot}}\lesssim
m_{\mathrm{p}}\lesssim{86}\>\mathrm{M_{\odot}}$,
${52}\>\mathrm{M_{\odot}}\lesssim
m_{\mathrm{s}}\lesssim{59}\>\mathrm{M_{\odot}}$) during the first $t=50$ Myr
of the simulations agree with the GWTC-3 sources very well.

The most numerous gravitational-wave merger category during the first $t=50$
Myr of the cluster assembly simulations in this study are IMBH-BH mergers.
Most IMBH-BH mergers have their total masses
$m_{\mathrm{b}}=m_{\mathrm{p}}+m_{\mathrm{s}}$ and mass ratios
$q=m_{\mathrm{s}}/m_{\mathrm{p}}$ ranging from
${1300}\>\mathrm{M_{\odot}}\lesssim
m_{\mathrm{b}}\lesssim{2220}\>\mathrm{M_{\odot}}$ and $1.5\times
10^{-2}\lesssim q\lesssim 4.5\times 10^{-2}$. One IMBH-BH merger has a lower
total mass of $m_{\mathrm{b}}\sim{600}\>\mathrm{M_{\odot}}$ and $q\sim 0.1$.
Finally, the only true IMBH-IMBH merger of this study has a total mass of
$m_{\mathrm{b}}\sim{1500}\>\mathrm{M_{\odot}}$ and $q\sim 0.3$. Based on the
literature work (e.g. Rizzuto et al. 2021; Rizzuto et al. 2022), we expect
that the inclusion of primordial binary stars in the initial cluster setups in
our forthcoming study FROST-CLUSTERS II will populate the BH-BH
LIGO/Virgo/KAGRA parameter space in Fig. 16. Beside this regime, our
hierarchical star cluster assembly simulations also indicate a distinct second
GW population in the IMBH mass range around
$m_{\mathrm{b}}\sim{10^{3}}\>\mathrm{M_{\odot}}$ with typical mass ratios of
$1.5\times 10^{-2}\lesssim q\lesssim 4.5\times 10^{-2}$.

Signatures of BH mergers in the IMBH mass range with various mass ratios have
been searched from the data of Advanced LIGO and and Advanced Virgo detectors
from the first, second and third observing runs (Abbott et al., 2017, 2019,
2022). The current ground-based detectors are sensitive to BH-BH mergers at
the lower end of the IMBH mass range as demonstrated by the detection of the
event GW190521 from the formation of an IMBH with
$m_{\bullet}\sim{150}\>\mathrm{M_{\odot}}$ (Abbott et al., 2020). So far, no
signatures of gravitational waves from BH-BH mergers in the IMBH mass range
have been found, and only upper limits on the IMBH merger rates can have been
placed. The searches do not typically extend well beyond total (IM)BH binary
masses of
$m_{\mathrm{b}}=m_{\mathrm{p}}+m_{\mathrm{s}}\sim{500}\>\mathrm{M_{\odot}}$ as
the current ground-based detectors become less sensitive at higher total
binary masses, corresponding to lower frequencies below $\sim 10$ Hz. In
general, more massive merging binaries produce a shorter but louder GW signal
with little inspiral but more ringdown signal at the current observable
frequency range (Abbott et al., 2022).

In general, IMBH mergers with a total mass of
${500}\>\mathrm{M_{\odot}}\lesssim
m_{\mathrm{b}}\lesssim{1000}\>\mathrm{M_{\odot}}$ could be observable with the
current LIGO/Virgo/KAGRA detectors with a sufficient signal-to-noise ratio of
S/N$\sim 8$ to cosmological distances up to $\sim$ a few Gpc, corresponding to
redshifts $z\lesssim 1$ (e.g. Mazzolo et al. 2014; Haster et al. 2016b; Mehta
et al. 2022; Fragione et al. 2022b). For a fixed total mass binary mass
$m_{\mathrm{b}}$, binaries with less massive secondaries $m_{\mathrm{s}}$ have
smaller signal-to-noise ratios and are thus more difficult to detect from
large distances than equal-mass binaries. For IMBH-BH mergers
($m_{\mathrm{p}}\gg m_{\mathrm{s}}$), the detection horizon distance scales
approximately proportional to as
$\left(m_{\mathrm{s}}/m_{\mathrm{p}}\right)^{1/2}$ when
$m_{\mathrm{b}}<{200}\>\mathrm{M_{\odot}}$ (inspiral domination) and with
$m_{\mathrm{s}}/m_{\mathrm{p}}$ at higher binary masses when the ringdown
signal dominates (Haster et al., 2016a, b), assuming Advanced LIGO design
sensitivity. Thus, we expect that most of our IMBH mass scale mergers would
not be detectable with the current GW instruments, as our
$z=0.01\>\mathrm{Z_{\odot}}$ simulations mimic the low-metallicity environment
at $z\gtrsim 6$. However, such mergers from low metallicity pockets at lower
redshifts in the Universe are possible, yet likely rare.

Future ground-based (Einstein Telescope, Cosmic Explorer) and space-borne
detectors such as LISA will be able to detect IMBH mergers with
$m_{\mathrm{b}}\sim{10^{3}}\>\mathrm{M_{\odot}}$ well up to $z\sim 10$ (Jani
et al., 2020). For multi-band detections, the Einstein Telescope + LISA
network will be the most sensitive for binaries with
$m_{\mathrm{b}}\sim{1.5\times 10^{3}}\>\mathrm{M_{\odot}}$–${2.0\times
10^{3}}\>\mathrm{M_{\odot}}$, independent of the mass ratio (Jani et al.,
2020). This is exactly the mass range of the most massive IMBHs of IMBH-M
category in our hierarchical assembly simulations.

##  7 Summary and conclusions

We have used the updated BIFROST N-body simulation code now coupled to stellar
evolution to study the formation of VMSs and IMBHs in hierarchically
assembling massive star clusters. We assume the universal power-law mass
function for the individual young sub-clusters with a slope of $\alpha=-2$
supported both by observations (e.g. Elmegreen & Efremov 1996; Zhang & Fall
1999; Adamo et al. 2020) and state-of-the-art solar-mass-resolution
hydrodynamical simulations of star cluster formation (Lahén et al., 2020).

Besides the observationally motivated shallow power-law slope for the star
cluster mass-size relation ($\alpha=0.18$), we choose small sub-pc birth radii
$r_{\mathrm{h}}$ for our initial sub-cluster populations based on the observed
sizes of the still embedded nascent star clusters (e.g. Marks & Kroupa 2012).
An additional motivation for the small sub-cluster sizes are the recent JWST
observations of several close, dense, parsec-size, massive
($M_{\star}\sim{10^{6}}\>\mathrm{M_{\odot}}$) star clusters in the Cosmic Gems
arc at $z=10.2$ (Adamo et al., 2024). The structure of the initial star
cluster assembly regions is adopted from representative cluster formation
regions of the hydrodynamical low-metallicity dwarf galaxy starburst
simulations of Lahén et al. (2020).

Our main findings are the following. First, IMBHs form in each three
hierarchical cluster assembly simulations during the first $t=50$ Myr with
maximum IMBH masses $m_{\bullet}\sim{2200}\>\mathrm{M_{\odot}}$. The primary
formation channel for IMBH formation is the collapse of a VMS which formed
through a sequence of collisions with massive stars, the primary progenitor
having an initial mass exceeding $m_{\star}\gtrsim{80}\>\mathrm{M_{\odot}}$.
Tidal disruption of stars, BH-BH, IMBH-BH and IMBH-IMBH (in one case) mergers
also contribute to the final mass of the IMBHs, but their contribution is sub-
dominant during the first $50$ Myr. Typically only a single IMBH forms in a
star cluster, but multiple IMBHs will end up in the most massive star cluster
through the cluster assembly process. Multiple IMBHs per final star cluster is
essentially an unique prediction of the model as in isolated setups typically
no more than a single object can considerably grow by stellar mergers.

Second, no IMBHs originate from stars in the initially most massive star
clusters in the hierarchical simulations. Although the initially most massive
star clusters have the highest central stellar densities, they also have
longest mass segregation timescales and highest central velocity dispersions,
impeding the collisions of massive stars and IMBH formation. In addition, when
simulating the evolution of isolated massive star clusters with density
profiles comparable to the hierarchical assembly end products, no IMBHs form.
The key to the presence of IMBHs in assembled massive star clusters with
$M_{\star}>{5\times 10^{5}}\>\mathrm{M_{\odot}}$ appears to be the
hierarchical assembly process itself.

The hierarchical star cluster assembly process has additional important
implications for the VMS growth and IMBH formation. Mass loss in stellar
collisions, though not modeled in this study, may be less severe in
hierarchical formation of star clusters due to lower velocity dispersions of
the star clusters in which the VMSs formed. Moreover, the IMBHs end up in
higher escape velocity environments than their original host star clusters,
making the IMBHs more resilient to gravitational-wave recoil kicks when
inevitably merging with stellar-mass black holes at later times. However, the
fact that massive star clusters in our simulation may host multiple IMBHs
poses an additional challenge to their survival in the clusters, as IMBH-IMBH
mergers have closer to equal mass ratios compared to IMBH-BH mergers, making
strong GW recoil kicks more probable.

The various GW-driven BH-BH and IMBH-BH mergers and the single IMBH-IMBH
merger are detectable with either current LIGO/Virgo/KAGRA GW observatories or
the next-generation ground-based and space-borne detectors (e.g. Jani et al.
2020). Besides the binary mergers consistent with the current LIGO/Virgo/KAGRA
observations, our simulations show another distinct population of GW sources:
IMBHs in with masses $m_{\bullet}\sim{10^{3}}\>\mathrm{M_{\odot}}$–${2\times
10^{3}}\>\mathrm{M_{\odot}}$ merging with stellar-mass BHs yielding mass
ratios of $1.5\times 10^{-2}\lesssim q\lesssim 4.5\times 10^{-2}$.

We have further studied the VMS and IMBH growth in dense, isolated, massive
star clusters in detail in a set of $100$ star cluster simulations containing
up to a million stars. Star clusters with masses $M_{\star}\gtrsim{2\times
10^{4}}\>\mathrm{M_{\odot}}$ form IMBHs above the (P)PISN mass gap. Initially,
the maximum reached VMS masses through collisions scale linearly with the
cluster mass, $\max{(m_{\star})}\approx 0.02\times M_{\star}$, yielding a BH
formation efficiency of $\epsilon_{\bullet}\sim 0.02$. In low-mass clusters,
VMS growth is inhibited by the relatively low central densities and the low
number of massive stars. The maximum IMBH mass in the isolated cluster
simulations during the first $t=5$ Myr is reached at cluster masses
$M_{\star}\sim{1.5\times 10^{5}}\>\mathrm{M_{\odot}}$ with
$m_{\bullet}\sim{1630}\>\mathrm{M_{\odot}}$, somewhat lower than in the
structured cluster assembly runs.

The key results of the isolated star cluster simulations are the following.
First, more massive star clusters have more stellar mergers, and the first
mergers occur on average earlier in more massive clusters. Second, most
mergers occur from bound orbits at low cluster masses
($M_{\star}\lesssim{6\times 10^{4}}\>\mathrm{M_{\odot}}$) while unbound
mergers make up to over $90\%$ of the total mergers at high cluster masses.
Next, most of the VMS mass growth originates from a relatively small number of
mergers of VMS primaries with massive secondaries: stars less massive than
$m_{\star}\sim{50}\>\mathrm{M_{\odot}}$ only contribute less than $\sim 10\%$
to the total mass accreted by the VMSs. VMS-VMS mergers are rare. Mergers
between massive stars ($m_{\star}\gtrsim{10}\>\mathrm{M_{\odot}}$) originate
mostly from bound orbits, and occur later in more massive clusters. The
massive-massive bound mergers become increasingly rare at high cluster masses
above $M_{\star}\gtrsim{2\times 10^{5}}\>\mathrm{M_{\odot}}$. This efficiently
suppresses the VMS growth and thus IMBH formation in isolated massive star
clusters.

The suppression of VMS growth and rapid IMBH formation in massive, isolated
star clusters caused by a combination of two factors. In star clusters with
masses in the range of ${2\times 10^{4}}\>\mathrm{M_{\odot}}\lesssim
M_{\star}\lesssim{2\times 10^{5}}\>\mathrm{M_{\odot}}$ binaries with massive
stars can form through the 3BBF channel and subsequent dynamical three-body
exchanges. As the velocity dispersion of these clusters is relatively low
($\lesssim 40$ km/s), these binaries can harden according to Heggie-Hills law,
and combined with their high eccentricities and low pericenter distances also
eventually merge with each other. In the case of more massive clusters the
velocity dispersion is high and the 3BBF channel is suppressed. If binaries
with massive stars form, they are too wide to harden with further interactions
with the cluster stars. Thus, massive stellar mergers and the VMS growth is
inhibited. Second, the relatively short life-times of massive stars ($t=2.82$
Myr for $m_{\star}={150}\>\mathrm{M_{\odot}}$) ensure that the onset of
massive star mergers is truly prevented and not just postponed as the massive
stars reach the ends of their lives.

Finally, we emphasize that primordial binaries, ignored in this study for
simplicity, may help to counter the suppression of VMS growth observed in star
cluster setups with high velocity by providing a source of massive binary
stars independent of the 3BBF and binary-single exchange channels. We will
explore this in the forthcoming second study of the FROST-CLUSTERS project.

## Data availability statement

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

## Acknowledgments

The authors thank Francesco Paolo Rizzuto and Jacob Kosowski from the
Theoretical Extragalactic Group of the University of Helsinki, Finland. The
numerical simulations were performed using facilities hosted by the Max Planck
Computing and Data Facility (MPCDF) and the Leibniz Supercomputing Centre
(LRZ), Germany. TN acknowledges support from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094 - 390783311 from the DFG Cluster of Excellence
"ORIGINS".

## References

  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations. Cambridge University Press
  * Aarseth & Hills (1972) Aarseth S. J., Hills J. G., 1972, A&A, 21, 255
  * Abbott et al. (2017) Abbott B. P., et al., 2017, Phys. Rev. D, 96, 022001
  * Abbott et al. (2019) Abbott B. P., et al., 2019, Phys. Rev. D, 100, 064064
  * Abbott et al. (2020) Abbott R., et al., 2020, Phys. Rev. Lett., 125, 101102
  * Abbott et al. (2022) Abbott R., et al., 2022, A&A, 659, A84
  * Abbott et al. (2023) Abbott R., et al., 2023, Physical Review X, 13, 041039
  * Adamo et al. (2017) Adamo A., et al., 2017, ApJ, 841, 131
  * Adamo et al. (2020) Adamo A., et al., 2020, Space Sci. Rev., 216, 69
  * Adamo et al. (2024) Adamo A., et al., 2024, arXiv e-prints, p. arXiv:2401.03224
  * Agrawal et al. (2020) Agrawal P., Hurley J., Stevenson S., Szécsi D., Flynn C., 2020, MNRAS, 497, 4549
  * Allison et al. (2009) Allison R. J., Goodwin S. P., Parker R. J., de Grijs R., Portegies Zwart S. F., Kouwenhoven M. B. N., 2009, ApJ, 700, L99
  * Antonini et al. (2019) Antonini F., Gieles M., Gualandris A., 2019, MNRAS, 486, 5008
  * Arca-Sedda et al. (2021) Arca-Sedda M., Rizzuto F. P., Naab T., Ostriker J., Giersz M., Spurzem R., 2021, ApJ, 920, 128
  * Arca Sedda et al. (2023a) Arca Sedda M., et al., 2023a, arXiv e-prints, p. arXiv:2307.04805
  * Arca Sedda et al. (2023b) Arca Sedda M., Kamlah A. W. H., Spurzem R., Rizzuto F. P., Giersz M., Naab T., Berczik P., 2023b, arXiv e-prints, p. arXiv:2307.04807
  * Arca Sedda et al. (2023c) Arca Sedda M., Kamlah A. W. H., Spurzem R., Rizzuto F. P., Naab T., Giersz M., Berczik P., 2023c, MNRAS, 526, 429
  * Ardi et al. (2008) Ardi E., Baumgardt H., Mineshige S., 2008, ApJ, 682, 1195
  * Askar et al. (2017) Askar A., Szkudlarek M., Gondek-Rosińska D., Giersz M., Bulik T., 2017, MNRAS, 464, L36
  * Atallah et al. (2024) Atallah D., Weatherford N. C., Trani A. A., Rasio F., 2024, On Binary Formation from Three Initially Unbound Bodies (arXiv:2402.12429)
  * Banerjee (2018) Banerjee S., 2018, MNRAS, 473, 909
  * Banerjee (2021) Banerjee S., 2021, MNRAS, 500, 3002
  * Banerjee & Kroupa (2017) Banerjee S., Kroupa P., 2017, A&A, 597, A28
  * Bastian et al. (2005) Bastian N., Gieles M., Efremov Y. N., Lamers H. J. G. L. M., 2005, A&A, 443, 79
  * Baumgardt (2017) Baumgardt H., 2017, MNRAS, 464, 2174
  * Baumgardt & Klessen (2011) Baumgardt H., Klessen R. S., 2011, MNRAS, 413, 1810
  * Baumgardt et al. (2003) Baumgardt H., Makino J., Hut P., McMillan S., Portegies Zwart S., 2003, ApJ, 589, L25
  * Baumgardt et al. (2019) Baumgardt H., et al., 2019, MNRAS, 488, 5340
  * Begelman & Rees (1978) Begelman M. C., Rees M. J., 1978, MNRAS, 185, 847
  * Bekenstein (1973) Bekenstein J. D., 1973, ApJ, 183, 657
  * Belczynski & Banerjee (2020) Belczynski K., Banerjee S., 2020, A&A, 640, L20
  * Belczynski et al. (2020) Belczynski K., et al., 2020, A&A, 636, A104
  * Belkus et al. (2007) Belkus H., Van Bever J., Vanbeveren D., 2007, ApJ, 659, 1576
  * Bestenlehner et al. (2011) Bestenlehner J. M., et al., 2011, A&A, 530, L14
  * Bettwieser & Sugimoto (1984) Bettwieser E., Sugimoto D., 1984, MNRAS, 208, 493
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press, Princeton, NJ USA
  * Boekholt et al. (2018) Boekholt T. C. N., Schleicher D. R. G., Fellhauer M., Klessen R. S., Reinoso B., Stutz A. M., Haemmerlé L., 2018, MNRAS, 476, 366
  * Bonnell et al. (1998) Bonnell I. A., Bate M. R., Zinnecker H., 1998, MNRAS, 298, 93
  * Breen & Heggie (2013) Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779
  * Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  * Brown & Gnedin (2021) Brown G., Gnedin O. Y., 2021, MNRAS, 508, 5935
  * Bulirsch & Stoer (1966) Bulirsch R., Stoer J., 1966, Numerische Mathematik, 8, 1
  * Bunker et al. (2023) Bunker A. J., et al., 2023, A&A, 677, A88
  * Campanelli et al. (2007a) Campanelli M., Lousto C. O., Zlochower Y., Merritt D., 2007a, Phys. Rev. Lett., 98, 231102
  * Campanelli et al. (2007b) Campanelli M., Lousto C., Zlochower Y., Merritt D., 2007b, ApJ, 659, L5
  * Casertano & Hut (1985) Casertano S., Hut P., 1985, ApJ, 298, 80
  * Charbonnel et al. (2023) Charbonnel C., Schaerer D., Prantzos N., Ramírez-Galeano L., Fragos T., Kuruvanthodi A., Marques-Chaves R., Gieles M., 2023, A&A, 673, L7
  * Chen et al. (2015) Chen Y., Bressan A., Girardi L., Marigo P., Kong X., Lanza A., 2015, MNRAS, 452, 1068
  * Chilingarian et al. (2018) Chilingarian I. V., Katkov I. Y., Zolotukhin I. Y., Grishin K. A., Beletsky Y., Boutsia K., Osip D. J., 2018, ApJ, 863, 1
  * Chin (1997) Chin S. A., 1997, Physics Letters A, 226, 344
  * Chin (2007) Chin S. A., 2007, arXiv e-prints, p. arXiv:0704.3273
  * Chin & Chen (2005) Chin S. A., Chen C. R., 2005, Celestial Mechanics and Dynamical Astronomy, 91, 301
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Cignoni et al. (2016) Cignoni M., et al., 2016, ApJ, 833, 154
  * Clarke & Bonnell (2008) Clarke C. J., Bonnell I. A., 2008, MNRAS, 388, 1171
  * Cohn (1980) Cohn H., 1980, ApJ, 242, 765
  * Colbert & Mushotzky (1999) Colbert E. J. M., Mushotzky R. F., 1999, ApJ, 519, 89
  * Cook et al. (2019) Cook D. O., et al., 2019, MNRAS, 484, 4897
  * Costa et al. (2021) Costa G., Bressan A., Mapelli M., Marigo P., Iorio G., Spera M., 2021, MNRAS, 501, 4514
  * Crowther et al. (2010) Crowther P. A., Schnurr O., Hirschi R., Yusof N., Parker R. J., Goodwin S. P., Kassim H. A., 2010, MNRAS, 408, 731
  * Dale et al. (2015a) Dale J. E., Ercolano B., Bonnell I. A., 2015a, MNRAS, 451, 987
  * Dale et al. (2015b) Dale J. E., Ercolano B., Bonnell I. A., 2015b, MNRAS, 451, 987
  * Davis et al. (2020) Davis T. A., et al., 2020, MNRAS, 496, 4061
  * Dehnen & Hernandez (2017) Dehnen W., Hernandez D. M., 2017, MNRAS, 465, 1201
  * Dehnen & Read (2011) Dehnen W., Read J. I., 2011, European Physical Journal Plus, 126, 55
  * Della Croce et al. (2024) Della Croce A., Pascale R., Giunchi E., Nipoti C., Cignoni M., Dalessandro E., 2024, A&A, 682, A22
  * Denissenkov & Hartwick (2014) Denissenkov P. A., Hartwick F. D. A., 2014, MNRAS, 437, L21
  * Di Carlo et al. (2021) Di Carlo U. N., et al., 2021, MNRAS, 507, 5132
  * Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
  * Dong et al. (2012) Dong R., Greene J. E., Ho L. C., 2012, ApJ, 761, 73
  * Elmegreen & Efremov (1996) Elmegreen B. G., Efremov Y. N., 1996, ApJ, 466, 802
  * Elson et al. (1987) Elson R. A. W., Fall S. M., Freeman K. C., 1987, ApJ, 323, 54
  * Escala (2021) Escala A., 2021, ApJ, 908, 57
  * Farias et al. (2024) Farias J. P., Offner S. S. R., Grudić M. Y., Guszejnov D., Rosen A. L., 2024, MNRAS, 527, 6732
  * Farrell et al. (2009) Farrell S. A., Webb N. A., Barret D., Godet O., Rodrigues J. M., 2009, Nature, 460, 73
  * Ferrarese & Merritt (2000) Ferrarese L., Merritt D., 2000, ApJ, 539, L9
  * Fowler & Hoyle (1964) Fowler W. A., Hoyle F., 1964, ApJS, 9, 201
  * Fragione & Kocsis (2018) Fragione G., Kocsis B., 2018, Phys. Rev. Lett., 121, 161103
  * Fragione & Silk (2020) Fragione G., Silk J., 2020, MNRAS, 498, 4591
  * Fragione et al. (2018) Fragione G., Leigh N. W. C., Ginsburg I., Kocsis B., 2018, ApJ, 867, 119
  * Fragione et al. (2022a) Fragione G., Kocsis B., Rasio F. A., Silk J., 2022a, ApJ, 927, 231
  * Fragione et al. (2022b) Fragione G., Loeb A., Kocsis B., Rasio F. A., 2022b, ApJ, 933, 170
  * Fragos et al. (2023) Fragos T., et al., 2023, ApJS, 264, 45
  * Fregeau et al. (2004) Fregeau J. M., Cheung P., Portegies Zwart S. F., Rasio F. A., 2004, MNRAS, 352, 1
  * Freitag et al. (2006a) Freitag M., Rasio F. A., Baumgardt H., 2006a, MNRAS, 368, 121
  * Freitag et al. (2006b) Freitag M., Gürkan M. A., Rasio F. A., 2006b, MNRAS, 368, 141
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Fujii & Portegies Zwart (2013) Fujii M. S., Portegies Zwart S., 2013, MNRAS, 430, 1018
  * Fujii & Portegies Zwart (2014) Fujii M. S., Portegies Zwart S., 2014, MNRAS, 439, 1003
  * Fujii et al. (2012) Fujii M. S., Saitoh T. R., Portegies Zwart S. F., 2012, ApJ, 753, 85
  * Gaburov et al. (2008) Gaburov E., Gualandris A., Portegies Zwart S., 2008, MNRAS, 384, 376
  * Gatto et al. (2021) Gatto M., et al., 2021, MNRAS, 507, 3312
  * Gebhardt et al. (2000) Gebhardt K., et al., 2000, ApJ, 539, L13
  * Gebhardt et al. (2002) Gebhardt K., Rich R. M., Ho L. C., 2002, ApJ, 578, L41
  * Giacobbo & Mapelli (2020) Giacobbo N., Mapelli M., 2020, ApJ, 891, 141
  * Gieles et al. (2018) Gieles M., et al., 2018, MNRAS, 478, 2461
  * Giersz & Heggie (1994) Giersz M., Heggie D. C., 1994, MNRAS, 268, 257
  * Glebbeek et al. (2009) Glebbeek E., Gaburov E., de Mink S. E., Pols O. R., Portegies Zwart S. F., 2009, A&A, 497, 255
  * Gold et al. (1965) Gold T., Axford W. I., Ray E. C., 1965, in Robinson I., Schild A., Schucking E. L., eds, Quasi-Stellar Sources and Gravitational Collapse. p. 93
  * González Prieto et al. (2022) González Prieto E., Kremer K., Fragione G., Martinez M. A. S., Weatherford N. C., Zevin M., Rasio F. A., 2022, ApJ, 940, 131
  * González et al. (2007) González J. A., Sperhake U., Brügmann B., Hannam M., Husa S., 2007, Phys. Rev. Lett., 98, 091101
  * Goodman & Hut (1993) Goodman J., Hut P., 1993, ApJ, 403, 271
  * Goswami et al. (2012) Goswami S., Umbreit S., Bierbaum M., Rasio F. A., 2012, ApJ, 752, 43
  * Gragg (1965) Gragg W. B., 1965, SIAM Journal on Numerical Analysis, 2, 384
  * Grasha et al. (2017) Grasha K., et al., 2017, ApJ, 840, 113
  * Greene & Ho (2007) Greene J. E., Ho L. C., 2007, ApJ, 656, 84
  * Greene et al. (2020) Greene J. E., Strader J., Ho L. C., 2020, ARA&A, 58, 257
  * Grudić et al. (2018) Grudić M. Y., Guszejnov D., Hopkins P. F., Lamberts A., Boylan-Kolchin M., Murray N., Schmitz D., 2018, MNRAS, 481, 688
  * Gürkan et al. (2004) Gürkan M. A., Freitag M., Rasio F. A., 2004, ApJ, 604, 632
  * Gürkan et al. (2006) Gürkan M. A., Fregeau J. M., Rasio F. A., 2006, ApJ, 640, L39
  * Gvaramadze & Gualandris (2011) Gvaramadze V. V., Gualandris A., 2011, MNRAS, 410, 304
  * Haster et al. (2016a) Haster C.-J., Wang Z., Berry C. P. L., Stevenson S., Veitch J., Mandel I., 2016a, MNRAS, 457, 4499
  * Haster et al. (2016b) Haster C.-J., Antonini F., Kalogera V., Mandel I., 2016b, ApJ, 832, 192
  * Hausammann et al. (2022) Hausammann L., Gonnet P., Schaller M., 2022, Astronomy and Computing, 41, 100659
  * Heggie (1975) Heggie D. C., 1975, MNRAS, 173, 729
  * Heggie (1993) Heggie D. C., 1993, arXiv e-prints, pp astro–ph/9312018
  * Heggie & Hut (2003) Heggie D., Hut P., 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics. Cambridge University Press
  * Hellström & Mikkola (2010) Hellström C., Mikkola S., 2010, Celestial Mechanics and Dynamical Astronomy, 106, 143
  * Hénault-Brunet et al. (2020) Hénault-Brunet V., Gieles M., Strader J., Peuten M., Balbinot E., Douglas K. E. K., 2020, MNRAS, 491, 113
  * Hill & Zaritsky (2006) Hill A., Zaritsky D., 2006, AJ, 131, 414
  * Hills (1975) Hills J. G., 1975, AJ, 80, 1075
  * Howard et al. (2018) Howard C. S., Pudritz R. E., Harris W. E., 2018, Nature Astronomy, 2, 725
  * Hunter et al. (2003) Hunter D. A., Elmegreen B. G., Dupuy T. J., Mortonson M., 2003, AJ, 126, 1836
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Ibata et al. (2009) Ibata R., et al., 2009, ApJ, 699, L169
  * Iorio et al. (2023) Iorio G., et al., 2023, MNRAS, 524, 426
  * Jalali et al. (2012) Jalali B., Baumgardt H., Kissler-Patig M., Gebhardt K., Noyola E., Lützgendorf N., de Zeeuw P. T., 2012, A&A, 538, A19
  * Jani et al. (2020) Jani K., Shoemaker D., Cutler C., 2020, Nature Astronomy, 4, 260
  * Jiménez-Forteza et al. (2017) Jiménez-Forteza X., Keitel D., Husa S., Hannam M., Khan S., Pürrer M., 2017, Phys. Rev. D, 95, 064024
  * Kaaret et al. (2001) Kaaret P., Prestwich A. H., Zezas A., Murray S. S., Kim D. W., Kilgard R. E., Schlegel E. M., Ward M. J., 2001, MNRAS, 321, L29
  * Kalari et al. (2018) Kalari V. M., Rubio M., Elmegreen B. G., Guzmán V. V., Zinnecker H., Herrera C. N., 2018, ApJ, 852, 71
  * Kamann et al. (2016) Kamann S., et al., 2016, A&A, 588, A149
  * Kamlah et al. (2022) Kamlah A. W. H., et al., 2022, arXiv e-prints, p. arXiv:2205.04470
  * Khorrami et al. (2019) Khorrami Z., et al., 2019, MNRAS, 485, 3124
  * Kim et al. (2020) Kim M., López K. M., Jonker P. G., Ho L. C., Im M., 2020, MNRAS, 493, L76
  * King (1962) King I., 1962, AJ, 67, 471
  * Kochanek (1992) Kochanek C. S., 1992, ApJ, 385, 604
  * Krause et al. (2016) Krause M. G. H., Charbonnel C., Bastian N., Diehl R., 2016, A&A, 587, A53
  * Kremer et al. (2020) Kremer K., Ye C. S., Chatterjee S., Rodriguez C. L., Rasio F. A., 2020, in Bragaglia A., Davies M., Sills A., Vesperini E., eds, Vol. 351, Star Clusters: From the Milky Way to the Early Universe. pp 357–366 (arXiv:1907.12564), doi:10.1017/S1743921319007269
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kruckow et al. (2018) Kruckow M. U., Tauris T. M., Langer N., Kramer M., Izzard R. G., 2018, MNRAS, 481, 1908
  * Küpper et al. (2011) Küpper A. H. W., Maschberger T., Kroupa P., Baumgardt H., 2011, MNRAS, 417, 2300
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA&A, 41, 57
  * Lahén et al. (2019) Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., 2019, ApJ, 879, L18
  * Lahén et al. (2020) Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., Steinwandel U. P., Moster B. P., 2020, ApJ, 891, 2
  * Lanzoni et al. (2013) Lanzoni B., et al., 2013, ApJ, 769, 107
  * Larsen (2004) Larsen S. S., 2004, A&A, 416, 537
  * Lattimer & Yahil (1989) Lattimer J. M., Yahil A., 1989, ApJ, 340, 426
  * Lauer et al. (1995) Lauer T. R., et al., 1995, AJ, 110, 2622
  * Lee (1987) Lee H. M., 1987, ApJ, 319, 801
  * Lin et al. (2018) Lin D., et al., 2018, Nature Astronomy, 2, 656
  * Livermore et al. (2017) Livermore R. C., Finkelstein S. L., Lotz J. M., 2017, ApJ, 835, 113
  * Lützgendorf et al. (2011) Lützgendorf N., Kissler-Patig M., Noyola E., Jalali B., de Zeeuw P. T., Gebhardt K., Baumgardt H., 2011, A&A, 533, A36
  * Lynden-Bell & Eggleton (1980) Lynden-Bell D., Eggleton P. P., 1980, MNRAS, 191, 483
  * Lynden-Bell & Wood (1968) Lynden-Bell D., Wood R., 1968, MNRAS, 138, 495
  * Ma et al. (2020) Ma X., et al., 2020, MNRAS, 493, 4315
  * Makino (1996) Makino J., 1996, ApJ, 471, 796
  * Maliszewski et al. (2022) Maliszewski K., Giersz M., Gondek-Rosinska D., Askar A., Hypki A., 2022, MNRAS, 514, 5879
  * Mann et al. (2019) Mann C. R., et al., 2019, ApJ, 875, 1
  * Mapelli (2016) Mapelli M., 2016, MNRAS, 459, 3432
  * Mapelli et al. (2020) Mapelli M., Spera M., Montanari E., Limongi M., Chieffi A., Giacobbo N., Bressan A., Bouffanais Y., 2020, ApJ, 888, 76
  * Mapelli et al. (2021) Mapelli M., et al., 2021, MNRAS, 505, 339
  * Mardling & Aarseth (2001) Mardling R. A., Aarseth S. J., 2001, MNRAS, 321, 398
  * Marks & Kroupa (2012) Marks M., Kroupa P., 2012, A&A, 543, A8
  * Marques-Chaves et al. (2024) Marques-Chaves R., et al., 2024, A&A, 681, A30
  * Massey & Hunter (1998) Massey P., Hunter D. A., 1998, ApJ, 493, 180
  * Matsumoto et al. (2001) Matsumoto H., Tsuru T. G., Koyama K., Awaki H., Canizares C. R., Kawai N., Matsushita S., Kawabe R., 2001, ApJ, 547, L25
  * Mazzolo et al. (2014) Mazzolo G., et al., 2014, Phys. Rev. D, 90, 063002
  * McMillan et al. (2007) McMillan S. L. W., Vesperini E., Portegies Zwart S. F., 2007, ApJ, 655, L45
  * Mehta et al. (2022) Mehta A. K., Buonanno A., Gair J., Miller M. C., Farag E., deBoer R. J., Wiescher M., Timmes F. X., 2022, ApJ, 924, 39
  * Menon et al. (2021) Menon S. H., et al., 2021, MNRAS, 507, 5542
  * Meurs & van den Heuvel (1989) Meurs E. J. A., van den Heuvel E. P. J., 1989, A&A, 226, 88
  * Mezcua (2017) Mezcua M., 2017, International Journal of Modern Physics D, 26, 1730021
  * Mezcua et al. (2013) Mezcua M., Roberts T. P., Sutton A. D., Lobanov A. P., 2013, MNRAS, 436, 3128
  * Mezcua et al. (2015) Mezcua M., Roberts T. P., Lobanov A. P., Sutton A. D., 2015, MNRAS, 448, 1893
  * Mezcua et al. (2018) Mezcua M., Civano F., Marchesi S., Suh H., Fabbiano G., Volonteri M., 2018, MNRAS, 478, 2576
  * Mikkola (2008) Mikkola S., 2008, in Vesperini E., Giersz M., Sills A., eds, IAU Symposium Vol. 246, Dynamical Evolution of Dense Stellar Systems. pp 218–227, doi:10.1017/S1743921308015639
  * Mikkola & Aarseth (1996) Mikkola S., Aarseth S. J., 1996, Celestial Mechanics and Dynamical Astronomy, 64, 197
  * Mikkola & Aarseth (2002) Mikkola S., Aarseth S., 2002, Celestial Mechanics and Dynamical Astronomy, 84, 343
  * Mikkola & Merritt (2006) Mikkola S., Merritt D., 2006, MNRAS, 372, 219
  * Mikkola & Tanikawa (1999) Mikkola S., Tanikawa K., 1999, MNRAS, 310, 745
  * Miller & Hamilton (2002) Miller M. C., Hamilton D. P., 2002, MNRAS, 330, 232
  * Moeckel & Bonnell (2009) Moeckel N., Bonnell I. A., 2009, MNRAS, 400, 657
  * Moeckel & Clarke (2011) Moeckel N., Clarke C. J., 2011, MNRAS, 410, 2799
  * Morscher et al. (2015) Morscher M., Pattabiraman B., Rodriguez C., Rasio F. A., Umbreit S., 2015, ApJ, 800, 9
  * Naab et al. (2009) Naab T., Johansson P. H., Ostriker J. P., 2009, ApJ, 699, L178
  * Nguyen et al. (2022) Nguyen C. T., et al., 2022, A&A, 665, A126
  * Noyola et al. (2010) Noyola E., Gebhardt K., Kissler-Patig M., Lützgendorf N., Jalali B., de Zeeuw P. T., Baumgardt H., 2010, ApJ, 719, L60
  * Omelyan (2006) Omelyan I. P., 2006, Phys. Rev. E, 74, 036703
  * Pasham et al. (2014) Pasham D. R., Strohmayer T. E., Mushotzky R. F., 2014, Nature, 513, 74
  * Patruno et al. (2006) Patruno A., Portegies Zwart S., Dewi J., Hopman C., 2006, MNRAS, 370, L6
  * Pauldrach et al. (2012) Pauldrach A. W. A., Vanbeveren D., Hoffmann T. L., 2012, A&A, 538, A75
  * Pechetti et al. (2022) Pechetti R., et al., 2022, ApJ, 924, 48
  * Peebles (1972) Peebles P. J. E., 1972, General Relativity and Gravitation, 3, 63
  * Plummer (1911) Plummer H. C., 1911, MNRAS, 71, 460
  * Poisson & Will (2014) Poisson E., Will C. M., 2014, Gravity. Cambridge University Press
  * Pols et al. (1998) Pols O. R., Schröder K.-P., Hurley J. R., Tout C. A., Eggleton P. P., 1998, MNRAS, 298, 525
  * Portegies Zwart & McMillan (2002) Portegies Zwart S. F., McMillan S. L. W., 2002, ApJ, 576, 899
  * Portegies Zwart et al. (1999) Portegies Zwart S. F., Makino J., McMillan S. L. W., Hut P., 1999, A&A, 348, 117
  * Portegies Zwart et al. (2004) Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillan S. L. W., 2004, Nature, 428, 724
  * Portegies Zwart et al. (2010) Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431
  * Preto & Tremaine (1999) Preto M., Tremaine S., 1999, AJ, 118, 2532
  * Quinlan & Shapiro (1990) Quinlan G. D., Shapiro S. L., 1990, ApJ, 356, 483
  * Rantala & Naab (2024) Rantala A., Naab T., 2024, MNRAS, 527, 11458
  * Rantala et al. (2017) Rantala A., Pihajoki P., Johansson P. H., Naab T., Lahén N., Sawala T., 2017, ApJ, 840, 53
  * Rantala et al. (2020) Rantala A., Pihajoki P., Mannerkoski M., Johansson P. H., Naab T., 2020, MNRAS, 492, 4131
  * Rantala et al. (2021) Rantala A., Naab T., Springel V., 2021, MNRAS, 502, 5546
  * Rantala et al. (2023) Rantala A., Naab T., Rizzuto F. P., Mannerkoski M., Partmann C., Lautenschütz K., 2023, MNRAS, 522, 5180
  * Rastello et al. (2021) Rastello S., Mapelli M., Di Carlo U. N., Iorio G., Ballone A., Giacobbo N., Santoliquido F., Torniamenti S., 2021, MNRAS, 507, 3612
  * Reines (2022) Reines A. E., 2022, Nature Astronomy, 6, 26
  * Reinoso et al. (2018) Reinoso B., Schleicher D. R. G., Fellhauer M., Klessen R. S., Boekholt T. C. N., 2018, A&A, 614, A14
  * Reinoso et al. (2020) Reinoso B., Schleicher D. R. G., Fellhauer M., Leigh N. W. C., Klessen R. S., 2020, A&A, 639, A92
  * Rizzuto et al. (2021) Rizzuto F. P., et al., 2021, MNRAS, 501, 5257
  * Rizzuto et al. (2022) Rizzuto F. P., Naab T., Spurzem R., Arca-Sedda M., Giersz M., Ostriker J. P., Banerjee S., 2022, MNRAS, 512, 884
  * Rizzuto et al. (2023) Rizzuto F. P., Naab T., Rantala A., Johansson P. H., Ostriker J. P., Stone N. C., Liao S., Irodotou D., 2023, MNRAS, 521, 2930
  * Rodriguez et al. (2016) Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, Phys. Rev. D, 93, 084029
  * Rodriguez et al. (2019) Rodriguez C. L., Zevin M., Amaro-Seoane P., Chatterjee S., Kremer K., Rasio F. A., Ye C. S., 2019, Phys. Rev. D, 100, 043027
  * Sakurai et al. (2017) Sakurai Y., Yoshida N., Fujii M. S., Hirano S., 2017, MNRAS, 472, 1677
  * Schaerer (2002) Schaerer D., 2002, A&A, 382, 28
  * Selman & Melnick (2013) Selman F. J., Melnick J., 2013, A&A, 552, A94
  * Sills et al. (2018) Sills A., Rieder S., Scora J., McCloskey J., Jaffa S., 2018, MNRAS, 477, 1903
  * Smith et al. (2011) Smith R., Slater R., Fellhauer M., Goodwin S., Assmann P., 2011, MNRAS, 416, 383
  * Spera & Mapelli (2017) Spera M., Mapelli M., 2017, MNRAS, 470, 4739
  * Spera et al. (2015) Spera M., Mapelli M., Bressan A., 2015, MNRAS, 451, 4086
  * Spera et al. (2019) Spera M., Mapelli M., Giacobbo N., Trani A. A., Bressan A., Costa G., 2019, MNRAS, 485, 889
  * Spitzer (1987) Spitzer L., 1987, Dynamical evolution of globular clusters. Princeton University Press, Princeton, N.J.
  * Spitzer & Hart (1971) Spitzer Lyman J., Hart M. H., 1971, ApJ, 164, 399
  * Spitzer & Shull (1975) Spitzer L. J., Shull J. M., 1975, ApJ, 200, 339
  * Spitzer & Stone (1967) Spitzer Lyman J., Stone M. E., 1967, ApJ, 147, 519
  * Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New Astron., 6, 79
  * Springel et al. (2021) Springel V., Pakmor R., Zier O., Reinecke M., 2021, MNRAS, 506, 2871
  * Spurzem & Aarseth (1996) Spurzem R., Aarseth S. J., 1996, MNRAS, 282, 19
  * Strohmayer & Mushotzky (2003) Strohmayer T. E., Mushotzky R. F., 2003, ApJ, 586, L61
  * Sugimoto & Bettwieser (1983) Sugimoto D., Bettwieser E., 1983, MNRAS, 204, 19P
  * Sun et al. (2017) Sun N.-C., et al., 2017, ApJ, 835, 171
  * Takahashi (1995) Takahashi K., 1995, PASJ, 47, 561
  * Torniamenti et al. (2022) Torniamenti S., Pasquato M., Di Cintio P., Ballone A., Iorio G., Artale M. C., Mapelli M., 2022, MNRAS, 510, 2097
  * Trani & Spera (2023) Trani A. A., Spera M., 2023, in Bisikalo D., Wiebe D., Boily C., eds, Vol. 362, The Predictive Power of Computational Astrophysics as a Discover Tool. Cambridge University Press, pp 404–409 (arXiv:2206.10583), doi:10.1017/S1743921322001818
  * Tremou et al. (2018) Tremou E., et al., 2018, ApJ, 862, 16
  * Valtonen & Karttunen (2006) Valtonen M., Karttunen H., 2006, The Three-Body Problem. Cambridge University Press
  * Vanzella et al. (2019) Vanzella E., et al., 2019, MNRAS, 483, 3618
  * Vanzella et al. (2023) Vanzella E., et al., 2023, ApJ, 945, 53
  * Vergara et al. (2023) Vergara M. C., Escala A., Schleicher D. R. G., Reinoso B., 2023, MNRAS, 522, 4224
  * Vink et al. (2015) Vink J. S., et al., 2015, Highlights of Astronomy, 16, 51
  * Wang et al. (2015) Wang L., Spurzem R., Aarseth S., Nitadori K., Berczik P., Kouwenhoven M. B. N., Naab T., 2015, MNRAS, 450, 4070
  * Wang et al. (2016) Wang L., et al., 2016, MNRAS, 458, 1450
  * Wang et al. (2020) Wang L., Nitadori K., Makino J., 2020, MNRAS, 493, 3398
  * Weidner & Kroupa (2006) Weidner C., Kroupa P., 2006, MNRAS, 365, 1333
  * Werchan & Zaritsky (2011) Werchan F., Zaritsky D., 2011, AJ, 142, 48
  * Woosley (2017) Woosley S. E., 2017, ApJ, 836, 244
  * Woosley et al. (2007) Woosley S. E., Blinnikov S., Heger A., 2007, Nature, 450, 390
  * Yan et al. (2023) Yan Z., Jerabkova T., Kroupa P., 2023, A&A, 670, A151
  * Yoshida et al. (2016) Yoshida T., Umeda H., Maeda K., Ishii T., 2016, MNRAS, 457, 351
  * Yu et al. (2011) Yu J., de Grijs R., Chen L., 2011, ApJ, 732, 16
  * Yungelson et al. (2008) Yungelson L. R., van den Heuvel E. P. J., Vink J. S., Portegies Zwart S. F., de Koter A., 2008, A&A, 477, 223
  * Zevin et al. (2020) Zevin M., Kelley L. Z., Nugent A., Fong W.-f., Berry C. P. L., Kalogera V., 2020, ApJ, 904, 190
  * Zhang & Fall (1999) Zhang Q., Fall S. M., 1999, ApJ, 527, L81
  * Zhang et al. (2001) Zhang Q., Fall S. M., Whitmore B. C., 2001, ApJ, 561, 727
  * Zlochower & Lousto (2015) Zlochower Y., Lousto C. O., 2015, Phys. Rev. D, 92, 024022
  * Zlochower et al. (2011) Zlochower Y., Campanelli M., Lousto C. O., 2011, Classical and Quantum Gravity, 28, 114015
  * van der Marel & Anderson (2010) van der Marel R. P., Anderson J., 2010, ApJ, 710, 1063

##  Appendix A The updated BIFROST code: gravitational dynamics

###  A.1 Simplified gradient term calculation for fourth-order forward
symplectic integrator

The fourth-order forward symplectic integrator (FSI, Chin 1997; Chin & Chen
2005; Chin 2007; Dehnen & Read 2011; Dehnen & Hernandez 2017) and its
hierarchical version (HHS-FSI, Rantala et al. 2021) need to calculate the so-
called gradient accelerations $\boldsymbol{g}_{\mathrm{i}}$ for a modified
mid-step kick operation in addition to the common Newtonian accelerations
$\boldsymbol{a}_{\mathrm{i}}=\boldsymbol{a}_{\mathrm{i}}(\boldsymbol{r}_{\mathrm{i}})$
in order to reach the fourth-order accuracy. We do not derive the integrator
here or review its properties, but instead refer the reader to Rantala et al.
(2021) and references therein for the details.

Briefly, the mid-step gradient kick operation requires the calculation of the
modified accelerations $\boldsymbol{\tilde{a}}_{\mathrm{i}}$ defined as

$\begin{split}\boldsymbol{\tilde{a}}_{\mathrm{i}}&=\boldsymbol{a}_{\mathrm{i}}+\frac{\epsilon^{2}}{48}\boldsymbol{g}_{\mathrm{i}}\\\
\boldsymbol{g}_{\mathrm{i}}&=2G\sum_{\mathrm{j\neq
i}}^{\mathrm{N}}\frac{m_{\mathrm{j}}}{r^{5}_{\mathrm{ji}}}\bigg{(}r_{\mathrm{ji}}^{2}\boldsymbol{a}_{\mathrm{ji}}-3(\boldsymbol{x}_{\mathrm{ji}}\cdot\boldsymbol{a}_{\mathrm{ji}})\boldsymbol{x}_{\mathrm{ji}}\bigg{)}\end{split}$
(16)

in which $\epsilon$ is the time-step, $m_{\mathrm{j}}$ the particle masses,
$\boldsymbol{x}_{\mathrm{ij}}$ the relative particle separations and
$\boldsymbol{a}_{\mathrm{ij}}$ the relative Newtonian accelerations. However,
this step can be simplified by using the so-called Omelyan’s approximation
(Omelyan, 2006) as noted by Chin (2007). The essence of the approximation is
the fact that the term $\boldsymbol{\tilde{a}}_{\mathrm{i}}$ can be very well
approximated in Keplerian potentials by ordinary Newtonian accelerations
calculated using displaced positions

$\boldsymbol{r}_{\mathrm{i}}^{\prime}=\boldsymbol{r}_{\mathrm{i}}+\frac{\epsilon^{2}}{24}\boldsymbol{a}_{\mathrm{i}}$
(17)

so that
$\boldsymbol{\tilde{a}}_{\mathrm{i}}(\boldsymbol{r}_{\mathrm{i}})\approx\boldsymbol{a}_{\mathrm{i}}(\boldsymbol{r}_{\mathrm{i}}^{\prime})$.
Thus, Omelyan’s approximation allows for using the same numerical procedures
for the gradient acceleration calculation as for the common Newtonian
calculations. While the computational expense of the approximation is still
comparable to the original method, the algorithmic complexity of the required
code is greatly reduced. We find this especially convenient for the CUDA
kernels of BIFROST.

###  A.2 Slow-down algorithmic regularization for subsystem integration

Instead using the standard logarithmic Hamiltonian (LogH) variant of the
algorithmically regularized few-body integration (e.g. Mikkola & Tanikawa
1999; Preto & Tremaine 1999; Mikkola & Aarseth 2002; Mikkola & Merritt 2006;
Mikkola 2008; Hellström & Mikkola 2010; Rantala et al. 2017, 2020; Trani &
Spera 2023) as in the previous BIFROST version (Rantala et al., 2023), we
implement a version of slow-down algorithmic regularization (SDAR, Wang et al.
2020) into the novel version of BIFROST. The technical details of the
implementation closely follow Wang et al. (2020). We very briefly outline the
main properties and advantages of the method here.

The key idea of slow-down regularization is to identify weakly perturbed
binary systems and slow down their relative orbital motion by hand by a factor
of $\kappa$ so that the new orbital period of the binary is $\kappa P$ instead
of the Keplerian period $P$ (Mikkola & Aarseth, 1996). If there are no weakly
perturbed binaries in the system, the SDAR reduces to the ordinary algorithmic
regularization technique. The increased orbital period obviously makes the
orbital integration of the system containing the binary computationally less
expensive. The slow-down factor $\kappa>1$ can be determined from the strength
of the tidal perturbation onto the binary and the fact that the perturbation
is only allowed to change a little during one slow-down period $\kappa P$.
From the point of view of the binary, the slow-down orbital motion corresponds
to an increased strength of the external perturbation. Formidably, the slow-
down orbital evolution of the binary corresponds to its evolution in the
secular approximation (Wang et al., 2020). This fact with the lowered
computational costs make the SDAR a powerful tool for integrating few-body
systems which include weakly perturbed binaries.

In BIFROST, we replace both the LogH algorithmic regularization and the
secular integration of three-body systems used in Rantala et al. (2023) with
our implementation of the SDAR method. The main difference of our
implementation of SDAR compared to the original SDAR of Wang et al. (2020) is
the use of the Gragg-Bulirsch-Stoer (GBS) extrapolation method (Gragg, 1965;
Bulirsch & Stoer, 1966) to ensure the user-desired numerical accuracy. The
original SDAR implementation instead uses a high-order symplectic integrator.
In this study we do not use the slow-down method (that is, we set $\kappa=1$)
for subsystems requiring post-Newtonian (PN) equations of motion, i.e.
subsystems including black holes. As in Rantala et al. (2023), we use PN terms
up to the order PN3.5.

Following Wang et al. (2020), the equations of motion our SDAR implementation
originate from a slow-down Hamiltonian $H_{\mathrm{sd}}$ defined as

$H_{\mathrm{sd}}=\sum_{\mathrm{i}}^{\mathrm{N_{b}}}\frac{1}{\kappa_{\mathrm{i}}}H_{\mathrm{b,i}}+\left(H-\sum_{\mathrm{i}}^{\mathrm{N_{b}}}H_{\mathrm{b,i}}\right),$
(18)

in which $H=T+U$ is the Hamiltonian of the original system, $N_{\mathrm{b}}$
the number of binary systems, $H_{\mathrm{b,i}}$ their two-body Hamiltonians
and $\kappa_{\mathrm{i}}$ their individual slow-down factors. The slow-down
Hamiltonian can be separated into kinetic and potential terms as
$H_{\mathrm{sd}}=T_{\mathrm{sd}}+U_{\mathrm{sd}}$ with $T_{\mathrm{sd}}$ and
$U_{\mathrm{sd}}$ defined as

$\begin{split}T_{\mathrm{sd}}&=\sum_{\mathrm{i}}^{\mathrm{N_{b}}}\frac{1}{\kappa_{\mathrm{i}}}T_{\mathrm{b,i}}+\sum_{\mathrm{i}}^{\mathrm{N_{b}}}T_{\mathrm{b,cm,i}}+\sum_{\mathrm{i=2N_{\mathrm{b}}+1}}^{\mathrm{N}}T_{\mathrm{i}}+B_{\mathrm{sd}}\\\
U_{\mathrm{sd}}&=\sum_{\mathrm{i}}^{\mathrm{N_{b}}}\frac{1}{\kappa_{\mathrm{i}}}U_{\mathrm{b,i}}+\left[U-\sum_{\mathrm{i}}^{\mathrm{N_{b}}}\frac{1}{\kappa_{\mathrm{i}}}U_{\mathrm{b,i}}\right].\end{split}$
(19)

In the definition of $T_{\mathrm{sd}}$ the first term is the kinetic energy of
the slow-down binaries, while the second and the third terms are the unaltered
kinetic energies of binary center-of-masses and single particles. The indexes
are ordered so that the binaries are first in the arrays, and then the
singles. The last term is the so-called momentum of time term
$p_{\mathrm{t}}=-H_{\mathrm{sd}}(t=0)=B_{\mathrm{sd}}$ in the extended phase-
space of algorithmic regularization, i.e. the initial binding energy. The
potential term $U_{\mathrm{sd}}$ is somewhat more simple, consisting of the
internal potential energy in the slow-down binaries (first term) and the
potential energy of the rest of the system (the term in the brackets).

Now, the time-transformed equations of motion can be constructed. First, the
equation of motion for time $t$ as a function of the new independent variable
$s$ follows directly from the chosen time transformation
$\mathrm{ds}=U_{\mathrm{sd}}\mathrm{dt}$ as

$\frac{\mathrm{d}t}{\mathrm{d}s}=\frac{1}{U_{\mathrm{sd}}}=\frac{1}{T_{\mathrm{sd}}}.$
(20)

Note that here $T_{\mathrm{sd}}$ includes the momentum-of-time binding energy
term $B_{\mathrm{sd}}$, which is typically displayed on its own when
discussing the time transformation (e.g. Rantala et al. 2017, 2020). Next, the
equations of motion can be written for the single particles as

$\begin{split}\frac{\mathrm{d}\boldsymbol{r}_{\mathrm{i}}}{\mathrm{d}s}&=\frac{1}{T_{\mathrm{sd}}}\boldsymbol{v}_{\mathrm{i}}\\\
\frac{\mathrm{d}\boldsymbol{v}_{\mathrm{i}}}{\mathrm{d}s}&=\frac{1}{U_{\mathrm{sd}}}\boldsymbol{a}_{\mathrm{i}}\end{split}$
(21)

and for the particles in binary systems as

$\begin{split}\frac{\mathrm{d}\boldsymbol{r}_{\mathrm{k,i}}}{\mathrm{d}s}&=\frac{1}{T_{\mathrm{sd}}}\left[\frac{1}{\kappa_{\mathrm{i}}}\left(\boldsymbol{v}_{\mathrm{k,i}}-\boldsymbol{v}_{\mathrm{cm,i}}\right)+\boldsymbol{v}_{\mathrm{cm,i}}\right]\\\
\frac{\mathrm{d}\boldsymbol{v}_{\mathrm{k,i}}}{\mathrm{d}s}&=\frac{1}{U_{\mathrm{sd}}}\left[\frac{1}{\kappa_{\mathrm{i}}}\boldsymbol{a}_{\mathrm{bin,i}}+\left(\boldsymbol{a}_{\mathrm{k,i}}-\boldsymbol{a}_{\mathrm{bin,i}}\right)\right].\end{split}$
(22)

Here $\boldsymbol{a}_{\mathrm{bin,i}}$ refer to the acceleration within each
individual binary.

We use the so-called jumping-$\kappa$ method closely following Wang et al.
(2020) to change the slow-down factors of individual binaries. The slow-down
factor is determined using the relative strength at the apocenter approach
already present in a simple form in Mikkola & Aarseth (1996). For binary with
a label $i$ and $N_{\mathrm{p}}$ perturbers the slow-down factor
$\kappa_{\mathrm{i}}$ is

$\kappa_{\mathrm{i}}=k_{\mathrm{ref}}\frac{m_{\mathrm{i,1}}m_{\mathrm{i,2}}}{(m_{\mathrm{i,1}}+m_{\mathrm{i,2}})r_{\mathrm{i,apo}}^{3}}\sum_{\mathrm{j}}^{\mathrm{N_{\mathrm{p}}}}\frac{\|\boldsymbol{r}_{\mathrm{j}}-\boldsymbol{r}_{\mathrm{i,cm}}\|^{3}}{m_{\mathrm{j}}}$
(23)

in which $k_{\mathrm{ref}}=10^{-6}$. The maximum slow-down factor is
calculated by estimating the how rapidly the perturbation itself is changing.
Denoting $\langle X\rangle_{\mathrm{m}}$ as the mass-weighted average of
quantity $X_{\mathrm{j}}$ over the perturbers $j$, we have

$\kappa_{\mathrm{i,max}}=\frac{c}{P}\frac{\|\langle\boldsymbol{r}_{\mathrm{j}}-\boldsymbol{r}_{\mathrm{i,com}}\rangle_{\mathrm{m}}\|}{\|\langle\boldsymbol{v}_{\mathrm{j}}-\boldsymbol{v}_{\mathrm{i,com}}\rangle_{\mathrm{m}}\|},$
(24)

in which $c=0.1$.

We enforce the step-size control and ensure the numerical accuracy of the SDAR
integration by using the Gragg-Bulirsch-Stoer (GBS) extrapolation method
(Gragg, 1965; Bulirsch & Stoer, 1966). We use at maximum $k_{\mathrm{max}}=8$
sub-step divisions, and the user given tolerance parameters are
$\eta_{\mathrm{gbs}}=5\times 10^{-8}$ for the GBS algorithm and $10^{-2}$ for
the end-time iteration procedure. For further details on practical
implementation of algorithmically regularized integrators and the GBS method
in regularized integration see Rantala et al. (2017, 2020, 2023). For further
details of the SDAR method we refer the reader to Wang et al. (2020).

##  Appendix B The updated BIFROST code: stellar evolution, collisions and
mergers

###  B.1 SEVN and this study

Fast stellar population synthesis codes for single and binary star evolution,
pioneered by SSE and BSE by Hurley et al. (2000, 2002), with a straightforward
direct N-body or cluster Monte Carlo interface are essential for realistic
simulations of star cluster evolution. SSE uses computationally efficient
fitting polynomials to the stellar evolution tracks of Pols et al. (1998) to
follow the evolution of stars. Instead of the pre-computed fitting
polynomials, recent fast population synthesis codes, such as combine (Kruckow
et al., 2018), metisse (Agrawal et al., 2020), posydon (Fragos et al., 2023)
and SEVN (Iorio et al., 2023) rely on multi-dimensional interpolation of the
stellar tracks. A major advantage of the multidimensional interpolation
compared to the fitting formulas is the ability to change the stellar tracks
once new ones become available.

SEVN (Iorio et al., 2023; Mapelli et al., 2020; Spera et al., 2015; Spera &
Mapelli, 2017) is a fast stellar population synthesis code, written in C++.
Binary evolution prescriptions are based on analytic and semi-analytic
formulas following (Hurley et al., 2002) plus updates for common envelope
treatments and stellar mergers (Spera et al., 2019; Iorio et al., 2023). The
reason for us choosing SEVN to be coupled to BIFROST over the other options
are the following. First, SEVN is written in C++ making the C/CUDA coupling to
BIFROST more straightforward. Second, SEVN is easily available, and the code
framework enables importing improved stellar tracks in the future. In this
work (FROST-CLUSTERS I), we model systems with no primordial binary stars, so
we focus here in single stellar evolution, and stars in dynamically formed
binaries evolve as they were single stars. The effect of primordial binary
dynamics and binary stellar evolution will be a topic for a forthcoming
follow-up study (FROST-CLUSTERS II).

SEVN can rapidly evolve stars with a wide range of masses and metallicities by
interpolating pre-calculated H-rich and pure-He stellar tracks selected by the
user. At the moment, SEVN has two PARSEC-based (Bressan et al., 2012; Chen et
al., 2015; Costa et al., 2021; Nguyen et al., 2022) track sets and one set of
MIST (Choi et al., 2016) tracks available. As we are interested in the
evolution of massive stars and their collisions, for this work we choose the
PARSEC-based SEVNtracks_parsec_ov04_AGB tracks ranging in stellar mass from
${2.2}\>\mathrm{M_{\odot}}\leq m_{\star}\leq{600}\>\mathrm{M_{\odot}}$ as they
reach higher stellar masses than the other available stellar tracks. The
tracks assume an overshooting parameter of $\lambda_{\mathrm{ov}}=0.4$, and we
use a fixed metallicity of $Z=0.0002=0.01\>Z_{\odot}$ in all of our runs. The
other H-rich PARSEC tracks not used in this work, SEVNtracks_parsec_ov05_AGB
tracks have an overshooting parameter of $\lambda_{\mathrm{ov}}=0.5$, extend
up to $m_{\star}={450}\>\mathrm{M_{\odot}}$ in stellar mass, but only for
specific metallicities the tracks extend up to
$m_{\star}={600}\>\mathrm{M_{\odot}}$. The available MIST (Choi et al., 2016)
SEVNtracks_MIST_AGB only extend up to $m_{\star}={150}\>\mathrm{M_{\odot}}$.

In practice, the main difference between the two PARSEC tracks sets for our
simulations would be the stellar life-times from the beginning of the main
sequence (MS) until the beginning of the remnant phase, and the remnant
masses. The life-times from the beginning of the MS to the beginning of the
remnant phase at the selected metallicity for the $\lambda_{\mathrm{ov}}=0.5$
tracks would be typically $4\%$–$8\%$ shorter compared to the
$\lambda_{\mathrm{ov}}=0.4$ tracks. The use of $\lambda_{\mathrm{ov}}=0.5$
tracks would lead to increased remnant BH masses by up to $\sim 20\%$–$25\%$
for stars with initial masses of
$m_{\star}\sim{15}\>\mathrm{M_{\odot}}$–${20}\>\mathrm{M_{\odot}}$, and
decreased remnant masses near the lower end of the (P)PNS gap by $\sim
17\%$–$23\%$. The (P)PISN mass gap is somewhat narrower for
$\lambda_{\mathrm{ov}}=0.5$ tracks at our metallicity of $Z=0.01\>Z_{\odot}$
compared to the $\lambda_{\mathrm{ov}}=0.4$ tracks.

SEVN stellar process | model | reference
---|---|---
Core-collapse SN | delayed | Fryer et al. (2012)
Neutrino mass loss | lattimer89 | Lattimer & Yahil (1989)
(P)PISN processes | mapelli20 | Woosley (2017)
|  | Spera & Mapelli (2017)
|  | Mapelli et al. (2020)
SN kicks | unified | Giacobbo & Mapelli (2020)
Table 8: SN-related stellar evolution models used in this work.

In addition to the track interpolation, SEVN models a variety of single
stellar evolution processes, which include supernova-related aspects of
stellar evolution not included in the pre-computed evolutionary tracks. The
most relevant such processes to this work are the core collapse supernovae
(SN), SN neutrino mass loss (Zevin et al., 2020; Lattimer & Yahil, 1989),
(P)PISN processes (Woosley et al., 2007; Yoshida et al., 2016; Woosley, 2017)
and SN kicks. The models used in this work are listed in Table 8. For in-depth
descriptions of the individual physical processes, see Iorio et al. (2023) and
references therein.

###  B.2 Interfacing SEVN with BIFROST

####  B.2.1 Code interface

We use SEVN as a library and have written a library interface in BIFROST using
C++ and standard C. The basic building blocks of the interface are the various
wrapper functions to access evolving stellar properties from SEVN in BIFROST,
functions for memory allocation, initialization of stars, evolving the stars,
and finally performing stellar collisions and mergers. As a fast population
synthesis code, SEVN is efficient enough to be used in a serial manner in
BIFROST. As the slow-fast hierarchical Hamiltonian splitting approach for
time-stepping in BIFROST requires regular pivot sorting of the data (Rantala
et al., 2021), we store the SEVN stellar evolution data into a separate
structure in BIFROST, and only sort pointer indexes of the particles to the
stellar evolution data alongside particle masses, positions and velocities.

The stellar evolution is performed after the FSI phase in the BIFROST HHS-FSI.
Every star in the SEVN mass range is not evolved each time they enter the FSI.
Instead, each star is assigned a time stamp for their next evolution time
$t_{\mathrm{next}}$ based on the simulation time and their stellar evolution
time-step $\epsilon_{\mathrm{se}}$ from SEVN. If the simulation time when
entering the FSI exceeds the next evolution time $t_{\mathrm{next}}$, the star
is evolved using SEVN, and $t_{\mathrm{next}}$ is updated. In addition, at the
end of each BIFROST integration interval, each star is evolved up to that
time.

####  B.2.2 Time-step assignment

The BIFROST time-step assignment for simulation particles is modified if
stellar evolution is enabled. For a star, given the dynamical time-step
$\epsilon_{\mathrm{dyn}}=\min{(\epsilon_{\mathrm{ff}},\epsilon_{\mathrm{fb}},\epsilon_{\mathrm{\nabla}})}$
and the stellar evolution time-step $\epsilon_{\mathrm{se}}$ obtained from
SEVN, the time-step is
$\epsilon=\min{(\epsilon_{\mathrm{dyn}},\epsilon_{\mathrm{se}})}$. For
additional details of the adaptive time-step algorithm in SEVN, see Iorio et
al. (2023). In addition, stars very close to the end of their lives
(phase_star=6) may have their time-steps reduced due to the possible SN kick
received at the beginning of the remnant phase phase_star=7 so that the
dynamical interactions with other stars after the SN kick are accurately
treated. Namely, we use
$\epsilon=\min{(\epsilon_{\mathrm{dyn}},\epsilon_{\mathrm{se}},\epsilon_{\mathrm{SNkick}})}$
with
$\epsilon_{\mathrm{SNkick}}=\eta_{\mathrm{fb}}r_{\mathrm{ngb}}/\tilde{v}_{\mathrm{kick}}$
in which $\eta_{\mathrm{fb}}\sim 0.1$ is the typical fly-by time-step factor
of BIFROST, $r_{\mathrm{ngb}}\sim 10^{-3}$ pc the typical spatial size of the
subsystem regions and $\tilde{v}_{\mathrm{kick}}\sim 10^{3}$ km/s a typical
large SN kick velocity. These order-of-magnitude estimates yield a SN kick
time-step of the order of $\epsilon_{\mathrm{SNkick}}\sim 10^{-7}$ Myr.

####  B.2.3 Low-mass stars and very massive stars

The stellar tracks SEVNtracks_parsec_ov04_AGB used in this work range from
${2.2}\>\mathrm{M_{\odot}}\leq m_{*}\leq{600}\>\mathrm{M_{\odot}}$ in initial
stellar mass. The longest simulations times in this work are $50$ Myr, and
stars with masses of $\lesssim{2.2}\>\mathrm{M_{\odot}}$ remain in the main
sequence experience little evolution during this time. Thus, in this work we
do not evolve such low-mass stars and treat them as they were in the beginning
of their main sequence phase. Low-mass stars may become evolving and enter the
SEVN mass range through stellar mergers with other low-mass stars.

As massive stars grow by mergers, they may exceed the maximum mass of
${600}\>\mathrm{M_{\odot}}$. These very massive stars follow simplified
stellar evolution by scaling the $m_{\star}={600}\>\mathrm{M_{\odot}}$
SEVNtracks_parsec_ov04_AGB stellar track. The VMS inherits its life-time and
spin magnitude from the $m_{\star}={600}\>\mathrm{M_{\odot}}$ track while the
wind mass loss and the VMS radius are scaled from the
$m_{\star}={600}\>\mathrm{M_{\odot}}$ using power-law scaling. The mass loss
rate scales proportional to $m_{\mathrm{VMS}}/{600}\>\mathrm{M_{\odot}}$ while
the VMS radius scales as $(m_{\mathrm{VMS}}/{600}\>\mathrm{M_{\odot}})^{1/3}$
to maintain a constant stellar density. We note that while our VMS treatment
is simplified, detailed VMS and SMS models have been developed (e.g. Gieles et
al. 2018). Including such a model in BIFROST remains a topic for future work.

###  B.3 Collisions and mergers

####  B.3.1 The collision detection procedure

We perform the procedures for checking for particle collisions within the
regularized subsystem integration routine. The current BIFROST code version
uses a two-fold procedure for particle collision detection based on the
separation of two particles
$r_{\mathrm{sep}}=\|\boldsymbol{r}_{\mathrm{j}}-\boldsymbol{r}_{\mathrm{i}}\|$
and their osculating pericenter distance $r_{\mathrm{p}}$ calculated from
their osculating two-body orbital elements. A two-body collision criterion
based on osculating two-body orbital elements alone may suffer from a boosted
collision rate due to false positive collisions: while the osculating elements
indicate a collision, it is not guaranteed to occur in a system with more than
two bodies.

In the current code version there are three main types of collisions: BH-BH
mergers, BH-star tidal disruption events, and star-star collisions. Each type
of collision is associated with a characteristic radial separation
$r_{\mathrm{coll}}$ below which a collision certainly occurs, i.e.
$r_{\mathrm{sep}}<r_{\mathrm{coll}}$. We term these true collisions. We also
nevertheless check the pericenter criterion
($r_{\mathrm{p}}<r_{\mathrm{coll}}$) using the osculating two-body orbital
elements. The particle pairs fulfilling the osculating pericenter collision
criterion while not fulfilling the separation criterion
($r_{\mathrm{sep}}>r_{\mathrm{coll}}$) are treated as collision candidates.

If a collision candidate is detected within the SDAR subsystem integration, we
first disable the slow-down treatment, i.e. set $\kappa=1.0$. Next, we
evaluate the time to the next pericenter (bound orbits) or to the only
pericenter passage (hyperbolic orbits) of the collision candidate particle
pair using their osculating two-body orbital elements. Now instead of letting
the GBS routine choose the regularized time-step, we force the regularized
integrator to take (typically less than a few) steps until the actual
pericenter is reached. At the pericenter we re-evaluate the collision criteria
to actually resolve whether there is a collision or not.

####  B.3.2 BH-BH mergers

A BH-BH occurs in the current version of BIFROST based on the radius of the
innermost stable circular orbit of a zero-spin BH defined as
$r_{\mathrm{ISCO}}=6Gm_{\bullet}/c^{2}$. The BH pair collision criterion is
thus $r_{\mathrm{sep}}<\max{(r_{\mathrm{ISCO,p}},r_{\mathrm{ISCO,s}})}$. The
mass, spin and recoil velocity of the remnant BH is calculated using the
fitting formulas to numerical relativity simulations by Zlochower & Lousto
(2015). Examples for recoil kick velocities for various BH binary mass ratios
and spins are elaborated in section 2.4. The gravitational-wave mass loss is
typically of the order of $\sim 5\%$ of the original total mass of the system
(e.g. Jiménez-Forteza et al. 2017).

####  B.3.3 Tidal disruption of stars by BHs

We use a simple order-of-magnitude estimation for the radius below which a
tidal disruption event occurs (TDE) as
$r_{\mathrm{sep}}<r_{\mathrm{TDE}}=1.3(m_{\bullet}/m_{\star})^{1/3}R_{\star}$
based on the standard cube root scaling of Kochanek (1992). Most TDEs in the
simulations in our runs occur from almost parabolic orbits with low
pericenters. We assume an accretion factor of $f_{\mathrm{c}}=0.5$ following
Rizzuto et al. (2023), so half of the stellar material is instantly accreted
into the BH while the remaining half is removed from the simulation. Most
tidal disruption events in this study occur between IMBHs and stars with
masses less than $m_{\star}\lesssim{100}\>\mathrm{M_{\odot}}$, and TDEs are a
subdominant channel in IMBH growth in the setups we have examined. Our simple
approach in this work does not include partial or repeated tidal disruption
events.

####  B.3.4 Stellar collisions and mergers

In the earlier BIFROST version of Rantala et al. (2023) no distinction was
made between a stellar collision (the radii of the two stars overlap) and a
stellar merger (two stars are mixed into one). If two non-remnant stars had
their pericenter distance smaller than the sum of their radii, i.e.
$r_{\mathrm{p}}<r_{\mathrm{\star,p}}+r_{\mathrm{\star,s}}$, the stars were
immediately merged.

In the current BIFROST version, we use SEVN to evaluate the outcome of each
stellar collision ($r_{\mathrm{p}}<r_{\mathrm{\star,p}}+r_{\mathrm{\star,s}}$)
within the SEVN mass range. Instead of merging and complete mixing of the two
stars, a stellar collision can lead to a common envelope (CE) phase from which
the two stars may survive without actually merging. This occurs especially in
the late life stages of stars. Binary systems surviving the CE phase have
reduced orbital separations and low eccentricities. Stars colliding from
hyperbolic orbits are in practice merged from almost parabolic bound orbits as
the current version of SEVN assumes bound binary orbits for colliding stars.

As in SEVN, in this study we assume no rejuvenation and no mass loss in
stellar mergers. In SEVN, when two stars merge, their total, He core and CO
core masses are simply summed. The merger product inherits the evolutionary
phase and percentage of life of the more evolved progenitor star, and the SEVN
interpolation algorithm finds the new post-merger track self-consistently. For
the somewhat more complex post-CE merger procedure, see Iorio et al. (2023).
Low-mass stellar mergers ($m_{\mathrm{p}}<{2.2}\>\mathrm{M_{\odot}}$,
$m_{\mathrm{s}}<{2.2}\>\mathrm{M_{\odot}}$) and any stellar mergers involving
stars above $m_{\star}>{600}\>\mathrm{M_{\odot}}$ are treated without SEVN
simply by summing up the stellar masses and updating the stellar radii.