The formation, evolution, and disruption of star clusters with improved gravitational dynamics in simulated dwarf galaxies

ABSTRACT

So far, even the highest resolution galaxy formation simulations with gravitational softening have failed to reproduce realistic life cycles of star clusters. We present the first star-by-star galaxy models of star cluster formation to account for hydrodynamics, star formation, stellar evolution, and collisional gravitational interactions between stars and compact remnants using the updated sphgal + ketju code, part of the griffin project. Gravitational dynamics in the vicinity of |$>3$| M|$_\odot$| stars and their remnants are solved with a regularized integrator (ketju) without gravitational softening. Comparisons of idealized star cluster evolution with sphgal + ketju and direct N-body show broad agreement and the failure of simulations that use gravitational softening. In the hydrodynamical simulations of idealized dwarf galaxies run with sphgal + ketju, clusters up to |$\sim 900$| M|$_\odot$| form compact (effective radii 0.1–1 pc) and their sizes increase by up to a factor of 10 in agreement with previous N-body simulations and the observed sizes of exposed star clusters. The sizes increase rapidly once the clusters become exposed due to photoionizing radiation. On average 63 per cent of the gravitationally bound clusters disrupt during the first 100 Myr of evolution in the galactic tidal field. The addition of collisional dynamics reduces the fraction of supernovae in bound clusters by a factor of |$\sim 1.7$|⁠; however, the global star formation and outflow histories change by less than 30 per cent. We demonstrate that the accurate treatment of gravitational encounters with massive stars enables more realistic star cluster life cycles from the earliest stages of cluster formation until disruption in simulated low-mass galaxies.

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

1 INTRODUCTION

Star clusters are collisional systems whose internal evolution, after gas-removal, is governed by stellar evolution and gravitational interactions between individual stars, their compact remnants and the galactic tidal field. The picture is especially complex at the early stages of cluster evolution due to ongoing star formation and the interactions between young (massive) stars and their gaseous surroundings. Observationally, the detailed structure of star clusters within the Local Group can be resolved down to individual stars (Kharchenko et al. 2013; Crowther et al. 2016). Multiwavelength studies of integrated light in young star clusters have been used to characterize the ages, masses, metallicities, and sizes of young clusters both in the local Universe (Whitmore et al. 1999; Hunter et al. 2003; Adamo et al. 2017, 2020b; Cook et al. 2023) and at increasingly high redshifts (Adamo et al. 2024; Mowla et al. 2024) thanks to the superb resolving capabilities of e.g. the Hubble Space Telescope (HST), Very Large Telescope (VLT), and the James Webb Space Telescope (JWST).

Young star clusters have mass functions |${\rm d}N/{\rm d}M\propto M^\alpha$| with a power-law index close to |$\alpha =-2$| across different star-forming environments (Elmegreen & Efremov 1996; Elmegreen & Falgarone 1996; Zhang & Fall 1999; Hunter et al. 2003; Fall & Chandar 2012; Mok, Chandar & Fall 2019). A truncation in the mass function at high cluster masses has been suggested (Adamo et al. 2015; Messa et al. 2018; Adamo et al. 2020b), but quantifying this cut-off mass is difficult due to low-number statistics. The mass–size relation of star clusters is observed to be shallow (half-mass or half-light radii |$\propto M^\beta$| with |$\beta \lesssim 1/3$|⁠, see the data collected in Krumholz, McKee & Bland-Hawthorn 2019 and Brown & Gnedin 2021). These measurements exclude the youngest, embedded phase of early cluster evolution that is only revealed at rest-frame infrared wavelengths. The unprecedented sensitivity of JWST enables the detection of embedded star formation in the form of individual proto-stars down to subsolar masses within the nearby galactic environments (De Marchi, Panagia & Beccari 2017; Jones et al. 2023), thus clear improvements to the measurements of embedded cluster sizes can be expected in near future.

Both the amount of star formation occurring in bound clusters and the amount of disruption that clusters undergo are still uncertain. The cluster formation efficiency (CFE or |$\Gamma$|⁠, Bastian 2008) defines the fraction of stellar mass forming in bound clusters. It has been studied within galaxies (Silva-Villa, Adamo & Bastian 2013; Johnson et al. 2016; Messa et al. 2018) and across galaxy samples (Goddard, Bastian & Kennicutt 2010; Adamo et al. 2020b; Chandar et al. 2023; Cook et al. 2023) but its correlation with the star-forming environment has not been firmly established (see the discussion e.g. in Chandar et al. 2017 and Adamo et al. 2020a). Analytic models (Kruijssen 2012) and numerical star formation simulations (Pfeffer et al. 2019; Li et al. 2022; Grudić et al. 2023) often find that |$\Gamma$| increases with increasing star formation rate (SFR) surface density in massive or starbursting galaxies, while the scatter in |$\Gamma$| increases towards lower galaxy mass and no environmental correlation is found (Hislop et al. 2022; Lahén et al. 2023; Andersson et al. 2024). Some observations indicate that the cluster disruption rate depends on the cluster mass (Lamers et al. 2005; Bastian et al. 2012), while others find a universal disruption rate (Chandar, Fall & Whitmore 2010; Fall & Chandar 2012). The disruption rate has been suggested to depend on the galactic environment (Bastian et al. 2012; Messa et al. 2018) with clusters in galactic outskirts or low-mass galaxies disrupting on a longer time-scale (Alvarez-Baena et al. 2024). Numerical work has shown that the response of the cluster to gas-removal and the overall mass-loss can depend on the initial structure of the cluster with respect to the tidal field (Baumgardt & Makino 2003; Baumgardt & Kroupa 2007; Lamers, Baumgardt & Gieles 2010; Smith et al. 2013; Shukirgaliyev et al. 2018). Quantifying the initial sizes and very early evolution of star clusters is thus important for understanding their long-term evolution.

The distribution of gas around young clusters has been used to estimate the time-scale of gas expulsion (Whitmore et al. 2011; Corbelli et al. 2017; Grasha et al. 2018; Messa et al. 2021; Hannon et al. 2022), in order to quantify the role of various stellar feedback processes in determining the initial cluster properties. Recently exposed clusters are often very young (⁠|$< 5$| Myr) and have been associated with gaseous outflows (Levy et al. 2021; Sirressi et al. 2024) conceptually similar to expanding superbubbles (Watkins et al. 2023). With the exposed clusters extending to ages as low as |$\sim 2$| Myr, pre-supernova feedback (H ii regions and stellar winds) has been proposed as an important driver of gas expulsion in young star clusters even before the destructive supernovae (SNe) occur. Counter-examples exist as well: Kim et al. (2023) and Calzetti et al. (2023) find massive, embedded clusters that are 5–6 Myr old, indicating that massive and/or compact enough clusters can resist gas expulsion. The majority of intermediate age and old massive globular clusters (GCs) exhibit chemical variations but only in their light-element abundances (Bastian & Lardo 2018; Gratton et al. 2019). This could be the smoking gun of gas retention in massive clusters that undergo star formation for an extended period of time while massive stars release chemically enriched stellar winds (Krause et al. 2016; Szécsi & Wünsch 2019; Lahén, Naab & Szécsi 2024). What role cluster mass, compactness and metallicity play in the ability of a star cluster to continue forming stars even after the ignition of the first massive stars still remain open questions. He, Ricotti & Geen (2019) have for instance found in molecular cloud simulations that regardless of the initial mass or the compactness of the progenitor gas clouds, they will continue forming stars for a few sound crossing times after which radiation disperses the remaining gas. However, galaxy scale simulations are needed to address the full life cycle of molecular clouds and star clusters (see e.g. the discussion in Jeffreson, Semenov & Krumholz 2024).

Resolving the initial properties of star clusters and the dominant sources of feedback energy on scales of individual stars is difficult in observations of clustered environments. High-resolution hydrodynamical simulations of galactic star formation can shed light on these complex processes that originate from small spatial scales (see e.g. Naab & Ostriker 2017). Individual (massive) stars and their stellar feedback processes have been included in galaxy-scale simulations to study e.g. galactic interstellar medium (ISM) structure and outflows (Hu et al. 2017; Steinwandel et al. 2020; Gutcke et al. 2021; Steinwandel et al. 2024b), star formation efficiency (Hislop et al. 2022), cosmological star formation history and Population III stellar feedback in low-mass galaxies (Gutcke et al. 2022; Sugimura et al. 2024; Andersson et al. 2025), clustering of stars and their feedback in idealised (Lahén et al. 2020; Smith et al. 2021; Deng et al. 2024; Lahén et al. 2024) and cosmological settings (Calura et al. 2022; Garcia et al. 2023; Calura et al. 2024; Gutcke 2024), and the impact of runaway stars on star formation and outflows (Andersson, Agertz & Renaud 2020; Andersson et al. 2023; Steinwandel et al. 2023). Lahén et al. (2023) introduced a numerical model (based on Hu et al. 2017 and Hu 2019 and references therein) within the Galaxy Realizations Including Feedback From INdividual massive stars¹ (griffin) project which considers every newly formed star as an individual particle. The inclusion of stellar evolution models describing the mass, energy, and metal output of massive and very massive stars throughout their lives enables the exploration of their role in the formation and evolution of star clusters star by star in their galactic environment. En route to building a self-consistent picture of how GCs form in a galactic context, Lahén et al. (2024) analysed the early chemical enrichment in proto-GCs and demonstrated that massive clusters can retain gas and recycle stellar wind-material efficiently compared to SN material.

However, one key ingredient missing in previous studies of galactic-scale star formation is the collisional treatment of the individually realized stars, which we address in this paper. Gravitational two-body and multibody interactions drive the internal evolution of star clusters through relaxation, mass-segregation, dynamical binary formation and core collapse (Spitzer & Hart 1971; Spitzer 1987; Goodman & Hut 1993), all of which are suppressed when gravitational softening is used to smooth out the gravitational potential in close encounters. Concurrent modelling of star formation, stellar feedback and collisional N-body dynamics has been successfully demonstrated in simulations of individual giant molecular clouds and star-forming regions (Dinnbier & Walch 2020; Wall et al. 2020; Fujii et al. 2021; Grudić et al. 2021; Rieder et al. 2022).

Collisional simulation codes typically use higher order integration methods such as the fourth-order Hermite (Makino 1991) or the fourth-order forward integrator (Rantala, Naab & Springel 2021) that provide better integration accuracy at the same computational cost compared to e.g. the standard leap-frog algorithm widely used in galactic-scale codes. These fourth-order algorithms are usually, though not always, coupled to regularized few-body solvers. Hydrodynamical simulation codes can be coupled with dedicated N-body solvers such as the direct summation code ph4 (McMillan et al. 2012) or the hybrid PeTaR code (Wang et al. 2020) implemented in the amuse framework (Portegies Zwart et al. 2009; Pelupessy et al. 2013).

Directly summing all two-body forces in a galaxy-scale simulation would however be prohibitively inefficient. Jo et al. (2024) recently implemented the N-body code nbody6++ gpu (Wang et al. 2015) in the enzo code, which they used to solve the collisional evolution of idealized star clusters in a hydrodynamical Milky Way-mass galaxy. In this study, we use the regularized tree code ketju (Rantala et al. 2017; Mannerkoski et al. 2023), which includes a number of numerical algorithms to enable fast and accurate integration of selected simulation regions using the mstar-integrator (Rantala et al. 2020). ketju was first developed in gadget-3 to be primarily used in larger-scale simulations to model the dynamics of supermassive black holes and supermassive black hole binaries in their stellar-dynamical environments (Rantala et al. 2017, 2018; Mannerkoski et al. 2021; Liao et al. 2023; Mannerkoski et al. 2023).

Partmann et al. (2025) introduced recently ketju in the sphgal version of gadget-3 in a new suite of star-by-star griffin-simulations where ketju was used to compute accurate interactions of central massive black holes with surrounding individual stars, while interactions between stars were still softened. Here, we expand the code-implementation of Partmann et al. (2025) by adding the ketju regularized integration regions around all forming massive stars in an entire galaxy. We execute comparison simulations with sphgal-ketju and the direct N-body code bifrost (Rantala et al. 2023; Rantala, Naab & Lahén 2024). Star clusters modelled using ketju expand in size and undergo dynamical mass-loss. We then run hydrodynamical simulations with star formation and an evolving galactic tidal field resembling that of the Wolf–Lundmark–Melotte galaxy as described in Hu et al. (2017). An accurate accounting of gravitational encounters with massive stars produces more realistic star clusters. The clusters form compact and expand rapidly once the first massive stars ignite and expel any remaining gas.

This article is structured as follows. Section 2 provides an overview of the sphgal simulation method and describes briefly the updated interface with the ketju-integrator module. We also give a brief introduction of the direct N-body code bifrost used to validate the sphgal + ketju code. A description of the code comparison tests and galaxy initial conditions is given here as well. Section 3 describes the results of the code comparison, where we investigate the evolution of the size, density, and velocity distribution of idealized star clusters run in isolation. Section 4 discusses the formation, evolution, and disruption of star clusters in hydrodynamical simulations of isolated low-metallicity dwarf galaxies run with the updated sphgal + ketju code, in comparison to simulations that use gravitational softening in all gravitational particle interactions. Our conclusions and final remarks are presented in Section 5.

2 SIMULATIONS

2.1 Hydrodynamical simulation code sphgal

The simulations for this study were run with the sphgal-version (Hu et al. 2014, 2016, 2017; Hu 2019) of the gadget-3 code (Springel 2005). We used the modern smoothed particle hydrodynamics (SPH) implementation in sphgal, which includes several improvements to the numerical accuracy of the SPH methodology as outlined e.g. in Hu et al. (2014, 2016). Recent upgrades (Lahén et al. 2023) to the code include a method for locally mass-conserving sampling of individual stars from a given stellar initial mass function (IMF) and stellar models that describe the evolution of massive stars described briefly below. Here, we expand the methodology for accurate treatment of collisional gravitational dynamics introduced in Partmann et al. (2025) to include high-accuracy ketju (Rantala et al. 2017) integration regions around all massive stars and their remnants. In the simulations, we assume an initial metallicity of |$Z\sim 0.01$|  |$\mathrm{Z}_\odot$| and adopt a fiducial resolution of 4 |$\mathrm{M}_\odot$| in gas and old stellar disc particles and |$6.8\times 10^3$|  |$\mathrm{M}_\odot$| in dark matter. The stars formed during the simulation have masses sampled from an IMF down to the hydrogen burning limit of 0.08 |$\mathrm{M}_\odot$|⁠.

2.1.1 Cooling, chemistry, and star formation

As outlined in detail in Hu et al. (2016) and Hu et al. (2017), cooling and chemistry of the ISM are treated in two temperature regimes. We use a chemical network at low temperatures (⁠|$< 3\times 10^4$| K) and tabulated metal-dependent cooling rates from Wiersma, Schaye & Smith (2009) at high temperatures (⁠|$>3\times 10^4$| K). The chemistry network includes six chemical species (H|$_2$|⁠, H|$^+$|⁠, H, CO, C|$^+$|⁠, O) and free electrons, based on the methods of Nelson & Langer (1997), Glover & Mac Low (2007) and Glover & Clark (2012) closely following the implementation in the silcc-project (Walch et al. 2015). The chemistry network takes into account the spatially resolved interstellar radiation field that is attenuated based on the gas and dust column density as outlined in Section 2.1.2. The code tracks the abundances of 13 individual elements in gas and stars: H, He, N, C, O, Si, Al, Na, Mg, Fe, S, Ca, and Ne.

For star formation, we consider a threshold according to the local Jeans-mass (⁠|$M_\mathrm{J}$|⁠) estimate within the SPH kernel, using a definition for |$M_\mathrm{J}$| as

$$\begin{eqnarray} M_\mathrm{J} = \frac{\pi ^{5/2}c_\mathrm{s}^3}{6G^{3/2}\rho ^{1/2}}, \end{eqnarray}$$

(1)

where |$c_\mathrm{s}$| is the sound speed, |$\rho$| is the SPH-averaged gas density and G is the gravity constant. When |$M_\mathrm{J}$| drops below half of the SPH kernel mass (⁠|$0.5\times 400$|  |$\mathrm{M}_\odot$|⁠), gas particles are turned into star particles. The star particles are first considered to be reservoir particles for a dynamical time according to the local gas density, |$t_\mathrm{dyn}=(4\pi G \rho)^{-1/2}$|⁠. This approximates further gravitational collapse of the parental dense gas phase. In this phase, the stellar reservoir particles are decoupled from the hydrodynamics and interact via gravitational forces only. After one dynamical time, the reservoir particle mass is sampled into individual stars along the Kroupa (2001) IMF between 0.08 |$\mathrm{M}_\odot$| and 500 |$\mathrm{M}_\odot$|⁠. To conserve mass on scales comparable to individual star-forming regions, we perform the IMF-sampling particle by particle considering only the combined mass of the reservoir particles within the Jeans-length |$R_\mathrm{J}$| that was measured at conversion, defined as

$$\begin{eqnarray} R_\mathrm{J} = \left(\frac{3}{4\pi } \frac{M_\mathrm{J}}{\rho } \right) ^{1/3}. \end{eqnarray}$$

(2)

In other words, when the stellar mass sampled per particle exceeds the fiducial resolution of 4 |$\mathrm{M}_\odot$|⁠, the overshoot mass is borrowed from other nearby reservoir particles within |$R_\mathrm{J}$|⁠. If there is not enough reservoir particle-mass nearby, the last stellar mass is re-sampled. This results in under-sampling of massive stars in regions where the local reservoir for sampling is small (low-SFR regions), while regions of intense star formation (e.g. in a starburst) will result in more fully populated input IMFs.

2.1.2 Stellar feedback

The individually sampled stars in the simulation release radiation, energy, momentum, and chemically enriched matter according to their initial stellar masses, metallicities, and evolutionary stages. All stars are assigned with fixed mass-dependent far-ultraviolet luminosities integrated in the range of 6–13.6 eV. For stars below 9 |$\mathrm{M}_\odot$|⁠, we combine the BaSeL spectral library at |$Z\sim 0.01$|  |$\mathrm{Z}_\odot$| (Lejeune, Cuisinier & Buser 1997, 1998; Westera et al. 2002) and the geneva stellar models at |$Z=0.0004\sim 0.02$|  |$\mathrm{Z}_\odot$| (Groh et al. 2019).² The resulting interstellar radiation field within the galaxy, produced by the stellar distribution, is attenuated at the location of each gas particle accounting for dust shielding and gas self-shielding, assuming optically thin conditions. We use the treecol-algorithm (Clark, Glover & Klessen 2012) to compute the attenuation around a gas particle, in 12 solid-angles divided equally using the healpix-algorithm (Górski & Hivon 2011).

Stars more massive than 9 |$\mathrm{M}_\odot$| are supplemented with radiation and stellar wind properties that evolve according to the stellar age. We adopt the Bonn Optimized Stellar Tracks (BoOST; Szécsi et al. 2022) that describe the evolution of stars between 9 and 500 |$\mathrm{M}_\odot$|⁠. We track the far-ultraviolet and photoionizing (⁠|$>13.6$| eV) luminosities, and the stellar winds characterized by the mass-loss rate, velocity and chemical composition. Photoionization is implemented using the Strömgren approximation by setting the gas to be fully ionized at |$10^4$| K within the Strömgren sphere of any photoionizing star, iteratively in case overlapping H ii regions occur (see section 2.5.2 in Hu et al. 2017). Stellar winds are released by injecting momentum and metal-enriched material into |$12\times 8(\pm 2)$| gas particles in 12 equally divided solid angles around the star (see details in Lahén et al. 2023).

In our model, stars with initial masses between 8–40 |$\mathrm{M}_\odot$| and 107.2–203.4 |$\mathrm{M}_\odot$| are assumed to explode as core-collapse and pair-instability SNe, respectively. The latter corresponds to stars with final helium-core masses of 65–133 |$\mathrm{M}_\odot$| that are expected to undergo thermonuclear runaway caused by the pair-production instability (Heger & Woosley 2002). In practice, the low-SFR galaxies considered in this study do not produce any pair-instability SN progenitors due to the locally mass-conserving IMF-sampling method. Core-collapse SNe are supplemented with thermal explosion energies of |$10^{51}$| erg, and metal yields and remnant masses from Chieffi & Limongi (2004), spatially distributed in a similar way as the stellar wind using healpix. Stars between the range of 40 |$\mathrm{M}_\odot$| and 107.2 |$\mathrm{M}_\odot$| are assumed to collapse to black holes without SNe.

2.1.3 Regularised integrator ketju

In order to accurately model the gravitational interactions at small separations between stars in the young star clusters without gravitational softening, we utilise the ketju integration technique (Rantala et al. 2017) implemented originally in the gadget-3 code. The most recent version of ketju (Mannerkoski et al. 2023) and its regularized integrator library mstar (Rantala et al. 2020) are publicly available,³ as a module for the gadget-4 code (Springel et al. 2021). The gadget-3-ketju code has been used in a wide range of studies examining the dynamics of supermassive black holes in isolated mergers of early-type galaxies (Rantala et al. 2018, 2019), dark matter haloes, and low-mass galaxies (Partmann et al. 2024), and gas-rich galaxies (Liao et al. 2023, 2024a, b) as well as in cosmological simulations of forming early-type galaxies (Mannerkoski et al. 2021, 2022). We use the gadget-3-version of ketju as described in Mannerkoski et al. (2021) as it is directly compatible with sphgal.

The key idea of ketju is to select regions (hereafter ketju-regions) with radius |$r_\mathrm{{\small KETJU}}$| up to few tens of parsecs in size within the gadget-3 simulation and to integrate the gravitational dynamics within the regions at high accuracy without the need for gravitational softening. The regularized integrator mstar used in the regions is based on three key numerical methods that enable efficient and accurate integration as follows:

• Algorithmic regularization (Mikkola & Tanikawa 1999; Preto & Tremaine 1999; Mikkola & Merritt 2006, 2008) using time-transformed equations of motion together with a custom leapfrog integrator for efficient non-softened integration of N-body systems avoiding the Newtonian coordinate singularity at small particle separations.

• Minimum spanning tree coordinate system (Rantala et al. 2020) to minimize floating-point round-off error and thus improve the integration accuracy.

• Gragg–Bulirsch–Stoer (Gragg 1965; Bulirsch & Stoer 1966) extrapolation method to reach extremely high user-desired integration accuracy.

mstar is efficiently parallelized and can treat up to thousands of simulation particles per regularized region. The library version of mstar also includes an option for post-Newtonian equations of motion and a method for sophisticated integration order control (Mannerkoski et al. 2021).

In previous studies the ketju-regions have been placed around intermediate-mass or supermassive black holes in simulations with a stellar particle mass resolution of |$\sim 0.08$|–|$10^5\, \mathrm{M}_\odot$|⁠. However, ketju allows for placing the regions at arbitrary locations and centring on very massive objects is not required for the integration technique to work. In this study for the first time, we use ketju regions around massive stars in a star-by-star mass resolution dwarf galaxy simulation. More specifically, we include ketju regions around every star above a certain mass threshold |$m_i$| to capture the collisional dynamics of these stars without gravitational softening. Outside the regions the code time integration is second-order accurate and all gravitational interactions are softened with the softening lengths |$\epsilon$| of the particles.

As the ketju-region size has to be at least 2.8 times the gravitational softening length of the particles, we adopt here a factor of 3, i.e. |$r_\mathrm{{\small KETJU}}=3\epsilon _*$| where |$\epsilon _*$| is the stellar softening length. ketju-regions that overlap are merged into one, thus the central regions of star cluster can host ketju-regions that are significantly larger than the single region size. The dynamics of low-mass stars below the threshold mass |$m_i$| remains softened unless the low-mass stars are within a regularized region of a nearby massive star. The reservoir stars (before IMF-sampling is done), old disc particles, gas particles, and dark matter particles are also not included in the ketju calculations. All interactions that are not between stars in a ketju-region are computed using the standard leapfrog integrator and tree force algorithm in gadget-3. We do not use the optional star–star softening inside the ketju-regions since our N-body particles represent actual physical stars instead of macro-particles as has been the case in most previous ketju studies. For the user-given accuracy parameters we use a Gragg–Bulirsch–Stoer tolerance of |$\eta _\mathrm{GBS}=10^{-7}$| and an end-time iteration tolerance of |$\eta _\mathrm{t}=10^{-3}$|⁠. The simulations in this study are fully Newtonian.

2.1.4 Code comparison: bifrost

In this study we use the direct N-body code bifrost for comparison star cluster simulations to verify the improved accuracy of our gadget-3-ketju star cluster dynamics. The benchmark code bifrost (Rantala et al. 2023; Rizzuto et al. 2023; Rantala & Naab 2024; Rantala et al. 2024) is a modern GPU-accelerated direct-summation N-body simulation code based on the hierarchical (Rantala et al. 2021) fourth-order forward symplectic integrator (Chin 1997; Chin & Chen 2005; Chin 2007; Dehnen & Hernandez 2017). Besides the forward integrator, bifrost uses a number of secular and regularized few-body integration techniques closely related to mstar (Rantala et al. 2020) to solve the gravitational dynamics of binary stars, close fly-bys, multiple systems and small clusters around massive black holes. In bifrost, no gravitational softening is used. For additional details of the code, see Rantala et al. (2023) and appendix A of Rantala et al. (2024).

2.2 Simulation set-ups and initial conditions

2.2.1 N-body simulations

In order to compare our sphgal + ketju code with the direct N-body code bifrost, we run two sets of isolated star cluster simulations without stellar evolution. The cluster initial conditions have Plummer density profiles and the half-mass radii |$r_{50~{{\rm per\ cent}}}$| are set according to the observed mass–size relation of Marks & Kroupa (2012). Individual stellar masses are sampled from a Kroupa IMF, adopting the cluster mass-dependent upper mass limit to stellar masses according to Yan, Jerabkova & Kroupa (2023), see also Rantala et al. (2024) for details. The first set follows the evolution of a star cluster consisting of |$\sim 10^3$| stars and a cluster mass of |$\sim 600$|  |$\mathrm{M}_\odot$| with |$r_{50~{{\ \rm per\ cent}}}=0.18$| pc. The second set of simulations considers the evolution of a cluster with |$\sim 10^4$| stars and a mass of |$\sim 6000$|  |$\mathrm{M}_\odot$| with |$r_{50~{{\ \rm per\ cent}}}=0.28$| pc. We created 10 random realizations of both initial conditions. The sample of runs is detailed in Table 1 and further information about the simulation parameters is given in Section 2.3. We label the simulations executed with bifrost as bifrost, simulations run while ketju is enabled as ketju and the simulations where ketju is disabled as gadget.⁴

Table 1.

Description of the N-body simulations. The columns give the simulation name, the number of stars subject to ketju-integration, the number of random initializations, the total simulation time, the code used, the stellar gravitational softening length, the size of the regularized region and the initial stellar mass of particles used as centres of ketju regions. Run hyd_KETJU_m3_e001 was limited to a total runtime of 400 Myr due to the increased computational cost.

Name	|$N_*$|	|$N_\mathrm{sample}$|	Time (Myr)	Code	|$\epsilon _*$|	|$r_\mathrm{{\small KETJU}}$| (pc)	ketju  |$m_i$| (⁠|$\mathrm{M}_\odot$|⁠)
nb_1k_BIFROST	1155	10	50	bifrost	–	–	–
nb_1k_KETJU_m0_e001	1155	10	50	sphgal + ketju	0.01	0.03	0.08
nb_1k_KETJU_m3_e001	1155	10	50	sphgal + ketju	0.01	0.03	3
nb_1k_KETJU_m3_e01	1155	10	50	sphgal + ketju	0.1	0.3	3
nb_1k_KETJU_m8_e001	1155	10	50	sphgal + ketju	0.01	0.03	8
nb_1k_GADGET_e001	1155	10	50	sphgal	0.01	–	–
nb_10k_BIFROST	10 660	10	50	bifrost	–	–	–
nb_10k_KETJU_m3_e001	10 660	10	50	sphgal + ketju	0.01	0.03	3
nb_10k_KETJU_m8_e001	10 660	10	50	sphgal + ketju	0.01	0.03	8
nb_10k_SPHGAL_e001	10 660	10	50	sphgal	0.01	–	–
hyd_KETJU_m3_e01	513 794	1	500	sphgal + ketju	0.1	0.3	3
hyd_SPHGAL_e01	489 226	1	500	sphgal	0.1	–	–
hyd_KETJU_m3_e001	478 024	1	400	sphgal + ketju	0.01	0.03	3
hyd_SPHGAL_e001	487 158	1	500	sphgal	0.01	–	–

2.2.2 Dwarf galaxy simulations

The low-metallicity (⁠|$Z=0.01$|  |$\mathrm{Z}_\odot$|⁠) dwarf galaxy initial conditions are adopted from Lahén et al. (2023), based on the dwarf galaxy models of Hu et al. (2016). We use the compact initial condition with virial mass of |$4\times 10^{10}$|  |$\mathrm{M}_\odot$|⁠, including a |$4\times 10^7$|  |$\mathrm{M}_\odot$| gas disc and a |$2\times 10^7$|  |$\mathrm{M}_\odot$| old stellar disc with disc scale lengths of 0.73 kpc. The initial mass resolution is 4 |$\mathrm{M}_\odot$| for gas and disc stars and |$\sim 6800$|  |$\mathrm{M}_\odot$| for dark matter. The gravitational force-softening length is set to 62 pc for dark matter, 0.1 pc for gas and 0.1 pc for pre-existing disc star particles. The smallest SPH smoothing lengths containing the 100 nearest gas particles are |$\sim 0.3$| pc, corresponding to densities of |$10^5$| cm|$^{-3}$| at the typical star formation threshold. For newly formed stars we use |$\epsilon _*$| of 0.1 pc or 0.01 pc (see Section 2.3 for details) except when for stars that are within ketju regions in the ketju simulations. The hydrodynamical simulations are named as hyd_KETJU when ketju is enabled and as hyd_SPHGAL when ketju is disabled.

2.3 Softening and regularization parameters

We can estimate the desired minimum size for the ketju region by calculating the typical strong encounter impact parameter or |$90^\circ$| deflection radius |$b_{90}$| (Binney & Tremaine 2008) for two equal mass stars |$m_*$| as

$$\begin{eqnarray} b_{90} = \frac{2Gm_*}{V^2} \end{eqnarray}$$

(3)

with the encounter velocity |$V\sim \sigma$|⁠, the stellar velocity dispersion. Another possible limit for the region size is the transition from a hard to a soft binary, defined as the limit where the potential energy of a binary equals the mean kinetic energy of a particle in the surrounding cluster, which for an equal-mass binary becomes

$$\begin{eqnarray} a_\mathrm{bin} \sim \frac{Gm_*^2}{\left< m_*\right>\sigma ^2}, \end{eqnarray}$$

(4)

where |$a_\mathrm{bin}$| is the binary semimajor axis. As we are interested mainly in the evolution of massive stars, we can assume that they segregate rapidly to the central region of the cluster resulting in |$\left< m_*\right>\sim m_*$| at the centre and yielding |$a_\mathrm{bin} \sim b_{90}$|⁠.

The 1D central velocity dispersions of the |$\sim 600$|  |$\mathrm{M}_\odot$| and |$\sim 6000$|  |$\mathrm{M}_\odot$| clusters with Plummer profiles are |$\sigma _\mathrm{0,1k}=1$| km s|$^{-1}$| and |$\sigma _\mathrm{0,10k}=2.6$| km s|$^{-1}$|⁠, respectively. The largest stellar masses are |$\sim 15$|  |$\mathrm{M}_\odot$| and |$\sim 55$|  |$\mathrm{M}_\odot$|⁠, respectively. Using equations (3) and (4), we get |$b_{90,1k}\sim 0.06$| pc and |$b_{90,10k}\sim 0.04$| pc. Compact (⁠|$r_{50~{{\ \rm per\ cent}}}<1$| pc) bound young star clusters that form in the isolated dwarf galaxies have typical masses in the range from a few tens to |$\sim 900$|  |$\mathrm{M}_\odot$| and typical 1D velocity dispersions of 0.1–2.5 km s|$^{-1}$|⁠, resulting in typical |$b_{90}$| values of less than |$\sim 0.1$| pc.

Optimal |$r_\mathrm{{\small KETJU}}$| values for strong gravitational interactions to be treated with the regularized ketju integrator should thus be larger than 0.01–0.1 pc for our simulations. In addition to defining which interactions are computed with the ketju-integrator, |$r_\mathrm{{\small KETJU}}$| also controls the computational cost of the integrator due to the scaling of the direct N-body problem with the particle number (N) as |$\sim N^2$|⁠. We do not therefore wish to make the ketju-regions larger than necessary. In an idealized environment of constant stellar density, the total N in a ketju-region scales as |$N\propto r_\mathrm{{\small KETJU}}^3$|⁠. In case dense systems form, a large value of |$r_\mathrm{{\small KETJU}}$| can lead to an increased computational load in ketju.⁵ A large value of |$r_\mathrm{{\small KETJU}}$| would additionally allow close ketju-regions to be combined as one, further increasing the number of particles integrated in one region. We thus test values of |$\epsilon _*$| of 0.01 and 0.1 pc in both N-body and hydrodynamical runs, which define |$r_\mathrm{{\small KETJU}}$| as 0.03 or 0.3 pc, respectively. These ketju-regions sizes bracket the optimal values estimated using equations (3) and (4).

We run a sample of ten random initializations of both idealized cluster initial conditions (⁠|$10^3$| and |$10^4$| stars) with bifrost; with sphgal (that is, without hydrodynamics or stellar evolution, the standard gadget-3) using a |$\epsilon _*=0.01$| pc; and with sphgal + ketju using |$\epsilon _*=0.01$| pc and initiating ketju integration within |$r_\mathrm{{\small KETJU}}=0.03$| pc either around all |$m_i>3$|  |$\mathrm{M}_\odot$| or around all |$m_i>8$|  |$\mathrm{M}_\odot$| stars. In addition, we run the 1k initial conditions with sphgal + ketju using |$r_\mathrm{{\small KETJU}}=0.3$| pc around |$m_i>3$|  |$\mathrm{M}_\odot$| stars and |$r_\mathrm{{\small KETJU}}=0.03$| pc around all stars. We do not execute the 10k-sample with |$r_\mathrm{{\small KETJU}}=0.3$| pc because it includes initially more than 5700 stars in the central 0.3 pc, and running it with the larger ketju-region size is prohibitively expensive for tens of Myrs (Rantala et al. 2020). The simulations are run for 50 Myr, which is significantly longer than the half-mass relaxation times of 1.4 and 5.7 Myr of the 1k and 10k initial conditions, respectively. The labels and parameter details of the simulations are collected in Table 1.

In the hydrodynamical runs, we test both |$\epsilon _*=0.1\, \mathrm{pc}$| and |$\epsilon _*=0.01\, \mathrm{pc}$|⁠, i.e. |$r_\mathrm{{\small KETJU}}=0.3$| and |$r_\mathrm{{\small KETJU}}=0.03$| pc, and run the simulations with and without ketju. As clusters more massive than |$\sim 100$|  |$\mathrm{M}_\odot$| are practically guaranteed in our simulations to form with stars that are more massive than 3 |$\mathrm{M}_\odot$| (see e.g. Lahén et al. 2023, fig. 14), we initiate the regions around |$m_i>3$|  |$\mathrm{M}_\odot$| stars to allow the majority of the massive clusters to have at least one ketju-region for part of their early evolution. ketju is triggered in total for |$\sim 7000$| stars during both regularized hydrodynamical simulations and the maximum number of stars within the individual regularized regions is 1000–2000 in hyd_KETJU_m3_e01 and 100–200 in hyd_KETJU_m3_e001.

2.4 Identification of star clusters

In the hydrodynamical dwarf galaxy simulations, we identify bound star clusters using friends-of-friends and subfind (Springel et al. 2001; Dolag et al. 2009), which are structure finding algorithms included in gadget-3. In order to follow the long-term evolution of clusters especially in the ketju-runs, we use a linking length of 1 pc to recover the extended distribution of stars. The subfind-algorithm has been modified to account for the non-softened gravitational potential within the ketju-regions in the runs where ketju is enabled. We require at least 50 bound stars to be counted as a cluster, and exclude clusters with a mean stellar age younger than 1 Myr when making direct comparisons to observations where deeply embedded, extremely young clusters are often excluded. The cluster finding is performed for each snapshot over a time-span of 500 in 1 Myr intervals.

3 COMPARISON OF KETJU WITH DIRECT N-BODY

The number of interactions integrated with ketju depends on the number and sizes of the regularized regions. Stars with initial masses of |$m_i>3$|  |$\mathrm{M}_\odot$| and |$m_i>8$|  |$\mathrm{M}_\odot$| account typically for |$\sim 2~{{\rm per\ cent}}$| and 0.5  per cent of all stars in the idealized cluster, respectively. With |$r_\mathrm{{\small KETJU}}=0.3$| pc and |$m_i=3$|  |$\mathrm{M}_\odot$| in nb_1k_KETJU_m3_e01, initially all interactions between stars inside |$r_{50~{{\ \rm per\ cent}}}$| are integrated with ketju as one combined ketju-region can cover the entire volume within |$r_{50~{{\ \rm per\ cent}}}$|⁠. The fraction of regularized gravity interactions reduces as mass-segregation proceeds and the clusters expand, dropping to a few tens of per cent during the first 10–20 Myr of evolution. When the size of the regularized region is reduced to |$r_\mathrm{{\small KETJU}}=0.03$| pc, the fraction of regularized interactions within |$r_{50~{{\ \rm per\ cent}}}$| is initially |$\sim 10~{{\ \rm per\ cent}}$| and |$\sim 20~{{\ \rm per\ cent}}$| in the nb_1k_KETJU_m3_e001 and nb_10k_KETJU_m3_e001 runs, respectively. The fraction reduces gradually to a few per cent as massive stars are segregated and/or ejected out of the cluster (since stellar evolution is switched off here). The fraction of stars within ketju regions is close to the fraction of |$m_i>3$|  |$\mathrm{M}_\odot$| stars, indicating that most of these stars have moved on to more distant orbits with low stellar densities or escaped the cluster and evolve in isolation. When |$m_i$| is increased to 8 |$\mathrm{M}_\odot$|⁠, the nb_1k_KETJU_m8_e001-run has only a handful of stars above the mass threshold. The initial fraction of stars within ketju-regions is (averaged over 10 runs) only a few per cent and plateaus at half a per cent during the simulations, close to the fractions of |$m_i>8\, \rm {\rm M}_{\odot }$| stars. The nb_10k_KETJU_m8_e001-run has a |$\sim 7~{{\ \rm per\ cent}}$| initial fraction of regularized interactions within |$r_{50~{{\ \rm per\ cent}}}$|⁠, which drops to 1–2  per cent after 10 Myr of evolution.

As the massive stars segregate to the cluster centre, they can form binaries, undergo strong gravitational two or few-body interactions and gradually get ejected from the clusters. Low-mass clusters simulated with ketju can thus end up being integrated with fully softened gravity if all of the relatively rare massive stars escape. In this case, the size evolution of the cluster halts entirely. This is essentially the cause of the plateau in the inner Lagrangian radii after 25 Myr in the nb_1k_KETJU_m8_e001-sample discussed in Section 3.1.

3.1 Lagrangian radii

Fig. 1 shows the evolution of the Lagrangian radii enclosing 5 per cent, 50 per cent, and 90 per cent of the total mass of the clusters for the N-body simulations during 50 Myr. The nb_1k_BIFROST and nb_10k_BIFROST runs are the benchmarks against which the sphgal and sphgal + ketju simulations are compared to.

Figure 1.

The 5 per cent (top), 50 per cent (middle), and 90 per cent (bottom) Lagrangian radii in the 1k (left) and the 10k (right) simulations, run in isolation without stellar evolution. The lines show the mean values of the sample of 10 runs per parameter combination, and the shaded regions show the standard deviations of the samples. The simulations on the left are coloured as follows: nb_1k_BIFROST in dashed black, nb_1k_KETJU_m0_e001 in solid black, nb_1k_KETJU_m3_e001 in blue, nb_1k_KETJU_m3_e01 in dashed blue, nb_1k_KETJU_m8_e001 in magenta and nb_1k_GADGET_e001 in orange, and equivalently for the 10k-runs on the right. The sizes evolve increasingly similar to the direct N-body simulations (bifrost) in simulations that solve increasing fractions of close gravitational interactions without softening. However, even a small fraction of accurately solved interactions (e.g. nb_1k_KETJU_m8_e001 and nb_10k_KETJU_m8_e001) already drastically improves the cluster evolution compared to fully softened simulations (nb_1k_GADGET_e001 and nb_10k_GADGET_e001).

The fully softened simulations nb_1k_GADGET_e001 and nb_10k_GADGET_e001 show practically no evolution in the Lagrangian radii, especially when compared to the (partially) collisional simulations. Clusters in the softened gravity would only lose mass if the clusters were exposed to an external tidal field. To check whether changing the softening length impacts the cluster size, we ran the sphgal-simulations without ketju using a larger |$\epsilon _*=0.1$| pc softening and saw no changes compared to the |$\epsilon _*=0.01$| pc run. We thus omit the simulation with |$\epsilon _*=0.1$| from further analysis.

The clusters in the simulations with regularization start compact and their centres contract even further during the first few Myr while mass-segregation proceeds. In response, the outer regions expand, up to a factor of 10 and 100 for the 50 per cent and 90 per cent radii, respectively. The inner regions of individual nb_10k_BIFROST-runs and nb_10k_KETJU-runs go through cycles of gradual contraction (see Fig. A1 in Appendix 5) that is halted by binary formation and rapid expansion of |$r_{5~{{\ \rm per\ cent}}}$|⁠. The expansion is then finished through the ejection of massive star(s) from the centre. However, this is not visible after the first core collapse cycle (around |$\sim 1$| Myr in the bifrost run) when the cluster size is averaged over 10 runs. This is because each random simulation realization has a different time-scale for the collapse-cycle depending on the masses of the most massive stars still present in the cluster. Overall, the mean 5 per cent Lagrangian radius increases after the first few Myr, as the clusters have too low masses to undergo further core collapse (see Fig. A1 and Heggie & Ramamani 1989).

The simulations that either integrate all close interactions with ketju (nb_1k_KETJU_m0_e001) or include a larger fraction of the cluster stars in the ketju integration (nb_1k_KETJU_m3_e01) evolve almost identically to the direct N-body runs. The results of the other simulations that use ketju but only include stars above a mass threshold as the ketju region centres (KETJU_m3, KETJU_m8) are in between the direct N-body results and the fully softened runs. This signifies that the dynamical processes important for the internal evolution of the clusters, such as mass-segregation and core collapse, are captured significantly better in all runs with ketju compared to fully softened runs. Comparing nb_10k_KETJU_m3_e001 and nb_10k_KETJU_m8_e001, we see that the dynamical evolution of the clusters is not very sensitive to the exact number of ketju regions. In other words, even only accounting for the collisional evolution of a few per cent of the stars in the cluster results in significantly more realistic cluster evolution.

3.2 Stellar density and velocity

Fig. 2 shows the density and velocity distributions of stars in the N-body simulations at 50 Myr, together with the initial conditions. We indicate the escape velocity at |$r_{50~{{\ \rm per\ cent}}}$|⁠,

$$\begin{eqnarray} v_\mathrm{esc}=\sqrt{2|\Phi _{50~{{\ \rm per\ cent}}}|}, \end{eqnarray}$$

(5)

where |$\Phi _{50~{{\ \rm per\ cent}}}$| is the gravitational potential at |$r_{50~{{\ \rm per\ cent}}}$|⁠, in the initial conditions and at 50 Myr. The sample-averaged |$v_\mathrm{esc}$| is between |$\sim 1$| and 5 km s|$^{-1}$| for the 1k runs and |$\sim 4$|–13 km s|$^{-1}$| for the 10k-runs. We also show separately stars that are isolated by at least 1 pc from any other star, which comprise walkaway (10 km s|$^{-1} <$|  v  |$< 30$| km s|$^{-1}$|⁠) and runaway (v  |$>30$| km s|$^{-1}$|⁠) stars with velocities larger than the typical escape velocity of star clusters. Some of the stars have large velocities as they reside in binaries. Thus, any binary members are excluded from the isolated star analysis to avoid confusing binary motion with escaping stars. As stellar evolution is not included in the N-body comparison runs, escapes e.g. due to exploding binary companion stars do not occur, and the stars reach their escape velocities purely due to few-body interactions.

Figure 2.

The radial stellar density (top), the stellar velocity distribution (middle), and the velocity of isolated stars (zero stars within a radius of 1 pc; bottom) in the final snapshot (50 Myr) of the 1k (left) and 10k (right) N-body runs. The line colours are as in Fig. 1 and the lines show the mean values of the sample of 10 runs per parameter combination. The density and velocity distribution in the initial conditions is shown in grey thick lines in the top and middle rows, and the initial escape velocity within |$r_{50~{{\ \rm per\ cent}}}$| is plotted with a thin grey vertical line in the middle panels. The escape velocity within |$r_{50~{{\ \rm per\ cent}}}$| in the final snapshot of each run is indicated with the vertical line in the same linestyle. Star clusters evolved with bifrost or ketju evolve towards energy equipartition, expand, and produce escaping stars, while clusters simulated with fully softened gravity evolve very little from the initial condition.

As with the Lagrangian radii, the nb_1k_GADGET_e001 and the nb_10k_GADGET_e001 runs do not evolve much from the initial condition, neither in the density nor the velocity distribution. The density distribution expands very little, and the stars maintain almost the same velocity distributions throughout the 50 Myr of evolution. As two-body interactions are softened in the dense central region of the cluster, energy equipartition cannot proceed, mass segregation is suppressed and no escaping stars are produced due to the lack of strong few-body interactions.

On the other hand, the simulations that do capture at least some of the close gravitational encounters show significant expansion during the evolution of the clusters (the top panels of Fig. 2). During 50 Myr of evolution, the central densities decrease and the outer regions expand. The effect is stronger in the runs that treat a larger fraction of gravitational interactions accurately without softening. As we saw in Fig. 1, the expansion of the clusters is slower in the ketju-runs compared to bifrost, because interactions between low-mass stars outside the regularized regions are still softened and their segregation is thus somewhat suppressed.

As mass-segregation proceeds, interactions lead to the redistribution of stars on to more distant (less bound) orbits with lower average orbital velocities while a handful of stars reach velocities above the escape velocity. Since the strongest interactions occur in the cluster centre and often include massive stars, we see more consistent numbers of escaping stars being produced in simulations that have at least some regularized interactions (bifrost-runs and ketju-runs). The runs with ketju-regions around all stars (nb_1k_KETJU_m0_e001) or in a larger volume (nb_1k_KETJU_m3_e01) show almost identical results to the nb_1k_BIFROST run, both in density and in the velocity distributions. The runs in the nb_10k_BIFROST-sample, on the other hand, have |$\sim 90$| stars with velocities larger than 10 km s|$^{-1}$| at 50 Myr, compared to the equivalent value of |$\sim 120$| and |$\sim 130$| in the nb_10k_KETJU_m3 and nb_10k_KETJU_m8 samples. The cause for the larger number of runaway stars in the ketju simulations are the higher central densities that enable an elevated number of strong gravitational interactions. Here it is worth stressing again that the clusters evolved entirely with softened gravity produce practically no escaping stars due to the suppression of strong gravitational encounters, even though they have the highest stellar densities for the majority of the simulation time.

4 STAR FORMATION AND CLUSTER EVOLUTION ON GALACTIC SCALES

Next we turn to the full hydrodynamical galaxy-simulations that follow the formation, evolution, and disruption of star clusters in the evolving tidal field of a low-metallicity dwarf galaxy. In Fig. 3, we show the surface density maps of stars and gas, as well as the gas temperature and thermal pressure, at the end of the hydrodynamical hyd_KETJU_m3_e01 and hyd_SPHGAL_m3_e01 simulations at 500 Myr. The gaseous distributions are fairly similar in the two runs, given that they are sensitive to the temporal stochasticity of the feedback produced by young stars mainly on small spatial scales. Bubbles of photoionization and SNe (e.g. the red blobs in the temperature panels) are visible in all gaseous quantities. The sphgal-panels show good examples of SNe that have exploded within photoionized bubbles that are visible as concentric expansion fronts especially in the pressure panel.

Figure 3.

The stellar and gaseous surface density (top left and right), the gas temperature (bottom left), and the thermal gas pressure (bottom right) in the hyd_KETJU_m3_e01 (four left panels) and hyd_SPHGAL_m3_e01 (four right panels) simulations at 500 Myr. The tidal tails of star clusters that are losing mass can be seen as leading and trailing arms around the concentrations of stellar mass on the left, while clusters on the right remain compact and only lose a little mass in tidal tails. SN-bubbles and photoionized regions are visible in the gaseous distribution in both simulations. The field of view is 2 kpc and the image resolution is |$\sim 4$| pc per pixel.

On the other hand a clear difference between the two simulations can be seen in the distribution of stars. The regularised ketju-simulation shows clear signs of cluster disruption through strong leading and trailing tidal tails associated with individual star clusters. The stellar streams with no clear concentrations of stellar mass are fully disrupted clusters. Some mass-loss is seen in the sphgal-run as well, enabled by the cluster interacting with the tidal field. The clusters, however, still remain very compact when all gravitational forces are softened (see Hislop et al. 2022, for a detailed discussion).

4.1 Star formation rate and cluster formation rate

We first inspect the impact of adding collisional stellar dynamics to hydrodynamical galaxy simulations by comparing the global properties of the simulated galaxies in the sphgal and sphgal + ketju runs. Fig. 4 shows the galaxy-wide SFRs, as well as the star cluster formation efficiencies |$\Gamma$| as defined by the ratio between the cluster formation rate (CFR) and the SFR measured over a similar timespan. The SFRs shown in the top panel are averages over the past 10 Myr, accompanied by the corresponding overall average values. CFRs are computed using the total mass in young (1–10 Myr, middle panel) or intermediate age (10–100 Myr, bottom panel) bound clusters more massive than 100 |$\mathrm{M}_\odot$|⁠, and are then divided by the average SFR over 10 and 100 Myr, respectively, to arrive at the value of |$\Gamma$|⁠. We show both the e01 and e001 versions of the hydrodynamical simulations in Fig. 4, and conclude that there are very little differences in the SFRs (top panel) between the runs with and without ketju or in the runs with varying |$\epsilon _*$|⁠, as indicated by the near identical time-averages over the 400–500 Myr of galaxy evolution.

Figure 4.

The SFR averaged in 10 Myr bins (top), and |$\Gamma$| computed as the ratio between the CFR and SFR for bound star clusters between ages of 1–10 Myr (middle) and 10–100 Myr (bottom). The thin lines show the time variation for hyd_KETJU_m3_e01 (black solid), hyd_KETJU_m3_e001 (black dashed), hyd_SPHGAL_m3_e01 (orange solid) and hyd_SPHGAL_m3_e001 (orange dashed), and the thick horizontal lines show the respective time-averages.

The value of |$\Gamma$| (middle panel), on the other hand, decreases with the inclusion of more accurate stellar dynamics. As already indicated in Fig. 2, collisional interactions cause a more rapid expansion and thus faster tidal disruption of clusters. Interactions with gas clouds and other clusters can also cause mass-loss events. Some of this evolution is captured in the sphgal simulations as well, as shown already in Lahén et al. (2023) where star clusters modelled with the standard sphgal lost tens of per cents of their mass over hundreds of Myrs even with a 0.1 pc gravitational softening.

In Fig. 4, the sphgal runs have time-averaged |$\Gamma$| (1–10 Myr) of 40–50  per cent, while the ketju-runs result in |$\Gamma$| (1–10 Myr) |$\sim 30$|–35 per cent. The global stellar mass and the number of gravitationally bound star clusters formed in the simulations are almost identical between the sphgal and ketju runs, but the ketju-clusters begin losing bound stars dynamically immediately at formation. This is in addition to the stellar wind mass-loss, which happens in the sphgal-run. The values of CFR decrease when averaged over 1–10 Myr, reducing also the 10 Myr averaged |$\Gamma$| when collisional dynamics are included. We will discuss the evolution of the individual clusters in more detail in Section 4.5. We note the persistent variation of |$\Gamma$| from snapshot to snapshot caused by stochasticity in star formation that is active in only a few regions of the galaxy at any given time. Observations of |$\Gamma$| in dwarf galaxies show similar variations from galaxy to galaxy, with sample mean or median values of 10–40 per cent and a standard deviation up to 20 per cent in young (1–10 Myr) clusters (Chandar et al. 2023; Cook et al. 2023). However, more important for the current study is the change of |$\Gamma$| with cluster age; observed values of |$\Gamma$| drop to a few per cent when cluster formation in the age range of 10–100 Myr is analysed. In Fig. 4, only the ketju-simulations recover 10 per cent-level values of |$\Gamma$| for older clusters, caused by more efficient mass-loss.

The clusters in the idealized nb_1k_GADGET_e001 and nb_10k_GADGET_e001 runs showed no size-evolution in Fig. 1. In the bottom panels of Fig. 4, we however see that the impact of the galactic tidal field in the hydrodynamical simulations on the dynamical mass-loss of the clusters is stronger when a smaller softening length is used in hyd_SPHGAL_m3_e001. The sphgal-run with a smaller |$\epsilon _*$| shows some reduction in |$\Gamma$| especially at older cluster ages caused by cluster expansion and mass-loss. The evolution is, however, still too slow to reach mass-loss rates equivalent to the ketju-runs. Based on Fig. 4, we conclude that the evolution of clusters in the softened galaxy-scale sphgal-simulations is sensitive to the specific adopted value of the softening length. The opposite is true for the ketju-simulations where an order of magnitude change in the ketju-region size leads to qualitatively similar results.

In the following, we will therefore concentrate on analysing the differences between the hyd_SPHGAL_m3_e01 and hyd_KETJU_m3_e01 runs.

4.2 Chemical enrichment and galactic outflows

The gas and stellar particles include information about how much each stellar feedback process contributed to the chemical composition of the particle. This enables us to trace the propagation of the wind and SN-material throughout the galaxy. In Fig. 5, we show the time-evolution of the outflow properties and mass-loading factors in the hyd_KETJU_m3_e01 and hyd_SPHGAL_m3_e01 simulations. The outflows have been measured at a height of 1 kpc above and below the galaxy mid-plane in 100 pc thick slabs.⁶ We show the total gas outflow rate (⁠|$\dot{M}_\mathrm{outflow}$|⁠); outflow rate of SN and wind-material (⁠|$\dot{M}_\mathrm{SN,outflow}$| and |$\dot{M}_\mathrm{Wind,outflow}$|⁠); the global SFR; and the locking rate of SN and wind-material into new stars (⁠|$\dot{M}_\mathrm{SN,*}$| and |$\dot{M}_\mathrm{Wind,*}$|⁠). The loading factors |$\eta$| have been computed as the ratio between the respective rate of outflow and star formation, which removes the material from the ISM and locks it into new stars: |$\eta _\mathrm{gas}=\dot{M}_\mathrm{outflow}/\mathrm{SFR}$|⁠, |$\eta _\mathrm{SN}=\dot{M}_\mathrm{SN,outflow}/\dot{M}_\mathrm{SN,*}$| and |$\eta _\mathrm{Wind}=\dot{M}_\mathrm{Wind,outflow}/\dot{M}_\mathrm{Wind,*}$|⁠. We also compute the cumulative values to compare the total masses in formed stars, the total chemical enrichment and the total mass-outflow out of the galaxy in the two simulations.

Figure 5.

Top: The galactic outflow rate of gas, SN-material and wind-material, and the rate at which wind and SN-material are locked into new stars (left) and the respective cumulative values (middle). The lines from top to bottom are |$\dot{M}_\mathrm{outflow}$|⁠, SFR, |$\dot{M}_\mathrm{SN,outflow}$|⁠, |$\dot{M}_\mathrm{SN,*}$|⁠, |$\dot{M}_\mathrm{Wind,outflow}$|⁠, and |$\dot{M}_\mathrm{Wind,*}$| (see the text for details). The hyd_KETJU_m3_e01 and hyd_SPHGAL_m3_e01 simulations are shown in solid and dashed line style. The right column shows the ratios between the cumulative masses in the ketju-run and the sphgal-run. The dashed horizontal line indicates the value of one. Bottom: The mass-loading factor of gas (black lines), SN-material (yellow lines), and wind-material (white lines). The panels show the rates and loading factors (left), the cumulative masses and loading factors (middle) based on the cumulative values from the top middle panel, and the ratios between the cumulative loading factors in the ketju and sphgal simulations (right).

Based on Fig. 5, the efficient cluster dissolution in such a low-SFR galaxy has only a minor impact on the galaxy-wide baryon cycle. The total mass in outflows, formed stars and locked metals show similar values in the ketju and sphgal runs to within a few tens of per cent. The massive stars in the ketju simulation are less clustered and deposit their wind and SN material in a larger ISM volume. This leads to periods of time when a lower fraction of wind and SN material is being locked in new stars. The more efficient disruption of star clusters thus results in a locally lower recycling efficiency of massive star ejecta. This is also reflected in the mass-loading of these two components, which is overall larger in the ketju-run. As time progresses, the enriched material injected outside of star clusters gets mixed in the galactic ambient ISM and other star-forming regions. The wind and SN-material locking rates and loading factors of the ketju-simulation thus approach the sphgal-run over a time-scale of several 100 Myr. The subresolution method for metal-diffusion introduced in Aumer et al. (2013) is not enabled here, thus the mixing of metals is the result of turbulent mixing of the particles (see the discussion in Steinwandel et al. 2024a on the impact of metal-diffusion in Lagrangian codes).

Our sphgal + ketju simulation method can be used to model the galactic-scale impact of runaway massive stars that are ejected from star clusters due to their dynamical evolution. Previous high-resolution galaxy simulations have also explored the impact of runaway stars on the star formation and outflow properties of galaxies (Andersson et al. 2020; Andersson et al. 2023; Steinwandel et al. 2023), however using ad-hoc subresolution models to give the stars their ‘excess’ velocities. Notably, the models give a velocity kick to massive stars immediately when they form. This essentially means that the stars are almost instantly removed from the star-forming regions and instead injected into the less dense ISM when their pre-SN feedback channels have just become active. With a typical escape velocity of 10 km |$s^{-1}$|⁠, stars travel at least 10 pc from the cluster during the first Myr,⁷ therefore their critical role in regulating the gas-content of the clusters (see Section 4.5.3) may be fundamentally impacted by such ‘primordial’ kicks. In collisional N-body models, the escape of the stars instead occurs gradually as the stars segregate, form binaries and receive kicks in dynamical interactions. Compared to the previous studies that found increased outflow rates and both increased and decreased star formation histories associated with runaway stars, our simulations show almost identical mass-outflow rates and star formation histories with and without efficient star cluster disruption. The only notable change here is seen in the recycling of the metal-enriched material in newly formed stars, which is reduced due to the less clustered feedback. Our galaxies, however, also have lower metallicites and lower global SFRs compared to the previous studies, therefore a direct comparison is not possible.

4.3 SN clustering

In Fig. 6, we show the ISM densities around young stars at formation, around photoionizing young stars (⁠|$<1$| Myr) and around stars that have exploded as SN. To inspect how the accurate stellar dynamics impacts clustering of stellar feedback, we also show separately feedback events that occur in bound star clusters (bottom panel) recovered in both simulations (hyd_KETJU_m3_e01 and hyd_SPHGAL_m3_e01). While photoionization and star formation occur at slightly lower densities in the ketju-run, the SN environmental densities are quite similar in both runs. This manifestly shows the impact of effective early stellar feedback in regulating the gaseous densities in star-forming regions. The only significant difference in the two runs is the number of SNe in clusters, which is consistent with Fig. 5 where we saw a lower recycling efficiency of wind and SN-material. The number of SNe in bound clusters reduces from 993 to 509, out of 1428 and 1273 SNe in total, respectively. The fraction of SNe in clusters therefore drops from 70 per cent in the sphgal simulation to 40 per cent in the ketju-run. The effect of lowered SN-clustering on the galaxy outflows was however very minor as shown in Fig. 5, with a difference of only 10  per cent in the total gas mass that escapes the galactic disc between the two simulations.

Figure 6.

ISM density around stars younger than 1 Myr (dashed) and at the location of SN events (solid), compared to the star formation density (dotted). The top panel shows the values for all stars and the bottom panel shows the values for stars in bound star clusters. The photoionization and SN densities have been stacked in 1 Myr steps. The hyd_KETJU_m3_e01 and hyd_SPHGAL_m3_e01 simulations are shown in black and orange, respectively.

4.4 Evolution of the star cluster population

The cluster mass function (CMF) of all gravitationally bound star clusters is shown as a cumulative distribution in Fig. 7. We show both young clusters (0–10 Myr) and evolved clusters (100–110 Myr). The clusters are stacked over 400 Myr of evolution (between 0–400 Myr for the 0–10 Myr age bin and 100–500 Myr for the 100–110 Myr age bin) in 10 Myr steps so that the same clusters should appear once in both age bins. We limit the mass range to clusters less than |$\sim 900$|  |$\mathrm{M}_\odot$|⁠, which is the upper mass limit of the clusters formed in the ketju simulations. The simulations without ketju allow star clusters to grow on average to larger masses, thus we exclude three clusters more massive than |$\sim 900$|  |$\mathrm{M}_\odot$| in the hyd_SPHGAL_m3_e01-simulation from Fig. 7.

Figure 7.

The cumulative CMF of clusters in the age range of 0–10 Myr (solid) and 100–110 Myr (dashed) in the hyd_KETJU_m3_e01 (black) and hyd_SPHGAL_m3_e01 (orange) simulations. The CMFs are shown for clusters less massive than 900 |$\mathrm{M}_\odot$|⁠, limited to include the same mass range in both simulations. The lines on the top right show power-law functions |${\rm d}N/{\rm d}M\propto M^{\alpha }$| with slopes |$\alpha =-2$| and |$\alpha =-3$|⁠.

The CMF in the ketju-runs differ from the fully softened simulations both at young and at old cluster ages. First, the ketju-run has a shallower power-law shape at high cluster masses already at young cluster age. Observationally, the power-law slope of the CMF is often seen to be close to |$-2$|⁠, thus the ketju simulations seem to be in better agreement with the observations. However, the observed cluster samples are often more massive (⁠|$>10^{3-4}$|  |$\mathrm{M}_\odot$|⁠), have higher metallicity or are older, and are rarely complete at such low masses as considered here. We do not therefore attempt a direct comparison here. Other simulations of star cluster formation (e.g. Hislop et al. 2022; Garcia et al. 2023; Andersson et al. 2024) have previously found that the shape of the CMF and mass of the most massive cluster can also depend strongly on the stellar feedback and star formation prescriptions.

Second, the inefficient cluster disruption in the softened sphgal-run results in very little evolution in the CMF over 100 Myr, while the evolved CMF in the ketju-run has significantly lower number of clusters in all mass-bins, highlighting again the efficient mass-loss. The total number of clusters formed in both runs is approximately the same. In total 63 per cent of the clusters by number that formed in the ketju-run disrupt during the first 100 Myr of evolution, which is almost twice the corresponding value of 38 per cent of disrupted clusters in the sphgal-run.

4.4.1 Expansion of young compact clusters

In Fig. 8, we show the effective radii r|$_\mathrm{eff}$| (2D in x–y plane) of clusters in three age bins of 1–5, 5–10, and 10–100 Myr. The age ranges have been selected to enable direct comparison with observed star clusters in dwarf galaxies. We compare our cluster sizes to cluster data from the LEGUS survey (Brown & Gnedin 2021), here limited to low-mass galaxies with |$M_*< 10^9$|  |$\mathrm{M}_\odot$|⁠, and clusters observed in the LMC and SMC from Hunter et al. (2003) and Gatto et al. (2021). The radii used in Hunter et al. (2003) are from Bica et al. (1999). These three samples were selected as they cover the full age and mass range considered here and provide an estimate of the half-mass or half-light radius, instead of e.g. a core-radius or a 90 per cent light radius as is often the case in Magellanic cloud cluster data. Fig. 8 shows also the best-fitting relation from Brown & Gnedin (2021) and this relation scaled down by a factor of 8. The latter illustrates the initial cluster sizes expected from N-body simulations that start with compact configurations and end on the observed size–mass relation (Banerjee & Kroupa 2017; Arca Sedda et al. 2024).

Figure 8.

The effective radius (2D) of bound star clusters (circles) in three age bins: 1–5 Myr (left), 5–10 Myr (middle), and 10–100 Myr (right). The top row shows the hyd_KETJU_m3_e01-simulation and the bottom row shows the hyd_SPHGAL_m3_e01-simulation. Clusters that are still embedded, with at least 10 per cent of the mass within |$r_\mathrm{eff}$| in gas, have been indicated with open circles. Open grey squares are observed results from Brown & Gnedin (2021), Hunter et al. (2003), and Gatto et al. (2021, see the text for details). The solid line is the best-fitting relation from Brown & Gnedin (2021) and the dotted line is the same relation scaled down by a factor of 8.

The young clusters, especially more massive than 100 |$\mathrm{M}_\odot$|⁠, have initial sizes that are similar or larger than the scaled Brown & Gnedin (2021) relation in both simulations. The sphgal-clusters remain compact throughout their evolution. On the other hand, the increase in size of clusters simulated with ketju can be seen already at very young ages. The youngest massive (⁠|$>100$|  |$\mathrm{M}_\odot$|⁠) ketju-clusters have already started to show signs of expansion during the first 5 Myr. The 5–10 Myr clusters have already expanded by a factor of a few compared to their compact initial size, and the still older clusters that remain bound have effective radii up to 10 times larger than initially. The cluster sizes in the ketju-simulation match very well with the observed sizes in their respective age bins, while the sphgal-clusters run with softened gravity have too compact sizes at old age. In particular the more massive clusters are too compact in all age bins when no collisional dynamics is included.

4.5 Cluster evolution in a galactic context

4.5.1 Mass-loss through stellar evolution and dynamics

A more in-depth look into the evolution of individual clusters in the ketju and sphgal-simulations is provided in Fig. 9. We have selected long-lived (at least 150 Myr) clusters with maximum masses in the range of 400–900 |$\mathrm{M}_\odot$| from both simulations (8 clusters in hyd_KETJU_m3_e01, 17 clusters in hyd_SPHGAL_m3_e01), and followed the evolution of their properties from birth until the final snapshot at 500 Myr. We show in the left panel of Fig. 9 (from top to bottom) the evolution of the bound mass of the cluster, the number of bound stars, the effective radius (2D, similar to Fig. 8) and the median stellar mass. In the right panel of Fig. 9, these quantities are normalized to their initial value when the clusters reach maximum mass.

Figure 9.

From top to bottom: The bound mass, the number of bound stars, the effective radius (2D), and the median stellar mass in star clusters with maximum masses between 400 and 900 |$\mathrm{M}_\odot$| in the hyd_KETJU_m3_e01 (black) hyd_SPHGAL_m3_e01 (orange) simulations. The evolution is shown starting at the time of maximum cluster mass |$t_\mathrm{max({\it M})}$|⁠. The left column shows the absolute values and the right column shows the values scaled by the equivalent value at |$t_\mathrm{max({\it M})}$|⁠. The vertical lines indicate epochs of 10 and 100 Myr.

The dynamical loss of bound stars in the ketju-run is accompanied by rapid size–evolution already at very early cluster ages. This rapidly emerging effect is not seen in the sphgal-run where the loss of stars is more gradual. Clusters in both simulations lose mass, especially through stellar evolution, but the softened simulation underestimates the dynamical mass-loss. It has been shown in N-body simulations that less extended clusters that have underfilled Roche lobes disrupt slower (Tanikawa & Fukushige 2005; Shukirgaliyev et al. 2018). Thus the secondary effect of suppressed two-body interactions and the compact softening-supported size of the sphgal-clusters is to lose mass at an even slower rate.

The total lifetimes of star clusters can be estimated by fitting a mass-loss rate parametrized in Baumgardt & Makino (2003) and Lamers et al. (2005) as

$$\begin{eqnarray} \frac{{\rm d}M_\mathrm{cluster}}{{\rm d}t}=-\frac{\mathrm{M_\odot }}{t_0} \left(\frac{M_\mathrm{cluster}}{\mathrm{M_\odot }}\right)^{1-\gamma }, \end{eqnarray}$$

(6)

with time-scale parameter |$t_0$| and |$\gamma =0.7$|⁠. Here, we set the initial mass of the fit |$M_i$| at |$t-t_\mathrm{max({\it M})}=100$| Myr, when the mass-loss has reached a gradual stage in the top panels of Fig. 9. The best-fitting total dissolution time |$t_\mathrm{dis}=t_0(M_i/M_\odot)^{\gamma }$| is in the range of 0.8–2.3 Gyr for these selected clusters in the ketju simulation. The equivalent values for the sphgal-run are between 2.0 and 4.2 Gyr. The total lifetimes of the clusters evolved with softened gravity in the tidal field of a dwarf galaxy are overestimated by at least factor of two at initial cluster mass of 400–900 |$\mathrm{M}_\odot$|⁠. The discrepancy between the runs with and without ketju increases with decreasing cluster mass, as the initial cluster mass is less clearly correlated with the total cluster lifetime. On the other hand, the smallest mass clusters in the ketju-simulation have always shorter estimated lifetimes compared to the more massive clusters.

The mass-segregation and escape of lower mass stars is quantified in the median stellar mass in the bottom panel of Fig. 9. The median stellar mass increases, in both simulations, which can only be the result of low-mass stars being removed from the stellar population. The median stellar masses of most of the ketju-clusters evolve faster compared to the softened simulation. There is, however, some preferential removal of low-mass stars in the sphgal-simulation as well. This indicates that the mass-loss process is partially but not entirely captured even when softened gravitational forces are used.

4.5.2 Evolution in the mass–size plane

In Fig. 10, we show the mass and size evolution of the individual clusters in the hyd_KETJU_m3_e01 and hyd_SPHGAL_m3_e01-simulations from Fig. 9, compared to the same observed mass-size data from Fig. 8. The clusters are shown through 100 Myr of evolution starting from the first snapshot in which they were identified in, as opposed to the snapshot of maximum mass as in Fig. 9. Here, we include the early embedded phase, defined in this analysis as the period of time when the gas mass fraction within 1 pc of the cluster centre of mass is larger than 10  per cent. We highlight these epochs to indicate when the clusters are embedded and when they later become fully exposed. Observed comparison data is the same as in Fig. 8, here stacked across all clusters with ages less than 100 Myr located in low-mass galaxies (⁠|$M_*< 10^9$|  |$\mathrm{M}_\odot$|⁠) from Brown & Gnedin (2021), Hunter et al. (2003) and Gatto et al. (2021).

Figure 10.

The evolution of the effective radius (2D, measured in the x–y plane) and the bound stellar mass of star clusters with maximum masses between 400 and 900 |$\mathrm{M}_\odot$| in the hyd_KETJU_m3_e01 (left) and hyd_SPHGAL_m3_e01 (right) simulations. The data points along each evolution track have been coloured according to the cluster age, and the black-outlined triangles indicate epochs when the cluster is embedded in at least 10 per cent of its current mass in gas (within |$r_\mathrm{eff}$|⁠). The observed data show young and intermediate age clusters (⁠|$< 100$| Myr) in dwarf galaxies from the LEGUS dataset (blue hexbins, Brown & Gnedin 2021) and the Magellanic clouds (green crosses from Gatto et al. 2021 and pluses from Hunter et al. 2003).

As was already shown in Fig. 9, the ketju-clusters evolve more rapidly towards the observed range of cluster sizes. The sizes first contract as the clusters grow in mass, but the radii begin increasing already when the clusters are still embedded. The sphgal-clusters, on the other hand, reach maximum mass at high densities and remain in a narrow region of parameter space because neither the size nor the total mass evolve much as shown in Fig. 9.

The young clusters in both simulations are located in a region of parameter space that is not well populated by observed clusters. The observed comparison data shown in Figs 8 and 10 extend down to young cluster ages but lacks the very youngest, embedded phase. Observations of the total sizes of embedded clusters in the LMC (Romita, Lada & Cioni 2016) span a range from |$\sim 0.3$| pc upwards, with typical values of 1 pc. Low-mass embedded clusters in the Milky Way also exhibit small, subparsec sizes (Allen et al. 2007). The effective sizes of observed embedded clusters are thus expected to be smaller than the exposed clusters included in the comparison sample. Our measured sizes in the embedded phase should thus be in agreement with the sizes of real embedded low-mass star clusters.

The size–measurement technique we use is also not fully consistent with the observational method of measuring the cluster sizes. The sizes of the observed comparison clusters have been obtained through fitting a light-profile to the cluster emission. This makes the sizes sensitive to background estimation and difficult to measure in crowded environments. We instead only consider gravitationally bound stars when computing the sizes and masses of the clusters. The ketju-clusters, in particular, have extended envelopes of unbound stars in the tidal tails that are excluded by the boundness criterion in the definition of our star clusters. A less strict clustering algorithm might thus result in even larger measured sizes at evolved cluster ages. The sizes measured through synthetic photometry might also differ from the measurements of direct particle data.

4.5.3 Gas expulsion

Exposed clusters in local galaxies are coincident with young stellar ages (Whitmore et al. 2014; Hollyhead et al. 2016; Grasha et al. 2018; Hannon et al. 2022), indicating that the embedded phase lasts only a few Myr. Fast gas expulsion is further supported by high-resolution simulations e.g. by Dinnbier & Walch (2020) and Farias et al. (2024). Based on the number of embedded data points per track in Fig. 10, we see the bound clusters remain embedded for up to 5 Myr before full gas removal.

We further inspect the intersection of embedded and exposed clusters in the hyd_KETJU_m3_e01  ketju-run in Fig. 11 by showing the time of the birth of the first photoionizing star and the first SN event within 20 pc of the cluster centre of mass superimposed on the mass–size evolution sequence from Fig. 10. The embedded phase comes to an end when the average stellar age is a couple Myr⁸ and often coincident with the ignition of the first photoionizing star. The first SNe occur only Myrs later, and do not contribute to the gas removal of the host cluster. We note that in the simulations of the present study, the shortest lifetimes of the SN-progenitor stars with masses up to 40 |$\mathrm{M}_\odot$| are |$\sim 5$| Myr in the absence of stars in the pair-instability mass range. In a more intensely star-forming galaxy (e.g. Lahén et al. 2024), the onset of SNe could be hastened through the formation of more massive star clusters that are able to host more massive stars with lifetimes as short as a couple Myrs.

Figure 11.

The mass–radius evolution of clusters in the hyd_KETJU_m3_e01 simulation as shown in Fig. 10. The time of the first phototionizing star (red circles) and the first SN (black cross) within 20 pc of the cluster centre of mass. The clusters become gas free right after the first massive stars appear, while SNe occur only Myrs later.

The small sizes and the gas expulsion induced by photoionization in our simulated young clusters therefore seem to be in agreement with observed star clusters, and a picture emerges wherein the clusters evolve rapidly towards lower central densities during and after gas removal. The driver of gas removal is photoionizing radiation, well before SNe in the vicinity of the clusters come into play. This result is in agreement with previous simulations that have shown that a sufficiently luminous source can significantly disrupt the gas cloud before the SNe occur (Dale, Ercolano & Bonnell 2013; Geen et al. 2016; Kim, Kim & Ostriker 2018; He et al. 2019; Fukushima et al. 2020; He et al. 2020; Guszejnov et al. 2022), and our results extend these previous works with a model that includes individually resolved radiation sources in a full galactic context. The present results are also in qualitative agreement with previous analysis of the simulations in the griffin-project on starburst systems (Fotopoulou et al. 2024) where we have shown that even the most massive cold star-forming ISM-clouds (⁠|$>10^5$|  |$\mathrm{M}_\odot$|⁠) can be destroyed already after the first SN. The pre-processing by early stellar feedback is thus a crucial driver of evolution of clustered star-forming regions across environments, though full gas removal by photoionization might be only possible in relatively low-mass clusters. The drastic difference in the size evolution of the collisional and softened clusters in Fig. 10, however, highlights the fact that the full impact of gas dispersal and removal on the evolution of star clusters can only be addressed in hydrodynamical simulations that account for their collisional nature. The internal evolution of more massive clusters with a prolonged embedded phase will be addressed in a future study.

5 SUMMARY AND CONCLUSIONS

We present new hydrodynamical + collisional N-body simulations of star formation in low-metallicity dwarf galaxies, concentrating on the formation and evolution of gravitationally bound star clusters. The simulations are run with the modern SPH code sphgal now also including regularized small-scale gravitational dynamics using the ketju integrator. We have updated the interface between sphgal and ketju which enables us to place regularised integration regions around all selected stars above a mass threshold and their remnants, as outlined in Table 1. This enables collisional dynamics within star clusters in the vicinity of massive stars to be resolved accurately without the need for gravitational softening.

We first benchmark the updated sphgal + ketju method with the direct N-body code bifrost by running idealized star clusters in isolation without stellar evolution. sphgal + ketju captures the mass-segregation and two/few-body interactions within the central regions of the clusters, thus producing similar numbers of unbound stars (to within a few tens of per cent) as the direct N-body code. These stars escape the cluster as runaway or walkaway stars even without an external tidal field. Because the majority of interactions between lower mass stars in the outskirts of clusters are still softened, the evolution towards energy equipartition is slowed down. This results in a somewhat slower expansion of the clusters in the sphgal + ketju-simulations compared to bifrost. When the number of ketju-regions is restricted to only stars more massive than 8 |$\mathrm{M}_\odot$|⁠, only a few per cent of the gravitational interactions are regularized. Even then, the clusters still expand significantly and produce runaway stars, owing to mass-segregation and strong few-body interactions in the central regions of the clusters.

We then analyse the long-term evolution of low-metallicity dwarf galaxies integrated with the sphgal + ketju code, in comparison to simulations that adopt gravitational softening everywhere. The global properties, such as the total new stellar mass, the outflow rate and mass-loading factors are quite similar, with differences in the cumulative properties only of the order of a few tens of per cent. SNe are less clustered in the ketju-run with regularization due to cluster dissolution, which leads to less efficient recycling of wind and SN-material in new stars. Most of the wind and SN-material still escapes the galaxy in metal-enriched outflows. Photoionization regulates the gas densities, thus the lowered SN-clustering has a little impact on the global star-forming environment.

The major impact of adding accurate gravitational interactions is seen in the formation and evolution of star clusters in the galactic tidal field. The average values of |$\Gamma$| for clusters in the age interval of 1–10 Myr are lowered from |$\sim 40~{{\ \rm per\ cent}}$|–50  per cent in the softened simulation to |$\sim 30~{{\ \rm per\ cent}}$|–35  per cent in the ketju-run due to the rapid loss of stars from the clusters. Evolved clusters in the range 10–100 Myr contain only |$\sim 13~{{\ \rm per\ cent}}$| of all stellar mass in the system in the ketju-run. This is in agreement with observations where |$\Gamma$| at older cluster age can be used to probe cluster dissolution when cluster mass-loss cannot be directly followed. Approximately, 63 per cent of the clusters by number in the ketju-simulation dissolve completely during the first |$\sim 100$| Myr of evolution, compared to 38 per cent in the run where softening suppresses cluster disruption.

We also took a more detailed look at the evolution of a limited sample of individual long-lived clusters formed in the hydrodynamical simulations, with initial masses of a few hundred solar masses. The clusters in both simulations suffer mass-loss through stellar winds and SNe. The ketju-runs have a large fraction (up to 40  per cent by number) of initially bound stars that become unbound during the first 10–100 Myr. The |$r_\mathrm{eff}$| of selected clusters in the initial mass range of 400–900 |$\mathrm{M}_\odot$| are shown to increase by a factor of 2–6 during the first 10 Myr of evolution, and by a factor of 5–9 during the first 100 Myr. Also the fully softened clusters expand somewhat, by a factor of 2–3, and they lose bound stars (up to |$\sim 10~{{\ \rm per\ cent}}$|⁠) due to the interactions with the galactic tidal field. This is in contrast with the pure N-body simulations of isolated star clusters where the structure of the fully softened clusters did not evolve nearly at all. The mass-loss in the softened clusters is, however, too slow by at least by a factor of two.

As the ketju-clusters evolve in the mass–size plane, they end up in the observed range of parameter space that is occupied by exposed relatively young star clusters in low-mass galaxies. We showed that the sizes first contract while clusters are still forming stars, but the sizes can start increasing already when the clusters are still embedded. The sizes then increase rapidly once the first photoionizing stars are ignited and the clusters become fully exposed. The gas expulsion is governed by the pre-SN stellar feedback while SNe occur only Myrs after the gas has already been removed.

In this study we discuss clusters less massive than 1000 |$\mathrm{M}_\odot$|⁠. The star formation time-scale is longer (He et al. 2019; Lahén et al. 2020) and the gas retention may be more efficient in more massive clusters as indicated in our previous work on star-by-star starburst environments in Lahén et al. (2024). In more massive clusters, where the gas can be retained in the cluster centre for a longer period of time, stellar winds and SN feedback may become more important for chemical enrichment and gas removal. Investigation of more massive clusters will be left for a future study, where we will apply the new hydrodynamical + N-body method in a low-metallicity starburst. An implementation for stellar interactions and mergers would enable us to study the collisional growth of very massive stars in dense star clusters (Portegies Zwart et al. 2004). Such dense cluster-forming regions may also have given rise to globular clusters in the early Universe (Adamo et al. 2024), which may have been important sources of ionising radiation (He et al. 2020). The rapid (Myr-scale) collisional growth of massive stars in these dense clusters may provide a channel for enhanced chemical enrichment (Gieles et al. 2018) and black hole formation in the mass range of |$\sim 10^3$|  |$\mathrm{M}_\odot$|–|$10^4$|  |$\mathrm{M}_\odot$| (Fujii et al. 2024; Rantala et al. 2024). With this simulation methodology we may begin to address the origin of the chemical enhancements (Bunker et al. 2023) and the seeds of the supermassive black holes (Bogdán et al. 2024) recently uncovered with JWST at unexpectedly high redshifts.

ACKNOWLEDGEMENTS

NL thanks Deidre Hunter and Bruce Elmegreen for providing the data reported in Hunter et al. (2003). 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 (LRZ, www.lrz.de) under project numbers pn49qi and pn72bu. 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’. PHJ acknowledges the support by the European Research Council via ERC Consolidator Grant KETJU (no. 818930) and the support of the Academy of Finland grant 339127. The computations were carried out at SuperMUC-NG hosted by the LRZ and the COBRA and FREYA clusters hosted by The Max Planck Computing and Data Facility (MPCDF) in Garching, Germany.

This research made use of python packages scipy (Virtanen et al. 2020), numpy (Harris et al. 2020), matplotlib (Hunter 2007), pygad (Röttgers et al. 2020), h5py (Collette 2013) and Astropy⁹ (Astropy Collaboration et al. 2022).

DATA AVAILABILITY

The data will be made available on reasonable request to the corresponding author.

Footnotes

1

https://www.mpa-garching.mpg.de/~naab/griffin-project

2

Stars below the lower mass limit of 1.7 |$\mathrm{M}_\odot$| in the |$Z=0.0004$||$\sim 0.02$|  |$\mathrm{Z}_\odot$| Geneva tables are supplemented with fluxes from the |$Z=0.002$||$\sim 0.1$|  |$\mathrm{Z}_\odot$| tables of Georgy et al. (2013) between 0.8 and 1.7 |$\mathrm{M}_\odot$| by scaling up the mass-dependent fluxes by a factor of 2, which results in approximately continuous mass–flux relation across the mass range. Below 0.8 |$\mathrm{M}_\odot$| the fluxes are extrapolated as |$L\propto M^{3.5}$|⁠.

3

https://www.mv.helsinki.fi/home/phjohans/ketju

4

The latter is named as such because there are no SPH particles or stellar evolution in the pure N-body runs, and the computation of gravitational dynamics is done in sphgal with the standard leap-frog integrator and tree force algorithm of gadget-3.

5

For a detailed description of the scaling properties of regularized integrators, see section 5.2 of Rantala et al. (2017). But please note that the discussion considers an older version of ketju, before updates to its perturber treatment (Mannerkoski et al. 2022) and parallelization (Rantala et al. 2020).

6

See section 5.4.1 in Lahén et al. (2023) for more details and further discussion regarding our definition of the mass-loading.

7

The minimum velocity kick imposed in the subresolution models is often 3 km s|$^{-1}$|⁠.

8

As opposed to the total stellar age spread or the total lifetime of the star cluster, which can be larger.

9

http://www.astropy.org

APPENDIX A: CENTRAL 5 PER CENT LAGRANGIAN RADII OF INDIVIDUAL CLUSTERS

Fig. A1 shows the inner 5 per cent Lagrangian radii in each individual simulation of the code comparison sample discussed in Section 3.1. The mean standard deviation of the simulations shown in each panel of Fig. A1 are shown in Fig. 1. The cycles of contraction and rapid expansion are most clearly seen in the BIFROST runs. The KETJU simulations do also show these cycles, albeit often with a smaller amplitude. The inner regions of the GADGET simulations show practically no evolution.

Figure A1.

The evolution of the 5 per cent Lagrangian radii of the 10 individual randomly generated clusters run in Section 3. The lines from top to bottom in each panel show the simulations started from the same initial condition. The lines have been shifted to ease their differentiation, therefore the y-axis is arbitrary. The simulation samples from left to right, top to bottom are nb_10k_BIFROST, nb_10k_KETJU_m3_e001, nb_10k_KETJU_m8_e001, and nb_10k_GADGET_e001.