© The Authors 2024

1. Introduction

Globular clusters (GCs), which are found in every massive galaxy, are some of the most ancient objects in the Universe. They serve as fossils that can reveal the elusive environment and physics of the early phases of galaxy assembly (Brodie & Strader 2006; Portegies Zwart et al. 2010; Renaud et al. 2017; Krumholz et al. 2019; Adamo et al. 2020). Yet because of their age, many aspects of cluster formation and evolution at high redshift are challenging to constrain through observation, and little is known about the efficiency and timescale at which gas is converted into stars to create such massive bound clusters.

Though the progenitors of GCs are too old to observe, there are younger star clusters that are as massive as GCs and currently forming in nearby galaxies. These young massive clusters (YMCs) have masses M ≥ 10⁴ M_⊙ and ages < 100 Myr (Portegies Zwart et al. 2010). The discovery of these objects has indicated that the mode of extreme star formation that forms massive star clusters still occurs today. Notably, even more of these clusters are being discovered with JWST, as many YMCs in starburst galaxies are too embedded to have been seen by Hubble (Whitmore et al. 2023). Although it has been suggested that YMCs are the present day analogs to young GCs, this is debated in the literature (see Renaud 2020).

Theory suggests that, despite the abundance of GCs, ≤1% of clusters survive to become GCs (Fall & Zhang 2001; Fall et al. 2005; Fall 2006). The conditions that lead to bound star clusters as massive as GCs remain a mystery, and observations of forming YMCs are sparse. Star formation must be fast and efficient enough to form bound stars that can survive the epoch of stellar feedback and the dispersal of the natal gas (Lada & Lada 2003). The plethora of GCs suggests these conditions were very common in the early Universe.

The process of star formation in a cluster begins with the global gravo-turbulent collapse of giant molecular clouds (GMCs; Larson 1981). As the collapse proceeds, fragmentation creates overdense clumps within the GMC that begin to form stars (Mac Low & Klessen 2004; McKee & Ostriker 2007; Klessen & Glover 2016). The feedback from these stars, in the form of stellar winds, jets, and radiation, begins to clear out dense gas in and around the forming sub-clusters, slowing down the local (sub-cluster scale) and global (cloud scale) star formation rate (SFR; e.g., Girichidis et al. 2020; Lewis et al. 2023). Eventually, massive stars explode as supernovae (SNe), further dispersing gas. However, it has been argued that the efficiency at which stellar feedback slows global star formation diminishes with higher gas surface density (Grudić et al. 2018). The sub-clusters eventually merge if they are mutually gravitationally bound, forming a final cluster cleared of all natal gas (Krause et al. 2020).

Many details of star cluster formation remain poorly understood due to the difficulty of modelling such a complex process. Stellar evolution and binary dynamics need to be resolved on timescales of years and distance scales of an AU, while the magnetohydrodynamics (MHD) of the collapsing gas covers regions several parsecs across, with crossing times of thousands to millions of years. Because of this, most computational star cluster formation models are limited and must make considerable approximations. Many simulations do not form individual stars: some apply stellar feedback as a combined source in the center of the cloud (Dale et al. 2005; Rahner et al. 2019), and others use sink particles representing sub-clusters (e.g., Bate et al. 1995; Federrath et al. 2010) or extract the properties and feedback of individual stars from the sink particles (e.g., Sormani et al. 2017; Howard et al. 2017; Kim et al. 2017; Grudić et al. 2018; Su et al. 2018). Other simulations do form single stars, but they do not resolve the stellar feedback of each individual star particle (Colín et al. 2013a; Li et al. 2019), instead including feedback from just the sink particles that created the stars. Simulations of dwarf galaxies can capture star cluster mass functions and formation times, but they do so without collisional dynamics of star particles and are therefore unable to accurately capture dynamical properties such as velocity dispersion and size (Lahén et al. 2019, 2024; Andersson et al. 2024).

Modelling individual stars is important, as this can change the efficiency and location of stellar feedback injection. Dynamical processes often eject high-mass stars (Fujii & Portegies Zwart 2011; Fujii et al. 2022a), and the location of massive stars directly affects how and when gas is dispersed. Gas dispersal stops star formation. Models of sub-cluster feedback may overestimate the strength of feedback, as they do not allow for spatial separation between the stars in the sub-cluster. This lack of separation also changes the morphology of the gas, affecting the number of low-density channels in the gas that can vent thermal energy from the sub-cluster. The degree to which the sub-cluster and star-by-star approaches differ must be constrained.

There are a few models that do evolve individual stars with both stellar feedback and higher order gravitational dynamics (Wall et al. 2020; Grudić et al. 2021; Fujii et al. 2021, 2022b,a; Cournoyer-Cloutier et al. 2021, 2023; Lewis et al. 2023; Wilhelm et al. 2023). While these models include most of the relevant physics, they lack the computational efficiency to simulate star clusters forming from clouds of masses > 10⁵ M_⊙, and instead the models focus on simulating star clusters forming from low-mass clouds ≤10⁵ M_⊙. This leaves a sizeable gap compared to the observed mass range of GMCs. While clusters with mass < 10⁵ M_⊙ are comparable to Local Group observations, YMC and GC formation is out of their reach. Furthermore, most star formation takes place in GMCs of mass ≥10⁵ M_⊙ (McKee & Williams 1997; Murray & Rahman 2009).

The goal of this work is to model the formation of massive clusters from their initial GMCs while following the formation of individual stars and their feedback. We aim to answer how and in what conditions YMCs can form while remaining bound throughout the onset of gas expulsion. We also seek to understand how efficient the process of star formation is in a cluster, what the timescale is over which star formation occurs, and whether the clusters formed from these massive clouds survive and remain bound or quickly disperse. We plan to compare our results to those that use a sub-cluster formation and feedback model.

To do this, we used the TORCH framework (Wall et al. 2019, 2020). TORCH employs the Astrophysical Multipurpose Software Environment (AMUSE) framework to couple separate physics codes that handle MHD, radiative transfer, stellar evolution, and N-body dynamics. TORCH uses the MHD code FLASH (Fryxell et al. 2000; Dubey et al. 2014), which accounts for the evolution of the gas and the formation of sink particles and stars. Stellar feedback in the form of winds and SNe is included, and the effect of ionizing and non-ionizing radiation is followed using a ray-tracing approach (Baczynski et al. 2015). The star formation model samples the Kroupa (2002) initial mass function (IMF) to form stars from sink mass reservoirs (Wall et al. 2019). SEBA (Portegies Zwart & Verbunt 1996) tracks stellar evolution from the zero-age main sequence, and, in the original version of TORCH, PH4 (McMillan et al. 2012) + MULTIPLES (Portegies Zwart & McMillan 2018) handled the stellar dynamics.

In that version (Wall et al. 2019, 2020), TORCH could not handle the hundreds of thousands of stars that form in massive GMCs > 10⁵ M_⊙. In this work, we solve this problem by making three updates: 1) We replace the combination of the N-body code PH4 and the higher-order interactions MULTIPLES with the code PETAR (Wang et al. 2020a); (2) we agglomerate stars with masses < 4 M_⊙ into summed-mass dynamic star particles with masses of ≥ 4 M_⊙; and 3) we mass load the stellar winds to reduce the peak temperatures beyond their termination shocks. These modifications enabled TORCH to then model clouds with an initial mass of up to 10⁶ M_⊙ that form hundreds of thousands of individual stars.

We present simulations of star clusters forming from turbulent spherical clouds with masses of 10⁴,  10⁵, and 10⁶ M_⊙. Each of these clouds is almost identical in terms of their initial properties, with only mass and density scaled between them. Our study investigates whether the formation of YMCs parallels that of low-mass clusters or if it varies significantly with initial cloud mass and density.

This paper is the first in a series exploring the results of these simulations. In this paper, we describe the TORCH code, the new features integrated into TORCH for handling massive GMCs, and the initial conditions of our three clouds in Sect. 2. We analyze the time evolution of global gas and stellar properties in Sect. 3. In Sect. 4, we discuss the results of our analysis, and in Sect. 5 we conclude with a summary of the most important results. We provide a data repository containing a sampling of the simulation data corresponding to the panels in Figure 1, the HTML file of the three-dimensional interactive plot shown in Figure 2, and the code used to generate the interactive plot.

Fig. 1. Slice plots of the three simulations in the x–y plane over time. The plane of the slices for a given cloud is the center of stellar mass in the final snapshot. Stellar positions are shown by white dots. The free-fall times tff are given in Table 1. The number of stars shows the amount of star particles in the domain, not the number sampled from the IMF. Due to our agglomeration of low-mass star particles, the number of stars sampled from the IMF is ∼10× greater.

Fig. 2. Still of the interactive plot of the embedded M4, M5, and M6 clusters (left to right) at 1 tff. The interactive plot file is available for download from the repository.

2. Methods

2.1. Standard TORCH

TORCH¹ is built upon the AMUSE framework, which couples multiple autonomous astrophysical codes. We chose codes that allowed efficient calculation of the disparate physical processes at work in star cluster formation.

The TORCH framework incorporates the adaptive mesh refinement MHD code FLASH v4.6.2 (Fryxell et al. 2000; Dubey et al. 2014) with a number of enhancements to the base version of FLASH. The base FLASH handles the MHD and sink particle formation and evolution. The modifications to FLASH presented in Wall et al. (2019, 2020) include heating and cooling, ionization, radiation transfer (using ray-tracing; see Baczynski et al. 2015), and feedback injection from stars. Stellar feedback is implemented in FLASH in the form of ionizing extreme ultraviolet (EUV) and non-ionizing far ultraviolet (FUV) radiation in the form of radiative heating and radiation pressure, as well as mechanical feedback from stellar winds and SNe. FUV rays are terminated when their flux drops below F_ray ≤ 16.9 G₀e^{−3.5 A_v}, where A_v is the visual extinction and G₀ is the Habing flux. This cutoff is 10× the applied background FUV field of F_ext = 1.69 G₀e^{−3.5 A_v} (Draine 1978). This limits the number of low energy rays on the grid for computational efficiency. We used the HLLD Riemann solver (Miyoshi & Kusano 2005) in FLASH paired with third-order piecewise parabolic method reconstruction (Colella & Woodward 1984).

To avoid artificial fragmentation, the Jeans length,

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{J}}=\sqrt{\pi c_{\mathrm{s}}^2/(G\rho )},\end{array}$$(1)

must be resolved by at least four cells (Truelove et al. 1997). We used a refinement criterion of 12 cells per Jeans length along with a derefinement criterion of 24 cells per Jeans length. As density increases during collapse, the Jeans length decreases until this criterion is no longer met at the highest level of AMR refinement. Sink particles were used to collect the gas that exceeds this density. The Truelove criterion sets the sink radius to R_sink = 2.5Δx_min and gives the sink threshold density during the entire run as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{sink}}=\frac{\pi c_{\mathrm{s}}^2}{G{\lambda }_{\mathrm{J}}^2}=\frac{\pi c_{\mathrm{s}}^2}{G{(5\mathrm{\Delta }{x}_{min})}^2},\end{array}$$(2)

where c_s was evaluated using the initial temperature of the gas. (If the gas heats during the run, the dense gas will be better resolved, making this a worst-case limit for the required density resolution.)

On each time step, the mass of gas within a distance R_sink of a sink particle that satisfies the criteria outlined in Federrath et al. (2010) is added to that sink’s mass reservoir for creating stars. When a sink forms, it randomly samples the Kroupa IMF (Kroupa 2002) and stores a long list of potential star masses to form (see also Sormani et al. 2017). Each time step, the sink forms as many stars from this mass list as possible until its current mass reservoir is depleted. It again forms one or more stars the next time it has accreted enough mass for at least the next star on the list. This is the standard stellar mass sampling method used in TORCH (Wall et al. 2019). Star positions are randomly sampled from a uniform spherical distribution within the sink’s accretion radius. Star velocities are set by the sink velocity added to an additional isotropic velocity dispersion with a Gaussian distribution having a standard deviation of the local sound speed.

Star particles are initially formed as zero-age main sequence stars, neglecting pre-main sequence evolution. Subsequent stellar evolution is tracked with SEBA (Portegies Zwart & Verbunt 1996), which passes the evolutionary properties informing stellar feedback to FLASH. The N-body dynamics of the stars are calculated using PETAR (Wang et al. 2020a), which is discussed further in the next section. Stars dynamically interact with the gas in FLASH through the AMUSE hierarchical coupling (Portegies Zwart et al. 2009) based on the gravity-bridge algorithm of Fujii et al. (2007).

2.2. PETAR N-body

TORCH was first designed to use the N-body code PH4 (McMillan et al. 2012) to handle direct stellar dynamics, paired with MULTIPLES (Portegies Zwart & McMillan 2018) to track binary orbital evolution and higher order perturbations. For TORCH runs using an initial cloud of 10⁴ M_⊙ and producing only a few thousand stars, this works well. However, the computational cost becomes unfeasible when pushing to higher initial cloud masses that produce far more than 10⁴ stars with over a few hundred binary systems. This is because MULTIPLES is a serial Python code, so with many interacting stars computation times become impractical. To speed up TORCH, we replaced PH4 and MULTIPLES with PETAR (Wang et al. 2020a).

PETAR is a state-of-the-art gravitational dynamics code optimized for solving the stellar dynamics of systems with millions of stars. It accomplishes this by dividing gravitational interactions into three regimes: distant interactions calculated with a Barnes & Hut (1986) tree and handled by the framework for developing parallel particle simulation codes (FDPS; Iwasawa et al. 2016, 2020), nearby interactions solved with a fourth-order Hermite direct N-body integrator (Makino & Aarseth 1992), and close interactions (binaries and higher order systems/perturbations) solved using the Slow Down Algorithmic Regularization (SDAR) technique (Wang et al. 2020b). For each particle, the force from neighboring particles is solved depending on what distance regime they are in, with a mass-dependent factor to increase the distance over which massive particles are considered close neighbors. The SDAR feature for handling higher-order dynamics is the novel component of PETAR, enabling it to handle large numbers of binaries and higher order systems in parallel.

In Figure 3, we plot the wall-clock time per evolution step for each of the N-body codes considered. For reproducibility, this test was done with the parameters $r_{\mathrm{out}}=0.001\mathrm{pc}$, $r_{\mathrm{bin}}=100\mathrm{AU}$ for PETAR and the stellar interaction radius $r_{\mathrm{int}}=15{\mathrm{R}}_{\odot }$ for MULTIPLES. PETAR is significantly faster and consistently performs well as the number of stars increases. The variability in the performance of PH4 and MULTIPLES is due to MULTIPLES taking longer if there are many third-body perturbations in a given step. We note that this test was done with single stars only; the scaling for a run with primordial binaries will be different.

Fig. 3. Wall-clock time for an evolution step for PETAR and PH4+MULTIPLES given the number of stars.

When running PETAR in TORCH, the time step of the long distance particle tree must be set (dtsoft), as well as the changeover radius between direct N-body and tree method for force calculations (rout). If the user sets these two parameters, all other parameters are set automatically. In TORCH, the MHD code sets the global time step for all worker codes based on the Courant condition. The tree time step was set as the nearest power of two in code units below the requested time step, as a power of two is required by PETAR (like most N-body codes). This sets dtsoft. For M5 and M6 we set the outer radius rout to 0.001 pc, the standard value used in PETAR simulations. We used a softening length of $\ell =15{\mathrm{R}}_{\odot }$ and a binary search radius of $r_{\mathrm{bin}}=16.5\mathrm{AU}$. For M4 we used a larger $r_{\mathrm{out}}=10r_{\mathrm{sink}}=7.8\mathrm{pc}$ to ensure accurate force calculations given the small number of stars and low stellar density. This corresponds to $r_{\mathrm{bin}}=0.63\mathrm{pc}$.

The code handling stellar mergers within PETAR is not active within the AMUSE framework, which results in star particles that approach within each others’ softening radius and should merge instead ending up with identical positions, leading to a halt in code execution. We have implemented code to check for particles in this state, and merge them. We intended to use SEBA to update the stellar mass of the surviving star, but later testing revealed that the surviving star’s mass remained unchanged. One star in the merger is removed meaning that stellar mass is unphysically lost. The effect of this error is negligible due to low merger rates: there are 0, 2, and 4 mergers in M4, M5, and M6, respectively. All of these mergers involve stars < 7 M_⊙. M5 and M6 lose only 8.4 M_⊙ and 22 M_⊙ of stellar mass due to unphysical mergers over the course of the simulations.

2.3. Stellar modifications

We made three alterations to the star formation and evolution procedures that vary from standard TORCH to accommodate the several orders of magnitude increase in number of stars formed when increasing the initial cloud mass from 10⁴ M_⊙ to 10⁶ M_⊙.

1. m We agglomerated low-mass star particles below M_agg = 4 M_⊙ as they formed until their summed mass is ≥M_agg. Then, a star particle is formed with a mass equal to the sum of the low-mass stars. This reduces the strain on the N-body calculations by reducing the number of dynamical star particles by 90%.

2. We mass-loaded stellar winds to raise the Courant time step by limiting the temperature of wind-blown bubbles to $T_w=3\times {10}^5\mathrm{K}$, which significantly sped up the simulations. This resulted in smaller, cooler, momentum-conserving bubbles instead of hot energy-conserving bubbles. The primary effect of wind feedback in cluster formation is to clear out extremely dense gas in order to allow ionizing radiation to form expanding H II regions. In this dense gas even hot stellar wind bubbles cool quickly, so there is little change in behavior in this regime.

3. We only injected feedback from stars above 20 M_⊙ to reduce the cost of ray-tracing. Massive stars output most of the ionizing radiation and mechanical wind energy in clusters: by setting this limit we lost less than 20% of the total feedback energy. Stars below the feedback cutoff mass did not go SN within the time frame of our simulations (≤10 Myr).

We further explain and examine the effects of these modifications in Appendix B, including providing a quantitative analysis of the amount of total energy lost by excluding feedback for stars < 20 M_⊙ in the M6 model.

2.4. Initial conditions

The initial properties of our three clouds are listed in Tables 1 and 2. We chose to keep the radius of all three clouds constant at R_cl = 11.7  pc. The radius was kept the same to have the same spatial distribution of star formation for each run. Constant radius allows the cell resolution and size of sink particles to be the same between the three simulations, and it facilitates directly comparing the morphology and dynamics of the forming clusters.

Consequently, the average initial densities of the clouds are 1.5, 15, and 150  M_⊙ pc⁻³, or 10⁻²²,  10⁻²¹, and 10⁻²⁰ g cm⁻³. The column densities of these clouds are 2.325 × 10^1, 2, 3 M_⊙ pc⁻², respectively. Assuming a 9:1 number ratio of H:He, resulting in a mean molecular weight of μ = 1.3, this gives total particle number densities of n = 46, 460, and 4600 cm⁻³. Each cloud has a column density consistent with observations. Observations show a strong positive correlation between the mass and density of GMCs in PHANGS galaxies (Sun et al. 2022), suggesting mass and density should be varied together.

The initial clouds must be in pressure equilibrium with their surroundings to avoid unphysical shocks from pressure imbalances at their surfaces. The M4 and M5 clouds are in the pressure regime where there is a stable two-phase medium at solar metallicity and Milky Way background UV field (Field et al. 1969; Wolfire et al. 2003), meaning there is a set of temperatures for the cloud and background for a given cloud density where the cold dense cloud and the warm ambient medium are both in thermal equilibrium at equal pressure. The cloud temperatures for the M4 and M5 clouds are T_cl = 103 K and 28 K, respectively, and the corresponding background temperatures and number densities are T_amb = 9000 and 4000 K, and n_amb = 3 and 1 cm⁻³. The M6 cloud, however, is at a high enough pressure that a two-phase medium no longer exists. Only the cold phase can be in thermal equilibrium. This means that the low-density envelope of the M6 cloud is inherently not in thermal equilibrium. To minimize the pressure imbalance with the core, we therefore raised the background density to n_amb = 100  cm⁻³. Both the cloud and background medium for M6 are at a temperature of T_cl = T_amb = 50 K.

The initial conditions described in the rest of this section and summarized in Table 2 apply to all three clouds. The clouds have a Gaussian density profile (Bate et al. 1995; Goodwin et al. 2004) with the standard deviation set such that the ratio of the cloud’s central to edge density is 3:1. The simulation domain is a cube of half-width $R_{\mathrm{box}}=20\mathrm{pc}$ with outflow boundary conditions. The outflow boundaries do allow gas flow onto the grid from ghost zones if the velocity at the edge of the grid is directed inward. The inflow of gas from the boundary is minimal: more gas exits the simulation than enters in all of our runs. Inflow was allowed, though, to prevent vacuums from forming at the boundaries. We used three refinement levels, yielding cell sizes that range from Δx_min = 0.3125 pc to Δx_max = 1.25 pc. Refinement and derefinement of the grid was determined by the Jeans criterion described in Sect. 2.1 and based on temperature and pressure gradients. The latter trigger refinement when the adapted Löhner (1987) estimator² of temperature or pressure exceeds 0.98 and trigger derefinement when the estimators drop below 0.6. We are interested in global formation properties of clusters rather than the fragmentation of the cloud or the origin of the IMF, so we find the chosen resolution to be sufficient.

We initially imposed a Kolmogorov (1941) turbulent velocity spectrum on all the gas in the domain. The peak Mach numbers for the turbulent spectrum are Ma = 30.3 for M6, Ma = 12.9 for M5, and Ma = 2.1 for M4. The same random seed was used to generate the turbulent velocity spectrum for all three clouds. This ensured the same turbulent collapse patterns, minimizing differences in the formation, location, and morphology of dense cores. From the edge of the cloud to the domain boundary, we linearly tapered the magnitude of the turbulent velocities from 100% to 25%. This tapering does not affect the low-density ambient background of the M4 and M5 cloud, but helps with stability in the M6 cloud by mixing the border of the cloud, where there is a small pressure jump.

The sink accretion radius and threshold density, derived in Sect. 2.1, are r_sink = 2.5Δx_min = 0.78 pc and ρ_sink = 8.35 × 10⁻²¹ g cm⁻³. This gives an initial sink mass resolution of m_sink = 245 M_⊙, meaning that when a sink initially forms it will accrete and form approximately m_sink worth of stellar mass, given the sink’s threshold density and accretion radius. The IMF sampling mass range is 0.08–100 M_⊙. The lower end is the hydrogen-burning limit, while the upper end is the most massive star thought to form in a star cluster with stellar mass ≈10⁴ M_⊙ (Weidner et al. 2009). This is the expected stellar mass limit for a cluster similar to M4, so we chose this value as a fixed parameter for consistency between the three clusters.

The critical virial ratio for stability is α_v = E_kin/|E_pot| = 0.5, below which collapse occurs. Massive clouds tend to be sub-virial, with clouds of 10⁶ M_⊙ observed to have virial parameters of α_v ≈ 0.05 − 0.35 (Kauffmann et al. 2013), though some surveys see super-virial massive clouds (see Fig. 2 of Chevance et al. 2023). We note that these values have been converted from the different virial parameter definition in Kauffmann et al. (2013). Therefore, we chose a fiducial virial parameter value of α_v = 0.15 for our models to promote rapid onset of collapse.

Magnetic fields are prevalent in the interstellar medium (Crutcher et al. 2003) and affect the collapse of GMCs and subsequent star formation. Although they are not the dominant factor in determining how star formation proceeds within a cloud, their presence has been shown to alter the fragmentation of cores (Price & Bate 2008; Peters et al. 2011) and slow down the global evolution of the cloud (Heitsch et al. 2001). With a strong enough field, clouds can be supported against gravitational collapse (Heiles 1976), although generally observed magnetic fields are not strong enough to inhibit collapse (Klessen & Glover 2016). The critical value of the mass-to-flux ratio for a cloud to be supported by magnetic fields against gravitational collapse is given by (Mouschovias & Spitzer 1976; Mouschovias 1991)

$$\begin{array}{c}\hfill [\frac{M}{\mathrm{\Phi }}{]}_c=\frac{\zeta }{3\pi }\sqrt{\frac{5}{G}}=490\frac{\mathrm{g}}{{\mathrm{Gauss}\mathrm{cm}}^2},\end{array}$$(3)

where G is the gravitational constant and a correction factor ζ = 0.53 for a uniform sphere is used (Strittmatter 1966).

In our simulations, each cloud’s initial magnetic field $\mathit{B}=B_z\widehat{z}$ is uniform in z and decreases radially in the x–y plane, following the mid-plane density ρ(x, y, z = 0) as

$$\begin{array}{c}\hfill B_z(x,y)=B_{0}exp[-(x^2+y^2)ln(3)/R_{\mathrm{cl}}^2],\end{array}$$(4)

with $B_0=0.185,1.85,18.5\mathrm{\mu }\mathrm{G}$ for the M4, M5, and M6 clouds, respectively. These values match observations for M5 and M6, while the field is a factor 10 weaker for M4 (Crutcher et al. 2010). The integrated magnetic flux Φ = 2πB₀R_cl²/(3ln(3)), so all clouds have an initial mass-to-flux ratio M_cl/Φ = 4.5 × 10⁴ g Gauss⁻¹ cm⁻² much larger than Equation (3). The initial magnetic fields are thus weak and do not inhibit collapse in any of our simulations.

3. Results

3.1. Cluster formation overview

At the onset of the simulation, each cloud begins to gravitationally collapse. Turbulent velocities fragment the cloud and create overdense hubs and filaments. Because the same random seed was used in all three clouds to generate the initial turbulent velocity spectrum, the web of dense gas is the same for each cloud. This means that the spatial distribution of star formation is similar for all three clouds. This can be seen in the time evolution of the three clouds in Figure 1. The first stars all form in the largest over-density in the middle of the bottom of the cloud. Then, more stars form along the filaments of the dense cloud forming a V shape. The M5 and M6 clusters in particular look very similar in terms of sub-clustering and merging. The M4 cluster forms significantly fewer stars and therefore fewer sub-clusters.

By a free-fall time t_ff, the sub-clusters in M4 and M5 have mostly merged, forming a single central spherical cluster. The M6 model is still forming stars in various sub-clusters and has not assembled its main central cluster yet. By looking at the spatial distribution of sub-clusters and the density of the gas, one can see that stellar feedback becomes most efficient once the sub-clusters have merged into a single cluster. Whether feedback is only strong enough to disperse gas when clustered or this is coincidental with the timing of feedback needs further examination, but this is outside of the scope of this introductory paper. Low density bubbles begin to occupy a significant fraction of the cloud volume once the central star cluster has been assembled.

Once most of the stars have formed, the efficiency of stellar feedback for dispersing the natal gas varies greatly for the three cloud masses. The final row in Figure 1 shows the M4, M5, and M6 systems at 1.2t_ff. At this point, only the M4 and M5 clusters have blown large bubbles. The feedback from the M6 cluster has hardly slowed the collapse of the densest gas, and rapid star formation continues. The M4 cloud has dispersed nearly all of the remaining gas, and star formation has halted completely.

3.2. Visualizing cluster morphology

The complex 3D structure of star clusters is hard to visualize using 2D plots. Figure 2 shows a still of an interactive plot of the M4, M5, and M6 simulations after one free-fall time generated with Plotly (Plotly Technologies Inc 2015). After downloading the HTML file available in the online version of this paper, readers can zoom, pan, and rotate for a complete look at the morphology of each cluster. The color is an isosurface of the gas density in log scale, and the points are stars with sizes scaled to their stellar radius. This tool makes it clear just how non-spherical these clusters are. Comparing the still of the interactive plot to the slice plots in Figure 1, one can already extract much more information on the system’s morphology.

Zooming into the core of each cluster shows the immense stellar density of the M6 cluster, whereas in the M4 cluster one can easily distinguish individual stars. The gas is also far less dense in the M4 system compared to M5 and M6.

The shape of the M6 cluster is highly irregular. Stemming from the largest cluster, one can see a row of sub-clusters forming along a filament. Branching perpendicularly off this filamentary cluster are two more star forming filaments in a configuration resembling the letter “F”. The M5 cluster has a shape congruous to the shape of M6, but with fewer stars bridging the gaps between clusters in the main filament. The M5 cluster also has only one finger perpendicular to the main filament, which contains many fewer stars than the fingers of the M6 cluster. The M4 cluster is much less dense, with its few stars outlining the same core filament cluster seen in M5 and M6. However, in M4 sub-clusters can no longer be distinguished. The M4 sub-clusters already merged into a singular central cluster.

3.3. Star formation history

The global properties of the star clusters that form from the M4, M5, and M6 clouds over the period of star formation are shown in Figure 4 as a function of time in units of the global free-fall time of the cloud t_ff (see Table 1). For reference, this same figure is shown as a function of physical time in Appendix A. We analyze these properties and discuss how they compare across the three clouds to assess the effect of the initial cloud mass and density on the resultant cluster properties. Table 3 highlights the final properties of the clusters formed in the three models.

Fig. 4. Global properties of the clusters and gas over time for models M4 (orange), M5 (maroon), and M6 (blue-violet) in units of free-fall time tff of the initial cloud (given in Table 1). From top left to bottom right: (a) SFR, where the transparent lines show the SFR at each star formation event, and the solid lines give the SFR smoothed using a Gaussian filter with σ = 0.005tff. (b) Instantaneous and integrated SFEs of the clouds, where ϵinst = M⋆/(Mgas + Msink + M⋆) and ϵint = M⋆/Mcloud = ϵ⋆. (c) Most massive star formed. (d) Number of formed stars. Dashed line: actual number of stars that would form from sampling the IMF given the amount of gas mass collected for star formation by sink particles. Solid line: number of stars followed in TORCH after the sampled stellar population below 4 M⊙ has been agglomerated. Dotted line: number of stars above 20 M⊙ on the grid that are generating feedback. The number of stars can drop due to SN, mass loss, or exiting the grid. (e) 3D stellar velocity dispersion. (f) Half-mass radius of the entire star cluster. (g) Total mass (dotted line), mass of stars (dashed line) and gas (solid line) on the grid. (h) Virial parameter of stars (dashed line) and gas (solid line), where αv = 0.5 is the equilibrium value. (i) Fraction of mass bound for stars (dashed line) and gas (solid line).

Our results suggest that star formation in a GMC is a fast and efficient process regardless of initial cloud mass and density, with all three clouds converting at least 30% of their initial gas into stars within an initial free-fall time. Star formation becomes faster and more efficient as the mass and density of the GMC increases.

Stars begin forming in the M5 and M6 clouds at t = 0.3 t_ff, while in the M4 run it is delayed until t = 0.5 t_ff. Because of the turbulent field, regions of the clouds have locally shorter free-fall times leading to star formation earlier than the global free-fall time. The duration of star formation is the shortest for the M4 cloud, lasting t_sf = 0.7 t_ff. The M5 cloud forms stars for a longer period in terms of its initial free-fall time, t_sf = 1.3 t_ff. M6 is still forming stars 1.3 t_ff after the onset of SF.

The SFR (Fig. 4a) increases with cloud mass. The peak SFR for the M4, M5, and M6 clouds are SFR_peak = 0.4, 5.5, and 392 M_⊙ yr⁻¹ respectively. The average SFRs also increase with mass, with values of SFR_ave = 0.02,  0.06, and 1.5 M_⊙ yr⁻¹.

The integrated star formation efficiency (SFE; Fig. 4b) we discuss here is given by the ratio of stellar mass formed to the initial gas mass of the cloud ϵ_⋆ = M_⋆/M_cloud. The instantaneous SFE is the ratio of stars to the total mass on the grid at that point in time ϵ_inst = M_⋆/(M_gas + M_sink + M_⋆).

For all three clouds, the instantaneous and integrated SFE closely follow each other at early times. In the M4 track, they diverge at t ∼ 1t_ff: ϵ_inst increases from 35 to 70%, while ϵ_⋆ levels out due to gas expulsion. When all gas is fully expelled from the domain all three values will converge to ϵ_inst = 100%. The instantaneous SFE for M5 is beginning to increase toward 100% as the integrated SFE levels out. The ϵ_inst for M6, however, levels out ≈10% lower than ϵ_⋆. This is due to the higher background density and some inflow from the boundary. Inflow is expected to occur for systems like M6, so this suggests observed SFEs of massive embedded clusters may be lower than the conversion ratio of gas to stars from the original cloud. For consistency, the best estimate for simulated SFE is the final integrated SFE value, and this is the value we use for all further comparisons to observations. We delve into the limitations of comparing observed and modelled SFEs in Sect. 4.1.

The M4 cloud is representative of typical GMCs at the solar circle, and its integrated SFE lies just over the top of this range at 36%. Typical SFE values of nearby clusters in the Milky Way lie between 0.1–0.3 (Lada & Lada 2003). In the higher mass clouds of M5 and M6, the SFE is much higher. The M5 cloud converted 65% of its gas into stars, and the M6 cloud converted 89% of gas into stars. This suggests that the SFE in high-mass, high-density environments can be much higher than seen in low-mass local clusters.

The free-fall time becomes so short in these high-mass clouds that the stellar feedback simply does not act quickly enough to stop star formation before most of the gas has formed into stars. Free-fall times of dense environments that are shorter than the development times for winds and SNe have indeed been shown to result in high SFE (Dekel et al. 2023). However, the more crucial factor may be that the total force budget of feedback from winds and radiation is not enough to surpass the gravity from gas and stars in dense, high-mass embedded clusters. This force-balance argument is supported by our results presented in Sect. 4.1 where the stellar feedback in a 1D model of M6 fails to expel gas from the embedded cluster’s potential well. Regardless of free-fall time, in high density clouds the total feedback energy won’t equal gravity until over half the cloud mass is converted to stars. However, if the free-fall time becomes so long ($t_{\mathrm{ff}}\ge 10\mathrm{Myr}$) that SNe become a dominant form of feedback during the primary epoch of star formation, the energy from frequent SNe may start to overpower gravity and affect the SFE. We have not explored this regime, as such a high free-fall time for a 10⁶ M_⊙ cloud is rare.

Figure 4c shows the mass of the most massive star that has formed from random draws from the IMF. By a free-fall time, each cloud has already formed the most massive star in its cluster. We find that the mass of the most massive star increases with cluster mass. For the M5 and M6 clouds, the most massive star is at the maximum sampling mass of 100 M_⊙, while the M4 cloud’s most massive star is around 70 M_⊙. This is a stochastic effect; as more stars are sampled from the IMF, you will eventually sample the most massive star in the distribution. This reproduces the effect suggested by Weidner et al. (2009) and Yan et al. (2023) that the cluster mass limits the most massive star mass. In each cloud, it is interesting to note that each instance of the formation of a very massive star, that is, above 40 M_⊙, correlates with a slowing of star formation indicated by a reduction in the SFE slope.

The M6 cloud forms from IMF sampling ∼10⁶ stars, M5 forms ∼10⁵ stars, and M4 forms ∼10⁴ stars, shown in Figure 4d. With agglomeration, the number of stars in the simulation are about 10% of these numbers, so the improved version of TORCH with PETAR can simulate clusters of > 10⁵ individual stars.

3.4. Cluster evolution

The evolution of the global properties of the formed star clusters occurs quite similarly for all three clouds, but the magnitude of their values depends greatly on the cloud’s initial mass.

The stellar velocity dispersion (Fig. 4e) generally increases with initial cloud mass. The velocities of stars increase at a slow pace before leveling out after 1 t_ff. For the M4 cluster, the velocity dispersion levels out at just 1.0 km s⁻¹. At late times, the velocities of the stars begin to slightly decrease in M4. This decrease correlates to the increasing half-mass radius of M4, indicating the star cluster is expanding and the stellar velocities are slowing. The M5 cluster reaches a velocity dispersion of 5 km s⁻¹ and the M6 cluster has a velocity dispersion of 20 km s⁻¹. The deeper potential wells of the higher mass clusters, going as the square root of the mass for these similar sized objects, drive the higher velocity dispersions, although the measured dispersion increases somewhat faster with mass than the potential well depth. In the case of M6, the potential well depth exceeds the sound speed of ionized gas, preventing gas from escaping even after ionization. For M6, the average sound speed in cells where the gas is fully ionized is c_ion= 16.4 km s⁻¹.

The evolution of the half-mass radius R_1/2 of all the stars in the cluster (Fig. 4f) seems to be split into a high and low-mass regime. The M5 and M6 clusters follow the same track closely. From 0.3 t_ff to 1.0 t_ff, R_1/2 increases to a peak of ∼3.75 pc at ∼0.6 t_ff then goes down to ∼1 pc (M5) and ∼0.5 pc (M6) and begins to level out. The similarities in the evolution of the two clusters are most likely due to the fact that both clouds at early times form enough stars for distinct sub-clusters to form and merge. The sub-clusters have formed in the same places so both clusters peak at roughly the same R_1/2. All of the stars in M5 and M6 remain bound, suggesting the clusters are relaxing into gravitational equilibrium. Longer runs following just the stars after gas dispersal will ultimately be needed to demonstrate this. M6 is a smaller and denser cluster than M5, likely due to the much deeper potential well of M6. Similar stellar densities are observed in the super-star clusters (SSCs) of starburst NGC 253 (Rico-Villas et al. 2020) with sizes of R < 1.7 pc and stellar masses of 10^4 − 6 M_⊙.

The M4 cluster grows slightly differently. It increases to a peak of R_1/2 = 2 pc at 1 t_ff, decreases to R_1/2 = 1 pc, then linearly increases to 2.5 pc by 1.3 t_ff. The M4 cluster is expanding, but 85% of its stars remain bound, so complete dissolution has not yet occurred (Fig. 4i).

The onset and duration of gas dispersal from the star clusters depends strongly on the initial mass and density of the cloud. Figure 4g–i shows the time evolution of the mass, virial parameter, and bound mass fraction of the gas and stars. With these three plots we can track the degree of gas dispersal. From the mass plot we see that by the end of star formation, the M4 and M5 clusters have expelled a significant fraction of the initial cloud. Only 10% of the original gas mass remains in M4 and < 10% in M5. Both M4 and M5 are well on their way to full gas expulsion, as in both clusters < 10% of the gas still on the grid is bound. The gas in the M4 and M5 systems become super-virial by a free-fall time. The gas in M5 takes ∼10% longer to become unbound, but progresses identically to the M4 gas. The gas in M6 differs significantly as it remains sub-virial even beyond one free-fall time. The potential well created by the massive cluster is enough to keep the remaining gas infalling, suggesting star formation is not yet quenched in the M6 cluster.

The only star cluster that is starting to disperse is the M4 cluster (Fig. 4i). The other two star clusters remain fully bound. The stars in M5 and M6 remain sub-virial, while the stars in M4 just reach virial equilibrium by the final simulation time (Fig. 4h). The dispersal time of the gas and stars increases with initial cloud mass as there is more gravity for the stellar feedback to counteract. Although massive clusters have more stars injecting feedback, the increasing gravity overpowers the feedback. At high densities, where the potential well depth exceeds the sound speed of ionized gas, ionization feedback cannot disperse gas, while the short free-fall time assembles dense gas more quickly than feedback can work against gravity.

4. Discussion

4.1. Limitations of comparing observed and modelled star formation efficiencies

We compare the SFEs of our modelled clusters to observations by using the integrated quantity ϵ_⋆, the stellar mass divided by the initial cloud mass. This is the total fraction of cloud mass that has been converted into stars. However, this value is impossible to calculate for observed clusters, as the only information available is how much gas and stars are present within a certain area. Thus, observations of star clusters only quote the instantaneous value ϵ_inst. There is a 3–4 order of magnitude spread in the observed SFE and SFR of Galactic GMCs (Lee et al. 2016).

A proper comparison between simulated and observed SFEs requires accounting for the amount of ongoing star formation, and determining whether the embedding gas is collapsing or dispersing. Comparisons done without accounting for the evolutionary stage of the cluster are misleading. The apparent SFE ranges from 0–100% over the lifetime of every star cluster that reaches a gas-free state, which may explain the spread in observed SFE. This issue is starting to be explored. Geen et al. (2017) suggest techniques for converting observed to integrated SFEs, although conclude that this is non-trivial and find errors up to a factor of 10. They find overall higher observed SFEs than integrated SFEs when they applied observational techniques to their simulations. On the other hand, Grudić et al. (2019) find lower observed than integrated SFEs in their models due to the inaccuracies of techniques for estimating stellar mass. One example of this comes from using only the young stellar population as a tracer for stellar mass, which underestimates the total stellar mass. They also discuss variability in observed SFEs due to the changing SFE over the course of a GMC lifetime, from first star formation to gas dispersal.

Further studies must be dedicated to outlining a systematic way to convert between observed SFEs and the final integrated SFEs we define in our models. Until then, direct comparisons should be interpreted with caution.

4.2. Observations

Galactic surveys of embedded clusters in the Milky Way typically find the SFE to be ≲30% (Lada & Lada 2003), with some studies finding lower values of ≲8% (Evans et al. 2009; Peters et al. 2011). The M4 cloud, which is a good representative of galactic GMCs³, agrees well with this SFE albeit at the high end of observed values. This could be due to the low virial parameter used, which is appropriate for M6 but lower than the average observed value of α_v = 1 seen in clouds similar to M4. The missing radiative feedback from stars < 20 M_⊙ could also be a factor causing the high SFE of M4 given the low-density of the cloud. Feedback contributions from low-mass stars could be a key factor in quenching star formation in Milky Way-like clouds.

The higher mass clouds similar to M5 and M6, however, have SFEs well above 30%. While there are no Milky-Way analogs to the M6 cloud, there are a few for M5. There is the W43 GMC with 1.32 × 10⁵ M_⊙ of gas within R ∼ 10 pc (Lin et al. 2016), similar to the M5 cloud with R = 11.7 pc. The W49 star forming region has a central YMC with stellar mass ≳ 5 × 10⁴ M_⊙ and gas mass ∼ 2 × 10⁵ M_⊙ and ∼ 1.1 × 10⁶ M_⊙ within 6 and 60 pc respectively (Galván-Madrid et al. 2013). This gives a current SFE of 20% in the inner region. W49 has ongoing star formation, and its embedded gas cloud is twice as massive and 15 times as dense as our M5 model. Based on our results for M5 and the fact that SFE increases with density, we predict W49 will exceed the SFE found for M5 of ϵ_⋆ = 65%.

Though conditions required to form the M6 cloud are not observed in the Milky Way, they are present in other galaxies. Starburst galaxies have been observed to host SSCs with SFEs of ϵ_inst ≈ 50 − 80%. These SSCs cover a size and mass range comparable to our models (see Table 4).

Additionally, the disks of gas-rich high redshift galaxies can be violently unstable and are thought to form clouds similar to M6 (see Tacconi et al. 2020). We can now directly observe the high redshift environment of forming GCs with JWST. (Li et al. 2024) predicts JWST will find feedback-free starbursts (Dekel et al. 2023), which are massive galaxies at z ≳ 10 with high SFEs due to dense gas with free-fall times ≤ 10 Myr forming stars effectively free of stellar feedback. Recent JWST observations uncovered “younger” populations of GCs in galaxies at redshift z = 0.38 (Harris & Reina-Campos 2023), and others are expected to observe GCs up to z = 1 without lensing (Reina-Campos & Harris 2024).

With lensing, clumps that are likely proto-GCs can be observed at redshift z > 1 (Adamo et al. 2023; Claeyssens et al. 2023). One such proto-GC candidate was found through lensing at z ∼ 6 with ≲ 10⁶ M_⊙ and a core radius of R_c < 13 pc (Vanzella et al. 2019). Another more massive bound YMC 3 Myr old was found at z = 2.37 with ∼ 10⁷ M_⊙, and R ∼ 8 pc (Vanzella et al. 2022a). A strongly lensed galaxy at z = 4 contains three bound YMCs each younger than < 30 Myr with masses between (0.7–4.0) × 10⁶ M_⊙ and radius estimates of 3–20 pc (Vanzella et al. 2022b). The Sunrise arc is a strongly lensed z ≈ 6 galaxy found to contain six YMCs with masses ∼ 10⁶⁻⁷ M_⊙, radii of ∼1–20 pc, and ages 1–30 Myr (Vanzella et al. 2023). Most of these recently discovered YMCs or proto-GCs are analogous in size and mass to the M6 cluster or larger. Dense regions of prolific star formation that form these objects seem pervasive in the early Universe, and more will surely be discovered as more JWST data arrives. Due to their similar properties, we argue that these clusters formed in the same manner as M6.

Another situation that can form GMCs similar to M6 is major galaxy mergers with small mass ratios. Tidal interactions of major galaxies are linked to bursts in star formation (Larson & Tinsley 1978; Lonsdale et al. 1984; Barton et al. 2000; Ellison et al. 2008; Renaud et al. 2019). Since most massive galaxies are believed to undergo at least one merger in their lifetime, this is not a rare occurrence. Galaxy mergers have been suggested as the progenitors of YMCs and younger GCs (Ashman & Zepf 1992; van den Bergh 2001). We note that this only applies to major galaxy interactions: minor galaxy interactions with large mass ratios produce little to no enhancement of the overall SFR (Cox et al. 2008; Tress et al. 2020).

In the interacting Antennae galaxies, the Firecracker cloud, which resembles M6, was observed by Whitmore et al. (2014). Finn et al. (2019) constrains its mass and characteristic radius to (1–9) × 10⁶ M_⊙ and 22 pc. The Firecracker cloud is in the very early stages of star formation, as it is estimated to have only formed M_⋆ ≲ 10⁴ M_⊙. This is less than 10% of the expected stellar mass of the final star cluster (Johnson et al. 2015). These observations show that progenitor clouds similar to M6 can form before any significant amount of star formation occurs.

A survey of the molecular clouds in the Antennae galaxies done by Wei et al. (2012) revealed two populations of MCs, with a distinct break in the differential mass function at $log(M_{\mathrm{cloud}}/{\mathrm{M}}_{\odot })=6.5$. Clouds above this mass were found in the regions of intense star formation, while the lower mass clouds were in more dormant regions. The large velocities seen in the high SF regions suggest compression by shocks, supporting the idea that galaxy mergers lead to high-mass GMCs that become sites of extreme star formation.

Finn et al. (2019) measured the velocity dispersion in the Firecracker cloud and found it to be neither in virial equilibrium nor free-fall. They conclude that there must be a high pressure background to contain the gas at such high densities in equilibrium.

We compare this velocity dispersion to those in our clouds over time to test whether we reach such high velocities through the addition of stellar feedback to free-fall collapse alone or whether a high pressure background is indeed needed. The results of this comparison are shown in Figure 5. We plot the size-line width coefficient σ_v²/R against surface density Σ = M/πR² for each of our clouds to compare with the observations from Finn et al. (2019). The observations are shown by the white diamonds in Figure 5 corresponding to 4 aperture sizes they used for R: 6.4, 15, 26, and 37 pc. We use four smaller aperture sizes: 5, 10, 15, and 20 pc, as our cloud is half the radius of the Firecracker. We plot four times for each of our simulations corresponding to t = 0.0, 0.5, 1.0, 1.25 t_ff, with a line connecting the points showing the set of apertures. The apertures for each set increase from right to left with decreasing surface density. The colors of the points indicate the amount of stellar mass formed, colored white when no SF has occurred yet. The sizes indicate the initial gas mass of the cloud, with the smallest being the M4 simulations and the largest being M6. The two labelled lines on the plot indicate the analytical conditions for free-fall and virial equilibrium (see the discussion of Fig. 2 in Ballesteros-Paredes et al. 2011). We include these lines for easier comparison to the observations presented by (Finn et al. 2019) rather than for inferring the state of the system. Unintuitively, velocity dispersions in a state of free-fall exceed those of virial equilibrium.

Fig. 5. Size-line width coefficient versus surface density for the Firecracker cloud (Finn et al. 2019) and the M4, M5, and M6 clouds. The four red diamonds are observations, done using apertures (1–4) with radii R = 6.4, 15, 26, and 37 pc. The other points are our simulations, with connected points corresponding to apertures (1–4) of radii R = 5, 10, 15 and 20 pc. The colors indicate how much stellar mass has been produced (with empty points indicating no star formation). The lines correspond to virial equilibrium and free-fall as labeled (see Fig. 2 of Ballesteros-Paredes et al. 2011).

Though the Firecracker cloud is as massive as M6, its surface density is lower: more comparable to M5. This is why we see an overlap between M5 and the Firecracker cloud when M5 has formed 4 × 10⁴ M_⊙ of stars, which is on the order of the stellar mass estimated to have already formed in the Firecracker cloud of ≤ 10⁴ M_⊙. This aligns with the idea that surface density is more influential than mass in the formation of star clusters. This also supports the possibility that the Firecracker gas velocities could be caused by the contribution of stellar feedback to the velocity dispersion in addition to free-fall collapse. The fast dynamical evolution of our models suggests that these are not equilibrium objects, making it unnecessary to invoke a high-pressure background to keep the cloud from expanding. This suggests objects like the Firecracker cloud can form from collapse with observed velocity dispersions without invoking a high pressure background medium.

4.3. Other simulations

For massive star clusters to form, they must survive the epoch of gas dispersal and remain bound. An analytical model predicts a positive correlation of SFE to SFR and initial cloud mass (Zamora-Avilés & Vázquez-Semadeni 2014). At high enough surface densities, stellar feedback cannot compete with star formation. Numerical studies done by Geyer & Burkert (2001) found that if the stars are initially in virial equilibrium with the remaining gas, only clusters with SFE ≥ 50% remain bound against the outflow of the gas. Li et al. (2018) finds in their cosmological galaxy formation models that though galactic properties are unaffected by varying ϵ_⋆, the properties of star clusters are affected. In particular, they found that the initial bound fraction of stars increases with ϵ_⋆ and cloud mass. Farias et al. (2023) ran cluster formation models from 2 × 10⁴ M_⊙ clouds and finds that SFE and gas expulsion time correlate with global bound fraction, with all SFEs ≤ 20% and all bound fractions ≤40%. Although it is still possible to form bound clusters with low SFE, these studies imply massive bound star clusters were most likely formed with high SFEs.

Menon et al. (2023) also finds high SFEs of ∼80% for 10⁶ M_⊙ clouds with feedback in the form of radiation pressure solved using a variable Eddington tensor approach as opposed to our ray-tracing method. In this density regime, radiation pressure is the dominant feedback mechanism⁴. They conclude that radiation pressure simply cannot regulate star formation for clouds with surface densities Σ ≳ 10³ M_⊙ pc⁻².

Our results are more constraining: We included more feedback physics, and we still achieved SFEs of ϵ_⋆ > 80%. Our M6 cloud is above this surface density with Σ = 2.3 × 10³ M_⊙ pc⁻². They also tested a larger 10⁶ M_⊙ cloud with roughly the same surface density as our M5 cloud, and find an SFE of ϵ_⋆ ≈ 60% comparable to the SFE of our M5 cluster of ϵ_⋆ = 65%.

Other simulations of massive star cluster formation with initial cloud mass of 10⁶ M_⊙ find high SFEs of ∼65% (Grudić et al. 2018) and 38% (Kim et al. 2018) for surface densities of $\mathrm{\Sigma }=1.27\times {10}^4{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}$ and $\mathrm{\Sigma }=500{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}$, respectively. Kim et al. (2018) finds an even higher SFE of 51% for a 10⁵ M_⊙ cloud but with a surface density of $\mathrm{\Sigma }=1.27\times {10}^3{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}$. Cluster models in Kimm et al. (2022) reached ϵ_⋆ = 50–72% from an initial cloud of 1.4 × 10⁶ M_⊙ and $\mathrm{\Sigma }=647{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}$ despite SN and radiation feedback. Another cluster, modeled with radiation feedback by Fukushima & Yajima (2021), reached ϵ_⋆ = 78% from a 10⁶ M_⊙ cloud with $\mathrm{\Sigma }=3.2\times {10}^3{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}$. This study also finds that bound cluster formation only occurs with ϵ_⋆ ≥ 30%. Recent models of ${10}^6{\mathrm{M}}_{\odot }$ clouds with R = 10 pc, α_v = 0.1, and Z = 0.2 Z_⊙ reached ϵ_⋆ = 50% (Fujii et al. 2024). The SFE in all of these studies increases strongly with surface density and slightly with initial cloud mass. Our results combined with those from previous models provide evidence that the formation of bound YMCs requires not only a high cloud mass but also, and more importantly, a high surface density.

Simulations of star clusters forming from clouds similar to our M5 cloud resulted in lower SFEs of 10–30%. The main difference between these models and ours is the use of sink particles to represent sub-clusters with combined feedback compared to our tracking of feedback from individual massive stars. A cloud modelled by He et al. (2019) with an initial mass of 10⁵ M_⊙, peak number density $n=1.8\times {10}^3{\mathrm{cm}}^{-3}$, metallicity $\mathrm{Z}={\mathrm{Z}}_{\odot }$, a higher virial parameter α_v = 0.4, and stellar feedback only through UV radiation reached a SFE of 13.7% by 6 t_ff. The cloud in Ali (2021) is the same mass and metallicity as M5, almost the same radius, $R_{\mathrm{cloud}}=11.9\mathrm{pc}$, but is initially super-virial, with α_v = 1. By 0.75 t_ff, the SFE reached only 10% while M5 reached ϵ_⋆ = 20% by this time. This difference may come from the different initial virial parameters, which prolongs the formation of the cluster in Ali (2021). Another possible cause is different feedback models. Injecting stellar feedback from entire sub-clusters rather than individual stars could artificially strengthen the effect of feedback resulting in a lower SFE.

Fujii et al. (2021) presents star-by-star cluster models with feedback in the form of radiation, radiation pressure and stellar outflows. For an initial cloud of ${10}^5{\mathrm{M}}_{\odot }$, $R=20\mathrm{pc}$, and α_v = 0.25 they find ϵ_⋆ = 40%. The same cloud with a larger virial parameter of α_v = 1.0 only reached ϵ_⋆ = 40%. Their sub-virial model agrees with our findings for M5. The super-virial model with a lower SFE further indicates that the high SFEs we find are possibly due to the low initial virial parameter, particularly for lower density clouds.

A colliding flow model of star formation in GMC environments described in Colín et al. (2013a) with individual star formation and ionizing radiation found SFEs of ϵ_⋆ = 10–30% depending on the degree of concentration by the flows. The two cylindrical streams were very large, with r = 64 pc and $\ell =112\mathrm{pc}$ and rarefied, with $n=1{\mathrm{cm}}^{-3}$, with the total mass in the two streams equalling 9 × 10⁵ M_⊙. The different initial conditions hinder a direct comparison to our SFE values, but reaching high SFEs from low density flows aligns with our results for M4.

Simulations have broadly found star formation to be suppressed with each additional form of stellar feedback included in the model (see Dale 2015). The exclusion of protostellar jets in our feedback model may artificially raise the SFR, as they contribute to the dispersal of gas around even low-mass stars at small scales (Chevance et al. 2023). Due to the quantity of low-mass stars and the collimated shape of the outflow, jets are drivers of turbulence at large scales in GMCs (e.g., Nakamura & Li 2007; Federrath 2015; Appel et al. 2022). These models do show that jets are an important factor in slowing the growth rate of the integrated SFE, though the final SFE is not known due to the duration of the simulations (Federrath 2015). Guszejnov et al. (2021) performed simulations of star-by-star cluster formation from 2 × 10⁴ M_⊙ clouds with stellar feedback, including protostellar jets as well as radiation, winds, and SNe. Simulations were repeated that isolated each form of feedback. They found jets to be important in regulating the growth of low-mass stars and constraining the IMF. Radiation and jets were the primary form of feedback that slowed star formation and dispersed the cloud. However, again the simulations were not run until the end of star formation, so the degree to which each effect changes the final SFE remains uncertain.

Despite their ubiquity, jets cannot prevent gas in high-density GMCs from forming stars eventually nor contain the power needed to disperse GMCs (see Chevance et al. 2023). The effect that jets would have on more massive clouds remains unclear. Although they may indeed slow star formation to observed values, as proposed by Chevance et al. (2023), analytic work by Matzner (2002) suggests that more massive clouds would be resistant to dispersal by jets, consistent with simulations by Guszejnov et al. (2022).

In order to verify the physical plausibility of the high SFE in M6 despite the large number of formed stars, we have directly compared our 3D results to a followup calculation using the one-dimensional (1D) code Winds And Radiation Pressure: Feedback Induced Expansion, colLapse and Dissolution (WARPFIELD; Rahner et al. 2019). This code models the effect of stellar feedback from young clusters on their natal gas cloud in spherical symmetry. WARPFIELD is designed to solve for the self-consistent motion of a 1D spherical gas shell evolving under the influence of feedback mechanisms including stellar winds, SNe, and radiation pressure, with consideration of gravity. We ran WARPFIELD using the same initial conditions as chosen for the M6 run (i.e., mass, density, temperature), with the addition that we varied the SFE from ϵ_⋆ = 0.1–0.9 in bins of 0.1, as shown in Figure 6. This varies the strength of the stellar feedback to test whether the M6 cloud would still be stable given amount of stars formed.

Fig. 6. WARPFIELD evolution of shell radius versus time with different initial SFEs with the same parameters described in M6. In all cases, stellar feedback is inefficient in dispersing the surrounding dense cloud, and the shell eventually undergoes re-collapse.

For all SFE values, the shocked gas eventually re-collapses. At this high density, the included feedback is not strong enough to completely disperse the cloud. The SFE of our M6 cluster most closely resembles the WARPFIELD runs with ϵ_⋆ = 80% & 90% in Figure 6. The WARPFIELD model reaches a maximum radius of $R\sim 55\mathrm{pc}$ at $t\sim 6.5\mathrm{Myr}$, and collapses back to $R=0\mathrm{pc}$ by $t\sim 11.25\mathrm{Myr}$. In the 1D model, all feedback occurs at a single point, so it is more effective than in our star-by-star 3D model, as in multiple dimensions channels that vent thermal energy can exist. Nevertheless, the gas still re-collapses, promoting further star formation. The expanding gas is not accelerated fast enough to escape the deep potential well of the massive cloud and the cluster that forms from it. This result supports the idea that the high SFE is due to the total feedback strength being weaker than gravity at these densities.

However, our results do suggest that more dispersed star formation leading to increased energy dissipation by radiative cooling may not even allow that much expansion. To resolve SFE well, feedback must be modelled for individual stars instead of for entire clusters. Approximating feedback as a sum for an entire cluster underestimates the SFE.

5. Conclusions

We performed numerical simulations of star cluster formation from gas clouds that run until star formation ceases or slows significantly due to stellar feedback dispersing any remaining gas. We tested initial cloud masses of 10⁴,  10⁵, and 10⁶ M_⊙ with radius R = 11.7 pc, holding all other characteristics of the initial cloud and simulation parameters the same. We analyzed the star formation histories and followed the evolution of the gas and forming star clusters. From this study, we conclude the following:

• Giant molecular clouds with surface density $\mathrm{\Sigma }\ge {10}^2{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}$ and mass $M_{\mathrm{cloud}}\ge {10}^5{\mathrm{M}}_{\odot }$ can form fully bound star clusters with stellar mass $M_{\star }\ge {10}^4{\mathrm{M}}_{\odot }$ with a high SFE ϵ_⋆ ≥ 65% over a short time t_sf ≈ 1t_ff, as seen by M5 and M6. The lower mass and density M4 cloud forms a cluster with a lower bound mass fraction of 60%.

• The Firecracker cloud in the Antennae galaxies, with a mass of 1–9$\times {10}^6{\mathrm{M}}_{\odot }$ and a radius of $22\mathrm{pc}$ (Finn et al. 2019), is a close analog to our M6 cloud, though with a surface density more closely matching our M5 cloud. From our results we can estimate that the Firecracker cloud will convert 65–85% of its mass into stars within a free-fall time and that it will form a YMC.

• It has been suggested that the Firecracker cloud must be surrounded by a high pressure medium to contain it because of its high surface density and size-line width coefficient σ_v²/R (Johnson et al. 2015; Finn et al. 2019). However, the M5 cluster reaches the same values by the time it forms $M_{\star }\approx {10}^4{\mathrm{M}}_{\odot }$ worth of stars, the same amount of stellar mass estimated to have formed in the Firecracker cloud. This suggests another possibility: Rather than being an equilibrium object confined in a high pressure environment, the Firecracker cloud is actually dynamically collapsing and forming stars, and the high velocity dispersion of the gas is from the combination of free-fall collapse and stellar feedback.

• Star formation from GMCs is capable of achieving up to 85% efficiency at high densities. Our M6 cloud is the most efficient of our models, converting ϵ_⋆ = 85% of its gas into stars. Even with hundreds of massive stars producing feedback, the short timescale of gravitational collapse for dense massive clouds renders the stellar feedback inefficient at slowing early star formation. However, even at much lower densities and masses, the M5 and M4 clouds achieved high SFEs of ϵ_⋆ = 65% and 36%, respectively. In dense, massive clouds, the total dispersing force of stellar feedback from winds and radiation cannot counteract the gravity from stars and gas until over half the cloud mass is converted into stars.

• The M4 cloud has a typical mass and size of Milky-Way GMCs. The SFE of M4 matches the maximum observed SFE values. This high SFE could be because of the low initial virial parameter of the cloud, or it could be due to the missing FUV radiation from stars $<20{\mathrm{M}}_{\odot }$. Alternatively, the exclusion of the protostellar jet feedback mechanism may be important for clouds similar to M4 clouds, as suggested, for example, by Chevance et al. (2023). Further studies must be done to constrain the effect of varying the virial parameter and including protostellar jets on integrated SFE.

• Star formation is fast in our models of clouds with low α_v. Regardless of the initial mass or density, the majority of star formation occurs within the first global free-fall time of the collapsing GMC. Collapse occurs and stars are produced so rapidly that stellar feedback is only prevalent and strong enough to clear dense gas from the cluster’s deep potential well after most of the cloud has formed into stars. The speed of star formation may depend strongly on the initial virial parameter and the inclusion of jets.

• A 1D stellar feedback model WARPFIELD was run using the same mass and density as the M6 simulation. In it, the gas re-collapses even for SFEs up to 90%. Even centralized feedback cannot expel the gas from the potential well of the massive cluster that forms. The WARPFIELD results indicate that the expanding gas shell for ϵ_⋆ = 85% collapses back to R = 0 by 11 Myr.

• Including feedback for individual stars rather than adding the total energy for the cluster at a single point is important for correctly constraining star formation histories. Modelling individual stellar feedback spreads the feedback energy enough to greatly reduce its effectiveness at clearing the natal gas because of the resulting enhanced radiative cooling. Models that add stellar feedback for the entire star cluster at a single point appear to overestimate the effect of the feedback on the gas and the star formation timescale and to underestimate the final SFE.

In conclusion, bound massive star clusters such as YMCs and GCs readily form from high-mass, dense GMCs. The GMCs can become this dense and massive naturally, even in the present day, as shown by the Firecracker cloud in the Antennae galaxies. In the early Universe, where galaxies were much more gravitationally unstable, these conditions would be much more common. The subsequent star formation from these dense high-mass clouds is highly efficient, converting ≥40% of the gas mass into stars within the first free-fall time of the initial cloud. The short timescales of star formation and/or the deep gravitational potential wells of dense, massive clouds render stellar feedback unable to significantly slow star formation, leading to integrated efficiencies as high as 85% for more massive clouds. After their formation, the clusters born in these environments remain bound after 90% of the gas is expelled.

Until recently, directly observing proto-GCs has been elusive. Now with JWST, observers have discovered five bound stellar clumps just $460\mathrm{Myr}$ after the Big Bang at $z\sim 10.2_{-0.2}^{+0.2}$ (Adamo et al. 2024). These clusters in the strongly lensed galaxy SPT0615-JD1 (alias the Cosmic Gems arc) have intrinsic masses of $\sim {10}^6{\mathrm{M}}_{\odot }$, half-light radii of $R_{\mathrm{eff}}\sim 1\mathrm{pc}$, and ages between 9 and 35 Myr. Roughly ∼60% of the total F150W flux of the galaxy comes from these five clusters, indicating that the predominant mode of star formation in these systems occurs in massive clusters. The resemblance between M6 and these objects is notable, indicating that the mode of star formation described in this study is a probable path for the formation of YMCs and proto-GCs in the present and early Universe, respectively.

Data availability

Simulation data and the interactive plot file is available for download at https://doi.org/10.5531/sd.astro.8.

Acknowledgments

B.P. would like to thank Gastón Escobar and Michela Mapelli for the helpful discussions that led to solving a persistent bug allowing continuation of the simulations. B.P. thanks Shyam Menon for many fruitful discussions of YMC formation, particularly regarding the battle between stellar feedback and gravity. We thank the referee for their useful comments and insights. B.P. was partly supported by a fellowship from the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). M.-M.M.L., B.P., E.P.A., and A.T. were partly supported by NSF grants AST18-15461 and AST23-07950. C.C.-C. is supported by a Canada Graduate Scholarship - Doctoral (CGS D) from the Natural Sciences and Engineering Research Council of Canada (NSERC). This work used Stampede 2 at TACC through allocation PHY220160 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants 21-38259, 21-38286, 21-38307, 21-37603, and 21-38296. The code development that facilitated this study was done on Snellius through the Dutch National Supercomputing Center SURF grants 15220 and 2023/ENW/01498863. S.A. acknowledges the support of NSF grant AST-2009679. M.W. acknowledges the support of NOVA project 10.2.5.12. R.S.K. and S.C.O.G. acknowledge financial support from the European Research Council via the ERC Synergy Grant “ECOGAL” (project ID 855130), from the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) “STRUCTURES”, funded by the German Excellence Strategy, and from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 50OO2206). The team in Heidelberg also thanks The Länd and the German Science Foundation (DFG) for computing resources provided in bwHPC supported by grant INST 35/1134-1 FUGG and for data storage at SDS@hd supported by grant INST 35/1314-1 FUGG. L.W. thanks the National Natural Science Foundation of China for support through grants 21BAA00619, 12073090 and 12233013, and the one-hundred-talent project of Sun Yat-sen University, the Fundamental Research Funds for the Central Universities, Sun Yat-sen University (22hytd09). The authors acknowledge Interstellar Institute’s program “II6” and the Paris-Saclay University’s Institut Pascal for hosting discussions that nourished the development of the ideas behind this work.

References

Appendix A: Stellar properties: Physical times

Figure A.1 reproduces Figure 4 using physical time rather than free-fall times to show global stellar properties over time. This demonstrates how much the duration of star formation shortens while its intensity increases as the cloud mass increases.

Fig. A.1. Global properties of the clusters and gas over time for models M4 (orange), M5 (maroon), and M6 (blue-violet) for comparison to Figure 4, where units of free-fall time (see Table 1) are used. From top left to bottom right: (a) SFR, where the transparent lines show the SFR at each star formation event, and the solid lines give the SFR smoothed using a Gaussian filter with σ = 0.005tff. (b) Instantaneous and integrated SFEs of the clouds, where ϵinst = M⋆/(Mgas + Msink + M⋆) and ϵint = M⋆/Mcloud = ϵ⋆. (c) Most massive star formed. (d) Number of formed stars. Dashed line: actual number of stars that would form from sampling the IMF given the amount of gas mass collected for star formation by sink particles. Solid line: number of stars followed in TORCH after the sampled stellar population below 4 M⊙ has been agglomerated. Dotted line: number of stars above 20 M⊙ on the grid that are generating feedback. The number of stars can drop due to SNe, mass loss, or exiting the grid. (e) 3D stellar velocity dispersion. (f) Half-mass radius of the entire star cluster. (g) Total mass (dotted line), mass of stars (dashed line) and gas (solid line) on the grid. (h) Virial parameter of stars (dashed line) and gas (solid line), where αv = 0.5 is the equilibrium value. (i) Fraction of mass bound for stars (dashed line) and gas (solid line).

Appendix B: Stellar modifications

B.1. Low-mass star agglomeration

Upwards of 10⁶ stars can be expected to form from a 10⁶ M_⊙ cloud with a peak number density of $n\approx 1000{\mathrm{cm}}^{-3}$. Even with the best modern N-body codes, evolving this many single stars and higher order stellar systems in such a dense stellar environment with a gravity bridge from each star to the gas in a separate code is immensely computationally taxing. To reduce the strain on the N-body portion of the calculations, we chose to agglomerate all stars under a given mass into gravitational super-star particles of equivalent mass to their sum. We refer to this mass cutoff as the agglomerate mass.

When a sink progresses through the list of stellar masses it will form, stars with masses under M_agg are put aside until the sum of their masses is above M_agg. Then a star particle is formed with the summed mass. Figure B.1 shows the reduction in number of stars formed in a cloud for a given agglomerate mass. For our choice of 4 M_⊙, we only had 10% of the stars undergoing gravitational interactions compared to the case with no agglomeration. This reduced our N-body execution time by a factor of somewhere between the 10log10 expected for the tree and 10³ expected for the direct N portion of the PETAR algorithm. We note that the feedback from these low-mass stars is shown in Appendix B.2 to be negligible compared to that of the higher-mass stars, and in any case the current TORCH version does not model feedback from stars $<8{\mathrm{M}}_{\odot }$ as we neglect jets, while the ionizing radiation from such low-mass stars is negligible. In this study, we further limited feedback to only come from stars $\ge 20{\mathrm{M}}_{\odot }$ as we discuss in the next subsection. The primary missing contribution from low-mass stars physically is their mutual gravitational interactions, which could potentially lead to the ejection of some fraction of them. However, the dynamics driven by those low-mass stars is also expected to be negligible in comparison to the effect of gas and more massive stars in the cluster. TORCH simulations with no mass agglomeration were done by Cournoyer-Cloutier et al. (2023), and in analyzing the morphology of clusters they find that the gravitational effects of the gas dominate over any stellar dynamics effect for the overall evolution of the cluster while it remains embedded.

B.2. Feedback mass limit

We limited all forms of stellar feedback—winds, radiation, and SNe—to stars above 20M_⊙ instead of the value of 8M_⊙ (lower bound for SN explosions) usually adopted in TORCH. This is necessary to significantly reduce the number of rays on the grid, which greatly decreases the calculation time and memory overhead for the ray-tracing algorithm. We quantify the effects of excluding radiation and winds from stars with masses below 20M_⊙ by comparing the power output in the form of winds and radiation from all stars above 8M_⊙ and above 20M_⊙. We only allowed stars above 20 M_⊙ to explode as SNe. Our simulations ran for $\le 10\mathrm{Myr}$, which is roughly the main-sequence lifetime of a 20 M_⊙ star. Stars under 20 M_⊙ do not explode as SNe in the timeframe of our simulations, so excluding their SNe feedback makes no practical difference for this comparison.

The power as a function of mass in the form of EUV radiation, non-ionizing FUV radiation, and stellar wind is shown in Figure B.2. We calculate these powers by taking stars from 8 to 100M_⊙ in 1M_⊙ increments, evolving them in SEBA for 1 Myr, and summing the energy output of each feedback channel. From this figure we can see that the power output of stellar winds and UV radiation is several orders of magnitude higher for stars above 20M_⊙ than for stars closer to 8M_⊙. Although stars in the 8–20 M_⊙ mass range still output a considerable amount of FUV radiation, stars above 20 M_⊙ account for over 80% of the total radiation power.

Although the feedback power is much stronger for stars above 20M_⊙, stars with masses 8–20M_⊙ greatly outnumber them. To find the ratio of feedback power for stars below and above 20 M_⊙, we convolve the number of stars of each mass with the power output for each stellar mass (Fig. B.3). In the top left histogram, we show the ratio of stars with mass 8–20 M_⊙ to stars with mass 20–100 M_⊙ in all three simulations, sampled at their respective initial free-fall times. All three runs have more stars in the lower-mass bin. In the top right plot, we show the ratio of total stellar feedback power P_FB (excluding SNe) for the stars in the two mass bins considered. We can see that although the lower-mass stars outnumber the higher-mass stars, the higher-mass stars still account for > 80% of the total stellar feedback energy. This shows that only including feedback from stars above 20 M_⊙ still retains almost all of the feedback energy produced after the formation of all three star clusters.

The bottom panel in Figure B.3 shows the feedback power per mass bin for each separate feedback process. For the EUV radiation and wind feedback, the low-mass stars contribute practically nothing to the feedback energy in comparison to the high-mass stars. The FUV feedback of low-mass stars is not negligible, but is still well below 20% of the total FUV feedback energy from all stars.

B.3. Mass-loading stellar winds

In TORCH, the stellar wind feedback implementation is inspired by Simpson et al. (2015), using a method of momentum injection, the details of which can be found in Wall et al. (2020). The energy of the cells within the wind injection radius of the star is increased based on the mechanical luminosity of the wind L_w = (1/2)̇Mv_w², where ̇M is the stellar mass loss rate (Vink et al. 2000) and v_w is the terminal wind velocity (Kudritzki & Puls 2000). The wind injection radius is set by comparing the cell width Δx to the wind termination shock radius (Weaver et al. 1977)

$$\begin{array}{c}\hfill R_{\mathrm{w}}=0.74{\left(\frac{M}{{\rho }_0}\right)}^{3/10}{v}_{\mathrm{w}}^{1/10}{t}_{\mathrm{w}}^{2/5},\end{array}$$(B.1)

where ρ₀ is the background density and t_w is the age of the wind-blowing star at the given time step. If R_w < Δx the injection radius is set to Δx, otherwise it is set to a maximum value of $6\sqrt{3}\mathrm{\Delta }x$, at which we have found that spherical winds are well resolved. Momentum and energy are conserved when injecting stellar winds.

Within a stellar wind bubble, in dense clumpy regions of star formation such as the ones in our simulations, material will be swept up into the flow of the hot bubble by mass loading processes such as photoevaporation and hydrodynamic ablation (Dyson & Hartquist 1992; Hartquist & Dyson 1996; Pittard et al. 2001; Lancaster et al. 2021). With enough mass loading, the density increase will result in much more efficient cooling and create momentum-driven rather than energy-driven bubbles. The amount of mass-loading in the case of hydrodynamic ablation depends on the prevalence of dense clumps within the wind region as well as the Mach number M of the flow around the clump. With a supersonic flow, the mass-loading rate saturates. With a subsonic flow, the mass-loading rate is proportional to M^4/3 (Smith et al. 1984; Hartquist et al. 1986). Accounting for mass loading in stellar wind models has been shown to successfully reproduce the kinematic properties of the observed stellar wind bubble of the Wolf-Rayet star RCW 58 (Arthur et al. 1993, 1996; Arthur 2007).

Simply injecting winds at v_w does not account for these mass-loading processes and results in unphysically hot bubbles. Therefore, we chose a lower temperature target for our bubbles and lowered the wind velocity v_w such that the final temperature of the wind bubble is the correct one. We conserved momentum and energy when injecting stellar winds, so while lowering the wind velocity, we also infused correspondingly more mass into the bubble than the stellar mass loss calculated. This mass is not taken off the grid elsewhere, meaning mass was not entirely conserved. At 1.25 t_ff, the total amount of mass that has been injected due to mass-loading as a fraction of initial cloud mass is 0.04, 0.02, and 0.01 for M4, M5, and M6.

Observed circumstellar bubbles cooled by suspected mass loading have been seen with temperatures as low as $T_{\mathrm{b}}\approx 1.1\times {10}^6\mathrm{K}$ in the S308 bubble (Chu et al. 2003). The spectra of the NGC 6888 bubble indicates a dominant component almost as cool, with $T_{\mathrm{b}}=1.5\times {10}^6\mathrm{K}$ (Bochkarev 1988; Wrigge et al. 1994; Wrigge 1999).

In the simulations presented here, we heavily mass loaded the stellar winds to achieve a lower than observed bubble temperature of T_b = 3 × 10⁵ K. This temperature is at the peak of the cooling curve, so the shocked wind rapidly cools, resulting in smaller, cooler, momentum-driven bubbles instead of hot bubbles filled with 10⁶ K gas. We chose to do so because the high sound speeds in hot wind bubbles lower the Courant time step significantly, making the computation impractical. Since we do not follow X-rays through ray-tracing, having cooler bubbles is adequate. Bubbles at this temperature also do not affect the ionization of the surrounding gas. The primary action of wind feedback during cluster formation is to clear out dense regions of gas so that radiatively ionized H II regions can expand. The only hot gas (≥10⁶ K) on the grid comes from SNe. Capping the temperature of gas on the grid at 3 × 10⁵ K until SNe occur significantly speeds up the simulations.

Fig. B.1. Fraction of the number of stars formed with agglomeration of stars below the mass on the x-axis over the total number of stars sampled by the IMF. We used $M_{\mathrm{agg}}=4{\mathrm{M}}_{\odot }$, which means the number of stars in the grid is 10% of the number that are formed by our IMF.

Fig. B.2. Power of stellar feedback in the form of winds and FUV and EUV radiation for different stellar masses. Left of the vertical line shows the amount of feedback power lost per star by excluding feedback from stars below 20 M⊙.

Fig. B.3. Power output of stellar feedback modes for each stellar mass regime. (Top left): Histogram showing the fractional stellar population of the three runs at one free-fall time, split into the mass regimes of 8–20 and 20–100 M⊙. (Top right): Fraction of feedback power in each mass regime. (Bottom): Histograms showing the fraction of feedback power for FUV, EUV, and winds in each mass regime. Although there are more lower-mass stars, the feedback produced by them is less than 20% of the total feedback energy for all stars.

B.4. Effect on star formation efficiency

Limiting the temperature of stellar winds and only modelling feedback for stars above 20 M_⊙ could potentially lead to un-physical runaway star formation. To test this, we re-ran the M6 model at early times to see if these two approximations are the cause for the extremely high SFE of 85%. For the first new M6 run we raised the wind temperature from $3\times {10}^5\mathrm{K}$ to $5\times {10}^6\mathrm{K}$. For the second test, we both raised the wind temperature and modelled feedback for all stars above 8 M_⊙. The SFE over time for the fiducial M6 run with our standard approximations and the new M6 models are shown in Figure B.4.

The two runs without the approximations that reduce the strength of the stellar feedback have similar SFEs as the M6 model with the aforementioned approximations. This validates our approximations and supports our argument that the high SFE in model M6 is not an artifact of under-estimating the strength of stellar feedback.

Fig. B.4. Star formation efficiency over time for the fiducial M6 cloud with Tw = 300,000 K and $M_{\mathrm{feedback}}\ge 20{\mathrm{M}}_{\odot }$ (cool winds), the M6 cloud with Tw = 5,000,000 K and $M_{\mathrm{feedback}}\ge 20{\mathrm{M}}_{\odot }$ (hot winds), and the M6 cloud with Tw = 5,000,000 K and $M_{\mathrm{feedback}}\ge 8{\mathrm{M}}_{\odot }$ (hot winds, $M_{\mathrm{fb}}=8{\mathrm{M}}_{\odot }$ ).

All Tables

All Figures

Fig. 1. Slice plots of the three simulations in the x–y plane over time. The plane of the slices for a given cloud is the center of stellar mass in the final snapshot. Stellar positions are shown by white dots. The free-fall times tff are given in Table 1. The number of stars shows the amount of star particles in the domain, not the number sampled from the IMF. Due to our agglomeration of low-mass star particles, the number of stars sampled from the IMF is ∼10× greater.

In the text

	Fig. 2. Still of the interactive plot of the embedded M4, M5, and M6 clusters (left to right) at 1 tff. The interactive plot file is available for download from the repository.
In the text

	Fig. 3. Wall-clock time for an evolution step for PETAR and PH4+MULTIPLES given the number of stars.
In the text

Fig. 4. Global properties of the clusters and gas over time for models M4 (orange), M5 (maroon), and M6 (blue-violet) in units of free-fall time tff of the initial cloud (given in Table 1). From top left to bottom right: (a) SFR, where the transparent lines show the SFR at each star formation event, and the solid lines give the SFR smoothed using a Gaussian filter with σ = 0.005tff. (b) Instantaneous and integrated SFEs of the clouds, where ϵinst = M⋆/(Mgas + Msink + M⋆) and ϵint = M⋆/Mcloud = ϵ⋆. (c) Most massive star formed. (d) Number of formed stars. Dashed line: actual number of stars that would form from sampling the IMF given the amount of gas mass collected for star formation by sink particles. Solid line: number of stars followed in TORCH after the sampled stellar population below 4 M⊙ has been agglomerated. Dotted line: number of stars above 20 M⊙ on the grid that are generating feedback. The number of stars can drop due to SN, mass loss, or exiting the grid. (e) 3D stellar velocity dispersion. (f) Half-mass radius of the entire star cluster. (g) Total mass (dotted line), mass of stars (dashed line) and gas (solid line) on the grid. (h) Virial parameter of stars (dashed line) and gas (solid line), where αv = 0.5 is the equilibrium value. (i) Fraction of mass bound for stars (dashed line) and gas (solid line).

In the text

Fig. 5. Size-line width coefficient versus surface density for the Firecracker cloud (Finn et al. 2019) and the M4, M5, and M6 clouds. The four red diamonds are observations, done using apertures (1–4) with radii R = 6.4, 15, 26, and 37 pc. The other points are our simulations, with connected points corresponding to apertures (1–4) of radii R = 5, 10, 15 and 20 pc. The colors indicate how much stellar mass has been produced (with empty points indicating no star formation). The lines correspond to virial equilibrium and free-fall as labeled (see Fig. 2 of Ballesteros-Paredes et al. 2011).

In the text

	Fig. 6. WARPFIELD evolution of shell radius versus time with different initial SFEs with the same parameters described in M6. In all cases, stellar feedback is inefficient in dispersing the surrounding dense cloud, and the shell eventually undergoes re-collapse.
In the text

Fig. A.1. Global properties of the clusters and gas over time for models M4 (orange), M5 (maroon), and M6 (blue-violet) for comparison to Figure 4, where units of free-fall time (see Table 1) are used. From top left to bottom right: (a) SFR, where the transparent lines show the SFR at each star formation event, and the solid lines give the SFR smoothed using a Gaussian filter with σ = 0.005tff. (b) Instantaneous and integrated SFEs of the clouds, where ϵinst = M⋆/(Mgas + Msink + M⋆) and ϵint = M⋆/Mcloud = ϵ⋆. (c) Most massive star formed. (d) Number of formed stars. Dashed line: actual number of stars that would form from sampling the IMF given the amount of gas mass collected for star formation by sink particles. Solid line: number of stars followed in TORCH after the sampled stellar population below 4 M⊙ has been agglomerated. Dotted line: number of stars above 20 M⊙ on the grid that are generating feedback. The number of stars can drop due to SNe, mass loss, or exiting the grid. (e) 3D stellar velocity dispersion. (f) Half-mass radius of the entire star cluster. (g) Total mass (dotted line), mass of stars (dashed line) and gas (solid line) on the grid. (h) Virial parameter of stars (dashed line) and gas (solid line), where αv = 0.5 is the equilibrium value. (i) Fraction of mass bound for stars (dashed line) and gas (solid line).

In the text

	Fig. B.1. Fraction of the number of stars formed with agglomeration of stars below the mass on the x-axis over the total number of stars sampled by the IMF. We used $M_{\mathrm{agg}}=4{\mathrm{M}}_{\odot }$, which means the number of stars in the grid is 10% of the number that are formed by our IMF.
In the text

	Fig. B.2. Power of stellar feedback in the form of winds and FUV and EUV radiation for different stellar masses. Left of the vertical line shows the amount of feedback power lost per star by excluding feedback from stars below 20 M⊙.
In the text

Fig. B.3. Power output of stellar feedback modes for each stellar mass regime. (Top left): Histogram showing the fractional stellar population of the three runs at one free-fall time, split into the mass regimes of 8–20 and 20–100 M⊙. (Top right): Fraction of feedback power in each mass regime. (Bottom): Histograms showing the fraction of feedback power for FUV, EUV, and winds in each mass regime. Although there are more lower-mass stars, the feedback produced by them is less than 20% of the total feedback energy for all stars.

In the text

Fig. B.4. Star formation efficiency over time for the fiducial M6 cloud with Tw = 300,000 K and $M_{\mathrm{feedback}}\ge 20{\mathrm{M}}_{\odot }$ (cool winds), the M6 cloud with Tw = 5,000,000 K and $M_{\mathrm{feedback}}\ge 20{\mathrm{M}}_{\odot }$ (hot winds), and the M6 cloud with Tw = 5,000,000 K and $M_{\mathrm{feedback}}\ge 8{\mathrm{M}}_{\odot }$ (hot winds, $M_{\mathrm{fb}}=8{\mathrm{M}}_{\odot }$ ).

In the text