Early evolution and three-dimensional structure of embedded star clusters

ABSTRACT

We perform simulations of star cluster formation to investigate the morphological evolution of embedded star clusters in the earliest stages of their evolution. We conduct our simulations with Torch, which uses the Amuse framework to couple state-of-the-art stellar dynamics to star formation, radiation, stellar winds, and hydrodynamics in Flash. We simulate a suite of 10⁴ M_⊙ clouds at 0.0683 pc resolution for ∼2 Myr after the onset of star formation, with virial parameters α = 0.8, 2.0, 4.0 and different random samplings of the stellar initial mass function and prescriptions for primordial binaries. Our simulations result in a population of embedded clusters with realistic morphologies (sizes, densities, and ellipticities) that reproduce the known trend of clouds with higher initial α having lower star formation efficiencies. Our key results are as follows: (1) Cluster mass growth is not monotonic, and clusters can lose up to half of their mass while they are embedded. (2) Cluster morphology is not correlated with cluster mass and changes over ∼0.01 Myr time-scales. (3) The morphology of an embedded cluster is not indicative of its long-term evolution but only of its recent history: radius and ellipticity increase sharply when a cluster accretes stars. (4) The dynamical evolution of very young embedded clusters with masses ≲1000 M_⊙ is dominated by the overall gravitational potential of the star-forming region rather than by internal dynamical processes such as two- or few-body relaxation.

1 INTRODUCTION

Most stars form within embedded clusters (Lada & Lada 2003; Portegies Zwart, McMillan & Gieles 2010). They remain shrouded in their natal gas for a few megayears after the onset of star formation (see e.g. Kim et al. 2022, for recent observations), while the cloud is still actively star-forming. Although most stars do not remain in bound star clusters for their whole lives, their formation and early evolution is shaped by the dense stellar environment in which they are born, which is in turn shaped by the interplay between gravity, turbulence, and stellar feedback. On smaller scales, stars also do not form in isolation: most stars form in multiple stellar systems (Offner et al. 2022, and references therein), most often in binaries. Binaries are known to be dynamically important for cluster long-term evolution (Heggie 1975; Hills 1975). Recent simulations by Torniamenti et al. (2021) further suggest that the presence of binaries impacts a cluster’s structure over time-scales of a few megayears after it has become free of gas. Despite their ubiquity, binaries in embedded clusters are seldom modelled numerically due the range of physical processes involved and the high numerical cost of modelling concurrently stellar dynamics on the scale of binaries and feedback processes impacting the gas in the embedded cluster.

Simulations of star cluster formation show that star clusters assemble through the merging of smaller embedded clusters over a few megayears (e.g. Fujii, Saitoh & Portegies Zwart 2012; Vázquez-Semadeni, González-Samaniego & Colín 2017; Grudić et al. 2018; Howard, Pudritz & Harris 2018; Chen, Li & Vogelsberger 2021). Karam & Sills (2022) have further shown that those mergers have an important impact on the boundedness of the stars and gas in the resultant cluster: some head-on collisions between clusters do not result in a single bound cluster, while there is mass loss and an increase in radius even in the successful mergers. The simulations conducted by Karam & Sills (2022) however do not account for the formation of new stars during cluster assembly. Recent work by Dobbs et al. (2022), which relies on star particles representing low-mass stellar populations or massive stars to model clusters, also reveals a more complex picture: clusters cannot only merge, but also split. They also trace the mass and size of their clusters throughout their simulations and find no clear correlation between mass and size. They however assume a spherical shape when measuring the size of their clusters, which is not the case for observed embedded clusters (e.g. Kuhn et al. 2014). Furthermore, neither of these recent suites of simulations include binaries, which are expected to influence stellar dynamics on the cluster scale, at least once the cluster becomes free of gas.

The virial parameter α of the star-forming cloud of gas, which describes the balance between the effects of self-gravity and turbulent support of the gas, is also important for cluster formation and evolution. For a spherical cloud, the virial parameter is defined as

where T is the kinetic energy of the cloud, U is its gravitational potential energy, σ is its velocity dispersion, R is its radius, M is its mass, and G is the gravitational constant (Bertoldi & McKee 1992). Thus clouds with smaller α are more strongly bound, and α = 1 corresponds to virial equilibrium. Observed clouds in galaxies cover a large range of virial parameters, from α ≲ 0.1 to α ≳ 100 (Kauffmann, Pillai & Goldsmith 2013). A cloud’s virial parameter systematically affects its star formation efficiency (SFE) and cluster formation efficiency (CFE; Kruijssen 2012; Howard, Pudritz & Harris 2016), with regions with higher α generally having lower SFE and CFE.

In this work, we use numerical simulations to investigate the effects of stellar dynamics and cloud-scale hydrodynamics on the structure and evolution of embedded star clusters. To test the relative importance of stellar dynamics, we explore the impact of forming (or not forming) primordial binaries with different underlying populations, as binaries are known to play an important role in setting cluster structure in systems dominated by stellar dynamics (e.g. Heggie 1975; Fujii & Portegies Zwart 2011; Torniamenti et al. 2021). To test the relative importance of cloud-scale hydrodynamics, we vary the cloud’s initial virial parameter α, which is known to have a strong effect on the CFE (e.g. Howard et al. 2016). We want to determine (1) whether cloud-scale hydrodynamics or stellar dynamics have the strongest impact on cluster structure (mass, size, and shape) and cluster formation efficiency and (2) how cluster structure evolves during the earliest stages of formation.

In Section 2, we describe our numerical framework and our simulations. In Section 3, we follow the evolution of the bulk properties of the stars in the simulation domain, we investigate the instantaneous properties of the clusters as a population, and we examine the assembly history of individual clusters; Section 3.3 contains the key results of the paper. In Section 4, we discuss the broader implications of our findings. We summarize our results in Section 5.

2 METHODS

2.1 Numerical framework

We use Torch (Wall et al. 2019, 2020; Cournoyer-Cloutier et al. 2021), which relies on the Amuse framework (Portegies Zwart et al. 2009; Pelupessy et al. 2013; Portegies Zwart et al. 2013; Portegies Zwart & McMillan 2019) to couple hydrodynamics to stellar dynamics, star and binary formation via sink particles, stellar evolution, and stellar feedback in the form of winds and radiation. Torch is optimized to investigate the effects of stellar and binary dynamics in young, gas-rich clusters, in particular stable multiple systems and dynamical short-range encounters between stars. We model the self-gravitating gas with the adaptive mesh refinement code Flash (Fryxell et al. 2000). We use simultaneously two types of refinement criteria for our adaptive grid. We first require that the Jeans length be resolved by at least four resolution elements in order to avoid numerical fragmentation (Truelove et al. 1997; Federrath et al. 2010). To improve stability, we also refine around sharp changes in pressure, temperature, total energy, and internal energy (i.e. at shocks and contact discontinuities, Löhner 1987; MacNeice et al. 2000). Although Flash can evolve magnetic fields, we do not include them in our simulations due to their high computational cost. We treat gas dynamics with a Harten-Lax-van Leer Riemann solver (Miyoshi & Kusano 2005) and an unsplit (magneto)-hydrodynamics solver (Lee 2013) with third-order piecewise parabolic method reconstruction (Colella & Woodward 1984). We handle the gas self-gravity with a multigrid solver (Ricker 2008) while we handle the gravitational attraction of the gas on the stars and vice-versa with a leapfrog scheme (Wall et al. 2019, based on Fujii et al. 2007).

On the stellar dynamics side, we handle long-range stellar dynamics with the direct N-body code Ph4 (McMillan et al. 2012), which uses a fourth-order Hermite predictor–corrector scheme (Makino & Aarseth 1992). For stable binary (and higher order) systems, resonant encounters, and scattering, we use the Amuse module Multiples (Portegies Zwart & McMillan 2019), which itself uses the codes smallN (Hut, Makino & McMillan 1995; McMillan & Hut 1996) and kepler (originally developed as part of Starlab, Portegies Zwart et al. 1999; Hut et al. 2010).

Star formation takes place within sink particles that are treated as star factories. The details of the sink implementation are presented in Wall et al. (2019) for single star formation and Cournoyer-Cloutier et al. (2021) for binary formation. Briefly, a sink particle is formed when the local gas density and convergence criteria outlined in Federrath et al. (2010) are satisfied. Once formed, it samples an initial mass function (Kroupa 2001) between 0.08 M_⊙ and 150 M_⊙ to generate a list of stars to be formed, using a Poisson sampling method first tested by Sormani et al. (2017) and implemented in Torch by Wall et al. (2019). Each star in the list is formed when the sink has accreted sufficient mass, in order to ensure quasi-local mass conservation. The sink must also sit in cold (<100 K) gas to form stars. Stars are formed with a gas-to-star conversion efficiency of 100  per cent. The additive properties of the Poisson distribution ensure that the sampling for the full simulation domain reproduces the IMF, despite possible stochastic variations within individual clusters. The decoupling allows the stars to be handled by the N-body solver Ph4, which is fourth-order accurate, instead of the second-order leapfrog scheme used for sink particles. Although the formation of individual stars is unresolved in our simulations, stellar dynamics are followed self-consistently after star formation.

Stars with masses above 7 M_⊙ inject radiative and momentum feedback on the grid. The details of the feedback implementation are presented in Wall et al. (2020). The far ultraviolet (between 5.6 eV and 13.6 eV) and ionizing (above 13.6 eV) radiative feedback is implemented within Flash as a modified version of the adaptive ray-tracing module Fervent (Baczynski, Glover & Klessen 2015). The total and average photon energy are calculated for each star from the surface temperature and mass obtained from stellar evolution, which is performed with SeBa (Portegies Zwart & Verbunt 1996). All radiative feedback heats the gas. Massive stars further provide feedback in the form of momentum-driven winds with mass loss rates based on Vink, de Koter & Lamers (2000). Radiative cooling of the gas from atomic and molecular lines and dust is included (Wall et al. 2019).

2.2 Simulations and star formation prescriptions

We conduct a total of 12 simulations, summarized in Table 1. All simulations are initialized from a spherical, turbulent cloud of neutral dense gas with a mass of 10⁴ M_⊙ and a radius of 7 pc in a cubic box of side 17.5 pc, following the model used in Cournoyer-Cloutier et al. (2021). The mean gas surface density is 50 M_⊙ pc⁻². Those values are consistent with a typical cloud in the Solar neighbourhood (Chen et al. 2020). The cloud follows a Gaussian density profile with a central density 8.75 × 10⁻²² g cm⁻³ and temperature 20.64 K, and sits in a warm neutral medium with density 2.18 × 10⁻²⁴ g cm⁻³ and temperature 6.11 × 10³ K. These values were chosen to ensure pressure and thermal equilibrium between the cloud and surrounding medium. The free-fall time for the cloud is 1.45 Myr. The gas follows an adiabatic equation of state with γ = 5/3, although radiative cooling maintains the dense neutral gas almost isothermal. We adopt the same gas spatial resolution of 0.0683 pc at the maximum refinement level and density threshold for the formation of sink particles of 3.82 × 10⁻²¹ g cm⁻³ as used in Cournoyer-Cloutier et al. (2021).

We consider four different prescriptions for binaries (described in Appendix  A), in addition to models without primordial binaries. Our models with primordial binaries span a range of mass-dependant binary fractions, mass ratios, and orbital periods. We stress that the details of those prescriptions are not the focus of this paper—rather, we test diverse models for primordial binaries to fully explore the impact that a change in stellar dynamics has on embedded cluster structure and evolution. Binaries can also form dynamically, and the properties of primordial and dynamically formed binaries will be modified by dynamics over the course of our simulations. We refer the interested reader to Cournoyer-Cloutier et al. (2021) for a detailed discussion of the effects of dynamical interactions on the initial population of binaries.

Eight of our 12 simulations (with names starting with B-P and S-R) are initialized with a virial parameter α = 2T/|U| = 0.8. The gas is initially gravitationally bound and its collapse is expected to result in abundant star formation. The gas initial conditions, including the random turbulent field, are identical for those 8 simulations. We also perform simulations with larger virial parameters, of α = 2.0 and α = 4.0. We perform pairs of simulations with our fiducial prescription for binaries and our single stars only prescription (both with the default random seed) for both these models, and label them B-V2, S-V2, B-V4 and S-V4. Those initial conditions are set up by scaling up the gas velocities in each cell of the initial conditions for the α = 0.8 runs. We therefore increase α but conserve the direction of motion of the gas in each cell.

Beyond the binary prescriptions and virial parameters, we also vary the random seed used to sample the initial mass function and to form binaries, which sets the masses of the stars and the order in which they form. Our simulations labelled with the random seed default all use the same random seed; the other simulations all use different random seeds. We use different random seeds to ensure our general conclusions are not affected by the stochastic formation of massive stars. We have shown in Lewis et al. (2023) that early-forming massive stars can promote the formation of smaller, isolated clusters, and prevent the formation of massive clusters. By using different random samplings of the IMF, we can verify that our conclusions are not drawn from a single, extreme case, in which the formation times and masses of the massive stars providing radiative and mechanical feedback would be atypical.

2.3 Cluster identification

Most of our simulations have reached 2 Myr after the onset of star formation, and snapshots are written every 0.01 Myr. We inspect all snapshots in our simulations for clusters, which we identify from a combination of spatial clustering and boundedness. We initially select clusters with DBSCAN (Ester et al. 1996; Pedregosa et al. 2011) based on the positions of the stars. We require each cluster star to have five neighbours (following Sander et al. 1998, for three-dimensional data), which are other stars within a user-determined distance. For our analysis, we fix this distance to the sink accretion radius, 0.17 pc. Following our initial identification of the clusters, we perform a boundedness check on the stars with respect to their associated cluster. For each star, we calculate the gravitational potential energy from the local gas gravitational potential (including the sink particles) and the potential from the cluster’s stars. We also calculate the stars’ kinetic energy in the cluster’s centre of mass frame. We remove stars with positive total energy (i.e. unbound stars) from the cluster. After this boundedness check, clusters that have at least 100 members are saved for subsequent analysis. An example of the clusters satisfying our clustering, boundedness and minimum membership criteria in a given snapshot is shown in the left-hand panel of Fig. 1.

2.4 Cluster structure

Once clusters are identified, it is useful to describe their size, which in turn requires us to measure their shape. Observational studies have used respectively ellipses (e.g. Kuhn et al. 2014; Zhai et al. 2017) and ellipsoids (e.g. Pang et al. 2021) to describe the 2D and 3D shapes of embedded or open clusters. We similarly use 3D ellipsoids to describe the shape of some inner fraction of the stellar distribution in an individual cluster—here, we use 50  per cent and 90  per cent mass ellipsoids, as proxies for the 50  per cent and 90  per cent Lagrangian radii. We use the fact that any distribution of points can be described by an inertial ellipsoid that shares its rotational properties about its principal axes (see e.g. Goldstein, Poole & Safko 2001). This technique has been used previously in astrophysics to describe the 3D shape (or projected 2D shape) of dark matter halos in cosmological simulations (see e.g. Velliscig et al. 2015; Thob et al. 2019; Hill et al. 2021; Reina-Campos et al. 2022). We show the 90  per cent mass ellipsoids for the last snapshot of S-R0 (the same example as for the clustering plot) in the right-hand panel of Fig. 1. We present the details of our fitting routine in Appendix  B. An example of the 50  per cent and 90  per cent mass ellipsoids for an individual cluster identified in our simulations is also provided in Fig. B1, and the spherical half-mass and 90  per cent Lagragian radii are provided for comparison.

We use our fitted ellipsoids to define a proxy for the radius, to compare our clusters to established mass–radius relations. We do so by taking the geometric mean of the semi-major, intermediate, and semi-minor axes a, b and c to define a characteristic radius

which is similar to what is done in observational studies (e.g. Kuhn et al. 2014) to define sizes for elliptical 2D clusters. We can define such a radius for any enclosed mass fraction, and therefore for any Lagrangian radius. To quantify how non-spherical a cluster is, we define an ellipticity (see e.g. Kuhn et al. 2014)

that depends on the ratio between the semi-minor and semi-major axes. A spherical cluster has an ellipticity ϵ = 0 while a very elongated cluster has an ellipticity ϵ → 1. With equations (2) and (3), we characterize the size and shape of individual clusters at each snapshot in our simulations.

2.5 Cluster history

We follow the evolution of individual clusters throughout the simulations. For each cluster identified in the last snapshot of a simulation, we trace back its main progenitor in earlier snapshots by identifying the cluster sharing the largest fraction of its stellar mass in the previous snapshot. We also look for clusters that are present at earlier times but are no longer present in the last snapshot. We allow for clusters to be missing in some checkpoints (e.g. if a cluster with 100 bound members loses one star, then forms one or more later on) but require a cluster to be present over at least 0.1 Myr to trace its history. In practice, this means that a cluster can be used in our analysis of cluster populations (e.g. Sections 3.1 and 3.2) without being used in our analysis of cluster histories (Section 3.3) if it survives for less than 0.1 Myr.

We use our results on cluster histories in three main ways. First, we track the evolution of the mass, size, and shape of individual clusters to investigate the presence of evolutionary trends. Second, we investigate the relative contributions of cluster mergers, the accretion of unclustered stars, and new star formation to the build-up of our clusters during the first ∼2 Myr after the onset of star formation. Third, we evaluate what proportion of cluster stellar mass is lost over the same time. Those relative contributions are not final, as the clusters are still growing in mass at the end of the simulations. They however give us a picture of the variations in cluster history during the early stages of embedded cluster evolution.

We rely on the tags given to star particles to follow the assembly of individual clusters. Between two subsequent snapshots in which a cluster is identified, we identify all new star particles and all star particles that left the cluster. For new star particles, we verify whether they were present in the previous snapshot (either in another cluster or as unclustered stars). If they were not present in the previous snapshot, we consider them to be newly formed stars, and treat them as having formed in the cluster. If they were present, we consider them as accreted stars. Star particles that have left the cluster are recorded as lost stars. For accreted and lost stars, we ensure that there is no double-counting, which could occur for example if a merger is unsuccessful or if a cluster splits. To evaluate cluster assembly, retained stars are therefore treated as formed in the cluster, accreted, or lost, if they respectively fulfill the following criteria:

• Stars are considered formed in the cluster if they were not present (as a clustered or unclustered star) in the snapshot before they are identified as a cluster member, and are present in the cluster in the last snapshot in the simulation. Some of the stars complying with these criteria may have been lost and then re-accreted.

• Stars are considered accreted if they were present in another cluster or as an unclustered star in the snapshot before they are identified as a cluster member, and are present in the cluster in the last snapshot in the simulation. Such stars may also have been lost and then re-accreted.

• Stars are considered lost if they were present in the cluster at any earlier snapshot, and are not present in the cluster in the last snapshot in the simulation. Such stars may also have been lost, re-accreted, and then lost again.

The stars that were cluster members when the cluster was first identified are treated separately to avoid artificially driving up the formed or accreted fractions in low mass clusters. We record the composition of the cluster at the end of our simulations (i.e. the mass in initially present, formed, and accreted stars), as well as the mass lost throughout the history of the cluster.

3 RESULTS

We structure our results in three subsections, corresponding to three different approaches to analysing our simulations. In Section 3.1, we summarize the evolution of the full simulation by tracking properties such as the star formation rate (SFR) and the clustered stellar mass. In Section 3.2, we track the mass, size, and morphology of the identified clusters as a population, and compare them to observations of Galactic clusters. In Section 3.3, we explicitly follow the evolution of the clusters throughout the simulations by tracking how they assemble their mass and how their morphologies change.

3.1 Overview: properties of the full simulation domain

We first look at the global properties of the simulations. The starkest differences are between simulations with different initial virial parameters α. This is already obvious from the plots of the gas column density presented in Fig. 2. The three simulations shown in the figure have the same star formation model (single stars only, default random seed) but are initialized with virial parameters of respectively α = 0.8, 2.0, and 4.0. Some features in the gas (such as the inverted Y shape made by the densest gas) persist across the three plots, but the gas morphology is nonetheless obviously different in the three simulations. In particular, the gas is less centrally concentrated and closer to the edges of the domain in the simulations with larger virial parameters.

Those morphological differences naturally give rise to differences in the SFR. The SFR and integrated SFE (mass of all formed stars divided by the initial gas mass) for the different simulations are plotted against time since the onset of star formation (SF) in Fig. 3. We use a Gaussian filter with a kernel width of 0.1 Myr to smooth both the SFR and the SFE, in order to remove instantaneous peaks in the SFR caused by the formation of individual massive stars. By the time we stop the simulations, the SFR and the SFE are both about half an order of magnitude larger in our simulations with the fiducial α = 0.8 than in the simulations with α = 2.0, and more than an order of magnitude larger than in the simulations with α = 4.0. The different prescriptions for binary formation and the choice of random seed for star formation do not systematically affect the SFR or the SFE. They however cause scatter, which is smaller than the systematic effects associated with variations in α.

Simulations with higher virial parameters begin forming stars earlier than simulations with lower virial parameters. The first stars form respectively at t_SF = 1.12, 1.04, 0.87 Myr in simulations with α = 0.8, 2.0, 4.0. This is consistent with our expectations: turbulence both promotes star formation—in leading to an earlier onset of SF—and prevents it—in lowering the SFR (Ballesteros-Paredes et al. 2007) but its net effect is to decrease the SFR (Mac Low & Klessen 2004). We also note that the clusters (and the stars) in our simulations with α = 0.8 tend to form along a linear chain of width ∼ 1 pc (see Fig. 1). This is similar to the complexes of embedded clusters in DR 21, NGC 2264, NGC 1893, NGC 6334, and the Carina Nebula observed by Kuhn et al. (2014) in the MYStIX survey.

3.2 Average properties of individual clusters

We now turn our attention to the properties of embedded clusters identified in our simulations as a population. For this section, we use all clusters identified at all times in our simulations and measure their masses, sizes, and ellipticities.

In Fig. 4, we present a mass–radius plot for all individual clusters in our simulations, where the characteristic radius r₅₀ for the 50  per cent ellipsoid is calculated with equation (2) and the cluster mass is obtained from the sum of the masses of all stars identified as cluster members. The diagonal lines denote lines of constant mass density or surface density. We also show in Fig. 4 the characteristic radii of six deeply embedded clusters with median X-ray energy in the 0.5–8.0 keV band above 2.0 keV and at least 100 members from the MYStIX survey (Kuhn et al. 2014). Median X-ray energy is a proxy for extinction, and therefore anti-correlated with cluster age; the six selected clusters are expected to be the best match to our simulated cluster population in age and mass. We calculate the characteristic radii of the observed clusters from

where a and b are respectively the semi-major and semi-minor axes of the projected ellipses. We estimate the mass M of the observed clusters from their star counts, using

where N denotes the star count, and the slope is obtained by fitting the mass against the star count for our simulated clusters. Most identified embedded clusters have masses around 100 M_⊙, characteristic half-mass radii around 0.05 pc, and therefore densities (calculated within the characteristic half-mass radius) between 10⁴ and 10⁵ M_⊙ pc⁻³. The most massive clusters have densities around 10⁵ M_⊙ pc⁻³. This is approximately the same density as the Arches cluster (Serabyn, Shupe & Figer 1998), which has a similar age of ∼2 Myr but a mass of a few 10⁴ M_⊙, about two orders of magnitude larger than our clusters. We therefore conclude that our clusters have densities comparable to the upper limit of observed densities in young Galactic clusters.

Using the two-sample Kolmogorov–Smirnov (KS) test, we find no statistically significant difference in the masses, radii, or densities of the clusters identified in simulations with and without primordial binaries. This supports our earlier conclusion that there are no structural differences in clusters with and without primordial binaries over the time-scales spanned by our simulations. Similarly, the clusters formed in simulations with different virial parameters cover similar regions in mass–radius space.

We also find clusters with unphysically high densities (above 10¹⁰ M_⊙ pc⁻³) within r₅₀, which is often due to a single star accounting for ≳30 pre cent of the cluster’s mass. One star may contribute to up to ∼50 per cent of the cluster’s mass: an extreme example is a ∼ 90 M_⊙ star in a ∼ 200 M_⊙ cluster, in S-R2. This skews the mass density to much higher than that of observed clusters but the number density remains reasonable. The clusters discussed here are still actively forming, and we expect them to grow via the formation of new stars and mergers with other clusters before star formation halts; the massive star discussed above is therefore expected to become part of a more massive cluster or to be lost as a runaway star.

We plot cluster ellipticity against characteristic radius r₅₀ in Fig. 5. In the top panel, the ellipticity shown is that of the ellipsoid enclosing 50  per cent of a cluster’s mass, calculated with equation (3). We compare the radius-ellipticity distribution to that of the same six deeply embedded clusters with at least 100 members from the MYStIX survey (Kuhn et al. 2014). For the observed embedded clusters, we calculate the characteristic radius from equation (4). Given the apparent mismatch between the simulated and observed clusters, we explore projection effects. To complete this more robust comparison to observations, we calculate the size and ellipticity from 2D projections of the simulated clusters’ shapes and present them in the bottom panel. Each 3D ellipsoid is projected along a randomly selected axis, and the semi-major and semi-minor axes |$\tilde{a}$| and |$\tilde{b}$| of the projected ellipse are used to calculate the characteristic radius and ellipticity respectively from equations (4) and (3). When accounting for projection effects, we find that our simulated clusters have ellipticities similar to those of the deeply embedded objects in the MYStIX sample, although our simulated clusters tend to have smaller radii. Kuhn et al. (2014) however find that the sizes of embedded clusters in their sample are positively correlated with cluster age. For their subsample of very deeply embedded objects—which are the most comparable in age to our simulated clusters but not limited in star count—they find an average projected radius of 0.04 pc, which is in good agreement with our simulated clusters.

Our simulated embedded clusters have realistic masses, sizes, densities, and ellipticities: conclusions drawn from the study of their evolution can therefore inform our understanding of observed embedded clusters. We further note that there are no systematic differences in the structural properties of our simulated embedded clusters as a population regardless of the presence of primordial binaries or the choice of initial virial parameter for the star-forming cloud.

3.3 Time evolution of individual clusters

We now investigate the evolution of individual clusters throughout the simulations. We find no individual cluster satisfying our membership and boundedness criteria that survived more than 0.1 Myr and then merged with another cluster. We however find that some clusters acquire more than ∼100 M_⊙ due to accretion. Two processes contribute to this accretion budget without being registered as mergers. First, clusters satisfying our minimum membership and boundedness criteria may be accreted less than 0.1 Myr after they are first detected, and so before their histories are tracked. Second, groups of stars with fewer than 100 bound members—that are not recorded as clusters—may be accreted. We further find six examples of clusters splitting from an already-formed cluster, four of which survived for more than 0.1 Myr (the other two are identified because they are present in the last snapshot of the simulation).

In Fig. 6, we compare the relative impact of accretion and star formation on the assembly of our simulated embedded clusters. As examples, we show the four most massive clusters from the B-P2, B-P3, S-R2, and S-R3 simulations, which were run respectively to 2.0, 2.1, 2.0, and 2.2 Myr after the onset of star formation. In all cases except the second most massive cluster in S-R2, the formation of new stars within the cluster contributes more mass to the cluster than the accretion of already-formed stars. In that example, accreted mass contributes 53  per cent of the total final mass of the cluster. The most extreme example of splitting is shown in the third most massive cluster in S-R3: the cluster split from the most massive cluster in the simulation 0.1 Myr before the last snapshot, and had a mass of 327 M_⊙ just after splitting (see initial mass of the third cluster in the fourth row of Fig. 6).

We present in Fig. 7 an overview of the relative contributions of mass loss, accretion, and in-cluster star formation to the history of the embedded clusters in our simulations. The lost mass is calculated from the ratio of the mass lost by the cluster to the total mass acquired by the cluster over its history—i.e. the final mass plus the lost mass. The accreted (formed) fraction is calculated as the fraction of the final stellar mass of the cluster that was accreted (formed) after the cluster was first identified. The accreted and formed fractions for a cluster therefore do not add up to 100  per cent, as the stellar mass present in the cluster when it is first identified also contributes. The first violin plot for the lost fractions does not include clusters that split into two clusters surviving for more than 0.1 Myr.

There are no statistically significant differences in the final compositions of clusters with and without primordial binaries, as verified by a series of two-sample KS tests comparing the fractions of the stellar mass lost, accreted, and formed within the cluster for simulations with and without primordial binaries. We also find no statistically significant difference for clusters formed in simulations with α = 0.8, α = 2.0, and α = 4.0. We however find a rich variety of relative contributions from accretion and star formation, at all cluster masses. In other words, we find that there is no single dominant growth mechanism for clusters while they are still deeply embedded and actively star-forming, although generally stars formed in situ outnumber accreted stars in a given cluster.

We have shown in Section 3.2 that cluster radius and cluster mass are uncorrelated for our full population of simulated embedded clusters. We now investigate how the radii of individual clusters change as they grow in mass. In Fig. 8, we present as examples the evolution in radius–mass space of the most massive clusters in the B-P2, B-P3, S-R2, and S-R3 simulations. We find once again no correlation between radius and mass. We however note that radius can change by up to one order of magnitude without significant changes to the mass (see e.g. S-R3).

We also investigate the evolution of the mass, radius, and ellipticity of individual clusters as a function of the time since they were first identified. In Fig. 9, we plot these quantities against the time since cluster formation for all clusters identified in B-P2, B-P3, S-R2, and S-R3. The first important result we glean from this plot is that the mass growth of our simulated clusters is not always monotonic. A clear example is the most massive cluster in S-R3, which loses ≳300 M_⊙ within 0.01 Myr as it splits, as discussed above. In most cases, however, the mass of the most massive cluster tends to grow exponentially with time. In all cases except for S-R1, the most massive cluster is also the longest-lived cluster. In most cases, it is however not the one with the highest growth rate, which suggests that another cluster could become more massive at later times. Overall, simulations with and without primordial binaries follow the same general trends.

In Fig. 9, we also explore the time evolution of the characteristic radius r₅₀ of the clusters in our simulations. We find no correlation between the characteristic radius of a cluster and the time since it was formed. This is an important result as it suggests that the evolution of our embedded simulated clusters is not yet dominated by their internal dynamics, which should cause expansion (see e.g. Torniamenti et al. 2021, for recent simulations). We further highlight that considerable changes in cluster radius, of half an order of magnitude, occur on time-scales shorter than 0.01 Myr (i.e. between two consecutive snapshots). We also plot the ellipticity of the distribution of cluster stars enclosed within their characteristic radius r₅₀ against time since cluster formation. We once again find no correlation, and find that considerable changes can take place over ∼0.01 Myr. Changes in cluster size and shape, while the clusters are actively forming, are therefore driven by physical processes more complex than simply growth in mass or effects from stellar dynamics.

Rapid changes in morphology are driven by accretion and splitting events. We compare the timing of changes in cluster mass, characteristic radius, and ellipticity in Fig. 10, and find that they occur at the same times. In particular, we find that the times for the local minima and maxima in characteristic radius and ellipticity match. This is not due to a general correlation between size and shape, as shown in Fig. 5: rather, it indicates that clusters grow more elliptical and grow in size at the same time, when they are actively accreting an infalling group of stars. ‘Failed’ accretion events, or events followed by a splitting of the cluster, result in a rapid increase of the radius and ellipticity followed by a rapid decrease of the radius and ellipticity as the cluster returns to its original state (see left-hand panel of Fig. 10). When the accretion event is successful, the cluster’s radius and ellipticity also grow rapidly but decrease more smoothly after the event (see right-hand panel of Fig. 10).

Together, our results indicate that the early structure evolution of embedded clusters is driven by processes arising from the larger cluster-forming region—such as a burst of star formation due to an inflow of gas, or the accretion of a group of stars—rather than by their internal dynamics. In particular, we find that the composition of the clusters is being modified by the formation of new stars or the accretion of already formed stars over time-scales much shorter than the clusters’ relaxation time. We can obtain back-of-the-envelope estimates of the lower and upper limit on the relaxation times for the clusters in our simulations with

where t_cross is the crossing time,

and ρ is the average density inside a particle’s orbit (Binney & Tremaine 2008). Although both equations are only exact for spherical systems, they nonetheless provide us with a simple estimate of the time-scales relevant for our simulated embedded clusters. For the lower limit on the relaxation time, we assume a cluster with 100 members, the minimum possible in our analysis framework, and use a density of 10⁴ M_⊙ pc⁻³, which is towards the low end of our density values but still common. For the upper limit, we assume 1000 members, which is at the high end for our simulated embedded clusters, and use a density of 10⁵ M_⊙ pc⁻³. These give us estimated relaxation times of ∼0.32 and ∼0.68 Myr, which are longer than the time-scales over which the mass—and therefore the number of stars—of the embedded clusters change. Indeed, sudden accretion events, like the attempted merger shown in the left-hand panel of Fig. 10, may change a cluster’s characteristic radius r₅₀ by up to an order of magnitude, and its mass by a factor of ∼ 1.5, over ∼ 0.01 Myr. Those time-scales are more than 10 times shorter than the estimated relaxation times for the clusters. We conclude that the embedded clusters present in our simulations are not relaxed, despite their small sizes, due to how frequently they form or accrete new stars. Their dynamical evolution is therefore still driven by the overall gravitational potential of the simulation domain, dominated by the gas, rather than by two- or few-body encounters within the cluster.

4 DISCUSSION

We now discuss the results presented in Section 3. We compare the efficiency of star formation in our simulations to recent observations and simulations of cluster-forming regions. We also discuss the broader implications of our results for observational and computational studies of embedded clusters. We end by outlining areas for future work.

4.1 Star formation efficiency

We have explored a range of realistic initial virial parameters α, ranging from a low virial parameter typical of more massive clouds in which young massive clusters (YMCs) form to a moderately high virial parameter typical of the 10⁴ M_⊙ clouds in the solar neighbourhood. The general trend is for our simulations with higher α to have lower SFE and a smaller mass for the most massive cluster, in agreement with observations (e.g. Schruba, Kruijssen & Leroy 2019) and GMC-scale hydrodynamics simulations (e.g. Howard et al. 2016). The physical quantities obtained for the full simulated domain and for individual and average clusters are generally in agreement with observations of Galactic star-forming regions and young clusters. The integrated SFE after one free-fall time is about 1  per cent in our simulations with α = 4.0 while it is about 3  per cent in our simulations with α = 2.0. Observations suggest a SFE per free-fall time of about 1  per cent on pc-scale clouds, with scatter up to about 3  per cent (Krumholz, McKee & Bland-Hawthorn 2019, and references therein); this is consistent with the results from our simulations.

The integrated SFEs after one free-fall time in our simulations are also consistent with the results from the starforge simulations conducted by Guszejnov et al. (2022) with a cloud mass of 2 × 10⁴ M_⊙ and virial parameters α = 1.0, 2.0, 4.0. Those simulations include models for protostellar outflows, radiation, stellar winds, and supernovae. Our simulations with α = 0.8 have SFEs per free-fall time ≲10  per cent, which is similar to what they obtain in their simulations with α = 1.0. It is further possible to compare our simulations with α = 0.8 to the work conducted by Howard et al. (2016) using flash with radiative feedback. The SFR after one free-fall time is ∼10⁻³ M_⊙ yr⁻¹ for our 10⁴ M_⊙ clouds with α = 0.8. This is consistent—after scaling for cloud mass—with the SFRs obtained by Howard et al. (2016) for their clouds with α = 0.5 and α = 1.0, where the SFR after one free-fall time is ∼10⁻¹ M_⊙ yr⁻¹ for a 10⁶ M_⊙ cloud.

4.2 Implications for observations

There are two key takeaways from our simulations that can be applied to observed embedded clusters. First, the structure of embedded clusters—such as shape and size—can change considerably over time-scales as short as 0.01 Myr, due to new star formation or accretion. Those changes are not monotonic, do not follow a general trend, and are not driven by internal dynamics. This contrasts with studies of the early evolution of gas-free young star clusters. Torniamenti et al. (2021), for example, find that gas-free young clusters expand faster in the presence of primordial binaries. The evolution of the morphology of our embedded clusters, however, is driven primarily by the acquisition of news stars via star formation or accretion. Both processes are themselves driven by gas dynamics: star formation takes place within dense, converging flows of gas, while the accretion of already-formed stars is driven by the gravitational dynamics of the gas, which still accounts for |$\gtrsim ~70~{{\ \rm per\ cent}}$| of the mass within the simulation domain at the end of the runs. The simulations presented here focus on the deeply embedded stages of cluster formation. Tentative conclusions about the behaviour of our clusters up to and following gas expulsion can be reached by considering the results from the simulations that we presented in Lewis et al. (2023), albeit with some caveats: the simulations presented in Lewis et al. (2023) have a spatial resolution four times coarser than the present work (0.27 pc versus 0.0683 pc), consider only one virial parameter α selected to promote abundant star formation, and do not include binaries.

We stress therefore that the observed state of Galactic embedded clusters in their first stages of formation is instantaneous. We argue that no conclusions about the future evolution of a very young embedded cluster can be drawn from its current size or ellipticity: the cluster’s current state gives no information about whether the radius or ellipticity will increase or decrease in the future. Environment and recent changes in stellar mass play a role at least as important as internal processes such as two- or few-body encounters in setting embedded clusters’ dynamical states. We further note that more information about the stellar content and kinematics of very young embedded clusters –including information about binaries—would not allow us to predict their evolution better. We thus also predict that observations of stellar positions and velocities—that could be used to verify boundedness, investigate cluster expansion, and measure cluster shape—cannot be used to infer the presence of a significant number of binaries in embedded clusters.

Second, we find that clusters tend to have large ellipticities and large characteristic radii when they are accreting new stars. Examples are shown in Fig. 10. A large ellipticity for an embedded cluster with a large radius, that persists despite projection effects, could be clear observational evidence that the embedded cluster is currently accreting—or has recently accreted—new stars without requiring any stellar velocity data. Torniamenti et al. (2021), in their simulations of the early evolution of gas-free young clusters, similarly find that clusters in the process of merging appear more elongated. We thus argue that the size and shape of observed embedded clusters can inform our understanding of their recent history but not of their future evolution.

4.3 Implications for larger-scale simulations

Simulations of YMC formation with hydrodynamics and stellar feedback require very high gas masses for the initial GMC (three orders of magnitude above what we consider here, around 10⁷ M_⊙) and thus often model subgrid clusters with sink particles that can grow in mass by merging with other sinks and accreting gas (e.g. Howard et al. 2016, 2018). Karam & Sills (2022) have highlighted some of the limitations of this model, by showing that collisions between clusters do not always result in a single, merged cluster and that even when they do, the bound mass of the resulting cluster is less than the sum of the bound masses of the progenitors. They also find that cluster radii grow following a merger. We reinforce here those conclusions, and further note that groups of recently formed stars identified as cluster members—that would form within a subgrid cluster sink—can escape a cluster and can even be identified as a new cluster later in the simulation if they escape together. In particular, clusters can lose up to |$\sim 50~{{\ \rm per\ cent}}$| of their stellar mass if they split, and up to |$\sim 30~{{\ \rm per\ cent}}$| without splitting. A significant fraction of the stellar mass formed or accreted by a cluster can be lost on pre-supernova time-scales, which is not accounted for in cluster sink models.

Our embedded clusters tend to build up their mass mostly by forming new stars within the cluster, although they can accrete up to |$\sim 50~{{\ \rm per\ cent}}$| of their mass in already formed stars. Both processes contribute to the clusters’ growth in mass on time-scales much shorter that the clusters’ relaxation times. The dynamical evolution of the clusters remains driven by gravitational processes on the scale of the full simulation domain, such as the collapse of the gas, rather than by internal processes. This is a plausible cause of the diversity of cluster histories within the same simulation, as each individual cluster forms in a different local environment. Howard et al. (2016) found a similar spread for their cluster sinks: for their clusters in the 10²–10³ M_⊙ mass range, similar to our simulated embedded clusters, they find that between 0  per cent and ∼ 60  per cent of the clusters’ stellar mass is accreted.

Approximating embedded clusters as relaxed, spherical collections of gas and stars does not give an accurate representation of the clusters’ dynamical state. Furthermore, using spheres as a proxy for the shape of embedded clusters—or as a tool to measure cluster size, e.g. from Lagragian radii—may not be appropriate, as our clusters generally have ellipticities around 0.5, which indicates a factor of 2 difference between the major and minor axes of the stellar distribution.

4.4 Directions for future work

The simulation time for which we can evolve our models is currently limited by the high computational cost associated with following the dynamics of a large number of close binaries concurrently with radiative transfer and hydrodynamics. Including a self-consistent treatment of binary dynamics in simulations of embedded clusters as they reach gas expulsion is however essential to advancing our understanding of how star clusters form in galaxies. Although our results here indicate that gas dynamics dominate in the deeply embedded phase of cluster formation, the effects of binaries on the dynamics of clusters during gas expulsion remain unknown, and binaries are known to have an important impact on the evolution of gas-free clusters (e.g. Heggie 1975; Hills 1975; Torniamenti et al. 2021). Pursuing similar simulations with a large number of close binaries over time-scales sufficient to reach gas expulsion is therefore our next goal. This will require the use of a different N-body and few-body solver, to replace Ph4 and Multiples.

Directions for future work also include improvements to the treatment of gas and the stellar feedback in our simulations. Magnetic fields are not used in the current work due to their high computational cost. Future work will include comparisons of simulations with and without magnetic fields, as they are known to participate in the regulation of star formation (Price & Bate 2008). We also note that the amount of mass injected by wind feedback in our simulations is an upper limit, since the Vink et al. (2000) prescription for mass loss rates is likely too high by a factor of ∼3 (Smith 2014) and our winds are mass-loaded to avoid extremely short time-steps (Wall et al. 2020). The shock fronts in compact colliding wind binaries are not resolved due to our gas spatial resolution, such that we underestimate the heating from the winds. Any modulation of the feedback coming from interacting binaries is neglected. Including feedback from binaries is non-trivial, but is something we hope to address in future work. Our simulations also currently do not include protostellar jets and outflows. We expect the caveats outlined above to affect the spatial distribution of the feedback in our simulations, but not to significantly under- or overestimate the overall feedback budget.

5 SUMMARY

We have conducted a suite of hydrodynamics simulations of star cluster formation with a state-of-the-art treatment of stellar dynamics down to the scale of individual binaries, as well as active star formation via sink particles, and stellar feedback. We have explored a range of realistic initial virial parameters α = 0.8, 2.0, 4.0, at a fixed initial cloud mass of 10⁴ M_⊙, five different models for the formation (or not) of primordial binaries, and seven different random seeds for stochastic star formation. Most of our simulations have progressed to 2.0 Myr after the onset of star formation, which is the same as the time-scales we considered in Cournoyer-Cloutier et al. (2021). This allows us to investigate the relative impacts of the cloud-scale gas environment and internal two- or few-body dynamics while gas dynamics are still dominated by the gravitational collapse of the gas and not yet by the effects of stellar feedback. We have used a combination of tools to identify and characterize clusters, and arranged our analysis around three main axes: the properties of the full simulation domain (Section 3.1), the properties of the identified clusters as a population (Section 3.2), and the time evolution of individual clusters (Section 3.3). We have verified that the SFE of our simulation domains as well as the sizes, densities, and ellipticities of our embedded clusters are consistent with observations.

We explored the relative impact of the cloud’s initial virial parameter α and stellar dynamics (using the presence of primordial binaries as a proxy) on cluster structure and evolution. We have found the following:

• The choice of initial virial parameter α has the largest systematic effect on the global properties of the simulation domain, such as the SFR and SFE.

• The presence of primordial binaries or individual massive stars causes scatter in the SFR and SFE, but no systematic effect. The scatter is smaller than the systematic effects caused by changes in α. Stochastic effects from individual stars are important due to the low cluster masses (≲1000 M_⊙) considered in our simulations.

Our simulated embedded clusters are not relaxed, as their mass changes due to accretion or star formation on time-scales significantly shorter than their relaxation times. We thus find that their dynamical evolution is driven by the local gravitational potential (from the gas and stars) rather than by two- or few-body encounters (and therefore the presence of binaries). We have also tracked how cluster structure evolves during the earliest stages of cluster formation. We find considerable variation in cluster histories; examples are shown in Fig. 6. We summarize our results on cluster evolution as follows:

• Cluster mass generally grows through star formation rather than accretion, although some individual clusters acquire up to half of their final mass by accretion.

• The mass of individual clusters generally grows exponentially, although this growth is not monotonic. Clusters can lose up to half of their mass while they assemble.

• The size, density, and ellipticity of clusters do not follow any particular trend as the cluster acquires more mass. Changes in size, density, and ellipticity can take place over time-scales as short as 0.01 Myr.

• Recent accretion coincides with simultaneous sharp increases in characteristic radius and ellipticity. We propose that observed embedded clusters with high ellipticities are in the process of accreting stars.

The earliest stages of star cluster formation, when stars are still embedded in their natal gas and stars are still actively forming, are driven by a variety of competing physical processes; the structure of embedded star clusters changes quickly. We caution observers that the state in which an embedded cluster is observed is instantaneous. Over the time-scales considered in this work, cluster dynamical evolution is driven by the overall gravitational potential of the star-forming region, as individual clusters acquire new stars on time-scales much shorter than their relaxation times.

Acknowledgement

We warmly thank Marta Reina-Campos and James Wadsley for useful discussion. CCC is supported by a Canada Graduate Scholarship—Doctoral from the Natural Sciences and Engineering Research Council of Canada (NSERC). CCC also acknowledges funding from a Queen Elizabeth II Graduate Scholarship in Science and Technology (QEII-GSST) for the 2021–2022 academic year. AS and WEH are supported by NSERC. SA is grateful for funding from NSF grant AST-2009679. BP is supported through a fellowship from the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). AT was partly supported by NASA FINESST 80NSSC21K1383. AT was also supported through a NASA Cooperative Agreement awarded to the New York Space Grant Consortium. AT, MMML, and SLWM were partly supported by NSF grant AST18-15461. SLWM was also supported by NSF grant AST18-14772. MW acknowledges funding from NOVA under project number 10.2.5.12.

The simulations were conducted on the Digital Research Alliance of Canada supercomputer Graham under the resource allocation RRG #4398: The Formation of Star Clusters in a Galactic Context, and on the supercomputer Cartesius under the Dutch National Supercomputing Center SURF grant 15520. In addition to the software cited in-text, we have made use of the matplotlib (Hunter 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), and yt (Turk et al. 2011) Python packages for plotting and analysis.

DATA AVAILABILITY

The data underlying this paper will be shared on reasonable request to the corresponding author.

References

APPENDIX A: BINARY PRESCRIPTIONS

• Field distribution.    This is our fiducial distribution, based on statistics for all companions to main-sequence stars in the Galactic field. It is presented in detail in Cournoyer-Cloutier et al. (2021) and is based on observations compiled by Moe & Di Stefano (2017) and Winters et al. (2019).

• 10 per cent random pairing.    This prescription is based on that used in Sills & Bailyn (1999); similar prescriptions continue to be used in current state-of-the-art N-body or Monte Carlo simulations of massive star clusters (see e.g. Kamlah et al. 2022; Wang, Tanikawa & Fujii 2022). It imposes a mass-independent binary fraction of 10 per cent, with a period drawn from a flat distribution in log P (between 0.5 and 7.5, in days), and an eccentricity drawn from a thermal distribution. This model tends to underproduce binaries compared to the Galactic field, but nonetheless contains low-mass binaries that do not form naturally in models without primordial binaries.

• 100 per cent random pairing.    This prescription is the same as the one described above, with a binary fraction of 100 per cent at all masses.

• Field distribution for M < 0.6M_⊙ and no close massive binaries.    This prescription is also based on the algorithm presented in Cournoyer-Cloutier et al. (2021), but shifts all the periods to higher values |$\tilde{P}$| for stars with masses above 0.6 M_⊙ following |$\tilde{P} = 10^P$|⁠, where P is the period drawn from the algorithm. The specific choice of period shift is motivated by a typo we found in the binary generation algorithm we used in Cournoyer-Cloutier et al. (2021), which caused us to draw log P instead of P from our observation-based distribution. This typo does not affect the conclusions of the previous paper, as those were drawn from comparisons to the distribution of formed binaries (and not from comparisons to observations).

APPENDIX B: ELLIPSOIDS FROM INERTIA TENSORS

We present in Fig. B1 an example of 3D ellipsoidal surfaces enclosing 50 per cent and 90 per cent of the cluster mass, compared to the 50 per cent and 90 per cent Lagragian radii for the same stellar distribution. We use the reduced inertia tensor (Thob et al. 2019),

where the individual elements I_ij are calculated from

and

where |$\vec{x}_{a}$| is the vector distance from star a to the cluster’s centre of mass. The reduced inertia tensor minimizes the impact of stars in the outskirts of the cluster on the calculated shape. We obtain the principal axes a, b, and c from the eigenvalues λ_i of the reduced inertia tensor, such that

Solving this system of equations, we recover initial guesses for the principal axes:

The initial guesses from equations (B4) and (B5) are then rescaled iteratively to enclose 50  per cent or 90  per cent of the stellar mass. We adopt the idea of iterative fitting from Thob et al. (2019), and adapt the 2D code from Hill et al. (2021) to handle 3D distributions. The main steps of the fitting algorithm are as follows:

• Identify the centre of mass of the cluster;

• Take the stars enclosed within a given Lagrangian radius and get the shape for this distribution from the reduced inertia tensor in equations (B1) and (B2);

• Increase or decrease the size of the ellipsoid until the required (50  per cent or 90  per cent) fraction of the mass is enclosed;

• Recalculate the shape from the inertia tensor associated with the stars now enclosed in the ellipsoid;

• Repeat the steps above until the change in shape is less than a given tolerance between two iterations.

The default tolerance is 1  per cent but we raise it to 2  per cent for systems with fewer than 500 (but at least 200) stars and to 5  per cent for systems with fewer than 200 stars. We do not fit an ellipsoid when the most massive star in a cluster accounts for more than 50  per cent of its mass. We encounter this situation for one cluster in a few consecutive snapshots in S-R2.