Disruption of giant molecular clouds and formation of bound star clusters under the influence of momentum stellar feedback

TL;DR

This study uses simulations to explore how stellar feedback influences giant molecular clouds and the formation of bound star clusters, revealing key dependencies on initial conditions and feedback strength.

Contribution

It introduces a detailed simulation framework that links cloud properties and feedback to star formation efficiency and cluster boundness, advancing understanding of cluster formation.

Findings

Star formation efficiency scales with initial surface density and feedback strength.

Stellar density profiles are steeper than gas profiles during peak star formation.

Clusters remain bound even with star formation efficiencies below 50%.

Abstract

Energetic feedback from star clusters plays a pivotal role in shaping the dynamical evolution of giant molecular clouds (GMCs). To study the effects of stellar feedback on the star formation efficiency of the clouds and the dynamical response of embedded star clusters, we perform a suite of isolated GMC simulations with star formation and momentum feedback subgrid models using the moving-mesh hydrodynamics code \textsc{Arepo}. The properties of our simulated GMCs span a wide range of initial mass, radius, and velocity configurations. We find that the ratio of the final stellar mass to the total cloud mass, $\epsilon_{\rm int}$, scales strongly with the initial cloud surface density and momentum feedback strength. This correlation is explained by an analytic model that considers force balancing between gravity and momentum feedback. For all simulated GMCs, the stellar density profiles are systematically steeper than that of the gas at the epochs of the peaks of star formation, suggesting a centrally concentrated stellar distribution. We also find that star clusters are always in a sub-virial state with a virial parameter $\sim0.6$ prior to gas expulsion. Both the sub-virial dynamical state and steeper stellar density profiles prevent clusters from dispersal during the gas removal phase of their evolution. The final cluster bound fraction is a continuously increasing function of $\epsilon_{\rm int}$. GMCs with star formation efficiency smaller than 0.5 are still able to form clusters with large bound fractions.

Tables

Table 1. Table 1: Model Parameters

Name	$M_{\mathrm{GMC}}\left(M_{\odot }\right)$	$R$ (pc)	${\alpha }_0^R$	${\tau }_{\mathrm{ff}}$ (Myr)	$m_{\mathrm{res}}$ ($M_{\odot }$)	$l_{\mathrm{soft}}$ (${10}^{-3}$pc)
(i)	(ii)	(iii)	(iv)	(v)	(vi)	(vii)
RHO5 T/R	$1.25\times {10}^4$	5	0.1/0.9	1.6	$6.0\times {10}^{-3}$	1
RHO10	${10}^5$	10	0.1/0.9	1.6	$4.8\times {10}^{-2}$	2
RHO20	$8\times {10}^5$	20	0.1/0.9	1.6	$0.38$	4
RHO40	$6.4\times {10}^6$	40	0.1/0.9	1.6	$3.1$	8
RHO80	$5.12\times {10}^7$	80	0.1/0.9	1.6	$24.4$	16
SIGMA5	$2.5\times {10}^4$	5	0.1/0.9	1.13	$1.2\times {10}^{-2}$	1
SIGMA20	$4\times {10}^5$	20	0.1/0.9	2.26	$0.19$	4
SIGMA40	$1.6\times {10}^6$	40	0.1/0.9	3.2	$0.76$	8

Equations

$$\frac{{\dot{m}}_{\mathrm{w}}}{\mathrm{Gy}{\mathrm{r}}^{-1}}=\{\begin{array}{ll}4.763(0.01+Z/Z_{\odot })& t_6<1\\ 4.763(0.01+Z/Z_{\odot })t_6^{1.45+0.8\ln (Z/Z_{\odot })}& 1<t_6<3.5\\ 29.4(t_6/3.5{)}^{-3.25}+0.0042& 3.5<t_6<100\end{array}$$

$$l_{\mathrm{w}}=\left[\frac{5.94\times 10^4}{1+(t_6/2.5{)}^{1.4}+(t_6/10{)}^5}+4.83\right]\times 10^{12}{\dot{m}}_{\mathrm{w}}erg/g,$$

$$v_x^R(x,y,z)=-{\Omega }_{i}y;v_y^R(x,y,z)={\Omega }_{i}x;v_z^R(x,y,z)=0,$$

$${\mathcal{H}}_{ij}(\mathbf{k})={\eta }_{\mathrm{comp}}{\mathcal{H}}_{ij}^{\parallel }(\mathbf{k})+(1-{\eta }_{\mathrm{comp}}){\mathcal{H}}_{ij}^{\perp }(\mathbf{k}),$$

$$\mathrm{LSF}=\frac{M_{\ast }(<r_{\mathrm{h}})}{M_{\ast }(<r_{\mathrm{h}})+M_{\mathrm{gas}}(<r_{\mathrm{h}})}.$$

$$F_{\mathrm{sh},\ast }=\frac{G{M}_{\ast }{M}_{\mathrm{sh}}}{r_s^2{A}_s}=\frac{G{M}_{\ast }{\Sigma }_{\mathrm{sh}}}{r_s^2},$$

$$F_{\mathrm{sh}}=\frac{\beta G{M}_{\mathrm{sh}}^2}{r_s^2{A}_s}=\frac{\beta G{M}_{\mathrm{sh}}{\Sigma }_{\mathrm{sh}}}{r_s^2},$$

$$F_p=\frac{M_{\ast }\dot{p}}{4\pi r_s^2},$$

$$\frac{M_{\ast }\dot{p}}{4\pi r_s^2}=\frac{G{M}_{\ast }{\Sigma }_{\mathrm{sh}}}{r_s^2}+\frac{\beta G{M}_{\mathrm{sh}}{\Sigma }_{\mathrm{sh}}}{r_s^2}.$$

$$(1-{ϵ}_{\mathrm{int}})\left(\frac{\beta }{{ϵ}_{\mathrm{int}}}-\beta +1\right)=\frac{1}{\Gamma }.$$

$${ϵ}_{\mathrm{int}}=\frac{\sqrt{{\Gamma }^2+(4\beta -2)\Gamma +1}-(2\beta -1)\Gamma -1}{2(1-\beta )\Gamma }.$$

$${ϵ}_{\mathrm{int}}=\frac{\sqrt{{\Gamma }^2+1}-1}{\Gamma }.$$

$${\tau }_{\mathrm{ff}}=\sqrt{\frac{3\pi }{32G{\rho }_0}}\approx 1.6\mathrm{Myr}{\left(\frac{M}{10^5{M}_{\odot }}\right)}^{-1/2}{\left(\frac{R_{\mathrm{GMC}}}{10\mathrm{pc}}\right)}^{3/2}.$$

$${\tau }_{\mathrm{spread}}\equiv \frac{M_{\ast }}{<\dot{M}>}=\frac{M_{\ast }}{\int {\dot{M}}_{\ast }^{2}dt/M_{\ast }},$$

$$f(v)dv\equiv f(x)dx=\sqrt{\frac{2}{\pi }}{x}^2\exp (-x^2/2)dx,$$

$$v_{\mathrm{esc}}=\sqrt{-\frac{2{\Omega }_{\ast ,1}}{M_{\ast }}}=\sqrt{\frac{2{ϵ}_{\mathrm{int}}}{{\alpha }_{\ast }}}{v}_{\mathrm{rms}}.$$

$$f_{\mathrm{bound}}$$

Full text

††pubyear: 2018††pagerange: Disruption of giant molecular clouds and formation of bound star clusters under the influence of momentum stellar feedback–Disruption of giant molecular clouds and formation of bound star clusters under the influence of momentum stellar feedback

1 Introduction

Most, if not all, stars are formed in clusters (Lada & Lada, 2003), which emerged from giant molecular clouds (GMC; Shu et al., 1987; Scoville & Good, 1989; McKee & Ostriker, 2007; Krumholz et al., 2018). Due to the complex interplay of gravity, supersonic turbulence and stellar feedback from massive stars, the dynamical evolution and cluster formation activities in GMCs are still highly debatable (e.g. Krumholz & McKee, 2005; Ballesteros-Paredes & Hartmann, 2007; Heitsch et al., 2009; Hennebelle & Chabrier, 2011; Padoan & Nordlund, 2011; Federrath & Klessen, 2012; Burkert & Hartmann, 2013; Traficante et al., 2018). One of the key observables that can be used to constrain various physical process is the star formation efficiency (SFE) in star-forming regions.

It is well known that star formation is inefficient on galactic scales. The observed linear correlation between molecular gas surface density and star formation rate (SFR) surface density in normal star-forming galaxies suggests a nearly constant gas depletion time-scale around $\sim$ 2 Gyr, much longer than the dynamical time-scale of galactic disks (Kennicutt, 1989; Bigiel et al., 2008; Saintonge et al., 2011; Leroy et al., 2013; Genzel et al., 2015; Tacconi et al., 2018). In contrast, the SFE on GMC scales shows a large variation ranging from less than a few percent to nearly unity (Zuckerman & Evans, 1974; Krumholz & Tan, 2007; Wu et al., 2010; Evans et al., 2014; Heyer et al., 2016; Lee et al., 2016; Vutisalchavakul et al., 2016). The origin of this large scatter is usually explained as a combination of the time variability of the SFR during the course of cloud evolution and intrinsic scatter of SFEs due to the diversity of GMC properties (Feldmann & Gnedin, 2011; Kruijssen & Longmore, 2014; Lee et al., 2016; Grudić et al., 2018a; Kruijssen et al., 2018). For example, recent theoretical models and high-resolution magneto-hydrodynamics simulations suggest that the SFE depends on the local virial parameters of the cloud controlled by large-scale turbulence (e.g. Krumholz & McKee, 2005; Padoan et al., 2012). However, it has recently been recognized that large-scale turbulence can only account for an $\sim$ 0.3 dex scatter, which is not sufficient to explain the observed SFE variations (Lee et al., 2016). Another source of variation comes from different stellar feedback channels that alter the dynamical states of the GMCs (Fall et al., 2010; Murray et al., 2010; Dale et al., 2014; Myers et al., 2014; Raskutti et al., 2016; Kim et al., 2017; Grudić et al., 2018b). Previous studies found that GMC simulations adopting different stellar feedback mechanisms (stellar winds, ionizing radiation, or supernovae) lead to dramatically different final SFEs. The problem has recently been recognized to be more subtle than previously thought, since even small differences in numerical treatments, such as different radiative transfer schemes, massive star sampling, and momentum and energy deposition algorithms, can lead to drastic changes for the final SFE (Dale et al., 2005; Roškar et al., 2014; Raskutti et al., 2016; Grudić et al., 2018b; Kim et al., 2018). Therefore, how the SFE depends on GMC properties and the strength of stellar feedback remains an open question.

Stellar feedback not only changes the efficiency of star formation within GMCs, but also alters the dynamical state of star clusters by dispersing the cloud. The process of star cluster disruption due to rapid gas expulsion shortly after the cluster emerges from its natal cloud is believed to be the main culprit of the “infant mortality” of star clusters – a sharp decrease in the number of young star clusters with the increase of cluster age in local star-forming regions (e.g. Lada & Lada, 2003). A simple virial analysis suggests that a star cluster that is initially in virial equilibrium will dissociate if more than half of the mass is instantaneously lost (Hills, 1980; Mathieu, 1983). However, this statement does not take into account the highly non-linear star formation and stellar feedback process in realistic self-gravitating turbulent environments. For example, embedded clusters in star-forming regions are not necessarily in virial equilibrium. Recent hydrodynamical simulations suggest that the stellar velocity dispersions are in general much smaller than that of the gas, suggesting a sub-virial dynamical state of star clusters within GMCs (e.g. Offner et al., 2009). Moreover, stars are usually not well mixed with gas but instead formed in the densest part of the cloud. The difference between the gas and stellar distribution can strongly affect the dynamical response of star clusters to gas dispersal (e.g. Kruijssen et al., 2012a; Shukirgaliyev et al., 2018). Most importantly, GMCs are highly substructured. Stars are formed at the intersections of gas filaments and assembled into different subclusters hierarchically. Previous works have explored some of the above complications using different physical and numerical methods, including analytical models (Hills, 1980; Mathieu, 1983; Adams, 2000; Boily & Kroupa, 2003a; Kruijssen, 2012; Parmentier & Pfalzner, 2013), pure N-body simulations (Tutukov, 1978; Lada et al., 1984; Boily & Kroupa, 2003b; Goodwin & Bastian, 2006; Baumgardt & Kroupa, 2007; Smith et al., 2011; Farias et al., 2018), and hydrodynamic simulations (Bonnell et al., 2011; Girichidis et al., 2012; Moeckel et al., 2012; Fujii & Portegies Zwart, 2016; Gavagnin et al., 2017). Recent efforts have been made to include various relevant physical processes in hydrodynamic simulations (e.g. Parker et al., 2015; Gavagnin et al., 2017), however, due to the high computational costs, they usually focus on a handful of less massive GMCs.

In this paper, we perform a suite of hydrodynamic simulations of turbulent GMCs employed with a simple star formation and stellar feedback models in the moving-mesh code AREPO. We survey GMCs with a broad range of mass, size, and velocity configurations to investigate the physical origin of the intrinsic variations of SFEs and the properties of surviving star clusters. The structure of this paper is as follows. In Section 2, we describe the simulation setup, star formation and momentum stellar feedback implementations, and the design of initial conditions. In Section 3, we examine the dependence of integrated SFE of GMCs on cloud mass, size, and momentum feedback intensity. In Section 4, we describe the subsequent dynamical evolution of star clusters after the residual gas is completely expelled and investigate the relationship between SFE and cluster bound fraction. We summarize our conclusions in Section 5.

2 Methods

2.1 Simulation setup

The simulations in this work are performed with Arepo (Springel, 2010), a moving-mesh, finite-volume hydrodynamic code employing a second-order unsplit Godunov scheme. The control volumes are discretized by a Voronoi tessellation, which is generated from its dual Delaunay tessellation determined by a set of mesh-generating points. These points can move freely within the simulation domain and follow gas flows in a quasi-Lagrangian fashion. Therefore, Arepo captures the advantages of both grid- and particle-based hydrodynamic methods and has already been applied to various astrophysical problems (e.g. Kereš et al., 2012; Torrey et al., 2012; Vogelsberger et al., 2012, 2014; Springel et al., 2018). Our simulations include hydrodynamics, self-gravity, radiative cooling, star formation, and momentum feedback from stellar winds.

We use an adaptive softening scheme for gas cells so that the gravitational forces are resolved all the way down to the size of each cell. We employ a quasi-Lagrangian refinement scheme that keeps the mass of gas cells close to a target mass determined by initial conditions. In addition, we refine a cell if its volume is more than 32 times larger than the minimum volume of all its face-touching neighbours. This volume-limited refinement scheme prevents large volume contrast between adjacent cells and helps to better resolve the regions that experience fast expansion due to stellar feedback.

Since radiative cooling is responsible for gas fragmentation and subsequent star formation in GMCs, it is important to follow the cooling process explicitly over a large range of temperatures. Instead of adopting isothermal or effective equations of state, which have been used in many previous GMC simulations (e.g. Dale et al., 2005; Raskutti et al., 2016; Kim et al., 2018), we explicitly include radiative cooling from multiple channels: a network implementing hydrogen and helium cooling and heating processes due to collisions, recombinations, free-free emission and photoionization from UV background radiation; high-temperature ($T>10^{4}\,{\rm K}$) metal-line cooling that is added to the hydrogen and helium network following Vogelsberger et al. (2013); low-temperature ($T<10^{4}\,{\rm K}$), metal-line, fine-structure and molecular cooling implemented as a fitting function, depending on temperature, density and gas metallicity, to CLOUDY calculations presented in Hopkins et al. (2018) and Marinacci et al. (in prep.). We set the metallicity of the GMCs to solar abundance when evaluating the cooling rates from metals. One caveat is that the adopted cooling curves are based on the CLOUDY model calculations under the spatially uniform UV background used for galaxy formation simulations. The cooling rate calculated based on this uniform UV background may not be accurate for GMC simulations, especially in close proximity of massive stars.

2.2 Star formation

During each simulation time-step, we identify all star-forming cells and convert them to stellar particles probabilistically. Star-forming cells are defined as gas cells that are cold ($T_{\rm cell}<100$ K), contracting ($\nabla\cdot\bm{v}<0$), and self-gravitating ($|\nabla\cdot\bm{v}|^{2}+|\nabla\times\bm{v}|^{2}<2G\rho$), where $\bm{v}$ and $\rho$ are the velocity and density of the cells, respectively. We also employ a density threshold for star formation, $n_{\rm cell}>10^{5}{\rm cm^{-3}}$, to avoid rare situations where some self-gravitating clumps are formed in the very low density outskirt of the cloud.

A given star-forming cell is converted to stellar particles with a constant probability $p=\Delta t/\tau_{\mathrm{ff}}(\rho)$ at a given time-step $\Delta t$, where $\tau_{\mathrm{ff}}=(3\pi/32G\rho)^{1/2}$ is the free-fall time of the cell. The cells that are converted to stellar particles are removed and the volume of these cells is claimed by their neighbours. The mass, position, and velocity of the newly formed stellar particles are inherited from their parent gas cells. Therefore, the mass distribution of stars is similar to that of the gas particles, which is around the target mass of the simulations. After the stellar particles are created, they are treated as collisionless particles with a Plummer-equivalent softening length fixed to $10^{-4}$ of the initial diameter of the GMC.

2.3 Momentum stellar feedback

The overall evolution of GMCs depends strongly on the strength of stellar feedback. Unfortunately, the exact amount of feedback that is associated with massive stars in the simulations is still debated. As has already been noticed in previous studies, GMC simulations with different stellar feedback sources (stellar winds, ionizing radiation, or supernovae) show dramatically different gas evolution and star formation efficiencies (Dale et al., 2005; Roškar et al., 2014; Raskutti et al., 2016; Grudić & Hopkins, 2018; Kim et al., 2018). Moreover, it has recently been recognized that, even some small changes in numerical implementations, such as radiation hydrodynamic methods (Kim et al., 2018), sampling of massive star formation (Grudić et al., 2018a), and momentum/energy deposition algorithms (Hopkins et al., 2018), can contribute noticeable variation to the star formation efficiencies of the clouds. Since exploring accurate feedback implementation from various sources is not the main focus of this paper, we simply treat stellar feedback by depositing mass and momentum fluxes from stellar particles to their neighbouring gas cells.

We set the fiducial mass-loss and momentum deposition rate to the initial mass function (IMF)-averaged values of stellar winds from a single stellar population with a Kroupa initial mass function (Kroupa, 2001). Following Hopkins et al. (2018), the mass-loss rate per unit stellar mass is

[TABLE]

where $Z$ is the metallicity and $t_{6}$ is the age of stellar particles in unit of Myr. The kinetic luminosity of winds per unit stellar mass is

[TABLE]

for $t_{6}<100$. Winds from stellar particles older than $t_{6}>100$ are irrelevant here since the dynamical time-scales of our model GMCs are much shorter than 100 Myr.

The mass-loss and wind momentum are deposited to the gas cells around each stellar particle in the following way: for a given stellar particle of mass $m_{*}$ at time-step $\Delta t$, the total mass-loss is $\Delta m=\dot{m}_{\rm w}m_{*}\Delta t$ and wind momentum is $\Delta p=\sqrt{2l_{\rm w}\dot{m}_{\rm w}}f_{\mathrm{boost}}m_{*}\Delta t$, where $f_{\mathrm{boost}}$ is a boosting factor to the fiducial wind momentum to mimic the feedback intensity from different feedback sources. The mass and momentum fluxes from a stellar particle are distributed to its nearest 32 neighbouring gas cells in a weighted fashion so that cell $i$ with weight $w_{i}$ receives mass $\Delta m_{i}=(w_{i}/\Sigma_{j}{w_{j}})\Delta m$ and momentum $\Delta\bm{p}_{i}=(w_{i}/\Sigma_{j}{w_{j}})\Delta p\bm{r}_{i}/|\bm{r}_{i}|$, where $\bm{r}_{i}$ is the vector from the position of the stellar particle to the mesh-generating point of cell $i$. The weight can be chosen to be any physical quantities of the cells, such as volume, mass, or solid angle opened to the stellar particle. To test the robustness of the feedback implementation, we perform a series of numerical tests of wind-blowing bubbles created by a stellar particle with constant mass-loss rate $\dot{m}=10^{-5}\,M_{\odot}/$yr and wind velocity $v_{\rm w}=500$ km/s (see Appendix A). We find that the expansion history and the internal structure of the bubble are consistent with analytical solutions in Weaver et al. (1977). We also test the sensitivity of the star formation history of one GMC using different weighting methods (see Appendix B). To make consistent investigation across all GMCs, in the rest of the paper, we use solid angle as the weight to deposit mass and momentum.

2.4 Initial conditions

We set up the initial condition of GMCs as gas spheres of uniform density with initial turbulent velocity fields. The mass, radius, and other physical parameters of the initial conditions are listed in Table 1. We choose the initial mass ($M_{\mathrm{GMC}}$) and radius ($R_{\mathrm{GMC}}$) of the GMCs so that all “RHO” runs have the same initial volume density ($\rho_{0}=M_{\mathrm{GMC}}/(\frac{4}{3}\pi R_{\mathrm{GMC}}^{3})\approx 24\,M_{\odot}{\rm pc^{-3}}$) and all “SIGMA” runs have the same initial surface density ($\Sigma_{0}=M_{\mathrm{GMC}}/\pi R_{\mathrm{GMC}}^{2}\approx 318\,M_{\odot}{\rm pc^{-2}}$). The goal of this experimental design is to determine whether SFE depends on volume density or surface density.

The initial velocity field is initialized as a combination of turbulent motions and rigid rotation along the $z$-axis. We assign the rotational velocity field as

[TABLE]

where $v_{x}^{R},v_{y}^{R},v_{z}^{R}$ are the three components of the rotation velocities with circular frequency $\Omega_{i}$. For each run listed in Table 1, we construct two separate initial conditions: rotation-supported (“R”) and turbulent-supported (“T”) runs. The virial parameter contributed from rotation $\alpha_{0}^{R}\equiv 2E_{\rm rot}/|E_{G}|=0.1/0.9$ is used to calculate $\Omega_{i}$ for the corresponding “T”/“R” runs, where $E_{\rm rot}$ and $E_{G}$ are the rotational energy and gravitational energy of the cloud.

The turbulent velocity field is first initialized as a Gaussian random field in the Fourier space with the variance of the field determined by a given power spectrum, $P(k)$. Each dimension of the turbulent velocity field is treated independently and the result is a natural mixture of solenoidal and compressive turbulence. In order to rearrange the turbulent field into arbitrary solenoidal and compressive components, we perform a Helmholtz decomposition in $k$-space by applying the projection operator to the field (Federrath et al., 2010)

[TABLE]

where $\bm{\mathcal{H}}_{ij}^{\parallel}(\mathbf{k})=k_{i}k_{j}/k^{2}$ and $\bm{\mathcal{H}}_{ij}^{\bot}(\mathbf{k})=\bm{\delta}_{ij}-k_{i}k_{j}/k^{2}$ are the compressive and solenoidal operators, respectively, and $\bm{\delta}_{ij}$ is the Kronecker delta function. $\eta_{\rm comp}$ is the contribution from the compressive mode ranging from 0 to 1. $\eta_{\rm comp}=1$ means a purely compressive turbulence field and $\eta_{\rm comp}=0$ means a purely solenoidal. Varying $\eta_{\rm comp}$ would lead to a different density structure of the clouds. For simplicity, we use $\eta_{\rm comp}=1/3$ so that the ratio of kinetic energy of the compressive mode to that of the solenoid mode is 2:1 to mimic the natural mixture of the two turbulent modes. After projection, the turbulence field in $k$ space is Fourier-transformed to real space and is then interpolated to the position of gas cells within the sphere. The field is renormalized so that the virial parameter due to turbulence is $\alpha_{0}^{T}\equiv 2E_{\rm turb}/|E_{G}|$. We adopt a power-law power spectrum $P(k)\propto k^{-4}$, which is similar to the turbulence properties of GMCs (e.g. Dobbs et al., 2014). We assume the cloud is initially in virial equilibrium so that the virial parameter of the cloud $\alpha_{0}=\alpha_{0}^{R}+\alpha_{0}^{T}=1$.

The gas temperature is initialized to 10 K, which is commonly used for GMC simulations (Dale et al., 2014; Raskutti et al., 2016). The choice of the initial temperature does not change the evolution of the gas and star of the GMCs because of the short cooling time-scale compared to the dynamical time-scale of the clouds. For each initial condition, we perform five simulations with different momentum boosting factor, $f_{\mathrm{boost}}=$ 0.5, 1, 2, 4, 10, to test the effects of the strength of momentum feedback on the global evolution of GMCs and star clusters. GMCs are initially resolved by $128^{3}$ equal-mass gas cells, which sets the target mass $m_{\rm res}\approx M_{\mathrm{GMC}}/128^{3}$. We perform a convergence test for the RHO20T run with $f_{\mathrm{boost}}=2$ by varying the number of resolution elements from $64^{3}$ to $256^{3}$ and find that the star formation histories are not sensitive to mass resolutions. The final star formation efficiencies in runs with different resolutions vary by only a few percent.

2.5 Caveats of the sub-grid models

The sub-grid model used in this paper samples star particles probabilistically with star formation criteria described in Section 2.2. This is different from the more realistic sink particle approach that follows the accretion history of individual stars. We adopt this simplified star formation algorithm since the goal of this paper is to investigate the global properties of GMCs and star clusters but not to study the origin of the IMF or the detailed formation of single stellar objects. We interpret each star particle as a single stellar population, whose feedback intensity is estimated in an IMF-averaged fashion, see Section 2.3. For the most massive GMCs in our simulations, this IMF-averaged approach captures the overall energy budget of stellar feedback (see also Grudić et al., 2018b). However, it unavoidably underestimates the large variation of star formation efficiency in low-mass clouds, where a few massive stars can dominate the feedback process. In addition, the star particles in our simulations with a fixed gravitational softening length only trace the overall mass distribution of star clusters. Therefore, the detailed dynamical evolution could in principle depend on the choice of softening length (e.g. Bate et al., 2003; Bate, 2012). To fully capture the collisional process between star particles requires simulations that resolve the formation of individual stars over the whole mass spectrum and a more accurate gravity integrator, such as NBODY6 (Aarseth, 1999). Recent efforts have been made towards this direction (e.g. Parker et al., 2015; Gavagnin et al., 2017; Wall et al., 2019), however, due to the high computational cost, these simulations mainly focused on less massive GMCs and cannot explore the dynamical evolution of GMCs over a large parameter space.

Another main caveat is that we only take into account the momentum feedback from stellar winds, whose intensity is controlled by $f_{\rm boost}$. A larger $f_{\rm boost}$ is used to mimic stellar feedback from multiple feedback sources, such as stellar winds, ionizing radiation, and supernovae. However, in reality, different feedback mechanisms operate on different time-scales, at different locations, and through different physical processes. For example, ionizing radiation from massive stars strongly alters the ionizing state of the gas around massive stars and deposits both internal and kinetic energy to the surrounding medium. Moreover, for some of our simulations with very long $t_{\rm ff}$, for example SIGMA40, at the late stage of the cloud evolution when the age of the massive stars is longer than their main-sequence lifetime, core-collapse supernovae can deposit an enormous energy to the cloud and violently disrupt it. The combined effects of various feedback mechanisms will be investigated in an upcoming paper.

3 Integrated star formation efficiency (SFE)

In total, we have performed 80 GMC simulations with different masses, radii, velocity configurations, and feedback boosting factors. For all runs, we stop the hydrodynamical simulations when 99% of the gas mass is expelled from the initial spherical regions by momentum feedback. Although different GMCs show quantitative different star formation histories and final efficiencies, the general evolutionary stages of the clouds are very similar. Here we use the RHO20T run with $f_{\mathrm{boost}}=2$ as an example to describe the general pattern of GMC evolution.

Figure 1 shows the gas density projection of this run at four different epochs. The cloud evolution is initially governed by the turbulent velocity field which creates complex filamentary structures (upper left). After $\sim 0.3\tau_{\mathrm{ff}}$, a roughly log-normal density distribution is established due to the supersonic turbulence. As turbulent energy dissipates by supersonic compression, the cloud starts to experience global contraction under self-gravity. Many subclusters are formed at the intersection of the filaments where dense gas clumps experience local runaway collapse (upper right). These subclusters move along the filaments, merge with each other frequently, and eventually form more massive subclusters. Due to momentum feedback, some gas mass is channelled outwards through low-density regions. In contrast, the high-density regions are compressed further and form young stars subsequently (lower left). When the central star cluster is massive enough so that its momentum feedback is able to counteract gravitational contraction, the majority of the gas mass is expelled from the cloud centre, causing the formation of giant wind-blowing bubbles (lower right).

3.1 Time evolution of cloud properties

In Figure 2, we quantify the time evolution of various physical quantities of the cloud for the same run shown in Figure 1. Panel (a) shows the star formation history of the cloud until it is fully disrupted. After the first group of stars forms at $\sim 0.3\tau_{\mathrm{ff}}$, the SFR rises dramatically and peaks at around $t_{50}\sim\tau_{\mathrm{ff}}$. As momentum feedback from stellar particles clears some gas mass from the cloud centre, the SFR drops gradually. Although the whole star formation activity spans over $\sim 2\tau_{\mathrm{ff}}$, the majority of the stellar mass is formed around the epoch of the star formation peak. The central 80% of the stellar mass is assembled within 0.6-1.4 $\tau_{\mathrm{ff}}$.

The cumulative version of panel (a), the stellar mass growth history, is shown in panel (b), together with the evolution of gas mass. We split all gas mass into bound and outflow components. The bound gas is defined as the total mass of gas cells with negative (kinetic+potential) energy, while the outflow is defined as the unbound gas cells that are outside twice the initial GMC radius. It is clear that the decrease of bound gas mass is partly due to gas consumption by star formation and partly due to the increase of feedback-driven outflows. The final stellar mass reaches $\sim 60$% of the initial cloud mass, which is defined as the integrated SFE, $\epsilon_{\mathrm{int}}$.

Panel (c) shows the evolution of virial parameters for both gas ($\alpha_{\rm gas}$) and stars ($\alpha_{\rm star}$). By construction, gas is initially in virial equilibrium with $\alpha_{0}=1$. As the turbulent energy dissipates, the kinetic energy of the cloud decreases and the system collapses, which leads to a slight decrease of $\alpha_{\rm gas}$. The momentum feedback from stars adds kinetic energy to the gas cells and helps to increase $\alpha_{\rm gas}$ after $t\sim 1.3\tau_{\mathrm{ff}}$. Eventually, the virial parameter comes back to unity. Yet momentum feedback cannot keep the cloud in a quasi-equilibrium state. $\alpha_{\rm gas}$ keeps rising and becomes much larger than unity very quickly until the majority of the gas mass is removed from the central region of the cloud. Interestingly, the virial parameter of stars, $\alpha_{\rm star}$, is always smaller than unity, suggesting that the model star cluster is sub-virial. As we will show later, this sub-virial dynamical state before gas expulsion has a dramatic effect on the formation of bound clusters in GMCs.

In panel (d), we present the evolution of the half-mass radius of the gas and stellar components of the GMC. The evolution of the half-mass radius of the gas tightly follows the evolution of its dynamical state as is described in panel (c). The size of the gas cloud first shrinks slightly due to the dissipation of the initial turbulent energy until stellar momentum feedback puffs it up. The evolution of the half-mass radius of the star cluster is more complicated. At first, stars are formed in dense gas clumps distributed over a large volume of the cloud, which leads to a relatively large initial stellar radius. As the cloud contracts, star formation activities concentrate more towards the central region and the stellar half-mass radius decreases until $t_{50}$. Later, as gas removal shallows the overall gravitational potential, the star cluster expands dynamically to reach a new equilibrium state.

Many theoretical works on gas expulsion and the formation of bound fraction suggest that, rather than $\epsilon_{\mathrm{int}}$, the local stellar fraction (LSF) is considered as an effective SFE to better probe the bound fraction of the cluster after gas expulsion (Goodwin, 2009; Smith et al., 2011, 2013; Farias et al., 2018). The LSF is defined as the mass fraction of stellar mass within the stellar half-mass radius right before the gas expulsion:

[TABLE]

In panel (e), we show the evolution of the LSF in the simulations. Since neither the formation of stars nor the dispersion of gas happens instantaneously, the LSF changes dramatically during the course of GMC evolution. By definition, the LSF is initially zero. As star formation continues, the gas is gradually consumed and expelled from the central region, causing the increase of stellar mass and decrease of gas mass. Therefore, the LSF increases monotonically as a function of time until it reaches unity.

3.2 Effects of $f_{\mathrm{boost}}$ on the star formation history

Figure 3 shows the star formation histories of the RHO20T run with different $f_{\mathrm{boost}}$. During the early stages, the total mass of young stars is so small that momentum feedback is not enough to affect the dynamical state of the cloud. Therefore, the effects of $f_{\mathrm{boost}}$ on the SFR is not visible and all lines overlap with each other until the SFR reaches its peak. During this period of time, the SFR presents a linear increase with time, ${\rm SFR}\propto t$. This linear time dependence is consistent with previous theoretical and numerical prediction of turbulent self-gravitating cloud with virial parameters close to unity (Lee et al., 2015; Murray & Chang, 2015; Murray et al., 2017).

When a sufficient fraction of gas mass is converted to stars, the momentum feedback is able to alter the overall dynamical state of the cloud (see Section 3.1), and eventually reverts the increasing trend of the SFR. The exact epoch of the turning point and, in turn, the final $\epsilon_{\mathrm{int}}$, is determined by the balance between the feedback intensity and gravitational contraction, which will be discussed in the next section.

3.3 Surface density-dependent SFE

For all 80 runs, we obtain the total stellar mass at the end of the hydrodynamical simulations and calculate the integrated SFE, $\epsilon_{\mathrm{int}}$. Figure 4 shows $\epsilon_{\mathrm{int}}$ as a function of the cloud initial surface density, $\Sigma_{0}$, for all 80 runs. We find a positive correlation between $\epsilon_{\mathrm{int}}$ and $\Sigma_{0}$, for a given value of $f_{\mathrm{boost}}$. In contrast, we do not find clear correlations of $\epsilon_{\mathrm{int}}$ with either the initial mass, radius, or volume density of the clouds. For the same GMC, runs with larger $f_{\mathrm{boost}}$ produce less stars and smaller $\epsilon_{\mathrm{int}}$, consistent with the results described in Section 3.2. Moreover, we find that rotation-supported (“R”) runs in general show a slightly higher $\epsilon_{\mathrm{int}}$ than the corresponding turbulence-supported (“T”) runs, especially when a large $f_{\mathrm{boost}}$ is employed. This can be explained as follows. Because of the initial rotational velocity, GMCs in the “R” runs first collapse to a disk-like structure whose scale height is typically smaller than the cloud radius. The formation of the thin disk allows momentum feedback to escape easier than that of the “T” runs, where the spherical shape of cloud is roughly maintained. This geometric effect leads to a difference of $\epsilon_{\mathrm{int}}$ by about 10-20%.

We next build an analytical model to explain the correlation between $\epsilon_{\mathrm{int}}$ and $\Sigma_{0}$ by considering the force balance between gravitational contraction and gas expulsion by momentum feedback. We assume, when the balance is achieved, the residual gas forms a thin spherical shell with a radius $r_{s}$. For gravitational forces, we consider the contribution from both the central star cluster of mass $M_{*}$ and self-gravity of the gas shell of mass $M_{\rm sh}=M_{\mathrm{GMC}}-M_{*}$. The gravitational force per unit area of the shell from the cluster is evaluated as

[TABLE]

while self-gravity of the gas shell is

[TABLE]

where $A_{s}=4\pi r_{s}^{2}$ is the surface area of the shell, $\Sigma_{\rm sh}=M_{\rm sh}/4\pi r_{s}^{2}$, and $\beta$ is the geometric factor that takes into account the anisotropic distribution of the gas shell. Note that $\Sigma_{\rm sh}$ is the surface density of the spherical shell seen from the central cluster, different from $\Sigma_{0}$, the cloud column density, by a factor of 4: $\Sigma_{0}=4\Sigma_{\rm sh}$. For a uniform density gas distribution, $\beta=0.5$. The expel force per unit area exerted onto the gas shell by momentum feedback is

[TABLE]

where $\dot{p}$ is the momentum deposition rate per unit stellar mass. As described in Section 2.3, we use an IMF-averaged wind injection as the default setup for feedback with a boosting factor, $f_{\mathrm{boost}}$. In Equation 2, the deposition rate evolves as the stellar population ages. For simplicity, in this analytical model we assume a constant momentum deposition rate per unit mass $\dot{p}=f_{\mathrm{boost}}\dot{p}_{\rm w}$ where $\dot{p}_{\rm w}$ is the IMF- and time-averaged value.

Force balancing between gravitational collapse and the momentum-driven wind, $F_{\rm sh,*}+F_{\rm sh}=F_{p}$, gives

[TABLE]

By defining variables $\Sigma_{\rm crit}\equiv\dot{p}/\pi G$ and $\Gamma\equiv\Sigma_{\rm sh}/\Sigma_{\rm crit}=\pi G\Sigma_{0}/4f_{\mathrm{boost}}\dot{p}_{\rm w}$ and assuming $r_{s}=R_{\mathrm{GMC}}$, the above equation can be simplified to

[TABLE]

Finally, $\epsilon_{\mathrm{int}}$ can be solved as

[TABLE]

For a uniform density gas shell, Equation 11 reduces to

[TABLE]

For clouds with high surface density, $\Gamma\gg 1$, Equation 11 is reduced to $\epsilon_{\mathrm{int}}\approx 1-1/\Gamma$, which suggests that almost all gas mass is converted into stars before the gravitational collapse is balanced by momentum feedback. Since the mass of the gas shell is much smaller than the mass of the central star cluster, self-gravity of the gas shell is negligible and therefore $\epsilon_{\mathrm{int}}$ is independent of $\beta$. For $\Gamma\ll 1$, on the other hand, the gravitational force is dominated by the self-gravity of the shell and Equation 11 can be simplified to $\epsilon_{\mathrm{int}}\approx\beta\Gamma$, which shows a clear $\beta$ dependence. In this case, $\epsilon_{\mathrm{int}}$ depends linearly on the cloud surface density divided by the momentum deposition rate, $\Sigma_{0}/f_{\mathrm{boost}}\dot{p}_{\rm w}$.

We fit the value of $\epsilon_{\mathrm{int}}$ for all 80 GMC simulations using Equation 11, and obtain the best-fit parameters with $1\sigma$ uncertainty: $\beta=1.83\pm 0.89$ and $\dot{p}_{\rm w}=(3.32\pm 0.64)\times 10^{-9}{\rm cm/s^{2}}$. As can be seen in Figure 4, the analytical model is in good agreement with the simulated $\epsilon_{\mathrm{int}}$ over a large range of $\Sigma_{0}$ and $f_{\mathrm{boost}}$. We notice that the analytical model overestimates $\epsilon_{\mathrm{int}}$ for runs with $f_{\mathrm{boost}}=10$. This is possibly because clouds are disrupted earlier in $f_{\mathrm{boost}}=10$ runs than other runs and the time-averaged $\dot{p}$ is systematically higher due to the decreasing wind kinetic luminosity used in the simulations, see Equation 2.

3.4 Star formation time-scales

As described in Section 3.2 and 3.3, stronger momentum feedback changes the epoch of the peak of star formation to earlier times and reduces the final SFE. How important is the strength of feedback to the overall star formation time-scales? Can feedback be the main energy source to support the cloud and maintain a quasi-equilibrium state? We investigate these questions here by defining several relevant time-scales that characterize the star formation activities for the simulated GMCs. First, we define the initial free-fall time of the cloud as

[TABLE]

where $\rho_{0}=3M_{\mathrm{GMC}}/4\pi R_{\mathrm{GMC}}^{3}$ is the initial volume density of the GMC. The free-fall times of all GMCs are listed in Table 1. We define the star formation duration as the time-scale during which the clouds form the central 80% of their stars: $\tau_{\mathrm{dur}}=t_{90}-t_{10}$. Following Li et al. (2018); Grudić et al. (2018b), we also define an age spread of star cluster as the ratio between the final stellar mass $M_{*}$ and the mass-weighted SFR:

[TABLE]

where $\dot{M}_{*}$ is the instantaneous SFR. For a Gaussian-like star formation history with standard deviation $\sigma_{*}$, the age spread is approximated $\tau_{\rm spread}\approx 2\sqrt{\pi}\sigma_{*}$.

In Figure 5, we show the central star formation epoch ($t_{50}$), star formation duration ($\tau_{\rm dur}$), and age spread ($\tau_{\rm spread}$) as a function of initial free-fall time ($\tau_{\mathrm{ff}}$) for all 80 GMCs. In the top panel, we find a clear linear correlation between $t_{50}$ and $\tau_{\mathrm{ff}}$. We know that $t_{50}$ roughly represents the epoch of the peak of star formation because of the Gaussian-like shape of the star formation history, see Section 3.1. In fact, $t_{50}$ is close to $\tau_{\mathrm{ff}}$ for clouds that are turbulence-supported (“T” runs). For “T” runs, the initial turbulent energy dissipates within the turbulence crossing time, which is shorter than the free-fall time of the cloud. The peak of star formation is determined by gravitational collapse of the whole cloud and is therefore similar to the free-fall time. The “R” runs, on the other hand, show a systematically larger $t_{50}$ than the corresponding “T” runs. This is because the rotation-supported cloud first collapses along the rotational axis and form a gaseous disk. The rotating disk contains more coherent motions whose kinetic energy dissipates over a longer time-scale than turbulent motions. Interestingly, as shown in the middle and bottom panels, the star formation duration and age spread for “T” and “R” runs do not show clear difference, which suggests that once the runaway collapse starts, the details of the initial configuration of the gas motion do not affect the subsequent star formation process.

Similar to $t_{50}$, $\tau_{\mathrm{dur}}$ and $\tau_{\mathrm{spread}}$ also correlate linearly with $\tau_{\mathrm{ff}}$. We perform a linear fit to the correlation between $\tau_{\mathrm{ff}}$ and $\tau_{\mathrm{dur}}$ and find that runs with different $f_{\mathrm{boost}}$ show similar scalings but with different normalizations. Although the normalization of the relations shows an anticorrelation to $f_{\mathrm{boost}}$, increasing $f_{\mathrm{boost}}$ by a factor of 20 from 0.5 to 10 only shortens the time-scale by a factor of 3 to 4. Quantitatively, this weak dependence of the momentum feedback intensity on the star formation duration can be understood as follows. In Section 3.2, we find a linearly increasing SFR from $t_{10}$ to $\sim t_{50}$ regardless the choice of $f_{\mathrm{boost}}$. Here we define ${\rm SFR}=A(t-t_{10})$, where $A$ is an arbitrary normalization. Therefore, the approximated total stellar mass is $M_{*,\rm{linear}}=5/2\int_{t_{10}}^{t_{50}}A(t-t_{10})dt=5/4A(t_{50}-t_{10})^{2}\sim 5/16A\tau_{\mathrm{dur}}^{2}\propto\tau_{\mathrm{dur}}^{2}$, assuming $\tau_{\rm dur}=2(t_{50}-t_{10})$. We also obtained a correlation between $\epsilon_{\mathrm{int}}$ and $\Sigma_{0}$ from Equation 11, which gives the final stellar mass as $M_{*}=\epsilon_{\mathrm{int}}M_{0}$. For clouds with $\Sigma_{0}<<\Sigma_{\rm crit}$, the above expression can be simplified as $M_{*}\approx\beta\Gamma M_{\mathrm{GMC}}=\pi G\beta\Sigma_{0}M_{\mathrm{GMC}}/8f_{\mathrm{boost}}\dot{p}_{\rm wind}\propto f_{\mathrm{boost}}^{-1}$. Equating $M_{*,\rm{linear}}$ and $M_{*}$ gives a scaling $\tau_{\mathrm{dur}}\propto f_{\mathrm{boost}}^{-1/2}$. For high surface density clouds when $\Sigma_{0}>>\Sigma_{\rm crit}$, $\epsilon_{\mathrm{int}}\approx 1-1/\Gamma$ and therefore $\tau_{\rm dur}$ is almost independent of $f_{\mathrm{boost}}$. Indeed, we find that the correlation between $\tau_{\rm dur}$ and $f_{\mathrm{boost}}$ in the simulations scales between $\tau_{\mathrm{dur}}\propto f_{\mathrm{boost}}^{-1/2}$ and $\tau_{\mathrm{dur}}\propto f_{\mathrm{boost}}^{0}$, with a median power-law slope around -1/4. The weak dependence of star formation duration to the strength of momentum feedback suggests that the cluster formation time-scale is mainly determined by gravitational runaway collapse. Keep in mind that the turbulent velocity fields used in our simulations are initialized at the very beginning of the simulations. No subsequent turbulence driving is applied to feed in the kinetic energy after turbulence dissipation. Understanding how turbulent motions cascade from large-scale environments to the local star-forming regions and affect the long-term star formation activities requires simulations of GMCs in realistic galactic environments, which will be investigated in a future work.

4 Bound fraction of model clusters

As discussed in previous sections, momentum feedback from young stars disrupts star-forming regions, reduces the SFE, and shortens star formation time-scales. The gas expulsion and cloud disruption flatten the gravitational potential and inevitably leave some imprint on the dynamical state of the star clusters formed at the centre of the clouds. Previous works have studied extensively the effects of gas removal on the dynamical evolution of star clusters. A simple virial analysis shows that if more than 50% of the mass is instantaneously removed from a virialized system, the remaining mass will become gravitationally unbound (Hills, 1980; Mathieu, 1983). This conclusion is based on several assumptions: (1) the system is initially in virial equilibrium; (2) the removed mass is initially well mixed with the residual mass; (3) the mass-loss time-scale is much shorter than the dynamical time-scale of the system. In realistic star-forming environments, all of the above assumptions are not strictly applicable. The star-forming regions are not necessarily in virial equilibrium (e.g. Offner et al., 2009). Stars are not randomly distributed within the GMCs, but are formed hierarchically within the densest molecular cores at the intersections of the gas filaments (e.g. Smith et al., 2011, 2013; Farias et al., 2015; Lee & Goodwin, 2016; Farias et al., 2018). The non-star-forming gas is expelled outward gradually rather than instantaneously (Geyer & Burkert, 2001; Smith et al., 2013) and is preferentially channelled through low-density holes and tunnels rather than being removed homogeneously. Here, we explore these factors in our simulations and investigate how gas expulsion affects the bound fraction of star clusters.

4.1 Stellar and gas distribution at $t_{50}$

Figure 6 shows the density profiles for both gas and stars at $t_{50}$ for all 80 runs. The profiles are centred at the location of the deepest gravitational potential, which is usually the centre of the central star clusters. Both the gas and stellar profiles are normalized to the central stellar density of the corresponding run. We find that, for most of the runs, the stellar central density is systematically higher than that of the gas, suggesting that a large fraction of gas mass in the central region of the GMCs has already been converted to stars at $t_{50}$. As a result, the gravitational potential is dominated by stars rather than gas within $\sim 0.1R_{\mathrm{GMC}}$. We fit both the gas and stellar density profiles with a power-law shape, $\rho(r)\propto r^{-\gamma}$, and show the distribution of the power-law slopes $\gamma$ in the inset of the figure. The gas density profiles shows a roughly isothermal profile with a median power-law slope $\gamma\sim 1.9$. It is interesting that the gas distribution for all runs converges to quasi-isothermal profiles regardless of the fact that the initial conditions are uniform density spheres. This isothermal gas density profile is consistent with the observed radial profiles of star-forming molecular clumps (Mueller et al., 2002; Palau et al., 2014; Wyrowski et al., 2016; Csengeri et al., 2017) and is thought to be a natural consequence of scale-free gravitational collapse (e.g. Larson, 1969; Penston, 1969; Naranjo-Romero et al., 2015; Donkov & Stefanov, 2018; Li, 2018). The stellar density profiles, on the other hand, are systematically steeper than that of the gas with slopes centred around 2.8 with a large variation. The steeper stellar density profiles imply more centrally concentrated star formation activities.

4.2 Virial state of star clusters at the peak of star formation

In panel (c) of Figure 2, we showed the evolution of the virial parameter of stars in the RHO20T run with $f_{\mathrm{boost}}=2$. We found that the model star cluster is in a sub-virial state ($\alpha_{\rm vir,*}$<1) during the course of gas expulsion from $t_{50}$ to $t_{\rm exp}$. Here we calculate $\alpha_{\rm vir,*}$ at $t_{50}$ for all 80 GMCs in order to quantify the dynamical state of the star clusters before gas expulsion. We find that $\alpha_{\rm vir,*}$ has a median value around 0.61 with a 25-75% interquartile range 0.55-0.65, which suggests a systematic sub-virial dynamical state. This finding is qualitatively consistent with previous simulations, such as Offner et al. (2009), who suggest a sub-virial stellar velocity dispersion even in virialized GMCs. This is possibly because stars are preferentially formed in the densest molecular cores which on average have less velocity dispersion than the rest of the cloud. Moreover, as discussed in Section 4.1, the stellar distribution is much more compact than that of the gas, which also helps the central cluster to remain gravitationally bound. As will be shown later, the sub-virial state has a dramatic effect on the final boundness of star clusters after gas expulsion.

4.3 Gas expulsion time-scales vs dynamical time-scales

The dynamical response of star clusters to gas expulsion is a competition between the flattening of the gas potential that happened over gas expulsion time-scale and the energy exchange among stars that happened over the dynamical time-scale of the clusters. Previous $N$-body simulations mimic the gas expulsion process by gradually reducing the background gas potential over a given period of time (e.g. Geyer & Burkert, 2001; Baumgardt & Kroupa, 2007). However, in realistic star-forming regions, star formation and gas expulsion happen at the same time and there is no clear separation between the two processes.

Here we define the gas expulsion time-scale as the duration from the peak of star formation to the epoch when the contribution of the gravitational potential energy from gas mass within twice the half-mass radius of the star cluster is less than 10%, $\tau_{\rm exp}=t_{\rm exp}-t_{50}$. We also modified the potential energy threshold from 1 to 10% and find the expulsion time-scale is not sensitive to the choice of this value. Similar to the star formation duration, we find that $t_{\rm exp}$ depends strongly on the initial free-fall time of the GMC, suggesting that gas expulsion associates well with the end of star formation. It also suggests that gas expulsion happens neither instantaneously nor much longer than the dynamical time-scale of the clouds.

4.4 Calculating the bound fraction of star clusters

The final efficiency of star formation and all relevant time-scales investigated above are determined once the majority of the gas mass is removed from the central region of the clouds. Right after gas removal, the dynamical state of star clusters will readjust according to the changes of gravitational potential for the next couple of dynamical time-scales. Therefore, after momentum feedback expels more than 99% of the gas mass out of twice the stellar half-mass radius, we stop the hydro runs, remove residual gas cells, and continue evolving the star clusters in a gravity-only mode with the same softening length of stellar particles as the corresponding hydro runs. The simulations keep running for another two free-fall times of the GMC, $\tau_{\mathrm{ff}}$. We analyze the cluster bound fraction from $N$-body snapshots at different epochs and find that the bound fraction usually becomes stable after only $\sim 0.5\tau_{\mathrm{ff}}$.

We adopt two methods to estimate the bound fraction from the last snapshot of the $N$-body runs. The first and simplest way is to calculate the mass fraction of all stellar particles with negative energies. This method gives accurate results for clusters that have large bound fractions. However, for clusters with low bound fraction, simply summing up stellar particles with negative energy overestimates the bound fraction. Stars with negative energy are not guaranteed to be bound to the cluster since removing all stars with positive energy shallows the gravitational potential. Therefore, we design a new method that removes stellar particles with positive energies and updates gravitational energies for the remaining stars iteratively. The iteration stops when all remaining stars have negative energies and the bound mass fraction of the clusters, $f_{\mathrm{bound}}$, is obtained.

The second method is to use the “Lagrangian radius”, $R_{\rm lag}$, defined as a series of radii within which the star cluster contains a sequence of fractions of stellar mass. Brinkmann et al. (2017) suggest to use the evolution of the Lagrangian radii to determine the structural changes of the star cluster during gas expulsion. The bound fraction after gas expulsion is determined by the outermost Lagrange radius that shows a core collapse. Figure 7 shows the time evolution of the Lagrangian radii for the RHO20T run with $f_{\mathrm{boost}}=2$. We find that the Lagrangian radii of all mass fractions decrease during the first couple of $\tau_{\mathrm{ff}}$ due to the same reason of the decrease of the stellar half-mass radius as mentioned in Section 3.1. After the majority of stars are formed after $\tau_{\mathrm{ff}}$, $R_{\rm lag}$ starts to increase in response to gas removal and the decrease of the total gravitational potential. The $R_{\rm lag}$ of small mass fractions shows a turnover when the central component recollapses to form the central bound clusters. The evolution of $R_{\rm lag}$ for mass fraction the same as the bound fraction determined by the iterative method is highlighted in the same figure. We find that the evolution of this $R_{\rm lag}$ is indeed approximately the outermost Lagrange radius that shows a turnover. We have performed the same analysis for all 80 runs and find that the iterative and “Lagrangian radius” methods give consistent results, confirming the eligibility of both methods. Because the outermost Lagrangian radius with turnover is determined somewhat subjectively while the iterative method always converges to an accurate result, we will only report the bound fraction that is determined by iterative methods in later sections.

4.5 Bound fraction as a function of $\epsilon_{\mathrm{int}}$

Figure 8 shows a compilation of bound fractions for all 80 GMCs as a function of $\epsilon_{\mathrm{int}}$. The bound fractions are calculated from the last output of the $N$-body runs using the iterative method described in Section 4.4. We find that there exists a broad range of $f_{\mathrm{bound}}$ from almost completely bound to almost completely disruptive. There is an increasing trend of $f_{\mathrm{bound}}$ as a function of $\epsilon_{\mathrm{int}}$. Interestingly, several runs with $\epsilon_{\mathrm{int}}<0.5$ show large $f_{\mathrm{bound}}$, deviating from the simple virial analysis. For some runs with $\epsilon_{\mathrm{int}}\sim 0.2$, they still form star clusters with large $f_{\mathrm{bound}}>0.5$. The emergence of bound clusters in low-SFE clouds is actually in line with the findings in previous hydrodynamic simulations of GMCs (e.g. Parker et al., 2015; Gavagnin et al., 2017; Farias et al., 2018), although the detailed relationship between $\epsilon_{\mathrm{int}}$ and $f_{\mathrm{bound}}$ shows subtle differences.

Here, we present a simple one-zone cluster model to explain the relationship between $\epsilon_{\mathrm{int}}$ and $f_{\mathrm{bound}}$. Assuming stars in clusters always follow a Maxwellian velocity distribution

[TABLE]

where $x\equiv\sqrt{3}v/v_{\mathrm{rms}}$, $v_{\mathrm{rms}}=\sqrt{3kT/m}$, and $m$ and $T$ are the average mass and “temperature” of the cluster. Before gas expulsion, a cluster with mass $M_{*}$ and radius $r_{*}$ has a virial parameter $\alpha_{*}=-2T_{*,0}/\Omega_{*,0}$, where the kinetic energy of stars $T_{*,0}=M_{*}v_{\mathrm{rms}}^{2}/2$ and $\Omega_{*,0}=-GM_{\mathrm{GMC}}M_{*}/r_{*}$. Instantaneous gas expulsion flattens the gravitational potential but does not change the kinetic energy of the stars. Therefore, the potential energy of stars after gas expulsion drops according to the SFE, $\Omega_{*,1}=-GM_{*}M_{*}/r_{*}=\epsilon_{\mathrm{int}}\Omega_{*,0}$, while the kinetic energy does not change, $T_{*,1}=T_{*,0}=M_{*}v_{\mathrm{rms}}^{2}/2$. The escape velocity of the cluster after gas expulsion is

[TABLE]

Assuming stars keep their Maxwellian distribution after gas expulsion, the bound fraction $f_{\mathrm{bound}}$ can be estimated as the fraction of stars that have velocities below $v_{\rm esc}$:

[TABLE]

where $f_{\rm sat}$ is the saturation of bound fractions. This saturation is probably due to some small fraction of substructures that are formed with low local star formation efficiencies and do not merge into the central cluster.

We fit the simulated $f_{\mathrm{bound}}$ using the above relationship with two parameters, $\alpha_{*}$ and $f_{\rm sat}$. The best-fit values and their corresponding 1-$\sigma$ confidence intervals are $\alpha_{*}=0.48\pm 0.02$ and $f_{\rm sat}=0.94\pm 0.03$. The best-fit $\alpha_{*}$ is consistent with the measured virial parameters at $t_{50}$ (see Section 4.2), suggesting that the sub-virial dynamical state of model clusters is the main driver to prevent clusters from being dissociated by gas expulsion.

We compare our best-fit relation with the semi-analytical model developed by Adams (2000), who evaluates the dynamical response of star clusters to instantaneous gas expulsion based on an equilibrium cluster model with different stellar and gas density profiles. The analytical fit to his isotropic velocity distribution model gives $f_{\mathrm{bound}}=(2\tilde{\epsilon}-\tilde{\epsilon}^{2})^{2/3}$, where $\tilde{\epsilon}\equiv(10\epsilon_{\mathrm{int}}-1)/9$. We find that this model is roughly consistent with our simulation result but slightly underestimates $f_{\mathrm{bound}}$ for $\epsilon_{\mathrm{int}}<0.5$. Moreover, we also show the $\epsilon_{\mathrm{int}}$-$f_{\mathrm{bound}}$ relation used in our previous cosmological hydrodynamic simulations of a Milky Way-sized galaxy. In these simulations, a continuous cluster formation prescription is used to model the formation of individual star clusters from dense clumps resolved by parsec-scale resolutions (Li et al., 2017). The adopted $\epsilon_{\mathrm{int}}$-$f_{\mathrm{bound}}$ relation, $f_{\mathrm{bound}}={\rm min}(2\epsilon_{\mathrm{int}},1)$, determines the final bound mass of model clusters and in turn affects the cluster initial mass function, cluster formation efficiency, and the properties of evolved cluster populations (Li et al., 2018; Li & Gnedin, 2018). The piecewise function used in those cosmological simulations is in general agreement with the relationship found here.

Note that the bound fractions obtained above take into account the total bound stellar mass from not only the central star cluster but also all other surrounding substructures, which are not necessarily bound to the central cluster. In the following section, we will identify individual subclusters in the simulations, estimate the bound mass of central star clusters, and quantify its relationship to the total bound mass.

4.6 Properties of substructures

Previous studies on the bound fraction of star clusters after gas expulsion usually assume an initial spherical gas and stellar distribution. However, recent observational and theoretical works suggest a hierarchical star formation scenario due to the turbulent nature of GMCs (e.g. Elmegreen & Elmegreen, 2001; Bonnell et al., 2003; Allen et al., 2007; Bate, 2009; Gutermuth et al., 2009; Bressert et al., 2010; Girichidis et al., 2011; Maury et al., 2011). In Figure 9, we show the stellar particle distributions in the $x$-$y$ plane for the RHO20T run with $f_{\mathrm{boost}}=2$ at four different epochs. At $t=\tau_{\mathrm{ff}}$, the stellar distribution follows well with the gas distribution and subclusters are distributed along the filamentary structures. At $t=2\tau_{\mathrm{ff}}$ when the majority of the gas mass has already been pushed out of the central region, we find that some subclusters spiral into the centre of the GMCs due to gravitational attraction and, at the same time, merge with each other, and form more massive subclusters. Some of the most massive subclusters show a non-spherical shape because of recent mergers. At $t=$ 3-4 $\tau_{\mathrm{ff}}$, dynamical evolution after gas expulsion and violent relaxation during mergers erase the memory of the turbulent configurations of the gas cloud and help circularize the central star clusters. Some small substructures with high bulk velocities escape from the central region and never return to the central cluster.

To fully analyze the behavior of subclusters and quantify the central cluster properties, we identify bound substructures in the $N$-body simulations by adopting the SUBFIND algorithm (Springel et al., 2001). The star particles are first linked into friends-of-friends (FOF) groups with separation less than 0.17 of the mean particle separation. For each FOF group, the SUBFIND algorithm is applied to identify all bound subclusters. We report bound substructures that are resolved by at least 200 stellar particles and assign the most massive bound subcluster as the central cluster. In most cases, the central cluster sits very close to the centre of the GMC and is surrounded by other less massive subclusters.

We find that the number of subclusters in the “R” runs is systematically larger than that in the corresponding “T” runs, although the total bound stellar masses are similar. In Figure 10, we examine the relationship between the bound fraction, $f_{\rm bound}$, and the ratio of the mass of the central cluster to the total stellar mass, $f_{\rm central}$. In most cases, the central clusters dominate the total bound mass of the system. We also find that $f_{\rm central}$ is systematically higher for “T” runs than for “R” runs. Statistically, more than 50% of the bound stellar mass is contributed by the central clusters for “T” runs. This fraction is smaller for “R” runs but it exhibits a large scatter. As we discussed in the previous sections, GMCs in the “R” runs first collapse along the z-axis and form a gas disk. Therefore, many subclusters, that are formed in the dense clumps within the disk, obtain a similar bulk rotational velocity and are less likely to merge with each other than those formed in the “T” runs. The large number of subclusters in the “R” runs also implies that the stellar mass is distributed more broadly across different subclusters. Therefore, the central clusters in the “R” runs are less massive than those in the corresponding “T” runs. The exact hierarchical structure of the subclusters in the “R” runs depends strongly on the initial setup of the turbulent velocity fields, making the prediction of the mass of central star clusters less promising. This difficulty is reflected in the very large scatter of $f_{\rm central}$ for clouds with $f_{\mathrm{bound}}>0.8$.

It should be noted that recent observations of cluster formation efficiency, defined as the fraction of stars formed in bound clusters, identify centrally-concentrated clusters of spherical shape as bound clusters (e.g. Adamo et al., 2015). This cluster selection method suggests that only the central clusters (and maybe other most massive subclusters) formed in GMCs are identified as bound star clusters when estimating cluster formation efficiencies from given galaxy patches. When interpreting the relevant observational correlations, such as the positive correlation between the SFR surface density and the cluster formation efficiency, through physics models (Kruijssen et al., 2012b; Li et al., 2018; Li & Gnedin, 2018), the effect of substructures that are not bound to the central star clusters needs to be considered and taken with caution.

5 Summary

We have performed a suite of three-dimensional hydrodynamic simulations of turbulent GMCs using the moving-mesh code Arepo with self-gravity, explicit cooling, star formation and momentum stellar feedback. We survey a large range of GMC masses and radii, and investigate the physical origin of the large variation of intrinsic SFE, $\epsilon_{\mathrm{int}}$. After the gas clouds are fully disrupted by stellar feedback, we follow the subsequent dynamical evolution of star clusters formed within the GMCs with $N$-body simulations. Below, we summarize our key conclusions:

• •

All simulated GMCs follow an initial linear increasing SFR before stellar momentum feedback disperses the whole cloud. The accelerating star formation activity leads to a superlinear stellar mass growth with time, $M_{*}\propto t^{2}$. This superlinearity is consistent with previous theoretical expectations of the gravitational runaway collapse of turbulent clouds.

• •

Momentum feedback from stellar particles adds kinetic energy to their ambient gas cells, inflates the virial parameters and radius of the clouds, drives outflows through low-density regions, and finally creates a large cavity at the centre of the clouds. The peak epoch and final efficiencies of star formation decrease with increasing strength of momentum feedback.

• •

$\epsilon_{\mathrm{int}}$ does not depend on the initial mass or radii of the clouds, but depends strongly on initial cloud surface density. This dependence is successfully explained by an analytical model that considers force balancing between gravitational collapse and momentum output from stellar particles. The model predicts $\epsilon_{\mathrm{int}}\approx 1-\Sigma_{\rm crit}/\Sigma_{0}$ for clouds with high surface density while $\epsilon_{\mathrm{int}}\propto\Sigma_{0}/f_{\mathrm{boost}}\dot{p}_{\rm w}$ for clouds with low surface density.

• •

The duration of star formation in simulated GMCs is close to the initial free-fall time of the clouds, suggesting that the cluster formation time-scale is mainly determined by gravitational runaway collapse. The duration decreases with increasing feedback intensity, although the dependence is weak: $\tau_{\rm dur}\propto f_{\mathrm{boost}}^{-1/4}$.

• •

The model star clusters are assembled hierarchically. Subclusters are formed at the many density peaks across the GMCs controlled by the initial turbulence configuration. The subclusters move along the filamentary structures and merge with each other frequently. About 50% of the mass of subclusters is merged to form the most massive central clusters, but there are always a small fraction of subclusters that are unbound to the system and fly apart from the central clusters.

• •

The gas density distribution rearranges from an initial uniform density sphere to an isothermal profile with a power-law slope $\gamma\sim 2$. At the peak of star formation, the stellar density profiles are systematically steeper than that of the gas with a power-law slope $\gamma\sim 2.8$. The steeper stellar profiles suggest a fast conversion of gas to stars and gas expulsion by stellar feedback at the centre of GMCs.

• •

The model star clusters are always in a sub-virial state with a gradually increasing virial parameter as star formation continues. Interestingly, right before gas expulsion, the virial parameters of all simulated star clusters show a consistent value around $\alpha_{\rm vir,*}\approx 0.6$.

• •

Due to the steep density profiles and sub-virial dynamical state of model clusters, clouds with low $\epsilon_{\mathrm{int}}$ (0.2-0.4) are still able to form clusters with relatively high bound fractions (0.3-0.8). The bound fraction of model clusters, $f_{\mathrm{bound}}$, is a continuously increasing function of the integrated SFE, $\epsilon_{\mathrm{int}}$. This relation is explained by a physical model that takes into account the mass fraction of stars with velocities below the escape velocity of a sub-virial system that obeys Maxwellian velocity distribution. The best-fit virial parameter of this model is around 0.5, consistent with the values obtained directly from the simulations.

Acknowledgements

We thank the anonymous referee for detailed comments and suggestions. We thank Volker Springel for giving us access to Arepo. We are grateful to Peter Behroozi, Andreas Burkert, John Forbes, Nick Gnedin, Mike Grudic, Lee Hartmann, Jens Kauffmann, and Vadim Semenov for insightful comments and suggestions. MV acknowledges support through an MIT RSC award, a Kavli Research Investment Fund, NASA ATP grant NNX17AG29G, and NSF grants AST-1814053 and AST-1814259. FM is supported by the Program “Rita Levi Montalcini” of the Italian MIUR. OG acknowledges supported through NSF grant 1412144. The simulations of this work were run on the Harvard Odyssey clusters and the Comet HPC resource at San Diego Supercomputer Center as part of XSEDE through TG-AST170042 and TG-AST180025.

Appendix A Tests of wind-blowing bubbles

We test the momentum deposition algorithm used in this paper by performing idealized simulations of wind-blowing bubbles. In this test, a stellar particle is located at the centre of a uniform density box with size 20 pc and total gas mass 320 $\,M_{\odot}$. The gas mass is initially resolved by $128^{3}$ gas cells. The central star deposits its wind material with a constant mass-loss rate $\dot{m}_{\rm w}=10^{-5}\,M_{\odot}/$yr at a fixed velocity $v_{\rm w}=500$ km/s. No self-gravity or cooling is used in this test in order to compare to analytical solutions of the evolution and internal structure of the wind-blowing bubble derived by Weaver et al. (1977).

The bubble first experiences free expansion with constant wind velocity until the mass of the swept-up material, $\frac{4}{3}\pi(v_{\rm w}t)^{3}\rho_{0}$, is comparable to the mass of the wind ejecta, $\dot{m}_{\rm w}t$. The initial free expansion phase only lasts for several hundred years and is followed by an adiabatic expansion phase. The time evolution of the bubble is $R_{\rm bubble}\approx 0.88(\dot{m}_{\rm v}v_{\rm w}^{2}/2\rho_{0})^{1/5}$, see Equation (5) in (Weaver et al., 1977). Figure 11 shows the time evolution of the wind-blowing bubble from the numerical test. The edge of the bubble is identified as the densest gas shell surrounding the stellar particles. The radius of the bubble is calculated as the mass-weighted shell radius. In addition to the time evolution of the gas shell, we also examine the internal structure of the bubble. Figure 12 shows the density, velocity, and pressure profiles of the shocked interstellar gas at $t=0.2$ Myr. All profiles are normalized to the values at the outer shock of radius $R_{2}$. The analytical solution of the self-similar adiabatic flow within the shocked medium is obtained by numerically solving Eq. (6)-(8) in Weaver et al. (1977). We confirm that the momentum deposition algorithm used in our simulations reproduces the analytical time evolution and internal structure of the wind-blowing bubble with high precision.

Appendix B Convergence to momentum deposition methods

In Section 2.3, we describe various algorithms for wind deposition to the ambient medium. Here, we test how different algorithms affect the star formation activities of the GMCs. The test runs use an identical simulation setup as the main GMC simulations described in Section 2 except wind feedback. We use three different weights to deposit momentum flux to neighbouring cells: volume, mass, and solid angle from star particles. We also test an alternative algorithm that injects stellar winds in the form of pure thermal energy, rather than momentum.

In the upper panel of Figure 13, we show the star formation histories of the GMCs using different wind deposition algorithms from the same initial condition RHO20T with $f_{\mathrm{boost}}=2$. We find that wind feedback in pure energy form is not able to disrupt the cloud and the star formation history is almost the same as that of the run without wind feedback. The failure of this energy deposition algorithm seems inconsistent with the result found in Rogers & Pittard (2013), who used thermal energy injection to simulate wind feedback and obtained significant gas outflows from the central turbulent cloud. Note that, in their simulations, self-gravity is not included and the turbulent cloud never collapses to higher density. The maximum number density reached in their simulations is about $\sim 10^{4}{\rm cm}^{-3}$. In contrast, our simulations track gravitational runaway collapse of dense clumps until they are converted to stars. In fact, our simulations always form much higher density clumps, $>10^{10}{\rm cm}^{-3}$. The thermal energies deposited into these dense environments suffer significant radiative cooling and therefore make the energy deposition inefficient. To capture the correct thermal dynamics of the adiabatic phase of the wind-blowing bubbles, the cooling radius needs to be resolved: $R_{\rm cool}\propto n^{-3/7}l_{\rm w}^{2/7}$ (Cioffi et al., 1988; Thornton et al., 1998). Since the highest density in our simulations is about six orders of magnitudes higher than that in Rogers & Pittard (2013), our simulations require about three orders of magnitude finer spatial resolution than in Rogers & Pittard (2013) to resolve the wind energy feedback, which is computationally prohibitive. We conclude that under the current simulation setup, wind feedback through thermal energy deposition is not appropriate. Alternatively, wind feedback through momentum deposition can efficiently shut off star formation activities within 2 $\tau_{\mathrm{ff}}$. We find that using volume- or solid angle-weighted schemes for momentum deposition leads to faster cloud disruption than using a mass-weighted scheme. The reason is that the majority of the momentum is deposited to the densest clumps in a mass-weighted scheme and is therefore not able to channel gas out through low-density regions (see also Smith et al., 2018). In reality, wind momentum should be deposited isotropically around star particles. It is more physically plausible to use volume and solid angle as weights.

Appendix C Convergence to numerical resolutions

Hydrodynamic simulations discretize continuous space into a finite number of resolution elements. The choice of optimal numerical resolution is to achieve a balance between scientific accuracy and computational costs. To study the convergence of the simulation outcomes to numerical resolutions, we perform GMC simulations of RHO20T runs with $f_{\mathrm{boost}}=2$ with different numbers of initial gas elements, $64^{3}$, $128^{3}$, and $256^{3}$, corresponding to target cell masses around $0.191$, $2.38\times 10^{-2}$, and $2.98\times 10^{-3}\,M_{\odot}$. The test runs are performed following the same physics as the production runs in the main text, see Section 2.

We find that the integral star formation efficiencies for the $64^{3}$, $128^{3}$, and $256^{3}$ runs are $\epsilon_{\mathrm{int}}=$ 0.612, 0.610, and 0.615, respectively. This consistency suggests that $\epsilon_{\mathrm{int}}$ is not sensitive to mass resolutions but is solely controlled by the force balance between gravitational collapse and gas expulsion by momentum stellar feedback. As shown in the lower panel of Figure 13, the star formation histories of the three runs are also in general agreement with each other. The star formation first rises dramatically and peaks at around the free-fall time of the cloud. The only noticeable difference comes from the $64^{3}$ run. This run shows a narrower star formation history than the other two runs, suggesting delayed star formation at the beginning and an earlier gas removal process after a majority of the stellar mass is formed. The $128^{3}$ and $256^{3}$ runs present a more consistent star formation history across the course of GMC evolution. Therefore, for all production runs present in the main text, we choose $128^{3}$ as the default number of resolution elements.