Grouped star formation: converting sink particles to stars in hydrodynamical simulations

ABSTRACT

1 INTRODUCTION

Sink particles, first introduced by Bate, Bonnell & Price (1995) and subsequently developed further by Hubber, Walch & Whitworth (2013) and Bleuler & Teyssier (2014), are used to replace high density regions in simulations, as the individual time-steps of the high density gas particles become too small to follow. Sink particles are as such commonly used to represent the sites of star formation in simulations from individual stars or dense cores (e.g. Bate, Bonnell & Bromm 2003; Price & Bate 2008; Bate 2009, 2012; Lomax, Whitworth & Hubber 2015; Boss 2019), small groups of stars (e.g. Girichidis et al. 2011; Federrath & Klessen 2012; Balfour et al. 2015; Bertelli Motta et al. 2016; Ali & Harries 2019; He, Ricotti & Geen 2019; Ntormousi & Hennebelle 2019; Dobbs, Liow & Rieder 2020; Liow & Dobbs 2020; Dobbs & Wurster 2021), up to whole clusters (e.g. Vázquez-Semadeni et al. 2007; Renaud, Bournaud & Duc 2015; Howard, Pudritz & Harris 2018; Bending, Dobbs & Bate 2020; Ali 2021). By replacing high density regions with sink particles, this method ensures a cheaper and more efficient computation. Sink particles are usually decoupled from the physics of hydrodynamics and only interact with each other and gas particles via gravity. They can be used to incorporate physics such as following accretion on to forming protostars (e.g. Bate et al. 1995; Bonnell et al. 1997; Federrath et al. 2010; Dale et al. 2011; Padoan & Nordlund 2011; Stamatellos, Whitworth & Hubber 2012; Latif & Volonteri 2015; Jones & Bate 2018; Kuznetsova, Hartmann & Heitsch 2020), and disc formation (e.g. Stacy, Greif & Bromm 2010; Greif et al. 2012; Hennebelle et al. 2020).

In some smaller length scale and high mass resolution (i.e. low mass per gas particle) simulations (e.g. Bate 2009, 2012; Lomax et al. 2015), sink particles are formed when fragmentation is resolved at the opacity limit (Masunaga & Inutsuka 2000), so these sink particles are well resolved as individual stars. None the less, most hydrodynamical simulations do not have such high mass resolution, so sink particles are often introduced when the Jeans mass criterion is still satisfied (Bate & Burkert 1997), which is usually at much lower densities than the opacity limit for fragmentation. This means that even though each sink is a self-collapsing star-forming region, they represent small groups of stars and not individual protostars. Thus, cluster evolution cannot be studied as well compared to using simulations with resolved star particles (e.g. Bate 2012; Fujii & Portegies Zwart 2015). This scenario also assumes that any stars that are formed are always bound within their parent sinks, which is commonly assumed in many stellar feedback simulations (e.g. Dale et al. 2014; Geen et al. 2016, 2018; Bending et al. 2020). This situation would not capture, for example, the mass segregation of massive stars towards the centre of the larger clusters, or expulsion of stars due to close interactions, which may affect the evolution of clusters.

Instead of resolving sinks as stars, a workaround is to form star particles either using a local star formation efficiency function (Fujii & Portegies Zwart 2015) or a pre-determined initial mass function (IMF; Hu et al. 2017; Lahén et al. 2020; Ballone et al. 2021; Hirai, Fujii & Saitoh 2021; Hislop et al. 2021; Smith 2021). This allows the introduction of star particles in lower resolution simulations. None the less, mass is unable to be conserved locally in Fujii & Portegies Zwart (2015) and Smith (2021), so this can affect the stellar distribution and thus the modelling of stellar feedback, if included, would be inaccurate (Dinnbier & Walch 2020). Even though the method by Ballone et al. (2021) is excellent in conserving global and local stellar properties, it is not developed to form stars dynamically in hydrodynamical simulations and therefore unable to model the dense gas distribution. The star formation method by Hirai et al. (2021) that uses density-independent formulation of smoothed particle hydrodynamics (SPH) for gas hydrodynamics (Saitoh & Makino 2013) can be applied to simulations of various length scales, however it can create gas with unequal masses, which is by design not permitted in some codes that implement standard SPH (e.g. Price & Monaghan 2007; Price et al. 2018). Another way to introduce stars in the simulation is via ‘sink-star’ hybrid particles as demonstrated in Grudić et al. (2021), whereby a protostar is embedded in each sink and only interacts with the simulation through accretion from its sink reservoir and stellar feedback, but this method does not resolve stellar dynamics within the sinks.

A novel way to create stars from sinks in low mass resolution simulations was introduced by Wall et al. (2019) and is hereafter referred to as single star formation. Sink particles are formed using the prescription by Federrath et al. (2010), then each sink particle forms stars that are sampled from the user-input IMF using Poisson sampling (Sormani et al. 2017). The stellar population from each sink is self-consistent, i.e. each stellar population follows the IMF. This method of star formation ensures that stellar dynamical interaction happens not only between stars formed within the same sink but also between stars formed in different sinks. Moreover, local mass conservation is obeyed and the star formation prescription causes minimum disruption to the gas evolution code. None the less, to sample a complete IMF, each sink is expected to have ≳10² M_⊙ in mass (Wall et al. 2019; Rieder et al. 2021; Smith 2021), which is usually assumed to be a ‘cluster sink’ and only achievable in larger cloud or (sub-)galactic simulations, or simulations with very low mass resolution. As we investigate in this paper, applying the single star formation method on parsec-scale or smaller simulations, in which sinks are usually assumed to be small groups of stars, results in oversampling of low mass stars and undersampling of high mass stars. These simulations usually have a mass per gas particle of about 10⁻¹−10⁻⁴ M_⊙. Consequently, no sink is massive enough to sample a complete IMF, while the sinks are not small enough to be resolved as individual star either.

In this paper, we extend the method of converting sinks to stars as presented in Wall et al. (2019) and Rieder et al. (2021) to apply to higher mass resolution simulations if sinks cannot be resolved as individual stars, but are individually not large enough to sample a complete IMF. We perform grouped star formation, i.e. we group the sink particles and use the group mass to sample the IMF. We first group the sink particles according to certain length, speed, and time-scales so that each group is approximately a consistent star-forming region. After sampling the IMF using the group mass, stars are placed within the sinks according to the sink masses to optimize the local conservation of mass. This method allows lower mass sinks to be grouped so that the group mass is large enough to sample the IMF robustly. This means that we can afford the formation of smaller size sinks so that the stellar distribution can mimic the dense gas distribution more accurately. The grouped star formation method also allows us to test the limit of the single star formation method.

This paper is sectioned as follows: in Section 2, we lay out the prescription for grouped star formation, and the simulation setup to study this method in different length scales. The results are shown in Section 3. We explore the effect of changing the random seed for star formation and varying the upper star mass limit of the IMF in Section 4. Lastly, we summarize the paper in Section 5.

2 NUMERICAL METHODS

2.1 Methods for group assignment and star formation

In general, the grouped star-forming prescription is divided into two steps: the grouping of sink particles and the conversion of sink particle mass into star particles. The group assignment occurs when new sink particles are formed, while the conversion of sink mass to stars happens whenever there are sink particles in the simulation. The newly formed sinks are first arranged according to their masses in descending order. When the first sink particle is formed in the simulation, it forms a group by itself, since no group existed prior to that. After the first sink, for the a-th newly formed sink particle, we loop through all existing groups and add the new sink to a pre-existing sink group g, if:

• sink particle a is within a distance d_g away from the group g’s centre-of-mass (COM), i.e. |$|\boldsymbol {r}_a - \boldsymbol {r}_{\text {COM},g}| \le d_g$|⁠, where |$\boldsymbol {r}$| is the position vector,

• sink particle a is within a speed v_g of the group g’s centre-of-mass velocity, i.e. |$|\boldsymbol {v}_a - \boldsymbol {v}_{\text {COM},g}| \le v_g$|⁠, where |$\boldsymbol {v}$| is the velocity vector,

• the creation time of sink particle a is at most τ_g greater than the oldest member in group g, and

• the sink particle a is the most bound to group g compared to the other groups.

For ease of reference, henceforth the first three conditions are called the distance, speed, and age criteria, respectively. These criteria are set so that each group is approximated as a self-consistent star-forming region, such that |$d_{\text{g}}$|⁠, |$v_{\text{g}}$|⁠, and |$\tau _{\text{g}}$| are the length, speed, and time-scales of a typical star-forming region, respectively. The speed criterion ensures that any sinks that are likely moving away from the group are not included within the group, while the age criterion makes sure that sinks formed at later times are not grouped together with sinks formed much earlier, even though they could coincidentally be located near each other. The fourth condition is not essential but necessary in case a sink particle satisfies the grouping criteria of more than one group. If the a-th sink particle fails to satisfy all of the above conditions, then it will create a new group by itself, and the process continues for the other sink particles. Each sink particle is assigned to only one group and will be in the same group in subsequent time-steps to keep track of the stellar population in each group.

For each group, the sink particle masses are reduced according to the total star mass located within them. However, although unlikely, the total stellar mass can be greater than the mass of the sink the stars are located in, due to the probabilistic nature of allocating stars. To overcome the issue of having negative sink mass, the sink mass is set to |$m_{\text{thres}}$| (here, |$m_{\text{thres}}$| is chosen as the mass of a gas particle) when the difference in total star mass and their parent sink is less than |$m_{\text{thres}}$|⁠. This mass difference is accumulated for all sinks in the group and is then reduced from all sinks in the group in proportion to their masses. The mass conservation of the group is guaranteed using this scheme. Lastly, the sink particles shrink in size to keep their creation density constant. In subsequent time-steps, when a new sink joins an existing group, its mass is added to the leftover mass of the sink group. The new total group mass is then used to sample the IMF and make new stars, and finally, the masses of the sinks in the group are reduced as described.

2.2 Simulation setup

We use EKSTER (Rieder & Liow 2021; Rieder et al. 2021) to perform our simulations. EKSTER is a multiphysical code that combines gas hydrodynamics, gravitational dynamics, and stellar evolution via the AMUSE interface (Portegies Zwart et al. 2018). The gravitational dynamics between gas and non-gas (sinks and stars) are coupled using BRIDGE (Fujii et al. 2007). For gas hydrodynamics, we use PHANTOM (Price et al. 2018), a smoothed particle hydrodynamical (SPH) code for astrophysics. Artificial viscosity is included to model shocks (Monaghan 1997). The standard values of artificial viscosity parameters |$\alpha_{\text{min}}^{\text {AV}} = 0.1$|⁠, |$\alpha_{\text{max}}^{\text {AV}} = 1$|⁠, and |$\beta ^{\text {AV}} = 2$| are used (Morris & Monaghan 1997) except for Section 2.3.2 where higher values are used. PETAR (Wang et al. 2020) is used for the fast calculation of gravitational dynamics between and among sink particles and stars. Lastly, we use SEBA (Portegies Zwart & Verbunt 1996), a parametric code to calculate stellar evolution. Sink particles are introduced to replace high density gas particles, following the prescription used in Bate et al. (1995) and Price et al. (2018), when the gas density exceeds the sink creation density |$\rho _{\text{sink}}$|⁠. Depending on the system studied, we adjust |$\rho _{\text{sink}}$| such that the minimum Jeans mass criterion is satisfied (Bate & Burkert 1997). Sinks are allowed to accrete mass.

2.3 Initial conditions

We test the grouped star formation prescription on simulations of different length scales and mass resolutions as shown in Fig. 1. The aim is to demonstrate that grouped star formation method is needed at the mass resolution regime in which each sink is usually considered as a small group of stars. In each comparison described below, we vary the grouping parameters as described in Section 2.1. The ‘no grouping’ case, i.e. |$d_{\text{g}} = v_{\text{g}} = \tau _{\text{g}} = 0$|⁠, is equivalent to the single star formation prescription used in Wall et al. (2019) and Rieder et al. (2021). On the other hand, the ‘all grouping’ case simply groups all sinks as one, and we achieve this by setting |$d_{\text{g}}$|⁠, |$v_{\text{g}}$|⁠, and |$\tau _{\text{g}}$| unphysically large, usually |$d_{\text{g}} = 1000$| pc, |$v_{\text{g}} = 1000$| km s⁻¹, and |$\tau _{\text{g}} = t_{\text{ff}}$| if the period of interest is within a free-fall time, or |$\tau _{\text{g}} = 1000 \ t_{\text{ff}}$| otherwise. The main other settings that we explore are the ‘standard’ case, i.e. |$d_{\text{g}} = 1$| pc, |$v_{\text{g}}$| = 1 km s⁻¹, and |$\tau _{\text{g}}=t_{\text{ff}}$|⁠, where |$t_{\text{ff}}$| is the free-fall time-scale of the system, and the ‘turbulent’ case, i.e. the same as the‘standard’ case except with |$v_{\text{g}}$| equals to the initial turbulent speed (coincidentally, it is approximately 3 km s⁻¹ in most of our different comparisons). The standard case is chosen as such because they are convenient and appropriate choices for the length and speed scales of a typical star-forming region. Besides, the mass estimate for the group, in this case, |${\sim}d_{\text{g}} v_{\text{g}}^2/G \approx 230$| M_⊙ turns out to be the same order of magnitude of the mass expected to sample a complete IMF (Wall et al. 2019; Smith 2021). We also test other grouping parameters specific to individual comparisons and they are explained in detail in the subsequent subsections. We list the different models we use and the different grouping parameters for each in Table 1.

Figure 1.

The relationship of length scale and mass resolution is shown for the original simulations and models compared in this paper (see the text for details of different simulations presented). The plot is roughly divided into three regions by the mass thresholds 2 × 10⁻⁵ M_⊙and 1 M_⊙ which show the nature of sink particle (bold text) in the simulations. Resolution is traditionally defined as the inverse of mass per gas particle.

Open in new tabDownload slide

2.3.1 Sub-parsec-scale isolated cluster simulation

We first test our sink method on a simulation of cluster formation in a 500 M_⊙ isolated cloud. We use the simulation of Bate (2012) who perform a 35 million particle (mass per gas particle |$m_{\text{SPH}}\approx 1.4 \times 10^{-5}$| M_⊙) simulation which is able to resolve star formation down to a mass of ∼10⁻² M_⊙. We choose this simulation, which includes radiative transfer to model stellar heating (Whitehouse, Bate & Monaghan 2005), rather than Bate (2009), even though we do not include this additional physics. This is partly because the data were readily available, and also because with the radiative feedback, the IMF produced is much closer to a Kroupa IMF which matches our choice of IMF. We describe the initial conditions below. The aim of our simulations is to reproduce the properties of the stars in the original simulation of Bate (2012), but using a gas resolution much lower than the original simulation. We then compare the results from Bate (2012) with our lower resolution simulations performed with EKSTER. Our initial conditions are the same except for the number of particles. The cloud has a mass of 500 M_⊙ and a radius of 0.404 pc initially, which gives an initial cloud density of 1.2 × 10⁻¹⁹ g cm⁻³ and a spherical free-fall time of 0.19 Myr. We apply exactly the same velocity field as Bate (2012), which has a velocity dispersion of ≈3 km s⁻¹.

Bate (2012) started with an initial temperature of 10.3 K and used SPH with radiative transfer in the flux-limited diffusion approximation (Whitehouse et al. 2005) in their simulation. As their sink particles are formed during the second collapse of protostellar cores (Masunaga & Inutsuka 2000), they are well resolved as stars and brown dwarfs. On the other hand, we simply adopt the isothermal equation of state at a temperature of 10 K as our chosen |$\rho _{\text{sink}} = 10^{-16}$| g cm⁻³ is within the first phase of protostellar collapse, which can be modelled as isothermal collapse (Masunaga, Miyama & Inutsuka 1998; Masunaga & Inutsuka 2000). Our resolution is |$N_{\text{SPH}} = 5 \times 10^6$| (⁠|$m_{\text{SPH}}=10^{-4}$| M_⊙) as opposed to |$N_{\text{SPH}} = 3.5 \times 10^7$| (⁠|$m_{\text{SPH}} \approx 1.4 \times 10^{-5}$| M_⊙) used in Bate (2012), however our resolution and |$\rho _{\text{sink}}$| allow us to resolve the Jeans mass criterion (Bate & Burkert 1997). Unless otherwise stated, our gravitational softening length ϵ in this comparison is set at 0.001 pc ≈200 AU. In Bate (2012), gravitational softening between sink particles is turned off. The initial accretion radius of sink particles |$r_{\text{acc}} = 0.001$| pc, the same order of magnitude of the Jeans length in this system.

Based on Bate (2012), we set the grouping parameters described in 2.1 as follows. For the ‘standard’ case, |$d_{\text{g}} = 1$| pc and |$v_{\text{g}} = 1$| km s⁻¹ are chosen as described. For the ‘turbulent’ case, the grouping parameters are the same as the ‘standard’ case except for |$v_{\text{g}} = 3$| km s⁻¹, which is approximately the turbulent velocity dispersion of this gas cloud. The single star formation method is exactly the same as setting both |$d_{\text{g}}$| and |$v_{\text{g}}$| to zero, i.e. the ‘no grouping’ case, as these settings prevent the formation of any groups at all. We also consider the opposite extreme, i.e. the ‘all grouping’ case by setting both |$d_{\text{g}}$| and |$v_{\text{g}}$| to unphysically large values of 1000 pc and 1000 km s⁻¹ respectively to group all the sinks into one single group.

We include the age criterion |$\tau _{\text{g}} = t_{\text{ff}}$| but in this comparison, this check is irrelevant as we only compare up to the free-fall time, so the sink particles and stars in all our simulations are definitely formed within that time-scale. Lastly, we also perform additional simulations with lower resolution and with no gravitational softening separately. The ‘low resolution’ model has |$N_{\text{SPH}}=10^6$| (⁠|$m_{\text{SPH}} = 5 \times 10^{-4}$| M_⊙). This model is used to check whether our star formation method works in a lower resolution of the same length scale. For the ‘no softening’ model, the gravitational softening length ϵ is turned off but the gas particle softening length remains at 0.001 pc, the same as in other simulations in this comparison. This setting of ϵ < gas smoothing length is cautioned by Bate & Burkert (1997) as fragmentation at scales less than the gas smoothing length can be artificially induced. However, we run this model to investigate the dynamics of stars formed from our star formation prescription. The details of the models for this comparison are listed in the upper section of Table 1.

For star formation, we sample stellar masses from the Kroupa IMF (Kroupa 2001). In Section 3.1, we take a mass range from 0.01 M_⊙ to 100 M_⊙ which extends down to the minimum mass of stars resolved in Bate (2012). However since 100 M_⊙ is large compared to the typical group size in the models in this comparison, in Section 4.2 we also test taking a mass range from 0.01 M_⊙ to the individual mass of the group.

2.3.2 Parsec-scale cloud–cloud collision simulation

A cloud–cloud collision simulation from Liow & Dobbs (2020) is chosen to test our star formation method on larger scales with higher sink masses. In this case, individual stars are not resolved, rather each sink is simply a small group of stars. Another difference compared to Bate (2012) is that in this larger scale simulation, sinks form at much wider distances from each other. The initial conditions are described in Liow & Dobbs (2020) (known as the model with low collision speed, standard density, and low turbulence) but we also describe them again below.

Two ellipsoidal clouds of mass 5 × 10⁴ M_⊙ and minor radii of 7 pc collide along the major radius of 16 pc. This gives an initial cloud density of 1.03 × 10⁻²¹ g cm⁻³ and a spherical free-fall time of 2.07 Myr. In the original simulation, the resolution is 5 × 10⁶ (⁠|$m_{\text{SPH}} = 0.02$| M_⊙), however the resolution we use is 10⁶ SPH particles (⁠|$m_{\text{SPH}} = 0.1$| M_⊙). The gas clouds are initially about 0.3 pc apart, and are subjected to two separate turbulent fields. In this setup, the turbulence has a velocity dispersion of 2.5 km s⁻¹, and the clouds are given a relative velocity of 10 km s⁻¹. We choose to compare our models at |$t_{10{{\ \rm per\ cent}}}$|⁠, the time when 10 per cent of the gas mass is converted to sink particles, as the sink particle distribution shows multiple star-forming regions, notably the central cluster at the collision site, and the filamentary structures that are in the parts of the clouds yet to collide. The ongoing collision also gives multiple velocity components to be considered. The external forces exerted by the colliding clouds cancel out at the collision site, giving a zero net momentum for the star-forming regions at the collision site, whilst sinks away from the area where the clouds are colliding are still moving at about 10 km s⁻¹ relative to each other.

In the original model by Liow & Dobbs (2020), sinks are formed at |$\rho _{\text{sink}} = 10^{-18}$| g cm⁻³ which was typically too low for them to be considered individual stars but enough to satisfy the Jeans mass criterion (Bate & Burkert 1997). Here, we use the same |$\rho _{\text{sink}}$| as our sink creation density. In this comparison, we also set the artificial viscosity parameter |$\beta ^{\text {AV}} = 4$| as for the original simulation to minimize the effect of particle penetration in high speed shocks (Price & Federrath 2010). Similar to the original model, the gravitational softening length and the initial accretion radius |$r_{\text{acc}}$| are set to 0.01 pc. We consider the ‘no grouping’, ‘standard’, ‘turbulent’ (with |$v_{\text{g}} = 3$| km s⁻¹, approximately the initial turbulent speed of the clouds), and ‘all grouping’ cases, similar to the comparison described in Section 2.3.1. The details of the models in this comparison are listed in the second part of Table 1.

In this comparison, we choose a mass range of 0.5–100 M_⊙, as it roughly reflects the mass distribution of sinks from the original simulation at |$t_{10{{\ \rm per\ cent}}}$|⁠, which is also our chosen time-scale for comparison. Unfortunately, unlike the models described in Section 2.3.1, we cannot resolve the full IMF with larger scale simulations, so we simply compare our models with lower resolution to the original simulation. We could, however, potentially choose a lower stellar mass limit for our model (see Buckner et al., in preparation, which has a lower stellar mass limit of 0.01 M_⊙).

2.3.3 Parsec-scale isolated cluster simulation

We adopt the grouped star formation method to a parsec-scale, gravitational collapse of an isolated cloud simulation originally performed by Jaffa et al. (in preparation). We do not rerun the simulation, rather here we post-process the results. The mass and radius of the cloud are 10⁴ M_⊙ and 10 pc, respectively, which gives an initial cloud density of 1.61 × 10⁻²² g cm⁻³ and a free-fall time of 5.24 Myr. The resolution used in the simulation is 10⁶ SPH particles, which corresponds to |$m_{\text{SPH}} = 0.01$| M_⊙.

The cloud has an initially uniform density and a turbulent velocity field, with a natural mix of compressive and solenoidal modes and a power spectrum of P(k) ∝ k⁻⁴. The turbulence is allowed to decay as the cloud collapses to form filamentary structures, which merge to feed a central cluster with some smaller structures around it. Sink particles are inserted when gas reaches a density of 10⁻¹⁸ g cm⁻³, with some further checks such as proximity to other sinks (see Hubber, Rosotti & Booth 2018). The gravitational softening applied to interactions between two sink particles is the mean of their sink radii, which are set to their SPH smoothing length, h ∼ 10⁻² pc. Sinks continue to accrete but are not allowed to merge, and by 10 Myr more than 60 per cent of the gas has been turned into sink particles. The simulation is terminated at 20 Myr, but we note that some important physical processes such as feedback are not modelled in this simulation so the later evolution of the cloud and cluster may be affected by this. The full details of the simulation can be found in Jaffa et al. (in preparation).

We use this simulation partly to test the grouping age criteria. In the simulations described so far in Sections 2.3.1 and 2.3.2, the age criteria is not relevant because the stars form in a short time-scale, within a free-fall time. This simulation however is run for 20 Myr |${\approx}3.5 t_{\text{ff}}$|⁠, and so provides a good test of the age criteria. Although the resolution is not that high, the 20 Myr time-scale means the calculation still required considerable computational resources to run, and so we do not repeat the calculation with different grouping parameters, rather we apply these post process. Therefore, instead of forming stars dynamically, we apply a modified grouped star formation prescription on the sinks from the original simulation at specific snapshots to form stars, keeping track of the group indices assigned to the groups at previous snapshots. Unfortunately, this means that we cannot analyse the velocity dispersion or the stellar distribution, since we do not follow the dynamics, but we can still consider the mass functions and properties of the groups. Similar to our previous comparisons, we choose models with different grouping parameters to represent ‘no grouping’, ‘standard’, ‘turbulent’ (⁠|$v_{\text{g}} = 3$| km s⁻¹, approximately the turbulent velocity dispersion of the cloud), and ‘all grouping’ cases. Here, the ‘all grouping’ case sets all three grouping parameters (⁠|$d_{\text{g}}, v_{\text{g}}$|⁠, and |$\tau _{\text{g}}$|⁠) unphysically large. To investigate the age criteria more robustly, we include the ‘no age check’ model, whereby the grouping parameters are the same as the ‘standard’ model except that |$\tau _{\text{g}}$| is set unphysically high to prevent triggering this condition check. We also include the ‘free-fall’ model, whereby the grouping parameters are the same as the ‘all grouping’ model except that |$\tau _{\text{g}} = t_{\text{ff}} = 5.24$| Myr, since our comparison is up to 20 Myr (⁠|${\approx}3.5 t_{\text{ff}}$|⁠). In this comparison, we use the Kroupa IMF with mass range 0.5–100 M_⊙. The parameters of the models are listed in the third part of Table 1.

2.3.4 Kiloparsec-scale spiral arm simulation

For completeness, we take the stars from the spiral arm simulation by Rieder et al. (2021) and study the mass functions. This simulation models star cluster formation along a section of spiral arm. The simulation (labelled standard-arm-1 in Rieder et al. 2021) models a region of dimensions 600 pc, which includes a number of molecular clouds. The mass resolution is |$m_{\text{SPH}}=1$| M_⊙.

We do not rerun this simulation because in the original simulation, Rieder et al. (2021) uses the single star formation method to form stars. If the IMF of the stars in the original simulation is already close to a Kroupa IMF, then we can assume that the single star formation method is suitable for simulations of this length scale and mass resolution, or larger. In the original simulation, the Kroupa IMF is used as the input IMF, and the chosen mass range is 0.1 to 100 M_⊙.

3 RESULTS AND DISCUSSION

3.1 Sub-parsec-scale isolated cluster simulation

Fig. 2 shows the IMF of the models to compare with Bate (2012), i.e. model B12, at the free-fall time of 0.19 Myr. The IMF of B12 is similar to the Kroupa IMF. The ‘no grouping’ model (orange) that uses the single star formation to form stars is unable to sample a complete Kroupa IMF. The average mass reservoir to form stars, i.e the average sink mass, is 21.59 M_⊙/1126 = 0.02 M_⊙, meaning not many stars above this threshold are expected to form, causing an oversampling of low-mass stars and undersampling of high mass stars. Moreover, as shown in Table 1, the star mass fraction, defined as the fraction of star mass over the total of sink and star masses, is f = 0.05, so most of the mass in each sink is unable to be converted to stars. The value of f is low because for a population of stars to be introduced for each group (in this model, each group is a sink), the group mass has to exceed a mass threshold sampled from the Kroupa IMF. This is difficult for this model because the probability of forming a star less than 0.02 M_⊙ is only 7 per cent. We discuss setting the upper star mass limit to be equal to the group mass as a potential improvement to produce a higher f in Section 4.2.

Figure 2.

The IMF of our models are compared with that of Bate (2012), labelled as B12, at the free-fall time of 0.19 Myr. The markers are the maximum mass bin in the respective models. The errorbars are the standard deviations from Poisson sampling. The black dashed line is the analytical Kroupa IMF.

Open in new tabDownload slide

With grouped star formation of varying degrees, the maximum star mass bin (markers in Fig. 2) is greatly increased. The ‘standard’ model (green) with |$d_{\text{g}}=1$| pc and |$v_{\text{g}}=1$| km s⁻¹ is still unable to sample the full Kroupa IMF. We find that by setting |$v_{\text{g}}$| to 3 km s⁻¹, about the initial root-mean-square turbulent velocity of the cloud, while keeping the same |$d_{\text{g}}$|⁠, the ‘turbulent’ model (red) can sample the full Kroupa IMF like the ‘all grouping’ model (purple) that collects all sinks as one giant group, simply because setting the grouping parameters that match the length scale and turbulence of the initial cloud approximately covers the whole parameter space of this system. In all the models, by summing the mass of sinks and stars, we see that the total mass differs between models. However, the stellar mass in the ‘all grouping’ case most closely matches that of the B12 model, whilst the ‘no grouping’ model (single star formation) massively underproduces stellar mass as compared to the original B12 calculation.

We included two more models to test the effect of changing gravitational softening length and resolution on the IMF. The ‘no softening’ model (brown) has similar IMF as the ‘all grouping’ model simply because the total sink mass of the system is similar regardless of the stellar dynamics. Similarly, the ‘low resolution’ model (pink) is also able to sample a complete Kroupa IMF like the ‘all grouping’ model, showing the robustness of the grouped star formation at different resolutions.

Fig. 3 shows the stellar distributions of the ‘no grouping’, ‘standard’, ‘turbulent,’ and ‘all grouping’ models. The stars in the ‘no grouping’ model (top left-hand panel), formed from the single star formation method, are less visible due to their low mass. The same can be said for the ‘standard’ model (top right-hand panel), even though we can observe that the low mass stars are generally located at the region stars from B12 are expected. The ‘turbulent’ model (bottom left-hand panel) contains massive stars and reproduces the spatial stellar mass distribution reasonably well. Lastly, we note the ‘all grouping’ model (bottom right-hand panel) that groups all sinks together can form stars that are the most massive among all models, but the stellar distribution of this model deviates from that of B12 the most. On scales of around ∼10⁻² pc, mass is not conserved locally in this model. However, it is not necessarily important to model the exact spatial distribution of stars in a cluster in larger scale simulations where the minimum resolvable length scale may be >>10⁻² pc. The ‘no softening’ and ‘low resolution’ models are not shown in Fig. 3, however in terms of the degree of grouping, they are similar to the ‘all grouping’ model but with different gravitational softening or resolution, so they are similar to the ‘all grouping’ model in terms of their stellar distribution. In summary, the models that group the most sinks like the ‘turbulent’ or ‘all grouping’ models produce IMFs which match the IMF of B12. This is not surprising since the whole system is so small that it takes the whole region to sample the IMF. Thus, grouping most (which happens in the ‘turbulent’ model where grouping parameters approximately equal to the scales of the cloud) or all the sinks together (which happens in the ‘all grouping’ model) works best, and the single star formation method is unsuitable. Grouping most or all the sinks together means that the spatial distribution matches that of B12 less well, but this is less concerning in simulations where the gas resolution is low as compared to the size of the clusters.

Figure 3.

The stars from the ‘no grouping’, ‘standard’, ‘turbulent’, and ‘all grouping’ models (white) are compared against the sinks from Bate (2012) (blue) at the free-fall time of 0.19 Myr. The marker radius is proportional to the particle mass. The gas distribution from the original simulation is set as background for reference, as the gas distribution from the individual models are visually identical to the original one.

Open in new tabDownload slide

Fig. 4 shows the cumulative distribution velocity function of the stars for all our models as compared to that of B12 (blue). We included a Gaussian distribution (black dashed line) with arbitrary mean of 1.6 km s⁻¹ and standard deviation of 0.4 km s⁻¹ to show that all our star formation models have velocity dispersion that is approximately Gaussian (with varying means and variances), as expected given how velocity is assigned to the stars (Wall et al. 2019). The models with a greater degree of grouping have velocity dispersion closer to B12, but even with models that group all sinks as one like the ‘all grouping’ model (purple), they are still visually different from the velocity dispersion of B12. We use the Kolmogorov−Smirnov (KS) test with permutation (Præstgaard 1995) to test the hypothesis of equal distribution between the velocity dispersion of our models with that of B12. With a significance level of 0.01, the result shows that all models except the ‘no softening’ model (brown) have velocity dispersions that are different from that of B12. Statistically, this means that the stellar dynamics of the ‘no softening’ model and the original simulation of B12 could be similar, which would likely be caused by the similar stellar mass distribution and the fact that both have gravitational softening turned off. With gravitational softening, dynamical interactions between the stars are not captured properly, and the velocities of the stars will be lower as a result. We performed an additional no softening model but with grouping parameters the same as the ‘standard’ model (not shown here) to confirm that the mass distribution also plays a role in determining the velocity dispersion as expected.

Figure 4.

The cumulative velocity distribution of the models are shown are compared with that of Bate (2012) at 0.19 Myr. The black dashed line is the Gaussian cumulative distribution function with mean 1.6 km s⁻¹ and standard deviation of 0.4 km s⁻¹. It is not a fit and is included solely for comparison.

Open in new tabDownload slide

Fig. 5 shows the mass and age spread of the star groups (recall that sinks are first grouped together, then stars form in each group) for the ‘no grouping’, ‘standard’, ‘turbulent’, and ‘all grouping’ models at the free-fall time of 0.19 Myr. The models are arranged in increasing degree of grouping, which can be roughly deduced from the star mass fraction f shown in Table 1. In general, as the degree of grouping increases, we see an increase in the median mass of star groups. For the ‘no grouping’ model, a great proportion of mass is ‘trapped’ in the sinks and is unable to form stars, which causes the low star group masses. We do not see a clear trend in the group age spread in our models in comparison of this length and time-scales. There is no age spread for the groups in the ‘no grouping’ model as they are single-member groups. On the other hand, for the ‘all grouping’ model, the age spread of the star group simply shows the duration between the first star formation and the current free-fall time. We also analyse the virial radius and the velocity dispersion of the star groups, and they simply reflect the length scale (∼10⁻² pc as shown in Fig. 3) and speed scale (≈3 km s⁻¹, same as the turbulence of the cloud) of the star groups.

Figure 5.

The group mass and age spread of the star groups for each model at the free-fall time of 0.19 Myr are shown. For each box, the orange line is the median group mass, the length of the box is the interquartile range, and the whiskers span the whole range.

Open in new tabDownload slide

3.2 Parsec-scale cloud–cloud collision simulation

Fig. 6 shows the IMFs of all the models that compare with the chosen simulation from Liow & Dobbs (2020), labelled L20 (blue) at 1.75 Myr, the time when 10 per cent of the gas mass is converted to sinks and stars, the same time-scale shown in the figures in Liow & Dobbs (2020). Although the L20 simulation has a much lower resolution than B12 (see Section 3.1), the form of the sink mass function is not dissimilar to the Kroupa IMF at high masses; at low masses the mass function departs strongly from the Kroupa IMF. First, the ‘no grouping’ model (orange) that uses single star formation is not able to form any high mass stars >10 M_⊙, and deviates significantly from the Kroupa IMF at the high mass range. With grouped star formation, the resultant IMFs are closer to the Kroupa IMF, however the ‘standard’ model (green) with |$d_{\text{g}}=1$| pc and |$v_{\text{g}}=1$| km s⁻¹ is still not enough to sample a complete Kroupa IMF. Similar to the comparison with Bate (2012), both the ‘turbulent’ model (red) that uses the initial turbulent velocity as the grouping speed, and the ‘all grouping’ model (purple) can sample the full Kroupa IMF. However the ‘turbulent’ model does not need to group all the sinks together like the ‘all grouping’ model, which as we shall see is beneficial for larger scale simulations. The star mass fraction f = 0.901 of the ‘turbulent’ model is also high compared to other models.

Figure 6.

The IMF of our models are compared with that of Liow & Dobbs (2020), labelled as L20, at 1.75 Myr, the time when 10 per cent of the gas mass is converted to sink mass in the original simulation. The markers are the maximum mass bin in the respective models. The errorbars are the standard deviations from Poisson sampling. The black dashed line is the analytical Kroupa IMF.

Open in new tabDownload slide

Fig. 7 shows the stellar distribution of the ‘no grouping’, ‘turbulent’, and ‘all grouping’ models plotted over the sink distribution of L20. Since each sink in L20 is best approximated as a small group of stars, the size of the sinks are exaggerated as compared to the stars from our models. For the ‘no grouping’ model (top plot), the stars are located about the dense gas filaments, however their masses are small as noted in Fig. 6. As the degree of grouping in grouped star formation increases, more higher mass stars that are comparable to the original simulations are formed, again already inferred in Fig. 6. However, at the same time, we note that the stellar distribution reflects the true dense gas filament distribution less as grouping increases. This means that local mass conservation, i.e. within scales of order parsecs, becomes more violated as the degree of grouping increases. In the case of the ‘all grouping’ model (bottom plot) where all sinks are grouped as one, mass is unphysically transferred to the clustered area from less clustered regions. Similar unphysical transfer of mass to some extent can occur using the star formation schemes described in Fujii & Portegies Zwart (2015) and Smith (2021). This may be unwanted if the length resolution of a simulation is parsec or sub-parsec, whereas in the comparison with Bate (2012) in Section 3.1, the size scale of the whole system is still small compared to, for example, whole galaxy simulations.

Figure 7.

The stars from the ‘no grouping’, ‘turbulent’, and ‘all grouping’ models (white) are compared against the sinks from Liow & Dobbs (2020) (blue) at the 1.75 Myr. The marker radius is proportional to the particle mass. The gas distribution from the original simulation is set as background for reference, as the gas distribution from the individual models are visually identical to the original one. The orange dashed box shows the region of interest for the discussion of local mass conservation.

Open in new tabDownload slide

To quantify mass conservation, we take a box of length 4 pc centred at the origin, and calculate the percentage of total sink and star mass from each model (⁠|$M_{\text{sinks,model,box}} + M_{\text{stars,model,box}}$|⁠) on the total sink mass from L20 (⁠|$M_{\text{sinks,L20,box}}$|⁠). The location of the box contains about half of the most massive cluster obtained from the original simulation of L20 (see Liow & Dobbs 2020 for more information on the identification of this massive cluster, labelled ‘L4’), ideal to investigate local mass conservation. We find that for the ‘no grouping’ model, 97 per cent of the mass is conserved locally, however only about 65 per cent of the locally conserved mass are stars as this model does not form many massive stars, and a lot of mass is still ‘trapped’ in the sinks without any star formation. On the other hand, even though the ‘all grouping’ model is extremely efficient in converting sink mass to star mass (f = 0.999, so |$M_{\text{sinks}}$| and in extension |$M_{\text{sinks,model,box}}$| are practically negligible in this model) and producing massive stars, only about 40 per cent of mass is conserved locally. For the ‘turbulent’ model (middle plot), about 90 per cent of the mass is conserved locally, and about 90 per cent of the locally conserved mass is contributed by stars, suggesting that this model is the most successful in terms of forming a similar mass of stars to the original model, reproducing the Kroupa IMF, and obeying local mass conservation.

Fig. 8 shows the cumulative velocity distribution function for all models compared to that of L20. The velocity dispersion of L20 is fitted with a Gaussian distribution with mean 6.2 km s⁻¹ and standard deviation of 2.7 km s⁻¹ (black dashed line). As large sample sizes cause bias towards rejecting the null hypothesis for equal distributions, KS tests are not suitable for hypothesis testing here (Gómez-de Mariscal et al. 2021). None the less, all velocity distributions are visually Gaussian-like and similar to that of L20, regardless of the grouping parameters, suggesting that the velocity distribution of the larger scale models is relatively independent of the details of the star formation prescription.

Figure 8.

The cumulative velocity distribution of the models shown are compared with that of Liow & Dobbs (2020) at 1.75 Myr. The black dashed line is the Gaussian cumulative distribution function, fitted on the velocity distribution of the original simulation, with mean 6.2 km s⁻¹ and standard deviation of 2.7 km s⁻¹.

Open in new tabDownload slide

Similar to Fig. 5, we analyse the properties of the star groups in this comparison, however in general we find similar conclusions to those of Section 3.1. The mass of the star groups increases with the degree of grouping (which is deduced from the star mass fraction f) increases. Similarly, the group age spread shows no clear correlation with the models. The virial radius of the groups ranges from ∼10⁻¹ to 10⁰ pc, reflecting the length scales of star-forming regions seen in Fig. 7, while the velocity dispersion of the star groups is ≈6 km s⁻¹, arising from the collision and initial turbulence.

3.3 Parsec-scale isolated cluster simulation

Fig. 9 shows the IMF of the different models as compared with the original simulation by Jaffa et al. (in preparation), labelled J+ at 5 Myr and at 20 Myr. In this comparison, we do not rerun the original simulation for practical reasons explained in Section 2.3.3. If we were to rerun the simulation, stars would form dynamically and a fraction of sink mass is constantly converted to stars, leaving the sinks smaller than the sinks in J+ at any time. This means that we would expect more low mass stars than the IMFs shown in Fig. 9 if stars are formed dynamically. None the less, at 5 Myr, the maximum star mass achievable by the models with less to no grouping, e.g. the ‘no grouping’ (orange) and ‘standard’ (green) models are generally smaller than those with greater degree of grouping, the extreme case being the ‘all grouping’ model (pink). Once again, this is simply because greater grouping allows for more massive mass reservoir to form higher mass stars. Similar to our inference from Section 3.1, the ‘turbulent’ model (red), i.e. setting |$d_{\text{g}} = 1$| pc, |$v_{\text{g}}=3$| km s⁻¹ (about the turbulence speed), and |$\tau _{\text{g}} = t_{\text{ff}} = 5.24$| Myr is able to form high mass stars. At 20 Myr, the maximum star mass achievable by all models converges towards 100 M_⊙, the upper star mass limit of the Kroupa IMF. For the solid lines in Fig. 9, 5 Myr |$\lt t_{\text{ff}} = 5.24$| Myr, so the age criterion is not used at this point.

Figure 9.

The IMF of our models are compared with that of the original simulation by Jaffa et al. (in preparation), labelled as J+, at 5 Myr (solid lines with circle ends) and at 20 Myr (dashed lines with triangle ends). The markers are the maximum mass bin in the respective models. The errorbars are the standard deviations from Poisson sampling. The black dashed line is the analytical Kroupa IMF.

Open in new tabDownload slide

Fig. 10 shows the mass and age spreads of the groups at 20 Myr, where the age criteria do become relevant. In the top subplot, even though the median star group mass of the ‘turbulent’ model is comparable to other models that do not group all sinks, the upper limit for the star group mass is greater, allowing the formation of larger mass stars due to the larger mass reservoir. The ‘no age check’ model, which differs from the ‘standard’ model in its grouping age parameter |$\tau _{\text{g}}$|⁠, shows a similar number of groups (Table 1) and group mass distribution to the ‘standard’ model. This suggests either the grouping age criteria generally do not affect the degree of grouping, or the system needs to evolve much longer to allow the age criteria to take effect. The ‘free-fall’ model, which differs just with the velocity criterion, shows a greater median mass than the ‘no age check’ and ‘standard’ models, showing that the grouping distance and speed parameters affect the grouping of sinks more significantly than the grouping age parameter.

Figure 10.

The mass and age spread of the groups for each model at 20 Myr. For each box, the orange line is the median group mass, the length of the box is the interquartile range, while the whiskers span the whole mass range.

Open in new tabDownload slide

The bottom subplot of Fig. 10 shows the age spread of the sink groups prior to star formation at 20 Myr. Our modified star formation scheme in this comparison does not form stars dynamically, so all stars form at 20 Myr with zero age spread. Similar to Fig. 5, there is no age spread for the sink groups in the ‘no grouping’ model simply because each group consists of only one sink by construction, and also the age spread of the sink group in ‘all grouping’ model is simply the time between first sink formation and the current snapshot of 20 Myr. For the ‘standard’ model, only 32 out of 160 (20 per cent) sink groups have non-zero age spreads, and among them only three groups have age spread |$\gt t_{\text{ff}}/2 = 2.62$| Myr. The age spread distribution of the sink groups in the ‘turbulent’ model is more skewed towards greater values. In both models, the maximum age spread is about |$t_{\text{ff}}=5.24$| Myr, the value of the grouping age parameter |$\tau _{\text{g}}$|⁠. For the ‘no age check’ model, 36 out of 155 (23 per cent) sink groups have non-zero age spreads, similar to the ‘standard’ model. However, because the grouping age criteria is neglected, the maximum age spread can go beyond |$t_{\text{ff}}$|⁠, and in this case five sink groups have age spread beyond |$t_{\text{ff}}$|⁠, the maximum value being 8.49 |${\approx}1.6 t_{\text{ff}}$|⁠. This is undesirable because each small star-forming region (i.e. group) is expected to form stars within one free-fall time-scale, i.e. the formation time-scale under gravitational collapse. Lastly, the ‘free-fall’ model creates three groups at 20 Myr, all of which have age spread around |$t_{\text{ff}} = 5.24$| Myr as expected since sinks can be grouped freely in position and velocity but restricted temporally, so this model can create at most three groups (⌊20/5.24⌋ = 3). In conclusion, even though the grouping age criteria does not seem to alter the group masses much, it is useful to ensure consistent age spreads among the groups. This criteria is essential in star-forming simulations that run much longer than the free-fall or any formation time-scales.

3.4 Kiloparsec-scale spiral arm simulation

Fig. 11 shows the IMF of the stars from a simulation by Rieder et al. (2021) at 1.80 Myr and at 2.40 Myr. Even though the stars are formed using the single star formation method, i.e. no group is assigned, the IMFs are already close to a complete Kroupa IMF. In this low mass resolution simulation, the mass reservoir for star formation, i.e. the average sink mass in this case, is ∼10² M_⊙, enough to sample a complete IMF. This confirms that if the average sink mass is large, then we can fully sample the IMF, which would tend to be the case for larger galaxy or even cosmological simulations. In these simulations, the single star formation method is appropriate.

Figure 11.

The IMF of the sub-galactic simulation by Rieder et al. (2021) at 1.80 Myr and 2.40 Myr are shown. The black dashed line is Kroupa IMF. The errorbars are the standard deviations from Poisson sampling. The black dashed line is the analytical Kroupa IMF.

Open in new tabDownload slide

4 DISCUSSION

4.1 Effect of changing the random seed on cluster properties

In our simulations, the stars’ initial positions are allocated within the accretion radius and the velocities according to the gas dispersions, but we do use a random seed to assign individual positions and velocities. Here we determine whether properties of the clusters formed have any dependence on the random positions and velocities. To test the effect of randomness on our results, we choose the L20Turbulent model (Section 3.2) and perform four more simulations with the exact same initial conditions, except we set a different random seed for star formation compared to the original ‘turbulent’ model. This leads to different positions of the stars, and a different velocity field. Fig. 12 shows the IMF of the ‘turbulent’ colliding clouds model and the reruns at 1.75 Myr. Even though there are slight deviations between the IMFs, changing the random seeds for star formation does not alter the overall trend of the mass function. Other results, such as the cumulative velocity distribution and the properties of the groups of the reruns are also similar compared to those of the ‘turbulent’ model. Using the KS tests as described in Section 3.2, we find that the velocity distributions may be considered statistically the same.

Figure 12.

The IMF of the reruns with different random seeds for star formation are compared with that of Liow & Dobbs (2020), labelled as L20, and the ‘turbulent’ model in Section 3.2, labelled as L20Turbulent, at 1.75 Myr, the time when 10 per cent of the gas mass is converted to sink mass in the original simulation. The markers are the maximum mass bin in the respective models. The errorbars are the standard deviations from Poisson sampling. The black dashed line is the analytical Kroupa IMF.

Open in new tabDownload slide

The stellar distributions of the reruns are visually indistinguishable from that of the ‘turbulent’ model (middle plot in Fig. 7), but we perform additional quantitative analysis to show the similarity between the stellar distributions. Fig. 13 shows the density histograms of the separation of stars for the reruns and the ‘turbulent’ model at 1.75 Myr, similar to the analysis by Torniamenti et al. (2021) to study the spatial distribution of the stars. Subtle differences can be observed among the five models, but generally they share similar overall trends and features such as the peaks at about 1 pc and 3 pc, signifying similar length scales of clustering. We use the same density-based clustering algorithm (DBSCAN; Ester et al. 1996) and parameters (i.e. maximum separation between members of 0.5 pc and minimum number of members of 5) used in Liow & Dobbs (2020) to identify the most massive cluster in each of the ‘turbulent’ model and the reruns. These clusters are shown in Fig. 14 which shows that by eye the distribution of stars and grouping into clusters is very similar. Except for the second rerun (labelled as L20Extra2), the clustering algorithm identifies the same most massive cluster in the reruns and in the ‘turbulent’ model. The most massive cluster in the second rerun is smaller simply because its stellar distribution is slightly more fragmented compared to the other models by chance. By increasing the maximum separation between members slightly to 0.6 pc, the most massive cluster identified in the second rerun becomes the same as the rest, showing that these clusters are statistically indifferent. Therefore, we find that cluster formation and evolution are relatively independent of the stochasticity introduced in our star formation prescription.

Figure 13.

The density histograms of the separation of stars for the reruns are compared with that of the ‘turbulent’ model in Section 3.2 at 1.75 Myr.

Open in new tabDownload slide

Figure 14.

The stellar distribution of the simulation by Liow & Dobbs (2020), labelled as L20, the ‘turbulent’ model in Section 3.2, labelled as L20Turbulent, and the reruns at 1.75 Myr are shown with the most massive clusters (red) identified by DBSCAN.

Open in new tabDownload slide

4.2 Varying the upper star mass limit

So far in this paper, we fix a constant value of 100 M_⊙ as the upper star mass limit to our input Kroupa IMF, which is used to sample the stellar population and decide the mass threshold to form further stars. In the case where the group masses are sufficiently large, most of the mass in the sinks can be converted to stars, for example for our comparisons in Sections 3.2, 3.3, and 3.4. However in Section 3.1, the star mass fraction can be very low, which is clearly not ideal. If we take the extreme case, the ‘no grouping’ model, the average group mass (each group is equivalent to a sink particle in this case) is about 0.02 M_⊙ < 0.37 M_⊙, the expected mass value calculated using the input IMF and mass range. Probabilistically, even relatively low mass stars cannot be assigned to the sinks, and consequently much of the stellar mass remains ‘trapped’ in the sinks rather than being converted to stars.

An alternative is to set the upper star mass limit equal to the individual group mass. Besides eliminating a parameter from our star formation prescription, this new setting ensures that the mass threshold and the star mass that can form within the group are always less than the group mass, avoiding the problem of having the mass threshold and stellar masses larger than the mass available locally. Consequently, the sink-to-star conversion is more efficient as shown in Table 2, where we test this setting by rerunning the models that compare with B12 (Section 3.1) and L20 (Section 3.2). The star mass fraction of the reruns in B12, especially those with lower degree of grouping, is boosted using the new setting. Using the group mass as the upper limit of the IMF has less impact on the star mass fraction of the reruns of L20, as the sink-to-star conversion is already relatively efficient in those models.

One major disadvantage of using the group mass as the upper mass limit is that more low mass stars, and subsequently less high mass stars, are created since the small individual upper star mass limits prevent the formation of higher mass stars. Figs 15 and 16 show the IMF of the reruns as compared to B12 and L20 using the group mass as the upper star mass limit. In a smaller scale system as shown in Fig. 15, all the reruns have more low mass stars as compared to the models in Fig. 2. The most massive stars in the reruns are also less massive compared to their respective counterparts. In the larger scale L20 system as shown in Fig. 16, all models except the ‘all grouping’ case cannot form massive stars > 5 M_⊙ (c.f. Fig. 6). In the previous models that used 100 M_⊙ as the upper star mass limit, the group masses in these large-scale models are generally more massive, but note that these sink groups are usually small upon creation. Therefore, taking the group mass as upper star mass limit instead, these sink groups would form stars right after sink creation and become less massive. Since the density of the sinks is kept constant, their accretion radii become shorter and therefore prevent the groups to grow in mass and form higher mass stars. Potentially, this could be alleviated by including a delay after sinks form before converting them to stars, so the sinks have more time to accrete gas and grow but we do not consider this further here. In summary, setting the upper star mass limit equal to the group mass increases the sink-to-star conversion, but the resultant mass distribution is heavily skewed towards the lower mass end.

Figure 15.

The IMF of the reruns when the upper star mass limit equals to individual group mass are compared with that of Bate (2012), labelled as B12, at the free-fall time of 0.19 Myr. The markers are the maximum mass bin in the respective models. The errorbars are the standard deviations from Poisson sampling. The black dashed line is the analytical Kroupa IMF.

Open in new tabDownload slide

Figure 16.

The IMF of the reruns when the upper star mass limit equals to individual group mass are compared with that of Liow & Dobbs (2020), labelled as L20, at 1.75 Myr, the time when 10 per cent of the gas mass is converted to sink mass in the original simulation. The markers are the maximum mass bin in the respective models. The errorbars are the standard deviations from Poisson sampling. The black dashed line is the analytical Kroupa IMF.

Open in new tabDownload slide

5 CONCLUSION

Sink particles are used in star-forming simulations of various length scales and resolutions to replace dense gas regions, however a method to convert sinks to stars is needed if the sinks are not well resolved as individual stars, which is essential in simulations that require stellar properties like stellar feedback simulations. Wall et al. (2019) introduced the single star formation prescription, where a population of stars that is sampled from the input IMF is introduced for each sink. This method works well with cluster sinks, but undersamples high mass stars if used on smaller mass sinks that are usually approximated as small groups of stars. In this paper, we introduce grouped star formation, a modification of the single star formation method whereby sinks are first grouped according to their positions, velocities, and ages, then the group masses are used to sample the IMF and form stars. Using EKSTER (Rieder et al. 2021), we test this method in simulations of various physical scales, from a sub-parsec isolated cloud simulation up to a kiloparsec spiral arm simulation. We show that the grouped star formation prescription is robust in simulations of different physical scales, and is essential in parsec-scale or smaller scale simulations, whereby their typical mass resolution of 10⁻⁴ − 10⁻¹ M_⊙ means that each sink is only approximated as a small star-forming region. This method would allow us to study the evolution of star clusters more accurately.

One of the main advantages of increasing the degree of grouping of sinks is that the IMF is more complete, i.e. higher mass stars are sampled, as the degree of grouping increases. The disadvantage of increasing the grouping is that local mass may be less well conserved, e.g. stars are placed more preferentially at the few larger sinks rather than spread out in the simulation. Furthermore, stars formed at different times could be grouped together but this is easily constrained by our age criterion. For the models we present here, generally the ‘turbulent’ grouping scenario is optimal. For this we set the grouping distance |$d_{\text{g}} = 1$| pc, grouping speed |$v_{\text{g}}$| equal to the turbulence of the system, and grouping age |$\tau _{\text{g}}$| equal to the free-fall time-scale. This tends to sample the IMF well, group sinks which form within a free-fall time, and reproduces most of the expected stellar mass in stars. This model is also more physically motivated, since the values of these parameters are approximately the physical scales expected in a typical star-forming region.

For our smallest scale region (∼0.5 pc), the ‘all grouping’ method which groups all sinks as one produces the optimal results. This is not really surprising since this simulation is simply of one small star-forming region. We also find that we can produce a closer match to the velocity distribution of the stars formed in Bate (2012) by not including softening for the star particles, which again is not that surprising. With our grouping method we can basically reproduce the IMF, velocity characteristic of the stars, whilst the detailed spatial distribution of the stars is not that important since it is unlikely to be resolved in larger scale simulations. Our ‘turbulent’ grouping model also works well, in this case this model is not very different to the ‘all grouping’ case – again we would expect the ‘turbulent’ grouping to group most of the sinks together if the ‘turbulent’ method is tuned to select individual star-forming regions (where stars located together form together on a free-fall time). The ‘no grouping’ model i.e. single star formation prescription, is a particularly poor choice since only a small mass of the sinks are converted to stars, and there is an absence of massive stars.

With low mass groups, which may arise in particularly small scale simulations, mass is not necessarily efficiently converted from sinks to stars, which is not satisfactory if we are treating the stars as the stellar component of the simulation. We explored changing the upper limit of the IMF from a fixed value to the individual group masses, so that the fraction of mass in stars is increased, but the latter method also has the effect of suppressing the stellar masses such that low mass stars are overpopulated and high mass stars underpopulated. Using the ‘all grouping’ or ‘turbulent’ modes seem to be much better choices than the ‘no grouping’ model here in any case.

In our intermediate scale simulations (∼10 pc), again the ‘all grouping’ model samples the IMF well, but we note that on these scales it is not necessarily appropriate, or desirable to group all the sinks together. Again the ‘turbulent’ grouping still samples the IMF reasonably well, and better groups together sinks based on location, and on these scales, age is also relevant. We note that potentially it may not actually be optimal to sample small clusters of stars according to a full IMF, as such regions may be less likely to contain massive stars, but we do not consider this further here. Again the ‘no grouping’ case does not produce the IMF well, and has a lower fraction of mass in stars compared to the other cases.

Finally we checked our grouping method for a larger scale simulation (∼1 kpc), and verified that the ‘no grouping’ case, i.e. single star formation is sufficient and there is no need to adopt grouping of sink particles once the mass of the sinks is sufficient to sample the IMF.

The type of IMF is a free parameter in our simulations. We choose the Kroupa IMF (Kroupa 2001) for convenience, but one can also choose other types of IMF to suit the problem in hand, e.g the simple broken power-law IMF by Salpeter (1955) or Miller & Scalo (1979), the IMF by Chabrier (2003) which models better the stars at lower mass end, and the IMF from the simulation by Susa, Hasegawa & Tominaga (2014), also investigated in Hirai et al. (2021), to form stars in the range of 1–300 M_⊙ in cosmological simulations. We can only compare the resultant IMF of the stars with our selected choice of IMF.

Our algorithm is simple and easy to add in any star formation code, but it introduces three additional parameters to the simulation. We briefly explored other implementations like using the virial theorem to group the sinks, or reassigning group indices at every time-step, but found these algorithms are unable to group the sinks effectively. Other recently developed clustering methods (e.g. Torniamenti et al. 2021) can potentially be adapted into our grouped star formation prescription to be applied in hydrodynamical simulations. One potential improvement which we did not test further is allowing a delay, so that sink groups can accumulate mass before being converted to stars, although this has the disadvantage of introducing a further free parameter. Essentially though, our method can be easily refined to suit the needs of the user.

ACKNOWLEDGEMENTS

The authors thank Tim Naylor for the helpful discussion on statistics, and Nicholas Herrington for the initial discussion on the analysis of B12 simulation comparison. SR acknowledges funding from STFC Consolidated Grant ST/R000395/1. CLD acknowledges funding from the European Research Council for the Horizon 2020 ERC consolidator grant project ICYBOB, grant number 818940. SEJ acknowledges support from the STFC grant ST/R000905/1 Simulations in this paper were performed on ISCA, University of Exeter. The column density plot of B12 cluster was produced using SPLASH (Price 2007), while those of cloud–cloud collision simulations were produced with the help of FI (Pelupessy, van der Werf & Icke 2004) in AMUSE. The underlying gas velocity fields used to check the stellar velocity were computed using GADGET2 (Springel 2005) in AMUSE. PYTHON packages, such as NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), PANDAS (Wes McKinney 2010), and SCIPY (Virtanen et al. 2020) were used generously to analyse the results in this paper.

DATA AVAILABILITY

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

REFERENCES