Impact of a star formation efficiency profile on the evolution of open clusters
^{1} Astronomisches RechenInstitut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstr. 1214, 69120 Heidelberg, Germany
email: beks@ari.uniheidelberg.de
^{2} Fesenkov Astrophysical Institute, Observatory 23, 050020 Almaty, Kazakhstan
^{3} The International Center of Future Science of the Jilin University, 2699 Qianjin St., 130012 Changchun City, PR China
^{4} National Astronomical Observatories of China and Key Laboratory for Computational Astrophysics, Chinese Academy of Sciences, 20A Datun Rd, Chaoyang District, 100012 Beijing, PR China
^{5} Main Astronomical Observatory, National Academy of Sciences of Ukraine, 27 Akademika Zabolotnoho St, 03680 Kyiv, Ukraine
Received: 13 February 2017
Accepted: 5 June 2017
Aims. We study the effect of the instantaneous expulsion of residual starforming gas on star clusters in which the residual gas has a density profile that is shallower than that of the embedded cluster. This configuration is expected if star formation proceeds with a given starformation efficiency per freefall time in a centrally concentrated molecular gas clump.
Methods. We performed direct Nbody simulations whose initial conditions were generated by the program “mkhalo” from the package “falcON”, adapted for our models. Our model clusters initially had a Plummer profile and are in virial equilibrium with the gravitational potential of the clusterforming clump. The residual gas contribution was computed based on a localdensity driven clustered star formation model. Our simulations included mass loss by stellar evolution and the tidal field of a host galaxy.
Results. We find that a star cluster with a minimum global star formation efficiency (SFE) of 15 percent is able to survive instantaneous gas expulsion and to produce a bound cluster. Its violent relaxation lasts no longer than 20 Myr, independently of its global SFE and initial stellar mass. At the end of violent relaxation, the bound fractions of the surviving clusters with the same global SFEs are similar, regardless of their initial stellar mass. Their subsequent lifetime in the gravitational field of the Galaxy depends on their bound stellar masses.
Conclusions. We therefore conclude that the critical SFE needed to produce a bound cluster is 15 percent, which is roughly half the earlier estimates of 33 percent. Thus we have improved the survival likelihood of young clusters after instantaneous gas expulsion. Young clusters can now survive instantaneous gas expulsion with a global SFEs as low as the SFEs observed for embedded clusters in the solar neighborhood (15–30 percent). The reason is that the star cluster density profile is steeper than that of the residual gas. However, in terms of the effective SFE, measured by the virial ratio of the cluster at gas expulsion, our results are in agreement with previous studies.
Key words: galaxies: star clusters: general / stars: kinematics and dynamics / methods: numerical / stars: formation
© ESO, 2017
1. Introduction
The formation of bound star clusters can be divided into three phases: 1) star formation (SF); 2) expulsion of the residual starforming gas; and 3) violent relaxation, that is, the cluster dynamical response to gas expulsion. The dynamics of stars in young clusters during their formation and after gas expulsion is not fully understood yet. It has been the object of intense scrutiny over the past years, with the applied methods ranging from pure Nbody simulations (Tutukov 1978; Lada et al. 1984; Boily & Kroupa 2003b; Baumgardt & Kroupa 2007; Goodwin 2009; Proszkow & Adams 2009; Smith et al. 2011; Lee & Goodwin 2016) and combined hydrodynamical and Nbody simulations (Bonnell et al. 2011; Girichidis et al. 2012; Moeckel et al. 2012; Fujii & Portegies Zwart 2016) to analytical and semianalytical models (Hills 1980; Adams 2000; Boily & Kroupa 2003a; Parmentier & Pfalzner 2013).
Star clusters form in dense (>5–10 × 10^{3} cm^{3}) clumps of gas inside giant molecular clouds (GMC; Lada & Lada 2003; Lada et al. 2010; Kainulainen et al. 2014). The global star formation efficiency (SFE_{gl}), the mass fraction of a starforming region converted into stars, is defined as (1)where M_{⋆} is the total stellar mass and M_{gas} is the mass of unprocessed gas.
The star formation efficiencies (SFEs) measured from observations vary from a few to 30 percent for the dense clumps of molecular clouds (Lada & Lada 2003; Higuchi et al. 2009), and from 0.1 percent to a few percent for their host GMCs (Evans et al. 2009; Murray 2011).
Several mechanisms (stellar winds, ionizing radiation, radiation pressure, Type II supernova explosions) interrupt the SF process and blow up the unprocessed gas out of the cluster (Krumholz & Matzner 2009; Murray et al. 2010; Dib et al. 2013; Hopkins et al. 2013). The observed open clusters older than 10 Myr are already gas free (Leisawitz et al. 1989; Lada & Lada 2003). The duration of SF is on order of 1 Myr, with observations of young star clusters revealing stellar age spreads ranging approximately between 0.3 Myr and 5.0 Myr (Reggiani et al. 2011; Kudryavtseva et al. 2012).
The combination of the SF duration with the SFE per freefall time determines the global SFE achieved by a clusterforming clump at the time of gas expulsion. The SFE per freefall time, the fraction of gas mass converted into stars over one freefall time, was estimated to be 0.01 by Krumholz & Tan (2007), while Murray (2011) suggested that it varies between 0.01 and 0.50 depending on the GMC mass.
Baumgardt & Kroupa (2007) performed a grid of simulations assuming an embedded cluster in virial equilibrium with the residual gas, where both mass density profiles have identical shape. From these Nbody simulations (see also Fig. 1 in Parmentier & Gilmore 2007, for an overview of earlier works) it was concluded that a global SFE of at least 33 percent is needed to form a bound cluster after instantaneous gas expulsion. This minimum global SFE is slightly higher than the SFEs observed for molecular clumps, which vary by up to 30 percent and are frequently estimated to be around 10 percent (Higuchi et al. 2009; Murray 2011; Kainulainen et al. 2014). To address this discrepancy between theoretical works and observations, different solutions exist: adiabatic gas expulsion(Lada et al. 1984; Geyer & Burkert 2001; Baumgardt & Kroupa 2007; Brinkmann et al. 2017), a subvirial cluster at the time of gas expulsion (Verschueren & David 1989; Goodwin 2009; Farias et al. 2015), or hierarchically formed clusters (Smith et al. 2011; Lee & Goodwin 2016).
Goodwin (2009) stressed that the critical factor for a cluster to survive gas expulsion is its dynamical state (as measured by its virial ratio) at the onset of gas expulsion, and not its global SFE. Because star clusters are not necessarily in equilibrium, he introduced an alternative SFE derived from the virial ratio of a star cluster measured immediately after instantaneous gas expulsion and called effective SFE (eSFE; see also Goodwin & Bastian 2006): (2)with the virial ratio defined by (3)Here T_{⋆} and Ω_{⋆} are the total kinetic and potential energies of a star cluster immediately after gas expulsion, and Q_{⋆} = 0.5 corresponds to virial equilibrium. The eSFE is equivalent to the global SFE for these models, where the SFE is constant with the distance to the center of the starforming clump (i.e., stars and gas follow the same density profiles shape), and stars are in virial equilibrium with the gravitational potential of residual gas. According to Goodwin (2009), clusters whose virial ratio is Q_{⋆}< 1.5 (equivalent to eSFE> 0.33) can survive the instantaneous gas expulsion.
Smith et al. (2011, 2013) and Farias et al. (2015) studied this problem by proposing a hierarchical merging scenario of substructured embedded clusters. They also concluded that the dynamical state of a cluster at the time of gas expulsion is important to cluster survival, while the global SFE is not. However, their distributions of stars and starforming gas are different, while not depending on each other.
The more recent work by Lee & Goodwin (2016) proposed different dynamical states for the subclusters, which are in virial equilibrium with each other within the clusterforming region. The authors also concluded that the total dynamical state of the whole cluster at gas expulsion onset is the most important factor in predicting whether the cluster survives gas expulsion. They did not link the formation of a bound cluster to its global SFE, however, and, considered only eSFEs. One could nevertheless map their virial ratio to a subcluster mean SFE, assuming that in each subcluster, stars and gas present the same density profile shape.
Parmentier & Pfalzner (2013), however, proposed a semianalytical model of cluster formation in which the density profile of the embedded cluster is steeper than that of the clusterforming gas. That is, the SFE varies locally, as it steadily increases from the clump outskirts to the clump inner regions (see Fig. 10 in Parmentier & Pfalzner 2013). The reason is that the clusterforming clump is denser and therefore experiences faster SF in its central regions than in its outskirts. The response to the gas expulsion of a cluster like this differs from the most often investigated case, that is, the case where the gas and stars density profiles have identical shapes. This was first investigated by Adams (2000) with a semianalytical method. His choice for different gas and star density profiles was not physically motivated, however. Adams (2000) showed that if the stellar mass is more concentrated in the cluster center than the gas mass, the cluster survival probability is significantly increased (see his Fig. 3). The reason is a gaspoor region in the cluster central regions, which promotes the formation of a bound cluster, even when the global SFE is low. In addition, Pfalzner et al. (2014) performed Nbody simulations whose initial conditions build on the model of Parmentier & Pfalzner (2013) for a global SFE of about 18 percent.
In this contribution, we revisit this problem by means of direct Nbody simulations. Specifically, we build on the model of Parmentier & Pfalzner (2013) to generate the initial conditions of our Nbody simulations, thereby physically motivating our choice for the stars and gas density profiles of the embedded cluster at the time of the gas expulsion. Our simulations also take into account stellar evolution and the tidal field of a Milky Waylike galaxy.
The outline of the paper is as follows. In Sect. 2 we build on the model of Parmentier & Pfalzner (2013) to relate the density profiles of the embedding gas and of the cluster at the onset of gas expulsion. Section 3 presents the technical setup and the initial conditions of our Nbody simulations.
In Sect. 4 we discuss the results of our simulations and compare them to earlier works. Finally, we present our conclusions in Sect. 5.
2. Starformation efficiency and density profiles
A semianalytical model of star cluster formation from centrally concentrated spherically symmetric gas clumps with a constant SFE per freefall time, ϵ_{ff}, was developed by Parmentier & Pfalzner (2013). Because the freefall time is shorter in the clump inner (denser) regions than in the clump outer (less dense) regions, the density profile of the formed star cluster is steeper than the density profiles of the initial and residual gas. The authors considered that the total density profile ρ_{0} of the system remains constant.
The total density profile is the sum of the density profiles of the embedded cluster, ρ_{⋆}, and the residual gas, ρ_{gas}, at any time t after SF onset: (4)here r is the distance to the clump center. The density profile of the unprocessed gas at time t is described with Eq. (19) from Parmentier & Pfalzner (2013), which we reproduce here for the sake of clarity: (5)G is the gravitational constant and ϵ_{ff} is the SFE per freefall time. The mass of the embedded cluster at time t is distributed according to Eq. (20) in Parmentier & Pfalzner (2013): (6)Our aim is to investigate the dynamical evolution of such star clusters after instantaneous gas expulsion and to estimate their survival likelihood. The instantaneous gas expulsion corresponds to the case when the timescale of the gas expulsion is significantly shorter than the dynamical timescale of the system. In terms of cluster survivability, it is the worst scenario we can envision. If a star cluster survives instantaneous gas expulsion, it is also able to survive a longer gas expulsion timescale.
We can adopt two different approaches to study this problem.

A.
Either the starting point is a molecular clump, with a givendensity profile ρ_{0}(r), and we obtain the density profile of the star cluster that formed after some SF time span based on Parmentier & Pfalzner (2013).

B.
Alternatively, the starting point is an embedded cluster with a wellknown profile (e.g., Plummer or King), and assuming an SF time span, we recover the initial gas density profile of the molecular clump out of which the cluster has formed.
In case (A) we start with the initial spherically symmetric gas clump that has a certain density profile ρ_{0}(r) at time t = 0. Then we assume that a star cluster forms within a time interval t = t_{SF} called SF duration with a constant SFE per freefall time ϵ_{ff} = const. Its density profile is then given by Eq. (6). Depending on how long the SF process lasts, star clusters with different SFE_{gl} are formed (see Fig. 9 in Parmentier & Pfalzner 2013). At time t = t_{SF} (corresponding to a certain value of SFE_{gl}), we set the instantaneous gas expulsion, that is, we remove the unprocessed gas from the system. A star cluster then becomes supervirial because it has lost part of the gravitational potential within which it was in virial equilibrium. This would define the initial conditions of our direct Nbody simulations to study the effect of instantaneous gas expulsion.
Because model (A) depends on the SF duration, it results in star clusters with different global SFEs, masses, and spatial and velocity distributions, which makes it difficult to compare the results of the Nbody simulations to each other. Additionally, generating the initial conditions of such models for Nbody simulations is not trivial.
Our main point is to study the effect of instantaneous gas expulsion on a star cluster in dependence of its global SFE, therefore we can simplify the problem and consider clusters with a fixed stellar mass and spatial distribution while varying the global SFE at gas expulsion. This leads us to our case (B), on which we focus in this paper.
In case (B), we thus assume a fixed density profile and a fixed stellar mass for the embedded clusters at gas expulsion. Then the initial clumps that formed such clusters have therefore different total masses and spatial distributions depending on the assumed global SFE.
As we wish to study the response of a star cluster to gas expulsion as a function of the global SFE, we have to find the cluster velocity distribution for any given global SFE assuming it is in virial equilibrium with the gravitational potential of the residual gas. To find this distribution, which depends on the SF duration, we need to solve the inverse problem to that presented in case (A). That is, having the density distribution of a star cluster ρ_{⋆}(r), we determine the density profile of the residual gas at gas expulsion ρ_{gas}(r,t_{SF}), and following from this, the density profile of the whole clusterforming clump, ρ_{0}(r,t_{SF}). Then we modify the cluster velocity distribution function so as to account for its virial equilibrium with the gravitational potential of the residual gas.
The model developed by Parmentier & Pfalzner (2013) can be applied to any clump density profile. So we can choose a wellknown density profile ρ_{⋆} for the embedded cluster and define the density profile of the residual gas corresponding to an SF duration t_{SF} by modifying Eqs. (4) and (5) and setting t = t_{SF}: In this equation, we note that the SF duration t_{SF} is given in units of Myr if densities ρ_{(⋆, gas)} are expressed in M_{⊙} pc^{3} and the gravitational constant . Introducing the parameter (9)which depends on the SFE per freefall time ϵ_{ff} and SF duration t_{SF}, we can rewrite Eq. (8) as (10)Solving this equation provides us with the residual gas density profile ρ_{gas} as a function of the stellar density profile ρ_{⋆}, SFE per freefall time, ϵ_{ff}, and SF duration t_{SF}. The stellar density profile, ρ_{⋆}, can be any centrally concentrated spherically symmetric density profile. We find that two of the four roots of Eq. (10) are complex. The other two roots are real, one increasing and one decreasing function of stellar density. Since we consider a centrally concentrated clump, where the stellar density decreases with increasing radius, the gas density should follow the stellar density and decrease as well. Thus we choose the root that is an increasing function of stellar density. The mathematical expressions and a short analysis of these roots are provided in Appendix A.
Parmentier & Pfalzner (2013) used a powerlaw density profile with a slope of −2 for their clusterforming clumps. This yields a powerlaw density profile with a slope of about −3 for a star cluster. Such initial conditions are not ideal for Nbody simulations because of their infinite stellar and gas masses. We need either to truncate these powerlaw profiles, or choose steeper profiles with finite masses. Thus we decide to use one of the wellknown spatial density distribution functions for an embedded cluster, that is, the Plummer profile (Plummer 1911), (11)where M_{⋆} is the cluster total mass and a_{⋆} is the Plummer radius, which corresponds to the projected halfmass radius of a star cluster.
Choosing a Plummer profile has many advantages, for instance, a finite mass for both gas and stars, and an analytical expression. It is also supported by almost all Nbody codes, which makes it possible to compare the results of different works.
3. Technical setup
3.1. Initial conditions
Since our aim is to perform Nbody simulations, we adopted dimensionless Nbody ([NB]) units. They are associated with the star cluster parameters, not with the clusterforming clump parameters: (12)The Nbody units can be converted into physical units when G,a_{⋆},M_{⋆} are assigned numerical values in their respective units. Here we note that our “Nbody” time unit depends on the cluster stellar mass M_{⋆} and Plummer radius a_{⋆} at the time of the instantaneous gas expulsion. It does not represent the dynamical timescale of the cluster (stars only) because the cluster is supervirial and its dynamics bears the imprint of the clusterforming clump mass at the time of the gas expulsion. Neither does it represent the dynamical timescale of the clump (stars+gas) since the total mass and halfmass radius of the clusterforming clump differ from the stellar mass and stellar halfmass radius because their density profiles have different shapes. We applied these Nbody units for (i) to recover the density profile of the residual clusterforming gas with Eq. (10); and (ii) to perform the subsequent Nbody integration of the gasfree cluster after instantaneous gas expulsion.
To generate the initial conditions of our Nbody simulations is not trivial because we need a star cluster in virial equilibrium with the gravitational potential of the residual gas, where the shapes of the gas and star distributions differ. In that respect, our case differs from most previous Nbody simulations, as radial variations of the SFE increase the degree of complexity of the problem.
We used the falcON program mkhalo by McMillan & Dehnen (2007), which produces a spherically symmetric star cluster in virial equilibrium with an external potential as the initial conditions of direct Nbody simulations. External potential means a gravitational potential produced by anything but the stars of the cluster. In this framework, the gravitational potential produced by the residual gas constitutes an external potential.
To use mkhalo for our purpose, we wrote an additional acceleration plugin, that is, an additional code, which takes into account the new external potential of our models. In this plugin we calculated the gravitational potential and the forces produced by the residual gas knowing its density profile. For this, we used Eq. (3.15′) from Duboshin (1968): where ρ_{gas} was obtained by solving Eq. (10) (see Sect. 2), and R_{gas} is the adopted outer edge of the clump of residual gas. We used R_{gas} = 32a_{⋆}, which is the smallest radius possible to use in mkhalo. Because the gas density profile, ρ_{gas}, is not a simple function of r, the distance to the clump center, Eqs. (13) and (14) were integrated numerically using the Simpson method.
In the framework of this study, we adopted a fixed SFE per freefall time of ϵ_{ff} = 0.05. Then the only parameter we varied in our acceleration plugin to produce the initial conditions is the SF duration t_{SF}. To allow a comparison of our results with other works, however, it is better to use as a main parameter the SFE_{gl} than t_{SF}. To develop a grid of models with a given SFE_{gl} (0.05, 0.10, ...), we still need to infer the corresponding t_{SF} and add them to the models.
We defined the global SFE as the ratio between the stellar and total (stellar + gas) masses residing inside a chosen outer limit, R_{cl}. Because a Plummer model has no finite outer limit, we adopted R_{cl} = 10a_{⋆}, which is the radius inside which about 98 per cent of the stellar mass resides. Here we note that because of the slope difference between the density profiles of the embedded cluster and the residual (as well as total) gas, a larger outer limit would imply a lower global SFE, and vice versa.
The relation between the SFE_{gl} and the corresponding SF duration t_{SF} calculated in Nbody units is presented in Fig. 1. In our study we concentrated on models with an SFE_{gl}< 0.50. The corresponding values of SFE_{gl} and t_{SF} for different models are presented in Table 1. Using these values, we produced the initial conditions of our simulations with mkhalo. Then, using the generated positions and velocities of the stars, we calculated the initial potential (Ω_{⋆}) and kinetic (T_{⋆}) energies of our model clusters at the moment of instantaneous gas expulsion, as well as their initial virial ratios and effective SFEs, using Eqs. (3) and (2). The corresponding values are also presented in Table 1.
Fig. 1 Relation between the SF duration (t_{SF}) (given in Nbody units according to Eq. (12)) and the global SFE (SFE_{gl}), measured inside R_{cl} = 10a_{⋆}. The lowest value of the global SFE here in this plot is SFE_{gl} = 0.007, which corresponds to an SF duration t_{SF} = 0.01 [NB]. 

Open with DEXTER 
We find that in our models the effective and global SFEs are different, unlike the models that used the same density profile for both the stars and gas (see Boily & Kroupa 2003a; Goodwin & Bastian 2006; Goodwin 2009). Varying the outer edge R_{cl} of the cluster, we infer that these eSFEs are almost the same as the global SFEs measured inside 1.5a_{⋆}, that is, similar within 2–3 percent to the global SFEs measured inside a cluster halfmass radius. The latter also measures the local stellar fraction (LSF) as defined by Smith et al. (2011). Goodwin (2009) noted that star clusters with an initial virial ratio lower than 1.5 are able to survive the instantaneous gas expulsion (in the case of Plummer profiles for both stars and gas), which corresponds to an effective SFE of 33 percent. Taking this into account, we can expect the minimum SFE needed to form a bound cluster to be SFE_{gl} = 0.15, as it corresponds to Q_{⋆} ≈ 1.5 for our models (see Table 1). This is indeed what we show in Sect. 4.
The density profiles of the embedded cluster, its residual gas before instantaneous gas expulsion, and the initial clump gas for different global SFEs are shown in Fig. 2. The models were scaled to physical units assuming a star cluster mass of M_{⋆} = 3000 M_{⊙} and a 3D halfmass radius of r_{h} = 1 pc. The density units are given in M_{⊙} pc^{3} (right yaxis) and molecules per cm^{3} (left yaxis). Gas densities on this scale vary within the observed range of dense clumps (≲10^{5} cm^{3}). The SF durations t_{SF} for these three models are 0.39 Myr, 2.21 Myr, and 12.77 Myr. As we see, some of our models are inconsistent with observations in the sense that the inferred SF duration is longer than what is observed.
Fig. 2 Density profiles of the star cluster (black dashed line), of the residual (solid lines), and initial (dashdotted lines) gas for different SFE_{gl} in scaled physical units. A total stellar mass M_{⋆} = 3000 M_{⊙} and a 3D halfmass radius r_{h} = 1 pc are assumed. Note that the stellar density profile is a Plummer profile. 

Open with DEXTER 
For our highresolution direct Nbody simulations we have chosen the φgrapegpu code developed by Berczik et al. (2011). As a check of the stability of the initial conditions generated by the code mkhalo, we added our newly created external potential to the φgrapegpu code and tested the dynamics of the embedded cluster within the residual gas potential. We ran a few simulations with external potentials corresponding to different t_{SF}. The test runs were performed for isolated clusters with N = 10k singlemass particles over a time interval of up to 1000 Nbody time units, which is slightly shorter than three relaxation times of the same cluster in virial equilibrium without external potential. We checked the evolution of the Lagrangian radii and cumulative mass profiles for two models with low and high values of t_{SF}. These exercises showed that our newly generated initial conditions are indeed in virial equilibrium with the external gas potential.
3.2. Setting up the simulations
We performed two types of Nbody simulations. First, we simulated isolated singlemass clusters without stellar evolution to observe the pure dynamical effect of an instantaneous gas expulsion. We ran these simulations with N = 10^{4} particles and covered global SFEs ranging from 5 to 50 percent.
Second, we studied a more realistic scenario of violent relaxation by considering star clusters consisting of multimass stars that evolve in a Galactic tidal field. We refer to these two sets of simulations as “isolated” and “nonisolated” models.
For the stellar initial mass function (IMF) of our nonisolated models we adopted the IMF of Kroupa (2001) with the lower and upper mass limits of M_{low} = 0.08 M_{⊙} and M_{up} = 100 M_{⊙}, respectively.
To account for the tidal field of the Galaxy, we used an axisymmetric threecomponent model that consists of a bulge, disk, and halo, as described by the PlummerKuzmin model (Miyamoto & Nagai 1975). This was added to the φgrapegpu code by Just et al. (2009), and we used it keeping their parameters.
We considered that our clusters move on a circular orbit in the plane of the Galactic disk, at a distance of R_{0} = 8 kpc from the Galactic center.
We normalized our Nbody units to the real physical units in order to assign the correct timescale to the stellar evolution routines (SSE; Hurley et al. 2000) implemented in the φgrapegpu code. This means that we assigned certain values to the cluster mass (M_{⋆}) and the initial Plummer radius (a_{⋆}).
Our simulations encompass five different initial cluster stellar masses: M_{⋆} = 3000, 6000, 10 000, 15 000, and 30 000 M_{⊙}. Then knowing the distance of the cluster to the Galactic center, we can calculate the cluster tidal (Jacobi) radius r_{t} for a given mass M_{⋆} using Eq. (13) from Just et al. (2009), which we reproduce here for the sake of clarity: (15)Here β = 1.37 is the normalized epicyclic frequency and Ω = V_{0}/R_{0} is the angular speed of a star cluster on a circular orbit at a distance R_{0} from the Galactic center. For R_{0} = 8000 pc, the orbital speed recovered from the rotation curve of the Galaxy model provided in Just et al. (2009) is V_{0} = 234.24 km s^{1}.
To make the models comparable with each other, we fixed the halfmass radius to the tidal radius ratio. We calculated it for a cluster mass M_{⋆} = 3000 M_{⊙} with a halfmass radius of r_{h} = 1 pc. This means that we considered clusters with different stellar masses, but the same mean stellar volume densities. The tidal radius of a cluster with M_{⋆} = 3000 M_{⊙} and r_{h} = 1 pc is r_{t} = 19.2 pc, and therefore (16)Thus we normalized the Nbody length unit into physical units of pc knowing that in a Plummer model r_{h} ≈ 1.3a_{⋆}, (17)where the tidal radius is calculated using Eq. (15). For the models with an initial stellar mass M_{⋆} = 3000 M_{⊙}, for instance, the normalized length is equal to r_{norm} = 0.77 pc. With our definition of the outer limit R_{cl} of our clusterforming clumps (see Sect. 3.1), star clusters initially fill their tidal radius up to 40 percent. This means that the total radius of a star cluster is initially smaller than its tidal radius: R_{cl} = 10a_{⋆} = 0.4r_{t}. The properties of a star cluster as a function of its stellar mass are presented in Table 2.
Set of models corresponding to different initial stellar masses.
The scale factor of time units, as shown in Eq. (12), can be found as (18)when r_{h}/r_{t} = 0.052. Given our assumption of a fixed ratio of the halfmass to tidal radius (see Eqs. (16) and (17)), r_{norm} ∝ r_{t} ∝ (GM_{⋆})^{1/3}, t_{norm} is the same for all models as well as the mean stellar volume densities.
The corresponding values of SF duration and total (stars + gas) volume densities averaged within the initial stellar 3D halfmass radius are provided in Table 3. This table shows that the models are consistent with the observations in terms of SF duration and clump mean densities. It is thought that star cluster formation takes between one half and roughly 5 Myr in the solar neighborhood. We therefore consider the models in these limits to make our simulations as consistent with reality as possible. This limits the SF duration of our models between 2.7 and 27 Nbody time units, which corresponds to 0.5 and 5 Myr for our chosen scale factor of r_{h}/r_{t}. Consequently, this also limits us in the range of achievable SFE_{gl}. According to these limits, we decided to calculate models with SFE_{gl} between 10 percent and 25 percent for ϵ_{ff} = 0.05, which corresponds to an SF duration t_{SF} between 1.15 and 5.25 Myr for all initial cluster masses M_{⋆}. To cover still higher SFEs, we ran two additional models with global SFEs of 30 and 35 percent for M_{⋆} = 6000 M_{⊙}. To make these runs consistent with the observations in terms of SF duration, we built on the following feature of Eq. (10): because the parameter k in Eq. (10) is proportional to ϵ_{ff}t_{SF}, our results stand for any model where the product ϵ_{ff}t_{SF} is conserved (e.g., a twice higher SFE per freefall time with a twice shorter SF duration). In our two additional runs the SF durations can therefore be considered as t_{SF} = 3.64 Myr and 4.88 Myr, respectively, if ϵ_{ff} = 0.1. For comparison, additional models with the same spatial distribution of stars initially in virial equilibrium within a Galactic tidal field, but without any residual starforming gas, (i.e., equivalent to SFE_{gl} = 1.0), were run for M_{⋆} = 3000, 6000, and 10 000 M_{⊙}.
The mean (total and stellar) volume densities of models are also consistent with observations of starforming molecular clumps, where an SF density threshold has been suggested (≳10^{4} cm^{3} in Lada et al. 2010 and >5 × 10^{3} cm^{3} in Kainulainen et al. 2014) and with stellar densities in embedded clusters, which vary from 100–200 to 1–2 × 10^{4} stars pc^{3} (Lada et al. 1991; Hillenbrand & Carpenter 2000).
Based on these simulations, we now study the evolution of the bound fraction of star clusters. We present the results of isolated models scaled to physical units such that M_{⋆} = 6000 M_{⊙} and r_{h} = 1.26 pc (see Table 2) to compare them with M_{⋆} = 6000 M_{⊙} nonisolated models, which consist of 10 455 stars.
4. Results and discussions
4.1. Clusterbound fraction evolution
In the classical way, the bound fraction of isolated clusters in virial equilibrium is defined as the fraction of stars whose total energy is negative, that is, the potential energy dominates the kinetic energy. Our model isolated clusters become supervirial and start to expand quickly after instantaneous gas expulsion, however. Thus it is not trivial to distinguish between bound and unbound stars in such systems during their expansion. This process can take quite a long time in isolated systems, especially if the final bound fraction is small. For instance, the cluster with SFE_{gl} = 0.15 needs 1 Gyr to approach the final bound fraction if this is defined in the “classical” way. We therefore developed the following technique to determine the final bound fraction. We recalculated the total energies of stars, removing the unbound stars (i.e., stars whose total energy is positive) from the cluster even if they were located in the center of the cluster. We iterated until no unbound star remained in the cluster. Doing so, we excluded the contribution of unbound stars to the cluster potential. The final bound fraction of our isolated models is defined as the final fraction of stars that remained bound to the cluster – if any – at the end of the iterations. We emphasize that removing the unbound stars does not mean that we removed them from the simulations, but that we removed them from our selected sample at each snapshot for the purpose of our analysis only. For each snapshot in time, we always started the iterative process with the total number of stars. More discussion about this technique and its relevance can be found in Appendix B.
For the nonisolated clusters, that is, those evolving within a Galactic tidal field, the bound fraction is defined as the stellar mass residing inside the instantaneous tidal radius normalized to the initial stellar mass, (19)In the top panel of Fig. 3 we present the time evolution of the bound fractions of our nonisolated (solid lines) and isolated (dashed lines) models.
Fig. 3 Evolution of the bound fraction of clusters obtained in our simulations. Different colors correspond to different global SFEs (see the key). The gray lines show the impact of mass loss caused by stellar evolution alone. The vertical dotted line corresponds to t = 20 Myr.Top panel: comparison of the bound fraction evolution of isolated (dashed lines) and nonisolated (solid lines) models for M_{⋆} = 6000 M_{⊙} with N ≈ 10^{4} stars. Note that the isolated models are the singlemass models without stellar evolution and scaled to the same physical units as the nonisolated models (i.e., M_{⋆} = 6000 M_{⊙} and r_{h} = 1.26 pc). Bottom panel: bound fraction (instantaneous tidal mass as a fraction of initial stellar mass) evolution for nonisolated models with different initial stellar masses (see the key for the linecoding). 

Open with DEXTER 
We note that for the isolated clusters we present the final bound fraction defined with the technique described above, and not the fraction of stars with negative total energy. This clearly provides a method for determining the bound mass very early (see Appendix B ). For our nonisolated models, we also show the imprint of stellarevolution mass loss as the gray lines. The model clusters with SFE_{gl} = 1.0, that is, those that are initially in virial equilibrium without any residual starforming gas, are depicted as skyblue lines. To better visualize both the very fast evolution shortly after gas expulsion and the cluster longterm evolution, the scale of the xaxis (time) is logarithmic. The plateau at the very beginning of the bound mass evolution results from our definition of the bound mass fraction, which does not account for the highvelocity unbound stars within the tidal radius. Because the cluster initial size is smaller than its tidal radius (r_{h}/r_{t} = 0.052), almost all stars reside within the tidal radius during the first few Myr of evolution even when the cluster starts to expand. The bound mass fraction starts to decrease as the escaping stars reach the tidal radius and become unbound by our definition.
From the bound fraction evolution of our nonisolated model clusters we can identify two regimes of mass loss (the solid lines in the top panel of Fig. 3), in addition to the mass loss driven by stellar evolution. During the first 20 Myr after gas expulsion, clusters intensively lose their mass (top panel of Fig. 3: the solid lines on the lefthand side of the vertical dotted line). During this time span, the cluster evolution is dominated by the consequences of gas expulsion, and their response is mostly determined by the cluster initial virial ratio.More or less flat plateaus can be seen in between two identified massloss regimes around 20 Myr after gas expulsion. This means that the surviving part of the cluster is not expanding anymore.
The bound fraction then decreases more slowly with time. It is now mostly affected by stellar evolution and the tidal field of the host galaxy (top panel of Fig. 3: the solid lines on the righthand side of the vertical dotted line).
The bottom panel of Fig. 3 presents results for the nonisolated models alone for the global SFEs and initial cluster stellar masses quoted in the key. The colorcoding is identical to the coding used in the top panel. We find that the models with identical global SFE show similar evolutionary tracks within the first 20 Myr (i.e., during the violent relaxation) independently of their initial stellar mass (bottom panel of Fig. 3). In particular, the models with a low global SFE dissolve on similar timescales independently of the initial star cluster mass, M_{⋆} (see the black curves). Here we recall that the ratio r_{h}/r_{t} is kept constant for now. The model with a global SFE of 13 percent and M_{⋆} = 6000 M_{⊙} does survive as a bound cluster, although with a very small bound fraction, around 2 percent (the blue lines in the top panel of Fig. 3). All other stellar mass models with a global SFE of 13 percent dissolve, however, except for the M_{⋆} = 30 000 M_{⊙} model, which barely survives with a 0.17 percent bound fraction, which corresponds to a bound mass of about 52 M_{⊙}. Therefore we adopt the models with SFE_{gl} = 0.13 as the limit between survival and dissolution following cluster gas expulsion for our adopted tidal field impact r_{h}/r_{t} = 0.052. Cluster models with a global SFE of 0.15 and higher can survive instantaneous gas expulsion, as expected from their initial virial ratios (see Table 1). We note that SFE_{gl} = 0.13 is about 2.5 times lower than the SFE threshold for cluster survival when (i) the density profiles of the stars and residual gas have the same shape; (ii) gas expulsion is instantaneous; and (iii) the tidal field impact is negligible.
We performed a few simulations for a given parameter set but different random seeds to explore the bound fraction variations that are due to random realizations. Figure 4 shows that the bound mass fraction at t = 20 Myr displays a range of variations of about 10 percent for a cluster with N = 10 455 stars (i.e., M_{⋆} = 6000 M_{⊙} cluster). For a cluster with a higher number of stars (N = 26 138, M_{⋆} = 15 000 M_{⊙}), the range of bound mass fraction variations is about 6 percent. The duration of cluster massloss in response to gas expulsion remains shorter than 20 Myr for different random realizations.
Fig. 4 Different random realizations (see line types in the key) of the model clusters with M_{⋆} = 6000 M_{⊙} (N = 10 455 stars) for r_{h}/r_{t} = 0.052. 

Open with DEXTER 
After violent relaxation, the cluster lifeexpectancy depends on its stellar mass, as expected (red lines in the bottom panel of Fig. 3). A higher stellar mass implies a higher probability to survive a longer time (but see Ernst et al. 2015). The models corresponding to the high global SFEs in this set of simulations show longterm evolutionary patterns similar to the pure Plummer models (i.e., initially in virial equilibrium without any external potential).
This study shows us that the mass loss of the cluster in response to gas expulsion is completed within 20 Myr, independently of its initial stellar mass and global SFE, and that its dynamical evolution is mostly affected by the tidal field of the host galaxy thereafter (Fig. 3). Since we focus on the cluster bound mass evolution, we consider the violent relaxation as the time span when the cluster loses its mass intensively in response to an instantaneous gas expulsion. We therefore used t = 20 Myr to measure the final bound fraction of clusters. We note, however, that the outer shells of surviving clusters need a longer timespan to return to virial equilibrium, as shown by Fig. 5 (see below; see also Brinkmann et al. 2017). We note therefore that the violent relaxation duration, as defined here, does not strictly equate with the cluster revirialization time.
An example of the Lagrange radius R_{f} evolution of a nonisolated cluster with SFE_{gl} = 0.20 and M_{⋆} = 6000 M_{⊙} is presented in Fig. 5.
Fig. 5 Evolution of a nonisolated cluster of the solar neighborhood with SFE_{gl} = 0.20, M_{⋆} = 6000 M_{⊙} and r_{h}/r_{t} = 0.052 over 50 Myr after gas expulsion. We show 49 Lagrange radii, ranging from 2 percent to 98 percent in intervals of 2 percent (solid and dashed lines). The instantaneous tidal radius r_{t} of the cluster is shown as the dotted line. Thick solid lines correspond to every 10 percent of the Lagrange radius, and the dashed line corresponds to 50 percent of the Lagrange radius. 

Open with DEXTER 
The Lagrange radii are defined based on the fraction of initial stellar mass, not on the number fraction of stars. The dotted line in this figure represents the instantaneous tidal radius. This figure shows that the inner parts of a cluster recede, form a bound cluster, and return to virial equilibrium within 20–30 Myr after gas expulsion. The inner shells of a cluster revirialize faster than the outer shells, as also found in Brinkmann et al. (2017), for example. The cluster tidal radius stays roughly constant after 20 Myr.
4.2. Influence of the cluster initial stellar density
To quantify the effect of the tidal field impact, three additional sets of simulations were performed for cluster masses M_{⋆} = 6000 M_{⊙} and 15 000 M_{⊙}: one with a weaker tidal field impact (r_{h}/r_{t} = 0.025), and two with stronger tidal field impact corresponding to r_{h}/r_{t} = 0.07 and 0.10. That is, we varied the halfmass radius of our model clusters, keeping them in the solar neighborhood (i.e., keeping the same tidal radius for a given stellar mass). Because we varied the initial density of the cluster, we varied the normalization factor of time units as shown in Table 4. For denser clusters (smaller r_{h}/r_{t}), a given physical timespan represents a higher number of Nbody time units.
Scaling factors to the physical time and duration of 20 Myr in Nbody units corresponding to different ratios of the cluster halfmass to tidal radius.
We present the bound mass fraction evolution of M_{⋆} = 15 000 M_{⊙} clusters with different initial densities in Fig. 6.
Fig. 6 Evolution of the bound fraction of M_{⋆} = 15 000 M_{⊙} clusters with different impact of the tidal field (see the key for the linecoding). Different colors correspond to different global SFEs (see the key for the colorcoding). 

Open with DEXTER 
We find that the dense clusters evolve quicker than the less dense clusters within the first 20 Myr after gas expulsion. The violent relaxation duration of clusters and their bound mass fraction at t = 20 Myr depend on the initial cluster mean densities. This is consistent with Parmentier & Baumgardt (2012) and Banerjee & Kroupa (2013), who showed that dense clusters have shorter revirialization times than less dense clusters. However, we find the cluster violent relaxation to depend fairly weakly on the initial stellar density. The difference during the first 10 Myr results from our definition of the bound mass fraction. That is, the most compact cluster (r_{h}/r_{t} = 0.025) loses mass in response to gas expulsion faster than the most diffuse cluster (r_{h}/r_{t} = 0.100). Its escaping stars reach the tidal radius twice faster because their velocity dispersions differ by a factor of 2. Although it is not obvious how to define the duration of violent relaxation accurately, Fig. 6 shows that it remains shorter than 20–25 Myr regardless of the cluster initial density. The key to understanding the violent relaxation duration may reside in the mean stellar density within the tidal radius, which is the same for all considered clusters since they all have the same orbit.
We note that stars that would be bound to the cluster if the cluster is isolated now become unbound once they cross the tidal radius of the cluster. Additionally, they are taken away from the cluster by the Galactic tidal field (see the different behaviors of the 30 percent Lagrange radius in Figs. 5 and 7).
In Fig. 7 we present an example of Lagrange radii evolution of a most compact cluster with M_{⋆} = 6000 M_{⊙},SFE_{gl} = 0.20. We can see that the inner shells of the cluster with r_{h}/r_{t} = 0.025 revirialize faster than those of a more diffuse cluster with r_{h}/r_{t} = 0.052.
Figure 6 shows the differences in bound mass fraction at t = 20 Myr to be around 10 percent between the most compact and the most diffuse clusters. More diffuse clusters have a lower bound fraction since their outer shells expand beyond their instantaneous tidal radius. We have checked that a higher number of particles (M_{⋆} = 15 000 M_{⊙}) does not affect our results.
Fig. 7 Evolution of the cluster with M_{⋆} = 6000 M_{⊙}, SFE_{gl} = 0.20 and r_{h}/r_{t} = 0.025 over the first 50 Myr after the instantaneous gas expulsion. 

Open with DEXTER 
4.3. Final bound mass fraction in dependence of SFE
We compare our results (represented by the mean values taken from the model clusters with r_{h}/r_{t} = 0.052 and different stellar masses) with previous works (Adams 2000; Lada et al. 1984; Geyer & Burkert 2001; Boily & Kroupa 2003b; Fellhauer & Kroupa 2005; Baumgardt & Kroupa 2007) in the top panel of Fig. 8, which shows the bound fraction as a function of global SFE.
Fig. 8 Bound fraction as a function of global or effective SFE. We compare our results (red lines) with previous works. The isolated models are depicted by the red diamonds, and nonisolated models by red crosses. In the top panel we use our global SFE (SFE_{gl}) and in the bottom panel, the effective SFE, eSFE = 1 /Q_{⋆} (which is almost the same as the LSF, the local stellar fraction of our models). 

Open with DEXTER 
We note here the improved survival likelihood of star clusters after instantaneous gas expulsion, and that although the tidal field is included in our models. Our model star clusters whose global SFE is lower than 30 percent, but higher than 15 percent, survive and retain a significant fraction of their stars after violent relaxation.
The global SFEs required by our model clusters to survive gas expulsion provide a good match to the SFEs observed in embedded clusters, which are lower than 30 percent (Higuchi et al. 2009; Murray 2011; Kainulainen et al. 2014). Star clusters are now able to survive instantaneous gas expulsion despite a global SFE as low as 15−20 percent. We now have four avenues to understand the presence of star clusters with ages older than few Myr in the solar neighborhood despite the low SFE observed: adiabatic gas expulsion (Lada et al. 1984; Geyer & Burkert 2001; Baumgardt & Kroupa 2007; Brinkmann et al. 2017), subvirial clusters (Verschueren & David 1989; Goodwin 2009; Farias et al. 2015), hierarchically formed clusters (Smith et al. 2011; Lee & Goodwin 2016), and centrally concentrated cluster formation (this contribution).
We argue that the improved survivability of our model clusters is mostly caused by the difference in density profiles between the embedded cluster and its residual gas, as postulated by Parmentier & Pfalzner (2013), namely, the stars have a density profile steeper than that of the residual and initial gas. This is the consequence of SF taking place with a constant SFE per freefall time in a centrally concentrated molecular clump.
We reproduce similar results to previous works when we plot the bound fraction as a function of the effective SFE (eSFE) instead of the global SFE (see the bottom panel of Fig. 8). We note here again that in our models the eSFE and the local stellar fraction (LSF) are almost the same. Our study thus agrees with all present works, including the most recent paper of Lee & Goodwin (2016), who concluded that the effective SFE is the most important parameter in predicting the cluster survivability. Our study allows us to compare model SFEs to those achieved in observed forming clusters, however.
5. Conclusions
We have performed Nbody simulations of violent relaxation and bound cluster formation after instantaneous gas expulsion. The key point of our study is that we used special initial conditions built on the model of Parmentier & Pfalzner (2013). This means that the density profile of our model star cluster is steeper than the density profile of the starforming gas at the time of instantaneous gas expulsion. If this cluster is in virial equilibrium, including the gravitational potential of residual gas, it should be able to survive gas expulsion despite low SFEs, as shown by Adams (2000) based on a semianalytical model.
Since our Nbody simulations start from the time of instantaneous gas expulsion and do not cover the SF phase, we started with a wellknown star cluster model, namely the Plummer model, instead of starting from the starfree molecular clump. For this we obtained the dependency of the residual gas density profile on stellar density profile at the time of instantaneous gas expulsion under the assumptions of Parmentier & Pfalzner (2013). In their clusterformation model, the density profile of a clump (i.e., stars + gas) is constant during SF, and SF takes place with a constant SFE per freefall time. Using this, we produced the initial conditions of embedded clusters for our Nbody simulations, which depend on the product of two parameters (SFE per freefall time and SF duration). With the equations we provide in the Appendix A, one can use any centrally concentrated spherically symmetric density profile for the embedded cluster and recover the initial and residual gas density profiles.
We adapted the falcON program mkhalo by McMillan & Dehnen (2007) to our problem and have written an additional acceleration plugin, which represents the gravitational potential of the residual gas in dependence on the SF duration and SFE per freefall time for a Plummer embedded cluster. Building on this adapted version of mkhalo, we produced the initial conditions of our simulations, that is, a Plummer star cluster in virial equilibrium with its residual gas with their respective density profiles obeying Eqs. (18) and (19) of Parmentier & Pfalzner (2013). We related the SF duration to the global SFE to make our set of simulations comparable to earlier works, in which the bound fraction is often presented in dependence on the global SFE.
We performed two types of cluster simulations, each time covering different global SFEs: 1) isolated singlemass clusters; and 2) nonisolated models, that is, star clusters with stellar evolution and dissolving within a Milky Waylike galaxy. We studied the effect of different initial cluster stellar masses as well as of different cluster densities on the evolution of our nonisolated models. The latter are implemented by varying the cluster halfmass to tidalradius ratio, r_{h}/r_{t} = 0.025, 0.052, 0.07, and0.10.
Based on the performed simulations, we quantified the bound fraction evolution and the violent relaxation duration of young clusters. We defined the violent relaxation duration as the time span of cluster massloss in response to instantaneous gas expulsion. We note that with our definition, the violent relaxation duration differs from the cluster revirialization time. Our models for isolated singlemass clusters allowed us to define an upper limit to the bound fraction as a function of the global SFE. For the models considered in our work with their specific parameters – the stellar density profile (Plummer model), the cluster orbit (with a circular velocity, in the solar neighborhood in the Galactic disk plane), the stellar evolutionary massloss from the SSE routine (Hurley et al. 2000) and for the cluster mean volume density range (r_{h}/r_{t} = [0.025:0.100]) – we conclude that the violent relaxation is not longer than 20 Myr, and its duration depends weakly on the initial stellar density of a cluster. We found that the violent relaxation duration of nonisolated model clusters depends neither on the cluster initial stellar mass nor on the global SFE, keeping the same initial stellar density. Varying the tidal field impact, that is, varying the cluster size while retaining the cluster mass, does not affect the cluster massloss in response to instantaneous gas expulsion much.
We also found that the minimum global SFE necessary to form a bound cluster after instantaneous gas expulsion is SFE_{gl} = 0.15 for a cluster with a circular orbit in the Galactic disk plane at a distance of R_{0} = 8 kpc from the Galactic center. If the tidal field is stronger, that is, the cluster is closer to the Galactic center, the minimum SFE_{gl} needed to survive instantaneous gas expulsion may be different. For the given r_{h}/r_{t} ratio, the bound fraction of surviving clusters that achieved the same global SFE does not depend on the cluster initial stellar mass. The bound mass fraction at the end of violent relaxation for clusters with r_{h}/r_{t} of 0.025 and 0.10 differs by only about 10 percent, with denser clusters retaining a higher bound fraction than more diffuse clusters. The evolution of bound clusters after violent relaxation is mostly driven by the tidal field of the host galaxy, and their life expectancy then depends on their stellar mass.
We compared our results with earlier works. Our final bound fractions are similar to those found in previous works only when the bound fraction is plotted in dependence of the effective SFE. Thus we agree with Goodwin (2009) that the virial ratio of a cluster at the time of gas expulsion is a key parameter for predicting whether it survives gas expulsion. However, when working in terms of the global SFE, that is, the SFE that can be measured by observers as the ratio between the stellar mass and the total (gas+star) mass of a starforming region, the models proposed in this paper improve the survival likelihood of star clusters after instantaneous gas expulsion. This is caused by the difference in density profiles between the embedded cluster and its residual gas, namely, the stellar density profile has a steeper slope than that of the residual gas, which is a consequence of SF taking place with a constant starformation efficiency per freefall time in a centrally concentrated molecular clump.
Acknowledgments
We thank the referee for the helpful comments. This work was supported by Sonderforschungsbereich SFB 881 “The Milky Way System” (subproject B2) of the German Research Foundation (DFG). The authors gratefully acknowledge Walter Dehnen for his support and discussions in connection with implementing the code mkhalo for our purposes. The authors gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JURECA at Jülich Supercomputing Centre (JSC), under application No. 10341 (hhd28). We also used the smaller GPU cluster KEPLER, funded under the grants I/80 041043 and I/81 396 of the Volkswagen Foundation, and the special GPU accelerated supercomputer LAOHU at the Center of Information and Computing at National Astronomical Observatories, Chinese Academy of Sciences, funded by the Ministry of Finance of the People’s Republic of China under the grant ZDYZ20082. B.S. gratefully acknowledges Rainer Spurzem for providing computing time on the highperformance computing clusters JURECA and LAOHU.
B.S. acknowledges the support within program N007310/PCF15MON by the Ministry of Education and Science of the Republic of Kazakhstan.
P.B. acknowledges the support by the Chinese Academy of Sciences through the Silk Road Project at NAOC, through the “Qianren” special foreign experts program, and the President’s International Fellowship for Visiting Scientists program of CAS and also the Strategic Priority Research Program (Pilot B) “Multiwavelength gravitational wave universe” of the Chinese Academy of Sciences (No. XDB23040100).
P.B. acknowledges the support of the Volkswagen Foundation under the Trilateral Partnerships grant No. 90411 and the special support by the NASU under the Main Astronomical Observatory GRID/GPU computing cluster project.
References
 Adams, F. C. 2000, ApJ, 542, 964 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S., & Kroupa, P. 2013, ApJ, 764, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Baumgardt, H., & Kroupa, P. 2007, MNRAS, 380, 1589 [NASA ADS] [CrossRef] [Google Scholar]
 Berczik, P., Nitadori, K., Zhong, S., et al. 2011, in International conference on High Performance Computing, Kyiv, Ukraine, October 8–10, 8 [Google Scholar]
 Boily, C. M., & Kroupa, P. 2003a, MNRAS, 338, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Boily, C. M., & Kroupa, P. 2003b, MNRAS, 338, 673 [NASA ADS] [CrossRef] [Google Scholar]
 Bonnell, I. A., Smith, R. J., Clark, P. C., & Bate, M. R. 2011, MNRAS, 410, 2339 [NASA ADS] [CrossRef] [Google Scholar]
 Brinkmann, N., Banerjee, S., Motwani, B., & Kroupa, P. 2017, A&A, 600, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dib, S., Gutkin, J., Brandner, W., & Basu, S. 2013, MNRAS, 436, 3727 [NASA ADS] [CrossRef] [Google Scholar]
 Duboshin, G. N. 1968, Nebesnaia mekhanika, Osnovnuye zadachi i metody (Moscow: Nauka) [Google Scholar]
 Ernst, A., Berczik, P., Just, A., & Noel, T. 2015, Astron. Nachr., 336, 577 [NASA ADS] [CrossRef] [Google Scholar]
 Evans, II, N. J., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Farias, J. P., Smith, R., Fellhauer, M., et al. 2015, MNRAS, 450, 2451 [NASA ADS] [CrossRef] [Google Scholar]
 Fellhauer, M., & Kroupa, P. 2005, ApJ, 630, 879 [NASA ADS] [CrossRef] [Google Scholar]
 Fujii, M. S., & Portegies Zwart, S. 2016, ApJ, 817, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Geyer, M. P., & Burkert, A. 2001, MNRAS, 323, 988 [NASA ADS] [CrossRef] [Google Scholar]
 Girichidis, P., Federrath, C., Allison, R., Banerjee, R., & Klessen, R. S. 2012, MNRAS, 420, 3264 [NASA ADS] [CrossRef] [Google Scholar]
 Goodwin, S. P. 2009, Ap&SS, 324, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Goodwin, S. P., & Bastian, N. 2006, MNRAS, 373, 752 [NASA ADS] [CrossRef] [Google Scholar]
 Higuchi, A. E., Kurono, Y., Saito, M., & Kawabe, R. 2009, ApJ, 705, 468 [NASA ADS] [CrossRef] [Google Scholar]
 Hillenbrand, L. A., & Carpenter, J. M. 2000, ApJ, 540, 236 [NASA ADS] [CrossRef] [Google Scholar]
 Hills, J. G. 1980, ApJ, 235, 986 [NASA ADS] [CrossRef] [Google Scholar]
 Hopkins, P. F., Narayanan, D., Murray, N., & Quataert, E. 2013, MNRAS, 433, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Just, A., Berczik, P., Petrov, M. I., & Ernst, A. 2009, MNRAS, 392, 969 [NASA ADS] [CrossRef] [Google Scholar]
 Kainulainen, J., Federrath, C., & Henning, T. 2014, Science, 344, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Krumholz, M. R., & Matzner, C. D. 2009, ApJ, 703, 1352 [NASA ADS] [CrossRef] [Google Scholar]
 Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304 [NASA ADS] [CrossRef] [Google Scholar]
 Kudryavtseva, N., Brandner, W., Gennaro, M., et al. 2012, ApJ, 750, L44 [NASA ADS] [CrossRef] [Google Scholar]
 Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Lada, C. J., Margulis, M., & Dearborn, D. 1984, ApJ, 285, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Lada, E. A., Depoy, D. L., Evans, II, N. J., & Gatley, I. 1991, ApJ, 371, 171 [NASA ADS] [CrossRef] [Google Scholar]
 Lada, C. J., Lombardi, M., & Alves, J. F. 2010, ApJ, 724, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, P. L., & Goodwin, S. P. 2016, MNRAS, 460, 2997 [NASA ADS] [CrossRef] [Google Scholar]
 Leisawitz, D., Bash, F. N., & Thaddeus, P. 1989, ApJS, 70, 731 [NASA ADS] [CrossRef] [Google Scholar]
 McMillan, P. J., & Dehnen, W. 2007, MNRAS, 378, 541 [NASA ADS] [CrossRef] [Google Scholar]
 Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
 Moeckel, N., Holland, C., Clarke, C. J., & Bonnell, I. A. 2012, MNRAS, 425, 450 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, N. 2011, ApJ, 729, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, N., Quataert, E., & Thompson, T. A. 2010, ApJ, 709, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Parmentier, G., & Baumgardt, H. 2012, MNRAS, 427, 1940 [NASA ADS] [CrossRef] [Google Scholar]
 Parmentier, G., & Gilmore, G. 2007, MNRAS, 377, 352 [NASA ADS] [CrossRef] [Google Scholar]
 Parmentier, G., & Pfalzner, S. 2013, A&A, 549, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pfalzner, S., Parmentier, G., Steinhausen, M., Vincke, K., & Menten, K. 2014, ApJ, 794, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [NASA ADS] [CrossRef] [Google Scholar]
 Proszkow, E.M., & Adams, F. C. 2009, ApJS, 185, 486 [NASA ADS] [CrossRef] [Google Scholar]
 Reggiani, M., Robberto, M., Da Rio, N., et al. 2011, A&A, 534, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smith, R., Fellhauer, M., Goodwin, S., & Assmann, P. 2011, MNRAS, 414, 3036 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R., Goodwin, S., Fellhauer, M., & Assmann, P. 2013, MNRAS, 428, 1303 [NASA ADS] [CrossRef] [Google Scholar]
 Tutukov, A. V. 1978, A&A, 70, 57 [NASA ADS] [Google Scholar]
 Verschueren, W., & David, M. 1989, A&A, 219, 105 [NASA ADS] [Google Scholar]
Appendix A: Expression of the density profile of unprocessed gas
Equation (10) can be easily solved using software such as mathematica. Since the roots of this equation obtained with mathematica are very long, we made them more compact by introducing the following intermediate terms: Then we can write the four roots of Eq. (10) as (A.5)The following relation is true for all real values of k and ρ_{⋆}(A.6)which gives complex numbers for two of the roots in case of (− K_{1}). The other two roots are real, with one decreasing and the other increasing with ρ_{⋆}.
In our case, we expect the residual gas density profile to be a decreasing function of the radius, as the stellar density profile. Thus we choose the root that increases together with stellar density toward the clump center with the following expression: (A.7)
Appendix B: Bound fraction of isolated models
For isolated models we first defined the bound fraction based on the total (i.e., kinetic + potential) energy of stars, that is, as the fraction of stars with a negative total energy (solid lines in Fig. B.1).
Fig. B.1 Time evolution of the bound fraction F_{b} of isolated models (N = 10^{4}) as defined by two methods: defined by the fraction of stars with a negative total energy (solid lines), and defined by recalculating the total energy of stars in an iterative process (see text for details; dashed lines). The vertical dotted line corresponds to t = 20 Myr when we scale the isolated models with the same scale factor as for a nonisolated model with M_{⋆} = 6000 M_{⊙}, which also has N ≈ 10^{4}. 

Open with DEXTER 
Fig. B.2 Top panel: bound and final bound fraction evolution of an isolated cluster with SFE_{gl} = 0.15 (solid and dashed lines, respectively). Bottom panel: Lagrange radius evolution of the same model. Lagrange radii are given in units of the initial stellar halfmass radius. 

Open with DEXTER 
Since our model clusters become supervirial after instantaneous gas expulsion, the unbound stars can be located anywhere inside the cluster and make a significant contribution to its gravitational potential depending on global SFE. Thus the bound fraction remains overestimated until bound and unbound stars are clearly spatially separated from each other. For instance, model clusters with global SFEs of 0.05 and 0.10 do not survive the instantaneous gas expulsion. In Fig. B.1, however, they retain a significant bound fraction for several Myr. For clusters with global SFEs of 0.13 and 0.15, we need to wait for a long time to reach the final bound fraction, that is, for the unbound stars to have evacuated the cluster region and to no longer contribute to the cluster gravitational potential.
These reasons motivated us to develop another technique to define the final bound fraction early on in the evolution of clusters. To do so, we recalculated the total energy of each star after removing the unbound stars. Some of these recalculated energies are now positive (i.e., some previously bound stars are now unbound) since the cluster gravitational potential is now shallower. We removed this new sample of unbound stars, and continued to iterate until only bound stars were left within our selection. The fraction of stars left in this selection determines the final bound fraction.
Figure B.2 presents an example of the evolution of the bound fraction and Lagrange radii of an isolated cluster with a global SFE of 15 percent. As we see when comparing the top and bottom panels, the bound fraction decreases during the cluster expansion and reaches its final value when its bound stars have collapsed back and form a bound cluster. With the new technique, however, we can predict the final bound fraction already after 2 Myr, when the cluster is still in the expansion phase. Figure B.1 and the top panel of Fig. B.2 show that the instantaneous bound fraction converges toward the final bound fraction determined with our technique by the end of the simulations. This shows that with our calculation method we can estimate the final bound fraction even before the inner part of the cluster starts to collapse back and return to virial equilibrium. We caution, however, that with this method we underestimate the final bound fraction of a cluster with a low global SFE during their early evolution after instantaneous gas expulsion. This is caused by removing all unbound stars, including the centrally concentrated ones, which contribute the most to the gravitational field of the cluster (see Fig. B.2). This is the reason for the unusual behavior of the final bound fractions of isolated clusters with a global SFE of 0.13 and 0.15, which is 0 at t ≲ 1 Myr, and why they rise at an early time in the evolution instead of decreasing.
All Tables
Scaling factors to the physical time and duration of 20 Myr in Nbody units corresponding to different ratios of the cluster halfmass to tidal radius.
All Figures
Fig. 1 Relation between the SF duration (t_{SF}) (given in Nbody units according to Eq. (12)) and the global SFE (SFE_{gl}), measured inside R_{cl} = 10a_{⋆}. The lowest value of the global SFE here in this plot is SFE_{gl} = 0.007, which corresponds to an SF duration t_{SF} = 0.01 [NB]. 

Open with DEXTER  
In the text 
Fig. 2 Density profiles of the star cluster (black dashed line), of the residual (solid lines), and initial (dashdotted lines) gas for different SFE_{gl} in scaled physical units. A total stellar mass M_{⋆} = 3000 M_{⊙} and a 3D halfmass radius r_{h} = 1 pc are assumed. Note that the stellar density profile is a Plummer profile. 

Open with DEXTER  
In the text 
Fig. 3 Evolution of the bound fraction of clusters obtained in our simulations. Different colors correspond to different global SFEs (see the key). The gray lines show the impact of mass loss caused by stellar evolution alone. The vertical dotted line corresponds to t = 20 Myr.Top panel: comparison of the bound fraction evolution of isolated (dashed lines) and nonisolated (solid lines) models for M_{⋆} = 6000 M_{⊙} with N ≈ 10^{4} stars. Note that the isolated models are the singlemass models without stellar evolution and scaled to the same physical units as the nonisolated models (i.e., M_{⋆} = 6000 M_{⊙} and r_{h} = 1.26 pc). Bottom panel: bound fraction (instantaneous tidal mass as a fraction of initial stellar mass) evolution for nonisolated models with different initial stellar masses (see the key for the linecoding). 

Open with DEXTER  
In the text 
Fig. 4 Different random realizations (see line types in the key) of the model clusters with M_{⋆} = 6000 M_{⊙} (N = 10 455 stars) for r_{h}/r_{t} = 0.052. 

Open with DEXTER  
In the text 
Fig. 5 Evolution of a nonisolated cluster of the solar neighborhood with SFE_{gl} = 0.20, M_{⋆} = 6000 M_{⊙} and r_{h}/r_{t} = 0.052 over 50 Myr after gas expulsion. We show 49 Lagrange radii, ranging from 2 percent to 98 percent in intervals of 2 percent (solid and dashed lines). The instantaneous tidal radius r_{t} of the cluster is shown as the dotted line. Thick solid lines correspond to every 10 percent of the Lagrange radius, and the dashed line corresponds to 50 percent of the Lagrange radius. 

Open with DEXTER  
In the text 
Fig. 6 Evolution of the bound fraction of M_{⋆} = 15 000 M_{⊙} clusters with different impact of the tidal field (see the key for the linecoding). Different colors correspond to different global SFEs (see the key for the colorcoding). 

Open with DEXTER  
In the text 
Fig. 7 Evolution of the cluster with M_{⋆} = 6000 M_{⊙}, SFE_{gl} = 0.20 and r_{h}/r_{t} = 0.025 over the first 50 Myr after the instantaneous gas expulsion. 

Open with DEXTER  
In the text 
Fig. 8 Bound fraction as a function of global or effective SFE. We compare our results (red lines) with previous works. The isolated models are depicted by the red diamonds, and nonisolated models by red crosses. In the top panel we use our global SFE (SFE_{gl}) and in the bottom panel, the effective SFE, eSFE = 1 /Q_{⋆} (which is almost the same as the LSF, the local stellar fraction of our models). 

Open with DEXTER  
In the text 
Fig. B.1 Time evolution of the bound fraction F_{b} of isolated models (N = 10^{4}) as defined by two methods: defined by the fraction of stars with a negative total energy (solid lines), and defined by recalculating the total energy of stars in an iterative process (see text for details; dashed lines). The vertical dotted line corresponds to t = 20 Myr when we scale the isolated models with the same scale factor as for a nonisolated model with M_{⋆} = 6000 M_{⊙}, which also has N ≈ 10^{4}. 

Open with DEXTER  
In the text 
Fig. B.2 Top panel: bound and final bound fraction evolution of an isolated cluster with SFE_{gl} = 0.15 (solid and dashed lines, respectively). Bottom panel: Lagrange radius evolution of the same model. Lagrange radii are given in units of the initial stellar halfmass radius. 

Open with DEXTER  
In the text 