EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 709 (May 2026) A&A, 709 (2026) A142 Full HTML
Open Access
Issue
A&A
Volume 709, May 2026
Article Number A142
Number of page(s) 20
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202558659
Published online 12 May 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

There are three main formation channels of seed black holes (BHs) for galactic nuclei that are typically considered today (see the reviews by Rees 1984; Greene et al. 2020, and sources therein). (i) Gravitational runaway mergers between stars and compact objects happen throughout cosmic time in dense star clusters. They can be separated into a ‘fast’ and a ‘slow’ regime following Greene et al. (2020). In the fast regime (a couple of million years from star cluster formation and natal gas expulsion), gravitational runaway mergers can happen during early star cluster evolution, when stars evolve and merge either through binary stellar evolution or dynamical collisions (Portegies Zwart & McMillan 2002; Sakurai et al. 2017; Gieles et al. 2018; Reinoso et al. 2018, 2022; Wang et al. 2022; Vergara et al. 2023, 2026; Rantala et al. 2025). The slow regime (around 100 Myr to a gigayear from star cluster formation and natal gas expulsion) is populated by gravitational runaway mergers that occur between BHs. Gravitational runaway mergers of both kinds are postulated to produce intermediate-mass BHs (IMBHs) with masses of order 102−104 M⊙. Both the fast and the slow regimes have been confirmed extensively by simulations of dense star clusters using various methods and both mechanisms are instrumental to growing IMBHs (e.g. Giersz et al. 2015; Rodriguez et al. 2019; Arca Sedda et al. 2019, 2020b, 2024; Di Carlo et al. 2020a,b, 2021; Rizzuto et al. 2021, 2022; Leveque et al. 2022a; Maliszewski et al. 2022; Vergara et al. 2025; Paiella et al. 2025). (ii) Population III (Pop III) stars above the pair-instability mass gap can collapse directly into BH seeds with masses of ~102 M⊙ (e.g. Bromm & Larson 2004; Bromm 2013; Woosley 2017; Haemmerlé et al. 2020; Mestichelli et al. 2024). Within the Pop III cluster, stars can collide and form an even more massive star, collapsing into a heavier BH seed (e.g. Katz et al. 2015; Sakurai et al. 2017; Reinoso et al. 2018, 2022; Vergara et al. 2021; Tanikawa et al. 2022a; Wang et al. 2022). Reinoso et al. (2025) shows that IMBHs with masses of up to 104 M⊙ are possible to form through stellar mergers and gas accretion. (iii) A pristine gas cloud can avoid fragmentation, collapsing directly into a massive star that eventually forms a heavy seed with masses of around 104−106 M⊙ (e.g. Bromm & Loeb 2003; Begelman et al. 2006, 2008; Begelman 2010; Chon & Omukai 2025; van Dokkum et al. 2025). Extremely massive stars with masses of 105−106 M⊙ can explode by the general relativistic instability supernova mechanism (e.g. Shibata & Shapiro 2002; Sakurai et al. 2015; Uchida et al. 2017).

As was noted above, metal-poor Population (Pop) III star clusters at high redshifts (z ≳ 10) with a top-heavy initial mass function (IMF; e.g. Sharda & Krumholz 2022) are prime candidates to produce BH seeds for galactic nuclei and nuclear star clusters (NSCs Neumayer et al. 2020; Greene et al. 2020; Askar et al. 2021, 2022; Schleicher et al. 2022). The IMBHs that form in these clusters can grow even more massive through tidal disruption events over long timescales (Stone et al. 2017; Sakurai et al. 2019; Wang et al. 2025b). The single- and multi-generation mergers from the gravitational runaway process of BHs and other compact objects originating from Pop III star clusters will be relevant gravitational wave (GW) detection events (e.g. Schneider et al. 2000, 2002, 2003; Kinugawa et al. 2014, 2016, 2021a,b; Hartwig et al. 2016; Belczynski et al. 2017; Ng et al. 2022; Mestichelli et al. 2024; Liu et al. 2024a,b; Wu et al. 2025; Wang et al. 2025a). Liu et al. (2025) shows that binary BH mergers in Pop III star clusters can produce GWs, similar to those observed in GW231123. Plunkett et al. (2025) develops observational strategies to constrain the stellar demographics of Pop III stars with next-generation telescopes, especially for the proposed third-generation ground-based GW detectors, the Cosmic Explorer (CE Reitze et al. 2019; Evans et al. 2021), Einstein Telescope (ET Punturo et al. 2010b,a; Sathyaprakash et al. 2012; Maggiore et al. 2020; Liu & Bromm 2020, 2021), and for future space-borne detectors operating in the millihertz and decihertz frequency bands, such as the laser interferometer space antenna (LISA Amaro-Seoane et al. 2013, 2017, 2023), (TianQin Luo et al. 2016; Liu et al. 2020), or the decihertz GW observatory (DECIGO Kawamura et al. 2011, 2021), where IMBHs are expected to be bright GW sources (e.g. Miller & Hamilton 2002; Amaro-Seoane et al. 2007; Amaro-Seoane 2018; Jani et al. 2020; Arca Sedda et al. 2020a, 2021).

Using detailed binary population synthesis (BPS), Tanikawa et al. (2021b) conducted studies on the merger rate density of Pop III binary BHs below, within, and above the pair-instability mass gap from isolated binary stellar evolution. They find that mergers between two low-mass BHs (in their models, low mass means M < 50 M⊙) independent of mass ratio and semi-major axis distributions of the primordial Pop III binaries could be identified from the BHBH mergers observed by the (Advanced) Laser Interferometer Gravitational-Wave Observatory ((a)LIGO Aasi et al. 2015; Abbott et al. 2018, 2019), (Advanced) Virgo Interferometer ((a)Virgo Acernese et al. 2015; Abbott et al. 2018, 2019) and by extension also the Kamioka Gravitational Wave Detector (KAGRA Abbott et al. 2018, 2020; Kagra Collaboration et al. 2019), although the predicted present-day (10 Gyr) merger rates would be comparatively low (~0.1 yr−1 Gpc−3) with regard to the merger rate density inferred by the GWTC-4 catalogue of 15.4−4.5+6.0Mathematical equation: $15.4^{+6.0}_{-4.5}$ yr−1 Gpc−3 (The LIGO Scientific Collaboration 2025). Similarly, mergers between a low-mass BH and a high-mass BH (in their models, high mass means M > 130 M⊙ due to the pair instability mass gap) or mergers between two high-mass BHs will be detectable using the aforementioned currently available GW detectors according to Tanikawa et al. (2021b). However, the authors caution that if two conditions hold simultaneously, namely a wide minimum semi-major axis in primordial Pop III binaries due stellar expansion in the protostellar phases and fast stellar rotation that causes excitation of (quasi-)chemically homogeneous evolution of these stars, during which Pop III stars stay more compact than their non-rotating counterparts due to mixed helium throughout the star (e.g. Maeder 1987; Yoon & Langer 2005; de Mink et al. 2009). Rapidly rotating Pop III stars (Nandal et al. 2024; Tsiatsiou et al. 2024) have been proposed as a possible source of the high nitrogen abundance observed in galaxies at high redshift (e.g. Bouwens et al. 2010; Tacchella et al. 2023; Marques-Chaves et al. 2024; Naidu et al. 2026). Rotation in Pop III stars induces strong internal mixing, which enlarges the convective core and results in brighter stars, improving the chances of detecting these primordial objects in the early Universe with JWST (Hassan et al. 2025).

Direct detections of Pop III stars, their remnants, or their host clusters remain lacking because they lie at very high redshift. Schauer et al. (2022) reported a candidate extremely massive Pop III star at z = 6.2, but this requires confirmation. Vanzella et al. (2020) claimed a Pop III complex at z = 6.629 with the MUSE Deep Lensed Field targeting the Hubble Frontier Field (HFF) galaxy cluster MACS J0416. Using JWST’s NIRCam and lensing, Fujimoto et al. (2025) identified a strong Pop III cluster candidate at z = 6.5. Another candidate is the high-redshift galaxy RX J2129-z8He, which exhibits He II emission at 1640À (Wang et al. 2024). Despite JWST, unambiguous direct observations remain difficult (Rydberg et al. 2013; Katz et al. 2022). de Souza et al. (2013) suggest that several hundred JWST supernova detections could constrain the Pop III IMF. Future wide-field surveys with Euclid (Laureijs et al. 2011; Tanikawa et al. 2022b) and Roman will likely outperform JWST in constraining directcollapse BHs above the pair-instability gap because they cover larger areas (Lazar & Bromm 2022; Vikaeus et al. 2022).

Because direct observations of Pop III stars and their environments are lacking, the parameter space for their formation is poorly constrained (e.g. Klessen 2019; Klessen & Glover 2023). Fraser et al. (2017) inferred a Pop III IMF from 29 extremely metal-poor Pop II stars that are expected to form in environments enriched by Pop III SNe, and found a Pop I/II-like IMF, but the small sample prevents strong conclusions from being drawn. Theory and simulations generally predict top-heavy Pop III IMFs; such IMFs favour XRB populations with brighter and more spatially irregular X-ray emission (Sartorio et al. 2023). Hydrodynamic models are therefore widely used to probe plausible Pop III IMFs and their environmental dependence (e.g. Hirano & Bromm 2017, 2018; Susa 2019; Chon et al. 2021; Sharda et al. 2021; Latif et al. 2022). The binary fraction and initial binary statistics remain highly uncertain (Stacy & Bromm 2013; Liu et al. 2021a) and, being environment-dependent, affect massive-star and compact-object demographics (Hirano et al. 2018; Sugimura et al. 2020; Kroupa 2025).

The observational and theoretical uncertainties are inherited by the initial conditions of Pop III star cluster simulations. One such parameter concerns the degree of initial Pop III star cluster rotation. A set of direct N-body simulations of initially rotating Pop III star clusters were performed by Vergara et al. (2021, N = 103−104 and evolution only up to 2 Myr). They used MiyamotoNagai models with flattening and rotation and found that not only the collision rate increases with increasing bulk rotation, but additionally the number of escapers is reduced the larger the initial rotation is. The low particle numbers, the short simulation time, and other simplifying assumptions, such as neglecting stellar evolution, make it difficult to generalize these results. In general, for simulations in collisional dynamics of rotating star clusters, distribution functions from, for example, Goodman (1983), Longaretti & Lagoute (1996), Einsel & Spurzem (1999), and Varri & Bertin (2012) are typically used. Kamlah et al. (2022b) recently ran simulations of Pop II star clusters with rotating King models from Einsel & Spurzem (1999, using the 2D Fokker-Planck code FOPAX; Einsel & Spurzem 1999; Kim et al. 2002, 2004, 2008). Apart from the formation and dissolution of a rotating bar of BHs that is related to the gravothermal-gravogyro catastrophe (Hachisu 1979, 1982; Akiyama & Sugimoto 1989; Einsel & Spurzem 1999; Hong et al. 2013; Kamlah et al. 2022b), they found a possible dependence of BHBH binary abundances on the initial star cluster bulk rotation; for a fast-rotating model (rotating King model W0 = 6.0, ω0 = 1.2), they found significantly more BHBH binaries during a phase of star cluster evolution when BHBH mergers in the aforementioned ‘slow’ regime are relevant (Greene et al. 2020). Similarly, Webb et al. (2019) found that increased initial star cluster rotation precipitates the formation of circularized BHBH binaries. Since Pop III are postulated to harbour many more (IM)BHs than Pop II star clusters, relative to their size, we tested whether earlier findings apply to Pop III clusters.

Furthermore, the lack of concrete observations of Pop III stars adds uncertainty to their stellar evolution. The most impactful parameter that differentiates stars of the same masses in Pop I and Pop III populations is the metallicity (therefore, Pop III stars are also referred to as extreme metal poor stars e.g. in Tanikawa et al. 2020). The lower the metallicity, the weaker the radiation-driven winds that affect massive stars. Generally, this statement also holds for pulsation-driven winds (Nakauchi et al. 2020). In the regime of Z/Z⊙ < 10−4 both winds become negligibly small. However, it has been suggested that Pop III stars form with very high rotation rates, when the magnetic fields are negligible, and thus there is negligible magnetic braking (Stacy et al. 2011, 2013; Hirano & Bromm 2018). The winds due to the weakened stellar magnetic field of such stars may be very powerful (Liu et al. 2021b; Jeena et al. 2023).

The lower the metal in a Pop III star, the more drastic its internal evolutionary changes, and the lower its opacity (Ekström et al. 2008). First, the lower the metallicity of a star, the more compact it is, because the line-driven radiation pressure diminishes. Furthermore, most massive Pop I/II stars are characterized by a red supergiant (RSG) evolutionary phase with convective envelopes, while most Pop III stars end with a blue supergiant (BSG) phase and radiative envelopes (Tanikawa et al. 2020). These properties affect binary stellar evolution. While Pop I/Pop II stars tend to undergo unstable mass transfer and common envelope evolution (CEE) in the RSG phase, BSGs from Pop III stars undergo stable mass transfer, so that less mass is ejected from the binary system. As a result, BHBH binaries from Pop III stellar populations can be more massive than Pop I/II counterparts, even ignoring wind mass loss, making Pop III star clusters a very attractive target of GW event detection and progenitor studies (e.g. Inayoshi et al. 2017; Kinugawa et al. 2021a,b; Tanikawa et al. 2021b,a, 2022c,a; Liu et al. 2024a). Dense stellar environments substantially increase the fraction of mergers involving BHs with primary masses exceeding the pair-instability mass limit (Mestichelli et al. 2024, 2025). Furthermore, shorter initial orbital periods lead to higher merger rates Santoliquido et al. (2023), which can boost the production of GW sources in the early Universe.

In this paper, we present eight direct N-body simulations of extremely massive, rotating Pop III clusters, run with and without self-consistent general relativistic (GR) merger recoil kicks. The main novelty is the inclusion of initial bulk rotation. We examine how initial bulk rotation, state-of-the-art Pop III stellarevolution models, GR merger recoil kicks (GRKs), primordial binaries, a continuous IMF, and tidal-field mass loss jointly shape cluster dynamics and IMBH formation.

The paper is structured as follows. Section 2 summarizes the methodology, the implementation of the GRKs, and the Pop III stellar evolution fitting formulae. Section 3 outlines the initial conditions. Section 4 presents the results. In Sect. 5 we summarize and discuss the work.

2 Methods

We used NBODY6++GPU (Spurzem 1999; Nitadori & Aarseth 2012; Wang et al. 2015) for the dynamical evolution of the Pop III clusters. The initial rotating conditions followed Kamlah et al. (2022b) and were generated with MCLUSTER (Kuepper et al. 2011; Leveque et al. 2022b) and FOPAX (Einsel & Spurzem 1999); details are given in Appendix C.2.

Population III stellar evolution was treated with the low-metallicity recipes of Tanikawa et al. (2020) at Z/Z⊙ = 10−8, without radiation-, pulsation- (e.g. Nakauchi et al. 2020), or rotation-driven (e.g. Liu et al. 2021b) wind mass loss; further details are in Appendix C.3. Fitting formulae for relativistic merger recoil kick velocity and remnant mass and spin followed Campanelli et al. (2007); Lousto et al. (2012); Jiménez-Forteza et al. (2017), as summarized in Appendix C.4. We refer to BHs with masses above 100 M⊙ as IMBHs.

3 Initial conditions

The mass function of Pop III star clusters is highly uncertain. Formation of Pop III star clusters has been simulated in cosmological mini-halos (cf. e.g. Stacy et al. 2016), and corresponding star-cluster simulations have covered masses from 103 to 105 M⊙ (Liu et al. 2025; Wu et al. 2025). However, Kimm et al. (2016) propose the formation of more massive Pop III clusters without dark matter directly from cooling gas clouds, with masses of up to 5 ∙ 105M⊙, and Liu et al. (2021a) find cluster masses reaching this value. Stacy et al. (2016) also discuss a possible Pop III cluster with a mass as high as 107 M⊙. In our simulations we explore the very high-mass end (over 105 objects, with 1% binaries and a top-heavy IMF; see Table 1) for two reasons: first, to extend previously published simulations to higher masses; and second, to probe very high initial densities that may promote rapid formation of massive stars and massive BHs (Vergara et al. 2025).

We started from an unsegregated, unfractal King model (King 1966; Goodwin & Whitworth 2004) with W0 = 6.0, initially in virial equilibrium. The half-mass radius was set to rh=1 pc. As is outlined in Appendix C.2, the initial model from MCLUSTER was then redistributed with a rotating King model, which is more compact than its non-rotating counterparts (Einsel & Spurzem 1999). Therefore, the internal structural parameters such as the rh and rc changed in this initialization step from their original MCLUSTER N-body distribution (see Fig. 3).

We used a flat IMF of between 8.0 and 300.0 M⊙ (Lazar & Bromm 2022), which gave an initial crossing time of 0.0074 Myr and half-mass relaxation time of 6.869 Myr. Primordial binaries were paired with uniform mass ratios (0.1 < q < 1.0)1 following previous studies (Kiminki et al. 2012; Sana & Evans 2011; Sana et al. 2012, 2013; Kobulnicky et al. 2014); their semimajor axes were uniformly distributed between 10 and 100 au and their eccentricities were thermal. The cluster metallicity was Z/Z⊙ = 10−8, following Wang et al. (2022) and matching the lowest currently available metallicity for the fitting formulae of Tanikawa et al. (2020).

We adopted a primordial binary fraction of fb = 1%. A higher fb at N ≃ 105 is computationally expensive. The Pop III binary fraction remains highly uncertain and likely depends on formation environment and feedback; hydrodynamic simulations can yield high multiplicities (tens of percent; e.g. Stacy & Bromm 2013), while other physical effects may reduce fragmentation, and thus the multiplicity (e.g. magnetic field; Sharda & Menon 2025; UV, Hosokawa et al. 2016). We therefore treated fb = 1% as a limiting assumption that isolates the role of rotation in a controlled pilot suite and provides a conservative lower bound on contributions from primordial binary evolution channels (Type-I and Type-III; Sect. 4.1). We discuss the expected impact of higher fb on IMBH formation in Sect. 5.

We placed the clusters on a circular orbit at 13.3 kpc following Kamlah et al. (2022b); this reproduced the mass-loss evolution of the eccentric NGC3201 orbit (8.60-29.25 kpc, e = 0.55, Gaia Collaboration 2018), assuming a point-mass MW of 1.78 × 1011 M⊙ and vc = 240.0 km s−1 (Gaia Collaboration 2018; Bobylev & Bajkova 2020), yielding an initial tidal radius of 264.50 pc. This set-up followed Kamlah et al. (2022a,b) to facilitate a direct comparison. our models adopted a static, Milky-Way-like tidal field as a controlled test; we emphasize that realistic high-z environments can produce time-dependent tidal tensors and repeated shocks - effects that cosmological studies show can materially increase mass loss and modify angular-momentum evolution (e.g. Renaud & Gieles 2015). A complementary study embedding our clusters in cosmological tidal histories is warranted.

We stopped the simulations at 500 Myr because (i) the gravothermal-gravogyro catastrophe has largely subsided and the cluster has undergone core-collapse; (ii) most ZAMS stars, particularly the most massive ones, have reached their final stages; and (iii) environmental effects (e.g. cluster mergers) make long-term isolated evolution of Pop III clusters unlikely (see e.g. Arca-Sedda & Gualandris 2018). Thus runs much longer than ~10 Gyr (e.g. Wang et al. 2022) are not always required.

Particles were removed once they reached twice the current tidal radius from the density centre; they were then counted as ‘escapers’. The tidal radius was recomputed from the current cluster mass. Apart from the common set-up above, we ran four initial rotation parameters, ω0, with and without GRKs, for eight simulations in total (Table 2).

Table 1

Initial parameters that are identical across all eight initial models for the NBODY6++GPU simulations.

Table 2

Model identifiers, rotation parameters, and whether simulations enable GR recoil kicks at relativistic mergers.

4 Results

4.1 Stars and compact objects

Figure A.1 tracks counts of different stellar types. Counts are similar during the first million years and then diverge as stellar evolution and core collapse proceed2. Larger ω0 speeds stellar loss, most clearly in (NO)K1.8, while same-ω0 K/NOK pairs remain nearly identical. The trend reflects the gravothermal-gravogyro catastrophe combined with tidal mass loss. By ~20 Myr all MS stars have evolved and none remain. Near 20 Myr all models produce > 1000 NSs (nns), which later decline mainly because of escapers (Sect. 4.5). The BHs appear earlier (~6 Myr) and nbh peaks at ~50 000 before 10 Myr.

We divided the ZAMS particles into four classes: low-mass core-collapse supernovae (lCCSNe): 8.0 ≤ mzams < 85.0 M⊙, pulsational pair-instability supernovae (PPISNe): 85.0 ≤ mzams < 120.0 M⊙, pair-instability supernovae (PISNe): 120.0 ≤ mzams < 265.0 M⊙, and high-mass core-collapse supernovae (hCCSNe): 265.0 ≤ mzams < 300.0 M⊙, where mzams is the ZAMS stellar mass of a single star3. Figures 1 and A.3 show the initial-final-mass relation (IFMR) for BHs retained in the cluster at 500 Myr, with mBH (M⊙) versus the progenitor ZAMS mass, mzams (M⊙). Progenitors produce BHs via lCCSNe (Fryer et al. 2012), PPISNe (Belczynski et al. 2016), or hCCSNe (Fryer et al. 2012); PISNe leave no remnant (Belczynski et al. 2016)4. Star-star and star-BH mergers produce IMBHs inside the gap; IMBHs are found above and below it.

We distinguish remnant-mass regimes of (IM)BHs produced by single versus binary evolution. Figure 1 shows these regimes for the NoK0.0 simulation5.

  • Type-I: the BHs from primordial binaries where H-envelopes of one or both members are fully stripped (by tides, RLOF, or CEE; see Appendix C.3), leaving naked HeMS cores that collapse into BHs. No mergers occur in this channel (c.f. Type-III); survivors may form a BH-BH binary depending on fallback and natal kicks.

  • Type-II: higher-mass BHs from single-star evolution using the delayed core-collapse SN model of Fryer et al. (2012).

  • Type-III: the BHs from primordial-binary mergers (CHeB+MS or CHeB+CHeB during RLOF); the merger product becomes a ShHeB star that later collapses to a BH.

  • Type-IV: dynamically formed IMBHs that always involve RLOF between an IMBH and a ShHeB star.

  • Type-V: the IMBHs from the PPISN channel (see Fig. A.3; e.g. K0.0).

  • Type-VI: two main dynamical channels for IMBH formation: (1) hyperbolic MS+MS collisions producing products that evolve to CHeB/ShHeB and then collapse to BHs; (2) hyperbolic coalescence of a CHeB+BH dynamical binary.

  • Type-VII: hyperbolic CE events in dynamical CHeB/ShHeB binaries whose H-envelopes are stripped off and expelled. The cores merge and the product is a CHeB star following Hurley et al. (2002). The product star then evolves into a ShHeB star and core-collapses into a BH.

  • Type-VIII: IMBHs form by single stellar evolution through core-collapse or the extension of the remnant mass functions of Fryer et al. (2012).

  • Type-IX: dynamically MS/CHeB binaries undergo hyperbolic collision or coalescence following a CE event. Technically in NBODY6++GPU, the product is labelled by the more massive member.

  • Type-X: rare MS+CHeB hyperbolic collisions with significant H-envelope loss, which yield CHeB, evolve into a ShHeB, and core-collapses into an IMBH.

We identified the formation channels present in the simulations using the Pop III fitting formulae of Tanikawa et al. (2020) and our initial conditions. Full N-body dynamics (instead of population synthesis) is crucial: many (IM)BHs result from interactions in dynamically formed binaries.

Within ≤20 Myr and during the pre-core collapse6, the Pop III clusters’ stellar population is radically altered. Compared to similar Pop II clusters (Kamlah et al. 2022b), this evolution is more extreme due to the extreme initial conditions. All Pop III cluster models produce (IM)BHs below, within, and above the pair-instability gap.

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Initial-final-mass relation (IFMR) for the NoK0.0 run: BH mass, mBH, versus progenitor ZAMS mass, mZAMS. The panels correspond to progenitors exploding as lCCSNe, PPISNe, PISNe, and hCCSNe. (IM)BH formation channels are shown in red.

4.2 Compact binary fractions and general properties

The temporal evolution of compact-binary abundances (Fig. A.2; Table B.1) traces the number of dynamical interactions. The (NO)K1.8 models have the lowest binary abundances. The (NO)K1.2 models show a comparable tendency, but their abundances lie close to those of the (NO)K0.0 and (NO)K0.6 models, so the trend is less clear. Notably, most abundances in the (NO)K0.6 models exceed those in the (NO)K0.0 models. See Sect. 4.5 for a detailed analysis based on escaper statistics.

Figure 2 summarizes binaries within the tidal radius at 500 Myr. The e2 distributions are bottom-heavy and similar across simulations. The cumulative a distributions diverge around a ≈ 50 au: most binaries are within a ≈ 50 au in NOK1.8 but not in NOK0.6 or NOK1.2, indicating stronger disruption of wide binaries in the fastest-rotating model mostly before core-collapse (see Fig. A.2).

The Epot distributions show fewer hard binaries in NOK1.8. This can be explained by (i) prior collisions/coalescences (Sect. 4.4, Figs. A.2 and A.5, Table B.1) and/or (ii) preferential escape of hard binaries in rapidly rotating runs (Fig. A.6). This effect ejects many would-be merging binaries early, lowering in situ merger rates of stars and BHs.

The 4th panel of Fig. 2 shows the radial distribution of binaries at 500 Myr (rdens); the vertical line marks rh of the NOK0.0 model. Binary heating dominates post core-collapse evolution and produces the self-similar rLagr evolution seen in Fig. A.4. Although higher initial rotation causes earlier core collapse via the gravothermal-gravogyro instability, by 500 Myr the binary radial distributions, if normalized to the retained binary counts, are similar across models, especially within rh. At 500 Myr, we therefore find no excess of harder binaries near the centres of faster-rotating clusters; long-term evolution largely erases the early structural (Figs. 3 and A.4) and angular-momentum differences, and equalizes the distributions of binaries that remain bound at 500 Myr.

The bottom panel of Fig. 2 shows the mass-ratio distribution q ≡ m2/m1 (m2 < m1). Its shape is insensitive to initial rotation and to whether GR recoil kicks are included.

In the (NO)K0.6 runs GRK reduces binaries by ~50 (mainly at a = 40-50 au, rdens = 10 pc, Epot = −1000 N-body units), while in (NO)K1.2 it increases binaries by ~30 (mainly at a = 50-60 au, rdens = 30 pc). The effects of GRKs in the nonrotating and fastest-rotating models are less visible. Overall, our results do not reveal a clear or systematic effect of GRKs on the binary distribution within the simulated clusters.

Note that the results are affected by small-number statistics, because each simulation contains only 1 000 primordial binaries. More simulations are needed to reach robust conclusions.

Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Binary distributions at 500 Myr. From top to bottom: e2; cumulative semi-major axis, a (au); cumulative binary potential energy, Epot (N-body units); cumulative distance to the density centre, rdens (parsecs); and mass ratio q ≡ m2/m1. In the rdens panel a, blue line shows the half-mass radius rh of the NOK0.0 model at 500 Myr (∼9.77 pc).

4.3 Global dynamical evolution

4.3.1 Structural parameter evolution

We ran each of the four initial models (ω0 = 0.0, 0.6, 1.2, 1.8) using NBODY6++GPU with GRKs switched on (K models) and without (NOK models). Models with the same ω0 evolve similarly whether or not GRKs are enabled, so the dominant control parameter of the global evolution is the initial bulk rotation. Figures 3 and A.4 show a four-phase evolution:

  • Phase I: before the stellar evolution mass loss begins to dominate the star cluster evolution, namely, 0-2 Myr. Mcl and rtid stay constant, because mass loss through escapers or stellar evolution is tiny. Depending on ω0, the other structural parameters in Fig. 3 show many more variations. Kamlah et al. (2022b) used the same rotating King model set-ups for Pop II clusters, and showed that angular-momentum transport is very rapid and usually occurs before core collapse; the larger ω0, the faster this phenomenon happens. In the Pop III clusters in this work, we also observe this fast angular momentum transport, for example in the panel showing the half-mass radius of Fig. 3. Within the first million years, the larger ω0, the more pronounced and earlier is the temporary increase in rh and decrease shortly after. This evolution is even more apparent in Fig. A.4, which shows the evolution of the Lagrangian radii rLagr within spheres that contain 1%, 5%, 10%, 30%, 50%, and 90% of the total cluster mass at the current simulation time step. For example, in the (NO)K1.8 models , the degree of initial bulk rotation has a pronounced effect: the outermost rLagr expands substantially, while the other Lagrangian radii generally contract up to 2 Myr. All rLagr show a small and synchronous bump (i.e. a brief increase followed by a decrease) at around 0.1 Myr. The (NO)K1.2 models show the same effect, but weaker and slightly delayed relative to the (NO)K1.8 models. Such a phenomenon is not shown in (NO)K0.0 and (NO)K0.6 models, where rh stays roughly constant until 2 Myr. This early evolution is also mirrored by rc and mc, where both the (NO)K1.2 and the (NO)K1.8 exhibit fast increases with subsequent decreases in the respective models. Again, this evolution is more pronounced and happens earlier in the (NO)K1.8 with respect to the (NO)k1.2 models. In summary, from 0 to 2 Myr, the evolution is governed only by the initial bulk rotation; tidal and stellar-evolution mass losses are negligible.

  • Phase II: rapid stellar evolution mass loss around 2-3 Myr, as shown in 1st panel of Fig. 3. This is mainly due to compact object formation, particularly PISN events, which produce no BHs and lose all the mass. Besides stellar evolution, mass loss is also caused by (a) escaping stars, through compact object natal kicks or GRKs, or (b) dynamical encounters, including strong encounters or by cumulative weak encounters. The cluster kinematics also change here owing to fallback-dependent natal kicks, which are generally larger for lower-mass BHs (see Kamlah et al. 2022a). As a result of the total mass drop, rtid drops, and rh with rc rises sharply in all models. Interestingly, the maxima of rc lie in the (NO)K1.2 model, while the cores of the (NO)K1.8 and (NO)K0.0 models are similarly compact. We compare this with the average mass in the core Mav within different Lagrangian radii rLagr, shown in Fig. A.4. Signified by Mav, we see that here all models show a similar segregation of masses into different spheres independent of initial bulk rotation; for example, the minimum Mav in the 10% rLagr has similar values occurring at similar times, although we note that this quantity is subject to large statistical fluctuations.

  • Phase III: from the end of phase II (−3 Myr) to the time of gravothermal core-collapse. As shown in Fig. A.4, the ending time of this phase varies between models, which is ~80 Myr for (NO)K0.0, ~70 Myr for (NO)K0.6, ~60 Myr for (NO)K1.2, and −50 Myr for (NO)K1.8. Therefore, the larger the degree of initial bulk rotation in our Pop III clusters, the higher the central density, and the shorter the gravothermal collapse timescale, as seen in previous studies (e.g. Akiyama & Sugimoto 1989; Einsel & Spurzem 1999; Kim et al. 2002, 2004, 2008; Fiestas et al. 2012; Kamlah et al. 2022b). Larger ω0 leads to higher Mav in the 5% and 10% spheres at core collapse. Thus, stronger initial rotation accelerates mass segregation. In (NO)K1.8 models, the cluster expansion ejects many stars7, producing a decrease in the 90% rLagr.

  • Phase IV: from core-collapse to 500 Myr (the simulations’ end). The Lagrangian radii evolve self-similarly, as shown in A.4, driven by binary heating in the core which causes global expansion. Models with larger ω0 expand more by 500 Myr because they undergo earlier core collapse; their subsequent rLagr evolution is otherwise similar. The clusters become fully mass-segregated (Mav ≈ 200 M⊙ centrally, ≈ 30 M⊙ in the halo). The same long-term trends appear in Fig. 3; longer runs may amplify differences. These results are consistent with energy-flow theorems across the half-mass radius rh (Hénon 1975; Breen & Heggie 2013).

The gravothermal-gravogyro catastrophe is therefore present in all our Pop III cluster models, but its timescale depends strongly on the initial rotation. The GRKs have little effect on the global structural evolution.

Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Total cluster mass, Mcl, tidal radius, rtid, half mass radius, rh, mass of the core, mc, and radius of the core, rc, in the four panels for all eight simulations with and without GR recoil kicks for different ω0. Time is shown on a logarithmic scale to highlight the cluster’s rapid early evolution. The models with K are plotted as solid curves and the models without (NOK models) are plotted as dash-dotted curves.

4.3.2 Angular momentum transport

We explored whether and how angular-momentum redistribution depends on stellar evolution, GRKs, and initial bulk rotation. We tracked particles from each group over the full evolution and computed their total angular momentum squared Lgroup2Mathematical equation: $L_{\mathrm{group}}^2$ as Kamlah et al. (2022b) normalized by Lω0=0.6,t=02Mathematical equation: $L^2_{\omega_{0}=0.6,t=0}$ for the (NO)K0.6 model. Figure A.7 shows these quantities together with the group masses, Mgroup.

Phase I: 0 to −2 Myr, before stellar-evolution mass loss dominates. Mgroup evolves similarly across groups, but Lgroup2/Lω0=0.6,t=02Mathematical equation: $L_{\mathrm{group}}^2/L^2_{\omega_{0}=0.6,t=0}$ does not. Higher ω0, notably (NO)K1.8, clearly shows larger Lgroup2Mathematical equation: $L_{\mathrm{group}}^2$ for all mass groups, reflecting rapid global expansion and outward migration of stars. This matches the outer-radius growth and inner-radius contraction visible in Figs. 3 and A.4 (Sect. 4.3.1).

Phase II: −2-3 Myr. Rapid stellar evolution drives a sharp fall in Mgroup, with higher-mass groups declining earlier owing to their shorter evolution timescales. The PISNe remove both mass and a large fraction of angular momentum: MPISNe carried the most angular momentum before these events and is effectively depleted afterwards, causing a large loss of rotational kinetic energy.

Phase III: −3 Myr to core collapse (around −50-80 Myr depending on ω0. See Sect. 4.3.1 and Fig. A.4). The quantity LhCCSNe2/Lω0=0.6,t=02Mathematical equation: $L_{\mathrm{hCCSNe}}^2/L^2_{\omega_{0}=0.6,t=0}$ initially falls after phase II but, in the (NO)K1.2 and (NO)K1.8 models, it instead rises sharply and peaks at −10-20 Myr. After the peak, the quantity of (NO)K1.2 shows a gradual decline, while that of (NO)K1.8 drops rapidly. Other mass groups follow the same pattern, except for MPISNe. In (NO)K1.8 the quantity attains a local minimum at core collapse, when the cluster is essentially fully mass-segregated (Fig. A.4).

Phase IV: core collapse to simulation end. In all models LhCCSNe2/Lω0=0.6,t=02Mathematical equation: $L_{\mathrm{hCCSNe}}^2/L^2_{\omega_{0}=0.6,t=0}$ decreases while LlCCSNe2/Lω0=0.6,t=02Mathematical equation: $L_{\mathrm{lCCSNe}}^2/L^2_{\omega_{0}=0.6,t=0}$ increases, indicating angular momentum transfer from high-mass to low-mass BHs. Low-mass BHs move to wider orbits and a compact IMBH subsystem forms at the centre (also see Fig. A.4). Combined with Sect. 4.3.1 and Fig. 1, this points to binary heating of the central IMBHs as the dominant driver of the cluster evolution in this phase.

Compared to rotating Pop II clusters (Kamlah et al. 2022b), the detailed transport history differs because of the IMF and stellar evolution, but the gravothermal-gravogyro catastrophe remains visible in both cases. The GRKs again produce no strong effect.

4.4 Collision and coalescence events

We defined a ‘coalescence’ as the merger of a nearly circular binary, while ‘collisions’ are a direct physical impact either on a hyperbolic or highly eccentric orbit; this is also how the NBODY6++GPU code separates them. Under our initial conditions and Pop III stellar evolution prescriptions, we find collisions of BHbHs, MSBHs, and MSmSs, while coalescences may involve MS, HG, CHeB, ShHeB, HeMS, NS, and BH progenitors. Figure A.8 shows the product masses, mcoal and mcoll, versus rdens. For coalescence, only MS, CHeB, ShHeB, and BH products appear, consistent with the merger-product treatment of Hurley et al. (2002) and Kamlah et al. (2022a). Event totals are given in Table B.2.

Our suite contains one realization per set-up (N = 1 per model ID). For integer event counts, N (specific coalescence and collision types, compact-object merger channels), the counting uncertainty is well approximated by Poisson statistics, with a reference 1σ uncertainty of σN≃NMathematical equation: $\sigma_N\simeq\sqrt{N}$ (relative uncertainty σN/N≃1/NMathematical equation: $\sigma_N/N\simeq 1/\sqrt{N}$). Consequently, differences of a few events are not, by themselves, strong evidence of a systematic dependence on ω0 when N is in the single-digit to few-tens regime8.

The coalescence events show no marked spatial differences: most occur within the central 2.5 pc, i.e. inside rh during the first ~40 Myr (Fig. 3). They are mostly completed by 4-5 Myr (precore-collapse; Fig. A.5). Many of these early mergers populate the PISN mass gap, indicating that stellar mergers in early star cluster evolution can seed BHs in and around the gap (see e.g. Costa et al. 2022; Ballone et al. 2022, for recent hydrodynamical simulations of stellar collisions on this topic). With our flat IMF, Pop III stellar evolution, and absent wind mass loss, producing BHs above the PISN gap is easier here than in typical young massive clusters (e.g. Rizzuto et al. 2022).

Most coalescences produce CHeB stars and the most massive products, typically occurring within a few million years in all runs. Faster initial rotation is consistent with more CHeB coalescences (Fig. A.5), but Poisson uncertainties and covariances between channels in a single realization make this indicative. The numbers of BH, ShHeB, and MS coalescences are similar across simulations. These mergers produce stars (mainly CHeB) with masses inside and above the pair-instability gap; some IMBHs end up in the pair-instability gap after accreting stellar material and grow further.

Collisions of MSMS, BHMS and BHBH offer useful diagnostics. Collision-produced IMBH masses are broadly similar across models (Fig. A.8), except for NOK1.2 (IMBH masses ≈ 1000 M⊙), which is not seen in K1.2. Figure A.5 suggests more MSMS and BHBH collisions at a higher initial rotation, but given the small absolute counts the implied trend should be read as suggestive and may be affected by stochastic fluctuations. Notably, K models show multi-generation IMBH growth.

Figure A.9 shows the primary versus secondary masses for MSMS, BHMS, and BHBH collisions across all runs; the distributions are similar regardless of rotation or GRKs. Contrary to theoretical expectations (e.g. Morawski et al. 2018, 2019)9, Fig. A.10 demonstrates that most GRK velocities, vGRk, exceed cluster escape speeds in the core, so merger remnants are often ejected even for q ≈ 1. Hence, IMBHs of > 1000M⊙ are unlikely to form via hierarchical merging in Pop III clusters unless additional processes retain the remnants10, such as strong gas inflow after cluster formation causing extreme contraction of the cluster (Kroupa et al. 2020).

To summarize, with increasing initial rotation, ω0, coalescences occur earlier and more frequently, and there is a tentative rise in MSMS and BHBH collision counts limited by low-number statistics. Enabling GRKs generates kick velocities that frequently exceed the cluster escape speed, favouring remnant ejection and thus suppressing collision-built IMBHs, even though multi-generation IMBH growth can still occur.

4.5 Escaper stars

Figure 4 (single escapers), Fig. A.6 (binary escapers), and Table B.3 show that a higher initial bulk rotation, ω0, is associated with earlier and stronger tidal mass loss in both NOK and K models, most clearly for (NO)K1.8 and (NO)K1.2. The (NO)K0.6 models show a non-monotonic behaviour, having nesc,singles, nesc,binarymembers and Mesc,singles comparable to or below the non-rotating case. A plausible physical interpretation is that two effects compete: (i) rotating King initial conditions are more compact than for their non-rotating counterparts (Sect. 1), which can reduce early tidal overflow, while (ii) stronger rotation accelerates early angular-momentum transport and structural readjustment (Sect. 4.3.1), which enhances mass loss once the system expands and interacts with the tidal boundary. In this framework, ω0 = 0.6 can fall in an intermediate regime where the increased initial compactness dominates the net mass loss, whereas for larger ω0 the rotation-driven early evolution dominates. This tidal field mass loss modifies cluster timescales and links directly to the structural and angular-momentum evolution discussed in Sects. 4.3.1 and 4.3.2.

We then focus on the escaper counts by stellar types. The third panel of Fig. 4 shows single MS escapers given by nescMS,s. The (NO)K1.8 runs eject hundreds more MSs than the (NO)K1.2 runs. Many of these MSs will form BHs via lCCSNe in the field (Sect. 4.1), and thus cannot fuel in-cluster hierarchical mergers. Although strong initial rotation speeds up the early collapse and gravogyro contraction (Fig. A.4; Sect. 4.3), the simultaneous11 loss of BH progenitors can hinder subsequent merger chains.

We next considered evolved phases, starting with escaping single CHeB stars (nescCHeB,s; Fig. 4 panel four). Their behaviour resembles that of escaping MSs, though differences between (NO)K1.2 and (NO)K1.8 are reduced. Many escaping MSs rapidly become CHeBs, and likewise many CHeBs shortly become ShHeBs (see Fig. A.1). After formation, most will remain in the ShHeB phase for a comparatively long time.

Escaping ShHeBs (nescShHeB,s) follow the overall ω0 trend, but show a pronounced increase in Phase IV (binary-heating dominated; Sect. 4.3), echoing the earlier rise seen in nescMS,s. Table B.3 indicates that the (NO)K0.6 runs yield the fewest single-star escapers across almost all types. For the rarest subtypes (e.g. nescCHeB,s ≲ O(1-10)), Poisson counting noise can be large and individual model-to-model differences should be interpreted cautiously.

For compact object escapers, nescNS,s shows no strong model-to-model differences, though (NO)K0.6 consistently yields the fewest single NS escapers. Escaping single BHs, nescBH,s, follow a similar behaviour: (NO)K1.8 produces the largest number at 500 Myr (Table B.3), while (NO)K0.0 exhibits the strongest late rise in BH losses (Fig. 4). Because our kicks are fallbackdependent (Appendix C.4), and given the flat IMF (8-300 M⊙) plus absent wind mass loss for BH progenitors, most BHs acquire such large fallback fractions that their natal kick velocities are ≈ 0. At 500 Myr, the (NO)K1.8 runs show the largest losses in compact objects; the (NO)K0.0 runs display accelerating latetime loss that may overtake them, while the (NO)K0.6 runs lose the least mass and the fewest compact objects.

About escaping binaries, Fig. A.6 and Table B.1 show that nesc,binarymembers depends on initial rotation. The (NO)K1.2 and (NO)K1.8 models eject binary members earlier and, by 500 Myr, have roughly twice as many escaped binary members as the (NO)K0.0 and (NO)K0.6 runs. The latter two remain similar and do not display the large long-term differences seen for single escapers (Fig. 4).

For escaping binary components excluding compact objects, a few MS escapers (counted by nescMS,b) occur only in the (NO)K1.2 and (NO)K1.8 runs; they are mostly from primordial binaries, low-mass, and unlikely to form IMBHs. Binary CHeB escapers (nescBH,b) are rarer still. Binary ShHeB losses (nescShHeB,b) are roughly twice as large and occur earlier in the ω0 = 1.2,1.8 models than in the ω0 = 0.0,0.6 models. Many late-time escaping BHBH binaries contain two IMBHs, and thus cannot fuel hierarchical in-cluster growth. The stronger-rotation models are therefore consistent with a reduced retained BH progenitor and seed population at late times.

Overall, enabling GRKs produces the charted differences: with GRKs enabled the (NO)K1.8 models show more CHeB single-star losses and fewer MS binary losses, while the (NO)K1.2 models show the opposite trend. We again note that this is limited by small-number statistics.

Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Escaping single stars in eight simulations: cumulative counts and masses. Top panel: cumulative number of escaping single stars, nesc,singles. Second panel: total mass in escaping single stars, Mesc,singles. Remaining panels: cumulative counts of escaping single MS, CHeB, ShHeB, NS, and BH stars, nescMS,s, nescCHeB,s, nescShHeB,s, nescNS,s, and nescBH,s. Solid curves show K models and dash-dotted curves show NOK models.

5 Summary and conclusions

We present a set of a few pilot direct N-body simulations that probe how bulk rotation and ultra-low metallicity with a top-heavy IMF shape Pop III cluster dynamics. Rotation follows generalized rotating King models (Einsel & Spurzem 1999) with rigid core rotation and differentially decreasing near the half-mass radius and into the halo. Kamlah et al. (2022b) analysed the general evolution and angular momentum transport in such models of Pop II star clusters; we find generally analogous dynamics, such as the gravogyro and gravothermal collapse, enhanced by the top-heavy IMF (8.0-300.0 M⊙, Lazar & Bromm 2022). We study the impact of switching on and off GRKs (Arca Sedda et al. 2024), and how the newly implemented stellar evolution for extremely metal-poor (EMP) stars (here used Z/Z⊙ = 10−8, Tanikawa et al. 2020) affects the formation and evolution of compact objects. While the mass function of Pop III star clusters is highly uncertain (see Sect. 1 and Stacy et al. 2016; Liu et al. 2021a), our pilot simulations model very massive Pop III star clusters (Table 1) to cover an unexplored parameter range and also because of stronger and faster dynamical effects, such as massive object formation, at high density (Vergara et al. 2025). Direct N-body simulations using the PETAR code (Wang et al. 2020) for Pop III star clusters (Liu et al. 2024a,b, 2025) have used smaller cluster masses and particle numbers, requiring many independent realizations to average N fluctuations. On the contrary, we used very large particle numbers (1.01 × 105 initially), which suffice to probe global dynamical effects of rotation and gravogyro-gravothermal collapse, but still provide a relatively small database for rare binaries and escapers (e.g. those involving BHs). We ran the models for 500 Myr, using one non-rotating King model and three rotating ones (Einsel & Spurzem 1999, ω0 = 0.0, 0.6, 1.2, 1.8), each with and without GR merger recoil kicks.

Primordial binary statistics are another major uncertainty for Pop III clusters. In this work we fixed the initial binary fraction to fb = 1% (Table 1), which should be regarded as a low- fb limiting case rather than a best-estimate value. This choice suppresses IMBH formation pathways tied to primordial binary evolution (in particular Type-I and Type-III; Sect. 4.1) and, by construction, increases the relative importance of dynamical formation channels in our event budget. If fb were higher, we would expect (i) a larger contribution of mergers from primordial binary channel, and (ii) stronger binary heating, which can delay gravothermal collapse and modify the balance between incluster hierarchical growth and dynamical ejections. Predicting the net change in the total IMBH production efficiency is therefore not straightforward without additional simulations, but the dominance of dynamical channels found here should be interpreted as an upper bound. A dedicated scan in fb is left for future work.

The global dynamical evolution can be divided into four phases (Figs. 3, A.4, and A.7):

  • Phase I (0-2 Myr): global mass (Mcl) and tidal radius (rtid) remain essentially constant because mass loss is negligible; however, initial bulk rotation (ω0) controls the structural response. Rapid angular-momentum transport occurs earlier and more strongly for larger ω0, producing a brief, synchronous bump in the Lagrangian radii near ~0.1 Myr, inner-radius contraction and outer-radius expansion, and temporary rises then falls in rh, rc, and mc in the faster-rotating models. Low and zero rotation models show little change in rh over this phase. Also, Mgroup evolves similarly across mass groups, while Lgroup2/Lω0=0.6,t=02Mathematical equation: $L_{\mathrm{group}}^2/L^2_{\omega_{0}=0.6,t=0}$ is larger for higher ω0, reflecting early outward angular-momentum migration;

  • Phase II (−2-3 Myr): rapid stellar-evolution mass loss (notably PISNe and compact-object formation) causes a sharp drop in total mass and rtid, and produces a rapid increase in rh and rc. This is due to PISNe both remove mass and carry away a large fraction of the total angular momentum, while escapers (natal or GR kicks) and dynamical encounters add to mass loss and alter kinematics through fallback-dependent kicks. Measures of central segregation (Mav within small rLagr) show a similar qualitative behaviour across ω0 at this stage;

  • Phase III (−3 Myr to −50-80 Myr; the latter is the core collapse time depending on ω0): larger initial rotation shortens the gravothermal collapse timescale, leading to a higher central density and faster mass segregation. Post-phase II angular-momentum diagnostics (e.g. LhCCSNe2/Lω0=0.6,t=02Mathematical equation: $L_{\mathrm{hCCSNe}}^2/L^2_{\omega_{0}=0.6,t=0}$) rise to a peak at ~10-20 Myr in the (NO)K1.2 and (NO)K1.8 runs then decline. The (NO)K1.8 models eject many stars;

  • Phase IV (core collapse to 500 Myr, the simulations’ end): post-collapse evolution is driven by binary heating in the core, producing self-similar Lagrangian-radius expansion and overall cluster expansion; models with larger ω0 collapse earlier and expand more. The system becomes fully mass-segregated. Angular momentum is redistributed from high-mass to low-mass remnants. Low-mass BHs migrate to wider orbits. A compact central IMBH subsystem forms.

We find a dependence of merger rates for stars and compact objects on the initial rotation of Pop III clusters (Fig. A.5), particularly visible in the K models. The largest initial rotation (ω0 = 1.8) corresponds to the largest number of coalescences and collisions (Table B.2), including the most BHBH collisions. Concerning the BHBH collisions, large initial rotation does also seem to precipitate these, with the largest number of collisions found in (NO)K1.8 simulations.

All simulations produce (IM)BHs below, inside, and above the PISN mass gap (Figs. A.3 and 1). The most massive IMBH has 943.19 M⊙ (from NOK1.2; Fig. A.9). Single-star core collapse can form IMBHs above the gap; elsewhere IMBH formation requires binary channels — primarily MS-BH hyperbolic collisions, CHeB or ShHeB coalescence with BHs, or MS-MS mergers in massive primordial binaries. Many BHs that cannot be produced by single star evolution instead form after a common-envelope phase and subsequent coalescence that produces a CHeB then ShHeB star, which then core-collapses into a BH.

With one realization per set-up, we separated our findings into (i) dynamically robust trends dominated by continuous quantities or large-number counts, including global structural evolution and timescales (core-collapse time and post-collapse self-similar expansion; Figs. 3 and A.4), angular-momentum transport signatures (Fig. A.7), and total mass loss as well as total escaper numbers (Table B.3); (ii) potentially stochastic trends dominated by small-number event statistics, such as specific collision or coalescence subtypes and compact-object channels (e.g. BHBH collision counts; Table B.2), rare escaper subtypes (e.g. nescCHeB,s), and single exceptional events (e.g. the lone BHNS merger).

Cluster evolution therefore reflects a competition: gravothermal-gravogyro contraction boosts early merger rates, while the loss of massive objects due to excess angular momentum suppresses subsequent IMBH growth and BH mergers. Other physical effects, which may increase (IM)BH retention even for rotating clusters are: (i) The presence of a dark matter mini-halo around Pop III clusters, as considered by Wang et al. (2022, we note that those simulations did not include rotation). However, Pop III star clusters may also form in larger dark matter halos, and would be dominated by their baryons after cooling; see discussion in Kimm et al. (2016). (ii) Gas infall could add more mass to the central regions, driving stronger contraction than predicted by gravothermal-gravogyro evolution. Such a contraction may trigger collapse of the central rotating IMBH subsystem and produce very massive seed BHs (see Kroupa 2020).

Our results establish a general understanding of how initial rotation, angular momentum transport, and numerous massive objects in the top-heavy IMF affect global evolution, counts of compact objects, and their merger rates. In the future, we aim to (i) continue the simulations for longer time, (ii) increase the number of simulations to improve the statistics of results, (iii) study the influence of an external potential (such as a dark matter mini-halo) on retention rates, and (iv) test a fractal initial structure, which may be present in the early phases.

Acknowledgements

We thank the anonymous referee for their constructive comments that helped to improve the quality of this paper. We thank Shuai Liu for helpful discussions and comments, and acknowledge computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS Booster. KW and RS acknowledge support by the German Science Foundation (DFG, project Sp 345/24-1). AT appreciates support by Grants-in-Aid for Scientific Research (17H06360, 19K03907, 24K07040, 25K01035) from the Japan Society for the Promotion of Science. This material is based upon work supported by Tamkeen under the NYU Abu Dhabi Research Institute grant CASS. MCV acknowledges funding through ANID (Doctorado acuerdo bilateral DAAD/62210038) and DAAD (funding program number 57600326), and by International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). BL acknowledges the funding of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1 - 390900948 (the Heidelberg STRUCTURES Excellence Cluster). MAS acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 101025436 (project GRACE-BH, PI Manuel Arca Sedda) and the support of the MERAC Foundation through the 2023 MERAC prize and support to research. NN is grateful for funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project-ID 138713538 - SFB 881 (“The Milky Way System”, subproject B08). RS has been supported by NAOC/CAS Beijing, International Cooperation Office (2023-2025), by National Natural Science Foundation of China (NSFC) under grant No. 12473017, and by the CAS President’s International Fellowship Initiative for Visiting Scientists (PIFI, 2026PVA0089).

References

  1. Aarseth, S. J. 2012, MNRAS, 422, 841 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aasi, J., Abadie, J., Abbott, B. P., et al. 2015, Class. Quant. Grav., 32, 115012 [NASA ADS] [CrossRef] [Google Scholar]
  3. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, LRR, 21, 3 [Google Scholar]
  4. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, ApJ, 882, L24 [Google Scholar]
  5. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020, LRR, 23, 3 [Google Scholar]
  6. Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Class. Quant. Grav., 32, 024001 [Google Scholar]
  7. Ahmad, A., & Cohen, L. 1973, J. Computat. Phys., 12, 389 [NASA ADS] [CrossRef] [Google Scholar]
  8. Akiyama, K., & Sugimoto, D. 1989, PASJ, 41, 991 [NASA ADS] [Google Scholar]
  9. Amaro-Seoane, P. 2018, PhRvD, 98, 063018 [Google Scholar]
  10. Amaro-Seoane, P., Gair, J. R., Freitag, M., et al. 2007, Class. Quant. Grav., 24, R113 [NASA ADS] [CrossRef] [Google Scholar]
  11. Amaro-Seoane, P., Aoudia, S., Babak, S., et al. 2013, GW Notes, 6, 4 [Google Scholar]
  12. Amaro-Seoane, P., Glaschke, P., & Spurzem, R. 2014, MNRAS, 445, 3755 [Google Scholar]
  13. Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv e-prints [arXiv:1702.00786] [Google Scholar]
  14. Amaro-Seoane, P., Andrews, J., Arca Sedda, M., et al. 2023, Liv. Rev. Relativ., 26, 2 [NASA ADS] [CrossRef] [Google Scholar]
  15. Arca Sedda, M. 2020, CmPhy, 3, 43 [Google Scholar]
  16. Arca-Sedda, M., & Gualandris, A. 2018, MNRAS, 477, 4423 [Google Scholar]
  17. Arca Sedda, M., Askar, A., & Giersz, M. 2019, arXiv e-prints [arXiv:1905.00902] [Google Scholar]
  18. Arca Sedda, M., Berry, C. P. L., Jani, K., et al. 2020a, CQGra, 37, 215011 [Google Scholar]
  19. Arca Sedda, M., Mapelli, M., Spera, M., Benacquista, M., & Giacobbo, N. 2020b, ApJ, 894, 133 [NASA ADS] [CrossRef] [Google Scholar]
  20. Arca Sedda, M., Amaro Seoane, P., & Chen, X. 2021, A&A, 652, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Arca-Sedda, M., Rizzuto, F. P., Naab, T., et al. 2021, ApJ, 920, 128 [NASA ADS] [CrossRef] [Google Scholar]
  22. Arca Sedda, M., Kamlah, A. W. H., Spurzem, R., et al. 2024, MNRAS, 528, 5119 [NASA ADS] [CrossRef] [Google Scholar]
  23. Askar, A., Davies, M. B., & Church, R. P. 2021, MNRAS, 502, 2682 [NASA ADS] [CrossRef] [Google Scholar]
  24. Askar, A., Davies, M. B., & Church, R. P. 2022, MNRAS, 511, 2631 [NASA ADS] [CrossRef] [Google Scholar]
  25. Baker, J. G., Boggs, W. D., Centrella, J., et al. 2008, PhRvD, 78, 044046 [Google Scholar]
  26. Ballone, A., Costa, G., Mapelli, M., & MacLeod, M. 2022, arXiv e-prints [arXiv:2204.03493] [Google Scholar]
  27. Banerjee, S. 2021, MNRAS, 500, 3002 [Google Scholar]
  28. Banerjee, S. 2022, A&A, 665, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Banerjee, S., Belczynski, K., Fryer, C. L., et al. 2020, A&A, 639, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Begelman, M. C. 2010, MNRAS, 402, 673 [NASA ADS] [CrossRef] [Google Scholar]
  31. Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, MNRAS, 370, 289 [NASA ADS] [CrossRef] [Google Scholar]
  32. Begelman, M. C., Rossi, E. M., & Armitage, P. J. 2008, MNRAS, 387, 1649 [NASA ADS] [CrossRef] [Google Scholar]
  33. Belczynski, K., & Banerjee, S. 2020, A&A, 640, L20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Belczynski, K., Heger, A., Gladysz, W., et al. 2016, A&A, 594, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Belczynski, K., Ryu, T., Perna, R., et al. 2017, MNRAS, 471, 4702 [Google Scholar]
  36. Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton University Press) [Google Scholar]
  37. Bobylev, V. V., & Bajkova, A. T. 2020, Astron. Rep., 64, 1042 [Google Scholar]
  38. Bouwens, R. J., Illingworth, G. D., González, V., et al. 2010, ApJ, 725, 1587 [NASA ADS] [CrossRef] [Google Scholar]
  39. Breen, P. G., & Heggie, D. C. 2013, MNRAS, 432, 2779 [NASA ADS] [CrossRef] [Google Scholar]
  40. Bromm, V. 2013, Rep. Progr. Phys., 76, 112901 [CrossRef] [Google Scholar]
  41. Bromm, V., & Loeb, A. 2003, ApJ, 596, 34 [Google Scholar]
  42. Bromm, V., & Larson, R. B. 2004, ARA&A, 42, 79 [Google Scholar]
  43. Campanelli, M., Lousto, C., Zlochower, Y., & Merritt, D. 2007, ApJ, 659, L5 [NASA ADS] [CrossRef] [Google Scholar]
  44. Chattopadhyay, D., Stevenson, S., Hurley, J. R., Bailes, M., & Broekgaarden, F. 2021, MNRAS, 504, 3682 [NASA ADS] [CrossRef] [Google Scholar]
  45. Chon, S., & Omukai, K. 2025, MNRAS, 539, 2561 [Google Scholar]
  46. Chon, S., Omukai, K., & Schneider, R. 2021, MNRAS, 508, 4175 [NASA ADS] [CrossRef] [Google Scholar]
  47. Cohn, H. 1979, ApJ, 234, 1036 [NASA ADS] [CrossRef] [Google Scholar]
  48. Costa, G., Ballone, A., Mapelli, M., & Bressan, A. 2022, arXiv e-prints [arXiv:2204.03492] [Google Scholar]
  49. de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. de Souza, R. S., Ciardi, B., Maio, U., & Ferrara, A. 2013, MNRAS, 428, 2109 [NASA ADS] [CrossRef] [Google Scholar]
  51. Di Carlo, U. N., Giacobbo, N., Mapelli, M., et al. 2019, MNRAS, 487, 2947 [NASA ADS] [CrossRef] [Google Scholar]
  52. Di Carlo, U. N., Mapelli, M., Bouffanais, Y., et al. 2020a, MNRAS, 497, 1043 [NASA ADS] [CrossRef] [Google Scholar]
  53. Di Carlo, U. N., Mapelli, M., Giacobbo, N., et al. 2020b, MNRAS, 498, 495 [NASA ADS] [CrossRef] [Google Scholar]
  54. Di Carlo, U. N., Mapelli, M., Pasquato, M., et al. 2021, MNRAS, 507, 5132 [NASA ADS] [CrossRef] [Google Scholar]
  55. Einsel, C., & Spurzem, R. 1999, MNRAS, 302, 81 [NASA ADS] [CrossRef] [Google Scholar]
  56. Ekström, S., Meynet, G., Chiappini, C., Hirschi, R., & Maeder, A. 2008, A&A, 489, 685 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Evans, M., Adhikari, R. X., Afle, C., et al. 2021, arXiv e-prints [arXiv:2109.09882] [Google Scholar]
  58. Fiestas, J., Porth, O., Berczik, P., & Spurzem, R. 2012, MNRAS, 419, 57 [Google Scholar]
  59. Fraser, M., Casey, A. R., Gilmore, G., Heger, A., & Chan, C. 2017, MNRAS, 468, 418 [NASA ADS] [CrossRef] [Google Scholar]
  60. Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [Google Scholar]
  61. Fujimoto, S., Naidu, R. P., Chisholm, J., et al. 2025, ApJ, 989, 46 [Google Scholar]
  62. Gaia Collaboration (Helmi, A., et al.) 2018, A&A, 616, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Gieles, M., Charbonnel, C., Krause, M. G. H., et al. 2018, MNRAS, 478, 2461 [Google Scholar]
  64. Giersz, M., Leigh, N., Hypki, A., Lützgendorf, N., & Askar, A. 2015, MNRAS, 454, 3150 [Google Scholar]
  65. Glaschke, P., Amaro-Seoane, P., & Spurzem, R. 2014, MNRAS, 445, 3620 [Google Scholar]
  66. Goodman, J. J. 1983, PhD thesis, Princeton University, NJ [Google Scholar]
  67. Goodwin, S. P., & Whitworth, A. P. 2004, A&A, 413, 929 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 [Google Scholar]
  69. Hachisu, I. 1979, PASJ, 31, 523 [Google Scholar]
  70. Hachisu, I. 1982, PASJ, 34, 313 [NASA ADS] [Google Scholar]
  71. Haemmerlé, L., Mayer, L., Klessen, R. S., et al. 2020, SSR, 216, 48 [Google Scholar]
  72. Hartwig, T., Volonteri, M., Bromm, V., et al. 2016, MNRAS, 460, L74 [NASA ADS] [CrossRef] [Google Scholar]
  73. Hassan, J. B., Perna, R., Cantiello, M., et al. 2025, ApJ, 992, 68 [Google Scholar]
  74. Hénon, M. 1975, in IAU Symposium, 69, Dynamics of the Solar Systems, ed. A. Hayli, 133 [Google Scholar]
  75. Henyey, L. G., Wilets, L., Böhm, K. H., Lelevier, R., & Levee, R. D. 1959, ApJ, 129, 628 [Google Scholar]
  76. Hijikawa, K., Tanikawa, A., Kinugawa, T., Yoshida, T., & Umeda, H. 2021, MNRAS [arXiv:2104.13384] [Google Scholar]
  77. Hirano, S., & Bromm, V. 2017, MNRAS, 470, 898 [CrossRef] [Google Scholar]
  78. Hirano, S., & Bromm, V. 2018, MNRAS, 476, 3964 [CrossRef] [Google Scholar]
  79. Hirano, S., Yoshida, N., Sakurai, Y., & Fujii, M. S. 2018, ApJ, 855, 17 [Google Scholar]
  80. Hong, J., Kim, E., Lee, H. M., & Spurzem, R. 2013, MNRAS, 430, 2960 [NASA ADS] [CrossRef] [Google Scholar]
  81. Hosokawa, T., Hirano, S., Kuiper, R., et al. 2016, ApJ, 824, 119 [NASA ADS] [CrossRef] [Google Scholar]
  82. Hughes, S. A. 2009, ARA&A, 47, 107 [Google Scholar]
  83. Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [Google Scholar]
  84. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
  85. Hut, P., Makino, J., & McMillan, S. 1995, ApJ, 443, L93 [NASA ADS] [CrossRef] [Google Scholar]
  86. Inayoshi, K., Hirai, R., Kinugawa, T., & Hotokezaka, K. 2017, MNRAS, 468, 5020 [Google Scholar]
  87. Jani, K., Shoemaker, D., & Cutler, C. 2020, Nat. Astron., 4, 260 [NASA ADS] [CrossRef] [Google Scholar]
  88. Jeena, S. K., Banerjee, P., Chiaki, G., & Heger, A. 2023, MNRAS, 526, 4467 [Google Scholar]
  89. Jiménez-Forteza, X., Keitel, D., Husa, S., et al. 2017, PhRvD, 95, 064024 [Google Scholar]
  90. Kagra Collaboration, (Akutsu, T., et al.) 2019, NatAs, 3, 35 [Google Scholar]
  91. Kamlah, A. W. H., Leveque, A., Spurzem, R., et al. 2022a, MNRAS, 511, 4060 [NASA ADS] [CrossRef] [Google Scholar]
  92. Kamlah, A. W. H., Spurzem, R., Berczik, P., et al. 2022b, MNRAS, 516, 3266 [CrossRef] [Google Scholar]
  93. Katz, H., Sijacki, D., & Haehnelt, M. G. 2015, MNRAS, 451, 2352 [Google Scholar]
  94. Katz, H., Kimm, T., Ellis, R. S., Devriendt, J., & Slyz, A. 2022, arXiv [arXiv:2207.04751] [Google Scholar]
  95. Kawamura, S., Ando, M., Seto, N., et al. 2011, Class. Quant. Grav., 28, 094011 [CrossRef] [Google Scholar]
  96. Kawamura, S., Ando, M., Seto, N., et al. 2021, Progr. Theor. Exp. Phys., 2021, 05A105 [Google Scholar]
  97. Kim, E., Einsel, C., Lee, H. M., Spurzem, R., & Lee, M. G. 2002, MNRAS, 334, 310 [NASA ADS] [CrossRef] [Google Scholar]
  98. Kim, E., Lee, H. M., & Spurzem, R. 2004, MNRAS, 351, 220 [NASA ADS] [CrossRef] [Google Scholar]
  99. Kim, E., Yoon, I., Lee, H. M., & Spurzem, R. 2008, MNRAS, 383, 2 [NASA ADS] [CrossRef] [Google Scholar]
  100. Kiminki, D. C., Kobulnicky, H. A., Ewing, I., et al. 2012, ApJ, 747, 41 [Google Scholar]
  101. Kimm, T., Cen, R., Rosdahl, J., & Yi, S. K. 2016, ApJ, 823, 52 [NASA ADS] [CrossRef] [Google Scholar]
  102. King, I. 1962, AJ, 67, 471 [Google Scholar]
  103. King, I. R. 1966, AJ, 71, 276 [NASA ADS] [CrossRef] [Google Scholar]
  104. Kinugawa, T., Inayoshi, K., Hotokezaka, K., Nakauchi, D., & Nakamura, T. 2014, MNRAS, 442, 2963 [NASA ADS] [CrossRef] [Google Scholar]
  105. Kinugawa, T., Miyamoto, A., Kanda, N., & Nakamura, T. 2016, MNRAS, 456, 1093 [CrossRef] [Google Scholar]
  106. Kinugawa, T., Nakamura, T., & Nakano, H. 2021a, MNRAS, 501, L49 [Google Scholar]
  107. Kinugawa, T., Nakamura, T., & Nakano, H. 2021b, MNRAS, 504, L28 [Google Scholar]
  108. Klessen, R. 2019, in Formation of the First Black Holes, eds. M. Latif, & D. Schleicher (World Scientific Publishing Co.), 67 [Google Scholar]
  109. Klessen, R. S., & Glover, S. C. O. 2023, ARA&A, 61, 65 [Google Scholar]
  110. Kobulnicky, H. A., Kiminki, D. C., Lundquist, M. J., et al. 2014, ApJS, 213, 34 [Google Scholar]
  111. Kokubo, E., Yoshinaga, K., & Makino, J. 1998, MNRAS, 297, 1067 [NASA ADS] [CrossRef] [Google Scholar]
  112. Kroupa, P. 2020, in IAU Symposium, 351, Star Clusters: From the Milky Way to the Early Universe, eds. A. Bragaglia, M. Davies, A. Sills, & E. Vesperini, 117 [Google Scholar]
  113. Kroupa, P. 2025, Contrib. Astron. Observ. Skalnate Pleso, 55, 217 [Google Scholar]
  114. Kroupa, P., Subr, L., Jerabkova, T., & Wang, L. 2020, MNRAS, 498, 5652 [Google Scholar]
  115. Kuepper, A. H. W., Maschberger, T., Kroupa, P., & Baumgardt, H. 2011, McLuster: A Tool to Make a Star Cluster [Google Scholar]
  116. Latif, M. A., Whalen, D., & Khochfar, S. 2022, ApJ, 925, 28 [NASA ADS] [CrossRef] [Google Scholar]
  117. Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints [arXiv:1110.3193] [Google Scholar]
  118. Lazar, A., & Bromm, V. 2022, MNRAS, 511, 2505 [NASA ADS] [CrossRef] [Google Scholar]
  119. Leveque, A., Giersz, M., Arca-Sedda, M., & Askar, A. 2022a, arXiv e-prints [arXiv:2206.03967] [Google Scholar]
  120. Leveque, A., Giersz, M., Banerjee, S., et al. 2022b, arXiv e-prints [arXiv:2206.03404] [Google Scholar]
  121. Liu, B., & Bromm, V. 2020, APJ, 903, L40 [Google Scholar]
  122. Liu, B., & Bromm, V. 2021, MNRAS, 506, 5451 [Google Scholar]
  123. Liu, S., Hu, Y.-M., Zhang, J.-d., & Mei, J. 2020, Phys. Rev. D, 101, 103027 [Google Scholar]
  124. Liu, B., Meynet, G., & Bromm, V. 2021a, MNRAS, 501, 643 [Google Scholar]
  125. Liu, B., Sibony, Y., Meynet, G., & Bromm, V. 2021b, MNRAS, 506, 5247 [NASA ADS] [CrossRef] [Google Scholar]
  126. Liu, B., Hartwig, T., Sartorio, N. S., et al. 2024a, MNRAS, 534, 1634 [NASA ADS] [CrossRef] [Google Scholar]
  127. Liu, S., Wang, L., Hu, Y.-M., Tanikawa, A., & Trani, A. A. 2024b, MNRAS, 533, 2262 [NASA ADS] [CrossRef] [Google Scholar]
  128. Liu, S., Wang, L., Tanikawa, A., Wu, W., & Fujii, M. S. 2025, ApJ, 993, L30 [Google Scholar]
  129. Longaretti, P. Y., & Lagoute, C. 1996, AAP, 308, 453 [Google Scholar]
  130. Lousto, C. O., Zlochower, Y., Dotti, M., & Volonteri, M. 2012, PhRvD, 85, 084015 [Google Scholar]
  131. Luo, J., et al. 2016, Class. Quant. Grav., 33, 035010 [NASA ADS] [CrossRef] [Google Scholar]
  132. Maeder, A. 1987, A&A, 173, 247 [NASA ADS] [Google Scholar]
  133. Maggiore, M., Van Den Broeck, C., Bartolo, N., et al. 2020, JCAP, 2020, 050 [CrossRef] [Google Scholar]
  134. Makino, J. 1991, ApJ, 369, 200 [Google Scholar]
  135. Makino, J. 1999, J. Computat. Phys., 151, 910 [Google Scholar]
  136. Maliszewski, K., Giersz, M., Gondek-Rosinska, D., Askar, A., & Hypki, A. 2022, MNRAS [arXiv:2111.09223] [Google Scholar]
  137. Marques-Chaves, R., Schaerer, D., Kuruvanthodi, A., et al. 2024, A&A, 681, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  138. McMillan, S. L. W. 1986, The Vectorization of Small-N Integrators, 267, eds. P. Hut, & S. L. W. McMillan, (Springer Berlin, Heidelberg), 156 [Google Scholar]
  139. Mestichelli, B., Mapelli, M., Torniamenti, S., et al. 2024, A&A, 690, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  140. Mestichelli, B., Mapelli, M., Santoliquido, F., et al. 2025, A&A, 704, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  141. Mikkola, S., & Aarseth, S. J. 1998, NA, 3, 309 [Google Scholar]
  142. Mikkola, S., & Tanikawa, K. 1999a, MNRAS, 310, 745 [NASA ADS] [CrossRef] [Google Scholar]
  143. Mikkola, S., & Tanikawa, K. 1999b, Celest. Mech. Dyn. Astron., 74, 287 [NASA ADS] [CrossRef] [Google Scholar]
  144. Miller, M. C., & Hamilton, D. P. 2002, MNRAS, 330, 232 [Google Scholar]
  145. Morawski, J., Giersz, M., Askar, A., & Belczynski, K. 2018, MNRAS, 481, 2168 [Google Scholar]
  146. Morawski, J., Giersz, M., Askar, A., & Belczynski, K. 2019, MNRAS, 486, 3402 [Google Scholar]
  147. Naidu, R. P., Oesch, P. A., Brammer, G., et al. 2026, Open J. Astrophys., 9, 56033 [Google Scholar]
  148. Nakauchi, D., Inayoshi, K., & Omukai, K. 2020, ApJ, 902, 81 [Google Scholar]
  149. Nandal, D., Sibony, Y., & Tsiatsiou, S. 2024, A&A, 688, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  150. Neumayer, N., Seth, A., & Böker, T. 2020, A&ARv, 28, 4 [Google Scholar]
  151. Ng, K. K. Y., Franciolini, G., Berti, E., et al. 2022, arXiv e-prints [arXiv:2204.11864] [Google Scholar]
  152. Nitadori, K., & Aarseth, S. J. 2012, MNRAS, 424, 545 [NASA ADS] [CrossRef] [Google Scholar]
  153. Paiella, L., Arca Sedda, M., Mestichelli, B., & Ugolini, C. 2026, A&A, 708, A200 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  154. Peters, P. C. 1964, PhRv, 136, 1224 [Google Scholar]
  155. Peters, P. C., & Mathews, J. 1963, PhRv, 131, 435 [Google Scholar]
  156. Plunkett, C., Mould, M., & Vitale, S. 2025, Phys. Rev. D, 112, 023039 [Google Scholar]
  157. Portegies Zwart, S. F., & McMillan, S. L. W. 2002, ApJ, 576, 899 [Google Scholar]
  158. Punturo, M., Abernathy, M., Acernese, F., et al. 2010a, CQGra, 27, 194002 [Google Scholar]
  159. Punturo, M., Abernathy, M., Acernese, F., et al. 2010b, CQGra, 27, 084007 [Google Scholar]
  160. Rantala, A., Lahén, N., Naab, T., Escobar, G. J., & Iorio, G. 2025, MNRAS, 543, 2130 [Google Scholar]
  161. Rees, M. J. 1984, ARA&A, 22, 471 [Google Scholar]
  162. Reinoso, B., Schleicher, D. R. G., Fellhauer, M., Klessen, R. S., & Boekholt, T. C. N. 2018, AAP, 614, A14 [Google Scholar]
  163. Reinoso, B., Leigh, N. W. C., Barrera-Retamal, C. M., et al. 2022, MNRAS, 509, 3724 [Google Scholar]
  164. Reinoso, B., Latif, M. A., & Schleicher, D. R. G. 2025, A&A, 700, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  165. Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019, in Bulletin of the American Astronomical Society, 51, 35 [Google Scholar]
  166. Renaud, F., & Gieles, M. 2015, MNRAS, 448, 3416 [NASA ADS] [CrossRef] [Google Scholar]
  167. Rezzolla, L., Barausse, E., Dorband, E. N., et al. 2008, PhRvD, 78, 044002 [Google Scholar]
  168. Rizzuto, F. P., Naab, T., Spurzem, R., et al. 2021, MNRAS, 501, 5257 [NASA ADS] [CrossRef] [Google Scholar]
  169. Rizzuto, F. P., Naab, T., Spurzem, R., et al. 2022, MNRAS, 512, 884 [NASA ADS] [CrossRef] [Google Scholar]
  170. Rodriguez, C. L., Zevin, M., Amaro-Seoane, P., et al. 2019, PhRvD, 100, 043027 [Google Scholar]
  171. Rydberg, C.-E., Zackrisson, E., Lundqvist, P., & Scott, P. 2013, MNRAS, 429, 3658 [NASA ADS] [CrossRef] [Google Scholar]
  172. Sakurai, Y., Hosokawa, T., Yoshida, N., & Yorke, H. W. 2015, MNRAS, 452, 755 [Google Scholar]
  173. Sakurai, Y., Yoshida, N., Fujii, M. S., & Hirano, S. 2017, MNRAS, 472, 1677 [NASA ADS] [CrossRef] [Google Scholar]
  174. Sakurai, Y., Yoshida, N., & Fujii, M. S. 2019, MNRAS, 484, 4665 [Google Scholar]
  175. Sana, H., & Evans, C. J. 2011, in IAU Symposium, 272, Active OB Stars: Structure, Evolution, Mass Loss, and Critical Limits, eds. C. Neiner, G. Wade, G. Meynet, & G. Peters, 474 [Google Scholar]
  176. Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
  177. Sana, H., de Mink, S. E., de Koter, A., et al. 2013, in Astronomical Society of the Pacific Conference Series, 470, 370 Years of Astronomy in Utrecht, eds. G. Pugliese, A. de Koter, & M. Wijburg, 141 [Google Scholar]
  178. Santoliquido, F., Mapelli, M., Iorio, G., et al. 2023, MNRAS, 524, 307 [NASA ADS] [CrossRef] [Google Scholar]
  179. Sartorio, N. S., Fialkov, A., Hartwig, T., et al. 2023, MNRAS, 521, 4039 [Google Scholar]
  180. Sathyaprakash, B., Abernathy, M., Acernese, F., et al. 2012, CQGra, 29, 124013 [Google Scholar]
  181. Schauer, A. T. P., Bromm, V., Drory, N., & Boylan-Kolchin, M. 2022, arXiv e-prints [arXiv:2207.02863] [Google Scholar]
  182. Schneider, R., Ferrara, A., Ciardi, B., Ferrari, V., & Matarrese, S. 2000, MNRAS, 317, 385 [Google Scholar]
  183. Schneider, R., Ferrara, A., Natarajan, P., & Omukai, K. 2002, ApJ, 571, 30 [CrossRef] [Google Scholar]
  184. Schneider, R., Ferrara, A., Salvaterra, R., Omukai, K., & Bromm, V. 2003, Nature, 422, 869 [NASA ADS] [CrossRef] [Google Scholar]
  185. Schleicher, D. R. G., Reinoso, B., Latif, M., et al. 2022, MNRAS, 512, 6192 [NASA ADS] [CrossRef] [Google Scholar]
  186. Sharda, P., & Krumholz, M. R. 2022, MNRAS, 509, 1959 [Google Scholar]
  187. Sharda, P., & Menon, S. H. 2025, MNRAS, 540, 1745 [Google Scholar]
  188. Sharda, P., Federrath, C., Krumholz, M. R., & Schleicher, D. R. G. 2021, MNRAS, 503, 2014 [NASA ADS] [CrossRef] [Google Scholar]
  189. Shibata, M., & Shapiro, S. L. 2002, ApJ, 572, L39 [Google Scholar]
  190. Spurzem, R. 1994, Gravothermal Oscillations, 430, eds. V. G. Gurzadyan, & D. Pfenniger (Springer Berlin, Heidelberg), 170 [Google Scholar]
  191. Spurzem, R. 1996, in IAU Symposium, 174, Dynamical Evolution of Star Clusters: Confrontation of Theory and Observations, eds. P. Hut, & J. Makino, 111 [Google Scholar]
  192. Spurzem, R. 1999, JCoAM, 109, 407 [Google Scholar]
  193. Spurzem, R., & Kamlah, A. 2023, Liv. Rev. Computat. Astrophys., 9, 3 [Google Scholar]
  194. Stacy, A., & Bromm, V. 2013, MNRAS, 433, 1094 [NASA ADS] [CrossRef] [Google Scholar]
  195. Stacy, A., Bromm, V., & Loeb, A. 2011, MNRAS, 413, 543 [NASA ADS] [CrossRef] [Google Scholar]
  196. Stacy, A., Greif, T. H., Klessen, R. S., Bromm, V., & Loeb, A. 2013, MNRAS, 431, 1470 [NASA ADS] [CrossRef] [Google Scholar]
  197. Stacy, A., Bromm, V., & Lee, A. T. 2016, MNRAS, 462, 1307 [NASA ADS] [CrossRef] [Google Scholar]
  198. Stiefel, E., & Kustaanheimo, P. 1965, J. Reine Angew. Math., 218, 204 [Google Scholar]
  199. Stone, N. C., Küpper, A. H. W., & Ostriker, J. P. 2017, MNRAS, 467, 4180 [NASA ADS] [Google Scholar]
  200. Sugimura, K., Matsumoto, T., Hosokawa, T., Hirano, S., & Omukai, K. 2020, ApJ, 892, L14 [NASA ADS] [CrossRef] [Google Scholar]
  201. Susa, H. 2019, ApJ, 877, 99 [NASA ADS] [CrossRef] [Google Scholar]
  202. Tacchella, S., Eisenstein, D. J., Hainline, K., et al. 2023, ApJ, 952, 74 [NASA ADS] [CrossRef] [Google Scholar]
  203. Takahashi, K., Yoshida, T., Umeda, H., Sumiyoshi, K., & Yamada, S. 2016, MNRAS, 456, 1320 [NASA ADS] [CrossRef] [Google Scholar]
  204. Takahashi, K., Yoshida, T., & Umeda, H. 2018, ApJ, 857, 111 [NASA ADS] [CrossRef] [Google Scholar]
  205. Takahashi, K., Sumiyoshi, K., Yamada, S., Umeda, H., & Yoshida, T. 2019, ApJ, 871, 153 [Google Scholar]
  206. Tanikawa, A., Yoshida, T., Kinugawa, T., Takahashi, K., & Umeda, H. 2020, MNRAS, 495, 4170 [NASA ADS] [CrossRef] [Google Scholar]
  207. Tanikawa, A., Kinugawa, T., Yoshida, T., Hijikawa, K., & Umeda, H. 2021a, MNRAS [arXiv:2010.07616] [Google Scholar]
  208. Tanikawa, A., Susa, H., Yoshida, T., Trani, A. A., & Kinugawa, T. 2021b, ApJ, 910, 30 [NASA ADS] [CrossRef] [Google Scholar]
  209. Tanikawa, A., Chiaki, G., Kinugawa, T., Suwa, Y., & Tominaga, N. 2022a, PASJ, 74, 521 [NASA ADS] [CrossRef] [Google Scholar]
  210. Tanikawa, A., Moriya, T. J., Tominaga, N., & Yoshida, N. 2022b, arXiv e-prints [arXiv:2204.09402] [Google Scholar]
  211. Tanikawa, A., Yoshida, T., Kinugawa, T., et al. 2022c, ApJ, 926, 83 [NASA ADS] [CrossRef] [Google Scholar]
  212. The LIGO Scientific Collaboration (the Virgo Collaboration, the KAGRA Collaboration). 2025, arXiv e-prints [arXiv:2508.20721] [Google Scholar]
  213. Tsiatsiou, S., Sibony, Y., Nandal, D., et al. 2024, A&A, 687, A307 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  214. Uchida, H., Shibata, M., Yoshida, T., Sekiguchi, Y., & Umeda, H. 2017, PhRvD, 96, 083016 [Google Scholar]
  215. van Dokkum, P., Brammer, G., Jennings, C., Pasha, I., & Baggen, J. F. W. 2025, ApJ, 990, L48 [Google Scholar]
  216. van Meter, J. R., Miller, M. C., Baker, J. G., Boggs, W. D., & Kelly, B. J. 2010, ApJ, 719, 1427 [Google Scholar]
  217. Vanzella, E., Meneghetti, M., Caminha, G. B., et al. 2020, MNRAS, 494, L81 [Google Scholar]
  218. Varri, A. L., & Bertin, G. 2012, AAP, 540, A94 [Google Scholar]
  219. Vergara, M. Z. C., Schleicher, D. R. G., Boekholt, T. C. N., et al. 2021, A&A, 649, A160 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  220. Vergara, M. C., Escala, A., Schleicher, D. R. G., & Reinoso, B. 2023, MNRAS, 522, 4224 [NASA ADS] [CrossRef] [Google Scholar]
  221. Vergara, M. C., Askar, A., Kamlah, A. W. H., et al. 2025, A&A, 704, A321 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  222. Vergara, M. C., Askar, A., Flammini Dotti, F., et al. 2026, A&A, 707, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  223. Vikaeus, A., Whalen, D. J., & Zackrisson, E. 2022, ApJ, 933, L8 [CrossRef] [Google Scholar]
  224. Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070 [NASA ADS] [CrossRef] [Google Scholar]
  225. Wang, L., Iwasawa, M., Nitadori, K., & Makino, J. 2020, MNRAS, 497, 536 [NASA ADS] [CrossRef] [Google Scholar]
  226. Wang, L., Tanikawa, A., & Fujii, M. 2022, arXiv e-prints [arXiv:2207.09621] [Google Scholar]
  227. Wang, X., Cheng, C., Ge, J., et al. 2024, ApJ, 967, L42 [NASA ADS] [CrossRef] [Google Scholar]
  228. Wang, H., Liu, S., Wang, H., et al. 2025a, A&A, 697, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  229. Wang, Z., Ma, Y., Li, Y., et al. 2025b, ApJ, 990, 160 [Google Scholar]
  230. Webb, J. J., Leigh, N. W. C., Serrano, R., et al. 2019, MNRAS, 488, 3055 [NASA ADS] [CrossRef] [Google Scholar]
  231. Woosley, S. E. 2017, ApJ, 836, 244 [Google Scholar]
  232. Woosley, S. E., & Heger, A. 2021, ApJ, 912, L31 [NASA ADS] [CrossRef] [Google Scholar]
  233. Wu, W., Wang, L., Liu, S., et al. 2025, ApJ, 986, 163 [Google Scholar]
  234. Yoon, S. C., & Langer, N. 2005, A&A, 443, 643 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  235. Yoshida, T., Takiwaki, T., Kotake, K., et al. 2019, ApJ, 881, 16 [NASA ADS] [CrossRef] [Google Scholar]
1

The Kobulnicky et al. (2014) distribution adopts a random pairing for stars with m ≤ 5 M⊙, but due to our IMF set-up, all stars have m > 5 M⊙.

2

Core collapse time is ∼50-80 Myr depending on ω0; see Sect. 4.3.1.

3

A primordial binary may have members of two different mass groups.

5

Stellar type acronyms follow Hurley et al. (2000); also see Sect. C.3.

6

Core collapse is at 50-80 Myr; see Fig. A.4.

7

See the next paper in this series.

8

By contrast, trends inferred from continuous structural quantities (e.g. Figs. 3 and A.4, and the angular-momentum transport in Fig. A.7) and from large-number counts (e.g. total escaper numbers; Table B.3) are dynamically more robust in our dataset.

9

From these theories, we would expect that mergers have a very low GRK velocity, which would make them promising candidates for growing seed BHs for galactic nuclei.

10

A surrounding dark matter halo is unlikely to help to retain merger remnants. Even assuming a NFW big-halo of 100 times of the cluster total mass with 10 times of the cluster tidal radius will only contribute to a few km s−1 of the escape speed in the core.

11

Escapers are counted after they pass twice the current tidal radius, so the actual escape process starts earlier than the figure suggests.

Appendix A Additional figures

Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Six vertically stacked panels display the evolution of total star count and the counts of MS (main sequence star), CHeB (core-helium burning star), ShHeB (shell-helium burning star), NS (neutron star), and BH (black hole), for eight simulations.
Thumbnail: Fig. A.2 Refer to the following caption and surrounding text. Fig. A.2

Top to bottom: Number of total binaries (nb), compact binaries (ncb), BHBHs (nbhbh), BH stars (nBHstar), and NS binaries (nNSb) for eight simulations.

Appendix B Tables of object counts

All tables in this section report raw integer counts from a single realization per set-up. For any event or object-count N, the intrinsic counting uncertainty is approximately Poisson, with a reference 68% (≈ 1σ) error σN≃NMathematical equation: $\sigma_N\simeq\sqrt{N}$ (relative error ∼eq1/NMathematical equation: ${\sim}eq 1/\sqrt{N}$). For small N, confidence intervals are non-Gaussian and can be asymmetric; therefore, differences of a few events between models should be interpreted cautiously when N is in the single-digit or few-tens regime.

Thumbnail: Fig. A.3 Refer to the following caption and surrounding text. Fig. A.3

Initial-final-mass relation (IFMR): the BH mass mBH versus the progenitor ZAMS star mass mZAMS for all eight simulations. Columns from left to right show increasing rotation ω0 = 0.0, 0.6, 1.2, 1.8. K runs use a white background; NOK runs are shaded light grey. Blue lines separate the progenitor-mass regimes, as annotated in red. Orange lines enclose the pair-instability mass gap (∼64-161 M⊙, Woosley & Heger 2021). IMBHs appear below, within and above the PISN gap.
Table B.1

Remaining binary populations at 500 Myr.

Table B.2

Counts of coalescences and collisions in the eight simulations.

Appendix C Methodology details

C.1 NBODY6++GPU

NBODY6++GPU12 is optimized for large-scale computing clusters by using MPI (Spurzem 1999), SIMD, OpenMP and GPU (Nitadori & Aarseth 2012; Wang et al. 2015) parallelization techniques. In combination with the Kustaanheimo-Stiefel (KS) regularization (Stiefel & Kustaanheimo 1965) for binary stars and close encounters, the Hermite scheme with hierarchical block time steps (McMillan 1986; Hut et al. 1995; Makino 1991, 1999) and the Ahmad-Cohen (AC) neighbour scheme (Ahmad & Cohen 1973), two-body and chain regularization (Mikkola & Tanikawa 1999b,a; Mikkola & Aarseth 1998) the code allows for highly accurate simulations to be made of star clusters of a realistic size, including binaries.

Thumbnail: Fig. A.4 Refer to the following caption and surrounding text. Fig. A.4

Lagrange radii (rLagr) and average mass (Mav) within spheres that contain 1%, 5%, 10%, 30%, 50%, and 90% of the total cluster mass at the current simulation time step for up to 500 Myr. Time is shown on a logarithmic scale to highlight the cluster’s rapid early evolution. Each column represents one rotational parameter ω0 of the rotating King model (ω0 = 0.0, 0.6, 1.2, 1.8). Models with K are plotted on a white background; models without GR kicks (NOK) are shaded light grey.
Table B.3

Escaping-object counts at 500 Myr in the eight simulations.

Thumbnail: Fig. A.5 Refer to the following caption and surrounding text. Fig. A.5

Time evolution of cumulative number of coalescences ncoal (top two rows) and collisions ncoll (bottom two rows) for the eight runs. Columns give increasing ω0 from left to right. K models (GR kicks on) use a white background; NOK models are shaded light grey. The horizontal axis is logarithmic in time (Myr) to clarify the pre-core-collapse phase. Colours indicate outcomes: MS (blue), CHeB (purple), ShHeB (orange), BH (black); collisions BH-BH (black) and MS-MS (cyan). Counts are in Table B.2

The AC scheme divide the gravitational forces acting on it into the regular component, originating from distant stars, and an irregular part, originating from nearby stars ("neighbours"). Regular forces, efficiently accelerated on the GPU, are updated in larger regular time steps, while neighbour forces are much more fluctuating and need updates in shorter time intervals. Since neighbour numbers are usually small compared to the total particle number, their implementation on the CPU using OpenMP (Wang et al. 2015) provides the best overall performance. PostNewtonian dynamics of relativistic binaries is currently still using the orbit-averaged Peters & Matthews formalism (Peters & Mathews 1963; Peters 1964), as described e.g. in Di Carlo et al. (2019, 2020a,b, 2021); Rizzuto et al. (2021, 2022); Arca-Sedda et al. (2021), but many upgrades were made in Arca Sedda et al. (2024) and this work.

We refer interested readers to Spurzem & Kamlah (2023) for a detailed review of NBODY6++GPU.

Thumbnail: Fig. A.6 Refer to the following caption and surrounding text. Fig. A.6

Cumulative counts of escaping binary members nesc,binarymembers and of binary components: nescMS,b, nescCHeB,b, nescShHeB,b, and nescBH,b. K: solid curves; NOK: dash-dotted curves.

C.2 MCLUSTER and Fopax: Initialization of rotating star clusters

  1. The star clusters models are first initialized with MCLUSTER (Kuepper et al. 2011; Kamlah et al. 2022a; Leveque et al. 2022b). The MCLUSTER output models are ready as an input for NBODY6++GPU as initial models. The input parameters are given in Sect. 3 and Table 1.

  2. Then, we generate a 2-D Fokker-Planck initial models using FOPAX (Einsel & Spurzem 1999; Kim et al. 2002, 2004, 2008). The code produces a 2-D mesh based output of density ρ and velocity dispersions σ as a function of r and z (as cylindrical coordinates) based on the rotating King model f(E, Jz) that are characterized by a pair of parameters (W0,ω0): frk∝(eβE−1)×e−βΩ0Jz,Mathematical equation: $f_{\mathrm{rk}}~\propto~\left(\mathrm{e}^{\beta E}-1\right)\times~\mathrm{e}^{-\beta\Omega_{0}J_{\mathrm{z}}},$ where β=1/(mσc2)Mathematical equation: $\beta = 1/(m\sigma_{\mathrm{c}}^2)$ and the dimensionless rotation parameter is given by ω0=9/4×πGnc×Ω0Mathematical equation: $\omega_{0}=\sqrt{9/4\times \pi Gn_{\mathrm{c}}}\times \Omega_{0}$. σc nc are central onedimensional velocity dispersion and the central density, respectively. Potential-density pairs (e.g. Binney & Tremaine 2008) for these models are created by relating β to the King parameter W0 via W0 = βm(ψ - ψt), where ψ and ψt are the central King potential and the King potential at the truncation radius, which is normally characterized by rtid, the tidal radius (for numerical and computational methods of this step, see also Henyey et al. 1959; Cohn 1979; Spurzem 1994, 1996).

  3. We use a Monte Carlo rejection technique to generate a discrete system of N particles, following the known distributions of ρ and σ. The output, ready for NBODY6++GPU, consist of mass, the 3-D position, velocity data of each particle. This creates a cluster with a differential rotation based on the rejection technique (it includes different rotating shells, which separate a certain number of particles, based on the used King model, ω0 and the number of particles). Finally, all data is scaled to standard Hénon units. As a result, we have an initial star cluster model with rotating King model N-body distribution, the chosen IMF, and relevant binary orbital parameter distributions conserved from MCLUSTER.

C.3 SSE, BSE, and Pop III stellar evolution

The role of stellar mass is extremely important, especially at larger masses. For stars less massive than 8 M⊙, we use the classical stellar evolution recipes of Hurley et al. (2000, 2002) with recent upgrades in Kamlah et al. (2022a). Adjustments for these recipes are applied where extremely metal-poor (EMP; Pop III) stars evolve differently due to the absence of heavy elements. Therefore, we assume EMP stars have negligible mass loss from winds and pulsations (Nakauchi et al. 2020) at extremely low metallicities. Since we use top-heavy IMFs, thus using numerous stars with M/M⊙ > 8, the classical stellar evolution recipes fail to model Pop III stars. In this paper, we present the implementation into NBODY6++GPU of the extended fitting formulae derived from the fitting to 1-D HOSHI stellar evolution models (Takahashi et al. 2016, 2018, 2019; Yoshida et al. 2019) of extremely metal-poor stars (Tanikawa et al. 2020, 2021b,a; Hijikawa et al. 2021). At the metallicity that we adopt in this work, Z/Z⊙ = 10−8, the fitting formulae by Tanikawa et al. (2020) are generally valid from 8 to 1280 M⊙. We adopt the extrapolation procedure of Wang et al. (2022) for stars outside this range. Extrapolation to stellar mass beyond 1280 M⊙ till several 103 M⊙ can be done with these formulae, because no abrupt changes in the fitting are expected for these large masses (e.g. stellar luminosity is almost proportional to stellar mass and the star lifetime is quite similar at these masses) Extrapolation is more difficult below 8 M⊙, where stellar mass has significantly greater effects on stellar evolution. For stars with masses below 8 M⊙, we adopt the evolutionary prescriptions of Hurley et al. (2000) at the lowest available metallicity of Z/Z⊙ = 10−2 Since only a small fraction of stars fall into this mass range under our initial conditions (see Table 1), this choice does not significantly affect the accuracy of our models. To summarize, besides Wang et al. (2022), we are among the first groups to combine direct N-body simulation with full Pop III stellar evolution from Tanikawa et al. (2020).

Thumbnail: Fig. A.7 Refer to the following caption and surrounding text. Fig. A.7

Top two rows: Total mass of the four mass groups (MICCSNE in red, MPPISNE in green, MPISNe in orange, MhCCSNe in blue). Bottom two rows: Each group’s angular-momentum squared, Lgroup2Mathematical equation: $L_{\mathrm{group}}^2$, normalized by Lω0=0.6,t=02Mathematical equation: $L^2_{\omega_{0}=0.6,t=0}$ (i.e. Lgroup2/Lω0=0.6,t=02Mathematical equation: $L_{\mathrm{group}}^2/L^2_{\omega_{0}=0.6,t=0}$). Columns show increasing ω0 from left to right (ω0 = 0.0, 0.6, 1.2, 1.8). K runs are on a white background; NOK runs are shaded light grey.

We briefly highlight the implementations in the evolution of the stars for Z/Z⊙ = 10−8, imported from Tanikawa et al. (2020). Normally, the stellar evolution is divided into five distinct phases in chronological order: Main Sequence (MS), Hertzsprung gap (HG), Core-Helium Burning (CHeB), Shell-Helium Burning (ShHeB) and the remnant phases of either NS or BH. In addition, Tanikawa et al. (2020) define the BSG and RSG phases that relate to the stellar surface, i.e. the effective temperature Teff, of the stars. Stars with log10(Teff/K) ≥ 3.65 are classified as BSGs. Likewise, stars with log10(Teff/K) < 3.65 are classified with RSGs.

The BSG phase begins at ZAMS and ends when log10(Teff/K) < 3.65. Subsequently, the RSG phase begins, which ends at the time of carbon ignition, which is the end point of the 1-D simulations used for the fitting formulae. For example, depending on Teff a CHeB star can become a BSG or a RSG.

For metallicity of Z/Z⊙ = 10−8, as shown in Tab. 2 of Tanikawa et al. 2020:

  • Stars with 8 < M/M⊙ < 13 have entered into their ShHeB phases by the time they become RSG stars.

  • Stars with masses of 13 ≤ M/M⊙ < 50 end with BSG stars.

  • Stars with M/M⊙ ≥ 50 are still CHeB stars when they become RSG stars.

  • In the case of binary stars, it is possible that the H envelope of the CHeB star is fully stripped by tides, Roche Lobe overflow (RLOF) or CEE changing the type of the star to a naked Helium MS (HeMS) star. For these stars, available fitting formulae of lowest metallicity of Z/Z⊙ = 5 × 10−3 from Hurley et al. (2002) are used instead.

We refer readers to Tanikawa et al. (2020) for detailed stellar evolution parameters, or check our code publicly available on GitHub13. Using a time step parameter η = 0.01 (Spurzem & Kamlah 2023), the local error scales as √η. The relative energy error is ~10−8 per N-body time unit (~104 yr). Due to the uncorrelated nature of these errors, the global accumulation is slow, reaching only ~10−6 after 100 Myr. Given the system’s chaotic nature, this precision is sufficient. Moreover, iterative Hermite schemes can further improve energy conservation by achieving near time-symmetry (e.g. Kokubo et al. 1998; Glaschke et al. 2014; Amaro-Seoane et al. 2014; Spurzem & Kamlah 2023).

Thumbnail: Fig. A.8 Refer to the following caption and surrounding text. Fig. A.8

All coalescence (top two rows) and all collision events (bottom two rows) from the eight N-body simulations. Plotted are product masses (mcoal or mcoll) versus distance to the cluster density centre rdens. Columns from left to right show increasing rotation ω0 = 0.0, 0.6, 1.2, 1.8. K models use a white background; NOK models are shaded light grey. Coalescence outcomes: MS (blue), CHeB (purple), ShHeB (orange), BH (black). Collision markers: BH-BH (black cross), MS-MS (cyan cross), BH-MS (violet cross). Counts in Table B.2.

C.4 General relativistic merger recoil kicks (GRKs)

The latest studies of IMBH growth (e.g. Di Carlo et al. 2019, 2020a,b, 2021; Rizzuto et al. 2021, 2022) and star cluster dynamics with NBODY6++GPU (e.g. Kamlah et al. 2022a,b) do not include a GR merger recoil treatment (in addition to missing PN terms). Arca-Sedda et al. (2021) included GRKs via a posteriori analysis, but this does not replace a fully self-consistent treatment during the simulation. The updated NBODY6++GPU presented in this paper and NBODY7 (Aarseth 2012; Banerjee et al. 2020; Banerjee 2021) contain a proper treatment of such kicks (also see Arca Sedda et al. 2024). These velocity kicks depend on spins and mass ratio, due to asymmetric GW radiation during the final inspiral and merger. Numerical relativity (NR) models (Campanelli et al. 2007; Rezzolla et al. 2008; Hughes 2009; van Meter et al. 2010; Jiménez-Forteza et al. 2017) have been used to formulate semi-analytic descriptions for MOCCA and NBODY code family (Morawski et al. 2018, 2019; Banerjee 2021; Belczynski & Banerjee 2020; Arca-Sedda et al. 2021; Banerjee 2022). The smaller the spin (e.g. FULLER model, (nearly) non-spinning BHs), the larger the kick velocity. The larger the mass ratio, the smaller the kick velocity (Morawski et al. 2018, 2019). In the other extreme case with non-aligned natal spins and small mass ratios, the asymmetry in the GW may produce GR merger recoils that reach of order 1000 km s−1 (Baker et al. 2008; van Meter et al. 2010).

Typically, the BHBH orbital angular momentum sets the angular momentum budget for the post-merger BH, so the final spin tends to align with the orbital angular momentum (Banerjee 2021). By contrast, in physical collisions or mergers during binary-single interactions, the orbital angular momentum does not dominate, and therefore the final BH spin can be low. Banerjee (2021) also includes a treatment for random isotropic spin alignment of dynamically formed BHs. Additionally, Banerjee (2021) sets the GRK to zero for NSNS and BHNS mergers (Arca Sedda 2020; Chattopadhyay et al. 2021), but applies numerical-relativity fitting formulae from van Meter et al. (2010, updated by Banerjee 2022) to assign kicks to BHBH merger products. The final spin of the merger product is then evaluated in the same way as a BHBH merger.

Thumbnail: Fig. A.9 Refer to the following caption and surrounding text. Fig. A.9

Primary vs secondary mass (m1 > m2) for collisions from all eight simulations: MS-MS (orange), BH-MS (purple), BH-BH (black). Columns show ω0 = 0.0, 0.6, 1.2, 1.8 left-to-right. K models are shaded white; NOK models are shaded grey. The red line denotes equal mass ratio (q = 1). Counts in Table B.2.
Thumbnail: Fig. A.10 Refer to the following caption and surrounding text. Fig. A.10

GRK velocity, vGRk, versus mass ratio, q, of the participating compact objects. All but one event are BHBH mergers; the exception is a BHNS merger, marked by an extreme q and low vGRk. The shaded part of the plot corresponds to the interval 26.9 - 159.4 km s−1, the escape speeds range in cluster cores among all models.

The GRK-enabled models of this work, for NS-NS, BH-BH, and NS-BH mergers, we draw the remnant spin from a Maxwellian distribution with σ = 0.2 and apply a mass loss due to GW emission: the remnant mass is 98.5% of the summed pre-merger compact object masses.

For further details of the implementation, we refer readers to Arca Sedda et al. (2020b) and Arca Sedda et al. (2024).

All Tables

Table 1

Initial parameters that are identical across all eight initial models for the NBODY6++GPU simulations.

In the text
Table 2

Model identifiers, rotation parameters, and whether simulations enable GR recoil kicks at relativistic mergers.

In the text
Table B.1

Remaining binary populations at 500 Myr.

In the text
Table B.2

Counts of coalescences and collisions in the eight simulations.

In the text
Table B.3

Escaping-object counts at 500 Myr in the eight simulations.

In the text

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Initial-final-mass relation (IFMR) for the NoK0.0 run: BH mass, mBH, versus progenitor ZAMS mass, mZAMS. The panels correspond to progenitors exploding as lCCSNe, PPISNe, PISNe, and hCCSNe. (IM)BH formation channels are shown in red.
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Binary distributions at 500 Myr. From top to bottom: e2; cumulative semi-major axis, a (au); cumulative binary potential energy, Epot (N-body units); cumulative distance to the density centre, rdens (parsecs); and mass ratio q ≡ m2/m1. In the rdens panel a, blue line shows the half-mass radius rh of the NOK0.0 model at 500 Myr (∼9.77 pc).
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Total cluster mass, Mcl, tidal radius, rtid, half mass radius, rh, mass of the core, mc, and radius of the core, rc, in the four panels for all eight simulations with and without GR recoil kicks for different ω0. Time is shown on a logarithmic scale to highlight the cluster’s rapid early evolution. The models with K are plotted as solid curves and the models without (NOK models) are plotted as dash-dotted curves.
In the text
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Escaping single stars in eight simulations: cumulative counts and masses. Top panel: cumulative number of escaping single stars, nesc,singles. Second panel: total mass in escaping single stars, Mesc,singles. Remaining panels: cumulative counts of escaping single MS, CHeB, ShHeB, NS, and BH stars, nescMS,s, nescCHeB,s, nescShHeB,s, nescNS,s, and nescBH,s. Solid curves show K models and dash-dotted curves show NOK models.
In the text
Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Six vertically stacked panels display the evolution of total star count and the counts of MS (main sequence star), CHeB (core-helium burning star), ShHeB (shell-helium burning star), NS (neutron star), and BH (black hole), for eight simulations.
In the text
Thumbnail: Fig. A.2 Refer to the following caption and surrounding text. Fig. A.2

Top to bottom: Number of total binaries (nb), compact binaries (ncb), BHBHs (nbhbh), BH stars (nBHstar), and NS binaries (nNSb) for eight simulations.
In the text
Thumbnail: Fig. A.3 Refer to the following caption and surrounding text. Fig. A.3

Initial-final-mass relation (IFMR): the BH mass mBH versus the progenitor ZAMS star mass mZAMS for all eight simulations. Columns from left to right show increasing rotation ω0 = 0.0, 0.6, 1.2, 1.8. K runs use a white background; NOK runs are shaded light grey. Blue lines separate the progenitor-mass regimes, as annotated in red. Orange lines enclose the pair-instability mass gap (∼64-161 M⊙, Woosley & Heger 2021). IMBHs appear below, within and above the PISN gap.
In the text
Thumbnail: Fig. A.4 Refer to the following caption and surrounding text. Fig. A.4

Lagrange radii (rLagr) and average mass (Mav) within spheres that contain 1%, 5%, 10%, 30%, 50%, and 90% of the total cluster mass at the current simulation time step for up to 500 Myr. Time is shown on a logarithmic scale to highlight the cluster’s rapid early evolution. Each column represents one rotational parameter ω0 of the rotating King model (ω0 = 0.0, 0.6, 1.2, 1.8). Models with K are plotted on a white background; models without GR kicks (NOK) are shaded light grey.
In the text
Thumbnail: Fig. A.5 Refer to the following caption and surrounding text. Fig. A.5

Time evolution of cumulative number of coalescences ncoal (top two rows) and collisions ncoll (bottom two rows) for the eight runs. Columns give increasing ω0 from left to right. K models (GR kicks on) use a white background; NOK models are shaded light grey. The horizontal axis is logarithmic in time (Myr) to clarify the pre-core-collapse phase. Colours indicate outcomes: MS (blue), CHeB (purple), ShHeB (orange), BH (black); collisions BH-BH (black) and MS-MS (cyan). Counts are in Table B.2
In the text
Thumbnail: Fig. A.6 Refer to the following caption and surrounding text. Fig. A.6

Cumulative counts of escaping binary members nesc,binarymembers and of binary components: nescMS,b, nescCHeB,b, nescShHeB,b, and nescBH,b. K: solid curves; NOK: dash-dotted curves.
In the text
Thumbnail: Fig. A.7 Refer to the following caption and surrounding text. Fig. A.7

Top two rows: Total mass of the four mass groups (MICCSNE in red, MPPISNE in green, MPISNe in orange, MhCCSNe in blue). Bottom two rows: Each group’s angular-momentum squared, Lgroup2Mathematical equation: $L_{\mathrm{group}}^2$, normalized by Lω0=0.6,t=02Mathematical equation: $L^2_{\omega_{0}=0.6,t=0}$ (i.e. Lgroup2/Lω0=0.6,t=02Mathematical equation: $L_{\mathrm{group}}^2/L^2_{\omega_{0}=0.6,t=0}$). Columns show increasing ω0 from left to right (ω0 = 0.0, 0.6, 1.2, 1.8). K runs are on a white background; NOK runs are shaded light grey.
In the text
Thumbnail: Fig. A.8 Refer to the following caption and surrounding text. Fig. A.8

All coalescence (top two rows) and all collision events (bottom two rows) from the eight N-body simulations. Plotted are product masses (mcoal or mcoll) versus distance to the cluster density centre rdens. Columns from left to right show increasing rotation ω0 = 0.0, 0.6, 1.2, 1.8. K models use a white background; NOK models are shaded light grey. Coalescence outcomes: MS (blue), CHeB (purple), ShHeB (orange), BH (black). Collision markers: BH-BH (black cross), MS-MS (cyan cross), BH-MS (violet cross). Counts in Table B.2.
In the text
Thumbnail: Fig. A.9 Refer to the following caption and surrounding text. Fig. A.9

Primary vs secondary mass (m1 > m2) for collisions from all eight simulations: MS-MS (orange), BH-MS (purple), BH-BH (black). Columns show ω0 = 0.0, 0.6, 1.2, 1.8 left-to-right. K models are shaded white; NOK models are shaded grey. The red line denotes equal mass ratio (q = 1). Counts in Table B.2.
In the text
Thumbnail: Fig. A.10 Refer to the following caption and surrounding text. Fig. A.10

GRK velocity, vGRk, versus mass ratio, q, of the participating compact objects. All but one event are BHBH mergers; the exception is a BHNS merger, marked by an extreme q and low vGRk. The shaded part of the plot corresponds to the interval 26.9 - 159.4 km s−1, the escape speeds range in cluster cores among all models.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.