Modeling UV Radiation Feedback from Massive Stars. II. Dispersal of Star-forming Giant Molecular Clouds by Photoionization and Radiation Pressure

For hydrodynamics, we employ Athena's van Leer-type time integrator (Stone & Gardiner 2009), the HLLC Riemann solver, and a piecewise linear spatial reconstruction scheme. We use the fast Fourier transformation Poisson solver with open boundary conditions (Skinner & Ostriker 2015) to calculate the gravitational forces from gas and stars; the stellar contribution is handled with the particle-mesh method using the triangular-shaped-cloud interpolation scheme (Hockney & Eastwood 1981).

We apply the sink particle technique of Gong & Ostriker (2013) to model the formation and growth of star clusters. We create a star particle if a cell with density above the Larson–Penston threshold density (${\rho }_{\mathrm{crit}}\equiv 8.86{c}_{{\rm{s}},\mathrm{neu}}^{2}/(\pi G{\rm{\Delta }}{x}^{2})$ for the sound speed of neutral gas ${c}_{{\rm{s}},\mathrm{neu}}$ and the grid spacing ${\rm{\Delta }}x$) has a converging velocity field and corresponds to the local minimum of the gravitational potential. The accretion rates of mass and momentum onto each sink are computed from the flux through the surfaces of the 3³-cell control volume centered at the sink particle; this control volume acts like internal ghost zones within the simulation domain. The positions and velocities of sink particles are updated using a leapfrog integrator. We allow sink particles to merge if their control volumes overlap.

Due to limited resolution, the sink particles in our simulations represent subclusters rather than individual stars,³ and may not fully sample the initial mass function. Using the mass-dependent light-to-mass ratios from Kim et al. (2016), we determine the total UV luminosity of cluster particles for two frequency bins, ${L}_{{\rm{i}}}$ and ${L}_{{\rm{n}}}$, representing Lyman continuum photons and far-UV photons, such as $L\equiv {L}_{{\rm{i}}}+{L}_{{\rm{n}}}={\rm{\Psi }}{M}_{* }$ and ${Q}_{{\rm{i}}}\equiv {L}_{{\rm{i}}}/(h{\nu }_{{\rm{i}}})={\rm{\Xi }}{M}_{* }$, where M_* is the total cluster mass, ${\rm{\Psi }}={10}^{2.98{{ \mathcal X }}^{6}/(29.0+{{ \mathcal X }}^{6})}\,{L}_{\odot }\,{M}_{\odot }^{-1}$, ${\rm{\Xi }}={10}^{46.7{{ \mathcal X }}^{7}/(7.28+{{ \mathcal X }}^{7})}\,{{\rm{s}}}^{-1}\,{M}_{\odot }^{-1}$ with ${ \mathcal X }={\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })$, and $h{\nu }_{{\rm{i}}}=18\,\mathrm{eV}$ is the mean energy of ionizing photons. Note that ${\rm{\Psi }}\simeq 943\,{L}_{\odot }\,{M}_{\odot }^{-1}$ and ${\rm{\Xi }}\simeq 5.05\times {10}^{46}\,{{\rm{s}}}^{-1}\,{M}_{\odot }^{-1}$ are almost constant for ${M}_{* }\,\gtrsim {10}^{3}\,{M}_{\odot }$, while varying steeply for ${M}_{* }\lt {10}^{3}\,{M}_{\odot }$ (see Figure 14 of Kim et al. 2016). The luminosity of individual cluster particles is assumed to be proportional to their mass. The use of mass-dependent conversion factors approximately captures the effect of undersampling at the massive end of the initial mass function when the total cluster mass is below $\sim {10}^{3}\,{M}_{\odot }$.

In the present paper, we do not consider the evolution of the stellar light-to-mass ratio. The UV luminosity of a coeval stellar population that well samples the initial mass function remains approximately constant for the main sequence lifetime of the most massive star, implying that ionizing and non-ionizing luminosities would in practice decay rapidly after timescales of $3\,\mathrm{Myr}$ and $8\,\mathrm{Myr}$, respectively (e.g., Parravano et al. 2003). Caveats related to our adoption of age-independent luminosity are discussed in Section 6.2.

The radiative transfer of ionizing and non-ionizing photons is handled by the adaptive ray tracing method (Abel & Wandelt 2002). For each radiating star particle, we generate 12 × 4⁴ = 3072 initial rays and allow photon packets to propagate radially outward, computing the cell-ray intersection length and the corresponding optical depth on a cell-by-cell basis. The volume-averaged radiation energy and flux densities of a cell are computed from the sum of the contributions from all rays that pass through the cell (see Section 2.1 of Paper I). These averages are used to compute the radiation force from the combined non-ionizing and ionizing radiation fields, and the ionization rate of neutral gas. Because the fluid variables and gravity in the 3³ cells surrounding each sink/source particle are unresolved, we do not allow photon packets to interact with the gas within these control volumes. The radiation field directions are angularly well resolved, and the rays contain all of the photon packets from each source, when they emerge from the control volume around each source particle. This eliminates the potential problem of momentum cancellation caused by volume-averaging in the cell containing a source, noted recently by Hopkins & Grudic (2018), without the need to introduce assumptions regarding the sub-cell distributions of gas density, gravity, or the radiation field within individual zones that contain sink/source particles.

The rays are discretized and split based on the HEALPix framework (Górski et al. 2005), ensuring that each grid cell is sampled by at least four rays per source. We rotate ray propagation directions randomly at every hydrodynamic time step to reduce the numerical errors due to angle discretizations (Krumholz et al. 2007). For simplicity, we use constant photoionization cross-section for neutral hydrogen atom ${\sigma }_{{{\rm{H}}}^{0}}=6.3\times {10}^{-18}\,{\mathrm{cm}}^{-2}$ (Krumholz et al. 2007)⁴ and constant dust absorption (and pressure) cross-section per H nucleon ${\sigma }_{{\rm{d}}}=1.17\times {10}^{-21}\,{\mathrm{cm}}^{-2}\,{{\rm{H}}}^{-1}$ both for ionizing and non-ionizing radiation (e.g., Draine 2011). We determine the degree of hydrogen ionization by solving the time-dependent rate equation for hydrogen photoionization and recombination (under the on-the-spot approximation). The ray-trace and ionization update are subcycled relative to the hydrodynamic update. The time step size for each substep is determined to ensure that the maximum change in the neutral fraction ${x}_{{\rm{n}}}={n}_{{{\rm{H}}}^{0}}/{n}_{{\rm{H}}}$ is less than 0.1. The source term updates for the ionization fraction and radiative force are explicit and performed at every substep in an operator-split fashion. The gas temperature is set to vary according to ${x}_{{\rm{n}}}$ between ${T}_{\mathrm{neu}}=20\,{\rm{K}}$ and ${T}_{\mathrm{ion}}=8000\,{\rm{K}}$—the temperature of the fully neutral and ionized gas, respectively. The corresponding isothermal sound speeds are c_s,neu = 0.26 km s⁻¹ and c_s,ion = 10.0 km s⁻¹.

We note that determining gas temperature solely based on the neutral fraction, together with adoption of constant cross-sections for photoionization and dust absorption, will certainly simplify temperature distributions inside H ii regions compared to what would be more complex in realistic environments. To properly model ionization and temperature structure, as well as small-scale dynamical instabilities of IFs that may operate (e.g., Kim & Kim 2014; Baczynski 2015), one has to accurately solve the energy equation after considering various cooling and heating processes, as well as the frequency dependence of the cross-sections.

Our initialization and boundary treatment are the same as in Skinner & Ostriker (2015) and Raskutti et al. (2016). The initial conditions for each model consist of an isolated, uniform density sphere of neutral gas with mass ${M}_{0}$ and radius ${R}_{0}$ surrounded by a tenuous external medium situated in a cubic box with sides ${L}_{\mathrm{box}}=4{R}_{0}$.⁵ The cloud is initially supplied with turbulent energy with velocity power $| \delta {{\boldsymbol{v}}}_{k}{| }^{2}\propto {k}^{-4}$ for $2\leqslant {{kL}}_{\mathrm{box}}/(2\pi )\leqslant 64$ and zero power otherwise. We adjust the amplitude of the velocity perturbations to make the cloud marginally bound gravitationally, with the initial virial parameter of ${\alpha }_{\mathrm{vir},0}=5{\sigma }_{0}^{2}{R}_{0}/(3{{GM}}_{0})=2$, where ${\sigma }_{0}$ denotes the initial (three-dimensional) velocity dispersion. There is no subsequent artificial turbulent driving. The standard resolution for the simulation domain is N_cell = 256³ cells.

We apply diode-like outflow boundary conditions both at the boundaries of the simulation box and at the boundary faces of the control volume of each sink particle. The diode-like conditions allow mass and momentum to flow into ghost zones but not out from them. When sink particles accrete while moving across grid zones, we ensure that the total mass and momentum of gas and stars are conserved by taking into account differences in the mass and momentum of cells entering and leaving the control volume. Over the course of evolution, we separately monitor the mass outflow rates of neutral and ionized gas from the computational domain, as well as escape fractions of ionizing and non-ionizing photons.

Our aim is to explore the effects of radiation feedback on cloud dispersal across a range of cloud masses and sizes, corresponding to GMCs observed in the Milky Way and nearby galaxies. These clouds are usually optically thick to UV radiation (${\rm{\Sigma }}\gtrsim {\kappa }_{\mathrm{UV}}^{-1}\sim {10}^{1}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$) but transparent to dust-reprocessed infrared radiation (${\rm{\Sigma }}\lesssim {\kappa }_{\mathrm{IR}}^{-1}\sim {10}^{3}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$). The momentum deposition by multiple scatterings of infrared photons may play an important role in the dynamics of extremely dense and massive clouds, such as the progenitors of super star clusters, in dust-rich environments (e.g., Skinner & Ostriker 2015; Tsang & Milosavljevic 2017).

We consider 14 clouds with ${M}_{0}$ in the range of ${10}^{4}\,{M}_{\odot }$ to ${10}^{6}\,{M}_{\odot }$ and ${R}_{0}$ from 2 to $80\,\mathrm{pc}$: the resulting initial surface density ${{\rm{\Sigma }}}_{0}={M}_{0}/(\pi {R}_{0}^{2})$ is in the range from $12.7\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ to $1.27\times {10}^{3}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$. Figure 1 plots the locations (red squares) of our model clouds in the M–Σ domain, compared to observed GMCs (circles) compiled from the literature (note that a range of different methods has been adopted in the literature to estimate observed cloud surface density and mass). Table 1 lists the initial physical parameters of our model clouds. Column 1 gives the name of each model. Columns 2–4 list ${M}_{0}$, ${R}_{0}$, and ${{\rm{\Sigma }}}_{0}$, respectively. Columns 5 and 6 give the number density of hydrogen atoms ${n}_{{\rm{H}},0}$ and the freefall time ${t}_{\mathrm{ff},0}=\sqrt{3\pi /(32G{\rho }_{0})}$. Finally, Column 7 gives the escape velocity ${v}_{\mathrm{esc},0}=\sqrt{2{{GM}}_{0}/{R}_{0}}$ at the cloud surface. The initial turbulent velocity dispersion in each cloud is ${\sigma }_{0}=0.77{v}_{\mathrm{esc},0}$. We take model M1E5R20 with ${M}_{0}={10}^{5}\,{M}_{\odot }$, ${R}_{0}=20\,\mathrm{pc}$, and ${t}_{\mathrm{ff},0}=4.7\,\mathrm{Myr}$ as our fiducial model, which is typical of Galactic GMCs and comparable in mass and size to the Orion A molecular cloud (e.g., Wilson et al. 2005).

Zoom In Zoom Out Reset image size

Table 1.  Model Cloud Parameters

Model	${M}_{0}$	${R}_{0}$	${{\rm{\Sigma }}}_{0}$	${n}_{{\rm{H}},0}$	${t}_{\mathrm{ff},0}$	${v}_{\mathrm{esc},0}$
(1)	(2)	(3)	(4)	(5)	(6)	(7)
M1E5R50	1 × 105	50.0	12.7	5.5	18.5	4.1
M1E5R40	1 × 105	40.0	19.9	10.8	13.2	4.6
M1E5R30	1 × 105	30.0	35.4	25.5	8.6	5.4
M1E4R08	1 × 104	8.0	49.7	134.7	3.7	3.3
M1E6R80	1 × 106	80.0	49.7	13.5	11.8	10.4
M5E4R15	5 × 104	15.0	70.7	102.2	4.3	5.4
M1E5R20	1 × 105	20.0	79.6	86.2	4.7	6.6
M1E4R05	1 × 104	5.0	127.3	551.8	1.9	4.1
M1E6R45	1 × 106	45.0	157.2	75.7	5.0	13.8
M1E5R10	1 × 105	10.0	318.3	689.7	1.7	9.3
M1E4R03	1 × 104	3.0	353.7	2554.6	0.9	5.4
M1E6R25	1 × 106	25.0	509.3	441.4	2.1	18.6
M1E4R02	1 × 104	2.0	795.8	8621.6	0.5	6.6
M1E5R05	1 × 105	5.0	1273.2	5517.8	0.6	13.1

Note. Parameters listed are for initial conditions of spherical clouds. Column 1: model name. Column 2: cloud mass (M_⊙). Column 3: cloud radius (pc). Column 4: gas surface density (M_⊙ pc⁻²). Column 5: number density of neutral gas (cm⁻³). Column 6: freefall time (Myr). Column 7: escape velocity at the cloud surface (km s⁻¹). The fiducial model M1E5R20 is shown in bold.

Download table as:  ASCIITypeset image

In addition to the standard runs including both photoionization and radiation pressure (PH+RP), we run two control simulations for each cloud, in which either only photoionization (PH-only) or only radiation pressure (RP-only) is turned on. This enables us to isolate the effect of each feedback mechanism and thus to indirectly assess their relative importance in cloud dispersal. For the fiducial cloud, we additionally run low- and high-resolution models (M1E5R20_N128 and M1E5R20_N512) with N_cell = 128³ and 512³, as well as a no-feedback run (M1E5R20_nofb), in the last of which stellar feedback is turned off. In Appendix A, we compare the results from the models with different resolution. Although star formation completes slightly later in simulations with higher resolution, the overall evolutionary behavior is qualitatively the same, and key quantitative outcomes (such as the SFE) are quite close at different resolution, following a converging trend. All simulations are run over four to seven freefall times, long enough for star formation to complete and for radiation feedback to evacuate the remaining gas from the simulation domain.

Similar to the simulations of Raskutti et al. (2016), the early evolution of our fiducial model is governed by turbulence and self-gravity. The density structure becomes increasingly filamentary as a result of shock compression. The densest regions—which may be at junctions of filaments, or simply overdense clumps within filaments—undergo gravitational collapse, leading to the formation of sink particles. In the fiducial model, the first collapse occurs at ${t}_{* ,0}/{t}_{\mathrm{ff},0}=0.4$.

Due to star formation and the ensuing feedback, the gas in our simulations is ultimately channeled into three distinct states. One portion of the gas gravitationally collapses and is accumulated in sink particles as a total stellar mass ${M}_{* }(t)$. Strong ionizing radiation from newly formed stars ionizes a portion of the initially neutral gas, with the total mass of photoevaporated gas increasing in time as ${M}_{\mathrm{ion}}(t)$. This includes the ionized gas currently in the simulation box and the cumulative ionized gas that has left the simulation domain. Over time, all of the neutral gas in the domain either collapses to make stars, is photoevaporated and ejected, or is ejected while still neutral. We denote the cumulative ejected mass of the neutral gas by ${M}_{\mathrm{ej},\mathrm{neu}}(t)$.⁶ Of course, gas that is photoionized may recombine before it flows out of the box, but at the end of the simulation when no gas remains in the domain, ${M}_{\mathrm{ion},\mathrm{final}}={M}_{\mathrm{ion}}({t}_{\mathrm{final}})$. The condition of mass conservation requires ${M}_{0}={M}_{* ,\mathrm{final}}+{M}_{\mathrm{ion},\mathrm{final}}+{M}_{\mathrm{ej},\mathrm{neu},\mathrm{final}}$ at the end of the run.

Figure 2(a) plots the temporal changes of the total gas mass ${M}_{\mathrm{gas}}$ in the simulation volume, as well as M_*, ${M}_{\mathrm{ion}}$, and ${M}_{\mathrm{ej},\mathrm{neu}}$, all normalized by the initial cloud mass ${M}_{0}$. Figure 2(b) plots the temporal evolution of the volume filling factor of the ionized gas ${\overline{x}}_{{\rm{i}}}=\int ({n}_{{{\rm{H}}}^{+}}/{n}_{{\rm{H}}}){dV}/\int {dV}$, the escape fractions of ionizing (${f}_{\mathrm{esc},{\rm{i}}}$) and non-ionizing (${f}_{\mathrm{esc},{\rm{n}}}$) radiation, and the fraction of ionizing photons absorbed by dust (${f}_{\mathrm{dust},{\rm{i}}}$). These escape fractions and dust absorption fraction here are the instantaneous probabilities of escape or absorption from a single ray-trace at a given time.

Zoom In Zoom Out Reset image size

Figure 3 plots selected snapshots of the fiducial model in the x–y plane. From top to bottom, the rows show column density of neutral gas projected along the z-axis, slices through the stellar center of mass of the neutral hydrogen number density ${n}_{{{\rm{H}}}^{0}}$, electron number density ${n}_{{\rm{e}}}$, ionizing radiation energy density ${{ \mathcal E }}_{{\rm{i}}}$, and non-ionizing radiation energy density ${{ \mathcal E }}_{{\rm{n}}}$. The projected positions of star particles are marked by circles colored by their ages. In the top row, all star particles in the simulation volume are shown, while in the bottom four rows only those within ${\rm{\Delta }}z=\pm 5\,\mathrm{pc}$ of the slice are shown. The selected times (columns from left to right) are $t/{t}_{\mathrm{ff},0}=$ 0.57, 0.80, 1.07, and 1.40 when 10%, 25%, 50%, and 90% of the final stellar mass has been assembled, respectively: here and hereafter, these epochs will be referred to as t_*,10%, t_*,25%, t_*,50%, and t_*,90%, respectively. Figure 4 displays volume rendered images of the fiducial model at these epochs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because of the initial turbulence, the background medium is highly clumpy and filamentary. Young H ii regions formed around star particles quickly break out of their dense natal clumps, and in the larger-scale turbulent cloud, some lines of sight have quite low optical depth. As a result, a substantial fraction of both ionizing and non-ionizing photons escape from the simulation box even at quite early stages of evolution. For example, at t_*,10% (only $0.17{t}_{\mathrm{ff},0}=0.8\,\mathrm{Myr}$ after the first collapse), five star particles have been created and more than 50% of the computational box is filled with ionized gas, although it accounts for only 3% of the total gas mass in the domain. By this time, 13% of the radiation has escaped from the domain overall (with an instantaneous escape probability of 16%). Even though the second quadrant in the x–y plane is rapidly ionized, at early times (persisting through t_*,25%), there are still regions (see lower-left panels of Figure 3) that are fully neutral because they are shadowed by dense gas that absorbs all the ionizing photons on rays emerging from the central star-forming cluster.

At t_*,25%, star formation is concentrated in the central $10\,\mathrm{pc}$ region. Multiple H ii regions expand simultaneously and merge with each other, accompanying a rapid increase of ${M}_{\mathrm{ion}}$ by intense photoevaporation. Due to thermal pressure gradients, the ionized gas is pushed to the outer low-density regions, developing a roughly exponential radial density profile outside the central cavity. Ionized gas starts to escape from the box at mildly supersonic outflow velocities of ∼20 km s⁻¹. Simultaneously, the remaining neutral gas filaments are pushed away from the central collection of stars by a combination of radiation pressure forces and the rocket effect as gas photoevaporates from their surfaces (see top two rows of Figure 3). At later times, star formation occurs mainly in dense clumps within structures that have been pushed to the periphery of the combined H ii region. We note that the expansion over time of the loci of star formation is not primarily due to shock-induced collapse as the H ii region expands (collapse in filaments would have occurred anyway), but because dense gas is progressively evacuated from a larger and larger volume.

Figure 5 plots the maps of the emission measure (EM) $\int {n}_{e}^{2}d{\ell }$ integrated along the z-direction of the fiducial model at t_*,10%, t_*,25%, t_*,50%, and t_*,90%. In the early phase of star formation, the EM map is dominated by bright, compact H ii regions around newly formed stars, while at later times it becomes rich in substructures resembling observed pillars and bright-rimmed globules (e.g., Koenig et al. 2012; Hartigan et al. 2015); these features are also apparent in the n_e slices of Figure 3. Similar structures have also been seen in previous numerical simulations of expanding H ii regions in a turbulent medium (e.g., Mellema et al. 2006; Gritschneder et al. 2009, 2010; Walch et al. 2012). The globules exhibit ionization-driven ablation flows and distinctive elongated tails pointing away from the ionizing sources. Anisotropic mass loss, preferentially from the sides of globules facing the star particles, exerts a recoil force on them; the globules rocket away from the stellar sources at a typical radial velocity of v_ej,neu ∼ 5 km s⁻¹.

Zoom In Zoom Out Reset image size

By t_*,90%, not only has the H ii region engulfed the entire cloud, but the ionized region has expanded to twice the original cloud radius. Furthermore, the central $10\,\mathrm{pc}$ has mostly been cleared of dense gas. The escape fractions of ionizing and non-ionizing photons keep increasing with time as neutral gas surrounding ionizing sources is pushed away (see also, e.g., Rogers & Pittard 2013; Raskutti et al. 2017).

The cloud has been completely destroyed by $t/{t}_{\mathrm{ff},0}\sim 2$, with the fraction of the remaining neutral gas <2% of the initial cloud mass. The final SFE is ${\varepsilon }_{* }\equiv {M}_{* ,\mathrm{final}}/{M}_{0}=0.13$ in the fiducial model. Only 7% of the initial cloud mass is ejected in the neutral phase, while 81% is lost to photoevaporation. The remaining cluster of star particles is loosely bound gravitationally and dissolves rapidly; 30% of the final stellar mass leaves the computational domain from $t/{t}_{\mathrm{ff},0}=2.09$ to $t/{t}_{\mathrm{ff},0}=4$.

In our simulations, the SFE of a cloud rises until the associated feedback photoevaporates and dynamically ejects the remaining gas, completely halting further star formation. We define the final net SFE as ${\varepsilon }_{* }\equiv {M}_{* ,\mathrm{final}}/{M}_{0}$. In addition to the net SFE, we similarly define photoevaporation and ejection efficiencies as ${\varepsilon }_{\mathrm{ion}}\equiv {M}_{\mathrm{ion},\mathrm{final}}/{M}_{0}$ and ${\varepsilon }_{\mathrm{ej}}\equiv {M}_{\mathrm{ej},\mathrm{final}}/{M}_{0}$, respectively, the latter of which counts both neutral and ionized gases. Note that these quantities are related by mass conservation through ${\varepsilon }_{\mathrm{ej}}=1-{\varepsilon }_{* }={\varepsilon }_{\mathrm{ion}}+{\varepsilon }_{\mathrm{ej},\mathrm{neu}}$, where ${\varepsilon }_{\mathrm{ej},\mathrm{neu}}\equiv {M}_{\mathrm{ej},\mathrm{neu},\mathrm{final}}/{M}_{0}$ and ${\varepsilon }_{\mathrm{ion}}={\varepsilon }_{\mathrm{ej},\mathrm{ion}}$, since no ionized gas is left inside the box at the end of the runs.

Figure 6 plots as star symbols (a) the net SFE ${\varepsilon }_{* }$, (b) the ejection efficiency ${\varepsilon }_{\mathrm{ej}}$, (c) the photoevaporation efficiency ${\varepsilon }_{\mathrm{ion}}$, and (d) the neutral ejection efficiency ${\varepsilon }_{\mathrm{ej},\mathrm{neu}}$ for all of our PH+RP runs as functions of the initial surface density; ${\varepsilon }_{* }$, ${\varepsilon }_{\mathrm{ion}}$, and ${\varepsilon }_{\mathrm{ej},\mathrm{neu}}$ are also tabulated in columns 2–4 of Table 2. The size of the symbols represents the initial cloud mass, a convention that will be adopted throughout this paper. The net SFE depends weakly on ${M}_{0}$ but strongly on ${{\rm{\Sigma }}}_{0}$, increasing from 0.04 for M1E5R50 to 0.51 for M1E5R05. For clouds with low ${{\rm{\Sigma }}}_{0}$ and intermediate-to-high ${M}_{0}$ (M5E4R15, M1E5R20, M1E5R30, M1E5R40, M1E5R50, M1E6R45, M1E6R80), photoevaporation is effective in destroying clouds: ${\varepsilon }_{\mathrm{ion}}\gtrsim 0.7$ and the resulting SFE is ${\varepsilon }_{* }\lesssim 0.2$. By contrast, the low mass clouds with ${M}_{0}={10}^{4}\,{M}_{\odot }$ (M1E4R02, M1E4R03, M1E4R05, M1E4R08) and/or high surface density (M1E5R10, M1E5R05, M1E6R25) clouds have ${\varepsilon }_{\mathrm{ej},\mathrm{neu}}\gtrsim 0.1$, implying that dynamical ejection of neutral gas is non-negligible for disruption and quenching of further star formation in these clouds.

Zoom In Zoom Out Reset image size

Table 2.  Simulation Results

Model	${\varepsilon }_{* }$	${\varepsilon }_{\mathrm{ion}}$	${\varepsilon }_{\mathrm{ej},\mathrm{neu}}$	t*,0	tSF	${t}_{\mathrm{dest}}$	${t}_{\mathrm{ion}}$	${v}_{\mathrm{ej},\mathrm{neu}}$	${v}_{\mathrm{ej},\mathrm{ion}}$	$\langle {f}_{\mathrm{ion}}\rangle$	$\langle q\rangle$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
M1E5R50	0.04	0.95	0.02	0.52	0.41	0.56	0.60	8.3	24.4	0.19	61
M1E5R40	0.05	0.93	0.03	0.46	1.03	0.78	0.75	9.4	24.3	0.22	81
M1E5R30	0.08	0.89	0.04	0.44	0.72	0.95	0.98	7.8	23.9	0.24	116
M1E4R08	0.04	0.52	0.45	0.56	1.24	1.04	0.96	6.3	18.1	0.42	69
M1E6R80	0.10	0.89	0.02	0.33	0.62	0.69	0.74	10.7	26.3	0.37	320
M5E4R15	0.12	0.79	0.09	0.44	1.14	1.45	1.31	7.1	23.1	0.26	148
M1E5R20	0.13	0.81	0.07	0.40	1.15	1.39	1.32	7.9	24.1	0.26	203
M1E4R05	0.13	0.56	0.32	0.52	1.19	2.07	1.61	6.8	20.2	0.36	199
M1E6R45	0.22	0.73	0.06	0.27	1.12	1.37	1.44	13.2	28.6	0.37	753
M1E5R10	0.28	0.56	0.17	0.33	1.79	2.26	2.08	8.3	25.6	0.42	807
M1E4R03	0.20	0.47	0.34	0.44	2.13	3.16	2.48	6.0	20.9	0.36	402
M1E6R25	0.38	0.52	0.11	0.24	1.70	2.41	2.38	14.7	35.6	0.43	2032
M1E4R02	0.40	0.32	0.29	0.40	2.32	4.10	2.97	6.7	24.1	0.46	1630
M1E5R05	0.51	0.31	0.19	0.27	2.35	3.59	3.46	9.7	30.2	0.50	2950
M1E5R20_N128	0.15	0.79	0.07	0.42	0.92	1.13	1.26	8.7	23.3	0.25	261
M1E5R20_N512	0.12	0.82	0.07	0.38	1.35	1.53	1.39	7.6	24.2	0.26	188
M1E5R20_nofb	0.88	⋯	0.13	0.40	5.85	⋯	⋯	3.9	⋯	⋯	⋯

Note. Column 1: model name. Column 2: net star formation efficiency. Column 3: photoevaporation efficiency. Column 4: neutral ejection efficiency. Column 5: time of first star formation (in units of ${t}_{\mathrm{ff},0}$). Column 6: star formation duration t_SF = t_*,95%–t_*,0 (units ${t}_{\mathrm{ff},0}$). Column 7: timescale for cloud destruction t_dest = t_neu,5%–t_*,0 (units ${t}_{\mathrm{ff},0}$). Column 8: timescale for photoevaporation (units ${t}_{\mathrm{ff},0};$ Equation (12)). Column 9: time-averaged outflow velocity of the neutral gas (km s⁻¹). Column 10: time-averaged outflow velocity of the ionized gas (km s⁻¹). Column 11: time-averaged hydrogen absorption fraction (see Equation (8)). Column 12: time-averaged shielding factor. The fiducial model M1E5R20 is shown in bold.

Download table as:  ASCIITypeset image

The effect of radiation feedback on cloud disruption is limited in the highest-${{\rm{\Sigma }}}_{0}$ clouds, such as M1E4R02, M1E5R05, and M1E6R25. These reach ${\varepsilon }_{* }\gtrsim 0.4$ before the remaining gas is dispersed. This is due to a deep gravitational potential well, which outflowing gases must overcome to escape from the system (e.g., Dale et al. 2012). These clouds also have a comparatively dense and thick layer of recombining gas and dust that absorbs most ionizing photons, resulting in relatively inefficient photoevaporative mass loss (see Section 4.1). Thus radiation pressure plays a relatively more important role for the highest surface density clouds.

One can assess the relative (negative) impact of photoionization and radiation pressure feedback on star formation by comparing the net SFEs of the PH-only or RP-only runs with those of PH+RH models. Figure 7 plots ${\varepsilon }_{* }$ of PH-only (red circles), RP-only (blue squares), and PH+RP (dark stars) models as a function of ${{\rm{\Sigma }}}_{0}$. Overall, the increasing trend of ${\varepsilon }_{* }$ with ${{\rm{\Sigma }}}_{0}$ is similar for all models. For most of our models, the net SFE of the PH-only run is smaller than that of the RP-only counterpart and closer to that of the PH+RP run, indicating that photoionization feedback plays a more important role than radiation pressure in disrupting the parent clouds. The two exceptions are the massive, high surface density clouds (M1E5R05 and M1E6R25), for which the net SFE of the RP-only run is smaller than that of the PH-only run.

Zoom In Zoom Out Reset image size

The net SFEs of the RP-only runs lie approximately on a single sequence as a function of ${{\rm{\Sigma }}}_{0}$, regardless of ${M}_{0}$, also in accordance with the recent theoretical predictions of the regulation of star formation by the radiation pressure feedback alone (Raskutti et al. 2016; Thompson & Krumholz 2016). We quantitatively show in Appendix B that our RP-only results are indeed in good agreement with the analytic predictions. Our net SFE is slightly lower than those obtained in the numerical simulations of Raskutti et al. (2016), which employed a grid-based moment method for radiation. Paper I showed that this method does not adequately resolve the radiation field near individual star particles (because radiation sources have finite size), leading to more gas accretion and thus a higher SFE.

Figure 7 also shows that model M1E5R20_nofb (green triangle) in the absence of feedback converts 87% of the gas into stars over the course of evolution. The reason that even this no-feedback run has a nonzero ejection efficiency of ${\varepsilon }_{\mathrm{ej}}=0.13$ is that initial turbulence makes a small fraction of the gas gravitationally unbound. Raskutti et al. (2016) found that the initial "turbulent outflow ejection" efficiency is ε_ej,turb = 10%–15% for all models (see their Figure 17(a)). In principle, this can be used to derive an "adjusted" feedback-induced SFE of ${\varepsilon }_{* ,\mathrm{adj}}={\varepsilon }_{* }/(1-{\varepsilon }_{\mathrm{ej},\mathrm{turb}})$.

Finally, we note that the feedback effects of photoionizing radiation and non-ionizing radiation on quenching star formation are not simply additive. If feedback effects were additive, the naive expectation would be $1/{\varepsilon }_{* ,\mathrm{PH}+\mathrm{RP}}=1/{\varepsilon }_{* ,\mathrm{PH}}+1/{\varepsilon }_{* ,\mathrm{RP}}$. However, we find that the true value of ε_*,PH+RP is 20%–30% larger than this naive expectation.

The cloud lifetime is important to understanding the life cycle of GMCs within the context of other ISM phases, yet observational constraints on it remain quite uncertain (Heyer & Dame 2015). We use our simulations to directly measure the relevant timescales as follows. We define the star formation timescale in our models as the time taken to assemble 95% of the final stellar mass after the onset of radiation feedback (i.e., t_SF ≡ t_*,95%–t_*,0). The approximate gas depletion timescale by star formation is then ${t}_{\mathrm{dep}}={t}_{\mathrm{SF}}/{\varepsilon }_{* }$. We regard a cloud as being effectively destroyed at t = t_neu,5% when only 5% of the initial cloud mass is left over as the neutral phase in the simulation box, with the other 95% transformed into stars, photoionized, or ejected. We define the destruction timescale as ${t}_{\mathrm{dest}}={t}_{\mathrm{neu},5 \% }-{t}_{* ,0}$.⁷ In most cases, ${t}_{\mathrm{dest}}\gt {t}_{\mathrm{SF}}$,⁸ but generally the two timescales are very similar. Columns 5–7 of Table 2 list t_*,0, t_SF, and ${t}_{\mathrm{dest}}$, respectively.

Figure 8 plots t_SF (triangles), ${t}_{\mathrm{dest}}$ (squares), and ${t}_{\mathrm{dep}}$ (circles) of our PH+RP runs in units of Myr (left panel) and in units of the initial freefall time (right panel). In physical units, t_SF and ${t}_{\mathrm{dest}}$ increase from ∼1–$3\,\mathrm{Myr}$ in the densest clouds (${n}_{{\rm{H}},0}\gtrsim 500\ {\mathrm{cm}}^{-3}$ and ${{\rm{\Sigma }}}_{0}\gtrsim 300\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$) to ∼8–$14\,\mathrm{Myr}$ in the lowest-density clouds (${n}_{{\rm{H}},0}\lesssim 10\ {\mathrm{cm}}^{-3}$ and ${{\rm{\Sigma }}}_{0}\lesssim 20\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$). Relative to the freefall time, t_SF (${t}_{\mathrm{dest}}$) decreases from $\sim 2{t}_{\mathrm{ff},0}$ ($\sim 4{t}_{\mathrm{ff},0}$) in the dense clouds to $\sim 0.4{t}_{\mathrm{ff},0}$ ($\sim 0.6{t}_{\mathrm{ff},0}$) in the low-density clouds. Our result that ${t}_{\mathrm{SF}}\lesssim 2{t}_{\mathrm{ff},0}$ is consistent with the picture of rapid star formation envisaged by Elmegreen (2000). The depletion timescale decreases considerably with ${{\rm{\Sigma }}}_{0}$ from $250\,\mathrm{Myr}$ for low-Σ₀ model (M1E5R40) to $2.7\,\mathrm{Myr}$ for high-Σ₀ models (M1E4R02 and M1E5R05).

Zoom In Zoom Out Reset image size

In the presence of ionizing radiation, the outflows in our simulations consist of both neutral and ionized gas, and it is interesting to measure the dependence of ejection velocities on cloud properties. We calculate the mass-weighted ejection velocities of neutral/ionized gas as

$$\begin{eqnarray}&&{v}_{\mathrm{ej},\mathrm{neu}/\mathrm{ion}}=\displaystyle \frac{{p}_{\mathrm{ej},\mathrm{neu}/\mathrm{ion},\mathrm{final}}}{{M}_{\mathrm{ej},\mathrm{neu}/\mathrm{ion},\mathrm{final}}},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{p}_{\mathrm{ej},\mathrm{neu}/\mathrm{ion},\mathrm{final}}=\int {dt}{\int }_{\partial V}{dA}\,{\rho }_{\mathrm{neu}/\mathrm{ion}}({\boldsymbol{v}}\cdot \hat{{\boldsymbol{n}}})({\boldsymbol{v}}\cdot \hat{{\boldsymbol{r}}})\end{eqnarray} \tag{ 2 }$$

is the time-integrated total radial momentum of neutral/ionized outflowing gas, with $\hat{{\boldsymbol{r}}}$ being the unit radial vector with respect to the stellar center of mass and $\hat{{\boldsymbol{n}}}$ being the unit vector normal to the surface area dA at the outer boundary of the box. These values are given in columns 9 and 10 of Table 2.

Figure 9 plots ${v}_{\mathrm{ej},\mathrm{neu}}$ (triangles) and ${v}_{\mathrm{ej},\mathrm{ion}}$ (circles) in units of $\mathrm{km}\,{{\rm{s}}}^{-1}$ (left panel) and in units of the initial escape velocity ${v}_{\mathrm{esc},0}$ (right panel) for all the PH+RP runs. The outflow velocity of the ionized gas is ∼18–36 km s⁻¹ and larger than ${v}_{\mathrm{esc},0}$ across the whole range of ${{\rm{\Sigma }}}_{0}$. In low surface density clouds, these supersonic outflows of the ionized gas are driven primarily by thermal pressure. Since ionized gas has low optical depth and hence high Eddington ratios, its ejection in high surface density clouds is further helped by radiation pressure, while strong gravity tends to reduce the outflow velocity.

Zoom In Zoom Out Reset image size

The neutral gas is ejected at a typical velocity of ∼6–15 km s⁻¹, about two to four times smaller than the ejection velocity of the ionized gas. This is roughly consistent with the characteristic rocket velocity ∼5–10 km s⁻¹ of cometary globules in a photoevaporation-dominated medium (e.g., Bertoldi & McKee 1990), although the outflow velocity is also affected by gravity as well as radiation pressure forces, especially for high surface density clouds. Note that ${v}_{\mathrm{ej},\mathrm{neu}}/{v}_{\mathrm{esc},0}\sim 0.7$–2.0 decreases slowly with increasing ${{\rm{\Sigma }}}_{0}$, owing to gravity.

We first consider the case in which star formation is quenched solely by photoevaporation. Ionizing photons emitted by star particles arrive at IFs after passing through an intervening volume in which the gas is nearly fully ionized, with a balance between ionization and recombination. Ionized gas created at IFs does not participate in star formation (provided it remains sufficiently exposed to radiation) and will be eventually pushed out of the cloud and simulation volume by thermal and radiation pressures. Here we use physical scaling arguments to estimate the net photoevaporation rate in a star-forming cloud.

4.1.1. Ionization Analysis

H ii regions formed in our simulations are highly irregular in shape. The EM map is dominated by pillars, ridges, and cometary globules near the H ii region peripheries, as Figure 5 illustrates. Each of these bright features marks the IF on the surface of an optically thick, neutral gas structure that is facing the ionizing sources. These neutral structures are photoevaporated as the IF advances. Since the newly ionized gas can expand into the surroundings, the IFs are typically D-critical or of strong-D type, and the ionized gas streams off the IF at approximately the sound speed and forms an ionized boundary layer (e.g., Oort & Spitzer 1955; Kahn 1969; Elmergreen 1976; Bertoldi 1989). In addition to the bright IF regions, Figures 3–5 show that there is substantial diffuse ionized gas in the volume, filling the space between and beyond the filaments and clumps of neutral gas. Along some lines of sight from stellar sources, ionizing photons moving toward optically thin regions produce a weak-R type IF that quickly propagates out of the simulation domain, such that in those directions the H ii region is effectively density bounded. In other directions, lines of sight from ionizing sources propagate through diffuse ionized gas until they end at an IF on the surface of dense neutral structure. As we will see below, almost all of the ionizations and recombinations in our simulations occur in the ionization-bounded Strömgren volume (IBSV), and ionization and recombination are in global equilibrium.

The equation for net electron production integrated over the entire domain reads

$$\begin{eqnarray}&&{\dot{N}}_{{\rm{e}}}=\int { \mathcal I }\,{dV}-\int { \mathcal R }\,{dV},\end{eqnarray} \tag{ 3 }$$

where ${ \mathcal I }={n}_{{{\rm{H}}}^{0}}(c{{ \mathcal E }}_{{\rm{i}}}/h{\nu }_{{\rm{i}}}){\sigma }_{{{\rm{H}}}^{0}}$ is the photoionization rate per unit volume, with c being the speed of light, and ${ \mathcal R }={\alpha }_{{\rm{B}}}{n}_{{\rm{e}}}{n}_{{{\rm{H}}}^{+}}$ is the recombination rate per unit volume, with ${\alpha }_{{\rm{B}}}=3.03\,\times {10}^{-13}{(T/8000{\rm{K}})}^{-0.7}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}$ being the case B recombination coefficient. The rate at which ionized gas is newly created within the domain (from ionizations exceeding recombinations) is ${\dot{M}}_{\mathrm{ion}}={\mu }_{{\rm{H}}}{\dot{N}}_{{\rm{e}}}$, where ${\mu }_{{\rm{H}}}=1.4{m}_{{\rm{H}}}$ is the mean atomic mass per hydrogen. Part of this may go into an increase in the mass of ionized gas within the box, and part may go into escape of ionized gas from the domain: ${\dot{M}}_{\mathrm{ion}}={\dot{M}}_{\mathrm{ion},\mathrm{box}}+{\dot{M}}_{\mathrm{ion},\mathrm{ej}}$.

In our simulations, some fraction of the ionizing photons emitted by the stars are absorbed by neutral hydrogen, while others are absorbed by dust or escape from the simulation domain. Let ${f}_{\mathrm{ion}}=1-{f}_{\mathrm{esc},{\rm{i}}}-{f}_{\mathrm{dust},{\rm{i}}}$ denote the fraction of ionizing photons absorbed by neutral hydrogen. Then, the total photoionization rate of neutral hydrogen inside the H ii region can be written as

$$\begin{eqnarray}&&\int { \mathcal I }\,{dV}={f}_{\mathrm{ion}}{Q}_{{\rm{i}}}\equiv {Q}_{{\rm{i}},\mathrm{eff}}.\end{eqnarray} \tag{ 4 }$$

Newly ionized gas is created at IFs, and we denote the total effective area of these IFs as ${A}_{{\rm{i}}}$. Gas streams away from the IFs with typical number density ${n}_{{\rm{i}}}$ and characteristic velocity ${c}_{{\rm{i}}}={(2.1{k}_{{\rm{B}}}{T}_{\mathrm{ion}}/{\mu }_{{\rm{H}}})}^{1/2}=10.0\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (e.g., Bertoldi 1989). We thus write

$$\begin{eqnarray}&&{\dot{N}}_{{\rm{e}}}\equiv {c}_{{\rm{i}}}{n}_{{\rm{i}}}{A}_{{\rm{i}}}\,\end{eqnarray} \tag{ 5 }$$

for the rate of election production at the IFs. Equation (5) may be thought of as defining the product ${n}_{{\rm{i}}}{A}_{{\rm{i}}}$ in terms of ${\dot{N}}_{{\rm{e}}}$.

To reach the IFs, photons traverse a volume of nearly fully ionized gas that is undergoing recombination (which is approximately balanced by ionization; see below). If we define the effective path length through this ionized volume as ${H}_{{\rm{i}}}$, we can write the total recombination rate as

$$\begin{eqnarray}&&\int { \mathcal R }\,{dV}\equiv {\alpha }_{{\rm{B}}}{n}_{{\rm{i}}}^{2}{A}_{{\rm{i}}}{H}_{{\rm{i}}}.\end{eqnarray} \tag{ 6 }$$

With ${n}_{{\rm{i}}}{A}_{{\rm{i}}}$ defined via Equation (5), we can think of Equation (6) as defining the ratio ${A}_{{\rm{i}}}/{H}_{{\rm{i}}}$, which has units of length.

The volume of highly ionized gas between the source and the IF acts as an insulator, absorbing photons that would otherwise reach the IF in order to balance recombination. We define the ratio of available ionizing photons to photons that actually reach the IF as a shielding factor:

$$\begin{eqnarray}&&q\equiv \displaystyle \frac{{Q}_{{\rm{i}},\mathrm{eff}}}{{\dot{N}}_{{\rm{e}}}}=1+\displaystyle \frac{\int { \mathcal R }\,{dV}}{\int { \mathcal I }\,{dV}-\int { \mathcal R }\,{dV}}=1+\displaystyle \frac{{\alpha }_{{\rm{B}}}{n}_{{\rm{i}}}{H}_{{\rm{i}}}}{{c}_{{\rm{i}}}}.\end{eqnarray} \tag{ 7 }$$

With this definition, q = 1 would imply all the ionizing photons are used to ionize neutrals at the IF, while $q\gg 1$ when most of the ionizing photons are shielded by the IBSV. The second term on the right-hand side of Equation (7), ${\alpha }_{{\rm{B}}}{n}_{{\rm{i}}}{H}_{{\rm{i}}}/{c}_{{\rm{i}}}$, represents the ratio of the characteristic flow timescale to the recombination timescale, or the number of recombinations occurring over the time it would take the flow to cross the IBSV.

Figure 10 plots the temporal histories of (a) the specific evaporation rate ${\mu }_{{\rm{H}}}{\dot{N}}_{{\rm{e}}}/{M}_{0}$, (b) the total photoionization rate ${Q}_{{\rm{i}},\mathrm{eff}}$, (c) the hydrogen absorption fraction of ionizing photons ${f}_{\mathrm{ion}}$, and (d) the shielding factor q, for all of our PH+RP runs. Both the evaporation rate and photoionization rate reach peak values and then decline with time as the cloud is destroyed.

Zoom In Zoom Out Reset image size

The hydrogen absorption fraction is largest in the early, embedded phase of star formation and decreases with time as optical depth drops and most of the ionizing photons escape the computational box. Column (11) of Table 2 gives

$$\begin{eqnarray}&&\langle {f}_{\mathrm{ion}}\rangle \equiv \displaystyle \frac{\int {\dot{N}}_{{\rm{e}}}{f}_{\mathrm{ion}}\,{dt}}{\int {\dot{N}}_{{\rm{e}}}\,{dt}},\end{eqnarray} \tag{ 8 }$$

where the angle brackets $\langle \,\rangle$ denote the ${\dot{N}}_{{\rm{e}}}$-weighted temporal average over the whole simulation period after t_*,0 (when ${\dot{N}}_{{\rm{e}}}$ is nonzero). The value of $\langle {f}_{\mathrm{ion}}\rangle$ ranges from 0.19 to 0.50. Overall, a larger fraction of ionizing radiation is absorbed by hydrogen in higher surface density clouds as H ii regions undergo a relatively longer embedded phase.

Figure 10 shows that the shielding factor is much larger than unity for all models at all times. Column (12) of Table 2 gives ${\dot{N}}_{{\rm{e}}}$-weighted time-averaged values $\langle q\rangle$, which range from 61 to 2950, increasing with mass and surface density. This is consistent with the expectation (cf. the right-hand side of Equation (7)) that a higher recombination rate in denser clouds, as well as longer path lengths in larger clouds, increases shielding and therefore reduces the efficiency of photoevaporation.

The fact that our clouds have $q\gg 1$ justifies the assumption of a global equilibrium between ionization and recombination:

$$\begin{eqnarray}&&{Q}_{{\rm{i}},\mathrm{eff}}=\int { \mathcal I }\,{dV}\approx \int { \mathcal R }\,{dV}={\alpha }_{{\rm{B}}}{n}_{{\rm{i}}}^{2}{A}_{{\rm{i}}}{H}_{{\rm{i}}}.\end{eqnarray} \tag{ 9 }$$

We solve this to obtain ${n}_{{\rm{i}}}\approx {[{Q}_{{\rm{i}},\mathrm{eff}}/({\alpha }_{{\rm{B}}}{A}_{{\rm{i}}}{H}_{{\rm{i}}})]}^{1/2}$, and from Equation (7) we obtain the shielding factor

$$\begin{eqnarray}&&q\approx \displaystyle \frac{{\alpha }_{{\rm{B}}}{n}_{{\rm{i}}}{H}_{{\rm{i}}}}{{c}_{{\rm{i}}}}\approx \displaystyle \frac{{\alpha }_{{\rm{B}}}^{1/2}}{{c}_{{\rm{i}}}}{\left(\displaystyle \frac{{Q}_{{\rm{i}},\mathrm{eff}}{H}_{{\rm{i}}}}{{A}_{{\rm{i}}}}\right)}^{1/2},\end{eqnarray} \tag{ 10 }$$

while Equation (5) yields an approximate expression for the evaporation rate:

$$\begin{eqnarray}&&{\dot{N}}_{{\rm{e}}}\approx \displaystyle \frac{{c}_{{\rm{i}}}}{{\alpha }_{{\rm{B}}}^{1/2}}{\left(\displaystyle \frac{{Q}_{{\rm{i}},\mathrm{eff}}{A}_{{\rm{i}}}}{{H}_{{\rm{i}}}}\right)}^{1/2}.\end{eqnarray} \tag{ 11 }$$

The timescale for the photoevaporative mass loss tends to be longer in high surface density clouds, similar to ${t}_{\mathrm{dest}}$ and t_SF. Despite the fact that they have a relatively high hydrogen absorption fraction and a small escape fraction of ionizing photons, the mass loss is inefficient because of large q. In addition, the motions of sink particles and the associated neutral envelopes cause a sudden change in the optical depth of ionizing photons. As a result, the size and shape of the H ii region and photoevaporation rate fluctuate with time until they completely exhaust accretion flows and break out to larger radii at late time (e.g., Peters et al. 2010; Dale et al. 2012).

We shall define the cloud photoevaporation timescale as

$$\begin{eqnarray}&&{t}_{\mathrm{ion}}\equiv \displaystyle \frac{{N}_{{\rm{e}}}}{\langle {\dot{N}}_{{\rm{e}}}\rangle }=\displaystyle \frac{{N}_{{\rm{e}}}^{2}}{\int {\dot{N}}_{{\rm{e}}}^{2}\,{dt}},\end{eqnarray} \tag{ 12 }$$

where ${N}_{{\rm{e}}}=\int {\dot{N}}_{{\rm{e}}}\,{dt}$ is the integrated number of electrons. Column (8) of Table 2 lists the value of ${t}_{\mathrm{ion}}$ for all models. Our numerical simulations show that this timescale is typically comparable to the initial freefall timescale for the cloud, ${t}_{\mathrm{ff},0}$.

4.1.2. Constraints from Simulations

Once ${Q}_{{\rm{i}},\mathrm{eff}}$, ${A}_{{\rm{i}}}$, and ${H}_{{\rm{i}}}$ are known, Equation (11) can be used to estimate the photoevaporation mass loss rate in star-forming clouds. It is natural to expect that the total area ${A}_{{\rm{i}}}$ of the IF scales with the surface area $4\pi {R}_{0}^{2}$ of the cloud, while the thickness ${H}_{{\rm{i}}}$ of the IBSV scales with the cloud radius R₀. Our numerical results, based on measured ${Q}_{{\rm{i}},\mathrm{eff}}$, ${\dot{N}}_{{\rm{e}}}$, and $\int {n}_{e}^{2}{dV}/\int {n}_{e}{dV}\to {n}_{{\rm{i}}}$, show that time-averaged values of ${A}_{{\rm{i}}}/(4\pi {R}_{0}^{2})={\dot{N}}_{{\rm{e}}}/(4\pi {n}_{{\rm{i}}}{c}_{{\rm{i}}}{R}_{0}^{2})$ and ${H}_{{\rm{i}}}/{R}_{0}\approx {c}_{{\rm{i}}}{Q}_{{\rm{i}},\mathrm{eff}}/({\alpha }_{{\rm{B}}}{n}_{{\rm{i}}}{\dot{N}}_{{\rm{e}}}{R}_{0})$ are indeed of order unity. In detail, the time-averaged values of ${A}_{{\rm{i}}}/(4\pi {R}_{0}^{2})$ and ${H}_{{\rm{i}}}/{R}_{0}$ are in the ranges 0.5–2.0 and 0.6–1.9, respectively. Motivated by Equation (11), together with the characteristic dynamical time of clouds, we introduce the dimensionless evaporation rate ${\phi }_{\mathrm{ion}}$ and evaporation timescale ${\phi }_{t}$, defined as

$$\begin{eqnarray}&&{\phi }_{\mathrm{ion}}\equiv \langle {\dot{N}}_{{\rm{e}}}\rangle \times \displaystyle \frac{{\alpha }_{{\rm{B}}}^{1/2}}{{c}_{{\rm{i}}}{Q}_{{\rm{i}},\max }^{1/2}{R}_{0}^{1/2}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\phi }_{t}\equiv \displaystyle \frac{{t}_{\mathrm{ion}}}{{t}_{\mathrm{ff},0}},\end{eqnarray} \tag{ 14 }$$

where ${Q}_{{\rm{i}},\max }$ is the maximum rate of the total ionizing photons over the lifetime of a star cluster (which is directly related to the SFE).

Figure 11 plots ${\phi }_{\mathrm{ion}}$ and ${\phi }_{t}$ from our PH+RP models as star symbols. Note that ${\phi }_{\mathrm{ion}}\sim 0.79$–1.34, roughly independent of the cloud mass and surface density. From the definition in Equation (13), an order-unity value of ${\phi }_{\mathrm{ion}}$ implies a mass loss rate from clouds in which the characteristic velocity is the ionized gas sound speed, the characteristic spatial scale is the size of the cloud, and the characteristic density is consistent with ionization-recombination equilibrium; comparison to Equation (11) implies that $\langle ({Q}_{{\rm{i}},\mathrm{eff}}/{Q}_{{\rm{i}},\max })$ $({A}_{{\rm{i}}}/{H}_{{\rm{i}}}){\rangle }^{1/2}/{R}_{0}^{1/2}\approx {\phi }_{\mathrm{ion}}\sim 1$.

Zoom In Zoom Out Reset image size

Figure 11 shows that ${\phi }_{t}$ ranges between 0.6 and 3.5 and increases with ${{\rm{\Sigma }}}_{0}$. A typical photoevaporation timescale is thus twice the freefall time in the cloud. Plotted as open squares is the normalized destruction timescale ${t}_{\mathrm{dest}}/{t}_{\mathrm{ff},0}$. This is very close to ${\phi }_{t}$ for low-density, massive clouds, suggesting that these clouds are destroyed primarily by photoevaporation. Low mass and/or high surface density clouds have a destruction timescale somewhat longer than the evaporation timescale, suggesting that cloud destruction involves ejection of neutrals (accelerated by a combination of rocket effect and radiation pressure) as well as ions.

The total number of electrons (and hence ions) created by photoevaporation over the lifetime of the cloud is ${N}_{{\rm{e}}}=\int {\dot{N}}_{{\rm{e}}}{dt}\equiv \langle {\dot{N}}_{{\rm{e}}}\rangle {t}_{\mathrm{ion}}$. The total mass photoevaporated from the cloud can therefore be expressed as

$$\begin{eqnarray}&&{M}_{\mathrm{ion}}={\mu }_{{\rm{H}}}{N}_{e}={\phi }_{t}{\phi }_{\mathrm{ion}}\displaystyle \frac{{\mu }_{{\rm{H}}}{c}_{{\rm{i}}}{Q}_{{\rm{i}},\max }^{1/2}{R}_{0}^{1/2}{t}_{\mathrm{ff},0}}{{\alpha }_{{\rm{B}}}^{1/2}}.\end{eqnarray} \tag{ 15 }$$

Since the right-hand side is proportional to ${M}_{* ,\mathrm{final}}^{1/2}$, this shows that the fraction of mass photoevaporated over the cloud lifetime, ${M}_{\mathrm{ion}}/{M}_{0}$, scales as the square root (rather than linearly) of the SFE. We discuss this result further in Section 5.1. Figure 11(c) plots the product ${\phi }_{t}{\phi }_{\mathrm{ion}}$ as a function of ${{\rm{\Sigma }}}_{0}$. We fit the numerical results as

$$\begin{eqnarray}&&{\phi }_{t}{\phi }_{\mathrm{ion}}={c}_{1}+{c}_{2}{\mathrm{log}}_{10}({S}_{0}+{c}_{3}),\end{eqnarray} \tag{ 16 }$$

where ${c}_{1}=-2.89$, c₂ = 2.11, c₃ = 25.3, and ${S}_{0}={{\rm{\Sigma }}}_{0}/({M}_{\odot }\,{\mathrm{pc}}^{-2})$, which is plotted as the red dashed line in Figure 11(c).

In our simulations, both neutral and ionized gas that does not collapse to make stars is eventually ejected, carrying both mass and momentum out of the cloud. In this section, we characterize the efficiency of the radial momentum injection by computing net momentum yields. We also separately assess the mean thermal and radiation pressure forces in the radial direction, where

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{\mathrm{thm}}=-{\rm{\nabla }}P,\quad {{\boldsymbol{f}}}_{\mathrm{rad}}=\displaystyle \frac{{n}_{{{\rm{H}}}^{0}}{\sigma }_{{{\rm{H}}}^{0}}}{c}{{\boldsymbol{F}}}_{{\rm{i}}}+\displaystyle \frac{{n}_{{\rm{H}}}{\sigma }_{{\rm{d}}}}{c}({{\boldsymbol{F}}}_{{\rm{i}}}+{{\boldsymbol{F}}}_{{\rm{n}}})\end{eqnarray} \tag{ 17 }$$

are the thermal and radiation pressure forces per unit volume, for ${{\boldsymbol{F}}}_{{\rm{i}}}$ and ${{\boldsymbol{F}}}_{{\rm{n}}}$ the ionizing and non-ionizing radiation fluxes, respectively.

4.2.1. Net Momentum Yield

When clouds are disrupted by feedback, the resulting outflow is roughly spherical, so it is useful to measure the total radial momentum of the ejected material induced by feedback. In Section 5, we will use the momentum ejected per unit mass of stars formed, ${p}_{* }/{m}_{* }={p}_{\mathrm{ej},\mathrm{final}}/{M}_{* ,\mathrm{final}}$, which here we term the net "momentum yield."

Figure 12 plots the ${p}_{* }/{m}_{* }$ in the PH-only (circles), RP-only (squares), and PH+RP (stars) models as a function of ${{\rm{\Sigma }}}_{0}$.⁹ We measure this by integrating the radial component of the momentum flux over the outer boundary of the simulation and over all time, $\int {dt}{\int }_{\partial V}{dA}\rho {\boldsymbol{v}}\cdot \hat{{\boldsymbol{n}}}{\boldsymbol{v}}\cdot \hat{{\boldsymbol{r}}}$, where $\hat{{\boldsymbol{n}}}$ is the normal to area dA. The PH-only runs have momentum yields higher than the RP-only results except for models M1E5R05 and M1E6R25, and close to the PH+RP results. Similar to Figure 7, this implies that thermal pressure induced by photoionization controls the dynamics of H ii regions in low surface density clouds, while radiation pressure is more important in dense, massive clouds. For the PH+RP and PH-only models, most of the outflow momentum is deposited in the ionized rather than neutral gas. We find that the total momentum yields in the PH+RP runs are well described by a power-law relationship:

$$\begin{eqnarray}&&{p}_{* }/{m}_{* }=\displaystyle \frac{{p}_{\mathrm{ej},\mathrm{final}}}{{M}_{* ,\mathrm{final}}}=135\,\mathrm{km}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{0}}{{10}^{2}{M}_{\odot }{\mathrm{pc}}^{-2}}\right)}^{-0.74},\end{eqnarray} \tag{ 18 }$$

shown as the dashed line in Figure 12.

Zoom In Zoom Out Reset image size

We note that by adopting a light-to-mass ratio that is time-independent, our numerical results are likely to overestimate the true feedback strength, especially for clouds with the destruction timescale longer than ${t}_{\mathrm{SN}}=3\,\mathrm{Myr}$ or ${t}_{\mathrm{UV}}=8\,\mathrm{Myr}$. Even with this overestimation, our results show that the momentum injection to the large-scale ISM from both ionizing and non-ionizing radiation feedback is substantially smaller than that from supernovae feedback, p_*/m_* ∼ 10³ km s⁻¹ (e.g., Kim & Ostriker 2015; Kim et al. 2017). Thus, while momentum injection from radiation feedback is quite important on the spatial scale of individual clouds and early life of clusters, it is subdominant in the large-scale ISM compared to other forms of massive star momentum injection, and is therefore not expected to be important in regulating either large-scale ISM pressure and star formation rates (e.g., Ostriker & Shetty 2011; Kim et al. 2013) or galaxy formation (e.g., Naab & Ostriker 2017).

4.2.2. Efficiency of Radial Momentum Injection

Dynamical expansion of a spherical, embedded H ii region in a uniform medium is well described by the thin shell approximation. Force balance ${{\boldsymbol{f}}}_{\mathrm{thm}}+\,{{\boldsymbol{f}}}_{\mathrm{rad}}=0$ holds in the interior (Draine 2011), and most of the swept-up gas is compressed into a thin shell that expands due to the contact forces on it. The effective momentum injection rate from gas pressure and direct radiation pressure are ${\rho }_{{\rm{i}}}{c}_{{\rm{i}}}^{2}{A}_{{\rm{i}}}$ and L/c, respectively (Krumholz & Matzner 2009; Kim et al. 2016). The former expression multiplied by a factor two has also been applied for blister H ii regions, assuming that in addition to thermal pressure, the cloud back-reacts to photoevaporation from a D-critical IF with ${\rho }_{{\rm{i}}}{c}_{{\rm{i}}}^{2}{A}_{{\rm{i}}}={\dot{M}}_{\mathrm{ion}}{c}_{{\rm{i}}}$ (e.g., Matzner 2002; Krumholz et al. 2006). In analytic solutions, the ionized gas density ${\rho }_{i}$ estimate assumes uniform density and ionization equilibrium. For either embedded or blister H ii regions, this leads to an expression for momentum injection that is given by Equation (11) for ${\dot{N}}_{{\rm{e}}}$ multiplied by ${\mu }_{{\rm{H}}}{c}_{{\rm{i}}}$ times an order-unity factor. For an idealized, dust-free, spherical, embedded H ii region, the theoretical point of comparison for the gas pressure force is ${\mu }_{{\rm{H}}}{c}_{i}^{2}{(12\pi {Q}_{{\rm{i}}}{R}_{0}/{\alpha }_{{\rm{B}}})}^{1/2}$.

In reality, H ii regions possess many holes and ionizing sources that are widely distributed in space. The momentum injection is then expected to be less efficient than in the idealized spherical case for the following reasons. First, the non-spherical distribution of radiation sources leads to momentum injection cancellation. Second, for applying back-reaction forces, normal vectors to irregular-shaped cloud surfaces have non-radial components with respect to the cluster center. Third, radiation escapes through holes (see also Dale 2017; Raskutti et al. 2017). Nevertheless, the results in Section 4.1.2 show that the numerically computed mass loss rates from photoevaporation are similar to the predictions from dimensional-analysis scaling arguments, so we may analogously expect the radial momentum injection rate from thermal pressure to be comparable to ${\dot{M}}_{\mathrm{ion}}{c}_{{\rm{i}}}$. This will be lower than the idealized spherical momentum injection, because dust absorption and escape of radiation reduce ${Q}_{{\rm{i}},\mathrm{eff}}$ below ${Q}_{{\rm{i}}}$. Similarly, the radiation pressure applied to the cloud would be reduced to $(1-{f}_{\mathrm{esc}})L/c$ by the escape of radiation, where ${f}_{\mathrm{esc}}=({L}_{{\rm{i}},\mathrm{esc}}+{L}_{{\rm{n}},\mathrm{esc}})/L$ is the overall escape fraction of UV radiation; the total radiation force may be further reduced by lack of spherical symmetry.

To determine the degree of reduction in the radial momentum injection, we calculate the time-averaged, normalized thermal pressure force

$$\begin{eqnarray}{\phi }_{\mathrm{thm}} & \equiv & \displaystyle \frac{{\displaystyle \int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}\displaystyle \int {{\boldsymbol{f}}}_{\mathrm{thm}}\cdot \hat{{\boldsymbol{r}}}\,{dVdt}}{{\displaystyle \int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}{\mu }_{{\rm{H}}}{c}_{{\rm{i}}}^{2}{(12\pi {Q}_{{\rm{i}}}{R}_{0}/{\alpha }_{{\rm{B}}})}^{1/2}\,{dt}},\\ \mathrm{and}\quad {\phi }_{\mathrm{thm}}^{{\prime} } & \equiv & \displaystyle \frac{{\displaystyle \int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}\displaystyle \int {{\boldsymbol{f}}}_{\mathrm{thm}}\cdot \hat{{\boldsymbol{r}}}\,{dVdt}}{{\displaystyle \int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}{\dot{M}}_{\mathrm{ion}}{c}_{{\rm{i}}}\,{dt}},\end{eqnarray} \tag{ 19 }$$

and normalized radiation pressure force

$$\begin{eqnarray}{\phi }_{\mathrm{rad}} & \equiv & \displaystyle \frac{{\displaystyle \int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}\displaystyle \int {{\boldsymbol{f}}}_{\mathrm{rad}}\cdot \hat{{\boldsymbol{r}}}\,{dVdt}}{{\displaystyle \int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}L/c\,{dt}},\\ \mathrm{and}\quad {\phi }_{\mathrm{rad}}^{{\prime} } & \equiv & \displaystyle \frac{{\displaystyle \int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}\displaystyle \int {{\boldsymbol{f}}}_{\mathrm{rad}}\cdot \hat{{\boldsymbol{r}}}\,{dVdt}}{{\displaystyle \int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}(1-{f}_{\mathrm{esc}})L/c\,{dt}},\end{eqnarray} \tag{ 20 }$$

where the range of time integration is taken as $({t}_{{\rm{i}}},{t}_{{\rm{f}}})\,=({t}_{* ,0},{t}_{\mathrm{ej},95 \% })$.¹⁰

In the above, we consider both the radial momentum injection from thermal pressure and radiation pressure relative to the maximum for a spherical shell, and relative to the actual material and radiation momentum available.

Figure 13 plots the normalized forces ${\phi }_{\mathrm{thm}}$ (circles), ${\phi }_{\mathrm{thm}}^{{\prime} }$ (diamonds), ${\phi }_{\mathrm{rad}}$ (squares), and ${\phi }_{\mathrm{rad}}^{{\prime} }$ (triangles) averaged over the time interval $({t}_{* ,0},{t}_{\mathrm{ej},95 \% })$ as functions of ${{\rm{\Sigma }}}_{0}$. There is significant reduction in both the gas and radiation pressure forces compared to the prediction for the "spherical maximum"; the normalized gas pressure force ${\phi }_{\mathrm{thm}}$ is in the range 0.10–0.19 with a mean value 0.13, roughly independent of ${{\rm{\Sigma }}}_{0};$ ${\phi }_{\mathrm{rad}}$ is as small as 0.06 for model M1E5R50 and increases to 0.27 for high surface density clouds. The values of ${\phi }_{\mathrm{thm}}^{{\prime} }$ range from 0.67 to 1.26, suggesting that the time-averaged radial momentum injection from thermal pressure force is within ∼30% of the naive estimate ${M}_{\mathrm{ion}}{c}_{i}$. The values of ${\phi }_{\mathrm{rad}}^{{\prime} }$ are slightly larger than ${\phi }_{\mathrm{rad}}$ and in the range 0.3–0.5, indicating that the momentum injection by radiation pressure force is reduced by a factor of ∼2–3 by flux cancellation alone.¹¹

Zoom In Zoom Out Reset image size

Finally, we evaluate the relative importance of gas and radiation pressure forces in radial momentum injection for the PH+RP models. Figure 14 plots the total injected momentum from gas (circles) and radiation (squares) pressure forces per unit mass of stars formed:

$$\begin{eqnarray}\displaystyle \frac{{p}_{\mathrm{ej},\mathrm{thm}}}{{M}_{* ,\mathrm{final}}} & \equiv & \displaystyle \frac{\iint {{\boldsymbol{f}}}_{\mathrm{thm}}\cdot \hat{{\boldsymbol{r}}}\,{dVdt}}{{M}_{* ,\mathrm{final}}},\\ \mathrm{and}\quad \displaystyle \frac{{p}_{\mathrm{ej},\mathrm{rad}}}{{M}_{* ,\mathrm{final}}} & \equiv & \displaystyle \frac{\iint {{\boldsymbol{f}}}_{\mathrm{rad}}\cdot \hat{{\boldsymbol{r}}}\,{dVdt}}{{M}_{* ,\mathrm{final}}}.\end{eqnarray} \tag{ 21 }$$

The results show trends similar to those exhibited by the total radial momentum yield in Figure 12 for the PH-only versus RP-only models; while the momentum injection due to gas pressure dominates in low surface density clouds, radiation pressure takes over for clouds with ${{\rm{\Sigma }}}_{0}\gtrsim 500\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$. Our numerical results are consistent with the qualitative trend in the previous analytic models (Krumholz & Matzner 2009; Fall et al. 2010; Murray et al. 2010; Kim et al. 2016) in which radiation pressure becomes relatively more important compared to ionized gas pressure for more massive, high surface density clouds.

Zoom In Zoom Out Reset image size

Comparison of the dashed line (Equation (18)) and symbols in Figure 14 shows that the sum of injected momentum from gas and radiation pressure forces ${p}_{\mathrm{ej},\mathrm{thm}}+{p}_{\mathrm{ej},\mathrm{rad}}$ is not necessarily equal to the total momentum of outflowing gas ${p}_{\mathrm{ej},\mathrm{final}}$. This difference arises because of the additional centrifugal and gravitational acceleration that contributes significantly to the latter.

We first consider a situation where the gas left over from star formation is all evaporated by photoionization. In this case, mass conservation requires

$$\begin{eqnarray}&&1={\varepsilon }_{* }+{\varepsilon }_{\mathrm{ion}},\end{eqnarray} \tag{ 22 }$$

where ${\varepsilon }_{\mathrm{ion}}={M}_{\mathrm{ion}}/{M}_{0}={\mu }_{{\rm{H}}}{N}_{{\rm{e}}}/{M}_{0}$.

Using Equation (15), we express ${\varepsilon }_{\mathrm{ion}}$ as

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{ion}}={\phi }_{t}{\phi }_{\mathrm{ion}}\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{ion}}}{{{\rm{\Sigma }}}_{0}}\right){\varepsilon }_{* }^{1/2},\end{eqnarray} \tag{ 23 }$$

where

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{ion}}={\mu }_{{\rm{H}}}{c}_{{\rm{i}}}{\left(\displaystyle \frac{{\rm{\Xi }}}{8G{\alpha }_{{\rm{B}}}}\right)}^{1/2}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&=\,140{\left(\displaystyle \frac{{\rm{\Xi }}}{5.05\times {10}^{46}{{\rm{s}}}^{-1}{M}_{\odot }^{-1}}\right)}^{1/2}{\left(\displaystyle \frac{{T}_{\mathrm{ion}}}{8000{\rm{K}}}\right)}^{0.85}\,{M}_{\odot }\,{\mathrm{pc}}^{-2},\end{eqnarray} \tag{ 25 }$$

with ${\rm{\Xi }}={Q}_{{\rm{i}},\max }/{M}_{* ,\mathrm{final}}$ being the conversion factor between stellar mass and ionizing photon output. Since ${\phi }_{\mathrm{ion}}$ is roughly constant in our numerical results, Equation (23) implies that the evaporation efficiency depends linearly on the duration (relative to ${t}_{\mathrm{ff}}$) of photoevaporation, ${\phi }_{t}$ (Equation (14)), and the square root of the SFE, and inversely on the cloud surface density, ${{\rm{\Sigma }}}_{0}$. The characteristic surface density ${{\rm{\Sigma }}}_{\mathrm{ion}}$ can be regarded as a constant for most of the models since ${M}_{* ,\mathrm{final}}\gt {10}^{3}\,{M}_{\odot }$, while it varies sensitively with ${M}_{* ,\mathrm{final}}$ in models with ${M}_{* ,\mathrm{final}}\lt {10}^{3}\,{M}_{\odot }$, as described in Section 2.1. The dependence of ${\varepsilon }_{\mathrm{ion}}$ on ${\varepsilon }_{* }^{1/2}$ (rather than ${\varepsilon }_{* }$) is caused by shielding within the IBSV, which makes photoevaporation inefficient. Most photons are used up in offsetting recombination within the H ii region, rather than in eroding neutral gas to create new ions at IFs.

We solve Equations (22) and (23) for ${\varepsilon }_{* }$ to obtain

$$\begin{eqnarray}&&{\varepsilon }_{* }={\left(\displaystyle \frac{2\xi }{1+\sqrt{1+4{\xi }^{2}}}\right)}^{2},\end{eqnarray} \tag{ 26 }$$

where $\xi \equiv {{\rm{\Sigma }}}_{0}/({\phi }_{t}{\phi }_{\mathrm{ion}}{{\rm{\Sigma }}}_{\mathrm{ion}})$. In the limit $\xi \ll 1$, ${\varepsilon }_{* }\approx {\xi }^{2}\,={[{{\rm{\Sigma }}}_{0}/({\phi }_{t}{\phi }_{\mathrm{ion}}{{\rm{\Sigma }}}_{\mathrm{ion}})]}^{2}$, while in the limit $\xi \gg 1$, ${\varepsilon }_{* }\approx 1-1/\xi \,=1-{\phi }_{t}{\phi }_{\mathrm{ion}}{{\rm{\Sigma }}}_{\mathrm{ion}}/{{\rm{\Sigma }}}_{0}$. The limiting behavior shows that photoevaporation becomes very efficient for destroying clouds when ${{\rm{\Sigma }}}_{0}\lesssim {{\rm{\Sigma }}}_{\mathrm{ion}}$.

Inserting the fit for ${\phi }_{t}{\phi }_{\mathrm{ion}}$ of Equation (16) in (26) allows us to calculate the net SFE as a function of ${{\rm{\Sigma }}}_{0}$. The resulting ${\varepsilon }_{* }$ and ${\varepsilon }_{\mathrm{ion}}$ are compared to the numerical results as dashed lines in Figures 6(a) and (c). For low surface density and massive clouds with ${\varepsilon }_{\mathrm{ion}}\gtrsim 0.7$, the semi-analytic ${\varepsilon }_{* }$ agrees with the numerical results extremely well. The agreement is less good, however, for high surface density (${{\rm{\Sigma }}}_{0}\gtrsim 300\,{M}_{\odot }\ {\mathrm{pc}}^{-2}$) or low mass (${M}_{0}={10}^{4}\,{M}_{\odot }$) clouds for which neutral gas ejection cannot be ignored; we address this next.

We consider the general case in which radiative feedback from star formation exerts combined thermal and radiation pressure forces on the surrounding gas. The gas is ejected from the cloud, quenching further star formation. In our simulations, the momentum ejected from the computational domain also includes a small contribution from the initial turbulence, which ejects a mass fraction ${\varepsilon }_{\mathrm{ej},\mathrm{turb}}={M}_{\mathrm{ej},\mathrm{turb}}/{M}_{0}\sim 0.1$ when the cloud is initially marginally bound.

Let ${v}_{\mathrm{ej}}$ denote the characteristic ejection velocity of the outgoing gas, and let ${p}_{* }/{m}_{* }$ denote the momentum per stellar mass injected by feedback. The total momentum of the ejected gas can be written as

$$\begin{eqnarray}&&{p}_{\mathrm{ej},\mathrm{tot}}=(1-{\varepsilon }_{* }){M}_{0}{v}_{\mathrm{ej}}=({p}_{* }/{m}_{* }){\varepsilon }_{* }{M}_{0}+{\varepsilon }_{\mathrm{ej},\mathrm{turb}}{M}_{0}{v}_{\mathrm{ej}}.\end{eqnarray} \tag{ 27 }$$

We divide the both sides of Equation (27) by ${M}_{0}{v}_{\mathrm{ej}}$ to obtain

$$\begin{eqnarray}&&{\varepsilon }_{* }=\displaystyle \frac{1-{\varepsilon }_{\mathrm{ej},\mathrm{turb}}}{1+({p}_{* }/{m}_{* })/{v}_{\mathrm{ej}}}.\end{eqnarray} \tag{ 28 }$$

From the results of model M1E5R20_nofb (see Table 2), we take ${\varepsilon }_{\mathrm{ej},\mathrm{turb}}=0.13$.¹² For ${v}_{\mathrm{ej}}$, we take the mass-weighted outflow velocity $({\varepsilon }_{\mathrm{ej},\mathrm{neu}}{v}_{\mathrm{ej},\mathrm{neu}}+{\varepsilon }_{\mathrm{ej},\mathrm{ion}}{v}_{\mathrm{ej},\mathrm{ion}})/{\varepsilon }_{\mathrm{ej}}$, which is found to be roughly constant for fixed ${M}_{0}$ with mean values 15, 23, 29 km s⁻¹ for ${M}_{0}={10}^{4}$, 10⁵, ${10}^{6}\,{M}_{\odot }$, respectively. Equation (28) gives ${\varepsilon }_{* }$ as functions of ${M}_{0}$ and ${{\rm{\Sigma }}}_{0}$ once ${p}_{* }/{m}_{* }$, ${v}_{\mathrm{ej}}$, and ${\varepsilon }_{\mathrm{ej},\mathrm{turb}}$ are specified. A fit to the momentum yield ${p}_{* }/{m}_{* }$ from our numerical results as a function of ${{\rm{\Sigma }}}_{0}$ is given in Equation (18). Altogether, Equation (28) then gives ${\varepsilon }_{* }$ as a function of ${M}_{0}$ and ${{\rm{\Sigma }}}_{0}$. For ${p}_{* }/{m}_{* }\gg {v}_{\mathrm{ej}}$ and ${\varepsilon }_{\mathrm{ej},\mathrm{turb}}\ll 1$, the predicted scaling dependence is ${\varepsilon }_{* }\sim {v}_{\mathrm{ej}}/({p}_{* }/{m}_{* })\propto {v}_{\mathrm{ej}}{{\rm{\Sigma }}}_{0}^{3/4}$.

Figure 6 overplots as solid lines the resulting (a) net SFE ${\varepsilon }_{* }$, (b) ejection efficiency ${\varepsilon }_{\mathrm{ej}}=1-{\varepsilon }_{* }$, (c) photoevaporation efficiency ${\varepsilon }_{\mathrm{ion}}$ calculated from Equation (23), and (d) neutral ejection efficiency ${\varepsilon }_{\mathrm{ej},\mathrm{neu}}={\varepsilon }_{\mathrm{ej}}-{\varepsilon }_{\mathrm{ion}}$. The line thickness indicates the initial cloud mass. Overall, our semi-analytic model for the dynamical mass ejection reproduces the net SFE of the simulations fairly well, with the average difference of 0.03. The model also explains the increasing tendency of ${\varepsilon }_{* }$ with increasing ${M}_{0}$. This weak dependence of the semi-analytic ${\varepsilon }_{* }$ on ${M}_{0}$ comes directly from ${v}_{\mathrm{ej}}$, indicating that stronger feedback is required to unbind gas in a more massive cloud (e.g., Kim et al. 2016; Rahner et al. 2017). A drop in the photoevaporation efficiency with decreasing ${{\rm{\Sigma }}}_{0}$ for ${M}_{0}={10}^{4}\,{M}_{\odot }$ and ${{\rm{\Sigma }}}_{0}\lesssim {10}^{2}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ is caused by the steep dependence of Ξ on ${M}_{* ,\mathrm{final}}$ below ${10}^{3}\,{M}_{\odot }$, as mentioned above.

• 1.  
  Evolutionary stages. All clouds in our simulations go through the following evolutionary stages (Figure 3), with some overlap and varying duration. (1) The initial turbulence creates filaments, and within these, denser clumps condense. Some clumps undergo gravitational collapse and form star particles representing subclusters. (2) Individual H ii regions form around each subcluster, growing toward directions of low optical depth until they merge and break out of the cloud. Simultaneously, both low- and high-density gas is accelerated away from the radiation sources due to thermal and radiation pressures. (3) At late stages of evolution, the brightest EM features are small globules, pillars, and ridges marking individual IFs at the periphery of the H ii region, quite similar to observed clouds (Figure 5). (4) At the end of the simulation, all of the gas has either collapsed into stars or been dispersed, flowing out of the computational domain.
• 2.  
  Timescales. Cloud destruction takes ∼2–$10\,\mathrm{Myr}$ after the onset of massive star formation feedback (Figure 8(a)). In units of the freefall time, the destruction timescale is in the range $0.6\lt {t}_{\mathrm{dest}}/{t}_{\mathrm{ff},0}\lt 4.1$ and systematically increases with ${{\rm{\Sigma }}}_{0}$ (Figure 8(b)). The timescale for star formation is comparable to or somewhat smaller than the destruction timescale. The effective gas depletion timescale ranges from $250\,\mathrm{Myr}$ to $2.7\,\mathrm{Myr}$, sharply decreasing with increasing ${{\rm{\Sigma }}}_{0}$.
• 3.  
  Star formation and mass loss efficiencies. We examine the dependence on the cloud mass and surface density of the net SFE ${\varepsilon }_{* }$, photoevaporation efficiency ${\varepsilon }_{\mathrm{ion}}$, and mass ejection efficiency ${\varepsilon }_{\mathrm{ej}}$ (Figure 6). The SFE ranges over ${\varepsilon }_{* }=0.04\mbox{--}0.51$, increasing strongly with the initial surface density ${{\rm{\Sigma }}}_{0}$ while increasing weakly with the initial cloud mass ${M}_{0}$. Photoevaporation accounts for more than 70% of mass loss for clouds with ${M}_{0}\gtrsim {10}^{5}\,{M}_{\odot }$ and ${{\rm{\Sigma }}}_{0}\lesssim {10}^{2}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$. The ejection of neutral gas mass by thermal and radiation pressures also contributes in quenching further star formation in low mass and high surface density clouds. The comparison of the net SFE among the models in which we turn on and off photoionization and radiation pressure suggests that photoionization is of greater importance in destroying GMCs in normal disk galaxies, whereas radiation pressure is more effective in regulating star formation in dense, massive clouds (Figure 7; see also Figure 14). This controlled experiment also demonstrates that the combined effects of photoionization and radiation pressure do not work in a simply additive manner in suppressing star formation.
• 4.  
  Photoevaporation. The photoevaporation rate ${\dot{N}}_{{\rm{e}}}$ (the number of free electrons produced per unit time) at IFs within the cloud is much smaller than the total number of ionizing photons absorbed by hydrogen per unit time ${Q}_{{\rm{i}},\mathrm{eff}}$ throughout the cloud (Figure 10). Most of the ionizations instead offset recombinations in diffuse gas throughout the cloud. The time-averaged hydrogen absorption fraction $\langle {f}_{\mathrm{ion}}\rangle =\langle {Q}_{{\rm{i}},\mathrm{eff}}/{Q}_{{\rm{i}}}\rangle$ ranges from 0.19 to 0.50, while the time-averaged shielding factor $\langle q\rangle =\langle {Q}_{{\rm{i}},\mathrm{eff}}/{\dot{N}}_{{\rm{e}}}\rangle$ ranges from 61 to 2950 (Table 2). Although dense, compact clouds tend to have a high hydrogen absorption fraction, they have a high shielding factor and hence inefficient photoevaporative mass loss. Assuming that the area of the IFs and thickness of the shielding layer scale with the dimensions of the initial cloud, we derive and calibrate expressions for the photoevaporation rate (Equations(11), (13); Figure 11(a)).
• 5.  
  Outflow acceleration and properties. We perform a detailed analysis of the momentum injection processes that are responsible for cloud disruption and mass loss. The total radial momentum yield (momentum per stellar mass formed) of outflowing gas ranges over ∼20–400 km s⁻¹ (Figure 12), scaling as ${p}_{* }/{m}_{* }\propto {{\rm{\Sigma }}}_{0}^{-0.74}$ (Equation (18)). The time-averaged total radial gas pressure force is smaller than the dustless, spherical case by a factor of ∼5–10 due to momentum cancellation and escape of radiation, but is within 30% of ${\dot{M}}_{\mathrm{ion}}{c}_{{\rm{i}}}$, where ${\dot{M}}_{\mathrm{ion}}$ and ${c}_{{\rm{i}}}$ refer to the mass evaporation rate and the sound speed of the ionized gas, respectively (Figure 13). This is consistent with expectations for both internal pressure forces within the ionized medium and the combined thermal pressure and recoil forces on neutral gas at IFs, both of which scale as ${n}_{{\rm{i}}}{c}_{{\rm{i}}}^{2}{A}_{{\rm{i}}}\sim {\dot{M}}_{\mathrm{ion}}{c}_{{\rm{i}}}$. Similarly, the time-averaged radial force from radiation pressure is much reduced below the naive spherical expectation to ∼0.08–0.27 L/c because of flux cancellation and photon escape. Although the overall momentum injection by radiation feedback is less efficient in realistic turbulent clouds than suggested by analytic predictions based on spherical H ii region expansion in smooth clouds (e.g., Krumholz & Matzner 2009; Fall et al. 2010; Murray et al. 2010; Kim et al. 2016), the ratio of radiation pressure forces to gas pressure forces increases for massive, high surface density clouds (Figure 14), consistent with expectations from these previous studies. The mean outflow velocity of ionized gas is mildly supersonic with ${v}_{\mathrm{ej},\mathrm{ion}}\sim 18$–$36\,\mathrm{km}\,{{\rm{s}}}^{-1}$, while that of neutral gas is v_ej,neu ∼ 6–15 km s⁻¹, at about 0.8–1.8 times the escape velocity of the initial cloud (Figure 9, Table 2).
• 6.  
  Semi-analytic models. Based on our analyses of mass loss processes, we develop simple semi-analytic models for the net SFE, ${\varepsilon }_{* }$, and photoevaporation efficiency, ${\varepsilon }_{\mathrm{ion}}$, as limited by radiation feedback in cluster-forming clouds. The predicted ${\varepsilon }_{\mathrm{ion}}$ is proportional to ${\varepsilon }_{* }^{1/2}{{\rm{\Sigma }}}_{0}^{-1}$ (Equation (23)), with an order-unity dimensionless coefficient ${\phi }_{t}{\phi }_{\mathrm{ion}}$ that we calibrate from simulations (Equation (16), Figure 11). When photoevaporation is the primary agent of cloud destruction, the net SFE depends solely on the gas surface density (Equation (26)). The resulting predictions for ${\varepsilon }_{* }$ and ${\varepsilon }_{\mathrm{ion}}$ agree well with the numerical results for massive ($\geqslant {10}^{5}\,{M}_{\odot }$) clouds (Figures 6(a) and (c)). In low mass clouds (${10}^{4}\,{M}_{\odot }$), the back-reaction to ionized gas pressure is effective at IFs, and these clouds lose 30%–50% of their initial gas mass as neutrals. Allowing for both ionized and neutral gas in outflows and assuming that the ejection velocity is constant for fixed cloud mass, with total momentum injection calibrated from our simulations (Equation (18)), the predicted net SFE (Equation (28)) has a weak dependence on the cloud mass, consistent with our numerical results (Figure 6).

It is interesting to compare our results with those of previous theoretical studies on star-forming GMCs with UV radiation feedback. Our net SFE and relative role of photoionization to radiation pressure are in qualitative agreement with the predictions of Kim et al. (2016), which adopted the idealizations of instantaneous star formation and spherical shell expansion (see also Fall et al. 2010; Murray et al. 2010; Rahner et al. 2017). However, the minimum SFE required for cloud disruption by ionized gas pressure found by Kim et al. (2016; ∼0.002–0.1 for $20\,{M}_{\odot }\,{\mathrm{pc}}^{-2}\leqslant {{\rm{\Sigma }}}_{0}\leqslant {10}^{2}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$, see their Figure 11) is smaller than what we find here. Turbulence-induced structure can lead to higher SFE for several reasons. First, a fraction of the photons can easily escape through low-density channels, without either ionizing gas or being absorbed to impart momentum. Second, low-density but high pressure ionized gas can also vent through these channels, reducing the transfer of momentum from photoevaporated to neutral gas. Third, turbulence increases the mass-weighted mean density and therefore the recombination rate of ionized gas, so that a higher luminosity is required to photoevaporate gas.

In this paper, we also find that radiation pressure becomes significant at ${{\rm{\Sigma }}}_{0}\sim 200$–$800\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ (lower for higher-mass clouds). This transition value of surface density for low mass clouds is somewhat higher than ${{\rm{\Sigma }}}_{0}\sim 100\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ in simple spherical models (e.g., Krumholz & Matzner 2009; Fall et al. 2010; Murray et al. 2010; Kim et al. 2016). These spherical models predicted that the expansion of an H ii region is dominated by radiation pressure force at least in the early phase of expansion (e.g., Krumholz & Matzner 2009; Kim et al. 2016). Even for the highest-${{\rm{\Sigma }}}_{0}$ models, however, we find that the volume-integrated radiation pressure force ($\int {{\boldsymbol{f}}}_{\mathrm{rad}}\cdot \hat{{\boldsymbol{r}}}\,{dV}$) in the radial direction tends to be smaller than that of gas pressure force ($\int {{\boldsymbol{f}}}_{\mathrm{thm}}\cdot \hat{{\boldsymbol{r}}}\,{dV}$) in the early phase of star formation, and that this tendency is reversed only at late time. This is due to the strong cancellation of radiation on a global scale, and should thus not be interpreted as evidence for radiation pressure being subdominant. Rather, radiation pressure plays a greater role in controlling the dynamics of sub H ii regions surrounding individual sources.

The outflow momentum yield of ${p}_{* }/{m}_{* }\sim 20$–400 km s⁻¹ by radiation feedback found from our simulations can be compared with those produced by other types of feedback. For example, the momentum yield of protostellar outflows is estimated as ∼40 km s⁻¹ (e.g., Matzner & McKee 2000), while recent numerical work on SNe-driven flows found a momentum yield of ∼1000–3000 km s⁻¹, weakly dependent on the background density and stellar clustering (e.g., Kim & Ostriker 2015; Kim et al. 2017). These feedback mechanisms are expected to be important at different scales (e.g., Fall et al. 2010; Krumholz et al. 2014; Matzner & Jumper 2015). The momentum injection by protostellar outflows plays a key role in driving turbulence and regulating star formation in small-scale clumps before massive stars form (e.g., Nakamura & Li 2014), and SNe are capable of driving ISM turbulence to regulate galaxy-wide star formation and possibly winds from low mass galaxies (e.g., Kim et al. 2017). Momentum injection by UV radiation is important to controlling dynamics at intermediate, GMC scales.

In this paper we have focused exclusively on the effects of radiation feedback, and as we neglect aging of stellar populations over the whole cloud lifetime, our mass loss efficiencies place an upper limit to the actual damage UV radiation can cause. To estimate lower limits to the destructive effects of radiation feedback, we can consider what UV feedback would be able to accomplish over two shorter intervals: ${t}_{\mathrm{SN}}=3\,\mathrm{Myr}$, the minimum lifetime of most massive stars (i.e., the time when first supernova is expected to occur after the epoch of the first star formation); and ${t}_{\mathrm{UV}}=8\,\mathrm{Myr}$, the timescale on which UV luminosity decays. Figure 15 shows the cumulative (a) ejected ion efficiency $\mathrm{EIE}={M}_{\mathrm{ej},\mathrm{ion}}/{M}_{0}$ and (b) ejected neutral efficiency $\mathrm{ENE}={M}_{\mathrm{ej},\mathrm{neu}}/{M}_{0}$ as functions of the initial cloud radius ${R}_{0}$ at t_*,0+t_SN (magenta), t_*,0+t_UV (cyan), and ${t}_{\mathrm{final}}$ (gray). In the case of compact, high surface density clouds with short freefall times, the mass loss efficiencies at t_*,0+t_SN are close to the final values at ${t}_{\mathrm{final}}$, suggesting that UV radiation is expected to clear out most of gas prior to the first supernova. For large, diffuse clouds with long freefall times, the mass ejection efficiencies at t_*,0+t_SN and t_*,0+t_UV are significantly smaller than the final values. For these diffuse clouds with long freefall timescales, a complete assessment of SFE and cloud destruction will have to include the effects of supernovae as well as radiation feedback. While some of the supernova energy escapes easily through low-density channels (e.g., Rogers & Pittard 2013; Walch & Naab 2015), and supernovae may undermine radiation feedback by compressing overdense structures, supernova feedback likely aids in destroying star-forming clouds overall (e.g., Geen et al. 2016). In principle, shocked stellar winds may also be important over the same period as radiation feedback is active, although recent studies have raised doubts about its effectiveness, as this hot gas can easily leak through holes or mix with cool gas (Harper-Clark & Murray 2009; Rosen et al. 2014).

Zoom In Zoom Out Reset image size

In this work, we have considered only unmagnetized, marginally bound clouds with ${\alpha }_{\mathrm{vir},0}=2$, while observed GMCs may have a range of virial values (e.g., Heyer et al. 2009; Roman-Duval et al. 2010; Miville-Deschênes et al. 2017). Preliminary simulations we have done reveal that the net SFE decreases from 0.31 to 0.02 as ${\alpha }_{\mathrm{vir},0}$ is varied from 0.5 to 5 for the fiducial mass and size. Clouds with large initial ${\alpha }_{\mathrm{vir},0}$ tend to form stars less efficiently, since turbulence unbinds a larger fraction of gas and reduces the gas mass in collapsing structures, as has been reported by recent simulations (e.g., Dale et al. 2013; Bertram et al. 2015; Howard et al. 2016; Raskutti et al. 2016; Dale 2017). For clouds that are initially strongly bound, an initial adjustment period leads to smaller clouds with order-unity α_vir, in which the subsequent evolution is similar to clouds with ${\alpha }_{\mathrm{vir},0}\sim 1$ (Raskutti et al. 2016). Magnetization is likely to increase the effectiveness of radiation feedback, as it reduces the density inhomogeneity that limits radiation pressure effects, and may also help distribute energy into neutral gas (Gendelev & Krumholz 2012).

Finally, we comment on the resolution of our numerical models. Our standard choice of grid resolution (N_cell = 256³) and domain size (${L}_{\mathrm{box}}=4{R}_{0}$) is a compromise between computational time and accuracy, and this choice of moderate resolution enabled us to extensively explore parameter space. While the present simulations can capture the dynamics of cluster-forming gas and outflowing gas on cloud scales with reasonable accuracy, they do not properly resolve compressible flows in high-density regions where sink particles are created. Our models may overestimate star formation rates, as all the gas that has been accreted onto sink particles is assumed to be converted to stars (cf. Howard et al. 2016, 2017), although a converging trend of ${\varepsilon }_{* }$ with resolution (see Appendix A) suggests that the final SFE may be primarily controlled by photoevaporation and injection of momentum in moderate density gas at large scales. Adaptive mesh refinement simulations can be used to address this and other resolution-related questions.