Abstract

Recent observations indicate that the progenitors of globular clusters (GCs) at high redshifts had high average stellar surface densities above $10^5$ M$_\odot$pc$^{-2}$. The internal structure and kinematics of the clusters, however, remain out of reach. Numerical simulations are necessary to decipher the origin of spatio-kinematic features in present-day GCs. Here we study star cluster formation in a star-by-star hydrodynamical simulation of a low-metallicity starburst in a merger of two gas-rich dwarf galaxies. The simulation accounts for the multiphase interstellar medium, stellar radiation, winds and supernovae, and the accurate small-scale gravitational dynamics near massive stars. We also include prescriptions for stellar collisions and tidal disruption events by black holes. Gravitationally bound star clusters up to $\sim2\times10^5$ M$_\odot$ form dense with initial half-mass radii of $\sim0.1$–1 pc. The most massive cluster approaches the observed high-redshift surface densities throughout its hierarchical and dissipative assembly. The cluster also hosts a collisionally growing very massive star of $\sim1000$ M$_\odot$ that will eventually collapse, forming an intermediate mass black hole. The assembly leaves an imprint in the spatio-kinematic structure of the cluster. The youngest stars are more centrally concentrated, they show significant bulk rotation and have radially biased velocity components at outer radii. The older population is more round in shape, rotates slowly, its velocity distribution is isotropic and exhibits higher dispersion. If chemically enriched star formation proceeds mainly in the later stages of cluster assembly, these results provide a possible explanation for some of the multiple population features observed in dynamically young GCs.

galaxies: dwarf – galaxies: star clusters: general – galaxies: star formation – gravitation – methods: numerical – stars: massive

1 Introduction↩︎

Young massive star clusters are among the densest gravitationally bound stellar structures [1]–[3]. The most extreme examples are super star clusters in local starburst galaxies such as the M82 [4], [5], and young massive star clusters observed through gravitational lensing in the high-redshift Universe, such as in the Sunrise Arc and the Cosmic Gems [6], [7]. In both cases, the highest resolved effective surface densities ($\Sigma_\mathrm{eff}$) within the effective radius ($R_\mathrm{eff}$; projected half-mass or half-light radius) exceed $10^5$ $\mathrm{M}_\odot$ pc$^{-2}$. The smallest $R_\mathrm{eff}$ of these massive ($M_\mathrm{cluster}$ $>10^5$ $\mathrm{M}_\odot$) clusters are less than 1 pc. Runaway stellar collisions and black hole growth [8]–[10] as well as chemical self-enrichment through gas retention [11], [12] have been suggested to take place at such high densities. However, the details of the internal structure of the densest clusters still remain unresolved.

In nearby systems, where the internal structure of clusters can be resolved down to individual stars, clusters that are both young and massive are rare. The R136 cluster in the Large Magellanic Cloud (LMC) is one of the best examples with an age of $\sim 1.5$ Myr [13]. The total mass of the cluster has been estimated as $>2\times 10^4$ $\mathrm{M}_\odot$ within the central few pc or up to $\sim 10^5$ $\mathrm{M}_\odot$ in total [14]. The $\Sigma_\mathrm{eff}$ of R136 is, however, only a few hundred $\mathrm{M}_\odot$ pc$^{-2}$. Other notable young massive clusters in the Local Group are NGC 3603 [15], Westerlund 1 [16], Westerlund 2 [17] and the Arches cluster [18] which have mostly higher metallicities than R136 ($Z\sim$ $0.5$ $\mathrm{Z}_\odot$). None of these clusters reach $\Sigma_\mathrm{eff}$ as high as the starburst or high-redshift massive clusters.

Present-day globular clusters (GCs) in the Local Group can also be used to infer some details of the initial properties of their massive progenitors. This is however not straightforward, as GCs undergo mass segregation, mass loss and angular momentum transport (in case they rotate), through stellar evolution, relaxation and tidal effects during the several Gyrs of evolution. While the two-body relaxation time can be short in the central region of an initially dense cluster, the more sparse outer parts can still retain some of the spatio-kinematic information that was imprinted during the assembly of the cluster [19]–[21]. [22], [23] and [24] for instance find that a significant number of GCs rotate. The rotation speed over velocity dispersion ($V/\sigma$) is observed to correlate with the mass and the relaxation time of the clusters. If the less dynamically evolved systems rotate more, then the angular momentum may indeed be primordial rather than the result of, for instance, later tidal interactions. Hydrodynamical simulations of cluster formation (e.g. [25]) support these findings, showing that increasingly massive clusters inherit a larger specific angular momentum from the progenitor gas cloud [26]. Clusters forming in filamentary complexes can form with anisotropic velocity distributions, for instance showing radial bias in the outer parts [27].

The spatial and kinematic differences of the stellar populations in present-day GCs provide additional information of the initial stellar distribution. Multiple populations (MPs) are identified through stellar light-element variations that are thought to originate in the early stages of the GC evolution and may thus carry information of the GC formation process ([28]–[33]; see e.g. [34]–[38] for reviews). [39] found that the light-element enhanced second population (P2) in dynamically young GCs tends to be more centrally concentrated than the primordial, field-like population (P1). [40] confirmed the possibility of centrally concentrated P2 but also found a number of clusters for which P1 is instead centrally concentrated. The majority of their clusters, however, showed fully mixed populations. In addition, [41] found that P2 has predominantly a higher $V/\sigma$ compared to the P1 in the majority of their analysed clusters. This scenario is supported by hydrodynamical simulations of the formation of the P2 through gas accretion and/or winds of evolved intermediate mass stars in a pre-existing cluster of P1 stars [42]–[44]. [23] and [24], on the other hand, only found one or two such clusters in their samples of 25 and 30 GCs and found no differences in the kinematics of P1 and P2 in the majority of their GCs. Finally, the P1 and P2 sometimes exhibit different velocity anisotropies [33], [45], [46]. N-body simulations indicate that primordial spatial or kinematic differences may be lost during the evolution of the GC in a galactic tidal field except for the most dynamically young clusters [47], and the rate of this relaxation depends on the initial velocity and anisotropy properties [48].

A self-consistent numerical approach to the initial properties of the populations of stars in young massive clusters forming in a realistic galactic environment is therefore needed. An increasing variety of numerical simulations is now able to probe a wide range of interstellar environments that give rise to the formation of massive star clusters [49]–[62]. Both observations and simulations agree that galactic systems exhibiting high interstellar gas pressures and high star formation densities harbour the most massive star clusters [51], [63], [64]. However, because dense star clusters are inherently collisional systems, models of their formation and internal evolution should account for stellar two and few body interactions (e.g. [65], [66]). A number of recent hydrodynamical studies have explored star cluster formation and/or evolution at such high fidelity: [67], [68], [69], [70], [71] and [72] modelled star cluster formation in individual molecular clouds at very high spatial resolution and including models for stellar evolution. [70], [71] and [72] showed that massive star clusters can exhibit very high densities and star formation efficiencies, in agreement or even exceeding the high-$z$ surface densities, when the clusters are still embedded. Cloud-scale simulations are, however, limited to a low number of clusters. They do not consider the role of the galactic environment and thus the impact of the galactic tidal field, gaseous inflows and outflows and the complex time and space varying distribution of feedback from other nearby star forming regions. Toward this end, [73] introduced a method for modelling the dynamics of stars in a galaxy-scale hydrodynamical simulation but did not include a realistic stellar mass function or star formation.

In [74] we presented the first star-by-star hydrodynamical simulations of galaxy-scale star cluster formation that accounted for collisional stellar dynamics in gravitational interactions that occur near massive stars. We used the ketju integration module [75], [76] incorporated in the sphgal galaxy evolution code [77]–[82]. These simulations followed the formation of relatively low mass (up to $\sim 1000$ $\mathrm{M}_\odot$) and low density ($\Sigma_\mathrm{eff}\lesssim$ $10^4$ $\mathrm{M}_\odot$ pc$^{-2}$) clusters in a quiescent low-metallicity dwarf galaxy environment. The clusters formed with compact, sub-parsec half-mass radii and rapidly expanded onto the observed mass-size relation over the first 10 Myr. In the present study we update and apply the numerical methods presented in [74] in an extreme starburst environment of a gas-rich dwarf galaxy merger. We investigate the initial structure and kinematics of massive star clusters star-by-star in a self-consistent low-metallicity galactic environment.

The paper is structured as follows. Section 2 gives a brief overview of the hydrodynamical simulation code, newly implemented methods, initial conditions and post-processing techniques. Section 3 describes the results of the starburst simulation, including a description of the mass-size and size-surface density distribution of star clusters. We also detail the stellar collisions that occur in the cores of the dense clusters, and the growth of the most massive star in the simulation. We describe the internal density and kinematic properties of the most massive cluster and its stellar populations in Section 3.2. Section 4 provides conclusions and a discussion of the main results.

2 Simulations↩︎

The starburst simulation presented here is a part of the Galaxy Realizations Including Feedback From INdividual massive stars² (griffin) project. Details of the numerical methods have previously been described in [74], [78]–[83] and we provide a brief summary of the main features of our code here.

2.1 Hydrodynamical simulations↩︎

We use the smoothed particle hydrodynamics (SPH) code sphgal [78]–[80] that is an updated version of the gadget-3 code [77] with improvements to the numerical and algorithmic accuracy of the hydrodynamical methods. The code solves the non-equilibrium cooling and heating processes at low temperatures ($<3\times 10^4$ K) using a chemical network (H$_2$, H$^+$, H, CO, C$^+$, O and free electrons) following the methods of [84], [85] and [86]. At higher temperatures tabulated metallicity dependent cooling rates from [87] are used.

2.1.1 Star formation↩︎

Star formation is triggered when the gas conditions exceed a threshold given by the local Jeans mass $M_\mathrm{J}$ set as $M_\mathrm{J}<0.5 M_\mathrm{SPH}$ where $m_\mathrm{SPH}$ is the SPH kernel mass ($\sim 400$ $\mathrm{M}_\odot$) and \[M_\mathrm{J} = \frac{\pi^{5/2}c_\mathrm{s}^3}{6G^{3/2}\rho^{1/2}},\] $c_\mathrm{s}$ being the sound speed, $\rho$ the gas density and $G$ is the gravitational constant. Any particle exceeding the threshold is instantaneously turned into a star particle with a mass and metallicity of the gas particle. The new star particles then wait for the local dynamical time scale $t_\mathrm{dyn}=(4\pi G \rho)^{-1/2}$ after which they are sampled into individual stars along the [88] IMF between 0.08 $\mathrm{M}_\odot$ and 500 $\mathrm{M}_\odot$.

Conservation of mass is enforced in the sampling process within the Jeans length \[R_\mathrm{J} = \left( \frac{3}{4\pi} \frac{M_\mathrm{J}}{\rho} \right) ^{1/3},\] measured when the particle was turned into a star particle. The mass of the star particle (fiducial gas mass resolution $\sim 4$ $\mathrm{M}_\odot$) is sampled into several stars and the last sampled star exceeding the particle mass is taken from other such particles waiting to be sampled within $R_\mathrm{J}$. This ensures that very massive stars (VMSs; $\gtrsim100$ $\mathrm{M}_\odot$) can only form in relatively dense star-forming regions as demonstrated in [81] and [89]. The most massive star formed this way reaches $\sim330$ $\mathrm{M}_\odot$, and more massive stars therefore form only through stellar collisions (Section 2.1.4). New stars are placed around the parent particle in a Gaussian position and velocity distribution with standard deviation equal to 0.1 pc and 0.1 km s$^{-1}$ respectively in each cartesian direction. Each star is initialised as a single star, and we do not consider primordial binaries in this study. To prevent the creation of tightly bound initial binaries, we resample positions and velocities of stars that are bound to another new star within a semimajor axis of less than $5\times 10^{-3}$ pc. We also check at the end of the routine that no star is placed closer than $10^{-3}$ pc to another new star, and shift the positions of such stars by an additional Gaussian 0.01 pc. Random initialisation of such binaries is however negligibly rare.

2.1.2 Stellar feedback↩︎

We consider mass and energy release through radiation, stellar winds and supernovae. Events that we consider to occur instantaneously at the end of the stellar lifetime are core-collapse supernovae (zero age main sequence mass $8$–$40$ $\mathrm{M}_\odot$, [90]), pair-instability supernovae ($107.2$–$203.4$ $\mathrm{M}_\odot$; according to the exploding helium core masses from [91]) and asymptotic giant branch winds ($0.5$–$8$ $\mathrm{M}_\odot$, [92]). The stellar lifetimes are estimated according to the Geneva stellar models [93] for low-mass stars ($<9$ $\mathrm{M}_\odot$) and using the Bonn Optimized Stellar Tracks (BoOST; [94]) for massive stars ($9$–$500$ $\mathrm{M}_\odot$). Winds of massive stars ($>9$ $\mathrm{M}_\odot$) are treated with a continuous momentum-conserving approach. The wind mass loss rates and velocities depend on the initial stellar mass and age, interpolated using the BoOSTmodels at metallicity $Z=$ $0.01$ $\mathrm{Z}_\odot$. We further describe in Section 2.1.4 our treatment of stellar mergers and stellar masses that exceed the upper mass limit of the tables. Chemical enrichment is followed by tracking the mass of 13 chemical elements per particle: H, He, N, C, O, Si, Al, Na, Mg, Fe, S, Ca and Ne.

Each star contributes to the interstellar radiation field in the far-ultraviolet wavelength range of $6$–$13.6$ eV. The fixed zero age main sequence luminosities for stars below $9$ $\mathrm{M}_\odot$ are integrated from the BaSeL spectral library at $Z\sim 0.01$ $\mathrm{Z}_\odot$ [95], 1997A&AS..125..229L?, 1998A&AS..130...65L? using the geneva stellar models at $Z=0.0004\sim 0.02$ $\mathrm{Z}_\odot$[96]³. Massive stars ($>9$ $\mathrm{M}_\odot$) are given initial mass and age dependent rates of far-ultraviolet and hydrogen ionizing photons ($>13.6$ eV) from the BoOST tables.

The strength of the interstellar radiation field is computed at the location of each gas particle using the treecol-algorithm [97]. The incoming flux is attenuated by gas and dust column density assuming optically thin conditions in 12 solid angles divided using the healpix algorithm [98]. For injection of supernovae and winds we likewise use healpix to distribute the $1/12$th fraction of mass and energy equally into $8\pm2$ closest gas particles in each solid angle. Photoionization is implemented as an approximation wherein the photoionized gas within the Strömgren sphere is set fully ionized at $10^4$ K and the ionization equilibrium is solved iteratively when HII regions overlap. More details for the feedback implementation can be found in [79], [80] and [81].

2.1.3 Collisional stellar gravitational dynamics↩︎

In [82] and [74] we described how to combine sphgal with the ketju dynamics module [75], [76]. All stellar gravitational interactions within a radius $r_\mathrm{KETJU}$ of any star or a remnant above a chosen initial mass limit $m_i$ are treated with the algorithmically regularised integrator mstar [99]. The regularised integration in the ketju regions enables accurate calculation of the stellar trajectories without the use of gravitational softening. The accuracy parameters in ketju are the same as in [74]: Gragg-Bulirsch-Stoer tolerance is set to $\eta_\mathrm{GBS}=10^{-7}$ and the end-time iteration tolerance is $\eta_\mathrm{t}=10^{-3}$.

The computational cost of ketju scales with both the number of regions (linear scaling) and the number of particles in the regions (near quadratic scaling), therefore both $r_\mathrm{KETJU}$ and $m_i$ have to be selected optimally for the specific problem. The size of the ketju-integration region has to be set as at least $2.8$ times the gravitational softening length ($\epsilon_*$) used in the standard tree force calculation of gadget-3. Here we use $\epsilon_*=0.01$ pc and therefore $r_\mathrm{KETJU}=0.03$ pc. This value of $\epsilon_*$ is selected to give a manageable number of up to a few hundred stars in the ketju-regions within massive clusters that can reach more than $10^6$ stars per pc$^3$ in our simulations. In practice, overlapping ketju-regions are combined into one, therefore the cores of the young massive clusters will have somewhat larger regularised regions as long as massive stars are not expelled. In the most massive cluster, up to a few $10^4$ stars are within the central ketju-region that extends to a radius of $0.1$ pc within the centre of mass.

As the present study concentrates on the most massive clusters that are dense and host a large number of massive stars, we set $m_i=8$ $\mathrm{M}_\odot$. This guarantees that the dynamics of massive stars, their remnants and their neighbours in the cores of most of the young clusters more massive than a few hundred $\mathrm{M}_\odot$ will be regularised (see the results of IMF sampling in [81]). This results in a significantly improved dynamical evolution of the stellar component in the clusters as demonstrated in [74] via comparison to direct N-body simulations performed with the bifrost code [100].

2.1.4 Stellar collisions↩︎

At high stellar densities, stellar collisions can lead to the growth of massive objects in the cores of the clusters (see e.g. N-body and cloud-scale simulations of [10], [70], [101]–[103]). Following stellar evolution combined with interactions and mergers in a sub-resolution model would require understanding how such mergers proceed (e.g. mixing and mass loss, [104], [105]) and how the remnants of the interactions evolve [106]–[108]. In the absence of comprehensive models for stellar evolution under repeated collisions we thus implement a mass conserving stellar collision approximation following ideas that are often applied in rapid binary population synthesis (e.g. [109]; see [110] and [111] for modern approaches).

We assume that a collision occurs when the radii of two stars overlap. N-body simulations indicate that there is usually only one massive object growing in the cluster core through collisions and that the majority of the mass gain comes in the form of other massive stars [103]. We use the stellar radii provided in the BoOST models ($m_i>9$ $\mathrm{M}_\odot$) to track the collisions, therefore collisions can only occur when one of the progenitors is a massive star. When a collision with another star occurs, the current mass of the less massive star is added to the more massive one and the age over lifetime ratio $f_\mathrm{age}=\mathrm{age}/\mathrm{lifetime}$ for the more massive star is recorded. We then go through BoOST stellar models with increasing initial mass to find an epoch that first matches the total (current) mass of the star and the fractional age $f_\mathrm{age}$ on the specific track. As a conservative estimate of the total lifetime of the star, we only update the current age of the star to match the respective fractional age on the new higher mass stellar track. We thus do not implement any explicit rejuvenation methods unlike is often done in N-body models that consider mass transfer or collisions [112], [113]. We will further refine the assignment of the new stellar age in the future, to account for the evolutionary stage (e.g. to retain the core He burning phase) using the equivalent evolutionary phases provided in the BoOST models.

The collision product is supplemented with the metals of the donor star. After the collision, the collision product will follow the newly assigned stellar track until another collision occurs or the star dies.

To our knowledge, at the time of starting the simulation there were no stellar models exceeding $\sim500$–$600$ $\mathrm{M}_\odot$ that would have also provided us with the chemical properties (e.g. Al, Na abundances) on the surface of the star throughout the stellar evolution. Our choice for the current study was thus to remain conservative and use the 500 $\mathrm{M}_\odot$ BoOST model at $Z=0.01$ $\mathrm{M}_\odot$ to evolve any star that exceeds 500 $\mathrm{M}_\odot$. This results in a lower limit for mass loss and the stellar radius. Recent extensions to the range of tables with initial stellar masses beyond $\sim500$ $\mathrm{M}_\odot$ include the PARSEC tables that now reach 2000 $\mathrm{M}_\odot$ at low metallicities [114]. For future applications we may consider updating our models to such tables. However, given that the BoOST models are among the only ones that provide enough details of the chemical species related to GC MPs, we may alternatively explore the possibility of using tables such as provided by [114] to build informed extrapolation techniques to extend the upper limit of 500 $\mathrm{M}_\odot$ using the BoOST models.

2.1.5 Tidal disruption events↩︎

A simple prescription for tidal disruption of massive stars by stellar black holes is also included. We detect a tidal disruption event (TDE) using the tidal radius given by \[r_\mathrm{tidal} = 1.3 r_\mathrm{*} \left(\frac{m_\mathrm{BH}+m_\mathrm{*}}{2m_\mathrm{*}}\right)^{1/3}\] where $r_\mathrm{*}$ and $m_\mathrm{*}$ are the current radius and mass of the star and $m_\mathrm{BH}$ is the mass of the black hole [115]. Estimates of the exact amount of stellar debris accreted by the black hole depend on the internal structure of the star and the collision parameters and are still uncertain [116], ranging from a few per cent (e.g. [117]) to unity [118]. We follow the simple procedure often adopted in N-body and Monte Carlo simulations and deposit a fixed fraction (e.g. [119], [120]) of 50% of the stellar mass into the black hole [121] and remove the other half from the simulation. The model can be further refined in the future to inject the returned mass back into the ISM. TDEs are extremely rare compared to stellar collisions in the current model since we only run the most intense starburst-phase for a short while. We will therefore leave exploration of TDEs as a source of chemical enrichment for future work.

2.2 Initial conditions and simulation parameters↩︎

The initial setup is based on a previous dwarf galaxy merger simulation presented in [89]. We re-simulate the starburst phase with the updated code version including collisional dynamics starting at a time when the starburst associated with the first passage is just starting to ramp up in the original simulation (simulation time 60 Myr in Fig. 4 in [89]). Some small clusters have already formed by this time in the disks of the galaxies, however the main burst occurs in the overlap region between the galaxies and remains spatially separated from the pre-existing stellar disks for the duration of the burst. Those clusters that formed before the burst in the fully gravitationally softened simulation are here left out of the analysis.

The dwarf galaxies have virial and baryon masses of $4\times 10^{10}$ $\mathrm{M}_\odot$ and $6\times 10^{7}$ $\mathrm{M}_\odot$, respectively, set up with the methods described in [77]. The orbital configuration is the same as in [51] but with reduced approach velocity (factor of 4) to achieve the starburst phase already after the first passage between the galaxies. The galaxies have initial gaseous and inert stellar disks with 66% gas fraction and 0.73 kpc scale lengths. The initial gas and stellar particle resolution is 4 $\mathrm{M}_\odot$ and the dark matter particle resolution is $6.8\times10^3$ $\mathrm{M}_\odot$. The initial metallicity is $\sim0.016$ $\mathrm{Z}_\odot$.

Gravitational interactions of dark matter, gas and pre-existing stars are softened, and the time integration is performed using the second order leapfrog integrator in gadget-3. The gravitational softening length is set as 62 pc for dark matter and 0.1 pc for gas and pre-existing stars. Stars formed during the simulation have a gravitational softening length of 0.01 pc if they are not within a ketju integration region (see Section 2.1.3). Only the gravitational interactions of new stars and their remnants that occur within $r_\mathrm{KETJU}=0.03$ pc of any star with $m_i>8$ $\mathrm{M}_\odot$ are solved using ketju without gravitational softening.

We continue the starburst simulation for a total duration of $\sim7.5$ Myr until the largest stellar collision product exceeds $\sim 1000$ $\mathrm{M}_\odot$, i.e. when the most massive star considerably exceeds the upper limit of stellar evolution models ($500$ $\mathrm{M}_\odot$). By this time, the peak star formation rate of the starburst has passed, the most massive cluster has formed most of its stars, and the majority of the star clusters formed in the simulation have become gas free.

Figure 1: The mass-size (left) and size-surface density (right) distributions of the simulated and observed young star clusters. The pink open circles show the general young star cluster population in the final snapshot. The squares indicate the evolution of the most massive cluster in \sim0.1 Myr steps starting when it had a bound mass equal to 25\% of its final mass (see Fig. 2), coloured according to its mean stellar age from light to dark shade. The comparison data are observed young star clusters (<10 Myr) in the LEGUS survey (grey hexbins with the shade from light to dark indicating the number of clusters ranging from 1 to 38 clusters on the left and 1 to 28 clusters on the right; [122]) and various surveys of strongly lensed clusters at redshifts z\gtrsim2.4 (filled symbols; [6], [7], [123]–[126]), for which we have used the reported half-light radii to compute the effective surface densities. The solid line shows the compactness parameter C_5=(M_\mathrm{cluster}/10^5\, \mathrm{M}_\odot)/(R_\mathrm{eff}/\mathrm{pc})=1 [11]. The dotted lines on the left show sizes corresponding to a constant surface density for masses in the range 10^3–10^6 \mathrm{M}_\odot pc^{-2} from left to right and the dashed lines on the right show the surface densities for a constant mass in the range 10^3–10^6 \mathrm{M}_\odot from left to right, in steps of 1 dex. — Figure 1: The mass-size (*left*) and size-surface density (*right*) distributions of the simulated and observed young star clusters. The pink open circles show the general young star cluster population in the final snapshot. The squares indicate the evolution of the most massive cluster in $\sim0.1$ Myr steps starting when it had a bound mass equal to $25\%$ of its final mass (see Fig. 2), coloured according to its mean stellar age from light to dark shade. The comparison data are observed young star clusters ($<10$ Myr) in the LEGUS survey (grey hexbins with the shade from light to dark indicating the number of clusters ranging from 1 to 38 clusters on the left and 1 to 28 clusters on the right; [122]) and various surveys of strongly lensed clusters at redshifts $z\gtrsim2.4$ (filled symbols; [6], [7], [123]–[126]), for which we have used the reported half-light radii to compute the effective surface densities. The solid line shows the compactness parameter $C_5=(M_\mathrm{cluster}/10^5\, \mathrm{M}_\odot)/(R_\mathrm{eff}/\mathrm{pc})=1$ [11]. The dotted lines on the left show sizes corresponding to a constant surface density for masses in the range $10^3$–$10^6$ $\mathrm{M}_\odot$ pc$^{-2}$ from left to right and the dashed lines on the right show the surface densities for a constant mass in the range $10^3$–$10^6$ $\mathrm{M}_\odot$ from left to right, in steps of 1 dex.

2.3 Cluster selection↩︎

We identify gravitationally bound star clusters using the friends-of-friends (fof) and subfind [127], [128] algorithms in gadget-3. The linking length of fof is set as 0.2 pc to prevent the algorithm from including extended sparse structures outside of the concentrated clusters of stars. A bound structure has to have at least 50 members after the unbinding step of subfind to be considered a cluster. After the selection of bound clusters we also include any stars within the half-mass radius regardless of their velocity, to catch stars that are in binaries or escaping but have not yet exited the central region of the cluster. These stars are included in the mass, light and density calculations but not in the computations of the velocity distribution.

2.4 Radiative transfer post-processing↩︎

We produce photometric images of the most massive cluster using the dusty Monte Carlo radiative transfer code skirt 9 [129]. We follow the methods for setting up the input radiation sources and the construction of the dust grid described in [81], [82], [130]. Briefly, all stars more massive than 0.8 $\mathrm{M}_\odot$ that have formed during the simulation are used as radiation sources with spectra given by atmosphere models of [131]. As with stellar evolution, stars beyond $500$ $\mathrm{M}_\odot$ are treated as stars with initial mass of $500$ $\mathrm{M}_\odot$ for consistency. The dust distribution is computed as a (60 pc)$^3$ octree grid constructed according to the SPH particle data and a fixed $\sim80\%$ dust-to-metals mass ratio. The grid is allowed to refine at most down to $\sim 0.01$ pc. Only gas below a temperature of 8000 K is considered to contain dust. Photometric images are then made at $\sim0.04$ arcsec per pixel resolution that is close to the pixel resolution of Hubble Space Telescope (HST) and James Webb Space Telescope (JWST), placing the cluster at 50 kpc distance that is similar to the distance of R136 in LMC. We use three standard broadband Johnson-Cousins filters B, V and I built in skirt. We apply a Gaussian point spread function with a full width at half maximum of 2 pixels.

3 Results↩︎

3.1 Properties of young star clusters↩︎

Approximately 850 star clusters that fill our selection criteria form during the $\sim 7.5$ Myr of starburst. Star clusters above $\sim 300$ $\mathrm{M}_\odot$ exhibit a cluster mass function that has a power-law index between $-1.9$ and $-2$ throughout the starburst, fit using 20 bins with equal number of 40–50 clusters per bin. The slope of the mass function agrees with observed mass functions and luminosity functions across environments from low-metallicity galaxies such as the Magellanic Clouds to starburst systems [132]–[137].

By the final snapshot, the majority of the clusters are already gas-free while all clusters have mean stellar ages less than 4.3 Myr. The most massive cluster at a mean stellar age of 1.9 Myr still has a gas-embedded core that is forming stars. While the tail end of the starburst is still ongoing, we limit the analysis only until this epoch. The majority of the stellar mass has already formed, and the main stellar collision product (Section 3.1.2) in the most massive cluster has not yet reached a mass significantly above 1000 $\mathrm{M}_\odot$, i.e. more than a factor of two beyond the upper stellar mass limit adopted in our VMS models (see Section 2.1.4). Fig. 1 shows the projected half-mass radii ($R_\mathrm{eff}$) and the corresponding surface densities $\Sigma_\mathrm{eff}=M_\mathrm{cluster}/(2\pi R_\mathrm{eff}^2)$ in the $x$–$y$ plane for the simulated gravitationally bound star clusters in the final snapshot. The results are compared to masses and half-light radii of young ($<10$ Myr) star clusters in the LEGUS survey [122] and of resolved or marginally resolved star clusters measured through gravitational lensing at redshifts $z\gtrsim 2.4$ ([6], [7], [123]–[126]). The surface densities in the observed clusters have been computed using the reported half-light radii.

The majority of the simulated star clusters show compact sizes of $R_\mathrm{eff}\sim0.1$–$1$ pc, though only the most massive cluster reaches surface densities of close to $10^5$ $\mathrm{M}_\odot$ pc$^{-2}$. The range of initial sizes is in general agreement with compact sizes found in our previous studies in [51], [89] and [74] for young clusters. Cloud-scale simulations of massive star clusters often also report compact sizes, e.g. in [71], [138]. Clusters in the new simulation that have initial masses of a few $10^2$–$10^5$ $\mathrm{M}_\odot$ are more extended on average by a factor of 1.5–2 compared to clusters of similar mass and age in the previous setup of [89] without collisional dynamics near massive stars. This difference is also seen as a slightly steeper mass-size relation slope of 0.32 compared to 0.26 fit to the clusters in [89]. Similar to the results of [74] for low-mass clusters, the fully exposed clusters undergo rapid size evolution due to internal gravitational dynamics already during the very first Myrs.

Our clusters have sizes that correspond to the most compact and dense young star clusters in local observations. Observationally it remains unclear whether very young clusters in general begin their lives with smaller sizes and higher densities compared to evolved clusters as suggested by theoretical models (e.g. [139], [140]) and our simulations. The larger sizes and lower surface densities of present-day GCs ([1]) compared to their high-redshift counterparts imply some evolution at least over extended periods of time. The size-mass relation of large samples of low-redshift star clusters e.g. in the LEGUS survey (e.g. [141], [122]; see also [142] for a discussion) instead show on average very little evolution with cluster age. Interestingly, [143] for instance verified that selection effects due to point-source removal should not be the cause for missing centrally concentrated bright massive clusters in the local Universe. On the other hand, young clusters may reach peak densities while they are still embedded, remaining mostly invisible to observations in visual wavelengths [144]. High-resolution infrared observations e.g. with the JWST may bring more light to this aspect in near future [5], [145].

In support of the latter argument, our compact clusters in Fig. 1 are very young and the majority have just become fully gas-free. As was shown in [74], the star clusters first become more compact while they are still embedded and forming stars. The most massive clusters in the present study form hierarchically through accretion of gaseous and stellar sub-clumps, and their internal structure will undergo violent relaxation [146]. The time evolution of the most massive star cluster is shown here in Fig. 1 as a specific example, starting when the most massive progenitor cluster had a bound mass of 25% of the total mass of the final cluster. The contraction of the cluster can be seen in the time-sequence in Fig. 1. Similar contraction was also shown in cloud-scale simulations of massive star cluster formation by [71].

Finally, we also show the compactness parameter $C_5=(M_\mathrm{cluster}/10^5\, \mathrm{M}_\odot)/(R_\mathrm{eff}/\mathrm{pc})$ with a value of $C_5=1$ in Fig. 1. This was estimated in [11] to be the approximate limit above which gas expulsion is no longer possible due to failure of stellar feedback (see also e.g. [147], [148] and [36]). Our most massive star cluster begins its life close to this relation, and progresses toward larger values of $C_5$ as it grows in mass and contracts. All this happens while VMSs are forming and growing in its core and releasing large amounts of pre-SN feedback. The $C_5=1$ line divides the mass-size and size-surface density planes into two regimes, where our simulated most massive cluster falls in the same regime with the high-redshift dense clusters. Our lower mass simulated clusters and the local observed young star clusters populate the $C_5<1$ regime. As suggested in [11] and supported by the simulations in [89], $C_5>1$ may well be the dividing line above which self-enrichment and MPs start occurring. Overall, the nitrogen enhancements observed in some high-redshift compact galaxies (e.g. [149]–[152]), the high compactness of the corresponding high-redshift clusters and their high stellar masses are in line with the approximate initial GC mass limit of $\gtrsim 10^5$ $\mathrm{M}_\odot$ for Local Group GCs to host MPs [37].

3.1.1 Stellar collisions and two TDEs↩︎

In total $\sim12500$ stars more massive than $8$ $\mathrm{M}_\odot$ form throughout the starburst. 69 are more massive than $100$ $\mathrm{M}_\odot$. When stellar collisions occur in our simulation, only one star grows per cluster after segregating into the cluster core, similar to results in the literature [103], [113], [153], [154]. In two cases the merger product of two stars below the limit of 100 $\mathrm{M}_\odot$ exceeds the VMS limit. The rest of the events that grow VMSs start with a star originally above 100 $\mathrm{M}_\odot$. The most massive star cluster hosts $>2000$ stars more massive than $8$ $\mathrm{M}_\odot$, up to $\sim18$ of which are more massive than $100$ $\mathrm{M}_\odot$ and present in the cluster at the same time. 12 remain in the cluster in the final snapshot. 13 of the VMSs form after the main progenitor has reached 25% of its final mass and they form within $\sim1$ pc of the cluster centre. The rest arrive via sub-cluster mergers. By the last snapshot, $78$ stars in the entire galaxy more massive than $40$ $\mathrm{M}_\odot$ have collapsed as black holes, out of which $8$ remain in the most massive cluster. 20 stars have exploded as pair-instability supernovae.

For comparison, R136 includes tens of stars that are more massive than a few tens of $\mathrm{M}_\odot$. Of these $\sim10$ may have started initially with masses of more than 100 $\mathrm{M}_\odot$ [155], the most massive estimated to have a current mass of $\sim220$ $\mathrm{M}_\odot$ and an initial mass of $\gtrsim250$ $\mathrm{M}_\odot$. The most massive stars in R136 are currently located in the central $\sim 0.5$ pc of the cluster. The number of VMSs in our most massive cluster is thus approximately twice as many as in the R136, consistent with the higher total cluster mass.

We track the collisions of stars throughout the simulation for pairs that include a star or a black hole that had an initial stellar mass of at least $9$ $\mathrm{M}_\odot$. During the $\sim7$ Myr of starburst there are in total $90$ collisions throughout the simulation, out of which $16$ involve two components both of which are more massive than $40$ $\mathrm{M}_\odot$. Almost $30\%$ of the collisions are between fairly equal mass objects (mass ratio $>1/3$). We note that compared to the hierarchical direct N-body simulations of [103], we do not see the growth of VMSs through collisions in the in-falling sub-clusters during the assembly of the most massive cluster. This may be due to somewhat delayed mass segregation and sub-cluster core collapse in our collisional dynamics implementation, compared to full direct N-body (see Fig. 1 in [74]). The core collapse time scale, for instance, affects when VMS growth can occur in the sub-clusters [154]. Optionally the spread in stellar ages in each cluster due to the star formation history may delay the birth of VMSs and their expansion into the giant phase compared to the uniform age initial conditions adopted in N-body models. This shifts the collisional growth of VMSs to the later, more centrally concentrated stage of cluster assembly.

There are two (micro-)TDEs: one between a 58 $\mathrm{M}_\odot$ black hole and a 125 $\mathrm{M}_\odot$ VMS (age $\sim 2.3$ Myr and initial mass 132 $\mathrm{M}_\odot$), and the other between a 160 $\mathrm{M}_\odot$ black hole and a 74 $\mathrm{M}_\odot$ VMS (age $\sim 2.9$ Myr and initial mass 89 $\mathrm{M}_\odot$). In the first TDE the VMS had recently reached its supergiant phase. In the second one the star was just about to collapse as a black hole. The collisions occurred not in the most massive cluster but in the cores of a $2\times10^4$ $\mathrm{M}_\odot$ and a $3\times10^4$ $\mathrm{M}_\odot$ cluster that had evolved to a mean stellar age of 2.7 Myr and 2.9 Myr, respectively. Black hole mergers do not occur in the simulation essentially because the post-Newtonian terms in ketju were disabled for the run.

There are several factors that affect the collision rate obtained in the simulation. Because close gravitational interactions are only possible in systems that are integrated with ketju and mergers are only allowed to happen for pairs that contain a $m_i>8$ $\mathrm{M}_\odot$ star, the collision rate found here is a lower limit. On the other hand, the self-consistent star formation history that we obtain for each cluster adds complexity to the stellar evolution and collisions compared to pure N-body simulations that often assume uniform age clusters. The stars for instance evolve to their giant stage later if they form later in the cluster formation process. This will naturally reduce the collision rate during the assembly process due to the reduced collision cross-section. Meanwhile, the most massive stars mainly form in the central regions of the clusters and mass segregation of the most massive collision partners can thus happen on a shorter time scale. We omit the initialisation of primordial binaries, which would additionally enhance the collision rate [156], [157].

Figure 2: Top: Colour composite images of the central 20 pc at four epochs along the formation of the most massive cluster (see text) in B (blue), V (green) and I (red) Johnson-Cousins broad-bands. The image resolution corresponds to that of the HST at the distance of LMC (\sim 50 kpc, 0.01 pc per pixel) degraded with a \sim2 pix Gaussian point spread function. The insets show the central 2 pc in an adjusted value range.Middle: V-band surface brightness maps of the stellar emission without dust attenuation. The white arrows show the centres of light (see text for details). The brightest spot close to the centre of the rightmost panel corresponds to the most massive star (see Section 3.1.2) whose radiation properties are modelled as an initially 500 \mathrm{M}_\odot star at an age of 1.2 Myr. This star, among others, is not visible in the top and bottom rows due to extinction. The insets in the middle and bottom rows show the central 2 pc in the same value range.Bottom: V-band surface brightness maps with attenuation, scattering and dust emission. The red arrows indicate the centres of light of the attenuated images measured similarly to the unattenuated one to highlight the offset in the direct and attenuated light. — Figure 2: *Top*: Colour composite images of the central 20 pc at four epochs along the formation of the most massive cluster (see text) in B (blue), V (green) and I (red) Johnson-Cousins broad-bands. The image resolution corresponds to that of the HST at the distance of LMC ($\sim 50$ kpc, 0.01 pc per pixel) degraded with a $\sim2$ pix Gaussian point spread function. The insets show the central 2 pc in an adjusted value range.Middle: V-band surface brightness maps of the stellar emission without dust attenuation. The white arrows show the centres of light (see text for details). The brightest spot close to the centre of the rightmost panel corresponds to the most massive star (see Section 3.1.2) whose radiation properties are modelled as an initially 500 $\mathrm{M}_\odot$ star at an age of 1.2 Myr. This star, among others, is not visible in the top and bottom rows due to extinction. The insets in the middle and bottom rows show the central 2 pc in the same value range.Bottom: V-band surface brightness maps with attenuation, scattering and dust emission. The red arrows indicate the centres of light of the attenuated images measured similarly to the unattenuated one to highlight the offset in the direct and attenuated light.

3.1.2 The most massive star↩︎

In the most massive cluster, one star grows to a peak mass of 1046 $\mathrm{M}_\odot$ through repeated collisions with other massive stars. It began its life as a 150 $\mathrm{M}_\odot$ star, and consequently merged with stars of 115 $\mathrm{M}_\odot$, 170 $\mathrm{M}_\odot$, 191 $\mathrm{M}_\odot$, 232 $\mathrm{M}_\odot$, 153 $\mathrm{M}_\odot$ and 44 $\mathrm{M}_\odot$ during a period of 1 Myr. The collisions are spaced by 0.1 Myr, 0.15 Myr, 0.4 Myr, 0.3 Myr and 4.4 kyr. The most massive star also merged with 13 low-mass stars ($<0.4$ $\mathrm{M}_\odot$), increasingly frequently during the final 0.2 Myr as it grew to the supergiant phase. We limit the analysis to this epoch, however the star would still have the potential to collide with a variety of other massive stars, including the most massive normal (i.e. not collisionally formed) VMS of 330 $\mathrm{M}_\odot$. This star is at a distance of only 0.05 pc from the cluster core and the most massive star in the final snapshot.

The final VMS mass that far exceeds the upper stellar mass limit of the adopted IMF, and its growth being dominated by collisions with other VMSs, are consistent with hierarchical direct N-body models of star cluster assembly [103]. All of the individual collision partners are within the range of estimated initial masses of the most massive stars in R136, therefore the physical scenario for the VMS growth through collisions is very well motivated based on observations. Why no such massive stars that could be collision products are observed in R136 may be due to its too low central/effective density, or its higher metallicity (factor of tens) that allows stronger mass-loss and smaller stellar radii at evolved stages.

Our resulting growth of the most massive star is less extreme than what is found in idealised cloud-scale hydrodynamical simulations (e.g. [70], [158], [159]). Such simulations of dense star-forming clouds find VMSs in excess of $10^3$–$10^4$ $\mathrm{M}_\odot$ forming through collisions and/or gas accretion. Future work with improved treatment for the evolution of the most massive stars will help address whether this difference is a result of the galactic environment and more complex cloud structure in our model, the stellar evolution models, the adopted IMF, or perhaps something else.

Though we have refrained from making assumptions about the detailed evolution of stars beyond $500$ $\mathrm{M}_\odot$, an order of magnitude estimate for the final mass of the most massive star and its black hole remnant can be made based on simplified arguments about the mass and metallicity dependent wind mass loss e.g. by [12] and [160]. The most massive star undergoes its last major collision at an age of $\sim 1.2$ Myr, and it can be expected to still live for one million years. A $\sim1000$ $\mathrm{M}_\odot$ star at $\sim 0.01$ $\mathrm{Z}_\odot$ would initially lose mass in optically thick winds at a rate between $\sim 10^{-4.3}$–$10^{-3.7}$ $\mathrm{M}_\odot$ yr$^{-1}$ depending on the rate estimate. The star would lose in total between $\sim25$–$160$ $\mathrm{M}_\odot$ during the final 1 Myr of its life, resulting in a black hole with a mass of at least $\sim840$ $\mathrm{M}_\odot$. Even if the star would be rotating fast as a consequence of its merger history and therefore lose mass at a factor of 2–3 faster rate as argued by [161], it would still result in a final black hole mass of at least $640$ $\mathrm{M}_\odot$. Additionally, if the stellar collisions were implemented with a prescription for mass loss, we can approximate the star to end up with 10% lower maximum mass, $\sim900$ $\mathrm{M}_\odot$, after the last major collision. This would further decrease the resulting black hole mass to a minimum of 600 $\mathrm{M}_\odot$ (rotating VMS with [160] mass loss rate) or 770 $\mathrm{M}_\odot$ (without rotation). The result would therefore be an intermediate mass black hole in a proto-globular cluster-like system. Whether such a black hole would be retained in the cluster or the galaxy over several Gyrs would require a more detailed modelling of black hole–black hole interactions and mergers in the simulation (see previous studies e.g. by [162], [163], [103]).

Figure 3: The radial projected profiles of stellar mass surface density (a), V-band surface brightness (b), 3D mass density (c) and stellar number density (d) of the most massive cluster at the four epochs shown in Fig. 2 and in the last snapshot when the cluster has a mean stellar age of 1.9 Myr. The solid lines in panels a, c and d indicate the gravitationally bound stellar component and the dashed lines show the total stellar mass in the region. In panel b the solid and dashed lines show the attenuated and the unattenuated surface brightness profiles, respectively. The power-law slopes of -0.5, -1, -2 and -3 are indicated in the panels. The best-fit surface density and brightness profiles of R136 from [1] are shown in black in panels a and b, and the right-side y-axis in panel b shows the surface brightness scale in mag arcsec^{-2}. The vertical bars at the bottom of the panels show the projected half-mass radii of the bound component (panels a, c and d) and the projected half-light radii of the attenuated and unattenuated light profiles (b) with respective linestyles. — Figure 3: The radial projected profiles of stellar mass surface density (a), V-band surface brightness (b), 3D mass density (c) and stellar number density (d) of the most massive cluster at the four epochs shown in Fig. 2 and in the last snapshot when the cluster has a mean stellar age of 1.9 Myr. The solid lines in panels a, c and d indicate the gravitationally bound stellar component and the dashed lines show the total stellar mass in the region. In panel b the solid and dashed lines show the attenuated and the unattenuated surface brightness profiles, respectively. The power-law slopes of $-0.5$, $-1$, $-2$ and $-3$ are indicated in the panels. The best-fit surface density and brightness profiles of R136 from [1] are shown in black in panels a and b, and the right-side y-axis in panel b shows the surface brightness scale in mag arcsec$^{-2}$. The vertical bars at the bottom of the panels show the projected half-mass radii of the bound component (panels a, c and d) and the projected half-light radii of the attenuated and unattenuated light profiles (b) with respective linestyles.

3.2 Properties of the most massive cluster↩︎

Fig. 2 shows a time sequence of the accumulation of mass in the most massive star cluster. The panels show colour composites and surface brightness maps of the direct and attenuated stellar light at four epochs. The epochs shown correspond to when the gravitationally bound main progenitor is 25%, 50%, 75% and $\sim90\%$ of the final mass ($\sim 2.2\times 10^5$ $\mathrm{M}_\odot$) of the most massive cluster. The mean stellar age of the main progenitor cluster is $0.8$ Myr, $0.9$ Myr, $1.3$ Myr and $1.6$ Myr in snapshots that are spaced by $0.2$ Myr, $0.4$ Myr and $1.1$ Myr, the last being $\sim0.4$ Myr before the final snapshot. The visual appearance of the most massive cluster in the final snapshot does not show additional features except for having less gas and a more compact central concentration of massive stars, therefore we only show the final surface density profile later in Section 3.2.1. Throughout the assembly, the mean gas surface density within 1 pc is $\sim 3$–5 $\mathrm{M}_\odot$ pc$^{-2}$, i.e. at least an order of magnitude less than that of the stellar component.

The first two epochs show the phase of cluster assembly where stellar clumps are concurrently coalescing and forming stars. The stellar structure at the top of the panels in the two left columns is falling toward the cluster. The centre of light indicated in the V-band images is measured with a 2D shrinking sphere of minimum radius of 10 pixels, masking out all pixels brighter than $10^5$ $\mathrm{L}_\odot$ (i.e. massive stars) before applying the point spread function. This prevents the algorithm from centring on the brightest stars.

The stellar light of the main cluster is attenuated and the other massive clump in the images remains the brightest cluster in V-band for the duration of this coalescence, irrespective of the masking applied in the centring. The V-band attenuation $A_V$ is $\sim 0.9$ mag at these epochs, computed as the AB-magnitude difference between the unattenuated and attenuated images within a radius of 1 pc of the centre of light of the unattenuated image. This is a moderate value for observed star clusters and not as extreme as the most attenuated clusters observed in some dwarf galaxies (e.g. [144]). The reddening $E(B-V)$ within the same 1 pc region is 0.7 and 0.5 in the two first epochs, similar to what is typically found for young clusters in dwarf galaxies or starbursting systems ($E(B-V)$ up to 0.7 e.g. in [64], [164]).

The last two epochs in Fig. 2 show centrally concentrated star formation. The central 1 pc has a star-forming disky gas component until the last snapshot. The gas disk is being fed along a gaseous filament that extends toward the cluster from north-west in the rightmost panels. $A_V$ peaks at $\sim 1.4$ during the last epochs, and $E(B-V)$ drops to 0.3–0.4.

3.2.1 Density, surface density and surface brightness profiles↩︎

In Fig. 3 we show the build up of the radial density, surface density and V-band light profile of the most massive cluster at epochs corresponding to Fig. 2 as well as in the final snapshot. The particle data is centred with a shrinking sphere ignoring all stars more massive than 50 $\mathrm{M}_\odot$, similar to the brightness profiles.

The density and brightness profiles of the most massive cluster are consistent with power-law shapes especially in the outer parts and at later epochs when the stellar distribution is less clumpy (see Fig. 2). The infalling clusters in our simulation can be seen as an increase in density at intermediate radii ($1$–$3$ pc) at earlier epochs. The outer radial power-law slope in surface density is $-2.2$, as expected from models of hierarchical sub-clusters assembly that predict slopes that approach $-2$ [103], [165], [166]. The outer power-law has a respective density slope of slightly shallower than $-3$ and continues well into the sub-parsec regime. The final central 3D mass density, 3D number density and surface mass density (inner $\sim1000$ stars) are a factor of $\sim190$, 10 and 10 times higher than the half-mass values shown in Fig. 1, respectively. The unattenuated central surface brightness profile likewise has a steep inner slope, dominated by the most massive stars in the cluster core. Based on these simulated brightness and density profiles, we therefore infer that the central surface densities $\Sigma_0$ of the observed high-redshift clusters may be at least an order of magnitude higher than the effective values $\Sigma_\mathrm{eff}$, i.e. $\Sigma_0>10^6$ $\mathrm{M}_\odot$ pc$^{-2}$. This may be even before the clusters undergo core collapse, based on our simulated cluster that has only just formed. The fact that the cluster is still forming stars in a central disk is consistent with the idea that the central stellar density of $>5\times 10^5$ $\mathrm{M}_\odot$ pc$^{-2}$ is high enough to prevent pre-SN stellar feedback from expelling the gas [148].

The number density profile shows a shallow central profile that is an expected feature in N-body simulations of repeated star cluster mergers [103], [167]. The relaxation time of the core is approximately one Myr, therefore star formation and relaxation occur simultaneously, together with mass segregation. The spatially resolved scale of the stellar component in the simulation is a few $0.01$ pc or less, however we note that the star sampling routine places stars in a Gaussian kernel of $\sigma=0.1$ pc that impacts the central-most density profile of the youngest stars. We therefore refrain from making quantitative statements regarding the density profile in the core region.

We also show a comparison to the best-fit surface density and V-band surface brightness profiles of the young massive cluster R136. The inner radius is limited to the range of data typically used in the fits such as provided by [168] and [1]. As with young clusters in general [1], the radial surface brightness profile of R136 is well fit with an Elson-Fall-Freeman profile [169]. The profile of R136 has an outer power-law slope of the order of $-2$ within the central $\sim 5$ pc [1], [168] and possibly a small central core of $\sim0.3$ pc. The central R136 cluster is less massive than our simulated cluster by a factor of a few, nonetheless its surface density and light profiles are qualitatively very similar to our simulated cluster. The central mass density of R136 has been estimated to be between $10^4$–$10^7$ $\mathrm{M}_\odot$ pc$^{-3}$ [170], depending on the assumed size of the central core. Such values are supported by our hydrodynamical results throughout the cluster assembly.

$Figure 4: The radial velocity distribution of the most massive cluster in 0.1 Myr steps. The lines are coloured with the mean stellar age of the cluster, increasing from light to dark. The panels show the 3D velocity dispersion (a), the velocity anisotropy (b), the rotation velocity along a \pm 1 pc slit along an inclined plane of rotation (c) and the tangential velocity (d). Each line is a mean of 100 bootstrapped samples oriented in a random inclination. The dashed line in panel b shows the final anisotropy profile without the removal of the bulk motion, \beta_\mathrm{rot}. The vertical bars in panel a show the projected R_\mathrm{eff}.$

Figure 4: The radial velocity distribution of the most massive cluster in 0.1 Myr steps. The lines are coloured with the mean stellar age of the cluster, increasing from light to dark. The panels show the 3D velocity dispersion (a), the velocity anisotropy (b), the rotation velocity along a $\pm 1$ pc slit along an inclined plane of rotation (c) and the tangential velocity (d). Each line is a mean of 100 bootstrapped samples oriented in a random inclination. The dashed line in panel b shows the final anisotropy profile without the removal of the bulk motion, $\beta_\mathrm{rot}$. The vertical bars in panel a show the projected $R_\mathrm{eff}$.

3.2.2 Velocity distribution↩︎

Next we inspect the velocity dispersion and velocity anisotropy profiles, as well as whether the cluster rotates. We compute the radial and tangential velocity distributions in spherical coordinates $v_\mathrm{rad}$ and $v_\mathrm{tan}=\sqrt{v_\theta^2 + v_\phi^2}$ (where $v_\theta, v_\phi$ are the angular velocity components) as well as using the line of sight (LOS) velocity in a $\pm1$ pc slit perpendicular to the cluster rotation axis, $V_\mathrm{LOS}$. The slit is placed according to the angular momentum vector within the 3D half-mass radius. The velocities are binned radially in 0.1 dex bins and the innermost bin includes all stars within that radius. The corresponding velocity dispersions $\sigma_\mathrm{rad}, \sigma_\theta, \sigma_\phi$, and $\sigma_\mathrm{tan}=\sqrt{\sigma_\theta^2 + \sigma_\phi^2}$ and total $\sigma=\sqrt{\sigma_\mathrm{rad}^2+\sigma_\theta^2 + \sigma_\phi^2}$ are computed using the standard deviation of each velocity component in each bin. The cartesian coordinates are centred at the cluster centre of mass. Here we also mask out stars that have velocities greater than $\sim70$ km s$^{-1}$, in excess of the bulk velocity distribution in the cluster, to remove stars that are either escaping or in close binaries. Standard deviations for the profiles are computed by bootstrapping the stars 100 times, each time in a random orientation to include the effect of inclination on $V_\mathrm{LOS}$. We use the velocity anisotropy parameter $\beta$ to describe the ratio in tangential and radial dispersion [171] as \[\beta=1-\frac{\sigma_\mathrm{tan}^2}{2\sigma_\mathrm{rad}^2}.\]

Fig. 4 shows the evolution of the velocity distribution in the most massive cluster with time. The increase of the velocity dispersion toward the centre steepens with time as the central density of the cluster increases (see Fig. 3). The peak central value is $\sigma\sim 37$ km s$^{-1}$. Once the cluster settles with one dominant core (as opposed to several merging sub-clumps), its peak rotation velocity shifts inwards. By the final snapshot the peak value reaches $V_\mathrm{LOS}\sim7$ km s$^{-1}$ and $v_\mathrm{tan}\sim12$ km s$^{-1}$ at a radius of $\sim 0.3$–$0.4$ pc. For comparison, GCs have similar or lower rotation velocities measured both in 1D and 3D [22], [24], [172], consistent with their dynamically more evolved state. R136 also rotates at a rate of a few km s$^{-1}$ [173], which is slower than our simulated cluster and consistent with its lower mass and previous results on star cluster dynamics in the griffin project [26].

The velocity distribution remains fairly isotropic within the inner 1 pc, while the outer parts develop with a radially biased velocity distribution. Qualitatively similar results were found by [27] in cloud-scale simulations of hierarchical assembly, where the outer regions of the final cluster show radial anisotropy. The stellar population giving rise to the rotation signal and to the anisotropy in our cluster is discussed in more detail in Sections 3.2.4 and 3.2.3.

As described in [24], the anisotropy profile can be computed using the standard deviations or the root mean square velocities. The former gives the intrinsic anisotropy, while the latter includes the bulk velocity component in each radial bin. That is, anisotropy including bulk motion $\beta_\mathrm{rot}$ is computed as $\sigma^2=\frac{1}{N} \sum_{i=0}^{N}{v_i}^2$ instead of $\sigma^2=\frac{1}{N} \sum_{i=1}^{N}{(v_i-\Bar{v})^2}$ where $\Bar{v}$ is the mean velocity per bin. The bulk motion mostly composes of rotation if the cluster is rotating and not expanding or contracting. The final $\beta_\mathrm{rot}$ of the simulated cluster is shown in Fig. 4. The bulk rotation in the inner radii is reflected in a tangentially biased $\beta_\mathrm{rot}$. The radially dominated $\beta_\mathrm{rot}$ at outer radii is reduced compared to $\beta$, however the anisotropy remains mildly positive even when the bulk motion in each component is included.

The rotational support $(V/\sigma)_\mathrm{tot}$, computed as the ratio between peak rotation velocity within $3 R_\mathrm{eff}$ and the central velocity dispersion, is $\sim0.3$–$0.5$ regardless of using the LOS velocity or the tangential rotational velocity. The value of $(V/\sigma)_\mathrm{tot}$ remained at a similar value during the final 1.3 Myr of the cluster assembly. With outward angular momentum transport, $(V/\sigma)_\mathrm{tot}$ is then expected to reduce later with time. The obtained $(V/\sigma)_\mathrm{tot}$ is in agreement with observations of GCs where typical values range between 0.1 and 0.6 in GCs that show statistically significant amounts of ordered rotation and that are dynamically young (relaxation time $>\mathrm{Gyr}$; e.g. [22]). Other definitions for $V/\sigma$, e.g. using the rotation amplitude over the total velocity dispersion observed in GCs in [24], also show similar values and similar correlation with dynamical age as in [22].

Figure 5: The radial velocity distribution of the stars in the most massive cluster divided in to the younger (blue) and older (red) population in the final snapshot. The panels show the stellar number density (a); 3D velocity dispersion (b); radial and tangential velocity dispersion (c); velocity anisotropy excluding (\beta, solid) and including (\beta_\mathrm{rot}, dashed) the bulk velocities (d); the projected LOS rotation velocity in a \pm 1 pc slit and the tangential rotation velocity (solid and dashed lines, respectively; e); and the corresponding V/\sigma of the tangential and LOS velocity components (f). The vertical gray line in the panels indicates the projected R_\mathrm{eff}, and the blue and red vertical lines in panel a shows the values of R_\mathrm{eff} for the young and old populations, respectively. — Figure 5: The radial velocity distribution of the stars in the most massive cluster divided in to the younger (blue) and older (red) population in the final snapshot. The panels show the stellar number density (a); 3D velocity dispersion (b); radial and tangential velocity dispersion (c); velocity anisotropy excluding ($\beta$, solid) and including ($\beta_\mathrm{rot}$, dashed) the bulk velocities (d); the projected LOS rotation velocity in a $\pm 1$ pc slit and the tangential rotation velocity (solid and dashed lines, respectively; e); and the corresponding $V/\sigma$ of the tangential and LOS velocity components (f). The vertical gray line in the panels indicates the projected $R_\mathrm{eff}$, and the blue and red vertical lines in panel a shows the values of $R_\mathrm{eff}$ for the young and old populations, respectively.

Figure 6: The shape parameters of the most massive cluster in the final snapshot (black) divided additionally by age into two groups: stars younger (blue) and older (red) than the mean stellar age of the cluster. The panels show the axis ratios b/a and c/a (solid and dashed, respectively; panel a), the ellipticity (panel b) and the triaxiality (panel c). All radial profiles are computed for stars enclosed within the increasing radial extent. The vertical lines show the projected R_\mathrm{eff}. The dashed horizontal lines in panel c indicate the triaxiality values of 2/3 and 1/3 above and below which the shape is prolate and oblate, respectively. — Figure 6: The shape parameters of the most massive cluster in the final snapshot (black) divided additionally by age into two groups: stars younger (blue) and older (red) than the mean stellar age of the cluster. The panels show the axis ratios $b/a$ and $c/a$ (solid and dashed, respectively; panel a), the ellipticity (panel b) and the triaxiality (panel c). All radial profiles are computed for stars enclosed within the increasing radial extent. The vertical lines show the projected $R_\mathrm{eff}$. The dashed horizontal lines in panel c indicate the triaxiality values of $2/3$ and $1/3$ above and below which the shape is prolate and oblate, respectively.

3.2.3 Kinematics of the age groups↩︎

Finally we investigate whether the hierarchical cluster assembly leaves an imprint in the shape and kinematics of the stars. We divide the stars into two age groups either younger or older than the mean stellar age of the cluster in the final snapshot. The groups of stars have similar total masses within 7%. The older population formed mainly in sub-clumps that accreted hierarchically before the cluster transitioned into the centrally concentrated star formation mode. The younger population corresponds to the latter. To draw a qualitative parallel to the MPs of GCs, we may suppose that the formation of the chemically peculiar P2 occurs after the first VMSs have had time to enrich their surrounding, and thus the youngest stars in the cluster should also be the most chemically enriched. The age groups do not signify generations, but simply the two $\sim3$ Myr halves of the fairly Gaussian star formation history of this cluster that proceeded over a time span of $\sim6$ Myr (see Fig. 7 in [89] for a similar star formation history).

We compute the kinematic properties for these two populations in the same radial bins used for the whole stellar population. Fig. 5 shows the number density and velocity profiles of the two age groups. The young component is slightly more centrally concentrated with $R_\mathrm{eff}\sim 0.56$ pc, compared to $R_\mathrm{eff}\sim 0.96$ pc of the old component and $R_\mathrm{eff}\sim 0.75$ pc of the entire cluster. The young component has a significantly higher rotation velocity compared to the old component. The peak rotation rate of the old population of a few km s$^{-1}$ is consistent with the residual random angular momentum resulting from the hierarchical assembly process, based on initially non-rotating hierarchical N-body models of similar mass ([103], [174]). The rotation velocity seen in Fig. 4 (panels c and d) is therefore mainly dominated by the younger population of stars that has not yet relaxed unlike the old population. The radially biased anisotropy profile seen in Fig. 4 (panel b) at the outer parts is dominated by the younger component as well.

If we do not reduce the bulk velocity and only use the root mean square velocities to compute $\beta_\mathrm{rot}$, the young population shows at intermediate radii ($0.1$ pc $<r<2$ pc) a strong tangential bias (down to $\sim-1$) compared to the fairly isotropic profile at $r\lesssim R_\mathrm{eff}$ shown in Fig. 5. This is simply a reflection of the rotation, that is encoded in bulk tangential motion. The radial bias in the outer radii of the young population is still retained even with bulk rotation included, as already shown in Fig. 4. The old population, on the other hand, remains isotropic when bulk rotation is not removed. The only exception is in the outermost parts where the $\beta_\mathrm{rot}$ of the old population shows tangential bias, i.e. bulk rotation.

$(V/\sigma)_\mathrm{tot}$ for the young and old populations are 0.6 and 0.1, respectively. Fig 5 (panel f) shows the radial $V/\sigma$ profiles of the two components. The intermediate and inner radii show very different $V/\sigma$ shapes for the young and old populations in line with the differences in the rotation velocity profiles. In the outer radii, however, the populations show better agreement in their $V/\sigma$ profiles, reflecting their similar rotation rates at the outer radii. The only difference in the outer parts is the significant radial component in the velocity dispersion of the young population (as seen in $\beta$ and $\beta_\mathrm{rot}$), which is the result of young sub-clusters falling in on radial orbits.

In the absence of strong chemical variations as expected of MPs, another way of differentiating the sub-populations in addition to age would be based on their birth location. As the most massive stars are located in the core of the cluster, then this would be expected as the location where the most chemically enriched stars form. Our previous simulations indicate that the chemical enrichment does indeed occur concentrated in the central regions of massive clusters [89]. We tested performing the above analysis using a division into two populations in which one (here the second population) contains stars born within the half-mass radius of the most massive progenitor sub-cluster. The other population consists of stars that originate from the outer regions of the progenitor or from accreted sub-clumps (here the first population). The results for Fig. 5 remain qualitatively similar. The peak rotation rate of the inner population is reduced by 1–2 km s$^{-1}$ from that shown for the young population in panel $e$ of Fig. 5. The rotation rate of the outer population peaks with $\sim5$ km s$^{-1}$ at $\sim 1$ pc. Some outer radial anisotropy is introduced into the outer population compared to the mostly isotropic old stars in panel $d$ of Fig. 5. The most drastic change is, however, in the density profiles. The inner-formed population remains strongly centrally concentrated and dominates the density profile within $r_\mathrm{eff}$ even though the stars only constitute 1/3rd of the total cluster mass. To regain a less biased prediction of the density distributions, we chose to show only the populations divided by age in Fig. 5. Stellar age can be viewed as a proxy for the progress of chemical enrichment in the cluster, and only selecting by age we allow later-formed stars also in the accreted sub-clusters to be considered a part of the population that would be analogous to P2 in observed GCs.

With this, we conclude that our results confirm the picture drawn out in previous simulation of star formation in pre-existing massive star clusters ([42]–[44]) wherein the second formed population rotates significantly more than the first population. The previous simulations assumed a gradual initial rotation velocity for the first population, while here we self-consistently follow the formation and evolution of both two populations.

3.2.4 Shapes of the age groups↩︎

Finally we compute in 3D the axis ratios $b/a$ and $c/a$, ellipticity $\epsilon=1-c/a$ (i.e. maximal ellipticity given by the longest and shortest principal axes) and triaxiality parameter $T=\frac{a^2-b^2}{a^2-c^2}$ using the inertial tensor [171] of the enclosed mass in increasing radial bins. The results are shown in Fig. 6.

Fig. 6 shows the shape parameters for stars enclosed within increasing radii in the entire cluster and in the two age groups. The core of the cluster has a significant ellipticity (up to $\epsilon\sim0.6$) and a mildly prolate ($T>0.66$) or triaxial ($0.33<T<0.66$) shape reflecting the star formation out of the disky internal gas structure. The cluster becomes increasingly triaxial or oblate ($T<0.33$) when stars at outer radii are included, but note that the overall shape becomes more spherical e.g. with $\epsilon\sim 0.2$ within $R_\mathrm{eff}$. The older population drives down the ellipticity in the cluster as it is fairly triaxial/round, while the young population has a more elliptic shape throughout the cluster.

Our new results verify the mildly elliptic shapes found for simulated young massive clusters in [26], where collisional dynamics was not yet considered. Observed (evolved) GCs show similarly low ellipticities [175], hence we conclude that massive dense clusters may have such spherical shapes already at birth.

3.2.5 Comparison to NGC 104/47 Tuc↩︎

The most massive simulated young cluster has similar rotation characteristics with the Galactic GC NGC 104/47 Tuc that has been widely analysed in kinematic surveys (e.g. [22], [176], [177]). NGC 104 is already $\gtrsim11$ Gyr old [178], [179], however it is still among the dynamically younger GCs [180]. The differences in spatial and kinematic properties of our simulated young and old stellar populations (centrally more concentrated young stars; faster rotating young stars; radially biased outer young stars) are similar to what is found in some studies of NGC 104 that compare the chemically peculiar P2 and chemically normal P1 stars (P2 is more centrally concentrated; P2 rotates faster; P2 has a radial anisotropy in the outer parts; e.g. [24], [41], [181]; but see also e.g. [23] where no statistically significant difference between is found for P2 and P1). While the two age groups defined here are not directly analogous to the chemically identified multiple populations, the striking qualitative resemblance of the spatial and kinematic properties of the young and old population with the observed P2 and P1 in NGC 104 may signify a similar temporal sequence in their formation. In this scenario the P2 of NGC 104 would indeed have formed after P1 in a more centrally concentrated star formation phase. Previous cloud-scale simulations have also suggested that massive stars would form only late in the cluster assembly process [182] supporting the late formation of the chemically peculiar population.

In some of the other dynamically young GCs, such as NGC 3201 and NGC 6101, the outer radii are also radially anisotropic, seemingly driven by the P2 stars [24]. These clusters, however, show no clear differences in the rotation rates of the populations, and they have more centrally concentrated P1 populations as opposed to the more centrally concentrated P2 of NGC 104.

It remains to be tested whether the simulated populations retain their differences over an extended period of time. [183] and [20] for instance showed with N-body simulations that a disk-like young population can retain the input lower velocity dispersion and outer radial anisotropy for a long time. The specific internal dynamical evolution depends on the initial stellar distribution and kinematics, as well as on the evolution of the tidal field (e.g. [184]). The present results would thus offer more self-consistent initial conditions for follow-up N-body simulations.

4 Conclusions and discussion↩︎

We have presented the first galactic scale simulation of massive star cluster formation including collisional stellar dynamics in the vicinity of massive stars. The new model includes prescriptions for stellar collisions and TDEs that allow us to follow the collisional growth of VMSs in the dense cores of star clusters.

The general population of very young ($\lesssim4$ Myr) intermediate mass star clusters (between a few $10^2$ $\mathrm{M}_\odot$ and few $10^4$ $\mathrm{M}_\odot$) matches with the smallest sizes and highest surface densities of observed clusters younger than $10$ Myr in the LEGUS survey. The agreement is expected to further improve once the simulated clusters expand during their first few Myr of dynamical evolution (not included here). These clusters have on average twice the extent of similar clusters in a previous simulation in [89] that started with otherwise similar initial setup but did not implement collisional stellar dynamics.

The most massive star cluster has a gravitationally bound mass of $\sim2\times 10^5$ $\mathrm{M}_\odot$. This is similar to what was found in the previous collisionless simulation. The $\lesssim1$ pc effective radius and surface density of $\sim 7\times 10^4$ $\mathrm{M}_\odot$ pc$^{-2}$ in the most massive cluster approaches the extreme values of surface density and compactness ([$M_\mathrm{cluster}/10^5$ $\mathrm{M}_\odot$ $]/[R_\mathrm{eff}/$ pc$]>1$) of observed young massive star clusters found in recent high-redshift data (e.g. [7]). The density profiles increase toward the cluster centre, without a flat core. The central density (central $>1000$ stars) is more than two orders of magnitude higher and the central surface density an order of magnitude higher compared to the effective values. This result is consistent with the resolved best-fit profiles of [1] where it is predominantly the dynamically young massive clusters (relaxation time $>$Gyrs) that show central-to-effective surface density ratios ($\Sigma_0/\Sigma_\mathrm{eff}$) of more than 10. The implication is that the GC progenitors observed at high-redshifts probably had at least an order of magnitude higher central densities than the resolution limited mean values of $\Sigma_\mathrm{eff}\sim10^5$–$10^6$ $\mathrm{M}_\odot$ pc$^{-2}$ reported thus far.

The most massive star cluster rotates at a rate of $V_\mathrm{LOS}\sim$7 km s$^{-1}$ or $v_\mathrm{tan}\sim$12 km s$^{-1}$, with the radial rotation profile peaking inside the half-mass radius. To decipher which stars dominate the rotation, we further divided the stars of the cluster into a younger and an older half across the $\sim 6$ Myr of star formation in the cluster. We find that the younger population is more centrally concentrated, it has a significantly faster peak rotation velocity ($>12$ km s$^{-1}$ compared to 2–4 km s$^{-1}$ in the old population), and its velocity distribution is radially anisotropic in the outer parts. The first formed, older population is more isotropic and it has a larger velocity dispersion. The cluster has a significant elliptic and mildly prolate shape in the inner parts while the overall shape of the cluster is mildly elliptic ($\epsilon\sim 0.1$–$0.3$). The older population is more triaxial/round in shape throughout the cluster while the young population has a milddy prolate shape in the inner parts and an oblate shape in the outer parts. The first formed population has thus already relaxed throughout the hierarchical assembly process, while the younger population exhibits its initial shape and kinematics, including some radial motions in the outer parts from sub-clump infall.

The radial stellar surface density and surface brightness profiles of the most massive simulated cluster are fairly similar in shape with the most massive local young star cluster R136. The shape and $V/\sigma$ values agree with dynamically young GCs. The kinematic properties of the young and old populations are in qualitatively agreement with those of the chemically peculiar P2 and normal P1 stars, respectively, in the Galactic GC NGC 104. In other dynamically young clusters, such as NGC 3201 and NGC 6101, the differences in the spatio-kinematics of the stellar populations can be less pronounced or even reversed compared to the features found here (e.g. a centrally more concentrated P1). Explaining the variety of population characteristics in GCs will thus require further investigation, including long-timescale N-body simulations.

Our stellar collision product of $\sim 1000$ $\mathrm{M}_\odot$ would result in a black hole of at least $\sim600$ $\mathrm{M}_\odot$ at our adopted metallicity of $0.01$ $\mathrm{Z}_\odot$. The most massive star and the resulting black hole are therefore $\sim0.5\%$ of the total mass of the cluster. This is similar to the lower limits for the most massive stars found in N-body simulations where the collision products typically contain between a fraction of a per cent to $\sim 10\%$ of the cluster mass [10], [103], [174], [185]. Between a few per cent and $10\%$ is also the expected upper limit for extremely massive or supermassive stars forming through collisions or gas accretion based on analytic arguments [161], [186]. We are therefore starting to probe the very high mass end of the stellar mass function, often found theoretically to extend to $10^3$ $\mathrm{M}_\odot$ or $10^4$ $\mathrm{M}_\odot$ (e.g. [70], [159]), now in galactic scale populations of simulated star clusters.

Our future aim is to explore the chemo-dynamical origin of multiple populations in GC progenitor clusters. While strong nitrogen enhancements are easily produced even with small amounts of hot hydrogen burning products, the currently implemented enrichment sources have not been able to produce large numbers of stars with as strong enhancements of, for instance, Al and Mg as observed ([89]). This remains true even if we add the mass loss of up to $400$ $\mathrm{M}_\odot$ from collisions and VMS winds related to the growth and evolution of the most massive star in the presently most massive cluster. Additional mass loss, for instance, through disruption of stars [187] or the growth of the most massive star to the extreme or supermassive regime [161], [188] will be needed to supply the star forming gas rapidly with significant amounts of light elements to explain the abundance anomalies in GCs. The growth could be enhanced through increased collision rates via primordial binaries [156], [157], [189] or with sustained gas accretion and associated mass loss [186], [190]. The stellar evolution methodology can also be expanded in the future with more accurate criteria and non-mass conserving prescriptions for mass transfer and collisions (e.g. adopting the Roche lobe overflow criterion from [191]). The chemical enrichment due to TDEs [192] and binary evolution [193], [194] may be explored with a method that returns the lost mass back into the ISM.

Acknowledgments↩︎

The authors thank Angela Adamo, Ellen Leitinger, Michael Hilker and Holger Baumgardt for clarifying discussions about the observational data. NL and TN gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SUPERMUC-NG at Leibniz Supercomputing Centre (www.lrz.de) under project numbers pn49qi. TN acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311 from the DFG Cluster of Excellence "ORIGINS". The computations were carried out at SuperMUC-NG hosted by the LRZ and the FREYA cluster hosted by The Max Planck Computing and Data Facility (MPCDF) in Garching, Germany.

This research made use of python packages scipy [195], numpy [196], matplotlib [197], pygad [198], h5py [199], and astropy⁴ [200].

Data availability statement↩︎

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

References↩︎

[1]

McLaughlin D. E., van der Marel R. P., 2005, @doi []10.1086/497429, https://ui.adsabs.harvard.edu/abs/2005ApJS..161..304M.

[2]

Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, @doi []10.1146/annurev-astro-081309-130834, https://ui.adsabs.harvard.edu/abs/2010ARA&A..48..431P.

[3]

Misgeld I., Hilker M., 2011, @doi []10.1111/j.1365-2966.2011.18669.x, https://ui.adsabs.harvard.edu/abs/2011MNRAS.414.3699M.

[4]

McCrady N., Graham J. R., 2007, @doi []10.1086/518357, https://ui.adsabs.harvard.edu/abs/2007ApJ...663..844M.

[5]

Levy R. C., et al., 2024, @doi []10.3847/2041-8213/ad7af3, https://ui.adsabs.harvard.edu/abs/2024ApJ...973L..55L.

[6]

Vanzella E., et al., 2023, @doi []10.3847/1538-4357/acb59a, https://ui.adsabs.harvard.edu/abs/2023ApJ...945...53V.

[7]

Adamo A., et al., 2024, @doi []10.1038/s41586-024-07703-7, https://ui.adsabs.harvard.edu/abs/2024Natur.632..513A.

[8]

Gold T., Axford W. I., Ray E. C., 1965, in Robinson I., Schild A., Schucking E. L., eds, Quasi-Stellar Sources and Gravitational Collapse. p. 93.

[9]

Spitzer Lyman J., Saslaw W. C., 1966, @doi []10.1086/148523, https://ui.adsabs.harvard.edu/abs/1966ApJ...143..400S.

[10]

Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillan S. L. W., 2004, @doi []10.1038/nature02448, https://ui.adsabs.harvard.edu/abs/2004Natur.428..724P.

[11]

Krause M. G. H., Charbonnel C., Bastian N., Diehl R., 2016, @doi []10.1051/0004-6361/201526685, https://ui.adsabs.harvard.edu/abs/2016A&A...587A..53K.

[12]

Vink J. S., 2018, @doi []10.1051/0004-6361/201832773, https://ui.adsabs.harvard.edu/abs/2018A&A...615A.119V.

[13]

Crowther P. A., et al., 2016, @doi []10.1093/mnras/stw273, https://ui.adsabs.harvard.edu/abs/2016MNRAS.458..624C.

[14]

Andersen M., Zinnecker H., Moneti A., McCaughrean M. J., Brandl B., Brandner W., Meylan G., Hunter D., 2009, @doi []10.1088/0004-637X/707/2/1347, https://ui.adsabs.harvard.edu/abs/2009ApJ...707.1347A.

[15]

Harayama Y., Eisenhauer F., Martins F., 2008, @doi []10.1086/524650, https://ui.adsabs.harvard.edu/abs/2008ApJ...675.1319H.

[16]

Lim B., Chun M.-Y., Sung H., Park B.-G., Lee J.-J., Sohn S. T., Hur H., Bessell M. S., 2013, @doi []10.1088/0004-6256/145/2/46, https://ui.adsabs.harvard.edu/abs/2013AJ....145...46L.

[17]

Ascenso J., Alves J., Beletsky Y., Lago M. T. V. T., 2007, @doi []10.1051/0004-6361:20066433, https://ui.adsabs.harvard.edu/abs/2007A&A...466..137A.

[18]

Espinoza P., Selman F. J., Melnick J., 2009, @doi []10.1051/0004-6361/20078597, https://ui.adsabs.harvard.edu/abs/2009A&A...501..563E.

[19]

Vesperini E., McMillan S. L. W., D’Antona F., D’Ercole A., 2013, @doi []10.1093/mnras/sts434, https://ui.adsabs.harvard.edu/abs/2013MNRAS.429.1913V.

[20]

Hénault-Brunet V., Gieles M., Agertz O., Read J. I., 2015, @doi []10.1093/mnras/stv675, https://ui.adsabs.harvard.edu/abs/2015MNRAS.450.1164H.

[21]

Tiongco M. A., Vesperini E., Varri A. L., 2017, @doi []10.1093/mnras/stx853, https://ui.adsabs.harvard.edu/abs/2017MNRAS.469..683T.

[22]

Bianchini P., van der Marel R. P., del Pino A., Watkins L. L., Bellini A., Fardal M. A., Libralato M., Sills A., 2018, @doi []10.1093/mnras/sty2365, https://ui.adsabs.harvard.edu/abs/2018MNRAS.481.2125B.

[23]

Martens S., et al., 2023, @doi []10.1051/0004-6361/202244787, https://ui.adsabs.harvard.edu/abs/2023A&A...671A.106M.

[24]

Leitinger E. I., Baumgardt H., Cabrera-Ziri I., Hilker M., Carbajo-Hijarrubia J., Gieles M., Husser T. O., Kamann S., 2025, @doi []10.1051/0004-6361/202452477, https://ui.adsabs.harvard.edu/abs/2025A&A...694A.184L.

[25]

Mapelli M., 2017, @doi []10.1093/mnras/stx304, https://ui.adsabs.harvard.edu/abs/2017MNRAS.467.3255M.

[26]

Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., 2020b, @doi []10.3847/1538-4357/abc001, https://ui.adsabs.harvard.edu/abs/2020ApJ...904...71L.

[27]

Karam J., Fujii M. S., Sills A., 2025, @doi []10.3847/1538-4357/adc719, https://ui.adsabs.harvard.edu/abs/2025ApJ...984...75K.

[28]

Gratton R. G., et al., 2001, @doi []10.1051/0004-6361:20010144, https://ui.adsabs.harvard.edu/abs/2001A&A...369...87G.

[29]

Carretta E., Bragaglia A., Gratton R., Lucatello S., 2009, @doi []10.1051/0004-6361/200912097, https://ui.adsabs.harvard.edu/abs/2009A&A...505..139C.

[30]

Piotto G., et al., 2015, @doi []10.1088/0004-6256/149/3/91, https://ui.adsabs.harvard.edu/abs/2015AJ....149...91P.

[31]

Pancino E., et al., 2017, @doi []10.1051/0004-6361/201730474, https://ui.adsabs.harvard.edu/abs/2017A&A...601A.112P.

[32]

Marino A. F., et al., 2019, @doi []10.1093/mnras/stz1415, https://ui.adsabs.harvard.edu/abs/2019MNRAS.487.3815M.

[33]

Cordoni G., et al., 2025, @doi []10.1093/mnras/staf102, https://ui.adsabs.harvard.edu/abs/2025MNRAS.537.2342C.

[34]

Renzini A., et al., 2015, @doi []10.1093/mnras/stv2268, https://ui.adsabs.harvard.edu/abs/2015MNRAS.454.4197R.

[35]

Charbonnel C., 2016, in EAS Publications Series. pp 177–226 (@eprint arXiv1611.08855), @doi10.1051/eas/1680006.

[36]

Bastian N., Lardo C., 2018, @doi []10.1146/annurev-astro-081817-051839, https://ui.adsabs.harvard.edu/abs/2018ARA&A..56...83B.

[37]

Gratton R., Bragaglia A., Carretta E., D’Orazi V., Lucatello S., Sollima A., 2019, @doi []10.1007/s00159-019-0119-3, https://ui.adsabs.harvard.edu/abs/2019A&ARv..27....8G.

[38]

Milone A. P., Marino A. F., 2022, @doi [Universe]10.3390/universe8070359, https://ui.adsabs.harvard.edu/abs/2022Univ....8..359M.

[39]

Dalessandro E., et al., 2019, @doi []10.3847/2041-8213/ab45f7, https://ui.adsabs.harvard.edu/abs/2019ApJ...884L..24D.

[40]

Leitinger E., Baumgardt H., Cabrera-Ziri I., Hilker M., Pancino E., 2023, @doi []10.1093/mnras/stad093, https://ui.adsabs.harvard.edu/abs/2023MNRAS.520.1456L.

[41]

Dalessandro E., et al., 2024, @doi []10.1051/0004-6361/202451054, https://ui.adsabs.harvard.edu/abs/2024A&A...691A..94D.

[42]

Bekki K., 2010, @doi []10.1088/2041-8205/724/1/L99, https://ui.adsabs.harvard.edu/abs/2010ApJ...724L..99B.

[43]

McKenzie M., Bekki K., 2021, @doi []10.1093/mnras/staa3376, https://ui.adsabs.harvard.edu/abs/2021MNRAS.500.4578M.

[44]

Lacchin E., Calura F., Vesperini E., Mastrobuono-Battisti A., 2022, @doi []10.1093/mnras/stac2328, https://ui.adsabs.harvard.edu/abs/2022MNRAS.517.1171L.

[45]

Richer H. B., Heyl J., Anderson J., Kalirai J. S., Shara M. M., Dotter A., Fahlman G. G., Rich R. M., 2013, @doi []10.1088/2041-8205/771/1/L15, https://ui.adsabs.harvard.edu/abs/2013ApJ...771L..15R.

[46]

Bellini A., et al., 2015, @doi []10.1088/2041-8205/810/1/L13, https://ui.adsabs.harvard.edu/abs/2015ApJ...810L..13B.

[47]

Decressin T., Baumgardt H., Kroupa P., 2008, @doi []10.1051/0004-6361:200810275, https://ui.adsabs.harvard.edu/abs/2008A&A...492..101D.

[48]

Pavlı́k V., Heggie D. C., Varri A. L., Vesperini E., 2024, @doi []10.1051/0004-6361/202450270, https://ui.adsabs.harvard.edu/abs/2024A&A...689A.313P.

[49]

Renaud F., Bournaud F., Duc P.-A., 2015, @doi []10.1093/mnras/stu2208, https://ui.adsabs.harvard.edu/abs/2015MNRAS.446.2038R.

[50]

Ma X., et al., 2020, @doi []10.1093/mnras/staa527, https://ui.adsabs.harvard.edu/abs/2020MNRAS.493.4315M.

[51]

Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., Steinwandel U. P., Moster B. P., 2020a, @doi []10.3847/1538-4357/ab7190, https://ui.adsabs.harvard.edu/abs/2020ApJ...891....2L.

[52]

Hirai Y., Fujii M. S., Saitoh T. R., 2021, @doi []10.1093/pasj/psab038, https://ui.adsabs.harvard.edu/abs/2021PASJ...73.1036H.

[53]

Dobbs C. L., Bending T. J. R., Pettitt A. R., Bate M. R., 2022, @doi []10.1093/mnras/stab3036, https://ui.adsabs.harvard.edu/abs/2022MNRAS.509..954D.

[54]

Grudić M. Y., Hafen Z., Rodriguez C. L., Guszejnov D., Lamberts A., Wetzel A., Boylan-Kolchin M., Faucher-Giguère C.-A., 2023, @doi []10.1093/mnras/stac3573, https://ui.adsabs.harvard.edu/abs/2023MNRAS.519.1366G.

[55]

Sameie O., et al., 2023, @doi []10.1093/mnras/stad1071, https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.1800S.

[56]

Mayer L., van Donkelaar F., Messa M., Capelo P. R., Adamo A., 2025, @doi []10.3847/2041-8213/adadfe, https://ui.adsabs.harvard.edu/abs/2025ApJ...981L..28M.

[57]

Calura F., et al., 2025, @doi []10.1051/0004-6361/202452876, https://ui.adsabs.harvard.edu/abs/2025A&A...698A.207C.

[58]

Andersson E. P., Mac Low M.-M., Agertz O., Renaud F., Li H., 2024, @doi []10.1051/0004-6361/202347792, https://ui.adsabs.harvard.edu/abs/2024A&A...681A..28A.

[59]

Gutcke T. A., 2024, @doi []10.3847/1538-4357/ad5c62, https://ui.adsabs.harvard.edu/abs/2024ApJ...971..103G.

[60]

Reina-Campos M., Gnedin O. Y., Sills A., Li H., 2025, @doi []10.3847/1538-4357/ad909f, https://ui.adsabs.harvard.edu/abs/2025ApJ...978...15R.

[61]

Williams C. E., et al., 2025, @doi []10.3847/1538-4357/adf19d, https://ui.adsabs.harvard.edu/abs/2025ApJ...990..135W.

[62]

Matsui H., Bekki K., McKenzie M., Saitoh T. R., 2025, @doi []10.1093/pasj/psaf004, https://ui.adsabs.harvard.edu/abs/2025PASJ..tmp...12M.

[63]

Johnson L. C., et al., 2017, @doi []10.3847/1538-4357/aa6a1f, https://ui.adsabs.harvard.edu/abs/2017ApJ...839...78J.

[64]

Adamo A., et al., 2020, @doi []10.1093/mnras/staa2380, https://ui.adsabs.harvard.edu/abs/2020MNRAS.499.3267A.

[65]

Aarseth S. J., 2003, Gravitational N-Body Simulations. Cambridge University Press, Cambridge, UK.

[66]

Heggie D., Hut P., 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics. Cambridge University Press, Cambridge, UK.

[67]

Fujii M. S., Saitoh T. R., Hirai Y., Wang L., 2021, @doi []10.1093/pasj/psab061, https://ui.adsabs.harvard.edu/abs/2021PASJ...73.1074F.

[68]

Grudić M. Y., Guszejnov D., Offner S. S. R., Rosen A. L., Raju A. N., Faucher-Giguère C.-A., Hopkins P. F., 2022, @doi []10.1093/mnras/stac526, https://ui.adsabs.harvard.edu/abs/2022MNRAS.512..216G.

[69]

Farias J. P., Offner S. S. R., Grudić M. Y., Guszejnov D., Rosen A. L., 2024, @doi []10.1093/mnras/stad3609, https://ui.adsabs.harvard.edu/abs/2024MNRAS.527.6732F.

[70]

Fujii M. S., Wang L., Tanikawa A., Hirai Y., Saitoh T. R., 2024, @doi [Science]10.1126/science.adi4211, https://ui.adsabs.harvard.edu/abs/2024Sci...384.1488F.

[71]

Polak B., et al., 2024, @doi []10.1051/0004-6361/202348840, https://ui.adsabs.harvard.edu/abs/2024A&A...690A..94P.

[72]

Cournoyer-Cloutier C., et al., 2024, @doi []10.3847/1538-4357/ad90b3, https://ui.adsabs.harvard.edu/abs/2024ApJ...977..203C.

[73]

Jo Y., Kim S., Kim J.-h., Bryan G. L., 2024, @doi []10.3847/1538-4357/ad6b16, https://ui.adsabs.harvard.edu/abs/2024ApJ...974..193J.

[74]

Lahén N., Rantala A., Naab T., Partmann C., Johansson P. H., Hislop J. M., 2025, @doi []10.1093/mnras/staf350, https://ui.adsabs.harvard.edu/abs/2025MNRAS.538.2129L.

[75]

Rantala A., Pihajoki P., Johansson P. H., Naab T., Lahén N., Sawala T., 2017, @doi []10.3847/1538-4357/aa6d65, https://ui.adsabs.harvard.edu/abs/2017ApJ...840...53R.

[76]

Mannerkoski M., Rawlings A., Johansson P. H., Naab T., Rantala A., Springel V., Irodotou D., Liao S., 2023, @doi []10.1093/mnras/stad2139, https://ui.adsabs.harvard.edu/abs/2023MNRAS.524.4062M.

[77]

Springel V., 2005, @doi []10.1111/j.1365-2966.2005.09655.x, https://ui.adsabs.harvard.edu/abs/2005MNRAS.364.1105S.

[78]

Hu C.-Y., Naab T., Walch S., Moster B. P., Oser L., 2014, @doi []10.1093/mnras/stu1187, https://ui.adsabs.harvard.edu/abs/2014MNRAS.443.1173H.

[79]

Hu C.-Y., Naab T., Walch S., Glover S. C. O., Clark P. C., 2016, @doi []10.1093/mnras/stw544, https://ui.adsabs.harvard.edu/abs/2016MNRAS.458.3528H.

[80]

Hu C.-Y., Naab T., Glover S. C. O., Walch S., Clark P. C., 2017, @doi []10.1093/mnras/stx1773, https://ui.adsabs.harvard.edu/abs/2017MNRAS.471.2151H.

[81]

Lahén N., et al., 2023, @doi []10.1093/mnras/stad1147, https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.3092L.

[82]

Partmann C., Naab T., Lahén N., Rantala A., Hirschmann M., Hislop J. M., Petersson J., Johansson P. H., 2025, @doi []10.1093/mnras/staf002, https://ui.adsabs.harvard.edu/abs/2025MNRAS.537..956P.

[83]

Hu C.-Y., 2019, @doi []10.1093/mnras/sty3252, https://ui.adsabs.harvard.edu/abs/2019MNRAS.483.3363H.

[84]

Nelson R. P., Langer W. D., 1997, @doi []10.1086/304167, https://ui.adsabs.harvard.edu/abs/1997ApJ...482..796N.

[85]

Glover S. C. O., Mac Low M.-M., 2007, @doi []10.1086/512238, https://ui.adsabs.harvard.edu/abs/2007ApJS..169..239G.

[86]

Glover S. C. O., Clark P. C., 2012, @doi []10.1111/j.1365-2966.2011.20260.x, https://ui.adsabs.harvard.edu/abs/2012MNRAS.421..116G.

[87]

Wiersma R. P. C., Schaye J., Smith B. D., 2009, @doi []10.1111/j.1365-2966.2008.14191.x, https://ui.adsabs.harvard.edu/abs/2009MNRAS.393...99W.

[88]

Kroupa P., 2001, @doi []10.1046/j.1365-8711.2001.04022.x, https://ui.adsabs.harvard.edu/abs/2001MNRAS.322..231K.

[89]

Lahén N., Naab T., Szécsi D., 2024, @doi []10.1093/mnras/stae904, https://ui.adsabs.harvard.edu/abs/2024MNRAS.530..645L.

[90]

Chieffi A., Limongi M., 2004, @doi []10.1086/392523, https://ui.adsabs.harvard.edu/abs/2004ApJ...608..405C.

[91]

Heger A., Woosley S. E., 2002, @doi []10.1086/338487, https://ui.adsabs.harvard.edu/abs/2002ApJ...567..532H.

[92]

Karakas A. I., 2010, @doi []10.1111/j.1365-2966.2009.16198.x, https://ui.adsabs.harvard.edu/abs/2010MNRAS.403.1413K.

[93]

Georgy C., et al., 2013, @doi []10.1051/0004-6361/201322178, https://ui.adsabs.harvard.edu/abs/2013A&A...558A.103G.

[94]

Szécsi D., Agrawal P., Wünsch R., Langer N., 2022, @doi []10.1051/0004-6361/202141536, https://ui.adsabs.harvard.edu/abs/2022A&A...658A.125S.

[95]

Westera P., Lejeune T., Buser R., Cuisinier F., Bruzual G., 2002, @doi []10.1051/0004-6361:20011493, https://ui.adsabs.harvard.edu/abs/2002A&A...381..524W.

[96]

Groh J. H., et al., 2019, @doi []10.1051/0004-6361/201833720, https://ui.adsabs.harvard.edu/abs/2019A&A...627A..24G.

[97]

Clark P. C., Glover S. C. O., Klessen R. S., 2012, @doi []10.1111/j.1365-2966.2011.20087.x, https://ui.adsabs.harvard.edu/abs/2012MNRAS.420..745C.

[98]

Górski K. M., Hivon E., 2011, HEALPix: Hierarchical Equal Area isoLatitude Pixelization of a sphere, Astrophysics Source Code Library (@eprint ascl1107.018).

[99]

Rantala A., Pihajoki P., Mannerkoski M., Johansson P. H., Naab T., 2020, @doi []10.1093/mnras/staa084, https://ui.adsabs.harvard.edu/abs/2020MNRAS.492.4131R.

[100]

Rantala A., Naab T., Rizzuto F. P., Mannerkoski M., Partmann C., Lautenschütz K., 2023, @doi []10.1093/mnras/stad1360, https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.5180R.

[101]

Portegies Zwart S. F., McMillan S. L. W., 2002, @doi []10.1086/341798, https://ui.adsabs.harvard.edu/abs/2002ApJ...576..899P.

[102]

Rizzuto F. P., et al., 2021, @doi []10.1093/mnras/staa3634, https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.5257R.

[103]

Rantala A., Naab T., Lahén N., 2024, @doi []10.1093/mnras/stae1413, https://ui.adsabs.harvard.edu/abs/2024MNRAS.531.3770R.

[104]

Lombardi Jr. J. C., Warren J. S., Rasio F. A., Sills A., Warren A. R., 2002, @doi []10.1086/339060, https://ui.adsabs.harvard.edu/abs/2002ApJ...568..939L.

[105]

Glebbeek E., Gaburov E., Portegies Zwart S., Pols O. R., 2013, @doi []10.1093/mnras/stt1268, https://ui.adsabs.harvard.edu/abs/2013MNRAS.434.3497G.

[106]

Glebbeek E., Gaburov E., de Mink S. E., Pols O. R., Portegies Zwart S. F., 2009, @doi []10.1051/0004-6361/200810425, https://ui.adsabs.harvard.edu/abs/2009A&A...497..255G.

[107]

de Mink S. E., Langer N., Izzard R. G., Sana H., de Koter A., 2013, @doi []10.1088/0004-637X/764/2/166, https://ui.adsabs.harvard.edu/abs/2013ApJ...764..166D.

[108]

Menon A., et al., 2024, @doi []10.3847/2041-8213/ad2074, https://ui.adsabs.harvard.edu/abs/2024ApJ...963L..42M.

[109]

Hurley J. R., Pols O. R., Tout C. A., 2000, @doi []10.1046/j.1365-8711.2000.03426.x, https://ui.adsabs.harvard.edu/abs/2000MNRAS.315..543H.

[110]

Agrawal P., Hurley J., Stevenson S., Szécsi D., Flynn C., 2020, @doi []10.1093/mnras/staa2264, https://ui.adsabs.harvard.edu/abs/2020MNRAS.497.4549A.

[111]

Iorio G., et al., 2023, @doi []10.1093/mnras/stad1630, https://ui.adsabs.harvard.edu/abs/2023MNRAS.524..426I.

[112]

Tout C. A., Aarseth S. J., Pols O. R., Eggleton P. P., 1997, @doi []10.1093/mnras/291.4.732, https://ui.adsabs.harvard.edu/abs/1997MNRAS.291..732T.

[113]

Portegies Zwart S. F., Makino J., McMillan S. L. W., Hut P., 1999, @doi []10.48550/arXiv.astro-ph/9812006, https://ui.adsabs.harvard.edu/abs/1999A&A...348..117P.

[114]

Costa G., et al., 2025, @doi []10.1051/0004-6361/202452573, https://ui.adsabs.harvard.edu/abs/2025A&A...694A.193C.

[115]

Kochanek C. S., 1992, @doi []10.1086/170966, https://ui.adsabs.harvard.edu/abs/1992ApJ...385..604K.

[116]

Guillochon J., Ramirez-Ruiz E., 2013, @doi []10.1088/0004-637X/767/1/25, https://ui.adsabs.harvard.edu/abs/2013ApJ...767...25G.

[117]

Metzger B. D., Stone N. C., 2016, @doi []10.1093/mnras/stw1394, https://ui.adsabs.harvard.edu/abs/2016MNRAS.461..948M.

[118]

Baumgarte T. W., Shapiro S. L., 2025, @doi []10.1103/PhysRevD.111.063039, https://ui.adsabs.harvard.edu/abs/2025PhRvD.111f3039B.

[119]

Giersz M., Leigh N., Hypki A., Lützgendorf N., Askar A., 2015, @doi []10.1093/mnras/stv2162, https://ui.adsabs.harvard.edu/abs/2015MNRAS.454.3150G.

[120]

Rizzuto F. P., Naab T., Rantala A., Johansson P. H., Ostriker J. P., Stone N. C., Liao S., Irodotou D., 2023, @doi []10.1093/mnras/stad734, https://ui.adsabs.harvard.edu/abs/2023MNRAS.521.2930R.

[121]

Rees M. J., 1988, @doi []10.1038/333523a0, https://ui.adsabs.harvard.edu/abs/1988Natur.333..523R.

[122]

Brown G., Gnedin O. Y., 2021, @doi []10.1093/mnras/stab2907, https://ui.adsabs.harvard.edu/abs/2021MNRAS.508.5935B.

[123]

Vanzella E., et al., 2022a, @doi []10.1051/0004-6361/202141590, https://ui.adsabs.harvard.edu/abs/2022A&A...659A...2V.

[124]

Vanzella E., et al., 2022b, @doi []10.3847/2041-8213/ac8c2d, https://ui.adsabs.harvard.edu/abs/2022ApJ...940L..53V.

[125]

Mowla L., et al., 2024, @doi []10.1038/s41586-024-08293-0, https://ui.adsabs.harvard.edu/abs/2024Natur.636..332M.

[126]

Messa M., et al., 2025, @doi []10.1051/0004-6361/202451695, https://ui.adsabs.harvard.edu/abs/2025A&A...694A..59M.

[127]

Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, @doi []10.1046/j.1365-8711.2001.04912.x, https://ui.adsabs.harvard.edu/abs/2001MNRAS.328..726S.

[128]

Dolag K., Borgani S., Murante G., Springel V., 2009, @doi []10.1111/j.1365-2966.2009.15034.x, https://ui.adsabs.harvard.edu/abs/2009MNRAS.399..497D.

[129]

Camps P., Baes M., 2020, @doi [Astronomy and Computing]10.1016/j.ascom.2020.100381, https://ui.adsabs.harvard.edu/abs/2020A&C....3100381C.

[130]

Lahén N., Naab T., Kauffmann G., 2022, @doi []10.1093/mnras/stac1594, https://ui.adsabs.harvard.edu/abs/2022MNRAS.514.4560L.

[131]

Castelli F., Kurucz R. L., 2003, in Piskunov N., Weiss W. W., Gray D. F., eds, IAU Symposium Vol. 210, Modelling of Stellar Atmospheres. p. A20 (@eprint arXivastro-ph/0405087), @doi10.48550/arXiv.astro-ph/0405087.

[132]

Elmegreen B. G., Falgarone E., 1996, @doi []10.1086/178009, http://adsabs.harvard.edu/abs/1996ApJ...471..816E.

[133]

Elmegreen B. G., Efremov Y. N., 1997, @doi []10.1086/303966, https://ui.adsabs.harvard.edu/abs/1997ApJ...480..235E.

[134]

Zhang Q., Fall S. M., 1999, @doi []10.1086/312412, http://adsabs.harvard.edu/abs/1999ApJ...527L..81Z.

[135]

Hunter D. A., Elmegreen B. G., Dupuy T. J., Mortonson M., 2003, @doi []10.1086/378056, https://ui.adsabs.harvard.edu/abs/2003AJ....126.1836H.

[136]

Larsen S. S., 2009, @doi []10.1051/0004-6361:200811212, https://ui.adsabs.harvard.edu/abs/2009A&A...494..539L.

[137]

Randriamanakoto Z., Väisänen P., Ryder S., Kankare E., Kotilainen J., Mattila S., 2013, @doi []10.1093/mnras/stt185, https://ui.adsabs.harvard.edu/abs/2013MNRAS.431..554R.

[138]

Tsang B. T. H., Milosavljević M., 2018, @doi []10.1093/mnras/sty1217, https://ui.adsabs.harvard.edu/abs/2018MNRAS.478.4142T.

[139]

Marks M., Kroupa P., 2012, @doi []10.1051/0004-6361/201118231, https://ui.adsabs.harvard.edu/abs/2012A&A...543A...8M.

[140]

Banerjee S., Kroupa P., 2017, @doi []10.1051/0004-6361/201526928, https://ui.adsabs.harvard.edu/abs/2017A&A...597A..28B.

[141]

Adamo A., et al., 2017, @doi []10.3847/1538-4357/aa7132, https://ui.adsabs.harvard.edu/abs/2017ApJ...841..131A.

[142]

Krumholz M. R., McKee C. F., Bland-Hawthorn J., 2019, @doi []10.1146/annurev-astro-091918-104430, https://ui.adsabs.harvard.edu/abs/2019ARA&A..57..227K.

[143]

Ryon J. E., et al., 2017, @doi []10.3847/1538-4357/aa719e, https://ui.adsabs.harvard.edu/abs/2017ApJ...841...92R.

[144]

Calzetti D., et al., 2023, @doi []10.3847/1538-4357/acbeac, https://ui.adsabs.harvard.edu/abs/2023ApJ...946....1C.

[145]

Whitmore B. C., et al., 2023, @doi []10.3847/2041-8213/acae94, https://ui.adsabs.harvard.edu/abs/2023ApJ...944L..14W.

[146]

Lynden-Bell D., 1967, @doi []10.1093/mnras/136.1.101, https://ui.adsabs.harvard.edu/abs/1967MNRAS.136..101L.

[147]

Silich S., Tenorio-Tagle G., 2017, @doi []10.1093/mnras/stw2879, https://ui.adsabs.harvard.edu/abs/2017MNRAS.465.1375S.

[148]

Grudić M. Y., Hopkins P. F., Quataert E., Murray N., 2019, @doi []10.1093/mnras/sty3386, https://ui.adsabs.harvard.edu/abs/2019MNRAS.483.5548G.

[149]

Charbonnel C., Schaerer D., Prantzos N., Ramı́rez-Galeano L., Fragos T., Kuruvanthodi A., Marques-Chaves R., Gieles M., 2023, @doi []10.1051/0004-6361/202346410, https://ui.adsabs.harvard.edu/abs/2023A&A...673L...7C.

[150]

Isobe Y., et al., 2023, @doi []10.3847/1538-4357/ad09be, https://ui.adsabs.harvard.edu/abs/2023ApJ...959..100I.

[151]

Topping M. W., et al., 2024, @doi []10.1093/mnras/stae682, https://ui.adsabs.harvard.edu/abs/2024MNRAS.529.3301T.

[152]

Marques-Chaves R., et al., 2024, @doi []10.1051/0004-6361/202347411, https://ui.adsabs.harvard.edu/abs/2024A&A...681A..30M.

[153]

Baumgardt H., Klessen R. S., 2011, @doi []10.1111/j.1365-2966.2011.18258.x, https://ui.adsabs.harvard.edu/abs/2011MNRAS.413.1810B.

[154]

Fujii M. S., Portegies Zwart S., 2013, @doi []10.1093/mnras/sts673, https://ui.adsabs.harvard.edu/abs/2013MNRAS.430.1018F.

[155]

Brands S. A., et al., 2022, @doi []10.1051/0004-6361/202142742, https://ui.adsabs.harvard.edu/abs/2022A&A...663A..36B.

[156]

Fregeau J. M., Cheung P., Portegies Zwart S. F., Rasio F. A., 2004, @doi []10.1111/j.1365-2966.2004.07914.x, https://ui.adsabs.harvard.edu/abs/2004MNRAS.352....1F.

[157]

Rantala A., Naab T., 2025, @doi []10.1093/mnrasl/slaf064, https://ui.adsabs.harvard.edu/abs/2025MNRAS.542L..78R.

[158]

Reinoso B., Klessen R. S., Schleicher D., Glover S. C. O., Solar P., 2023, @doi []10.1093/mnras/stad790, https://ui.adsabs.harvard.edu/abs/2023MNRAS.521.3553R.

[159]

Chon S., Omukai K., 2025, @doi []10.1093/mnras/staf598, https://ui.adsabs.harvard.edu/abs/2025MNRAS.539.2561C.

[160]

Sabhahit G. N., Vink J. S., Sander A. A. C., Higgins E. R., 2023, @doi []10.1093/mnras/stad1888, https://ui.adsabs.harvard.edu/abs/2023MNRAS.524.1529S.

[161]

Gieles M., Padoan P., Charbonnel C., Vink J. S., Ramı́rez-Galeano L., 2025, @doi []10.1093/mnras/staf1314, https://ui.adsabs.harvard.edu/abs/2025MNRAS.tmp.1257G.

[162]

Amaro-Seoane P., Konstantinidis S., Dewi Freitag M., Miller M. C., Rasio F. A., 2014, @doi []10.1088/0004-637X/782/2/97, https://ui.adsabs.harvard.edu/abs/2014ApJ...782...97A.

[163]

Arca Sedda M., Amaro Seoane P., Chen X., 2021, @doi []10.1051/0004-6361/202037785, https://ui.adsabs.harvard.edu/abs/2021A&A...652A..54A.

[164]

Whitmore B. C., et al., 2020, @doi []10.3847/1538-4357/ab59e5, https://ui.adsabs.harvard.edu/abs/2020ApJ...889..154W.

[165]

Sills A., Rieder S., Scora J., McCloskey J., Jaffa S., 2018, @doi []10.1093/mnras/sty681, https://ui.adsabs.harvard.edu/abs/2018MNRAS.477.1903S.

[166]

Grudić M. Y., Guszejnov D., Hopkins P. F., Lamberts A., Boylan-Kolchin M., Murray N., Schmitz D., 2018, @doi []10.1093/mnras/sty2303, https://ui.adsabs.harvard.edu/abs/2018MNRAS.481..688G.

[167]

Antonini F., Capuzzo-Dolcetta R., Mastrobuono-Battisti A., Merritt D., 2012, @doi []10.1088/0004-637X/750/2/111, https://ui.adsabs.harvard.edu/abs/2012ApJ...750..111A.

[168]

Mackey A. D., Gilmore G. F., 2003, @doi []10.1046/j.1365-8711.2003.06021.x, https://ui.adsabs.harvard.edu/abs/2003MNRAS.338...85M.

[169]

Elson R. A. W., Fall S. M., Freeman K. C., 1987, @doi []10.1086/165807, https://ui.adsabs.harvard.edu/abs/1987ApJ...323...54E.

[170]

Selman F. J., Melnick J., 2013, @doi []10.1051/0004-6361/201220396, https://ui.adsabs.harvard.edu/abs/2013A&A...552A..94S.

[171]

Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press, Princeton, N.J.

[172]

Kamann S., et al., 2018, @doi []10.1093/mnras/stx2719, https://ui.adsabs.harvard.edu/abs/2018MNRAS.473.5591K.

[173]

Hénault-Brunet V., et al., 2012, @doi []10.1051/0004-6361/201219472, https://ui.adsabs.harvard.edu/abs/2012A&A...545L...1H.

[174]

Rantala A., Lahén N., Naab T., Escobar G. J., Iorio G., 2025, @doi [arXiv e-prints]10.48550/arXiv.2506.04330, https://ui.adsabs.harvard.edu/abs/2025arXiv250604330R.

[175]

Harris W. E., 1996, @doi []10.1086/118116, https://ui.adsabs.harvard.edu/abs/1996AJ....112.1487H.

[176]

Lane R. R., et al., 2010, @doi []10.1111/j.1365-2966.2010.16874.x, https://ui.adsabs.harvard.edu/abs/2010MNRAS.406.2732L.

[177]

Bellini A., Bianchini P., Varri A. L., Anderson J., Piotto G., van der Marel R. P., Vesperini E., Watkins L. L., 2017, @doi []10.3847/1538-4357/aa7c5f, https://ui.adsabs.harvard.edu/abs/2017ApJ...844..167B.

[178]

Gratton R. G., Bragaglia A., Carretta E., Clementini G., Desidera S., Grundahl F., Lucatello S., 2003, @doi []10.1051/0004-6361:20031003, https://ui.adsabs.harvard.edu/abs/2003A&A...408..529G.

[179]

Dotter A., et al., 2010, @doi []10.1088/0004-637X/708/1/698, https://ui.adsabs.harvard.edu/abs/2010ApJ...708..698D.

[180]

Baumgardt H., Hilker M., 2018, @doi []10.1093/mnras/sty1057, https://ui.adsabs.harvard.edu/abs/2018MNRAS.478.1520B.

[181]

Milone A. P., et al., 2012, @doi []10.1088/0004-637X/744/1/5810.1086/141918, https://ui.adsabs.harvard.edu/abs/2012ApJ...744...58M.

[182]

Vázquez-Semadeni E., Palau A., Ballesteros-Paredes J., Gómez G. C., Zamora-Avilés M., 2019, @doi []10.1093/mnras/stz2736, https://ui.adsabs.harvard.edu/abs/2019MNRAS.490.3061V.

[183]

Mastrobuono-Battisti A., Perets H. B., 2013, @doi []10.1088/0004-637X/779/1/85, https://ui.adsabs.harvard.edu/abs/2013ApJ...779...85M.

[184]

Berczik P., Panamarev T., Ishchenko M., Kocsis B., 2025, @doi []10.1051/0004-6361/202452848, https://ui.adsabs.harvard.edu/abs/2025A&A...694A.163B.

[185]

Arca Sedda M., Kamlah A. W. H., Spurzem R., Rizzuto F. P., Naab T., Giersz M., Berczik P., 2023, @doi []10.1093/mnras/stad2292, https://ui.adsabs.harvard.edu/abs/2023MNRAS.526..429A.

[186]

Gieles M., et al., 2018, @doi []10.1093/mnras/sty1059, https://ui.adsabs.harvard.edu/abs/2018MNRAS.478.2461G.

[187]

Elmegreen B. G., 2017, @doi []10.3847/1538-4357/836/1/80, https://ui.adsabs.harvard.edu/abs/2017ApJ...836...80E.

[188]

Denissenkov P. A., Hartwick F. D. A., 2014, @doi []10.1093/mnrasl/slt133, https://ui.adsabs.harvard.edu/abs/2014MNRAS.437L..21D.

[189]

González Prieto E., Weatherford N. C., Fragione G., Kremer K., Rasio F. A., 2024, @doi []10.3847/1538-4357/ad43d6, https://ui.adsabs.harvard.edu/abs/2024ApJ...969...29G.

[190]

Ramı́rez-Galeano L., Charbonnel C., Fragos T., Tazakkati Z., Roman-Garza J., Gieles M., 2025, @doi []10.1051/0004-6361/202453462, https://ui.adsabs.harvard.edu/abs/2025A&A...699A.223R.

[191]

Eggleton P. P., 1983, @doi []10.1086/160960, https://ui.adsabs.harvard.edu/abs/1983ApJ...268..368E.

[192]

Kochanek C. S., 2016, @doi []10.1093/mnras/stw267, https://ui.adsabs.harvard.edu/abs/2016MNRAS.458..127K.

[193]

de Mink S. E., Pols O. R., Langer N., Izzard R. G., 2009, @doi []10.1051/0004-6361/200913205, https://ui.adsabs.harvard.edu/abs/2009A&A...507L...1D.

[194]

Nguyen M., Sills A., 2024, @doi []10.3847/1538-4357/ad4a61, https://ui.adsabs.harvard.edu/abs/2024ApJ...969...18N.

[195]

Virtanen P., et al., 2020, @doi [Nature Methods]10.1038/s41592-019-0686-2, https://rdcu.be/b08Wh.

[196]

Harris C. R., et al., 2020, @doi [Nature]10.1038/s41586-020-2649-2, 585, 357–362.

[197]

Hunter J. D., 2007, @doi [Computing in Science & Engineering]10.1109/MCSE.2007.55, 9, 90.

[198]

Röttgers B., Naab T., Cernetic M., Davé R., Kauffmann G., Borthakur S., Foidl H., 2020, @doi []10.1093/mnras/staa1490, https://ui.adsabs.harvard.edu/abs/2020MNRAS.496..152R.

[199]

Collette A., 2013, Python and HDF5. O’Reilly Media, Inc., Sebostopol, USA.

[200]

Astropy Collaboration et al., 2022, @doi []10.3847/1538-4357/ac7c74, https://ui.adsabs.harvard.edu/abs/2022ApJ...935..167A.

E-mail: nlahen@mpa-garching.mpg.de↩︎
https://wwwmpa.mpa-garching.mpg.de/~naab/griffin-project ↩︎
For stars of $<1.7$ $\mathrm{M}_\odot$ that are not tabulated in the $Z=0.0004$ $\sim 0.02$ $\mathrm{Z}_\odot$ Geneva tables we scale the fluxes from the $Z=0.002$ $\sim 0.1$ $\mathrm{Z}_\odot$ tables of [93] between $0.8$–$1.7$ $\mathrm{M}_\odot$ by a factor of 2. For stars below $0.8$ $\mathrm{M}_\odot$ we extrapolate as $L\propto M^{3.5}$.↩︎
http://www.astropy.org↩︎

Mergers all the way down: stellar collisions and kinematics of a dense hierarchically forming massive star cluster in a dwarf starburst