# The dynamical state of stellar structure in star-forming regions

J. M. Diederik Kruijssen, Thomas Maschberger, Nickolas Moeckel, Cathie J.

Clarke, Nate Bastian, Ian A. Bonnell

► To cite this version:

J. M. Diederik Kruijssen, Thomas Maschberger, Nickolas Moeckel, Cathie J. Clarke, Nate Bastian, et al.. The dynamical state of stellar structure in star-forming regions. Monthly Notices of the Royal Astronomical Society, 2012, 419, pp.841-853. 10.1111/j.1365-2966.2011.19748.x. insu-03612454

HAL Id: insu-03612454

https://insu.hal.science/insu-03612454

Submitted on 18 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The dynamical state of stellar structure in star-forming regions

J. M. Diederik Kruijssen,<fnref n="1" />,<fnref n="2" />,<fnref n="3" />* Thomas Maschberger,<fnref n="4" />,<fnref n="5" />,<fnref n="6" /> Nickolas Moeckel,<fnref n="4" />
Cathie J. Clarke,<fnref n="4" /> Nate Bastian<fnref n="4" />,<fnref n="7" /> and Ian A. Bonnell<fnref n="8" />

Accepted 2011 September 1. Received 2011 August 17; in original form 2011 April 19

## ABSTRACT

The fraction of star formation that results in bound star clusters is influenced by the density spectrum in which stars are formed and by the response of the stellar structure to gas expulsion. We analyse hydrodynamical simulations of turbulent fragmentation in star-forming regions to assess the dynamical properties of the resulting population of stars and (sub)clusters. Stellar subclusters are identified using a minimum spanning tree algorithm. When considering only the gravitational potential of the stars and ignoring the gas, we find that the identified subclusters are close to virial equilibrium (the typical virial ratio $Q_{vir} \approx 0.59$, where virial equilibrium would be $Q_{vir} \sim 0.5$). This virial state is a consequence of the low gas fractions within the subclusters, caused by the accretion of gas on to the stars and the accretion-induced shrinkage of the subclusters. Because the subclusters are gas poor, up to a length-scale of 0.1–0.2 pc at the end of the simulation, they are only weakly affected by gas expulsion. The fraction of subclusters that reaches the high density required to evolve to a gas-poor state increases with the density of the star-forming region. We extend this argument to star cluster scales, and suggest that the absence of gas indicates that the early disruption of star clusters due to gas expulsion (infant mortality) plays a smaller role than anticipated, and is potentially restricted to star-forming regions with low ambient gas densities. We propose that in dense star-forming regions, the tidal shocking of young star clusters by the surrounding gas clouds could be responsible for the early disruption. This ‘cruel cradle effect’ would work in addition to disruption by gas expulsion. We suggest possible methods to quantify the relative contributions of both mechanisms.

**Key words:** stars: formation – stars: kinematics and dynamics – open clusters and associations: general – galaxies: star clusters: general – galaxies: stellar content.

## 1 INTRODUCTION

Over the past years, the implications of clustered star formation have touched a range of astrophysical disciplines, from the scales of the star formation process itself (see the review by McKee & Ostriker 2007) to the fundamental properties of young star clusters (e.g. McMillan, Vesperini & Portegies Zwart 2007; Allison et al. 2009a; Moeckel & Bonnell 2009), or possibly even the global stellar mass assembly of galaxies (see e.g. Pflamm-Altenburg, Weidner & Kroupa 2007; Bastian, Covey & Meyer 2010). While it seems

evident that most stars form in a clustered setting (e.g. Parker & Goodwin 2007), estimations of the exact fraction are hampered by the substantial dissociation of stellar structure that occurs during (but is not necessarily related to) the transition from the gas-embedded phase to classical, gas-poor star clusters (Lada & Lada 2003; Portegies Zwart, McMillan & Gieles 2010). The traditional interpretation that most, if not all stars form in clusters, with gas expulsion leading to their early disruption ('infant mortality'; see Lada & Lada 2003; Bastian & Goodwin 2006; Goodwin & Bastian 2006) has recently been challenged by observational studies suggesting that stars form with a continuous distribution of densities, of which only the high-density tail eventually leads to bound stellar clusters (Bressert et al. 2010; Gieles & Portegies Zwart 2011).

<footnote n="1">¹Astronomical Institute, Utrecht University, PO Box 80000, 3508 TA Utrecht, the Netherlands</footnote>

<footnote n="2">²Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, the Netherlands</footnote>

<footnote n="3">³Max-Planck Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85748 Garching, Germany</footnote>

<footnote n="4">⁴Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA</footnote>

<footnote n="5">⁵Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany</footnote>

<footnote n="6">⁶UJF-Grenoble 1/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble, UMR 5274, Grenoble 38041, France</footnote>

<footnote n="7">⁷Excellence Cluster Universe, Boltzmannstraße 2, 85748 Garching, Germany</footnote>

<footnote n="8">⁸Scottish Universities Physics Alliance (SUPA), School of Physics and Astronomy, University of St. Andrews, North Haugh, St. Andrews, Fife KY16 9SS</footnote>

<footnote>*E-mail: kruijssen@mpa-garching.mpg.de</footnote>
Current advancements in numerical calculations of turbulent fragmentation in star-forming regions enable the study of clustered star formation in increasing detail (e.g. Bonnell, Bate & Zinnecker 1998; Klessen & Burkert 2000; Bate, Bonnell & Bromm 2003; Bonnell, Clark & Bate 2008). However, theoretical investigations of the response of stellar structure to gas expulsion are still largely based on the assumption of either a static gas potential<fnref n="1" /> or initial equilibrium between the stars and gas (e.g. Tutukov 1978; Adams 2000; Geyer & Burkert 2001; Boily & Kroupa 2003a,b; Baumgardt & Kroupa 2007; Parmentier et al. 2008), which need not apply to star-forming regions in nature. A more realistic setting was recently explored by Offner, Hansen & Krumholz (2009), who find that the velocity dispersions of the stars in hydrodynamic simulations of star formation are smaller than that of the gas by about a factor of 5, suggesting that the assumption of equilibrium between both components indeed does not hold. The response of star clusters to gas expulsion has also been investigated by Moeckel & Bate (2010), who consider N-body simulations of star clusters using initial conditions from hydrodynamic simulations, and by Moeckel & Clarke (2011), who address the dynamical evolution of star clusters under the condition of ongoing gas accretion. They propose that the disruptive effect of gas expulsion is limited by the way in which gas and stars are redistributed by the accretion-induced shrinkage of clusters.

The hydrodynamical calculations of Bonnell et al. (2008) cover the hierarchical formation of several stellar (sub)clusters, which have been identified and analysed by Maschberger et al. (2010). The simulation is very suitable for investigating the properties of the (sub)cluster population due to the relatively large range of mass scales of the modelled structure (see Section 2). In this paper, we analyse the simulations reported in Bonnell et al. (2008) to probe the response of gas-embedded stellar structure to gas expulsion. We consider the dynamical state of the stars while ignoring the gas, which is equivalent to observing the stellar structure under the condition of instantaneous gas expulsion at any time in the simulation.

This paper starts with a discussion of the set-up of the simulations, the subcluster identification algorithm and the characteristics of the stellar structure in Section 2. An analysis of the dynamical state of the subclusters follows in Section 3, which covers quantities such as the virial ratio, bound mass fraction and gas content. The response of the clusters to gas expulsion and an extension to the lengthscales of actual star clusters is considered in Section 4. The paper is concluded with a summary and an outlook in Section 5, where we discuss the possible dependence of the results on the initial conditions of the simulations and the input physics, and suggest ways in which our analysis could be improved and extended.

## 2 SIMULATIONS AND CLUSTER SELECTION

In this work, we analyse the hydrodynamical/N-body simulations of Bonnell, Bate & Vine (2003) and Bonnell et al. (2008), extending the analysis of Maschberger et al. (2010) and Maschberger & Clarke (2011). These smoothed particle hydrodynamics (SPH) simulations follow the evolution of a initially marginally unbound, homogeneous gas sphere of $10^3 M_\odot$ with a diameter of 1 pc (Bonnell et al. 2003), and a cylinder of $3 \times 10^4 M_\odot$ gas, bound in the upper part and unbound in the lower (Bonnell et al. 2008). Initial turbulent motions are modelled with an initially

divergence-free random Gaussian velocity field with a power spectrum $P(k) \propto k^{-4}$. The gas is kept at a temperature of 10 K, staying isothermal throughout the $10^3 M_\odot$ simulation. In the $10^4 M_\odot$ simulation, the gas follows a modified Larson-type equation of state (Larson 2005) composed of three barotropic equations of state:

$$P \propto \rho^\gamma, \quad (1)$$

with P the pressure, ρ the density, and where

$$\begin{align*} \gamma &= 0.75, & \rho \le \rho_1, \\ \gamma &= 1.0, & \rho_1 < \rho \le \rho_2, \\ \gamma &= 1.4, & \rho_2 < \rho \le \rho_3, \\ \gamma &= 1.0, & \rho > \rho_3, \end{align*} \quad (2)$$

and $\rho_1 = 5.5 \times 10^{-19} \text{ g cm}^{-3}$, $\rho_2 = 5.5 \times 10^{-15} \text{ g cm}^{-3}$ and $\rho_3 = 2 \times 10^{-13} \text{ g cm}^{-3}$. For reference, a density of $1 M_\odot (0.01 \text{ pc})^{-3}$ equals $1.6 \times 10^{-17} \text{ g cm}^{-3}$.

Star formation is modelled via sink particles, which are formed if the densest gas particle and its ~50 neighbours are gravitationally bound (the critical density for sink particle formation is $1.5 \times 10^{-15} \text{ g cm}^{-3}$ for the $10^3 M_\odot$ simulation, and $6.8 \times 10^{-14} \text{ g cm}^{-3}$ for the $10^4 M_\odot$ simulation; see Bonnell et al. 2008 for details). The mass resolutions of the sink particles are ~0.1 and 0.0167 $M_\odot$, respectively. Accretion on to sink particles occurs if SPH particles move within the sink radius (200 au for both simulations) and are gravitationally bound, or if SPH particles move within the accretion radius (40 au for both simulations). Gravitational forces between sink particles are softened at smoothing lengths of 160 au ($10^3 M_\odot$) and 40 au ($10^4 M_\odot$).

Under the influence of gravity, the initially smooth gas distributions quickly form filaments in which the sink particles are formed. The sink particles themselves group together in subclusters that merge into larger structures, leading to the formation of one ‘star cluster’ in the $10^3 M_\odot$ simulation and about three ‘star clusters’ in the $10^4 M_\odot$ simulation. Throughout this paper, we will focus on the $10^4 M_\odot$ simulation, which contains about 10 times more subclusters than the $10^3 M_\odot$ simulation, and therefore enables us to consider a population of subclusters rather than a selected set of examples. We also ran our analyses for the $10^3 M_\odot$ simulation, which gave results that are consistent with those from the $10^4 M_\odot$ simulation.

For the identification of the subclusters, we employ a minimum spanning tree (MST)-based clustering technique, which has been used in the context of young star-forming regions (e.g. Maschberger et al. 2010; Kirk & Myers 2011). The MST, which has the advantageous property of not imposing any geometrical symmetry on the data set, has also been used to quantify the amount of substructure (e.g. Cartwright & Whitworth 2004; Schmeja, Kumar & Ferreira 2008; Maschberger et al. 2010) and mass segregation (e.g. Allison et al. 2009a,b; Maschberger et al. 2010; Parker et al. 2011 – see Moeckel & Bonnell 2009; Maschberger & Clarke 2011 for alternative methods) in star-forming regions. The MST is a concept from graph theory, which represents the unique connection of all points of a data set, so that there are no closed loops (a ‘tree’), and so that the total length of all edges between the points is minimal. Typically, two separated groups of points are connected with one single, long edge, whereas the points within the groups have short edges. By simply removing edges that are longer than a chosen break distance the tree can be split in subtrees, which connect the points of the subclusters in the data set (further information on MST-based clustering techniques can be found in Zahn 1971). To avoid

<footnote n="1">¹ Except for a normalization of the gas potential that decreases with time when the gas is expelled.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Spatial distribution of sink particles that are present at the end of the 10⁴ M☉ simulation (t = 0.641 Myr), projected on the x–y plane. Black particles constitute subclusters, and dark grey particles belong to the field population. Since the spatial extent of the simulation in the z-direction is larger than in the x–y plane, some of the apparent clustering is the result of the projection.</figcaption></figure>

spurious detections, we require that a subcluster contains a minimum number of 12 stars.

The MST technique utilizes one free parameter, which is the break distance. Automated methods to determine the break distance are ill-suited for the analysis of the simulations due to the highly varying properties of the stellar distribution. We therefore choose a break distance of $d_{break} = 0.035$ pc, which gives subclusters that are comparable to those identified by the human eye, although they do not include the stellar haloes surrounding them. This break distance is larger than the choice of Maschberger et al. (2010), but leads to comparable clusters because we analyse the stellar structure in three spatial dimensions and not in projection. Because the choice of a single break distance could introduce an artificial length-scale into the analysis (Bastian et al. 2007), we have also performed our calculations for a set of other break distances in the range $d_{break} = 0.020–0.100$ pc, of which the results are used when discussing the implications and applicability of our findings in Sections 4 and 5.

An example of the results of the subcluster identification method can be seen in Fig. 1, which shows the spatial distribution of sink particles and subclusters at the end of the 10⁴ M☉ simulation at t = 0.641 Myr, some 0.3 Myr after the onset of star formation. At this point in the simulation, after one free-fall time, about 60 per cent of all stellar mass is constituted by subclusters. The spatial distribution shown in Fig. 1 is the result of a complex tree of hierarchical merging of small subclusters (Maschberger et al. 2010). This process is still ongoing at the end of the simulation, which is illustrated by the close proximity of the subclusters towards the right in the plane of Fig. 1.

Following Casertano & Hut (1985), we use the projected distance to the Nth nearest neighbour to determine the surface density distribution of the sink particles. For a rank number N, the local surface density at the locations of each of the sink particles is $\Sigma_{sink} = (N - 1)/(\pi D_N^2)$, with $D_N$ the projected distance to the Nth nearest neighbour in the x–y plane. The resulting distribution of surface densities for N = 7 is shown in Fig. 2, which includes all sink particles at the end of the simulation (t = 0.641 Myr), as well as those from a snapshot at t = 0.442 Myr, not too long after the

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Surface density distributions of sink particles. The solid histogram includes all sink particles present at the end of the simulation (t = 0.641 Myr), with the shaded histogram marking the subset of particles that belong to subclusters. The dashed curve represents a lognormal fit to the distribution with a peak at $log(\Sigma_{sink}) = 3.75$ and a dispersion of $\sigma_{log \Sigma} = 1.13$. For comparison, the dotted histogram denotes the surface density distribution of sink particles at t = 0.442 Myr.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Top: time evolution of the number of subclusters. Bottom: time evolution of the mean subcluster mass. The SFE of the entire simulation is indicated along the top axis.</figcaption></figure>

onset of star formation (also see Fig. 3). The difference between the surface density distributions at both times shows that the stellar structure in the simulation typically evolves towards higher densities, even though the density range spanned by the sink particles does not change much. As should be expected, the high end of the surface density distribution is occupied by the sink particles that are residing in subclusters (shaded area), reaching densities of more than $10^5$ stars pc⁻². These surface densities are several orders of magnitude higher than those observed in nearby (<500 pc)
star-forming regions by Bressert et al. (2010), which is not surprising for two reasons. First, crowding obstructs the observation of the densest parts of star-forming regions, which are therefore not included in their sample. Secondly, the high densities that are achieved in the simulation are likely the result of the initial conditions, with a mean initial gas mass surface density of over $10^3 \text{ M}_\odot \text{ pc}^{-2}$ in the x-y plane. None the less, the surface densities do compare well to the high-density region in the Orion Nebula cluster (ONC), of which the central surface density coincides with the peak of the distribution in Fig. 2 (Hillenbrand & Hartmann 1998), although the system under consideration here is about one order of magnitude younger than the ONC.

The subcluster assembly history is considered in Fig. 3, which shows how the number of subclusters and the mean subcluster mass evolve as a function of time and the star formation efficiency (SFE, the fraction of the gas that has been converted to stars) of the entire simulation (see Bonnell et al. 2011 for a detailed spatial analysis of the SFE). The number of subclusters initially increases until it reaches a maximum, which occurs when the formation of new concentrations of sink particles is neutralized by the hierarchical merging of other subclusters. This is nicely illustrated by the mean stellar subcluster mass, which steeply increases around $t = 0.55-0.60$ Myr, when the merging of new-formed subclusters causes the number of clusters to decrease. The mean mass keeps increasing until the end of the simulation due to the ongoing accretion of gas and small subclusters on to sink particles and the merging of small subclusters with a few massive ones. This mass increase would eventually be halted on time-scales much longer than the duration of the simulation, when the available gas reservoir is depleted or, more likely, when the further inflow of gas onto the subclusters is obstructed by feedback from supernovae and stellar winds (neither of which are included in the simulation).

# 3 DYNAMICAL STATE OF STELLAR SUBCLUSTERS

Whether or not a gas-embedded stellar structure survives gas expulsion depends on its dynamical state. Excluding the gas from the dynamical analysis is equivalent to observing the stellar structure at the moment of instantaneous gas removal. This represents the extreme case of gas expulsion, since a more gradual expulsion could allow a subcluster to (adiabatically) respond to the potential change and thereby retain a larger number of stars. As a result, the analysis of the dynamical state of solely the stellar component in simulations of star formation provides a lower limit to the retention of stellar structure upon gas removal.

## 3.1 Dynamical quantities

We have followed the evolution of several (dynamical) properties of the identified subclusters over the course of the simulation of Bonnell et al. (2008), such as the stellar mass, the stellar half-mass radius,<fnref n="2" /> the fraction of the subcluster mass that is bound and the virial ratio. The properties of the stellar component are supplemented with information on the gas, including the gas mass fractions within the subclusters.

The gravitational boundedness and virial ratio of subclusters are fundamental measures for their dynamical state. Both quantities

are based on the potential energy and internal kinetic energy of a subcluster. For a sink particle i, the potential energy $V_i$ is defined as

$$V_i = -\sum_{j \neq i} \frac{Gm_i m_j}{r_{ij}}, \qquad (3)$$

where $m_i$ and $m_j$ are the sink particle masses and $r_{ij}$ their mutual distance. The kinetic energy $T_i$ of a sink particle is

$$T_i = \frac{1}{2} m_i |v_i - v_{cl}|^2, \qquad (4)$$

where $v_i$ and $v_{cl}$ are the respective velocity vectors of the sink particle and the centre of mass of the subcluster. A particle is gravitationally bound if $T_i + V_i < 0$. We define the virial ratio $Q_{vir}$ as

$$Q_{vir} = -\frac{2 \sum_i T_i}{\sum_i V_i}. \qquad (5)$$

The factor 2 reflects the correction for counting the potential energy twice for each particle pair when combining equations (3) and (5). A subcluster is in virial equilibrium if $Q_{vir} \sim 0.5$ and gravitationally bound if $Q_{vir} < 1$. Supervirial subclusters have $Q_{vir} > 0.5$, while subvirial subclusters have $Q_{vir} < 0.5$.

It is possible that a single, dynamically hard binary, triple or multiple system dominates the energy of a subcluster. We correct for this by searching the sink particle list for binaries<fnref n="3" /> and replacing them with a single centre of mass particle. We repeat this step two more times, thereby correcting for triples and higher order multiple systems. During the last iteration, the kinetic and potential energies of the subclusters generally remain unchanged, which indicates that a correction for higher order multiples is not required. We quantify the effect of binaries on the observables of interest below.

In previous studies, the SFE has been identified as a key parameter which determines the survival chances of stellar clusters upon gas expulsion (e.g. Goodwin & Bastian 2006). However, a more fundamental critical factor is the dynamical state of the stars when the gas is removed. The effective SFE, eSFE = $1/2Q_{vir}$, was therefore introduced as a measure for the survival probability of stellar structure at the moment of instantaneous gas expulsion (e.g. Verschueren 1990; Goodwin 2009). If the gas and stars were in virial equilibrium before gas expulsion, the eSFE is equivalent to the actual SFE. If they were not in virial equilibrium, the eSFE no longer reflects the actual SFE. In that case, the survival chance of stellar clusters is not related to the actual SFE, but is solely determined by their dynamical state. The eSFE is naturally higher in the identified subclusters than in the simulation as a whole.

## 3.2 Virial ratio

As star formation progresses, the population of subclusters grows in terms of its total mass. The hierarchical merging of the subclusters inhibits the increase of their number and causes it to level off towards the end of the simulation, when the formation of new subclusters is balanced by their accretion on to more massive ones (see Fig. 3 and Maschberger et al. 2010). Another consequence of this hierarchical build-up is that the properties of the subcluster population as a whole are not a direct representation of the evolution of individual subclusters, but also include ‘emergent’ properties of the population due to the interactions between the subclusters. This

<footnote n="3">³ Binaries are selected by identifying a most bound partner for each sink particle. If it exists and the semimajor axis is smaller than 1000 au, it is considered a binary.</footnote>

<footnote n="2">² The subclusters have predominantly small elongations (Maschberger et al. 2010, fig. 10), which enables the use of a half-mass radius.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Time evolution of the mean virial ratio of the stellar subclusters (black solid line), of the mean virial ratio weighted by subcluster mass (grey dash-dotted line) and of the total virial ratio of all sink particles in the simulation (dashed line). The horizontal dotted lines indicate the marginally gravitationally bound case ($Q_{vir} < 1$) and the virialized case ($Q_{vir} = 0.5$). After $t = 0.5$ Myr, about 40–60 per cent of the stars resides in subclusters for the break distance adopted here ($d_{break} = 0.035$ pc).</figcaption></figure>

is also relevant when considering the mean virial ratio of the subclusters as a function of time. Individual subclusters can be formed either subvirially or supervirially with respect to the total potential, and would eventually viralize with the total potential if kept in isolation. Deviations from this trend occur when subclusters merge or accrete smaller stellar aggregates, which temporarily increase the virial ratio of the merger product due to the relative velocity of the progenitors. Another thing to keep in mind is that we only include the stars in our dynamical analysis, implying that the obtained virial ratio is always higher than its actual value by an amount that depends on the gas fraction. This affects the mean virial ratio of the population of subclusters, in which there is a continuous formation of new, gas-rich subclusters, which are typically still supervirial and only reach $Q_{vir} = 0.5$ after some further evolution.

Despite the complex setting of hierarchical star formation, the evolution of the mean virial ratio can be used as a first indication of how the dynamical state of the subcluster population evolves over time. This is shown in Fig. 4, which also includes the time evolution of the virial ratio of the entire simulation. As indicated earlier, we ignore the contribution of the gas to the gravitational potential, in order to assess the dynamical state of the stellar structure under the assumption of instantaneous gas removal. Even without accounting for the gas potential, the population of subclusters evolves to a near-virialized state on a time-scale of only a few tenths of a Myr. This would suggest that the subclusters are typically gas poor on length-scales corresponding to their half-mass radii. The difference between the mean virial ratio and the mean virial ratio weighted by subcluster mass in Fig. 4 indicates that more massive subclusters are typically somewhat closer to virial equilibrium than low-mass subclusters. This is more of a trend than a relation: a simple linear regression of the virial ratio and subcluster mass *M* yields a best fit of $Q_{vir} = 0.86 - 0.16 \log M$, but with scatter larger than the fitted slope. Lastly, it is also shown by Fig. 4 that the entire stellar population in the simulation does not reach virial equilibrium, but does become marginally bound. This occurs because the simulation as a whole has a higher gas fraction than the subclusters (see Section 3.3), which illustrates that the SFE depends on the location and length-scale on which it is computed. The dynamical state of the entire simulation also bears some traces of the initial conditions, covering

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Dependence of the bound mass fraction (ratio of the total mass of the bound sink particles to the subcluster mass) on the virial ratio for the subclusters from all snapshots of the simulation. Each symbol represents a subcluster. The dashed line represents a linear fit to the data. Like in Fig. 4, the dotted lines indicate the marginally gravitationally bound case ($Q_{vir} < 1$) and the virialized case ($Q_{vir} = 0.5$).</figcaption></figure>

a cylinder that contains a bound upper half and an unbound lower half (see Section 2).

The virial ratios of individual subclusters do not show notable correlations with subcluster mass or half-mass radius.<fnref n="4" /> Instead, they depend more strongly on the recent mass evolution of the subclusters. The virial ratio temporarily increases whenever the subcluster mass increases, be it due to the merging with other subclusters or by individual sink particles moving inside the MST break distance. When sink particles move more than a break distance away and the subcluster mass decreases, the virial ratio decreases as well. Both are natural consequences of the inclusion or exclusion of transient substructure in the identification of the subclusters. The same dependence is found when using different MST break distances to identify the subclusters: larger break distances yield more extended subclusters and consequently the mean virial ratio is higher. When set to extreme values ($d_{break} > 0.050$ pc), close passages of subclusters are incorrectly picked up as merger products, causing a spurious increase of the virial ratio. For the largest break distance used in our analysis ($d_{break} = 0.100$ pc), these fluctuations can yield mean virial ratios briefly hitting $\langle Q_{vir} \rangle = 1$, in clear contrast with the result from Fig. 4.

The quantity that most strongly correlates with the virial ratio is the bound mass fraction of the subclusters, i.e. the fraction of their mass that is bound even without accounting for the potential of the gas (see the definition in Section 3.1, a sink particle is bound if $T_i + V_i < 0$). It is shown in Fig. 5 that subclusters with high virial ratios tend to have lower bound mass fractions, albeit with substantial scatter. This is not surprising, because the virial ratio is efficiently increased by fast, unbound sink particles that are included by the MST algorithm but would be left out with a physically motivated identification. For larger break distances, the correlation between bound mass fraction and virial ratio is stronger, due to the erroneous identification of kinematically hot structure as subclusters. None the less, most subclusters contain only very few unbound stars, with typical bound mass fractions of 0.95. Fig. 5 also confirms that most subclusters are close to virial equilibrium, which was already suggested by the evolution of the mean virial ratio in Fig. 4. For

<footnote n="4">⁴ Except for the unbound subclusters ($Q_{vir} > 1$), which are generally small ($r_h < 0.04$ pc) and low mass ($M < 30$ M$_\odot$).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Histogram of the virial ratios of the subclusters from all snapshots of the simulation (solid line). The shaded histogram represents the set of subclusters from the last snapshot at $t = 0.641$ Myr. The dashed line is a Gaussian fit to the data for all snapshots, with mean value $Q_{vir} = 0.59$ and standard deviation $\sigma_Q = 0.16$. The vertical dotted lines again indicate the marginally gravitationally bound case ($Q_{vir} < 1$) and the virialized case ($Q_{vir} = 0.5$).</figcaption></figure>

further illustration, Fig. 6 shows the distribution of the virial ratios of the subclusters from all snapshots, as well as those from the last snapshot, which are shown as the shaded region. Only 8 per cent of the subclusters are unbound when excluding the gas, while 25 per cent remains subvirial. When considering the subclusters from all snapshots, a Gaussian fit to the distribution of virial ratios gives a mean of $Q_{vir} = 0.59$ and a standard deviation of $\sigma_Q = 0.16$. As in Fig. 4, the gradual decrease of the mean virial ratio towards $Q_{vir} = 0.5$ is also visible in Fig. 6. A comparison of the two histograms shows that the subclusters in the last snapshot are closer to being virialized than the population of subclusters from all snapshots. These virial ratios imply that the eSFE is close to unity, i.e. the majority of subclusters will not be strongly affected by gas expulsion (see Section 4).

Replacing hard binaries and higher order multiples by their centre of mass particles is essential to obtain a reliable picture of the subcluster dynamics. The disruption of the subclusters during gas expulsion is controlled by the dynamical state of the binary centres of mass rather than the binaries themselves. Had we not corrected for binaries or higher order multiples, the measures for the dynamical state of the subclusters would fluctuate with the orbital phase of a few tightly bound and eccentric binaries.

The bound mass fraction of the subclusters is not strongly affected by the presence of binaries (unbound sink particles are generally single), but because binaries are in virial equilibrium<fnref n="5" /> or slightly subvirial, the mean virial ratio of the subclusters from all snapshots is decreased by 0.1–0.2 if it is not corrected for multiples. About two-thirds of this difference is due to binaries, while the remaining third is accounted for by triples and quadruples. This shift of the virial ratio means that without correcting for binaries, the subclusters could be incorrectly interpreted as being slightly subvirial<fnref n="6" />

($Q_{vir} \sim 0.45–0.50$), and the entire simulation would be close to virialized ($Q_{vir} \sim 0.60–0.65$) instead of the marginally bound state that is shown in Fig. 4. With respect to the binary-corrected results from Fig. 4, this rather modest difference arises because the finite gravitational smoothing length used in the simulation inhibits the formation of very hard binaries. None the less, the correction for binaries improves the accuracy of our analysis, and therefore all results shown in this paper are corrected for binaries and higher order multiple systems.

### 3.3 Gas content

The key question is why the subclusters are so close to virial equilibrium when neglecting the gas potential. An obvious answer would be that the subclusters are generally gas poor, which would imply that they are hardly affected by the gas potential in the first place. To assess the gas potential and its time evolution, we have analysed the distribution of the gas in two snapshots of the simulation, at times $t_1 = 0.442$ Myr (when star formation is ongoing) and $t_2 = 0.641$ Myr (the last snapshot of the simulation, after one free-fall time; also see Fig. 3). For each of the identified subclusters in these snapshots, we calculate the fraction of the total mass within the stellar half-mass radius of the stellar distribution that is constituted by gas. The distribution of these gas fractions is shown in Fig. 7, which confirms that the subclusters are indeed gas poor on their typical length-scales, with gas fractions of $\langle f_{gas} \rangle = 0–0.2$. Because the simulation does not include feedback, this means that the accretion of gas on to the sink particles can keep up with the overall gas inflow towards the subclusters. Another mechanism that naturally leads to gas-poor subclusters is their accretion-driven shrinking (Bonnell et al. 1998; Moeckel & Bate 2010; Moeckel & Clarke 2011), which increases the density contrast between the subclusters and the surrounding gas.

Gas accretion and the time evolution of the structural properties of the population of subclusters both further decrease the gas fraction as time progresses. This evolution is illustrated by comparing the data of the two snapshots in Fig. 7, corresponding to times $t_1 = 0.442$ Myr and $t_2 = 0.641$ Myr. During the enclosed time interval, the mean gas fraction of all detected subclusters decreases by 0.63 dex, from $\langle f_{gas}(t_1) \rangle = 0.238$ to $\langle f_{gas}(t_2) \rangle = 0.056$. The

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Histogram of the gas fraction within the stellar half-mass radius $r_h$ of each subcluster from the two snapshots at $t_1 = 0.442$ Myr and $t_2 = 0.641$ Myr (solid line). The shaded histogram only shows the gas fractions for the last snapshot at $t = 0.641$ Myr.</figcaption></figure>

<footnote n="5">⁵ Instantaneously, this only holds for binaries on circular orbits. Binaries on eccentric orbits exhibit a variation of the virial ratio, with a subvirial state near apocentre and a supervirial state near pericentre. When considering the time-averaged virial ratio, eccentric binaries are in virial equilibrium. However, because the phase velocity is lowest near apocentre, a sampled population of randomly oriented binaries will typically be slightly subvirial.</footnote>

<footnote n="6">⁶ As in all of Section 3.2, this statement excludes the gravitational potential of the gas.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Subcluster mass-weighted, mean cumulative mass fractions ($\mu(\xi)$) (see equation 6) of sink particles (dashed line), gas particles (dashed-dotted line) and the sum of both (solid line), as a function of the radial distance in units of the half-mass radius ($\xi \equiv r/r_h$). Top: mean cumulative distributions for the subclusters present at $t_1 = 0.442$ Myr. Bottom: mean cumulative distributions for the subclusters present at $t_2 = 0.641$ Myr. The shaded areas enclosed by the dotted lines mark the 16th and 84th percentiles and illustrate the typical spread of the enclosed mass fractions of sink particles (dark grey) and gas (light grey).</figcaption></figure>

mean half-mass radius of the subclusters<fnref n="7" /> decreases from $\langle r_h(t_1) \rangle = 0.020$ pc to $\langle r_h(t_2) \rangle = 0.013$ pc, which is a decrease of 0.18 dex. Even though the shrinking of subclusters is a second-order effect caused by gas accretion, it is interesting to ask which of both mechanisms contributes most to the decrease of the gas fraction. Is it mainly driven by the increasing mean stellar density of the subclusters or by the ongoing gas accretion on to the sink particles?

To assess the relative contributions to gas depletion by accretion and subcluster shrinking, we consider the spatial distribution of the sink particles and the gas. Because of the relatively small numbers of stars in individual subclusters, it is best to examine the mean density profiles of the populations of subclusters in the two snapshots at $t_1$ and $t_2$. Such a combination of the different density profiles decreases the influence of low-number statistics on the result. In Fig. 8, we show the subcluster mass-weighted, mean cumulative mass distributions of gas, sink particles and both combined. The distributions represent the enclosed mass fractions $\mu$, normalized to the sum of the subcluster mass $M_{cl}$ and the enclosed gas mass within three stellar half-mass radii $M_{gas}$:

$$ \mu_i(\xi) = \frac{M_i(\xi)}{M_{cl} + M_{gas}}, \qquad (6) $$

with $i = \{\text{stars, gas, all}\}$ and $M_i(\xi)$ the enclosed mass at $\xi \equiv r/r_h$, which is the radial distance in units of the stellar half-mass

radius. The mean distributions shown in Fig. 8 are weighted by subcluster mass to emphasize those subclusters with better statistics. A first comparison of both panels in Fig. 8 shows that the gas fraction indeed decreases between $t_1$ and $t_2$. The contribution to this decrease by subcluster shrinking can be estimated by a simple thought experiment, in which the gas distribution is kept fixed and the distribution of stellar mass is compressed by the appropriate amount. In the top panel of Fig. 8, the half-mass radii $\langle r_h(t_1) \rangle$ and $\langle r_h(t_2) \rangle$ correspond to $\xi = 1$ and 0.66, between which the enclosed gas mass differs by 0.26 dex. In other words, if the gas distribution was held fixed and the stellar distribution was shrunk appropriately, then the gas fraction within the new half-mass radius would have declined by 0.26 dex. This is a probe for the decrease of the gas fraction that is solely caused by the shrinking of the subclusters. Comparing it with the actual decrease of the mean gas fraction between $t_1$ and $t_2$ of 0.63 dex, we see that it covers about half of the decrease, with the remaining 0.37 dex covered by gas accretion itself – not only by adding to the mass in stars, but also by decreasing the gas mass. We conclude that the evacuation of the gas due to ongoing gas accretion is about equally important for the gas depletion as the shrinking of the subclusters.

Apart from enabling a quantitative comparison of the effect of gas accretion and cluster shrinking, Fig. 8 also demonstrates the spatial variation of the gas fraction in the subclusters. At early times, the gas is still prevalent in the outskirts of the subclusters, contributing 20–60 per cent of the enclosed mass at $\xi = 3$. At the end of the simulation this gas has mostly vanished, leaving only a few per cent of the mass within the stellar half-mass radius as gas, and typically 20 per cent at $\xi = 3$. The radial dependence of the enclosed gas fraction is qualitatively reminiscent of the model of Adams (2000). It is interesting to note that the relative increase of the enclosed gas mass fraction with respect to the enclosed sink particle mass fraction only occurs at radii where the latter flattens, i.e. the subclusters only become gas rich at radii where very little stellar mass is present. The influence of the gas on the subcluster dynamics is therefore best evaluated at radii smaller than where the flattening of $\mu_{stars}$ occurs. At $t_1$, the ratio between the enclosed stellar mass and gas mass just before the flattening is about 4:1, while at $t_2$ it has increased to 11:1. This suggests that if feedback starts at a time $t > t_2$, the resulting gas expulsion will not strongly affect the subcluster dynamics, and that their virialized state (see Section 3.2) will be largely retained.

# 4 RESPONSE TO GAS EXPULSION

Motivated by the low gas fractions found in Section 3, we now address the response of the subclusters to gas expulsion in more detail.

## 4.1 The expansion of subclusters

The long-term response of the subclusters to gas expulsion can be evaluated by once again omitting the gas from the simulations and considering only the identified stellar subclusters and their evolution towards virial equilibrium. Given a certain virial ratio and bound mass fraction, does a subcluster expand or shrink after gas removal? We combine the data from the simulations with a simple energy argument similar to Hills (1980) to estimate how the subcluster masses and half-mass radii evolve after the expulsion of the gas. It is insightful to consider the system at two key moments.

(i) *The time of instantaneous gas removal*, which is equivalent to the current system in the simulations while omitting the gas. This

<footnote n="7">⁷ The subcluster sizes are strongly correlated with the adopted MST break distance. See Section 4.1 for a discussion.</footnote>
can be done for each snapshot, thereby providing a larger sample of subclusters than when only the last snapshot was to be considered. Of course, including subclusters from different snapshots implies a correspondingly extended range of moments of gas expulsion.

(ii) The time at which each subcluster attains virial equilibrium after removing the gas. These times are different for each cluster per definition, but by considering the subclusters at their respective times of virialization the long-term impact of gas removal is most clearly isolated and identified.

The evolution of the subclusters between these two moments can be quantified by evaluating the conservation of energy. For any subcluster, we can express the kinetic energy as

$$T = \frac{1}{2} M_{\mathrm{cl}} \langle v^2 \rangle = -V Q_{\mathrm{vir}} \approx \frac{GM_{\mathrm{cl}}^2}{2r_v} Q_{\mathrm{vir}}, \quad (7)$$

where $r_v$ is the virial radius, and $\langle v^2 \rangle$ denotes the mean square velocity, which as a result can be written as

$$\langle v^2 \rangle = Q_{\mathrm{vir}} \frac{GM_{\mathrm{cl}}}{r_v}. \quad (8)$$

The total energy at the moment of instantaneous gas removal thus becomes

$$E_1 = (Q_{\mathrm{vir},1} - 1) \frac{GM_{\mathrm{cl},1}^2}{2r_{v,1}}, \quad (9)$$

where the relevant quantities have been marked with subscript ‘1’ to indicate the moment of gas expulsion.

Given a deviation from virial equilibrium, a subcluster will respond by changing its radius and/or mass. As can be verified from Fig. 5, most subclusters contain a certain number of unbound sink particles, which were either previously retained by the gas potential, or are randomly passing the subcluster close enough to be included by the cluster identification algorithm. These unbound sink particles will escape the subcluster upon gas expulsion and take away some of the kinetic energy. We now consider a second moment in time, at which the gas-rid subcluster has reached virial equilibrium ($Q_{vir,2} = 0.5$) and the unbound sink particles have successfully escaped. At this time, energy conservation dictates

$$E_1 = E_2 + E_{\mathrm{esc}} = - \frac{GM_{\mathrm{cl},2}^2}{4r_{v,2}} + (M_{\mathrm{cl},1} - M_{\mathrm{cl},2}) \frac{\beta \langle v_2^2 \rangle}{2}, \quad (10)$$

where $E_2$ is the total energy of the (virialized) subcluster, $E_{esc}$ is the total energy of the escaped stars and the relevant quantities have been marked with subscript ‘2’ to indicate the moment of virialization. The parameter $\beta$ denotes the surplus energy per unit mass of the escaped stars after they clear the potential of the subcluster, in units of its mean square velocity. The values of $\beta$ can be estimated from the simulation by computing $(T_i + V_i)/(GM_{cl,2}m_i/2r_{v,1})$ for each of the unbound sink particles.<fnref n="8" /> For this, we use the relation between the virial and half-mass radius corresponding to a Plummer (1911) potential, which is given by $r_v = 1.3r_h$ (e.g. Heggie & Hut 2003). The escaping sink particles are the tail of an approximately Maxwellian velocity distribution of the sink particles in the subcluster, and consequently the distribution of $\beta$ declines exponentially as $f(\beta) \propto \exp(-\beta/\beta_0)$, with $\beta_0$ around unity. The mean of such a distribution equals $\beta_0$ per definition, which illustrates that unbound

stars typically retain velocities similar to the mean square velocity in the subcluster after they escape.

Combining equations (8) and (10), one obtains an expression for the evolution of the gas-rid subcluster as it approaches virial equilibrium, which relates the half-mass radii, masses, initial virial ratio and $\beta$. It is given by

$$\frac{r_{h,2}}{r_{h,1}} \approx \frac{r_{v,2}}{r_{v,1}} = \frac{1}{1 - Q_{\mathrm{vir},1}} \left[ \frac{1 + \beta}{2} \left( \frac{M_{\mathrm{cl},2}}{M_{\mathrm{cl},1}} \right)^2 - \frac{\beta}{2} \frac{M_{\mathrm{cl},2}}{M_{\mathrm{cl},1}} \right], \quad (11)$$

where $r_{h,2}$ is the only unknown and all other variables given by the simulation. For $\beta = 0$, i.e. all escaping stars are only marginally unbound, the expression returns the basic result that when unbound stars escape from a virialized system ($Q_{vir,1} = 0.5$), it contracts to reattain virial equilibrium.<fnref n="9" /> Inserting typical values of $\beta = 1$ (see above), $Q_{vir,1} = 0.59$ and $M_{cl,2}/M_{cl,1} = 0.95$ (see Section 3.2) in equation (11) yields $r_{h,2}/r_{h,1} = 1.04$, which justifies the earlier assumption that the radius does not change much between instantaneous gas removal and virialization (see Footnote 8). This minor expansion is driven by the slightly supervirial state of the subclusters, and inhibited by the escape of unbound stars, which have velocities larger than the escape velocity.

As discussed in Section 3.2, the virial ratios of the subclusters give an indication of their survival fraction after gas removal. Out of all 140 subclusters identified in the simulation, only 10 have virial ratios $Q_{vir} > 1$ and are therefore unbound. In the last snapshot of the simulation, only one of the 21 subclusters will disperse after the removal of the gas.<fnref n="10" /> As a result, typically 90–95 per cent of all the identified subclusters survive gas expulsion. The fate of these survivors depends on whether they expand, and how their environment affects them. Expanded subclusters with lower densities are more susceptible to disruption by tidal shocks. The evolution of the half-mass radii after gas expulsion can be considered in more detail by evaluating equation (11) for each of the subclusters in the simulation. This enables a comparison of the distribution of half-mass radii of the current subclusters (at the moment of instantaneous gas removal) with the distribution of their half-mass radii when they have reached virial equilibrium, which is shown in Fig. 9. The distribution of half-mass radii changes remarkably little after gas removal, as the means of the lognormal functions that are fitted to both distributions differ by 0.035 dex. This implies $r_{h,2}/r_{h,1} = 1.08$, very similar to the earlier, simple estimate of $r_{h,2}/r_{h,1} = 1.04$. The subclusters in the last snapshot experience roughly 1.5 times this expansion after gas removal, which is of the same order of magnitude as the expansion of the other subclusters. As gas expulsion does not increase the cluster radii, they become no more susceptible to tidal perturbations once the gas is removed than they are at birth. However, this does not exclude any expansion of the subclusters, since mass loss due to stellar evolution and tidal perturbations would still cause them to expand (also see Section 4.3).

The characteristic size of the distribution in Fig. 9 strongly correlates with the break distance of the MST. As indicated in Section 2,

<footnote n="8">⁸ The denominator holds a slightly modified form of the mean square velocity $\langle v_2^2 \rangle$ (see equation 8) and assumes that the virial radius does not change much between gas expulsion and virialization (i.e. $r_{v,1} \approx r_{v,2}$). This is required since $r_{v,2}$ is not available in the simulation. The assumption will be verified below.</footnote>

<footnote n="9">⁹ This situation, in which the naturally unbound component of a system escapes, should not be confused with the response of a virialized system to mass loss due to stellar evolution, when energy is injected into the system to unbind mass. In such a case, the surplus energy of the escaping mass is supplied by the energy injection and not by the dynamical system itself, which therefore mainly loses potential energy. This does not apply to the case under consideration in equation (11), where no energy is injected and the unbound stars take away more kinetic energy than potential energy.</footnote>

<footnote n="10">¹⁰ These unbound subclusters are typically low-mass, compact systems, and are often newly formed or have just experienced a subcluster merger.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Histogram of the stellar half-mass radii of the subclusters from all snapshots (solid line). The shaded histogram represents the set of subclusters from the last snapshot at *t* = 0.641 Myr. Top: for the simulated radii, representing the moment of instantaneous gas removal. Bottom: for analytically computed radii, reflecting the subclusters at a later moment, when they have reached virial equilibrium. The dashed lines are lognormal fits to the data for the subclusters from all snapshots, with mean values log (*r*<sub>h</sub>/pc) = {-1.98, -1.95} and standard deviations σlog*r* = {0.73, 0.61} for the top and bottom panel, respectively, with median half-mass radii log (*r*<sub>h</sub>/pc) = {-1.96, -1.91}. The radii depend on the adopted break distance (see text), but the lognormal shape of the distribution persists.</figcaption></figure>

we varied the break distance over the range $d_{\text{break}} = 0.020-0.100$ pc. By fitting lognormal functions to each of the obtained size distributions, we find that the corresponding mean sizes $\hat{r}_h$ correlate with break distance as $\hat{r}_h/\text{pc} = (0.88 \pm 0.16)(d_{\text{break}}/\text{pc})^{1.33\pm0.05}$ for $0.020 \le d_{\text{break}}/\text{pc} \le 0.100$. This demonstrates that the *characteristic* size of the size distribution has no physical meaning. None the less, it is consistent with a lognormal distribution for all break distances, i.e. independently of the definition for the subcluster radius, and the smallest subcluster size does remain constant at about 0.005 pc, corresponding to the smallest cores in the simulation (Smith, Clark & Bonnell 2009). Regardless of the adopted break distance, the smallest subclusters are always those with the lowest numbers of stars. These length-scales are not very surprising, because we are not considering actual stellar clusters, but the compact (sub)systems that develop during the early formation of a cluster.

Indeed a length-scale of a few 0.01 pc is much smaller than that of typical embedded clusters (see e.g. Lada & Lada 2003) and of the order of the maximum size of Class 0 protostars (Launhardt et al. 2010). The formation of multiple stars in such a small volume could be related to the use of sink particles in the simulation, although it should be noted that small stellar clusters and multiple systems within star-forming globules are actually found on scales as small as 0.005–0.015 pc (Kraus & Hillenbrand 2008; Gutermuth et al. 2009; Launhardt et al. 2010). Such separations are consistent with simulated low-*N* subclusters of which the half-mass radii are dominated by a single massive star or binary. In Section 4.3, the units

of the simulation are rescaled to cover the typical mass, length and time-scales of low-mass, young star clusters, causing the half-mass radii of the identified clusters to reach up to 0.3 pc.

## 4.2 The cluster formation efficiency

The instantaneous gas removal discussed in this paper is an extreme form of the more gradual expulsion occurring in nature. As a result, the described weak effect of gas expulsion should be even weaker in real subclusters. It seems that gas expulsion plays a negligible role on the length-scales of the compact stellar aggregates in star-forming regions. However, the regions between subclusters may still be gas dominated, implying that feedback could prevent the further merging of subclusters and thereby inhibit their hierarchical growth.

The length-scale on which the subclusters will have merged and have become gas poor depends on the moment at which feedback starts. For the MST break distance and corresponding length-scale that is used in most of this paper, the subclusters are gas poor irrespective of time. However, there should be a break distance at which a notable time evolution of the gas fraction appears. By comparing the subcluster gas fractions for different break distances, we find that at the end of the simulation (after one free-fall time), the subclusters have become gas poor ($f_{\text{gas}} < 0.1$) on a length-scale of about 0.1–0.2 pc. This length-scale will increase further with the number of free-fall times that are completed before the onset of gas removal. In turn, this increases the spatial extent over which the subclusters are allowed to merge before gas expulsion, which implies that the most massive bound structure is inversely related to the free-fall time.

The free-fall time is related to the density as $t_{\text{ff}} \propto \rho^{-1/2}$, which implies that the time of the onset of gas removal by feedback $t_{\text{fb}}$ is associated with a density $\rho_{\text{fb}}$ that has a free-fall time equal to $t_{\text{fb}}$. For a given density spectrum of subclusters (see e.g. Bressert et al. 2010), only those subclusters with densities $\rho \gg \rho_{\text{fb}} \propto t_{\text{fb}}^{-2}$ have the opportunity to undergo the collapse and shrinkage that we find in the simulations. The cluster formation efficiency (CFE) increases with the fraction of subclusters that forms in these density peaks. As subclusters merge, accrete gas and shrink, the density of the stellar structure further increases (see Section 3). Each free-fall time, more subclusters evolve into the density regime where $\rho \gg \rho_{\text{fb}}$, also on larger length-scales. This means that the length-scales on which star-forming regions produce virialized stellar systems that are insensitive to gas expulsion are larger in dense sites of star formation than in sparse ones. The resulting dense clusters are also less susceptible to disruptive tidal effects from their environment, which potentially further increases their survival chances. As a result, the CFE should increase with density. Through the Schmidt–Kennicutt law (Schmidt 1959; Kennicutt 1998), this suggests a relation between the CFE and the star formation rate (SFR) per unit volume $\rho_{\text{SFR}}$ or per unit surface area $\Sigma_{\text{SFR}}$. Indeed, first observational indications for such a relation have been found by Larsen & Richtler (2000), Larsen (2004) and recently also by Goddard, Bastian & Kennicutt (2010), who obtain CFE $\propto \Sigma_{\text{SFR}}^{0.24}$. A relation between the CFE and the SFR density would also be consistent with the high cluster formation efficiencies that are found in starburst galaxies (e.g. Zepf et al. 1999). However, dense star-forming regions are generally also more disruptive to clustered structure due to the higher frequency and amplitude of tidal shocks (Lamers et al. 2005; Kruijssen et al. 2011). This implies that a relation between the CFE and the ambient density would be weakened or could even be cancelled (also see Section 4.3).
## 4.3 Infant mortality and the ‘cruel cradle effect’

It is often said that most stars form in stellar clusters. The concept of ‘infant mortality’ (Lada & Lada 2003), i.e. the rapid dispersal of stellar structure following the change of the gravitational potential due to gas expulsion, has been put forward in the literature to explain the observed rapid, mass-independent decline of the number of stellar clusters between ages of a few Myr and several tens of Myr (e.g. Bastian & Goodwin 2006). Because the majority of stars has been thought to form in clusters, infant mortality was also held responsible for the low number of clusters per unit field star mass. However, recent (observational) evidence is pointing towards a picture in which star clusters are the dense end of a continuous density spectrum of star formation (see Fig. 2 and Bressert et al. 2010; Gieles & Portegies Zwart 2011). This view challenges the need for infant mortality in the early disruption of stellar structure.

The results presented in this paper show that stellar substructure can evolve towards a virialized state before the gas is removed. This occurs because the dynamics of the stars and the gas are decoupled (also see Offner et al. 2009), as the accretion of gas on to the stars together with the subcluster shrinkage can compensate the overall gas inflow on to the subclusters. In time, this causes the subclusters to become gas poor, thereby diminishing the disruptive effect of gas expulsion. It depends on the length- and mass-scales on which the gas is evacuated whether this result can be extended from subclusters to actual star clusters. Towards the end of the simulation, after about 0.3 Myr of star formation, the subclusters have a mean mass of $40 M_\odot$ and are gas poor on length-scales of 0.1–0.2 pc. As a very crude first-order estimate, one can rescale the units of the simulation to have the subclusters match the typical properties of young star clusters.<fnref n="11" /> By multiplying the mass unit by a factor of 25 and the length unit by a factor of 5, we rescale the mean cluster mass to $10^3 M_\odot$, and the length-scale on which the stellar structure will be gas poor to 0.5–1 pc. By scaling the time unit accordingly, we see that such gas depletion is reached on a time-scale of 0.8 Myr, which is of the same order as the expected $t_{fb}$ due to stellar winds and ionization feedback. This order-of-magnitude estimate is of course far from conclusive, but it does show the relevance of pursuing this problem further.

If clusters reach a relatively gas-poor state before the onset of gas removal by feedback, the influence of gas expulsion on the dynamical state of the clusters will be smaller than previously expected. Rather than leading to the disruption of clusters ('infant mortality'), the different spatial distributions of gas and stars imply that gas expulsion could leave clusters marginally affected, unbinding only a certain fraction of their stars (e.g. Moeckel & Bate 2010). This is in clear contrast with earlier (theoretical) approaches in literature (e.g. Boily & Kroupa 2003a,b; Bastian & Goodwin 2006; Goodwin & Bastian 2006; Baumgardt & Kroupa 2007; Parmentier et al. 2008), which assumed a model where the gas and stars are in equilibrium during gas expulsion. Clearly, this is not the case in the simulation of Bonnell et al. (2008).

Because the density of the stellar structure determines whether or not gas expulsion affects the survival chances of star clusters, a continuous density spectrum of young stellar structure as in Bressert et al. (2010) and Fig. 2 naturally leads to the situation in which the dispersed part of the new-born stellar structure is af-

fected by gas expulsion, while the other, dense and clustered part is not. However, this does not imply that the survival chances of these clusters are necessarily higher. Recent papers have argued that the disruption of star clusters due to tidal shocks from the natural environment could be substantial (Elmegreen & Hunter 2010; Kruijssen et al. 2011; Kruijssen & Bastian, in preparation). Although the disruption rate due to tidal shocks decreases with cluster density, sufficiently strong<fnref n="12" /> shocks would still be able to disrupt dense clusters<fnref n="13" /> (Gieles et al. 2006). Such shocks could be prevalent in dense star-forming regions. As clusters move out of their primordial environment, the typical ambient gas density decreases (Elmegreen & Hunter 2010; Kruijssen et al. 2011), which lessens the disruptive effect of tidal shocks. Observationally, this mechanism affects the star cluster population in a way that is very similar to infant mortality: the fact that young clusters are more efficiently disrupted than older clusters gives rise to a strong decline of the number of clusters with age at young ages. This decline acts on the age scale corresponding to the time it takes to migrate out of the star-forming region. Rather than being an internal effect, like infant mortality is, the primordial disruption by tidal shocks is driven by the state of the environment in which the clusters are born. We will therefore refer to this mechanism as the ‘cruel cradle effect’.

It will be interesting to quantify what the relative contributions of gas expulsion and the cruel cradle effect are to the early disruption of young stellar clusters. It is possible that both effects co-exist, and that the relative importance changes with the environment. It was explained in Section 4.2 that the fraction of clusters that is affected by gas expulsion decreases with the ambient density of the star-forming region. The cruel cradle effect shows the opposite dependence, as the disruptive effect of tidal shocks increases with the ambient gas density. It would therefore not be unlikely that gas expulsion and the cruel cradle effect each dominate a different side of the gas density spectrum of star-forming regions. Their relative strength would then determine the relation between the CFE and the ambient gas density.

# 5 SUMMARY AND OUTLOOK

In this paper, we have assessed the dynamical state of stellar structure in star-forming regions and its response to gas expulsion by analysing the properties of the stellar structure in the SPH/sink particle simulations of Bonnell et al. (2003, 2008). Subclusters have been identified using a MST algorithm (following Maschberger et al. 2010), and binaries have been replaced by their centre-of-mass particles when computing the global dynamical properties of the subclusters. We have also discussed the long-term implications of gas expulsion for the properties of star cluster populations. The main results are as follows.

<footnote n="12">¹² The strength of a tidal shock corresponds to the amount of energy it injects into the cluster.</footnote>

<footnote n="13">¹³ Note that this early tidal disruption of stellar clusters is very different from the mechanism studied by Parmentier & Kroupa (2011), who consider the enhanced loss of stars due to the expansion of a young cluster to its tidal boundary that is induced by gas removal. The difference again lies in their assumption of equilibrium between the stars and the gas. If a young cluster has evacuated the surrounding gas and gas expulsion plays a minor role, then any gas expulsion-induced tidal stripping will be minor as well. None the less, young clusters would still expand due to stellar evolutionary mass loss and the tidal shocks that are causing them to be disrupted.</footnote>

<footnote n="11">¹¹ The simulations are scale-free except for the sink particle radius and accretion radius mentioned in Section 2. For the scaling used in this example, these respective radii are $10^3$ and 200 au.</footnote>
(i) The surface density distribution of sink particles follows an approximately lognormal distribution similar to that observed by Bressert et al. (2010). However, the surface density corresponding to the peak of the distribution is several orders of magnitudes higher than the observed one, because the subclusters considered in our study are part of a region that would represent only one or two clusters in the observations. The high-density end of the distribution is occupied by sink particles belonging to the subclusters that are identified with the MST.

(ii) When excluding the potential of the gas from the dynamical analysis and only considering the sink particles, we find that the simulation as a whole becomes marginally bound after one free-fall time, and the population of individual subclusters is close to virial equilibrium. The mean value of a Gaussian fit to the distribution of virial ratios from all snapshots is $Q_{vir} = 0.59$, where virial equilibrium would imply $Q_{vir} = 0.5$. The mean virial ratio of the population slowly decreases with time, from $Q_{vir} = 0.70-0.80$ early on to $Q_{vir} = 0.55-0.60$ towards the end of the simulation.

(iii) The virialization of the subclusters occurs due to their low gas fractions. We consider the spatial distributions of gas and sink particles at two characteristic moments during the simulation ($t_1 = 0.442$ Myr and $t_2 = 0.641$ Myr, reflecting the system early on and after one free-fall time), and find that the mean gas fraction within the stellar half-mass radii ($r_h \sim 0.01$ pc) of the subclusters decreases by 0.63 dex during the enclosed time interval, from $\langle f_{gas}(t_1) \rangle = 0.238$ to $\langle f_{gas}(t_2) \rangle = 0.056$. By comparing the density profiles of gas and sink particles, we conclude that this decrease is caused by gas accretion and subcluster shrinkage to approximately the same degree.

(iv) Because the subclusters are relatively gas poor, they are only weakly affected by gas expulsion and the subsequent evolution towards virial equilibrium. According to our analytical estimate, they expand by an average factor of 1.08 after gas removal. The length-scale on which the subclusters are gas poor ($\langle f_{gas} \rangle < 0.1$) is about 0.1 pc at the end of the simulation. By scaling up the units of the simulation from subcluster to star cluster scales, we find that after about 0.8 Myr of star formation, star clusters with a mean mass of $10^3$ M$_\odot$ would be gas poor on a length-scale of 0.5–1 pc.

(v) Only those (sub)clusters with densities much larger than the density that is associated with a free-fall time equal to the gas expulsion time can exhibit the shrinkage and accretion that causes them to become gas poor. The fraction of clusters that reaches the required density to become insensitive to gas expulsion before the onset of gas removal therefore increases with ambient gas density. This suggests a relation between the CFE and the ambient gas or star formation rate density, with a larger fraction of star formation resulting in bound clusters in dense regions.

(vi) A possible relation between the CFE and the ambient gas density is affected by a second mechanism. In dense regions, the survival chances of stellar structure are not determined by gas expulsion or 'infant mortality', but by the disruptive effect of tidal shocks from the surrounding gas (Elmegreen & Hunter 2010; Kruijssen et al. 2011). This disruption of young clusters by their environment is titled the 'cruel cradle effect' and is suggested to take over the disruptive role of gas expulsion in dense star-forming regions. The strength and relative contributions of infant mortality and the cruel cradle effect as a function of ambient gas density will be the decisive factor to assess the relation between the CFE and the ambient gas density. This needs to be quantified in future studies, because the time-scale covered by the simulation is too short to include the disruption of subclusters due to the cruel cradle effect.

Throughout the paper, we have mentioned directions in which further research should be undertaken to verify and expand our conclusions. It is essential to check to what extent these results depend on the properties of the simulations we analysed, such as their initial conditions and input physics. The three chief concerns would be whether the results are affected by (1) the turbulence spectrum and initial set-up of the simulation, (2) the inclusion or exclusion of feedback and magnetic fields and (3) the choice of sink particle radius.

(i) The turbulence spectrum and overall boundedness of the simulation primarily influence the time evolution of the overall star formation efficiency (McKee & Ostriker 2007; Dale & Bonnell 2008), and will only impact the evolution of subclusters if the inflow of gas becomes too high to be compensated by accretion and subcluster shrinkage. Judging the relative ease at which the subclusters in the current simulation become gas poor, it is unlikely that the situation of a saturating gas inflow would take place. However, a dynamical analysis of a larger set of simulations will be needed to obtain a conclusive picture, also to include the formation of stars and subclusters in environments with lower densities.

(ii) Feedback from accreting sink particles would inhibit the inflow of the gas, which in turn would lead to a lower gas fraction within the subclusters. There have been several efforts in literature to quantify the effect of (positive or negative) feedback on the star formation process, which generally consider effects on the length-scales of the giant molecular clouds in which the star formation takes place (see e.g. McKee & Ostriker 2007). While global effects could influence the gas-poor state of the (sub)clusters, the nature of feedback is such that it evacuates the gas on the spatial scales of the subclusters, which therefore should not lead to a fundamentally different conclusion than made in this paper. Magnetic fields could slow down the accretion process. If this decrease of the accretion rate would be smaller in the outskirts of the subclusters than within the stellar half-mass radius, it would increase the gas fractions and virial ratios of the stellar component. Therefore, the influence of magnetic fields on our results will need to be investigated.

(iii) If the accretion and/or sink radii of the sink particles were comparable to the typical interstellar separation, the gas-poor nature of the subclusters would be a trivial result of a high 'filling factor' of the subclusters by the sink particles, as the volume where the gas could reside without being accreted would be too small to achieve a stable configuration. We have addressed this to first order by computing the accretion and sink volumes taken up by sink particles and dividing it by the enclosed volume. This was done for each sink particle while taking the nearest neighbour,<fnref n="14" /> and also by calculating a mean radial 'filling factor' profile for each subcluster, analogous to Fig. 8. Both methods returned low filling factors, with typical values of $10^{-2}$ for the sink radius and $10^{-4}$ for the accretion radius. In other words, less than 1 per cent of the volume inside the subclusters is taken up by the sink particles. In order to assess to which extent this allows us to neglect spurious accretion, we ran a set of simple test simulations with different accretion and sink particle radii. These tests show that the gas accretion rate is hardly affected by the accretion and sink radii, which validates our results. The details of the tests are given in Appendix A.

Ideally, the next step would be to perform the same type of calculation as in Bonnell et al. (2008) for different initial conditions,

<footnote n="14">¹⁴ We used the particle list that was corrected for multiples, implying that bound neighbours were ignored.</footnote>
including descriptions for radiative and mechanical feedback, potentially accounting for magnetic fields and varying the sink particle radius. The dynamical analysis of such simulations would provide a good verification of our conclusions, and would improve the current understanding of the dependence on initial conditions and input physics.

The order-of-magnitude extension of our results from subcluster to actual star cluster scales should be investigated further. With the continuously improving computational facilities, it will be possible to simulate systems on the scales needed to cover the formation of star clusters. The key ingredients of such an effort will be larger particle numbers and smaller sink particle radii. Additionally, infrared or spectroscopic observations can be used to verify the length-scales on which star-forming regions are gas poor prior to the onset of gas removal. The current and upcoming generation of telescopes will provide excellent opportunities for this.

If gas expulsion indeed only weakly affects the survival chances of stellar structure, it will need to be verified in which regimes infant mortality still plays a role. In order to understand the relation between the CFE and the local environment, the relative contributions to early star cluster disruption of infant mortality and the cruel cradle effect will need to be quantified. The kinematic signatures of both effects should differ and would therefore be measurable in the velocity maps of young disrupted clusters (Kruijssen 2011). Possible ways in which this could be done observationally include searching for young clusters that are currently going through gas expulsion and mapping the radial velocities of the stars, or tracing the velocity dispersion profiles of young, gas-poor clusters in dense regions. To aid this effort, the differences between the kinematic signatures of energy injection into a star cluster by gas expulsion or tidal shocks have to be established theoretically. The combination of these approaches should provide a conclusive picture of the mechanisms that determine which fraction of star formation results in bound star clusters.

## ACKNOWLEDGMENTS

We thank the anonymous referee for thoughtful comments. JMDK is grateful to Eli Bressert for helpful discussions and acknowledges the kind hospitality of the Institute of Astronomy in Cambridge, where a large part of this work took place. TM acknowledges funding by Constellation, a European Commission FP6 Marie Curie Research Training Network. This research is supported by the Leids Kerkhoven-Bosscha Fonds (LKBF) and the Netherlands Organisation for Scientific Research (NWO), grant 021.001.038.

## REFERENCES

* Adams F. C., 2000, ApJ, 542, 964

* Allison R. J., Goodwin S. P., Parker R. J., de Grijs R., Portegies Zwart S. F., Kouwenhoven M. B. N., 2009a, ApJ, 700, L99

* Allison R. J., Goodwin S. P., Parker R. J., Portegies Zwart S. F., de Grijs R., Kouwenhoven M. B. N., 2009b, MNRAS, 395, 1449

* Bastian N., Goodwin S. P., 2006, MNRAS, 369, L9

* Bastian N., Ercolano B., Gieles M., Rosolowsky E., Scheepmaker R. A., Gutermuth R., Efremov Y., 2007, MNRAS, 379, 1302

* Bastian N., Covey K. R., Meyer M. R., 2010, ARA&A, 48, 339

* Bate M. R., Bonnell I. A., Bromm V., 2003, MNRAS, 339, 577

* Baumgardt H., Kroupa P., 2007, MNRAS, 380, 1589

* Boily C. M., Kroupa P., 2003a, MNRAS, 338, 665

* Boily C. M., Kroupa P., 2003b, MNRAS, 338, 673

* Bonnell I. A., Bate M. R., Zinnecker H., 1998, MNRAS, 298, 93

* Bonnell I. A., Bate M. R., Vine S. G., 2003, MNRAS, 343, 413

* Bonnell I. A., Clark P., Bate M. R., 2008, MNRAS, 389, 1556

* Bonnell I. A., Smith R. J., Clark P. C., Bate M. R., 2011, MNRAS, 410, 2339

* Bressert E. et al., 2010, MNRAS, 409, L54

* Cartwright A., Whitworth A. P., 2004, MNRAS, 348, 589

* Casartano S., Hut P., 1985, ApJ, 298, 80

* Dale J. E., Bonnell I. A., 2008, MNRAS, 391, 2

* Elmegreen B. G., Hunter D. A., 2010, ApJ, 712, 604

* Geyer M. P., Burkert A., 2001, MNRAS, 323, 988

* Gieles M., Portegies Zwart S. F., 2011, MNRAS, 410, L6

* Gieles M., Portegies Zwart S. F., Baumgardt H., Athanassoula E., Lamers H. J. G. L. M., Sipior M., Leenaarts J., 2006, MNRAS, 371, 793

* Goddard Q. E., Bastian N., Kennicutt R. C., 2010, MNRAS, 405, 857

* Goodwin S. P., 2009, Ap&SS, 324, 259

* Goodwin S. P., Bastian N., 2006, MNRAS, 373, 752

* Gutermuth R. A., Megeath S. T., Myers P. C., Allen L. E., Pipher J. L., Fazio G. G., 2009, ApJS, 184, 18

* Heggie D., Hut P., 2003, The Gravitational Million-Body Problem: A Multi-disciplinary Approach to Star Cluster Dynamics. Cambridge Univ. Press, Cambridge

* Hillenbrand L. A., Hartmann L. W., 1998, ApJ, 492, 540

* Hills J. G., 1980, ApJ, 235, 986

* Kennicutt R. C., Jr, 1998, ApJ, 498, 541

* Kirk H., Myers P. C., 2011, ApJ, 727, 64

* Klessen R. S., Burkert A., 2000, ApJS, 128, 287

* Kraus A. L., Hillenbrand L. A., 2008, ApJ, 686, L111

* Kruijssen J. M. D., 2011, preprint (arXiv:1107.2114)

* Kruijssen J. M. D., Pelupessy F. I., Lamers H. J. G. L. M., Portegies Zwart S. F., Icke V., 2011, MNRAS, 414, 1339

* Lada C. J., Lada E. A., 2003, ARA&A, 41, 57

* Lamers H. J. G. L. M., Gieles M., Portegies Zwart S. F., 2005, A&A, 429, 173

* Larsen S. S., 2004, A&A, 416, 537

* Larsen S. S., Richtler T., 2000, A&A, 354, 836

* Larson R. B., 2005, MNRAS, 359, 211

* Launhardt R. et al., 2010, ApJS, 188, 139

* McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565

* McMillan S. L. W., Vesperini E., Portegies Zwart S. F., 2007, ApJ, 655, L45

* Maschberger T., Clarke C. J., 2011, MNRAS, 416, 541

* Maschberger T., Clarke C. J., Bonnell I. A., Kroupa P., 2010, MNRAS, 404, 1061

* Moeckel N., Bate M. R., 2010, MNRAS, 404, 721

* Moeckel N., Bonnell I. A., 2009, MNRAS, 400, 657

* Moeckel N., Clarke C. J., 2011, MNRAS, 410, 2799

* Offner S. S. R., Hansen C. E., Krumholz M. R., 2009, ApJ, 704, L124

* Parker R. J., Goodwin S. P., 2007, MNRAS, 380, 1271

* Parker R. J., Bouvier J., Goodwin S. P., Moraux E., Allison R. J., Guieu S., Güdel M., 2011, MNRAS, 412, 2489

* Parmentier G., Kroupa P., 2011, MNRAS, 411, 1258

* Parmentier G., Goodwin S. P., Kroupa P., Baumgardt H., 2008, ApJ, 678, 347

* Pflamm-Altenburg J., Weidner C., Kroupa P., 2007, ApJ, 671, 1550

* Plummer H. C., 1911, MNRAS, 71, 460

* Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431

* Schmeja S., Kumar M. S. N., Ferreira B., 2008, MNRAS, 389, 1209

* Schmidt M., 1959, ApJ, 129, 243

* Smith R. J., Clark P. C., Bonnell I. A., 2009, MNRAS, 396, 830

* Tutukov A. V., 1978, A&A, 70, 57

* Verschueren W., 1990, A&A, 234, 156

* Zahn C. T., 1971, IEEE Trans. Comput., 20, 68

* Zepf S. E., Ashman K. M., English J., Freeman K. C., Sharples R. M., 1999, AJ, 118, 752

## APPENDIX A: INDEPENDENCE OF RESULTS ON SINK PARAMETERS

In this appendix we verify that the resolution of the SPH simulation is not playing an important role in the evolution of the stellar-to-gas
<figure><img src="image_11.png" /><figcaption>**Figure A1.** Remaining gas mass as a function of time for our tests of the SPH simulation sink parameters. Top panel: the runs with $r_{\text{gas}} = 3r_{\text{sinks}}$; bottom panel: the runs with $r_{\text{gas}} = 10r_{\text{sinks}}$. Grey lines have $r_{\text{acc}} = 0.125$, black lines have $r_{\text{acc}} = 0.0625$ and light grey lines have $r_{\text{acc}} = 0.03125$. Solid lines have $5 \times 10^4$ gas particles. The dashed lines have $2 \times 10^5$ particles, and the dotted lines have $1.25 \times 10^4$ particles.</figcaption></figure>

mass ratio of the subclusters. We accomplish this via a series of controlled, idealized tests in which a cluster of 10 sink particles accretes from an envelope of gas. The total mass of the system is 1, divided equally between the sinks and the gas. The sinks are initially

of equal mass, thus each has mass 0.05. They are placed randomly in a Plummer model of virial radius $r_{\text{sinks}} = 1$, and we use the same initial configuration of the stars in each test. The median nearest neighbour separation of the sinks is 0.43. The gas is likewise in a Plummer sphere spatially, although with a larger radius than the sinks. The gas has zero initial kinetic energy and minimal thermal support, so that the gas falls on to the sink cluster and is accreted. We run two sets of tests, one in which the gas sphere's virial radius is 10 times $r_{\text{sinks}}$, i.e. $r_{\text{gas}} = 10$, and one in which $r_{\text{gas}} = 3$.

The two numerical scales we are concerned with are the accretion radius of the sinks $r_{\text{acc}}$, and the smoothing length of the gas particles. For the sink radii, we use the set $r_{\text{acc}} = \{0.125, 0.0625, 0.03125\}$. The middle value yields approximately the ratio of the neighbour distance to the accretion radius seen in the clusters in the simulation. The smoothing length of the gas is determined by the number of gas particles. To roughly match the simulated value, suppose the sinks have masses $1 M_{\odot}$. The total gas mass is then $10 M_{\odot}$, and $5 \times 10^4$ gas particles approximates the resolution of the large-scale simulation. We run the 0.0625 accretion radius cases with four times more and fewer gas particles, i.e. $2 \times 10^5$ and $1.25 \times 10^4$.

In Fig. A1 we show the gas mass as a function of time for the test runs. In the top panel we show the results for the case with $r_{\text{gas}} = 3$. Time is measured dimensionlessly where we have taken the gravitational constant $G = 1$; the crossing time of the sink cluster is $\sim 2$. In this set-up, the gas free-fall time is $\sim 6$ and the gas accretes quickly compared to the time for the N-body dynamics to dissolve the small-N sink system. The agreement between all the test runs is excellent. The bottom panel shows the $r_{\text{gas}} = 10$ cases, and the gas falls on to the sink system over a longer time-scale, with a free-fall time $\sim 35$. At early times the agreement is quite good, with some disagreement between the runs appearing after $t \sim 20$. We attribute this to the fact that at this point the N-body dynamics of the different runs have set the clusters on clearly divergent paths; recall that the gravitational smoothing length of the sinks is proportional to their sink radius. By $t = 80$ the cluster of sinks has effectively dissolved. We conclude that gas accretion time-scale is not greatly affected by the choices of the sink radius.

This paper has been typeset from a TeX/Latex file prepared by the author.