NUMERICAL SIMULATIONS OF TURBULENT MOLECULAR CLOUDS REGULATED BY RADIATION FEEDBACK FORCES. I. STAR FORMATION RATE AND EFFICIENCY

We run three-dimensional RHD simulations on a Cartesian grid using the Hyperion (Skinner & Ostriker 2013) extension of the Athena code (Stone et al. 2008). For this application, we solve the following simplified mixed-frame equations of RHD:

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\partial }_{t}(\rho {\boldsymbol{v}})+{\rm{\nabla }}\cdot (\rho {\boldsymbol{vv}}+P{\mathbb{I}})=-\rho {\rm{\nabla }}{\rm{\Phi }}+\rho \kappa \displaystyle \frac{{\boldsymbol{F}}}{c},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{\hat{c}}\,{\partial }_{t}{ \mathcal E }+{\rm{\nabla }}\cdot \left(\displaystyle \frac{{\boldsymbol{F}}}{c}\right)=-\rho \kappa { \mathcal E }+{\mathbb{S}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{\hat{c}}\,{\partial }_{t}\left(\displaystyle \frac{{\boldsymbol{F}}}{c}\right)+{\rm{\nabla }}\cdot {\mathbb{P}}=-\rho \kappa \displaystyle \frac{{\boldsymbol{F}}}{c},\end{eqnarray} \tag{ 4 }$$

where ρ, ${\boldsymbol{v}}$, and P are the gas density, velocity, and pressure, and Φ is the gravitational potential, all evaluated in the lab frame. We adopt the simplifying assumption of an isothermal equation of state for the gas with $P={c}_{s}^{2}\rho$ (see discussion below). The variables ${ \mathcal E }$, ${\boldsymbol{F}}$, and ${\mathbb{P}}$ are the radiation energy density, flux vector, and pressure tensor, respectively, again evaluated in the lab frame, while κ is the frequency-weighted specific material opacity calculated in the gas rest frame.

Hyperion closes the two radiative moment equations above by adopting the M₁ relation (Levermore 1984). This expresses the pressure tensor in terms of ${ \mathcal E }$ and ${\boldsymbol{F}}$, with ${\mathbb{P}}\to (1/3){ \mathcal E }{\mathbb{I}}$ in the diffusion limit ($| {\boldsymbol{F}}| /{ \mathcal E }c\ll 1$) and ${\mathbb{P}}\to { \mathcal E }\hat{{\boldsymbol{n}}}\hat{{\boldsymbol{n}}}$ in the streaming limit ($| {\boldsymbol{F}}| /{ \mathcal E }c\to 1$), where $\hat{{\boldsymbol{n}}}={\boldsymbol{F}}/| {\boldsymbol{F}}|$. We omit radiative emission terms from gas and dust in Equation (3), as we are interested in the limit in which the effects from direct stellar radiation dominate over IR emission from the dust. This stellar emission is captured in the term ${\mathbb{S}}$ in Equation (3), which describes radiative emission from star particles (see below). Finally, $\hat{c}\ne c$ is a reduced radiation propagation speed, adopted within the Reduced Speed of Light Approximation (RSLA; Gnedin & Abel 2001) to ensure that timesteps updating the radiation field are not unfeasibly short.

Hyperion divides Equations (1)–(4) into gas and radiation subsystems because these variables are transported on very different timescales. We therefore separate the two systems and solve for the gas subsystem using Athena's unsplit Van Leer (VL) integrator (Stone & Gardiner 2009), a Godunov finite-volume method adapted from the MUSCL-Hancock scheme of Falle (1991). The hydrodynamic timestep is determined using a radiation-modified CFL condition with a Courant number of 0.4 (the typical value adopted in VL integration schemes) and a radiation-modified effective sound speed that accounts for the effect of interactions between the gas and radiation fields, ${c}_{{\rm{eff}}}\equiv \sqrt{(\gamma P+4/9{ \mathcal E }(1-{e}^{-\rho {\kappa }_{0}{\rm{\Delta }}x}))/\rho }$ (Krumholz et al. 2007).

Hyperion solves the radiation subsystem using an operator-split method that separates the radiation source terms into explicit and implicit terms, with the explicit terms updated together with the update from the divergence of ${\boldsymbol{F}}$ and ${\mathbb{P}}$. Here, both the radiation energy and the flux absorption updates are solved using a standard θ-scheme update with $\theta =0.51$, very close to second order implicit in time.

The radiation subsytem is similarly solved using a VL scheme, with a Harten–Lax–van Leer (HLL) Riemann solver (González et al. 2007) used to calculate the flux between cells. In this case, the timestep is set by a CFL condition with the radiation signal speed $\hat{c}$. This reduced speed of light is chosen so that radiation timesteps are as long as possible, while still ensuring that the RSLA does not improperly affect the gas dynamics. We may achieve this so long as $\hat{c}$ is sufficiently large that the radiation field approaches equilibrium much faster than characteristic gas timescales. For streaming radiation, a practical condition is $\hat{c}\sim 10{v}_{{\rm{\max }}}\gg {v}_{{\rm{\max }}}$ (see Skinner & Ostriker 2013). There are then roughly 10 radiation substeps for each update to the gas subsystem.

In order to avoid unrealistically short timesteps for low-density regions accelerated by a strong radiative flux, there is an artificially imposed density floor. Cells whose density falls below a prescribed value at any given timestep are reset to the density floor with zero momentum. This may add mass to the grid over the course of a simulation, but in practice the density floor is chosen such that less than ∼0.1% of the initial cloud mass is artificially added over the course of a simulation.

Stars are represented within the code by point-mass sink particles (Gong & Ostriker 2013). Sink particles are formed dynamically when cells exceed a density threshold ${\rho }_{{\rm{th}}}=8.86\,{c}_{s}^{2}/(\pi G{\rm{\Delta }}{x}^{2})$ motivated by the Larson (1969) and Penston (1969) solutions for self-gravitating isothermal collapse. Locations where sink particles form must also be potential minima (Banerjee et al. 2009; Federrath et al. 2010; Vázquez-Semadeni et al. 2011). If a cell satisfies these criteria, a sink particle is created at the center of a control volume of width $3{\rm{\Delta }}x$. The sink particle has initial mass and momentum set by the sum over all control volume cells. Subsequently, gas is accreted onto the sink particles based on the HLL flux at the interface between sink control volumes and the rest of the grid. Gas variables within the control volume are set by extrapolating values from the surrounding active zones in the grid.

Sink particles are evolved in time using a leapfrog kick-drift-kick method (Springel 2005), where the particles' positions and momenta are updated alternately. The position is updated using the current velocities, while the momentum is updated based on gravitational potential differences. The potential itself is computed using particle-mesh methods with a triangular-shaped cloud kernel applied to map each particle's mass onto the grid (Hockney & Eastwood 1981). The combined particle + gas potential is found using a Fourier transform method on a domain equal to eight times the computational volume in order to implement vacuum boundary conditions for Φ (Hockney & Eastwood 1981). Finally, when the control volumes of two sink particles overlap, they are merged and placed at the center of mass of the two old particles.

Monochromatic radiation from the sink particles is emitted isotropically, representing idealized luminous stellar clusters. The source function ${\mathbb{S}}={j}_{{\rm{\ast }}}/c$ of each particle of mass ${M}_{{\rm{\ast }}}$ takes a Gaussian shape with

$$\begin{eqnarray}&&{j}_{{\rm{\ast }}}(r)=\displaystyle \frac{{L}_{{\rm{\ast }}}}{{(2\pi {\sigma }_{* }^{2})}^{3/2}}\exp \left(-\displaystyle \frac{{r}^{2}}{2{\sigma }_{{\rm{\ast }}}^{2}}\right),\end{eqnarray} \tag{ 5 }$$

with (fixed) radius ${r}_{{\rm{\ast }}}=\sqrt{2{\rm{log}}2}{\sigma }_{{\rm{\ast }}}=1\,\mathrm{pc}$ and fixed luminosity per unit mass ${\rm{\Psi }}\equiv {L}_{* }/{M}_{* }$ typical of young, luminous clusters. We adopt a fiducial value of ${\rm{\Psi }}=2000\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$, characteristic of a fully sampled Kroupa IMF (Dopita et al. 2006).

We note that undersampling of the IMF can lead to both a large stochastic variation (da Silva et al. 2012) and a systematic overestimate in the luminosity per unit mass (Weidner & Kroupa 2006). However, as we discuss in Section 3.2, all of our clouds have high enough SFE that the IMF would be close to fully sampled by the time star formation is complete. We further discuss the effects of varying Ψ in Section 4.4.

We consider the evolution of self-gravitating star-forming clouds over a period of ∼4 initial free-fall times. Each cloud is initialized as a uniform density sphere, with ${\rho }_{0}=3{M}_{{\rm{cl,0}}}/(4\pi {r}_{0}^{3})$, where r₀ is the initial cloud radius. The clouds are centered inside cubic simulation volumes of length $L=4{r}_{0}$ with outflow boundary conditions, so that we may track the mass expelled from the cloud by radiation forces. The gas surrounding the cloud is initialized at a factor of 10³ lower than the cloud density, so that the total mass surrounding the cloud is $\sim 0.015\,{M}_{{\rm{cl,0}}}$. The density floor is lower than this again by a factor of 10, so that only a small amount of mass is added to the cloud over the course of a simulation run.

In order to cover a realistic range in cloud surface density, we adopt a range above and below fiducial values ${M}_{{\rm{cl,0}}}=5\times {10}^{4}$ and ${r}_{0}=15$ pc, which corresponds to a cloud surface density of ${{\rm{\Sigma }}}_{{\rm{cl,0}}}\equiv {M}_{{\rm{cl,0}}}/(\pi {r}_{{\rm{0}}}^{2})=70.7\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$. Our Σ-series consists of a subset of models with initial cloud masses ${M}_{{\rm{cl,0}}}\;=5\times {10}^{3},{10}^{4},2\times {10}^{4},5\times {10}^{4},{10}^{5},2\times {10}^{5}\,{M}_{\odot }$ and radii ${r}_{0}=5,8,10,15,20,25,35\,\mathrm{pc}$. In Figure 1 we show the masses, surface densities and velocity dispersions of the full set of Σ-series models compared to the same quantities derived for a set of 158 Milky Way GMCs measured by Heyer et al. (2009) and a more extended sample of 580 molecular clouds measured in Roman-Duval et al. (2010). Similar to Dale et al. (2012), we cover the high-mass end of the observed distribution, with our fiducial model being roughly characteristic of the median observed cloud. The gas motions in our clouds are initially seeded by a turbulent velocity field, with power spectrum ${v}^{2}(k)\propto {k}^{-4}$ as is observed within GMCs (e.g., Dobbs et al. 2014). The turbulence is initialized as described in Stone et al. (1998) and Skinner & Ostriker (2015). Briefly, we generate a Gaussian random field in Fourier space, such that over the range $k\in [2,64]\times {dk}$ where ${dk}=2\pi /L$, $\delta {v}_{k}$ is chosen from a Gaussian distribution with variance $P(k)\propto {k}^{-4}$. This field is then transformed back to real space and renormalized in terms of the virial parameter ${\alpha }_{{\rm{vir,0}}}\equiv 2{E}_{K}/| {E}_{G}|$ so that the variance of the velocity distribution obeys ${\sigma }^{2}=2{E}_{K}/{M}_{{\rm{cl,0}}}={\alpha }_{{\rm{vir,0}}}{E}_{G}/{M}_{{\rm{cl,0}}}$, where ${E}_{K}={M}_{{\rm{cl,0}}}{v}_{{\rm{rms}}}^{2}/2$ is the total initial turbulent gas kinetic energy, and ${E}_{G}=-3{{GM}}_{{\rm{cl,0}}}^{2}/(5{r}_{0})$ is the cloud's initial gravitational binding energy. Finally, the momentum field is forced to have zero mean by subtracting off the initial net momentum of the cloud. The initial turbulent power spectrum is a mixture of solenoidal and compressive modes. We further discuss effects of the specific initialization of turbulence in Sections 4.2.3 and 4.5. For the fiducial model, and other models in the Σ-series, we set ${\alpha }_{{\rm{vir,0}}}=2$. We also consider another series of models, the α-series, in which the initial ${\alpha }_{{\rm{vir,0}}}$ is in the range 0.1–10.0.

Zoom In Zoom Out Reset image size

Clouds with ${\alpha }_{{\rm{vir}}}=2$ are still marginally bound so long as the thermal energy does not contribute significantly to the total kinetic energy i.e., ${c}_{s}^{2}\ll 3{{GM}}_{{\rm{cl,0}}}/(5{r}_{0})$. In practice, for the lowest surface density clouds that we consider, this is satisfied provided ${c}_{s}\ll 1\,\mathrm{km}\,{{\rm{s}}}^{-1}$. We adopt a constant isothermal sound speed ${c}_{s}=0.2\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for all simulations, consistent with a temperature of T ∼ 10 K, as is characteristic of most of the mass in observed GMCs (Scoville et al. 1987). Of course, ionizing UV radiation will heat a very small fraction of the gas to a much higher temperature, and non-ionizing radiation can raise the temperature of gas within regions near stellar sources (within ${A}_{V}\sim 1$ where it is absorbed). However, given the high (FUV) optical depths of clouds, most of the gas is shielded from both internal and external radiation sources. For the regime of the present study, the optical depth to IR is small, so radiation that is absorbed and reprocessed to IR near sources subsequently escapes from the cloud without significant re-absorption by dust. To the extent that regions near sources can be heated above ∼10 K, the increase in pressure would limit fragmentation, with important consequences for the IMF (Krumholz et al. 2007; Myers et al. 2014). In addition, gas that is heated above the escape speed of the cloud could directly evaporate.

In Table 1, we list simulation inputs for our fiducial model, including the initial cloud mass, radius, and virial parameter, and the parameters ${c}_{s},\hat{c},{\rm{\Psi }}$, and κ. Our standard resolution is 256³, although we have employed higher and lower resolution grids to test convergence. We also adopt a fixed opacity $\kappa =1000\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, consistent with the radiation pressure cross sections per H derived from the Weingartner & Draine (2001) dust model (Draine 2011).

Table 1.  Fiducial Parameters

Parameter	Value
${\alpha }_{{\rm{vir,0}}}$	2.0
${r}_{{\rm{0}}}$	$15\,{\rm{pc}}$
${M}_{{\rm{cl,0}}}$	$5\times {10}^{4}\,{M}_{\odot }$
${{\rm{\Sigma }}}_{{\rm{cl,0}}}$	$70.74\,{M}_{\odot }\,{{\rm{pc}}}^{-2}$
${t}_{{\rm{ff,0}}}$	$4.29\,{\rm{Myr}}$
${v}_{{\rm{rms}}}$	$4.16\,\mathrm{km}\,{{\rm{s}}}^{-1}$
${v}_{{\rm{esc}}}$	$5.36\,\mathrm{km}\,{{\rm{s}}}^{-1}$
cs	$0.2\,{\rm{km}}\,{{\rm{s}}}^{-1}$
$\hat{c}$	$250\,{\rm{km}}\,{{\rm{s}}}^{-1}$
Ψ	$2000\,{\rm{erg}}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$
κ	$1000\,{{\rm{cm}}}^{2}\,{{\rm{g}}}^{-1}$

Download table as:  ASCIITypeset image

We extend this in Table 2, to show ${M}_{{\rm{cl,0}}}$, ${r}_{{\rm{0}}}$, and ${\alpha }_{{\rm{vir,0}}}$ for all models in our Σ and α-series. For each simulation, we show also the initial surface density ${{\rm{\Sigma }}}_{{\rm{cl,0}}}\equiv {M}_{{\rm{cl,0}}}/(\pi {r}_{0}^{2})$, the initial rms velocity dispersion ${v}_{{\rm{rms}}}={[3{\alpha }_{{\rm{vir,0}}}{{GM}}_{{\rm{cl,0}}}/(5{r}_{0})]}^{1/2}$, the escape speed at the edge of the cloud ${v}_{{\rm{esc}}}={(2{{GM}}_{{\rm{cl,0}}}/{r}_{0})}^{1/2}$, the initial gravitational free-fall time ${t}_{{\rm{ff,0}}}={[3\pi /(32G{\rho }_{0})]}^{1/2}$ and the initial hydrogen number density ${n}_{0}={\rho }_{0}/(1.4{m}_{p})$, where ${\rho }_{0}=3{M}_{{\rm{cl,0}}}/(4\pi {r}_{0}^{3})$ is the initial gas density and we allow for 40% helium by mass. Model Σ-M5E4-R15 (the same as α-A2.0), shown in bold, is the fiducial model.

Table 2.  Model Parameters

Model	${{\rm{\Sigma }}}_{{\rm{cl,0}}}$	${M}_{{\rm{cl,0}}}$	${r}_{{\rm{0}}}$	${n}_{{\rm{H,0}}}$	${t}_{{\rm{ff,0}}}$	${v}_{{\rm{rms}}}$	${v}_{{\rm{esc}}}$	${\alpha }_{{\rm{vir,0}}}$
	(${M}_{\odot }\,{{\rm{pc}}}^{-2}$)	(${M}_{\odot }$)	(${\rm{pc}}$)	(${{\rm{cm}}}^{-3}$)	(${\rm{Myr}}$)	($\mathrm{km}\,{{\rm{s}}}^{-1}$)	($\mathrm{km}\,{{\rm{s}}}^{-1}$)
Σ-M2E4-R25	10.19	$2\times {10}^{4}$	25	8.860	14.6	2.04	2.63	2.0
Σ-M5E4-R35	12.99	$5\times {10}^{4}$	35	8.072	15.3	2.72	3.51	2.0
Σ-M2E4-R20	15.92	$2\times {10}^{4}$	20	17.30	10.5	2.28	2.94	2.0
Σ-M5E4-R25	25.46	$5\times {10}^{4}$	25	22.15	9.24	3.22	4.16	2.0
Σ-M1E5-R35	25.98	$1\times {10}^{5}$	35	16.14	10.8	3.85	4.97	2.0
Σ-M2E4-R15	28.29	$2\times {10}^{4}$	15	41.02	6.79	2.63	3.39	2.0
Σ-M1E4-R10	31.83	$1\times {10}^{4}$	10	69.22	5.23	2.28	2.94	2.0
Σ-M5E4-R20	39.79	$5\times {10}^{4}$	20	43.26	6.61	3.60	4.65	2.0
Σ-M1E4-R08	49.74	$1\times {10}^{4}$	8	135.2	3.74	2.54	3.29	2.0
Σ-M1E5-R25	50.93	$1\times {10}^{5}$	25	44.30	6.53	4.55	5.88	2.0
Σ-M2E5-R35	51.97	$2\times {10}^{5}$	35	32.29	7.65	5.44	7.02	2.0
Σ-M5E3-R05	63.66	$5\times {10}^{3}$	5	276.9	2.61	2.28	2.94	2.0
Σ-M2E4-R10	63.66	$2\times {10}^{4}$	10	138.4	3.70	3.22	4.16	2.0
${\boldsymbol{\Sigma }}$-M5E4-R15	70.74	${\bf{5}}\times {{\bf{10}}}^{4}$	15	102.5	4.29	4.16	5.36	2.0
Σ-M1E5-R20	79.58	$1\times {10}^{5}$	20	86.52	4.67	5.09	6.57	2.0
Σ-M2E4-R08	99.47	$2\times {10}^{4}$	8	270.4	2.64	3.60	4.65	2.0
Σ-M2E5-R25	101.9	$2\times {10}^{5}$	25	88.60	4.62	6.44	8.31	2.0
Σ-M1E4-R05	127.3	$1\times {10}^{5}$	5	553.7	1.85	3.22	4.16	2.0
Σ-M5E5-R35	129.9	$5\times {10}^{5}$	35	80.72	4.84	8.60	11.11	2.0
Σ-M1E5-R15	141.5	$1\times {10}^{5}$	15	205.1	3.04	5.88	7.59	2.0
Σ-M5E4-R10	159.2	$5\times {10}^{4}$	10	346.1	2.34	5.09	6.57	2.0
Σ-M2E5-R20	159.2	$2\times {10}^{5}$	20	173.0	3.31	7.20	9.29	2.0
Σ-M5E4-R08	248.7	$5\times {10}^{4}$	8	676.0	1.67	5.69	7.35	2.0
Σ-M2E4-R05	254.6	$2\times {10}^{4}$	5	1107	1.31	4.55	5.88	2.0
Σ-M2E5-R15	282.9	$2\times {10}^{5}$	15	410.2	2.15	8.31	10.7	2.0
α-A0.1	70.74	$5\times {10}^{4}$	15	102.5	4.29	0.93	5.36	0.1
α-A0.2	70.74	$5\times {10}^{4}$	15	102.5	4.29	1.32	5.36	0.2
α-A0.4	70.74	$5\times {10}^{4}$	15	102.5	4.29	1.86	5.36	0.4
α-A0.8	70.74	$5\times {10}^{4}$	15	102.5	4.29	2.63	5.36	0.8
α-A1.5	70.74	$5\times {10}^{4}$	15	102.5	4.29	3.60	5.36	1.5
${\boldsymbol{\alpha }}$-A2.0	70.74	${\bf{5}}\times {{\bf{10}}}^{4}$	15	102.5	4.29	4.16	5.36	2.0
α-A3.0	70.74	$5\times {10}^{4}$	15	102.5	4.29	5.09	5.36	3.0
α-A6.0	70.74	$5\times {10}^{4}$	15	102.5	4.29	7.21	5.36	6.0
α-A10.0	70.74	$5\times {10}^{4}$	15	102.5	4.29	9.30	5.36	10.0

Note. Columns display the following information: (i) model name, (ii) initial cloud surface density, (iii) initial cloud mass, (iv) initial cloud radius, (v) initial cloud hydrogen number density, assuming a mean atomic weight of $\mu =1.4$, (vi) initial cloud free-fall time, (vi) initial turbulent velocity dispersion, (viii) cloud escape velocity from the initial cloud radius, and (ix) initial virial parameter. The fiducial model is shown in bold (Σ-M5E4-R15 and α-A2.0).

Download table as:  ASCIITypeset image

Table 3.  Model Results

Model	ε	${\varepsilon }_{{\rm{adj}}}$	${t}_{{\rm{\ast }}}$	${t}_{{\rm{50}}}$	${t}_{90}$	${t}_{{\rm{unb}}}$	β	${t}_{{\rm{break}}}$	${\varepsilon }_{\mathrm{ff},\bar{\rho }}$
			(${t}_{{\rm{ff}},0}$)	(${t}_{{\rm{ff}},0}$)	(${t}_{{\rm{ff}},0}$)	(${t}_{{\rm{ff}},0}$)		(${t}_{{\rm{ff}},0}$)
Σ-M2E4-R25	0.12	0.14	0.67	0.88	1.51	${1.14}_{-0.24}^{+0.18}$	⋯a	⋯a	⋯a
Σ-M5E4-R35	0.18	0.20	0.63	1.07	2.05	${1.17}_{-0.23}^{+0.22}$	⋯a	⋯a	⋯a
Σ-M2E4-R20	0.22	0.25	0.61	1.08	1.72	${1.35}_{0.12}^{+0.25}$	⋯a	⋯a	⋯a
Σ-M5E4-R25	0.25	0.28	0.52	1.01	1.71	${1.18}_{-0.20}^{+0.16}$	1.64	0.80	0.49
Σ-M1E5-R35	0.23	0.26	0.49	1.00	1.53	${1.22}_{-0.22}^{+0.19}$	0.77	0.85	0.24
Σ-M2E4-R15	0.30	0.34	0.47	1.07	1.65	${1.21}_{-0.19}^{+0.12}$	1.08	0.70	0.37
Σ-M1E4-R10	0.30	0.34	0.54	1.09	1.64	${1.29}_{-0.20}^{+0.16}$	1.11	0.78	0.34
Σ-M5E4-R20	0.32	0.36	0.46	1.02	1.51	${1.19}_{-0.16}^{+0.14}$	1.17	0.73	0.45
Σ-M1E4-R08	0.39	0.44	0.47	1.11	1.75	${1.34}_{-0.20}^{+0.16}$	0.85	0.92	0.41
Σ-M1E5-R25	0.32	0.37	0.46	1.05	1.48	${1.26}_{-0.19}^{+0.14}$	1.25	1.02	0.36
Σ-M2E5-R35	0.31	0.35	0.47	1.08	1.51	${1.32}_{-0.24}^{+0.16}$	0.87	0.74	0.29
Σ-M5E3-R05	0.42	0.49	0.54	1.17	1.95	${1.45}_{-0.21}^{+0.17}$	0.69	1.02	0.37
Σ-M2E4-R10	0.42	0.48	0.44	1.08	1.58	${1.33}_{-0.22}^{+0.16}$	0.93	0.76	0.43
Σ-M5E4-R15	0.42	0.48	0.37	1.06	1.57	${1.33}_{-0.21}^{+0.14}$	0.95	0.81	0.43
Σ-M1E5-R20	0.41	0.47	0.32	1.02	1.56	${1.31}_{-0.21}^{+0.18}$	0.88	0.84	0.45
Σ-M2E4-R08	0.49	0.56	0.42	1.10	1.71	${1.46}_{-0.24}^{+0.18}$	0.94	0.80	0.49
Σ-M2E5-R25	0.41	0.41	0.32	1.04	1.57	${1.43}_{-0.24}^{+0.20}$	0.87	0.92	0.42
Σ-M1E4-R05	0.52	0.60	0.44	1.16	1.80	${1.59}_{-0.25}^{+0.23}$	0.83	0.91	0.51
Σ-M1E5-R15	0.52	0.59	0.29	1.04	1.63	${1.40}_{-0.13}^{+0.22}$	1.31	1.14	0.45
Σ-M5E4-R10	0.58	0.65	0.31	1.09	1.67	${1.55}_{-0.28}^{+0.21}$	0.95	0.78	0.49
Σ-M2E5-R20	0.48	0.54	0.28	1.21	1.86	${1.43}_{-0.18}^{+0.31}$	1.05	0.63	0.51
Σ-M5E4-R08	0.62	0.70	0.31	1.11	1.77	${1.73}_{-0.30}^{+0.28}$	1.15	0.78	0.49
Σ-M2E4-R05	0.61	0.69	0.37	1.20	1.86	${1.83}_{-0.32}^{+0.34}$	1.11	0.78	0.49
Σ-M2E5-R15	0.61	0.69	0.26	1.08	1.63	${1.75}_{-0.25}^{+0.25}$	1.30	0.56	0.47
α-A0.1	0.91	0.91	0.67	1.02	1.13	${1.16}_{-0.13}^{+0.04}$	1.24	0.83	0.42
α-A0.2	0.87	0.87	0.60	1.04	1.22	${1.27}_{-0.19}^{+0.06}$	1.19	0.78	0.43
α-A0.4	0.67	0.67	0.51	0.99	1.23	${1.28}_{-0.17}^{+0.07}$	1.39	0.75	0.44
α-A0.8	0.58	0.59	0.44	1.00	1.31	${1.28}_{-0.15}^{+0.10}$	0.98	0.77	0.36
α-A1.5	0.47	0.54	0.43	1.03	1.54	${1.29}_{-0.19}^{+0.17}$	1.10	0.70	0.41
α-A2.0	0.42	0.48	0.37	1.06	1.57	${1.33}_{-0.21}^{+0.14}$	1.03	0.80	0.43
α-A3.0	0.28	0.35	0.36	0.99	1.47	${1.33}_{-1.33}^{+0.18}$	0.92	0.80	0.36
α-A6.0	0.12	0.35	0.27	1.03	1.66	⋯b	1.57	0.87	0.16
α-A10.0	0.05	0.32	0.14	1.00	1.39	⋯b	0.20	1.95	0.07

Notes. Columns display the following information: (i) final efficiency of star formation, (ii) final efficiency adjusted for the inital turbulence-driven mass outflow, (iii) time of first star formation, (iv) time at which 50% of the final stellar mass is assembled, (v) time at which $90 \%$ of the final stellar mass is assembled, (vi) time at which the cloud becomes unbound (reaches a virial parameter of ${\alpha }_{{\rm{vir}}}={5}_{-3}^{+5}$), (vii) post-break power-law exponent of the stellar mass evolution using a double power-law fit, (viii) break time in double power-law fit to the stellar mass evolution, and (ix) star formation efficiency per free-fall time as defined in Equation (25).

^aIn these models, the star formation rate was not fit as there were insufficient numbers of discrete star particles to perform a meaningful fit. ^bIn these models, ${t}_{{\rm{unb}}}$ was not calculated since the models were initially unbound.

Download table as:  ASCIITypeset image

We begin by considering the overall time evolution of our fiducial model. In Figures 2 and 3 we show evolving column density maps of this cloud in the y–z and x–y planes, respectively. We see that turbulence drives the cloud into collapse very rapidly, and that by a little before half a free-fall time, mass has gathered preferentially along two perpendicular filaments, roughly coincident with the x and z-axes. The first star formation event occurs around this time, as the density peaks in the cloud continue to collapse under self-gravity. For the fiducial model and all others, the time of the first star formation t_∗ is given in Table 3. The filamentary structure of gas in the cloud and the shape of its density distribution do not, however, change significantly. By ∼t₁₀, where t_x denotes the time at which x% of the stellar mass is assembled, the main difference from earlier is that the density contrast between the filaments and surrounding gas has increased. As star formation progresses, the UV radiation begins to drive gas in even the densest filaments away from the sites of star formation. By the time half the stars have formed, much of the gas, even at higher densities, is already flowing outwards from the center of mass. At late stages, shown for example in the surface density projections at $t=2\,{t}_{{\rm{ff,0}}}$ in Figure 4, all of the material is streaming away from the center, which has been cleared of gas.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the same picture as in Figure 3 but for a single slice through the x–y plane with the radiation field flux directions and energy density overplotted in vectors and contours, respectively. Because of the filamentary nature of the cloud, even though stars form near the center of the original GMC, the radiation from these stars very quickly blows a hole in the surrounding gas. Therefore, by t₁₀, as shown in the second snapshot, a significant fraction of the radiation already escapes the cloud through the second quadrant, where the gas density is small.

Zoom In Zoom Out Reset image size

In detail, the gas structure surrounding the most massive star clusters is far from smooth and presents a very different picture from the gas expansion seen in outflows driven by reprocessed radiation (Skinner & Ostriker 2015). At late stages in Figures 2–5 we see a central cavity surrounded by prominent high-density fingers of gas extending inward toward the most massive stars. Regions of low column density, as seen by the central sources, are evacuated first once the radiative force on them becomes super-Eddington. However, regions that are shielded by higher-density clumps of gas only begin to be driven out at late times, giving rise to the prominent columns of gas.

This picture of filamentary collapse followed by rapid star formation and subsequent gas expulsion can also be seen in the time histories shown in Figure 6. In Figure 6(a), we show the evolution of the total mass in stars, gas, and outflows from the simulation volume. While the cloud is initially collapsing there is relatively little star formation, although a small number of stars are formed through the effects of turbulence initiating compression. Meanwhile, the outflows driven by this initial turbulence only start leaving the simulation volume at $\sim 0.8{t}_{{\rm{ff,0}}}$, which is roughly the time taken for gas traveling at ∼2 times the escape velocity to reach the corner of the box. Table 3 lists the times t₅₀ and t₉₀ when, respectively, 50% and 90% of the stellar mass has been assembled.

Zoom In Zoom Out Reset image size

At around the same time ($\sim 0.8{t}_{{\rm{ff,0}}}$) there is a break in the stellar mass evolution as stars begin forming more rapidly; the majority of the stellar mass is assembled over the next free-fall time. By $\sim 1.8{t}_{{\rm{ff,0}}}$ the accretion onto star particles is essentially complete, while radiation from these stars continues to accelerate the remaining gas so that it becomes unbound from the central cluster. By $\sim 3{t}_{{\rm{ff,0}}}$ most of the outflowing gas has left the simulation box.

The evolution of the global gas distribution is more difficult to characterize. One simple measure is the global virial parameter, which is a rough proxy for the collapse and expansion of the cloud and is shown in Figure 6(b). Initially this decreases from its starting value of ${\alpha }_{{\rm{vir}}}=2$ as the turbulence decays slightly and the total potential energy increases due to contraction along local filaments. However, as for the stellar mass, there is a break at around $\sim 0.8{t}_{{\rm{ff,0}}}$, where collapse ceases and radiation from the first star particles begins to unbind gas and drive outflows.

The evolution in Mach number parallels that of the virial parameter very closely. In fact, the two are well correlated until at least around $\sim 2\,{t}_{{\rm{ff,0}}}$, with ${ \mathcal M }\propto {\alpha }_{{\rm{vir}}}^{0.4}$. This roughly corresponds to the evolution in Mach number for a shell expanding around a fixed central mass M. For a shell of fixed mass ${M}_{{\rm{sh}}}$, expanding with velocity v, the virial parameter is ${\alpha }_{{\rm{vir}}}={v}^{2}r/{GM}$, while the mass-weighted Mach number is roughly ${{ \mathcal M }}^{2}=({M}_{{\rm{sh}}}/M)({v}^{2}/{c}_{s}^{2})$ so that ${{ \mathcal M }}^{2}=({{GM}}_{{\rm{sh}}}/r)({\alpha }_{{\rm{vir}}}/{c}_{s}^{2})$. The evolution we find is slightly different from ${ \mathcal M }\propto {\alpha }_{{\rm{vir}}}^{1/2}$ since the radius of the shell expands. In general though, we see that the evolution of both Mach number and virial parameter beyond $\sim 0.8\,{t}_{{\rm{ff,0}}}$ is dominated by the cloud expansion and not turbulence.

It is for this reason that the increase in Mach number has no correlated increase in the width of the log-normal distribution (see below), unlike the case for driven turbulence when the width of the density distribution is entirely set by the Mach number (Padoan & Nordlund 2011; Molina et al. 2012; Hopkins 2013; Thompson & Krumholz 2016). While the lowest density regions, which are expanding away at high velocity, dominate the evolution of the virial parameter and the Mach number, they have little influence on the density distribution since, until late times, they only represent around 10%–20% of the mass. Therefore, the log-normal density distribution, which is fit to the majority of gas remaining in the cloud, shows no significant change over the bulk of star formation.

A measure quantifying the evolution of the gas is the shape of the gas density distribution (PDF). In Figure 7 we show the density distribution at two separate times (${t}_{10}=0.59\,{t}_{{\rm{ff,0}}}$ in black and ${t}_{50}=1.06\,{t}_{{\rm{ff,0}}}$ in red). We show only the mass distribution, since the distribution by volume comprises more than 90% empty space. We see that by $\sim 0.6\,{t}_{{\rm{ff,0}}}$, when the virial parameter begins to turn around, the gas density is roughly log-normal in shape, characteristic of supersonic turbulence (Vazquez-Semadeni 1994; Ostriker et al. 2001; Vázquez-Semadeni & García 2001; Federrath et al. 2009). Interestingly, even towards the end of star formation, the shape of the distribution is not significantly different at the high-density end. Certainly, there is an excess of low-density gas escaping the cloud, but the highest density portion, which represents the star-forming regions, remains essentially the same.

Zoom In Zoom Out Reset image size

Our fiducial simulation displays a number of key stages of evolution. Initially, there is rapid collapse and structure formation driven by turbulent compression and self-gravity, resulting in a filamentary gas distribution within around half a free-fall time. This is followed by ongoing collapse under self-gravity, during which the filamentarity remains while collapse around the density peaks begins to form stars at an accelerating rate. Starting at $t\sim 0.8{t}_{{\rm{ff,0}}}$ there is a transition from collapse dominated by self-gravity to cloud expansion driven by radiative feedback, where the lowest density regions are accelerated out of the cloud first, but where the star-forming regions remain largely unaffected. Star formation is largely complete by $\sim 1.5{t}_{{\rm{ff,0}}}$, and the remnant gas in the cloud is mostly removed by radiation-driven outflows within the next free-fall time.

In numerically simulating turbulent cloud evolution with radiation feedback, we want to ensure that (i) we accurately model the formation of stars driven by turbulence and gravitational collapse and (ii) we accurately model gas dynamics driven by radiative forces.

The first of these is related to the methodology of sink-particle creation, which is necessary because gas collapse becomes unresolved on the numerical grid. To have a sufficiently "fine-grained" representation of star formation, we would like to ensure that a single star particle cannot represent a substantial fraction of the cloud. The minimum sink-particle mass is defined in terms of the density threshold and the cell size ${\rm{\Delta }}x=L/N=4{r}_{0}/N$ as

$$\begin{eqnarray}&&{M}_{{\rm{sink,min}}}={\rho }_{{\rm{th}}}{({\rm{\Delta }}x)}^{3}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{8.86{c}_{s}^{2}{\rm{\Delta }}x}{G\pi }=\displaystyle \frac{35.4{c}_{s}^{2}{r}_{0}}{G\pi N}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&=24.5\,{M}_{\odot }{\left(\displaystyle \frac{{c}_{s}}{0.2\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{2}\left(\displaystyle \frac{{r}_{0}}{10\,{\rm{pc}}}\right)\left(\displaystyle \frac{256}{N}\right).\end{eqnarray} \tag{ 8 }$$

As a practical limit, we require that this mass is not more than $\sim 0.1 \%$ of the initial cloud mass, so that a reasonable number of star particles can be created, even at low efficiencies. This translates to around ${\rm{\Sigma }}\gtrsim 10\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ for the masses and radii that we consider, though we will examine this limit in more detail below.

We also need to ensure that we can sufficiently resolve the highest density regions, which are converted into star particles. High-resolution simulations of compressible turbulence modeling molecular cloud conditions generally show both a log-normal component around the mean density, characteristic of supersonic turbulence (Vazquez-Semadeni 1994; Ostriker et al. 2001; Vázquez-Semadeni & García 2001; Federrath et al. 2009), and an extended power-law tail that arises when self-gravity is included (Klessen et al. 2000; Dib & Burkert 2005; Federrath et al. 2008; Vázquez-Semadeni et al. 2008; Ballesteros-Paredes et al. 2011; Kritsuk et al. 2011; Collins et al. 2012; Federrath & Klessen 2013). There is evidence that the power-law component is specifically associated with dense prestellar cores (see below). High-resolution observations of star-forming clouds also generally show power-law tails in the surface density whenever there is star formation (Kainulainen et al. 2009; Schneider et al. 2013, 2015), although there is some debate as to whether or not the power law arises from self-gravity (Brunt 2015) and what is the exact relationship between the surface density and density distributions (Brunt et al. 2010).

By contrast, when we consider the mass-weighted density distribution in our simulations, we find no strong evidence for a power-law tail. In other simulations, this tail is associated with self-gravity and, in particular, with collapsing cores, which themselves have power-law density profiles (Kritsuk et al. 2011; Lee et al. 2015). Thus, its absence could indicate that we do not capture the details of core collapse. In general, these power-law tails are seen to begin at densities ∼100 times the mean density (Dib & Burkert 2005; Vázquez-Semadeni et al. 2008; Kritsuk et al. 2011; Collins et al. 2012; Federrath & Klessen 2013; Lee et al. 2015). In our fiducial simulation, the threshold for star particle formation is at $\rho \sim 150{\rho }_{0}$ (or ${n}_{H}\sim 1.5\times {10}^{4}\,{\mathrm{cm}}^{-3}$). With the peak of the log-normal in Figure 7 at ${\rm{ln}}({\rho }_{{\rm{peak}}}/{\rho }_{0})\sim 2.2$ (or ${n}_{H}\sim {10}^{3}\,{\mathrm{cm}}^{-3}$) and a variance of ${\sigma }_{{\rm{ln}}\rho }\sim 2$, the $1\sigma$ density will be at ${n}_{H}\sim {10}^{4}\,{\mathrm{cm}}^{-3};$ hence there is only a limited range of density over which to sample the power law. Effectively, the majority of what would make up the power-law core is hidden in these star particles. Nevertheless, tests at higher resolution (see below) appear to indicate that all our models are converged in the SFE.

In order to test how well our code captures the physics of gas expansion, we conducted a number of tests involving a spherical gas shell surrounding a single central star particle, and the results of these are shown in the Appendix. These tests suggest that the primary numerical limit in simulating the interaction between radiation and gas with high accuracy is in ensuring that ${\tau }_{{\rm{cell}}}=\rho \kappa {\rm{\Delta }}x\lesssim 2$ for the optical depth within individual cells. This corresponds to a maximum density that is resolution dependent:

$$\begin{eqnarray}&&{n}_{{\rm{H,}}{\rm{\max }}}=2.8\times {10}^{3}\,{\mathrm{cm}}^{-3}{\left(\displaystyle \frac{{\rm{\Delta }}x}{0.1{\rm{pc}}}\right)}^{-1}\\ &&\times {\left(\displaystyle \frac{\kappa }{1000{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{-1}\left(\displaystyle \frac{{\tau }_{{\rm{cell,max}}}}{2}\right).\end{eqnarray} \tag{ 9 }$$

To gain an idea of how the optical depth evolves in the fiducial model, we show in Figure 8 snapshots of the spatial distribution of ${\tau }_{{\rm{cell}}}$ at the same four times as Figures 2–5. We see that even in this model, which has an initial surface density of ${{\rm{\Sigma }}}_{{\rm{cl,0}}}=71\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ and ${\tau }_{{\rm{cell,init}}}=3{{\rm{\Sigma }}}_{0}\kappa /N=0.17$, over time the densest regions are enhanced by at least two orders of magnitude and a portion of the cells have an optical depth ${\tau }_{{\rm{cell}}}\gtrsim 2$.

Zoom In Zoom Out Reset image size

Even if some localized regions have large single-cell optical depth and do not resolve the radiation field well, the overall resolution may still be acceptable. In particular, by around a free-fall time, when stars are beginning to drive gas away from the sites of star formation, only very small regions have ${\tau }_{{\rm{cell}}}\gtrsim 2$. While we may underestimate the radiation force applied to a small number of cells, for the bulk of the mass we still follow the physics governing the expulsion of gas and the propagation of the radiation field.

Quantitatively, the distribution of ${\tau }_{{\rm{cell}}}$ has a long tail, such that close to 90% of the mass is found an order of magnitude below the maximum density. In Figure 9(a) we show the $50{\rm{th}}$ percentile, 90th percentile, and maximum of the ${\tau }_{{\rm{cell}}}$ distribution as a function of time for the fiducial model. Even though the maximum ${\tau }_{{\rm{cell}}}$ reaches ∼20 at quite early times, ${\tau }_{{\rm{cell,90}}}\approx 1$ for the majority of the simulation, which is only a factor of 10 higher than the initial optical depth. A similar ratio of ${\tau }_{{\rm{cell,90}}}/{\tau }_{{\rm{cell,init}}}$ appears to hold in higher surface density models. We see in Figure 9(b) that a model with ${{\rm{\Sigma }}}_{{\rm{cl,0}}}\approx 280\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$, reaches ${\tau }_{{\rm{cell,90}}}\approx 8$. This suggests that, up to cloud surface densities ${{\rm{\Sigma }}}_{{\rm{cl,0}}}\approx 100\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$, ${\tau }_{{\rm{cell,90}}}$ will remain below ∼2, while for higher surface density clouds ${\tau }_{{\rm{cell,90}}}$ may be higher.

Zoom In Zoom Out Reset image size

The above considerations suggest that we can satisfactorily simulate clouds over a range of initial surface densities $10\,{M}_{\odot }\,{\mathrm{pc}}^{-2}\lesssim {{\rm{\Sigma }}}_{{\rm{cl,0}}}\lesssim 100\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$. However, in reality we may be able to do better than that, since as discussed in the Appendix, even at higher cell optical depths, we only underestimate the velocity of gas expulsion by around ∼10%. We are primarily interested in ensuring that we are able to accurately capture the net SFE and evolution in response to radiation. By looking at how these properties vary as we change the numerical resolution, we may gain a better idea of what range of parameters provide converged results.

Figure 10(a) shows evolution of the global efficiency of star formation in the fiducial model as a function of time with varying numerical resolution for N = 128, 256, and 512. We see that in the N = 128 case, the fiducial simulation is clearly not converged. It shows substantial differences from the higher resolution models, including "stair-stepping" in the stellar mass history. This makes sense, since the minimum sink mass for this simulation is around ${M}_{{\rm{sink}}}\approx 80\,{M}_{\odot }$, which is close to 0.2% of the global cloud mass. Importantly, however, the N = 256 and N = 512 models are converged, and even though the N = 128 model underestimates the SFE at intermediate times, its final SFE is the same as for the higher resolution models.

Zoom In Zoom Out Reset image size

We show also in Figure 11 the density distribution at t₅₀, when 50% of the final stellar mass has been assembled. With resolution reduced by a factor of two in the 128³ models, the star particle density threshold is at $\sim 40{\rho }_{0}$ or ${n}_{H}\sim 4\times {10}^{3}\,{\mathrm{cm}}^{-3}$ compared to $\rho \sim 150{\rho }_{0}$ or ${n}_{H}\sim 1.5\times {10}^{4}\,{\mathrm{cm}}^{-3}$ for $N={256}^{3}$. We believe that the lack of convergence evident in the low-resolution models (see Figure 10), which leads to an underestimation of the early SFE, is because this threshold is not sufficiently high compared to the log-normal distribution, as evident from the cutoff in Figure 11 here. However, the convergence of the 256³ and 512³ models in terms of stellar efficiency (see Figure 10) would suggest that in our highest resolution models, we capture the distribution of gas dominated by turbulence, with only collapsing regions assigned to sink particles. Thus, while insufficient resolution can lead to an underestimate of the SFE, we believe that increasing resolution beyond a certain point primarily provides improved resolution of collapsing cores without changing the SFE. This speculation could be tested with AMR simulations.

Zoom In Zoom Out Reset image size

A similar picture holds for models at the highest surface densities we simulate. In Figure 12(a) we again show evolution of the global stellar efficiency as a function of resolution, but in this case for a cloud of mass ${M}_{{\rm{cl,0}}}=2\times {10}^{5}\,{M}_{\odot }$. For this model, the 90th percentile in optical depth lies at ${\tau }_{{\rm{cell,90}}}\sim 8$, which is a factor of four larger than the limit found in the Appendix for best accuracy in capturing radiation forces. However, the larger values of ${\tau }_{{\rm{cell}}}$ do not have a major impact on the evolution of the global efficiency, presumably because errors are only ∼10% at larger ${\tau }_{{\rm{cell}}}$ (see the Appendix), and because only a small fraction of the gas has large ${\tau }_{{\rm{cell}}}$. As for the fiducial cloud, stars do begin forming slightly earlier in the higher resolution run and they do form at a slightly faster rate, but the final efficiency is the same.

Zoom In Zoom Out Reset image size

In Figures 10(b) and 12(b), we show the gas virial parameter for our fiducial and highest-Σ models at different resolution. As for the stellar efficiency, our simulations are not well converged at N = 128, but there are few differences between the N = 256 and N = 512 runs, at least until the bulk of star formation is already complete.

We begin by presenting results from a set of models with the initial turbulent energy a factor of 20 lower than in the fiducial case (i.e., with ${\alpha }_{{\rm{vir,0}}}=0.1$). In Figure 13, we show snapshots of the column density in the x–y plane for a model cloud with the fiducial initial surface density, but with reduced initial turbulence. The lack of turbulent support means that for the first $\sim 0.6{t}_{{\rm{ff,0}}}$, the cloud undergoes nearly free-fall collapse, contracting to around 40% of its initial radius and converting potential to kinetic energy. Star formation only begins once this contraction is complete and the virial parameter for the cloud is near unity.

Zoom In Zoom Out Reset image size

After the initial contraction and virialization, the evolution is quite similiar to the fiducial case as the gas density distribution becomes similarly filamentary, with preferential directions set by the seed turbulent field. Low initial turbulence simplifies the picture in several respects. First, stars tend to form very close to the center of the cloud, with very little initial momentum away from the center of mass. This can be seen in the lower two panels of Figure 13, where star particles close to the cloud center eventually converge to a single, massive, central star particle. Second, star formation occurs on a relatively short timescale compared to the initial free-fall time due to the high densities and small physical length scales involved. Therefore, the majority of stellar mass is assembled in less than 25% of a free-fall time after the initial collapse.

As a very simple model of star formation, one could imagine an initial roughly spherical collapse, until the cloud reaches a characteristic radius xr₀ at which the virial parameter is of the order of unity. Once stars have begun to form close to the center and have evacuated their local environment, we might expect the remnant cloud's effective radius to evolve under the competing effects of the inward force of gravity and the outward force of radiation. Star formation, or at least accretion onto the young central star cluster would then continue so long as gravity dominates over the effects of stellar feedback, and it would stop once enough stars form for radiation forces to disperse the remaining gas in the cloud.

This picture is essentially that of a radiation-dominated H ii region, with a single, central stellar cluster that stops accreting mass once the luminosity is sufficient to drive away all surrounding gas. Such systems, in different limits, have been analyzed previously in Elmegreen (1983), Scoville et al. (2001), Krumholz & Matzner (2009), Murray et al. (2010), Fall et al. (2010), and Kim et al. (2016), and considering effects of inhomogeneity by Thompson & Krumholz (2016). We review the ideas involved, connecting to the results of our simulations.

We consider a system consisting of a spherical shell of gas with mass ${M}_{{\rm{sh}}}=(1-\varepsilon ){M}_{{\rm{cl,0}}}$ surrounding a point mass representing a stellar cluster of mass ${M}_{* }=\varepsilon {M}_{{\rm{cl,0}}}$. At any point in time, the circumcluster shell will have a distribution of surface densities around the mean surface density $\langle {{\rm{\Sigma }}}^{c}(\varepsilon ,x,{{\rm{\Sigma }}}_{{\rm{cl,0}}})\rangle$, where ${{\rm{\Sigma }}}_{{\rm{cl,0}}}={M}_{{\rm{cl,0}}}/\pi {r}_{0}^{2}$ is the initial, observed cloud column density. The mean surface density would decrease as gas is converted to stars, but increase under global collapse, so the outcome depends on the initial cloud reference surface density ${{\rm{\Sigma }}}_{{\rm{cl,0}}}$, the time-dependent stellar efficiency ε, and the radius to which the cloud collapses xr₀. We note also that here $\langle {{\rm{\Sigma }}}^{c}\rangle ={M}_{{\rm{sh}}}/(4\pi {r}^{2})$ represents the mean surface density seen by the central star cluster, as distinct from the cloud column density seen by an external observer ${M}_{{\rm{sh}}}/(\pi {r}^{2})=4\langle {{\rm{\Sigma }}}^{c}\rangle$.

In this simple model, assuming all the flux is absorbed, the radiative force per unit area in the shell due to the central luminous source is ${F}_{{\rm{rad}}}=L/(4\pi {r}^{2}c)={\rm{\Psi }}{M}_{* }/(4\pi {r}^{2}c)$. Meanwhile, there are two components to the gravitational force applied to a local patch of surface density ${{\rm{\Sigma }}}^{c}$. First, there is the force per unit area due to the central point mass ${F}_{1}={{GM}}_{* }{{\rm{\Sigma }}}^{c}/{r}^{2}$. Second, unlike previous studies (Fall et al. 2010), we also include a term for the self-gravity of the shell per unit area, which is ${F}_{2}={{GM}}_{{\rm{sh}}}{{\rm{\Sigma }}}^{c}/(2{r}^{2})$. The exact normalization of this latter force is derived for a uniform density shell and may change with surface density variations. However, with the basic argument that gas at the radial center of the shell will feel the force of all interior mass, or half the total shell mass, this is at least approximately correct even when the shell surface density is not uniform.

For a region with surface density ${{\rm{\Sigma }}}^{c}$, the Eddington ratio ${F}_{{\rm{rad}}}/{F}_{{\rm{grav}}}={F}_{{\rm{rad}}}/({F}_{1}+{F}_{2})$ is equal to ${\rm{\Psi }}{M}_{* }/[2\pi {cG}{{\rm{\Sigma }}}^{c}(2{M}_{* }+{M}_{{\rm{sh}}})]$. At a time when the net SFE is ε, we can evaluate the Eddington surface density ${{\rm{\Sigma }}}_{E}$ such that the inward and outward forces balance and the Eddington ratio is unity:

$$\begin{eqnarray}{{\rm{\Sigma }}}_{E} & = & \displaystyle \frac{\varepsilon }{\varepsilon +1}\displaystyle \frac{{\rm{\Psi }}}{2\pi {cG}}\\ & = & 759\,{M}_{\odot }\,{\mathrm{pc}}^{-2}\displaystyle \frac{\varepsilon }{1+\varepsilon }\left(\displaystyle \frac{{\rm{\Psi }}}{2000\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}}\right);\end{eqnarray} \tag{ 10 }$$

here we have used ${M}_{* }=\varepsilon {M}_{{\rm{cl,0}}}$ and ${M}_{{\rm{sh}}}=(1-\varepsilon ){M}_{{\rm{cl,0}}}$. In local patches where surface densities within the shell exceed the Eddington level ${{\rm{\Sigma }}}^{c}\gt {{\rm{\Sigma }}}_{E}$, gas will be able to continue collapsing, and where ${{\rm{\Sigma }}}^{c}\lt {{\rm{\Sigma }}}_{E}$, the gas can be driven outward. We note that even for $\varepsilon \to 1$, ${{\rm{\Sigma }}}_{E}\lt 400{M}_{\odot }\,{\mathrm{pc}}^{-2}$, implying that UV radiation feedback by itself is not expected to be effective in expelling gas from very high surface density GMCs (see Skinner & Ostriker 2015, for a study of reprocessed radiation effects in this regime).

We now consider a very simple hypothesis, similar in spirit to Fall et al. (2010), in which star formation and accretion halts completely once the mean circumcluster surface density reaches the Eddington value, $\langle {{\rm{\Sigma }}}^{c}\rangle ={{\rm{\Sigma }}}_{E}$. For a given $\varepsilon (t)$, the mean circumcluster surface density can be related to the initial cloud mass and surface density by

$$\begin{eqnarray}\langle {{\rm{\Sigma }}}^{c}\rangle & = & \displaystyle \frac{(1-\varepsilon ){M}_{{\rm{cl,0}}}}{4\pi {({{xr}}_{0})}^{2}}\\ & = & \displaystyle \frac{(1-\varepsilon )}{4{x}^{2}}{{\rm{\Sigma }}}_{{\rm{cl,0}}},\end{eqnarray} \tag{ 11 }$$

where the $(1-\varepsilon )$ comes from the conversion of gas to stars, the factor of x from cloud contraction, and the additional factor of four from the fact that ${{\rm{\Sigma }}}_{{\rm{cl,0}}}={M}_{{\rm{cl,0}}}/(\pi {r}_{0}^{2})$ is an observed column density. If we substitute $\langle {{\rm{\Sigma }}}^{c}\rangle$ from Equation (11) for ${{\rm{\Sigma }}}_{E}$ in Equation (10), we obtain a simple relationship between the SFE and the initial cloud surface density:

$$\begin{eqnarray}\displaystyle \frac{1}{\varepsilon }-\varepsilon & = & \displaystyle \frac{2{\rm{\Psi }}{x}^{2}}{\pi {cG}{{\rm{\Sigma }}}_{{\rm{cl,0}}}}\\ & = & 30.5{x}^{2}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{cl,0}}}}{100{M}_{\odot }{\mathrm{pc}}^{-2}}\right)}^{-1}\left(\displaystyle \frac{{\rm{\Psi }}}{2000\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}}\right).\end{eqnarray} \tag{ 12 }$$

Similar to Equation (6) of Fall et al. (2010), this predicts a low efficiency ($\varepsilon \sim 3 \%$) of star formation in GMCs with ${\rm{\Sigma }}\sim 100\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$, typical of the Milky Way.

To test whether or not this simple radiative force balance argument can capture simulated cloud behavior, we first consider low-turbulence models, ${\alpha }_{{\rm{vir,0}}}=0.1$, for a range of surface densities. In Figure 14 we show the final SFE for this set of simulations, which have initial masses and radii that are the same as for our Σ-series.

Zoom In Zoom Out Reset image size

In Figure 14 we also show the fit of the simple model represented by Equation (12) to our cloud simulations. In the low-turbulence case, the shape of Equation (12) fits relatively well, matching the drop in efficiencies in our low surface density clouds. However, the best fit requires a value of ${x}^{2}=0.005$, or linear contraction of the cloud radius by $x\sim 0.07$. Modeling the cloud contraction is not straightforward due to the non-spherical gas distribution. Even so, it is clear from looking at Figure 13 that the low-turbulence clouds collapse to around $x\sim 1/2$ of their initial radii. This is a factor of more than five larger than the nominal best fit. We conclude that although Equation (12) can fit the numerical results for the reduced-turbulence models, the parameters required are not consistent with the true evolution and structure of the cloud.

The source of this discrepancy is likely that the surface density seen by the radiation field is not uniform, but is instead distributed over at least an order of magnitude. As a consequence, at a given ε, even if the mean surface density in a cloud is comparable to ${{\rm{\Sigma }}}_{E}$ as given in Equation (10), there will be many higher surface density regions that can still collapse, leading to an increase in ε. Conversely, as noted by Thompson & Krumholz (2016), this also means that low surface density regions can be driven away even when the stellar luminosity is not sufficient to drive off the bulk of a cloud's mass. To assess the impact of this, we now turn to the effect that turbulence has on changing the gas density distribution and therefore also this simple picture of star formation regulation.

Qualitatively, we expect turbulence to affect cloud evolution in two separate ways. First, since our Σ-series clouds are only marginally bound, the initial turbulence will unbind a certain amount of the gas on the edges of the cloud immediately. This effectively lowers the initial surface density, and sets a maximum on the SFE. Second, and more fundamentally, turbulence dramatically alters the density and surface density distributions.

We may gain a more concrete idea of what turbulence means for star formation by considering a simplified version of our Σ- and α-series in which radiative feedback is turned off, so that only turbulence and self-gravity set the SFR. In these simulations, star formation continues until all the gas is consumed, and the only gas that is not converted to stars is that which is initially unbound by the turbulence.

Figure 15 shows the resulting column density evolution for our fiducial model with no radiation. Until more than a free-fall time the no-feedback and feedback cases are very similar. Gas collapses along the same filaments, with the same initial density structure. Moreover, stars appear to form at roughly the same rate despite the lack of radiative feedback to stir up additional turbulence. This is likely because even though the no-feedback model does not drive out low-density gas, and so does not show the same prominent columns seen in Figures 3–5, the high-density sites of star formation remain similar. Differences between the two models appear only at late times, since in the model without feedback, star formation continues at the same rate for close to an extra free-fall time. Eventually, the majority of the no-feedback cloud is converted into stars. Only the gas driven away by the initial turbulence does not contribute to star formation.

Zoom In Zoom Out Reset image size

In fact, if we compare the evolution of stellar mass in all of our Σ and α-series models (discussed in detail below) with and without feedback, we find that they are essentially identical up until the point that radiation begins to disperse the gas. We may define a time, ${t}_{{\rm{fb}}}$, at which there is a 10% difference in stellar mass between simulations with and without radiation feedback. We then define ${f}_{{\rm{fb}}}$, the fraction of the final stellar mass in the feedback model that has been assembled by time ${t}_{{\rm{fb}}}$. In Figure 16 we show ${f}_{{\rm{fb}}}$ as a function of both initial surface density and virial parameter. We see that for almost all models, except for those at low surface density or high virial parameter (where there is little stellar mass formed), between 70% and 90% of the stellar mass is assembled before radiative feedback can significantly affect star formation.

Zoom In Zoom Out Reset image size

This suggests that the SFR in our clouds is determined primarily by the initial conditions, via the density structure that is imposed by the interplay between gravity and turbulence. The majority of stellar mass growth is through stars formed early on before radiative feedback becomes important. As will be discussed further in Section 4.6, radiation does little to suppress this early star formation. Instead, radiative feedback is primarily important in rapidly truncating star formation by driving gas from the cloud, once a sufficient number of stars have formed in the turbulent density field.

4.2.1. Turbulent Outflows

The most obvious initial effect of turbulence is the unbinding of small mass fractions with high velocities in the tail of the distribution. We may calculate the total mass in outflows due to the initial turbulence alone by running simulations with no radiative feedback until all of the gas in the box is converted to stars. In these models, any gas that leaves the cloud is unbound by the initial turbulence.

In Figure 17(a), we show the corresponding outflow mass fractions ${\varepsilon }_{{\rm{of}}}$ for our no-feedback Σ-series models. The outflowing mass amounts to between 10% and 15% of the initial cloud mass, with little variation between low and high surface density. Without radiation, the primary dimensionless parameter that varies for different values of R and M at fixed ${\alpha }_{{\rm{vir,0}}}=2.0$ is the Mach number ${ \mathcal M }$. However, there is little direct dependence of the outflow mass on ${ \mathcal M }$ at fixed ${\alpha }_{{\rm{vir,0}}}$, since the escape velocity from the edge of the cloud increases proportionally to the Mach number at fixed virial parameter.

Zoom In Zoom Out Reset image size

As shown in Figure 17(b), the outflow mass depends strongly only on the initial virial parameter ${\alpha }_{{\rm{vir,0}}}\propto {{ \mathcal M }}^{2}/({GM}/r)$, i.e., the strength of the turbulence relative to gravity. Below ${\alpha }_{{\rm{vir,0}}}\sim 1$, there is effectively no outflow, as seen for the case of ${\alpha }_{{\rm{vir,0}}}=0.1$. Meanwhile, above ${\alpha }_{{\rm{vir,0}}}\sim 5$, the outflow mass fraction approaches unity.

The dependence on ${\alpha }_{{\rm{vir,0}}}$ can be understood in terms of the fraction of gas that escapes due to an initial turbulent velocity above the local escape velocity. For a uniform density cloud, the escape velocity at radius r is ${v}_{e}^{2}(r)=({GM}/{r}_{0})(3-{r}^{2}/{r}_{0}^{2})$, while in our models the initial velocity field is Gaussian with a dispersion of ${v}_{\mathrm{rms}}^{2}=3{GM}{\alpha }_{{\rm{vir,0}}}/(5{r}_{0})$. Therefore, the total outflow mass fraction due to the initial turbulence is just the fraction of mass, at any given radius, above v_e(r):

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{of,init}}}={\int }_{0}^{1}{dx}\,3{x}^{2}\left[1-{\rm{erf}}\sqrt{\displaystyle \frac{5(3-{x}^{2})}{6{\alpha }_{{\rm{vir}}}}}\right].\end{eqnarray} \tag{ 13 }$$

This quantity is indicated by the black curve in Figure 17(b). For ${\alpha }_{{\rm{vir,0}}}\leqslant 3$, Equation (13) somewhat overestimates the total outflow mass, since the initial turbulence is damped, and gas in the cloud interior, which may have an initial velocity higher than the escape velocity, will nevertheless collapse to form stars. However, it does represent a reasonable estimate, since the majority of outflowing mass is on the outer edges of the cloud, where the escape velocity is lowest.

Above ${\alpha }_{{\rm{vir,0}}}\sim 3$, we no longer predict the outflow mass well since we cannot distinguish between gas that leaves the box and gas that becomes truly unbound. Even accounting for an escape velocity to 2r₀ rather than infinity, ${v}_{2}{(r)}^{2}=({GM}/{r}_{0})(2-{r}^{2}/{r}_{0}^{2})$, indicated by the red line in Figure 17(b), we still underestimate the outflow mass at large ${\alpha }_{{\rm{vir,0}}}$, since the remaining gas mass is at too low a density to form a significant mass in stars, and so may drift outside the box while still bound. However, without dramatically increasing the box size, with associated high computational cost, this degeneracy is unavoidable. Therefore, at values of ${\alpha }_{{\rm{vir,0}}}\gtrsim 3$, Equation (13) may underestimate the net SFE if, in reality, some gas that is expelled from our box were to ultimately recollapse.

4.2.2. Surface Density PDF

In addition to driving outflows, turbulence also has an effect on the gas density structure. In particular, the non-uniformity of the cloud means that it has significant variations in surface density. To quantify this, we calculate the surface densities projected in the x–y, x–z, and y–z planes, and we consider the probability distribution function over all three.

For our fiducial ${\alpha }_{{\rm{vir,0}}}=2$ model, the resultant surface density distributions by both area ${P}_{A}({\rm{\Sigma }})$ and mass ${P}_{M}({\rm{\Sigma }})$ are shown in Figures 18(a) and (c) at t₁₀. By this time, both P_A and P_M are roughly log-normal in shape at high Σ, with similar variances of ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}\sim 1$. The area distribution has a long tail at low surface density due to the fact that we are sampling over the whole simulation volume rather than just the cloud volume. Meanwhile, as there is little cumulative mass in these low-density regions, the mass distribution is much more obviously log-normal, consistent with previous simulations of GMCs with supersonic turbulence (see e.g., Ostriker et al. 2001; Vázquez-Semadeni & García 2001; Federrath et al. 2010).

Zoom In Zoom Out Reset image size

We may therefore parameterize the "cloud" portion of both ${P}_{A}({\rm{\Sigma }})$ and ${P}_{M}({\rm{\Sigma }})$ as log-normal distributions with mean ${\mu }_{(A/M)}$ and standard deviation ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$:

$$\begin{eqnarray}&&{P}_{(A/M)}({\rm{\Sigma }})\,d\mathrm{ln}{\rm{\Sigma }}=\displaystyle \frac{1}{{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}\sqrt{2\pi }}\times \\ &&{\rm{\exp }}\left[-\displaystyle \frac{{({\rm{ln}}{\rm{\Sigma }}-{\mu }_{(A/M)})}^{2}}{2{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2}}\right]\,d\mathrm{ln}{\rm{\Sigma }}.\end{eqnarray} \tag{ 14 }$$

At any given time, the total mass in the system is conserved so that ${P}_{M}({\rm{\Sigma }})\propto {\rm{\Sigma }}{P}_{A}({\rm{\Sigma }})$. This implies that both distributions must have the same standard deviation and gives a relation between the mean surface density and the means of the two log-normals:

$$\begin{eqnarray}&&{\rm{ln}}\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}={\mu }_{A}+\displaystyle \frac{1}{2}{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2}={\mu }_{M}-\displaystyle \frac{1}{2}{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2}.\end{eqnarray} \tag{ 15 }$$

Thus, half of the cloud's mass is at surface densities above/below $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}\exp ({\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2}/2)$, and half of its area is at surface densities above/below $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}\exp (-{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2}/2)$. Here, we emphasize the difference between the mean cloud surface density $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}\sim {M}_{{\rm{cl}}}/{A}_{{\rm{cl}}}$, which can be obtained using Equation (15) after fitting log-normals to obtain ${\mu }_{(A/M)}$ and ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$, and the mean surface density over the whole simulation volume $\langle {\rm{\Sigma }}{\rangle }_{{\rm{box}}}\sim {M}_{{\rm{cl}}}/{A}_{{\rm{box}}}$.

Comparing the fitted log-normal functional form to the PDFs in Figure 18 we see that the log-normal approximation is reasonable, particularly at the high-density end for P_A and around the peak for P_M. The fits we show are only calculated between the $10{\rm{th}}$ and 90th percentile in mass, so as to avoid sampling the low surface density regions outside the cloud. We also assume measurement errors proportional to $\sqrt{N}$ where N is the number of grid cells in each density bin. At the time shown in Figure 18(c), the median or peak surface density of the mass distribution for the fiducial model is around ${\rm{\Sigma }}\sim 100\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$, but part of the mass lies beyond even ${\rm{\Sigma }}\sim 300\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$. The peak of the distribution by area shown in Figure 18(a) is lower, ${\rm{\Sigma }}\sim 25\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$. As we shall discuss below, the wide variation in Σ implies that at early times, radiation forces can be much more effective in some regions than in others.

The distribution is not just a broad log-normal at the onset of star formation, but remains broad as it progresses. For example, Figures 18(b) and (d) show the area and mass PDFs at ${t}_{50}=1.06{t}_{{\rm{ff,0}}}$. At this time, the best-fit parameters are ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}=1.38$ (where this is taken as the average of fits to the area and mass distributions), ${\mu }_{A}=2.53$, and ${\mu }_{M}=4.45$, so that the mass conservation relations are approximately followed and $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}\sim 33\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$.

In Figures 19(a) and (b), we show the best-fit mean and standard deviation parameters as a function of time for log-normal fits to the surface density distribution of our fiducial cloud. A drop in reduced ${\chi }^{2}$ values in Figure 19(c) only really develops after $\sim 0.4{t}_{{\rm{ff,0}}}$, indicating that it is only by this time that turbulence has erased the initial conditions. This is expected, since the cloud is initially uniform and the development of density structure has a timescale set by the turbulent crossing time ${t}_{{\rm{cross}}}={r}_{0}/{v}_{{\rm{rms}}}=1.2{t}_{{\rm{ff,0}}}$. While the surface density distribution gradually broadens over the first ${t}_{{\rm{ff,0}}}$, after this time it effectively reaches a steady state in both the feedback and no-feedback models, with ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}\sim 1\,-\,1.5$. Even though the surface density distributions are visibly different at $t=1.06{t}_{{\rm{ff,0}}}$ for the feedback (Figure 3) and no-feedback (Figure 15) models, they are statistically similar, since it is only the lower column density regions that are affected by radiative feedback. This holds true until around $1.5{t}_{{\rm{ff,0}}}$, by which time the majority of stars have formed and the majority of remaining gas in the fiducial model is outflowing. As gas is accreted by star particles or is driven from the cloud, the mean and peak of the distribution naturally shift downward. However, the width of the distribution remains similar throughout star formation.

Zoom In Zoom Out Reset image size

As noted above, we cannot obtain $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}$ directly from our simulations, since our simulation volume comprises both the "cloud" and the surrounding empty region (or larger-scale lower-density ISM in a real system). However, if we assume that the "cloud" portion of the distribution is a log-normal, then we can fit to the PDFs and use Equation (15) to obtain $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}$. Our fits to obtain ${\mu }_{A}$, ${\mu }_{M}$, and ${\sigma }_{{\rm{ln}}{\rm{\Sigma }}}$ extend from the $10{\rm{th}}$ to the 90th percentile by mass, which effectively excises the low surface density regions external to the cloud.

We can define the mass-weighted mean surface density over the whole box as $\langle {\rm{\Sigma }}{\rangle }_{M,\mathrm{box}}=\langle {{\rm{\Sigma }}}^{2}{\rangle }_{{\rm{box}}}/\langle {\rm{\Sigma }}{\rangle }_{{\rm{box}}}$, where $\langle {\rm{\Sigma }}{\rangle }_{{\rm{box}}}$ is the area-weighted mean. Under the assumption that the "cloud" portion is a log-normal, $\langle {\rm{\Sigma }}{\rangle }_{M,\mathrm{cloud}}=\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}\exp ({\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2})$, where $\langle {\rm{\Sigma }}{\rangle }_{A,\mathrm{cloud}}\equiv \langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}$. Figure 20 shows, for the fiducial model, the evolution of the area- and mass-weighted surface densities computed in two different ways: taking direct averages over the box, and using the fitted log-normals to identify just the "cloud" material.

Zoom In Zoom Out Reset image size

From Figure 20, $\langle {\rm{\Sigma }}{\rangle }_{{\rm{box}}}$ starts a factor ∼5 below $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}$, because the projected surface area of the cloud is $\sim 1/5$ of the box surface area. Over time, as the cloud disperses and fills more of the simulation area and volume, these measures become more similar. In contrast, $\langle {\rm{\Sigma }}{\rangle }_{M,\mathrm{box}}$ and $\langle {\rm{\Sigma }}{\rangle }_{M,\mathrm{cloud}}$ are much closer to each other over the whole evolution, and are each an order of magnitude larger than $\langle {\rm{\Sigma }}{\rangle }_{{\rm{box}}}$. We conclude that both the definitions $\langle {\rm{\Sigma }}{\rangle }_{A,\mathrm{cloud}}\equiv {\rm{\exp }}({\mu }_{M}-(1/2){\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2})$ and $\langle {\rm{\Sigma }}{\rangle }_{M,\mathrm{cloud}}\equiv {\rm{\exp }}({\mu }_{M}+(1/2){\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2})$ well represent the area- and mass-weighted mean values for the cloud material, and that the latter is also similar to the mass-weighted mean over the whole box.

We note also that there is no significant change in the width of the surface density distribution with differing initial cloud mass and radius. For a number of different models in our Σ-series, we compute best-fit values of ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$ at times ${t}_{10},{t}_{50},$ and t₉₀, with the results shown in Figure 21(a). In most models, there appears to be a modest increase in the log-normal width with time from the width at t₁₀, since this occurs at around $\sim 0.6{t}_{{\rm{ff,0}}}$, and as discussed earlier, the log-normal has not quite reached a steady state by this stage. However, after $t\sim 1\mbox{--}1.2{t}_{{\rm{ff,0}}}$, the distribution widths remain roughly constant for the remainder of star formation. Moreover, the distribution width is relatively independent of cloud surface density at ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}\sim 1.3\mbox{--}1.5$ (slightly decreasing towards the low-Σ end, which has lower Mach number).

Zoom In Zoom Out Reset image size

This makes it clear why the simple hypothesis that stars will form until $\langle {{\rm{\Sigma }}}^{c}{\rangle }_{{\rm{cloud}}}={{\rm{\Sigma }}}_{E}$ is seriously flawed: even the peak by surface density's mass, $\langle {\rm{\Sigma }}{\rangle }_{M,\mathrm{cloud}}$, is a factor of ${\rm{\exp }}({\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2}/2)\sim 3.1$ higher than $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}$. Thus, half of the gas is at surface density more than three times $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}$, and a large fraction is in regions at even higher surface density. Forcing these high-Σ regions out of the cloud would demand ${{\rm{\Sigma }}}_{E}\gg \langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}$, requiring a higher luminosity and hence (from Equation (10)) larger ε than predicted by Equation (12).

In the above, we have analyzed the distribution of surface densities as would be seen by an external observer. However, for the purposes of gauging the effects of radiation forces from a cluster on the surrounding gas in a cloud, what matters is actually the circumcluster distribution of densities. As we shall show in a separate publication (Raskutti et al. 2016; Paper II), this distribution is in fact quite similar to the log-normal PDFs shown and discussed above. For example, for the fiducial model, we find $\langle {{\rm{\Sigma }}}^{c}{\rangle }_{{\rm{cloud}}}\sim 12\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ and ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}=1.42$ at $t=1.06\,{t}_{{\rm{ff,0}}}$, which can be compared to the variance and mean values seen by an external observer ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}=1.38$ and $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}\sim 33\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ (see also Figures 19 and 20).

The mean value of the cloud surface density at any time is related to its initial mean surface density via the efficiency and a radial contraction factor x (cf. Equation (11)), albeit adjusted for the initial turbulent outflows

$$\begin{eqnarray}&&\langle {{\rm{\Sigma }}}^{c}\rangle =\displaystyle \frac{(1-\varepsilon )}{4{x}^{2}}{{\rm{\Sigma }}}_{{\rm{adj}}},\end{eqnarray} \tag{ 16 }$$

where ${{\rm{\Sigma }}}_{{\rm{adj}}}={{\rm{\Sigma }}}_{{\rm{cl,0}}}(1.0-{\varepsilon }_{{\rm{of,init}}})$. In addition, fitting a log-normal form gives a value of $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}$ from Equation (15). Putting these relations together yields x. We shall show in Paper II that for models with ${\alpha }_{{\rm{vir,0}}}\sim 1\mbox{--}2$, this yields $x\approx 1$. That is, the overall size of clouds does not vary much over the star-forming period.

4.2.3. Dependence of the PDF on Seed Field

Before we turn to how the surface density field affects the final stellar efficiency, we briefly consider what factors are important in determining its shape. To do this, we run a variation of the fiducial model in which we initially allow our cloud—with the same seed turbulent velocity field—to evolve without self-gravity. Without star formation and feedback, the gas flows steadily out from the simulation volume due to the initial turbulence, while the Mach number decreases as the turbulence decays and the highest velocity material leaves the box.

However, after half a free-fall time, when stars are beginning to form in the fiducial model, the density distribution has reached a rough steady state (as shown from the variance of the column density PDF in Figure 22). From this point forward, the column density distribution remains roughly log-normal, with a fairly steady, though slowly declining width, and a steadily declining mean as the gas flows out of the box. It is also relatively similar to the fiducial model with gravity, suggesting that the initial density field is set almost entirely by the initial turbulence, with gravity only being important on the smallest scales (or highest densities). The only differences are a reduced width, due to the absence of self-gravity, which flattens the distribution at the high-density end, and a reduced mean density due to the increased gas outflows. At the onset of star formation, therefore, the bulk of gas mass (i.e., up to the highest densities) has a log-normal shape set almost entirely by the initial turbulent field.

Zoom In Zoom Out Reset image size

Figure 21(b) shows the dependence of ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$ on the initial virial parameter. Before the majority of star formation begins at t₁₀, the lower ${\alpha }_{{\rm{vir,0}}}$ models have slightly narrower distributions (lower ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$) due to their lower Mach number. However, this trend is reversed over time, with the lower virial parameter models tending to have much broader distributions by the end of star formation. Largely this is because these clouds have very high efficiencies, so that only a small fraction of the cloud mass remains once 90% of stars have formed, hence the log-normal fit is considerably worse.

It must also be kept in mind that ${\alpha }_{{\rm{vir}}}$ changes in time for most of the α series. At high initial virial parameter, ${\alpha }_{{\rm{vir,0}}}\gtrsim 5$, the initial turbulence rapidly drives the highest velocity regions from the box, so that the virial parameter rapidly decays. Similarly, at low initial ${\alpha }_{{\rm{vir,0}}}$, including the reduced-turbulence models considered in Section 4.1, the cloud contracts until the virial parameter is close to unity in all cases. This is evident if we consider the variation of ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$ with instantaneous rather than initial virial parameter, shown in Figure 23(b). Clouds in the range ${\alpha }_{{\rm{vir,0}}}=0.1\mbox{--}3.0$ all converge to a much smaller range of ${\alpha }_{{\rm{vir}}}=0.5\mbox{--}1.0$ by the time star formation begins. Then, as radiative feedback becomes important, ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$ increases slightly with increasing virial parameter.

Zoom In Zoom Out Reset image size

In this subsection, we present and analyze results for the full set of turbulent cloud collapse models with radiation feedback. In Figure 24 we show the final star formation efficiencies for our Σ- and α-series cloud models; all values are also listed in Table 3. We have normalized by both the initial cloud mass and the cloud mass corrected for turbulence-driven outflows. Interestingly, the efficiency appears to show a logarithmic dependence on Σ across almost two orders of magnitude of variation in the initial surface density. In fact, Figure 24 shows a remarkably good fit to the relation $\varepsilon =0.37{\rm{log}}{\rm{\Sigma }}-0.26$ (black line). This is true even if we account for the mass loss due to turbulent outflows from the simulation box by setting $\varepsilon \to {\varepsilon }_{{\rm{adj}}}\equiv \varepsilon /(1.0-{\varepsilon }_{{\rm{of,init}}})$. The normalization of the relation changes slightly, but we still obtain a logarithmic relationship between surface density and final efficiency, given by $\varepsilon =0.41{\rm{log}}{\rm{\Sigma }}-0.26$ (red line).

Zoom In Zoom Out Reset image size

We can also fit to a logarithmic dependence on the virial parameter, with $\varepsilon =-0.45{\rm{log}}\alpha +0.51$, although in this case it is less clear whether this is significant, as there are systematic errors at both low and high virial parameters. Moreover, our estimate of the net efficiency is uncertain at high virial parameter, since we overestimate turbulent outflows and correspondingly underestimate star formation.

Although the logarithmic form fits the Σ series well, it is quite different from simple predictions. For example, we can compare the numerical results to the uniform shell force balance model (with x = 1) of Equation (12), as shown with the red solid line in Figure 25. Evidently, both the shape and magnitude of the curve from Equation (12) compare poorly with the simulation results.

Zoom In Zoom Out Reset image size

As discussed earlier, the key difference between the simulated clouds and the simple Eddington limit prediction of Equation (12) is that the cloud does not have uniform surface density. Stars do not keep forming inside a uniform shell and then instantaneously drive that shell away once the efficiency ε and luminosity are large enough for ${{\rm{\Sigma }}}_{E}$ to match $\langle {{\rm{\Sigma }}}^{c}{\rangle }_{{\rm{cloud}}}$. Instead, as will be discussed in more detail in Raskutti et al. (2016), early star formation begins to drive away the lowest surface density regions (those that have have ${{\rm{\Sigma }}}^{c}\lt {{\rm{\Sigma }}}_{E}$). Stars continue forming from the remaining mass until the radiative force is sufficient to drive away much higher surface density gas (as, from Equation (10), an increase in ε raises ${{\rm{\Sigma }}}_{E}$). For example, in our fiducial model, close to the end of star formation at $t=1.45{t}_{{\rm{ff,0}}}$ when $\varepsilon \sim 0.35$, the mean surface density is $\langle {\rm{\Sigma }}{\rangle }_{{\rm{cloud}}}\sim 12\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$, with the circumcluster surface density around a factor of two or three lower again, while the Eddington surface density in Equation (10) is ${{\rm{\Sigma }}}_{E}\sim 200\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$, implying that enough stars have formed to drive away gas at close to 50 times the mean surface density of the cloud. Figure 19 shows that at this time ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}=1.49$, which yields ${\mu }_{M}=2.9$ such that ${\rm{ln}}{{\rm{\Sigma }}}_{E}=5.3$ is roughly 1.5-σ above the peak surface density. This would suggest that star formation is halted when enough stars form to drive away not the mean surface density, but instead something closer to the 90th percentile of surface density.

Recently, Thompson & Krumholz (2016) have argued that understanding how radiative feedback limits star formation in GMCs requires an accounting of the full surface density distribution set by turbulence. They propose a model in which the instantaneous mass-loss rate from a cloud is set by the cloud's free-fall time and fraction of mass in a log-normal PDF at surface densities below a critical value. At the same time, they assume the stellar mass and luminosity increases as $\dot{{M}_{* }}={\varepsilon }_{{\rm{ff}}}{M}_{{\rm{gas}}}/{t}_{{\rm{ff}}}$. Here, we develop a related model by considering the mass eligible to be expelled from the cloud at any time, for a given log-normal PDF and SFE.

From Equation (10), surface densities above the Eddington surface density ${{\rm{\Sigma }}}_{E}$ have a net inward force and can continue to either collapse and form stars, or be accreted onto the growing star clusters. Correspondingly, surface densities below ${{\rm{\Sigma }}}_{E}$ have a net outward force, and we assume here that such regions become instantaneously unbound without further mixing with other gas. Although this is clearly an oversimplification, we have found (see Raskutti et al. 2016 for details) that the distribution of outflowing velocities is consistent with essentially ballistic outflow of super-Eddington structures.

If the cloud has a circumcluster surface density distribution ${P}_{M}({{\rm{\Sigma }}}^{c})$ by mass, the fraction of mass eligible to be expelled will be ${\int }_{-\infty }^{\mathrm{ln}{{\rm{\Sigma }}}_{E}}{P}_{M}({{\rm{\Sigma }}}^{c})d\mathrm{ln}{{\rm{\Sigma }}}^{c}$. If the stellar mass in the cloud at a given instant during its evolution is ${M}_{* }=\varepsilon {M}_{{\rm{cl,0}}}$, then the gas mass remaining is $(1-\varepsilon ){M}_{{\rm{cl,0}}}$, and the fraction ${\varepsilon }_{{\rm{of}}}$ of the original gas that is eligible for outflow is

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{of}}}=(1-\varepsilon ){\int }_{-\infty }^{\mathrm{ln}{{\rm{\Sigma }}}_{E}}{P}_{M}({{\rm{\Sigma }}}^{c})d\mathrm{ln}{{\rm{\Sigma }}}^{c}.\end{eqnarray} \tag{ 17 }$$

In Section 4.2.2, we showed that the column density distribution is well approximated by a log-normal distribution, and here, we extend that approximation to the circumcluster surface density (see Raskutti et al. 2016 for details). In this case, Equation (17) may be evaluated as

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{of}}}=\displaystyle \frac{1}{2}(1-\varepsilon )(1+{\rm{erf}}({y}_{E})),\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}&&{y}_{E}\equiv \displaystyle \frac{{\rm{ln}}{{\rm{\Sigma }}}_{E}-{\mu }_{M}}{\sqrt{2}{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}}.\end{eqnarray} \tag{ 19 }$$

The quantity ${\mu }_{M}=\mathrm{ln}\langle {{\rm{\Sigma }}}^{c}\rangle +{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2}/2$ (see Equation (15)) is the mean of the mass distribution of ln${{\rm{\Sigma }}}^{c}$ and ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$ is the variance. Using Equation (11),

$$\begin{eqnarray}&&{\mu }_{M}={\rm{ln}}[{{\rm{\Sigma }}}_{{\rm{cl,0}}}(1-\varepsilon )]+\displaystyle \frac{{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2}}{2}-{\rm{ln}}(4{x}^{2})\end{eqnarray} \tag{ 20 }$$

if the circumcluster gas is concentrated in a thin shell. Substituting into Equation (19) and using Equation (10), we obtain

$$\begin{eqnarray}{y}_{E} & = & \displaystyle \frac{1}{\sqrt{2}{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}}{\rm{ln}}\left(\displaystyle \frac{4{{\rm{\Sigma }}}_{E}{x}^{2}}{{{\rm{\Sigma }}}_{{\rm{cl,0}}}(1-\varepsilon )}\right)-\displaystyle \frac{{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}}{\sqrt{8}}\\ & = & \displaystyle \frac{1}{\sqrt{2}{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}}\left[{\rm{ln}}\left[\displaystyle \frac{2{\rm{\Psi }}}{\pi {cG}{{\rm{\Sigma }}}_{{\rm{cl,0}}}}\right]\right.\\ & & \left.+2{\rm{ln}}x-\displaystyle \frac{{\sigma }_{\mathrm{ln}{\rm{\Sigma }}}^{2}}{2}+{\rm{ln}}\displaystyle \frac{\varepsilon }{1-{\varepsilon }^{2}}\right].\end{eqnarray} \tag{ 21 }$$

With Equation (21), Equation (18) gives the unbound or outflowing fraction ${\varepsilon }_{{\rm{of}}}$ of a cloud with initial surface density ${{\rm{\Sigma }}}_{{\rm{cl,0}}}\equiv {M}_{{\rm{cl,0}}}/(\pi {r}_{{\rm{cl,0}}}^{2})$ in terms of the current efficiency ε and the two free parameters x and ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$.

Figure 26 shows the behavior of ${\varepsilon }_{{\rm{of}}}$ from Equation (18) as a function of input stellar efficiency ε for varying ${{\rm{\Sigma }}}_{{\rm{cl,0}}}$ and ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$. When $\varepsilon =0$, ${\varepsilon }_{{\rm{of}}}=0$ since with no stars, the Eddington surface density is zero. As the mass in stars increases, ${\varepsilon }_{{\rm{of}}}$ initially also increases under two competing influences. More stars increase the radiative force and ${{\rm{\Sigma }}}_{E}$ so that a larger fraction of the cloud is eligible to become unbound; this increases the factor $1+{\rm{erf}}({y}_{E})$ in Equation (18). However, a larger stellar mass also decreases the gas mass since a larger fraction of the cloud is already bound up in stars; this decreases the factor $1-\varepsilon$. In the limit $\varepsilon \to 1$, there is no gas reservoir and therefore ${\varepsilon }_{{\rm{of}}}\to 0$. For some intermediate value of ε between 0 and 1, an infinitesimal increase in ε would decrease ${\varepsilon }_{{\rm{of}}}$ in Equation (18), since the overall decrease in the gas mass reservoir (lower $1-\varepsilon$) outweighs the increase in the fraction of gas that is super-Eddington (higher y_E). This point represents the maximum value of ${\varepsilon }_{{\rm{of}}}$ for any given set of initial cloud parameters and variance in the PDF.

Zoom In Zoom Out Reset image size

We consider a cloud such that, for a given PDF variance, the maximum in ${\varepsilon }_{{\rm{of}}}$ has been reached; its efficiency is then $\varepsilon =\text{arg max}{\varepsilon }_{{\rm{of}}}$. At this point, the outflowing gas mass cannot decrease because this material has already become super-Eddington. The remaining gas reservoir that is neither stars nor outflowing gas is a fraction $(1-{\rm{\max }}\,{\rm{}}\,{\varepsilon }_{{\rm{of}}}-\text{arg max}{\varepsilon }_{{\rm{of}}})$ of the original cloud. If this material collapses faster than the PDF in larger-scale cloud can adjust, it will all be added to the existing stars and the final SFE will be ${\varepsilon }_{{\rm{final}}}=1-{\rm{\max }}\,{\rm{}}\,{\varepsilon }_{{\rm{of}}}$. The luminosity from these additional stars would also increase the radiation pressure on the outflowing gas. Alternatively, if the collapse is slower, there may be time for the log-normal density distribution of the remaining gas to adjust, and the final SFE may rise to a level between $\text{arg max}{\varepsilon }_{{\rm{of}}}$ and $1-{\rm{\max }}\,{\rm{}}\,{\varepsilon }_{{\rm{of}}}$. This suggests that for a given PDF variance, there are upper and lower bounds on ${\varepsilon }_{{\rm{final}}}$:

$$\begin{eqnarray}{\varepsilon }_{{\rm{\max }}} & = & (1-\mathop{\max }\limits_{0\lt \varepsilon \lt 1}{\varepsilon }_{{\rm{of}}})\\ & = & \displaystyle \frac{1}{2}\mathop{\min }\limits_{0\lt \varepsilon \lt 1}[1+\varepsilon +(\varepsilon -1){\rm{erf}}({y}_{E})]\end{eqnarray} \tag{ 22a }$$

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{\min }}}={\rm{\arg }}\mathop{\max }\limits_{0\lt \varepsilon \lt 1}{\varepsilon }_{{\rm{of}}}.\,\end{eqnarray} \tag{ 22b }$$

From Figure 26, the maximum of the function ${\varepsilon }_{{\rm{of}}}$ increases with decreasing ${{\rm{\Sigma }}}_{{\rm{cl,0}}}$, since at lower surface densities the gas reservoir can be driven away more easily. Also, broader surface density distributions (larger ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$) tend to decrease ${\rm{\max }}\,{\varepsilon }_{{\rm{of}}}$ since more gas is at the highest density, which is more difficult to unbind.

In Figure 25, we compare the predictions of Equation (22) to our numerical results for net star formation efficiencies as a function of initial surface density. We note that we are primarily interested in the efficiency relative to the cloud mass adjusted for initial turbulent outflows ${\varepsilon }_{{\rm{adj}}}=\varepsilon /(1-{\varepsilon }_{{\rm{of,init}}})$, so we substitute the adjusted surface density ${{\rm{\Sigma }}}_{{\rm{cl,0}}}\to {{\rm{\Sigma }}}_{{\rm{adj}}}={{\rm{\Sigma }}}_{{\rm{cl,0}}}(1-{\varepsilon }_{{\rm{of,init}}})$ in Equation (21). For the parameters entering Equation (21), we use x = 0.84 and ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}=1.42$, based on the time average of the best-fit values for the circumcluster surface density in the fiducial model (see Paper II for details). We use these same values of x and ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$ at all ${{\rm{\Sigma }}}_{{\rm{cl,0}}}$. Additionally, since we are interested in comparing to the physical observed cloud, rather than the artificial initial conditions, we show the adjusted efficiency as a function of the adjusted cloud surface density $\langle {{\rm{\Sigma }}}_{{\rm{cloud}}}\rangle ={{\rm{\Sigma }}}_{{\rm{adj}}}/{x}^{2}$ rather than ${{\rm{\Sigma }}}_{{\rm{cl,0}}}$.

Figure 25 shows that ${\varepsilon }_{{\rm{\max }}}$ from Equation (22a) represents a reasonable estimate of the actual SFE found in the simulations, both in normalization and in the shape of the dependence on ${{\rm{\Sigma }}}_{{\rm{cloud}}}$. In principle, however, the final efficiency might be closer to ${\varepsilon }_{{\rm{\min }}}$ if conditions were such that the SFR were lower.

Equation (22) predicts much higher efficiencies than the fiducial model of Thompson & Krumholz (2016), shown for comparison in Figure 25. Their model has similar ingredients to ours (a log-normal surface density distribution, with both sub-Eddington and super-Eddington regions), with the principal difference being that the Thompson & Krumholz (2016) formalism assumes a fixed SFR per unit gas mass, whereas there is no assumption about the SFR in our model. Specifically, they assume that stars form at a rate $\dot{{M}_{* }}\,={\varepsilon }_{{\rm{ff}}}{M}_{g}/{t}_{{\rm{ff}}}$ from the total remaining mass (including super-Eddington ${\rm{\Sigma }}\lt {{\rm{\Sigma }}}_{E}$ regions), and that winds are driven out at a rate proportional to $1/{t}_{{\rm{ff}}}$ from the super-Eddington gas. It is their adoption of a very small fiducial value ${\varepsilon }_{{\rm{ff}}}=0.01$ (implying a vast discrepancy between star formation and wind expulsion rates) that leads to a highly suppressed final efficiency in their fiducial model. If instead we adopt ${\varepsilon }_{{\rm{ff}}}=0.44$ (similar to our simulation results in Section 4.6) and apply their formula, their model prediction is somewhat closer to ours, albeit with a lower normalization and shallower dependence on ${{\rm{\Sigma }}}_{{\rm{cloud}}}$ (see Figure 25).

The correspondence between the prediction of ${\varepsilon }_{{\rm{\max }}}$ in Equation (22a) and our Σ-series numerical model results is close enough to suggest that the dominant effect in suppressing star formation is radiative feedback driving out structures at successively higher surface densities until a maximum mass of outflowing material is reached. However, there are a number of issues, or at least questions, surrounding this model. First among these is whether the correlation timescale of the log-normal surface density distribution in the cloud is long enough to allow persistent acceleration by the central stars. This is because the dynamical evolution of any given fluid element depends on the coherence time of the (Lagrangian) evolution for the surface density region surrounding it. In principle, the surface density distribution could remain statistically log-normal at all times, while individual regions fluctuate rapidly. If these fluctuations in time are much shorter than the time taken for radiation to accelerate gas from the cloud, then any fluid element would fully sample the distribution of surface densities, and only the mean cloud surface density would be relevant.

The question of the column density correlation timescales in comparison to the cloud destruction timescale has already been discussed in a heuristic manner in Thompson & Krumholz (2016). They compare the turbulence crossing time to the acceleration timescale ${t}_{{\rm{acc}}}\propto {r}_{0}/{v}_{{\rm{esc}}}({r}_{0})$ and argue that so long as the radiation force is several times stronger than the force of gravity, densities will fluctuate on longer timescales than it takes for the cloud to unbind.

In our simulations, we can measure the temporal correlations of the column in a given area of the sky. We find a correlation time $\sim 0.5{t}_{{\rm{ff,0}}}$ once star formation has begun. This is roughly comparable to the timescales on which gas is accelerated out of the cloud, suggesting that outflowing low surface density regions might merge with collapsing higher surface density regions before they have a chance to escape the cloud. However, this overall Eulerian correlation time is not necessarily representative for low-density regions. In addition, we do find an outflowing velocity distribution consistent with low-density regions remaining correlated until they escape the cloud (see Raskutti et al. 2016 for details) and so conclude that this interpretation is not unreasonable. In future tests, to answer this question realistically, we would need to use tracer particles in the gas to track the flow of individual fluid elements.

We note also that the values of ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$ from our simulations (see Figures 21 and 23) are somewhat larger than current estimates from observations, which typically find ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}\lt 1$. This may in part be due to line-of-sight contamination, which tends to reduce the observed ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$ (e.g., Schneider et al. 2015), and in part to the absence of magnetic fields in the present models, as magnetization reduces shock compression and therefore density variance (e.g., Ostriker et al. 2001; Molina et al. 2012). Figure 27(a) shows the prediction of Equation (22a) for ${\varepsilon }_{\mathrm{final}}$ as a function of ${{\rm{\Sigma }}}_{{\rm{cl}}}$ for a range of ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$, demonstrating that the predicted net SFE in a cloud could be considerably lower at low ${\sigma }_{\mathrm{ln}{\rm{\Sigma }}}$.

Zoom In Zoom Out Reset image size

As a final point, it is worth considering the effects on our model of a varying Ψ due to undersampling of the IMF. Kim et al. (2016) recently studied the effects of an undersampled Chabrier (2003) IMF on the value of Ψ, both stochastically and on average, by using the SLUG code (Krumholz et al. 2015) to simulate the photon output from stellar clusters as a function of mass. They found that although $100\,{M}_{\odot }$ clusters produce significantly fewer photons, with ${\rm{\Psi }}\sim 200\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$, by the time ${M}_{\star }\sim 2\times {10}^{3}\,{M}_{\odot }$, the cluster luminosity per unit mass has settled to the final fully sampled value of ${\rm{\Psi }}\sim 2000\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$, with a variation of around ±0.3 dex around this median. For even larger clusters of ${M}_{\star }\sim {10}^{4}\,{M}_{\odot }$, this variance drops to only ±0.1 dex. We may therefore expect the stellar IMF to be fully sampled for ${M}_{\star }\gtrsim 2\times {10}^{3}\,{M}_{\odot }$.

In all of our simulations, the stellar mass exceeds $2\times {10}^{3}\,{M}_{\odot }$ by the end of the star formation epoch. Our lowest mass cloud has a mass of only $5\times {10}^{3}\,{M}_{\odot }$ but an efficiency of more than 40%, so that more than $2\times {10}^{3}\,{M}_{\odot }$ is in stars. When star particles first form, realistically the IMF would be undersampled. By setting Ψ to a constant, we overestimate the radiative force at early times, while correspondingly underestimating the SFR. However, the final efficiencies are largely set by the balance between the radiative force and the gas self-gravity at late times, so these are not changed significantly by our overestimation of Ψ at early times.

The smallest of our star particles are less than ${M}_{\star }\sim {10}^{2}\,{M}_{\odot }$, which would imply correspondingly low median values of ${\rm{\Psi }}\sim 200\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$. Therefore, when these star particles first form, we are overestimating the radiative force by a factor of almost 10. We may bracket the effect of this by simulating our fiducial model over a range of possible Ψ values. In Figure 28 we show the stellar efficiency for three different models with Ψ varying between 200 and $2000\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$. The models with lower Ψ form stars at a slightly faster rate and end up with a higher efficiency, since less gas is driven from the cloud by the effects of UV radiation.

Zoom In Zoom Out Reset image size

However, we note that the variation in efficiency is only ∼50% despite an order-of-magnitude variation in the luminosity per unit mass. In simple models of radiative feedback that assume a uniform shell of gas driven away by radiative pressure, $1/\varepsilon -\varepsilon$ is linear with Ψ. Therefore, variation by an order of magnitude in Ψ would correspond to almost an order of magnitude variation in ε, or at least a saturation at unity. By contrast, the model we develop, which takes into account the log-normal surface density distribution with best-fit values of x and ${\sigma }_{{\rm{ln}}{\rm{\Sigma }}}$, does a much better job of matching the final stellar efficiency, as shown in Figure 27(b).

The above shows that the final SFE is much less sensitive to undersampling of the IMF than simple models would predict. More importantly, by the time star formation stops, all of our clouds, even the lowest mass ones, have enough stellar mass to sample the IMF. Furthermore, the gas surface density structure of our clouds does not change significantly for varying Ψ, since it is largely set by the initial turbulence. Therefore, by the time the IMF becomes fully sampled, the cloud appears roughly the same regardless of the strength of radiative feedback prior to that time. We might then expect the final efficiency to only depend on the final value of Ψ, which is ${\rm{\Psi }}\sim 2000\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$ in all cases. In fact, when we run simulations in which the value of Ψ changes from ${\rm{\Psi }}\sim 200\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$ at early times to ${\rm{\Psi }}\sim 2000\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$ after a free-fall time, as shown in Figure 28, the final stellar efficiency is largely unchanged.

So far, we have only been considering the final efficiencies of star-forming clouds as a means of assessing when radiative feedback becomes important in unbinding GMCs. However, of just as much interest is the rate of star formation, which is often parameterized by the SFR per free-fall time:

$$\begin{eqnarray}&&{\dot{M}}_{* }(t)\equiv {\varepsilon }_{{\rm{ff}}}(t)\displaystyle \frac{{M}_{g}(t)}{{t}_{{\rm{ff}}}(t)}.\end{eqnarray} \tag{ 23 }$$

We now turn to an analysis of how this SFR varies in our models.

The majority of previous studies have tended to focus only on the mean SFR per free-fall time $\langle {\varepsilon }_{{\rm{ff}}}(t)\rangle$ averaged over the whole epoch of star formation and with a fixed ${t}_{{\rm{ff}}}={t}_{{\rm{ff,0}}}$ set by the cloud's initial mean density (e.g., Wang et al. 2010; Padoan & Nordlund 2011; Bate 2012; Federrath & Klessen 2012; Krumholz et al. 2012; Padoan et al. 2012; Myers et al. 2014). Here, we adopt the methodology of Lee et al. (2015) and fit the stellar mass history of our clouds to determine if there is any systematic evolution in time. We fit with a power law defined in terms of $t-{t}_{* }$, where t_* is the time at which the first star particle is formed. We also initially tested the fitting region of Lee et al. (2015) ranging from ${M}_{* }=0.015{M}_{{\rm{cl,0}}}$ to ${M}_{* }=0.3{M}_{{\rm{cl,0}}}$. Adopting this choice for our fiducial model, similar to Lee et al. (2015) we can confirm a super-linear, though slightly less than quadratic power-law evolution of the SFE, ${M}_{* }\propto {t}^{\beta }$ with $\beta \sim 1.5$, as shown in Figure 29.

Zoom In Zoom Out Reset image size

However, there are a number of issues with fitting a power law to the SFR. First, there is no reason a priori to assume that the gas density distribution or the SFR have reached some sort of steady state at t_*. This is certainly a concern for our simulations, since the cloud begins with uniform density, which is very far from a self-consistent quasi-steady state. It will undergo turbulent collapse, and the highest density regions may form stars while the shape of the density distribution is still changing considerably. Therefore, a power law starting at t_* will not necessarily capture the physics of star formation in realistic clouds as it may still be affected by the artificial initial conditions.

This can be seen very clearly if we consider the stellar mass evolution in our fiducial model as a function of $t-{t}_{* }$, as shown in Figure 29. There is an obvious break in the power law at around ${M}_{* }\approx 0.1{M}_{{\rm{cl,0}}}$ so that a single power law underestimates the stellar mass at both early and late times. This break is not the result of feedback; the model without radiative feedback also shows a similiar break, and an almost identical evolution up to at least $\sim 0.3{M}_{{\rm{cl,0}}}$.

A much better fit can be achieved if we allow for a broken power law ${M}_{* }\propto {(t-{t}_{{\rm{break}}})}^{\beta }$, with a transition time ${t}_{{\rm{break}}}$. For the fiducial model, ${t}_{{\rm{break}}}/{t}_{{\rm{ff,0}}}=0.80$ when ${M}_{* }/{M}_{{\rm{cl,0}}}=0.094$. In this case, as shown in Figure 30(b), we have a much smaller least-squares error for the broken power-law model (${\chi }_{{\rm{broken}}}^{2}$) as opposed to the simple power law (${\chi }_{{\rm{pl}}}^{2}$) with ${\chi }_{{\rm{broken}}}^{2}/{\chi }_{{\rm{pl}}}^{2}=0.042$ over the whole fitting range and ${\chi }_{{\rm{broken}}}^{2}/{\chi }_{{\rm{pl}}}^{2}=0.24$ just above the break mass. More importantly, the broken power law does not systematically depart from the evolution at either low or high stellar mass, and thus more accurately captures the cloud behavior below the transition.

Zoom In Zoom Out Reset image size

Interestingly, when we fit to a broken power law, both regimes show roughly linear growth in the stellar mass with time, albeit at significantly different rates. In fact, if we restrict both regimes to have a constant SFR (as shown in Figure 30(c)) we still get a much better fit than the simple power law above the break time, with ${\chi }_{{\rm{lin}}}^{2}/{\chi }_{{\rm{pl}}}^{2}=0.38$. This suggests two things. First, our fiducial cloud simulation exhibits some transient behavior even after star formation has begun; until close to a free-fall time it seems to be adjusting from the artificial initial state. Second, the cloud appears to emerge from this transient state at approximately t₁₀. It then exhibits roughly linear growth in stellar mass until radiative feedback becomes important.

We may test the presence of this transient state by considering a "pre-relaxed" model, in which the cloud is allowed to evolve without self-gravity for the first half of a free-fall time and then evolves with gravity beyond this point. When this is done, the stellar and outflowing mass evolve as shown in Figure 31. We see that the pre-relaxed model has both a lower final efficiency and lower SFR. It is, however, difficult to compare this result to our fiducial model, since it is at lower surface density due to both expansion of the cloud and the much larger fraction of gas unbound by the initial turbulence (around 30% compared to 10%). Moreover, it is at a lower virial parameter since the turbulence decays away over the first half of a free-fall time of relaxation. What we can see is that there is no evidence of a break in the SFR for this model, suggesting that this break is an artificial one born of our initial conditions.

Zoom In Zoom Out Reset image size

The origin of a constant SFR in later stages is not trivial to explain. The simplest star formation law, and one assumed in a number of previous studies, adopts the form

$$\begin{eqnarray}&&{\dot{M}}_{* }(t)={\varepsilon }_{{\rm{ff}}}\displaystyle \frac{{M}_{g}(t)}{{t}_{{\rm{ff}}}(t)},\end{eqnarray} \tag{ 24 }$$

where ${\varepsilon }_{{\rm{ff}}}$ is a constant and ${t}_{{\rm{ff}}}\propto {\rho }_{g}^{-1/2}$ is the instantaneous free-fall time for ${\rho }_{g}(t)$, the (time-dependent) volume-averaged gas density. A focus of both numerical and observational studies has been to fit to this form and estimate ${\varepsilon }_{{\rm{ff}}}$. In this picture, star formation halts through a combination of gas depletion by star formation and gas expulsion through cloud expansion and a subsequent increase in ${t}_{{\rm{ff}}}$. However, the form of Equation (24) does not seem to apply in our simulations. In our turbulent clouds, the cloud radius remains roughly constant while the majority of stars are forming so that ${\rho }_{g}\propto {M}_{g}$ and ${t}_{{\rm{ff}}}\propto {M}_{g}^{-1/2}$, which would yield $\dot{{M}_{* }}\propto {M}_{g}^{3/2}$ in Equation (24). Therefore, as the gas mass is depleted, Equation (24) would predict a decrease in the SFR until cloud expansion from feedback drives a rapid increase in the free-fall time. In fact, our power-law or broken power-law fits show that the SFR is constant or increasing until gas is expelled by feedback.

Fits to the simple star formation law of Equation (24) are shown in Figure 30(d). The gas mass M_g(t) is taken directly from each simulation, and ${t}_{{\rm{ff}}}(t)$ is found by fitting for the mean cloud density through the density PDF. If we fit this form with star formation beginning at $t={t}_{* }$, the best-fit value is ${\varepsilon }_{{\rm{ff}}}=0.45$, and there is a huge discrepancy between the numerical and analytic results. Even if we fit only after $t={t}_{{\rm{break}}}$, the fit is not good; the best-fit efficiency ${\varepsilon }_{{\rm{ff}}}=1.21$ is high since the early mass growth must be large to compensate for the steady decrease in SFR.

We may test the generality of these conclusions by applying the same methodology to other members of our Σ-series simulations. As a caveat, we saw in Section 3.2 that lower resolution simulations will capture the final SFE, but may underestimate the SFR if they are not converged. This is particularly evident in our low surface density and high virial parameter simulations in which the number of individual star particles formed is small. With this in mind, we omit simulations with $\varepsilon \lt 0.15$ or ${{\rm{\Sigma }}}_{{\rm{cl,0}}}\lt 20\,{M}_{\odot }\,{{\rm{pc}}}^{-2}$ from our studies of the SFR.

In Figure 32 we show the best-fit characteristics for both the single power-law fit and the high-mass portion of our broken power-law fit, as a function of surface density. Values for β and t_break are also listed in Table 3. For the Σ-series, we observe little variation in either the power-law exponent or the break time and mass across the sequence. As for our fiducial model, for the broken power law, the fitted exponent is close to $\beta =1$ for most cases. Meanwhile, a single power law almost always shows super-linear behavior, with $\beta \sim 1.5$.

Zoom In Zoom Out Reset image size

At the break time, roughly $0.1\,{M}_{{\rm{cl,0}}}$ in stellar mass has formed across all models. The break times also show little variation with surface density, being close to ${t}_{{\rm{break}}}\sim 0.8{t}_{{\rm{ff,0}}}$ for the broken power law, with ${t}_{* }\sim 0.5{t}_{{\rm{ff,0}}}$. The most massive clouds begin star formation a little earlier, since they are already denser, hence regions reach critical densities high enough to undergo local collapse at earlier times. This slight surface density dependence also persists in the break time for the broken power law, suggesting that it may be as arbitrary as t_*. Potentially, the initial conditions still affect the state of the system at ${t}_{{\rm{break}}}$, although our simulations do not show any strong evidence for a continuous acceleration in the SFR. Nevertheless, simulations with more realistic cloud initial conditions are needed before any definitive statements can be made about varying SFRs in turbulent clouds.

The fact that our star formation efficiencies grow roughly linearly with time allows us to quantify these SFRs quite easily. We replicate the method of Padoan & Nordlund (2011), Krumholz et al. (2012), and Myers et al. (2014) and fit a straight line to the stellar mass versus time above the transition mass discussed in Section 4.5. We then calculate the efficiency per free-fall time defined as

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{ff},\bar{\rho }}\equiv \displaystyle \frac{\langle \dot{{M}_{* }}\rangle {t}_{\mathrm{ff},\bar{\rho }}}{{M}_{{\rm{cl,0}}}}\end{eqnarray} \tag{ 25 }$$

where the mean density used in ${t}_{\mathrm{ff},\bar{\rho }}$ is calculated from directly fitting a log-normal form to the density PDF at the start of star formation, defined here as ${t}_{{\rm{break}}}$, and extracting the mean density $\langle \rho \rangle$. For the Σ-series, this is very close to the initial cloud density since only around $10 \%$ of the mass is lost in early turbulent outflows and there is very little global cloud contraction or expansion. However, the mean density differs significantly from the initial cloud value for the models with low virial parameter, where the cloud contracts significantly.

In Figure 33 (see also Table 3) we show the resulting SFR coefficients as a function of both surface density and virial parameter. Almost irrespective of both virial parameter and surface density, we find a rate coefficient ${\varepsilon }_{\mathrm{ff},\bar{\rho }}\sim 0.25\mbox{--}0.5$. Notably, the inclusion of radiative feedback only mildly decreases the SFR, with no systematic surface density dependence. This is not surprising as we had already noted that turbulence dominates the star formation at early times. These results again suggest that UV feedback in clouds can be effective as a means of limiting star formation and unbinding clouds, but does little to suppress the instantaneous SFR.

Zoom In Zoom Out Reset image size

The consequence of such high SFRs is that our model clouds either convert most of their gas mass to stars, or become unbound on short timescales. While it is difficult to characterize exactly when the bulk of gas mass from a filamentary cloud becomes unbound, the proxy we use is the virial parameter. Similar to Colín et al. (2013), we find that in all our Σ-series simulations, the virial parameter remains close to unity when the dynamics are set by gravitational collapse and turbulence, but then quickly expand to ${\alpha }_{{\rm{vir}}}\sim 50$ once radiative feedback begins to dominate. Therefore, we arbitrarily take ${\alpha }_{{\rm{vir}}}=5$ as our criterion for unboundedness, noting that the expansion from being formally unbound at ${\alpha }_{{\rm{vir}}}=2$ to ${\alpha }_{{\rm{vir}}}\gtrsim 10$ takes $\sim 0.3{t}_{{\rm{ff,0}}}$.

In Figure 34(a), we show the time ${t}_{{\rm{unb}}}$ when clouds become unbound (${\alpha }_{{\rm{vir}}}=5$) as a function of surface density. We find that clouds last between 1.2 and $1.9\,{t}_{{\rm{ff,0}}}$ at most and that radiative feedback acts very rapidly to unbind the clouds, with the transition from ${\alpha }_{{\rm{vir}}}=2$ to ${\alpha }_{{\rm{vir}}}=10$ never taking longer than half a free-fall time (see also Table 3). Star formation still continues slowly as clouds continue to expand beyond ${\alpha }_{{\rm{vir}}}=5$, but if we look at the time t₈₀ as shown in Figure 34(a), it is generally similar to ${t}_{{\rm{unb}}}$, while t₉₀ roughly corresponds to clouds reaching a virial parameter of ∼10.

Zoom In Zoom Out Reset image size

Given that the first stars only begin forming at around $\sim 0.5{t}_{{\rm{ff,0}}}$, and clouds only reach a state where the artificial initial conditions are erased at $\sim 0.8{t}_{{\rm{ff,0}}}$ (where this number is taken from both the minimum in the virial parameter, and the break times in fits to the SFR), this means that the majority of clouds form stars over a period shorter than a global free-fall time.

This can be seen if we consider the cloud star-forming time, defined as ${t}_{{\rm{unb}}}-{t}_{* }$, in Figure 34(b). Since ${\varepsilon }_{\mathrm{ff},\bar{\rho }}$ varies mildly but ${\varepsilon }_{{\rm{final}}}$ increases more strongly with Σ_cl,0 there is an increase in cloud lifetime in units of ${t}_{{\rm{ff,0}}}$ with Σ_cl,0 as the more massive clouds will convert more of their gas to stars before the effects of radiative feedback dominate and, additionally, begin forming stars at a slightly earlier time. However, as high-Σ_cl,0 clouds also have much shorter free-fall times, we find that, generically, cloud star formation times (${t}_{{\rm{unb}}}-{t}_{* }$ or ${t}_{{\rm{90}}}-{t}_{* }$) are between 2 and $8\,$ Myr (see Figure 34(c)), with more massive clouds being slightly shorter lived than their low surface density counterparts. Since SN feedback can only act more than 3 Myr after $\sim {t}_{* }$, it seems likely that if there were no other source of star formation suppression, direct radiation pressure would be able to disperse clouds across a wide range of surface densities before SNe begin to impact cloud dynamics. However, this cloud dispersal would be at the expense of a larger net SFE than inferred from observations in many Milky Way clouds.

This apparent conflict with observation may potentially be explained in several ways. One possibility is that GMCs are strongly affected by additional sources of internal feedback not modeled here (such as ionizing radiation¹ ) or external feedback (including supernova blast waves from stars formed in other GMCs). Another is that GMCs may have effective virial parameters exceeding 2, e.g., if the outer parts are still condensing even as the inner parts begin vigorous star formation.

To test our code behavior in highly idealized conditions, we consider the expansion of a spherical shell of gas due to the absorption of radiation momentum from a central source, similar to that described in Ostriker & Shetty (2011) and Skinner & Ostriker (2013). In this problem, we imagine an idealized spherical GMC of mass ${M}_{{\rm{cl,0}}}$ that forms stars of total mass M_* with efficiency $\varepsilon ={M}_{* }/{M}_{{\rm{cl,0}}}$. The remaining gas of mass ${M}_{{\rm{sh}}}\equiv (1-\varepsilon ){M}_{{\rm{cl,0}}}$ is ejected as an expanding, spherical shell of radius r due to the radiation force from the stellar component. The stars are modeled here as a centrally located cluster.

We initialize a shell at an initial radius r₀ with zero velocity at time t = 0. In reality, the shell would have some initial velocity, but this just provides an additive constant. Assuming that the ejected shell is thin and of uniform surface density ${\rm{\Sigma }}(r)={M}_{{\rm{sh}}}/(4\pi {r}^{2})$, the optical depth across the shell is ${\tau }_{{\rm{sh}}}(r)\approx {\rm{\Sigma }}(r)\kappa$, where κ is the absorption opacity of the gas to UV photons. If, as discussed earlier, we are in the limit where we consider only UV radiation, then the flux at r is given by

$$\begin{eqnarray}&&F(r)=\displaystyle \frac{{L}_{* }}{4\pi {r}^{2}}{\rm{\exp }}\left[-{\int }^{r}\rho (r^{\prime} )\kappa {dr}^{\prime} \right]\equiv \displaystyle \frac{{L}_{* }{e}^{-{\tau }_{{\rm{sh}}}(r)}}{4\pi {r}^{2}}.\end{eqnarray} \tag{ 26 }$$

The total radiation force on the shell is then

$$\begin{eqnarray}&&\int F(r)\rho (r)\displaystyle \frac{\kappa }{c}4\pi {r}^{2}{dr}=\displaystyle \frac{{L}_{* }}{c}\int {e}^{-\tau }d\tau =\displaystyle \frac{{L}_{* }}{c}(1-{e}^{-{\tau }_{{\rm{sh}}}({r}_{{\rm{\max }}})}).\end{eqnarray} \tag{ 27 }$$

Neglecting gravitational and internal pressure forces, the outward acceleration of the shell becomes

$$\begin{eqnarray}&&\ddot{r}=\displaystyle \frac{{L}_{* }[1-{\rm{\exp }}(-{\rm{\Sigma }}(r)\kappa )]}{{M}_{{\rm{sh}}}c},\end{eqnarray} \tag{ 28 }$$

which is independent of the shell's thickness. Substituting for the surface density and the cluster luminosity, the shell acceleration reduces to

$$\begin{eqnarray}&&\ddot{r}=\displaystyle \frac{{\rm{\Psi }}\varepsilon (1-{e}^{-{\tau }_{0}{({r}_{0}/r)}^{2}})}{c(1-\varepsilon )},\end{eqnarray} \tag{ 29 }$$

where we have introduced the shell optical depth at $r={r}_{0}$, given by

$$\begin{eqnarray}{\tau }_{0} & \equiv & {M}_{{\rm{sh}}}\kappa /(4\pi {r}_{0}^{2})\\ & = & 1.67{\left(\displaystyle \frac{{r}_{0}}{10{\rm{pc}}}\right)}^{-2}\left(\displaystyle \frac{\kappa }{1000\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}}\right)\left(\displaystyle \frac{{M}_{{\rm{sh}}}}{{10}^{4}\,{M}_{\odot }}\right).\end{eqnarray} \tag{ 30 }$$

We may simplify further by rewriting Equation (29) in terms of the dimensionless variables $\tilde{r}\equiv r/{r}_{0}$ and $\tilde{t}\equiv t/{t}_{0}$, so that

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}\tilde{r}}{d{\tilde{t}}^{2}}=\displaystyle \frac{\varepsilon (1-{e}^{-{\tau }_{0}/{\tilde{r}}^{2}})}{(1-\varepsilon )},\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}{t}_{0} & \equiv & \sqrt{\displaystyle \frac{{r}_{0}c}{{\rm{\Psi }}}}\\ & = & 0.68\,{\rm{Myr}}{\left(\displaystyle \frac{{r}_{0}}{10{\rm{pc}}}\right)}^{1/2}{\left(\displaystyle \frac{{\rm{\Psi }}}{2000\mathrm{erg}{{\rm{s}}}^{-1}{{\rm{g}}}^{-1}}\right)}^{-1/2}.\end{eqnarray} \tag{ 32 }$$

The presence of the exponential term in inverse radius precludes a general analytic solution to this problem (but see Raskutti et al 2016 for a velocity-radius relation). However, in the optically thick (${\tau }_{{\rm{sh}}}\gg 1$) limit, or close to the initial shell radius with ${\tau }_{{\rm{sh}}}(r)\approx {\tau }_{0}$, all explicit dependence on radius drops out of Equation (31). For the $\tau \approx {\tau }_{0}$ case it may then be solved trivially to give a quadratic expansion in time

$$\begin{eqnarray}&&\tilde{r}=\displaystyle \frac{\varepsilon (1-{e}^{-{\tau }_{0}})}{2(1-\varepsilon )}{\tilde{t}}^{2}+1.0;\end{eqnarray} \tag{ 33 }$$

in the optically thick case, we instead have $1-{e}^{-{\tau }_{0}}\to 1$. We note that if the shell is optically thick, then to first order the solution depends only on the inital shell radius, SFE ε, and luminosity per unit mass Ψ (see Equation (29)). Otherwise, dependence on the mass and opacity only enters through ${\tau }_{0}$ (see Equation (30)) and is only significant for relatively low optical depths.

A.1.1. Convergence Tests

In this test, we wish to explore the sensitivity of our code to changes in the key physical parameters κ and ${M}_{{\rm{cl,0}}}$, as well as numerical parameters N, and $\hat{c}$. We consider a central luminous cluster defined by a sink particle at the origin, with ${r}_{* }=1$ pc for the radiation source function given in Equation (5). This central source illuminates a thin shell with Gaussian density profile given by

$$\begin{eqnarray}&&{\rho }_{{\rm{sh}}}(r)=\displaystyle \frac{{M}_{{\rm{sh}}}}{4\pi {r}_{0}^{2}\,\sqrt{2\pi {\sigma }_{{\rm{sh}}}^{2}}}\exp \left(-\displaystyle \frac{{(r-{r}_{0})}^{2}}{2{\sigma }_{{\rm{sh}}}^{2}}\right),\end{eqnarray} \tag{ 34 }$$

where ${H}_{0}=2\sqrt{2\mathrm{ln}2}{\sigma }_{{\rm{sh}}}$ is the FWHM of the shell.

For this test, we initialize the radiation field using the result of a prior simulation in which we turn on the cluster and evolve the radiation field without evolving the gas hydrodynamics. The final radiation field after ∼10 radiation crossing times can then be used to initialize the flux and energy density for the case where we consider shell expansion.

To prevent the gas timesteps from becoming prohibitively small, we enforce a density floor of ${\rho }_{{\rm{\min }}}\equiv {10}^{-8}{\rho }_{{\rm{sh}}}(r={r}_{0})$ initially, as well as after each gas integration step. We employ an N³ grid on the domain $(x,y,z)\in ([0,3\,{r}_{0}],[0,3\,{r}_{0}],[0,3\,{r}_{0}])$ and enforce outflow boundary conditions along all faces of the box that do not touch the origin of the cloud. We note that due to spherical symmetry, we only run the test in a single octant of the sphere. We run each simulation for a time $t\approx 2{t}_{0}$ so that the planar shell reaches an outermost radius of $r\approx 3{r}_{0}$.

For the RSLA, it is necessary to choose a value of $\hat{c}$ such that ${v}_{{\rm{\max }}}\ll \hat{c}$ at all times. For each of our simulations, the shell reaches a maximum velocity of ${v}_{{\rm{\max }}}=\varepsilon {v}_{0}{\tilde{t}}_{{\rm{\max }}}/(1-\varepsilon )$, where

$$\begin{eqnarray}{v}_{0} & \equiv & \sqrt{\displaystyle \frac{{\rm{\Psi }}{r}_{0}}{c}}\\ & = & 14.4\,\mathrm{km}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{{r}_{0}}{10{\rm{pc}}}\right)}^{1/2}{\left(\displaystyle \frac{{\rm{\Psi }}}{2000\mathrm{erg}{{\rm{s}}}^{-1}{{\rm{g}}}^{-1}}\right)}^{1/2}.\end{eqnarray} \tag{ 35 }$$

As we use $\varepsilon =0.5$ and ${\tilde{t}}_{{\rm{\max }}}=2$, a choice of $\hat{c}=250\,\mathrm{km}\,{{\rm{s}}}^{-1}$ should be satisfactory for the RSLA in most situations. Of course, increased values of Ψ or the SFE ε will provide greater accelerations and so require higher values of $\hat{c}$. The robustness of our results to variations in $\hat{c}$ is therefore also verified below.

We begin by adopting a set of fiducial parameters roughly characteristic of Milky Way clouds. In addition to ${\rm{\Psi }}\,=2000\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$ and $\kappa =1000\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, we choose an initial shell radius of ${r}_{0}=10$ pc, a relatively thin initial shell width of ${H}_{0}=1.0$ pc, and a cloud mass of ${M}_{{\rm{cl,0}}}={10}^{4}\,{M}_{\odot }$. This mass provides a marginally optically thick shell ${\tau }_{0}\approx 2$ (see Equation (30)). Finally, we adopt a SFE of $\varepsilon =0.5;$ although this is high for Milky Way GMCs, it is similar to the upper range of what we find in our full cloud simulations.

For this test we adopt a sound speed of ${c}_{s}=0.5\,\mathrm{km}\,{{\rm{s}}}^{-1}$, which lies somewhere between the turbulent velocities in our full cloud simulations and the true sound speed (a factor of 2 lower). This helps ensure that when resolution is adequate, the shell does not develop thin-shell instability as it is accelerated outward, which would compromise our ability to compare to the analytic spherical solution.

We are initially interested in the numerical parameters N and $\hat{c}$ required to match the analytic thin-shell solution for a typical GMC. In Figure 35 we show the results for several fiducial simulations run at different numerical resolutions. Pictured are snapshots of density in slices through the x–y plane at $t=1.2{t}_{0}$ when the mean radius is $r\sim 1.7{r}_{0}$ for resolutions N = 64, 128, 256, and 512. We note that the shell width remains roughly constant or only slightly expands to a final shell width of $H\sim 2\,\mathrm{pc}$. However, in the lowest resolution simulation, N = 64, the spherical shell is disturbed by grid scale noise. By the time the shell has expanded to close to twice its initial radius, it is no longer spherical, but instead, has large-scale perturbations with angle caused by the initial difficulty of resolving a spherical shell on a square grid.

Zoom In Zoom Out Reset image size

Figure 36(a) shows the evolution of mass-weighted shell radius with time as compared to the analytic solution at varying resolution for a cloud of mass ${10}^{4}\,{M}_{\odot }$. Since the optical depth has a radial dependence that appears through the surface density, we may only find the analytic solution by numerically integrating Equation (31). Evidently, the numerical solution follows the analytic prediction fairly well even at low resolution, although for higher curvature structures, the accuracy would be reduced at each given resolution. The main conclusion from this test is that at the typical resolution of our simulations, we are able to satisfactorily capture the predicted expansion driven by radiation forces.

Zoom In Zoom Out Reset image size

The primary other numerical parameter that may affect our results is the reduced speed of light $\hat{c}$. In Hyperion simulations measuring the effect of reprocessed radiation, results can be quite sensitive to this parameter, since the RSLA static diffusion criterion requires that the effective radiation diffusion speed remain large compared to dynamical speeds, i.e., $\hat{c}/{\tau }_{{\rm{\max }}}\gg {v}_{{\rm{\max }}}$, where ${\tau }_{{\rm{\max }}}$ is the maximum optical depth across all cells. However, for direct radiation the RSLA criterion is simply $\hat{c}\gg {v}_{{\rm{\max }}}$, where the highest surface density clouds we consider typically have ${v}_{{\rm{\max }}}\sim {v}_{{\rm{esc}}}\sim 15\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The shell expansion test with the fiducial cloud parameters described above has ${v}_{{\rm{\max }}}\sim \sqrt{2}{v}_{0}=20\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Figure 36(b), showing results for varying $\hat{c}$, demonstrates that there is only a small error with respect to the analytic solution for $\hat{c}\sim 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and we recover it exactly for $\hat{c}=250\,\mathrm{km}\,{{\rm{s}}}^{-1}$. As discussed earlier, we have conservatively adopted the latter value for all our cloud simulations.

Figure 37(a) shows the same evolution of shell radius as Figure 36(a), but for a more massive cloud with ${M}_{{\rm{cl,0}}}=3\,\times {10}^{5}\,{M}_{\odot }$, i.e., thirty times larger than in Figure 36(a). The mass of the central cluster is again half this at $M=1.5\times {10}^{5}\,{M}_{\odot }$. Evidently, the numerical tests no longer agree closely with the analytic solutions for higher mass shells, for numerical resolution N in the range shown. In all cases, the numerical solution systematically underestimates the shell radius. While increasing the numerical resolution may marginally improve the agreement between analytic and numerical solutions, the error is still close to ∼10% once the shell has reached the edge of the simulation volume. Meanwhile, if we consider the same test for varying $\hat{c}$ (Figure 37(b)), we see that increasing the reduced speed of light beyond our adopted value of $\hat{c}=250\,\mathrm{km}\,{{\rm{s}}}^{-1}$ does not remedy this discrepancy either. As further tests show (see below), this discrepancy can in fact be traced to inadequate resolution of the flux.

Zoom In Zoom Out Reset image size

A.1.2. Flux Resolution

The optical depth across individual cells in our grid is ${\tau }_{{\rm{cell}}}=\rho \kappa {\rm{\Delta }}x$. If this optical depth exceeds unity by a considerable amount, then the flux across cells is not well resolved spatially and the impulse provided to the gas is not captured accurately. For the shell problem, with ${\rho }_{{\rm{sh}}}(r)\approx {\rm{\Sigma }}/H$ and ${\rm{\Delta }}x=4{r}_{0}/N$,

$$\begin{eqnarray}&&{\tau }_{{\rm{cell}}}\to \displaystyle \frac{4{\rm{\Sigma }}\kappa {r}_{0}}{{{HN}}_{x}};\end{eqnarray} \tag{ 36 }$$

this depends on the opacity, the shell surface density (and hence mass and radius), and the numerical resolution.

We may systematically test the dependence of the solution accuracy on ${\tau }_{{\rm{cell}}}$ by using the shell problem with varying κ and Σ (through ${M}_{{\rm{cl,0}}}$). Figure 38(a) shows the result from varying opacity over two orders of magnitude between 10³ and ${10}^{5}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. Above $\kappa \sim {10}^{4}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$, the shell is completely optically thick and there is no variation in the analytic solution with opacity. The numerical results match the analytic solutions well for opacities $\kappa =2000\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ and below, but then they deteriorate as the opacity increases to $\kappa =5000\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. Interestingly, beyond this opacity, even up to $\kappa ={10}^{5}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ there is not a significant further increase in the error relative to the analytic solution, and there is only a 10% deviation in radius out to two times the initial shell radius (where a real cloud could have become gravitationally unbound).

Zoom In Zoom Out Reset image size

The underlying reason for these deviations can be seen in the maximum optical depth across individual cells as shown in Figure 38(b). In all cases, this increases by up to a factor ∼10 from its starting value as the shell is compressed during expansion, and then decreases as the shell expands outwards and the mean shell density decreases. More importantly, the value of the maximum optical depth plays a key role in determining solution accuracy. We see that for $\kappa =2000$ and $\kappa =5000\,{\rm{g}}\,{\mathrm{cm}}^{-2}$, the optical depths peak at ${\tau }_{{\rm{cell}}}=1.3$ and ${\tau }_{{\rm{cell}}}=4.3$, respectively, and somewhere between these, there is a transition point at ${\tau }_{{\rm{cell}}}\sim 2-3$ beyond which the flux is not well resolved. For $\kappa =5000\,{\rm{g}}\,{\mathrm{cm}}^{-2}$, ${\tau }_{{\rm{cell}}}\gtrsim 2$ for around $0.5{t}_{0}$, and the numerical simulation differs from the analytic solution by around 5%. For larger $\kappa ={10}^{4}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$, ${\tau }_{{\rm{cell}}}\gtrsim 2$ for the majority of the shell evolution, which leads to the $10 \%$ errors discussed earlier. Beyond this point, increasing the opacity does not have a strong effect on the solution accuracy, since ${\tau }_{{\rm{cell}}}\gtrsim 2$ always.

Figure 39 shows results of similar tests, in which we vary the surface density (through ${M}_{{\rm{cl,0}}}$). Over a range of two orders of magnitude above the fiducial surface density, we see the same trends with cell optical depth. For ${M}_{{\rm{cl,0}}}\lesssim 3\times {10}^{4}\,{M}_{\odot }$, for which ${\tau }_{{\rm{cell}}}\leqslant 2$ at all times, we match the analytic solution reasonably well. However, at larger masses and correspondingly larger ${\tau }_{{\rm{cell}}}$, we again underestimate the shell expansion velocity at all times.

Zoom In Zoom Out Reset image size

We conclude that, provided ${\tau }_{{\rm{cell}}}$ remains below ∼2, we can obtain an accurate solution for this problem. Since ${\tau }_{{\rm{cell}}}$ is inversely proportional to N, we can in principle capture the behavior in increasingly high-density clouds by increasing the numerical resolution, although in practice this becomes numerically prohibitive for very dense systems. In any case, regardless of how high ${\tau }_{{\rm{cell}}}$ becomes, we never underpredict the shell radius by more than $10 \%$ even by the time gas in one of our turbulent clouds will have become unbound.

Finally, we note that the maximum values of ${\tau }_{{\rm{cell}}}$ depend strongly on gas compression and hence, on the detailed problem-specific evolution of turbulent clouds. Therefore, while the spherical tests show that a resolution of N = 256 is generally sufficient for this problem, we also need to directly test convergence in our full turbulent models.