# Feedback in Clouds II: UV photoionization and the first supernova in a massive cloud

Sam Geen, Patrick Hennebelle, Pascal Tremblin, Joakim A Rosdahl

► To cite this version:

Sam Geen, Patrick Hennebelle, Pascal Tremblin, Joakim A Rosdahl. Feedback in Clouds II: UV photoionization and the first supernova in a massive cloud. Monthly Notices of the Royal Astronomical Society, 2016, 463 (3), pp.3129 - 3142. 10.1093/mnras/stw2235. hal-01834131

HAL Id: hal-01834131

https://hal.science/hal-01834131

Submitted on 11 Jan 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Feedback in Clouds II: UV photoionization and the first supernova in a massive cloud

Sam Geen,<fnref n="1" />* Patrick Hennebelle,<fnref n="1" /> Pascal Tremblin<fnref n="2" /> and Joakim Rosdahl<fnref n="3" />

Accepted 2016 September 2. Received 2016 August 31; in original form 2016 June 20

## ABSTRACT

Molecular cloud structure is regulated by stellar feedback in various forms. Two of the most important feedback processes are UV photoionization and supernovae from massive stars. However, the precise response of the cloud to these processes, and the interaction between them, remains an open question. In particular, we wish to know under which conditions the cloud can be dispersed by feedback, which, in turn, can give us hints as to how feedback regulates the star formation inside the cloud. We perform a suite of radiative magnetohydrodynamical simulations of a $10^5$ solar mass cloud with embedded sources of ionizing radiation and supernovae, including multiple supernovae and a hypernova model. A UV source corresponding to 10 per cent of the mass of the cloud is required to disperse the cloud, suggesting that the star formation efficiency should be of the order of 10 per cent. A single supernova is unable to significantly affect the evolution of the cloud. However, energetic hypernovae and multiple supernovae are able to add significant quantities of momentum to the cloud, approximately $10^{43}$ g cm s$^{-1}$ of momentum per $10^{51}$ erg of supernova energy. We argue that supernovae alone are unable to regulate star formation in molecular clouds. We stress the importance of ram pressure from turbulence in regulating feedback in molecular clouds.

**Key words:** methods: analytical – methods: numerical – stars: massive – ISM: clouds – H II regions – ISM: supernova remnants.

## 1 INTRODUCTION

Massive stars release large quantities of energy into their environment. They produce protostellar jets, winds, radiation across a wide spectrum and supernovae. The first phase of stellar feedback occurs in dense molecular cloud environments in which the stars are born. In order for the energy from stars to propagate into the wider interstellar medium (ISM), it must first escape this cloud environment, either by destroying the cloud or creating sufficient channels through which the propagating shocks can escape.

In the previous paper, Geen et al. (2015b), we determined a limit at which ionizing radiation can escape molecular clouds using both numerical simulations and an analytic model. This model is based on the arguments made in Tremblin et al. (2014) and Didelon et al. (2015), which compare models of H II regions expanding into turbulent environments to observed H II regions. These models were constructed using previous analytic theory by Kahn (1954), Spitzer (1978), Whitworth (1979), Franco, Tenorio-Tagle & Bodenheimer

(1990), Williams & McKee (1997), Hosokawa & Inutsuka (2006) and Raga, Canto & Rodriguez (2012).

In Geen et al. (2015b), we proposed a limit at which ionizing photons are able to destroy their host cloud. This extends the argument of Dale, Ercolano & Bonnell (2012), who consider the case where the ionization front cannot expand beyond the initial Strömgren radius. We argue that if a calculated ‘stall’ radius (Draine & Woods 1991) is smaller than the radius of the cloud, the ionization front cannot escape the cloud. This, in turn, sets the ability for ionizing radiation to regulate the environments in which stars form, and determines whether ionizing radiation can suppress the star formation rate of molecular clouds.

Massive stars typically end their lives as supernovae (for estimates of which stars become supernovae, see e.g. Heger et al. 2003). The evolution of the supernova remnant depends on the environment into which it expands. Understanding the momentum deposition from supernovae in star-forming environments is crucial to understanding processes in galaxies as a whole. Sub-grid models by e.g. Hopkins et al. (2014) and Kimm et al. (2015) attempt to correct for a lack of numerical resolution by depositing a pre-calculated quantity of momentum around the supernova if the resolution is insufficient to resolve the blastwave properly (see

<footnote n="1">¹Laboratoire AIM, Paris-Saclay, CEA/IRFU/SAp – CNRS – Université Paris Diderot, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="2">²Maison de la Simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="3">³Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA, Leiden, the Netherlands</footnote>

<footnote>* E-mail: samuel.geen@cea.fr</footnote>
also Kim & Ostriker 2015, for a study of numerical limits on re-solving supernova blastwaves). Analytic and 1D simulation work by Chevalier (1974), Cioffi, McKee & Bertschinger (1988), Draine & Woods (1991), Thornton et al. (1998) and Haid et al. (2016) provides insights into this process, with simulations of supernova blastwaves by Iffrig & Hennebelle (2015), Kim & Ostriker (2015), Martizzi, Faucher-Giguere & Quataert (2015) and Körtgen et al. (2016) extending this to more complex environments using 3D numerical simulations. Supernovae shock against the surrounding medium, expanding adiabatically (Sedov 1946). Eventually, they reach a point where they begin to lose a significant fraction of their energy to radiative cooling (see the estimates by Cox 1972). After a longer period of time, the supernova remnant begins to merge with the surrounding medium (Cioffi et al. 1988).

Pre-supernova feedback, as either stellar winds or ionizing radiation, has been found in simulations to enhance the final energy and momentum of the supernova remnant by injecting additional momentum and reducing the density of the environment into which the supernova occurs (Dwarkadas 2007; Fierlinger et al. 2015; Geen et al. 2015a). Rogers & Pittard (2013) and Walch & Naab (2015) have had some success in driving outflows in simulations of molecular clouds with both supernova and pre-supernova stellar feedback. However, Draine & Woods (1991) suggest that if the medium is sufficiently turbulent, the H II region will re-collapse before the supernova occurs, depending on the mass of the progenitor and the density of the surrounding medium. Krause et al. (2016) find that stellar feedback in very massive extragalactic clouds is ineffective at reducing the star formation efficiency.

In this paper, we explore the competition between pre-supernova ionizing feedback and turbulence in molecular clouds, and the resulting evolution of the supernova remnant as it expands into the environment resulting from this competition. We simulate ionizing radiation and supernovae in a turbulent cloud using RAMSES-RT (Teyssier 2002; Fromang, Hennebelle & Teyssier 2006; Rosdahl et al. 2013). The cloud is $10^5 M_\odot$, 10 times more massive than the one studied in the previous paper. We choose this cloud mass because the slope of the cloud mass function ($dN/dM_c = M_c^{-1.7}$) means that more mass is expected to be found in clouds above $10^5 M_\odot$ (see the review by Hennebelle & Falgarone 2012). Therefore, most of the stars in our Galaxy are expected to form in these clouds. It is thus important to study these objects if we wish to understand how feedback from massive stars interacts with both the host cloud and the wider Galactic ISM.

We begin in Section 2 by presenting the simulations performed. We then extend the analysis of H II regions in the previous paper to a more massive cloud in Section 3, and produce simple models that describe the time-dependence of the evolution of the ionization front. In Section 4, we analyse the results of simulations that introduce a supernova into the cloud and H II region after the source of UV photons is extinguished, and how this compares to previous analytic and numerical theory. We extend the single supernova scenario to more energetic events in Section 5.

## 2 METHODS

We use the radiative magnetohydrodynamics code RAMSES-RT (Teyssier 2002; Fromang et al. 2006; Rosdahl et al. 2013). The system is described by an isolated turbulent, magnetized, self-gravitating initially spherical cloud placed at the centre of the simulation volume. After 2.53 Myr (one free-fall time for the cloud as a whole), we turn on a constant source of ionizing UV photons in the centre of the simulation volume. After 3 Myr, we turn off the

<table><caption>**Table 1.** Table of the simulations included in this paper. 'N' refers to the number of photons deposited per second by the source, $S_*$, in all frequency groups in photons per second as a power of 10 (with 'N00' referring to a zero photon emission rate). 'NSN' means that no supernova is included. 'SN' means that a single $10^{51}$ erg supernova is included. 'HN' means that a hypernova (modelled as a $10^{52}$ erg blast) is included. 'MSN' means that 10 supernovae are included, with $10^{51}$ erg every 0.1 Myr. See Section 2 for more details about the simulation setup.</caption><thead><tr><th>Name</th><th>log<sub>10</sub>(S<sub>*</sub>/s<sup>-1</sup>)</th><th>Supernova?</th></tr></thead><tbody><tr><td>N00-NSN</td><td>(no photons)</td><td>✗</td></tr><tr><td>N49-NSN</td><td>49</td><td>✗</td></tr><tr><td>N50-NSN</td><td>50</td><td>✗</td></tr><tr><td>N51-NSN</td><td>51</td><td>✗</td></tr><tr><td>N00-SN</td><td>(no photons)</td><td>✓</td></tr><tr><td>N49-SN</td><td>49</td><td>✓</td></tr><tr><td>N50-SN</td><td>50</td><td>✓</td></tr><tr><td>N51-SN</td><td>51</td><td>✓</td></tr><tr><td>N50-HN</td><td>50</td><td>Hypernova</td></tr><tr><td>N50-MSN</td><td>50</td><td>10 × SN</td></tr></tbody></table>

source of photons and inject a thermal blast representing a supernova. More details on each component are given in the following sub-sections. Table 1 lists all the simulations used in this paper.

### 2.1 Initial conditions

In this simulation, we consider one set of initial conditions only. For a theoretical description of the effect of varying cloud properties on the shape of H II regions, see Geen et al. (2015b). These initial conditions are similar to the setup given in Iffrig & Hennebelle (2015), Geen et al. (2015b) and Lee & Hennebelle (2016a).

The simulation volume is a cubic box of length 86.3 pc. This volume is divided into 256 cells on a side. We allow further two levels of refinement (i.e. sub-division of a single cell into 8 cells), giving 1024 cells on a side effective resolution. The simulation thus has a minimum spatial resolution of 0.33 pc everywhere and 0.084 pc in the most refined regions. At all times, if a cell is found to be Jeans-unstable, it is allowed to refine up to the maximum level. Additionally, shortly before the supernova is launched, we fully refine the central 1.5 pc of the simulation volume in order to capture the shock evolution properly. For N51-SN and N51-NSN, we run identical simulations with twice the box length but identical spatial physical resolution and refinement criteria in order to follow shocks that would otherwise escape the box.

The mass of the cloud studied in this paper is set to $10^5 M_\odot$. We impose a spherically symmetric density profile on to the simulation volume (see the left-hand panel of Fig. 1). This is given by $n(r, t = 0) = n_0/(1 + (r/r_c)^2)$ for hydrogen number density $n$ at radius $r$ and time $t = 0$ with peak density $n_0 = 2850$ atoms cm$^{-3}$ and characteristic radius $r_c = 3.6$ pc. We impose a cut-off at $3 r_c$ (where $n(3r_c, t = 0) = 0.1 n_0$) Outside this radius, a uniform density field of density $20$ atoms cm$^{-3}$ is imposed out to $21$ pc. Beyond this, the density field is a uniform $1$ atoms cm$^{-3}$. The cloud has a global free-fall time $t_{ff} = 2.53$ Myr, defined as $3\sqrt{\frac{3\pi}{32} \frac{G n_0 m_H}{X}}$, where $m_H$ is the mass of a hydrogen atom and $X = 0.76$ is the hydrogen mass fraction. The temperature inside the inner part of the cloud is set to $10$ K, while the temperature in the medium outside the cloud is set to $4000$ K. The magnetic field is initially $25$ µG in the centre of the cloud and $4.2$ µG outside, chosen such that the ratio between
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Density profile of the cloud. The plot on the left gives the profile at the start of the simulation, while the plot on the right is taken after one free-fall time, just before the H II source is turned on. We cast a set of rays from the source position to the edge of the simulation volume in evenly spaced directions as described in Appendix E of Geen et al. (2015b). Along each ray, we sample the density profile. The minimum and maximum values at each radius are given as a thick grey line. The median value at each radius is given as a solid black line, while the mean is a black dashed line. The interquartile range is bounded by black dotted lines. Each second percentile of probability from the median value is shaded from red to white, with the closest percentiles to the median shaded as red and the furthest as white.</figcaption></figure>

the free-fall time and the Alfvén crossing time is 0.2. Note that the magnetic field strength increases as the cloud evolves.

A turbulent velocity field is imposed over the grid, such that the kinetic energy in turbulence in the cloud is approximately equal to the gravitational energy of the cloud. The turbulence has a Kolmogorov power spectrum (i.e. $P(k) \propto k^{-5/3}$) with random phases.

## 2.2 Radiative transfer

RAMSES-RT (Rosdahl et al. 2013) uses a first-order moment method for the advection of photons, closing the set of equations with the local M1 expression for the radiation pressure tensor. It tracks the ionization states of hydrogen and helium in the gas, and couples the interactions between the photons and the gas on-the-fly. We split the radiation into three groups, bracketed by the ionization energies of H I, He I and He II (13.6, 24.6 and 54.2 eV for the lower bounds of each). In all of the simulations in this paper, we assume a Solar metallicity everywhere at all times. We do not include photons below the ionization energy of hydrogen, nor do we include radiation pressure (as in Rosdahl & Teyssier 2015). A reduced speed of light of $10^{-4}$ c (= 30 km s$^{-1}$, or 2.4 c) is used. This is in order to prevent the timestep becoming prohibitively short. We chose the minimum value such that the speed of ionization fronts in our simulations is the same as that for a larger speed of light.

## 2.3 Radiative cooling

In each cell, radiative cooling and heating of hydrogen and helium is performed by the radiative transfer module of RAMSES-RT as in Rosdahl et al. (2013). For metals, we have two cooling and heating functions. The first, ‘neutral’ function $N(T)$, is the prescription of Audit & Hennebelle (2005), which includes carbon, oxygen and dust grains as well as the ambient UV background in the ISM, transitioning to the prescription of Sutherland & Dopita (1993) above $10^4$ K.

A further ‘photoionized’ cooling function $P(T)$ is included, which is a piecewise fit to Ferland (2003) (see also Osterbrock & Ferland 2006). This fit has a constant $3 \times 10^{-24}$ erg cm$^3$ s$^{-1}$ below 9000 K

and $2.2 \times 10^{-22}$ erg cm$^3$ s$^{-1}$ above $10^5$ K, interpolating between these two points inside these temperatures. The metal cooling rate is set to $N(T)$ if $N(T) > P(T)$, where a positive value indicates cooling rather than heating. Otherwise, the cooling and heating rate is given by $xP(T) + (1-x)N(T)$, where x is the hydrogen ionization fraction.

## 2.4 Uv source properties

We implement UV radiation in a similar way to Geen et al. (2015b). The source of UV photons is modelled as a point source of ionizing photons in the three photon groups. For each source, we calculate a photon energy and emission rate assuming blackbody emission. Frequency-dependent cross-sections are taken from Verner et al. (1996) via Hui & Gnedin (1997). In principle, the spectrum from an OB star will differ from a blackbody spectrum, but, in practice, we find that the exact spectrum of ionizing photons is of secondary importance provided that the number of hydrogen-ionizing photons is the same (see also Haworth et al. 2015).

We use three sources of ionizing photons in this paper, as well as control simulations with no ionizing radiation. These sources have ionizing photon emission rates $S_*$ of $10^{49}$, $10^{50}$ and $10^{51}$ s$^{-1}$. In order to compare these values to physical sources, we use the results of Vacca, Garmany & Shull (1996) and STARBURST99 (Leitherer et al. 2014). The sources are taken to be, respectively, a 40 M$_\odot$ star, a 100 M$_\odot$ star, and a cluster of 10 100 M$_\odot$ stars. Alternatively, the hydrogen-ionizing photon emission rates correspond to clusters of masses 100, 1000 and 10 000 M$_\odot$ respectively, i.e. 0.1 per cent, 1 per cent and 10 per cent of the total mass of the cloud ($10^5$ M$_\odot$). See Appendix A for a calculation of these values. Each of these sources are turned on for 3 Myr, at which point the stars enter the horizontal giant branch phase and stop producing significant quantities of ionizing photons.

## 2.5 Supernova model

After 3 Myr, when the UV photons are extinguished, we launch the supernova. This is a simplification of the full stellar lifecycle, which is beyond the scope of this paper. Instead, the intent of this
paper is to use semirealistic prescriptions that allow us to explore the behaviour of feedback in molecular clouds using controlled conditions. We discuss more sophisticated models in Section 6.

We implement the supernova as a point injection of thermal energy into the centre of the simulation volume, at the same location as the source of UV radiation. Our spatial resolution (0.084 pc around the supernova) satisfies the resolution criterion given by Kim & Ostriker (2015) (0.14 pc at $10^4$ atoms cm$^{-3}$). In our primary supernova model (labelled ‘SN’), we inject $10^{51}$ erg of energy and 1 M$_{\odot}$ of mass, representing the ejecta from the supernova. We also include a hypernova model (labelled ‘HN’; Nomoto et al. 2005), which is identical to the ‘SN’ model except that we inject $10^{52}$ erg and 10 M$_{\odot}$. A further model, labelled ‘MSN’, includes 10 supernovae with energy $10^{51}$ erg and 1 M$_{\odot}$ of ejecta, launched 0.1 Myr apart (see Kim & Ostriker 2015, for a similar model).

# 3 THE H II REGION BEFORE THE SUPERNOVA

In this section, we present the results of our simulations containing H II regions before the supernova is launched. This includes both the phase before the source of ionizing photons is turned on and the period during which the source is shining. We then review predictions of the behaviour of the ionization front made in Geen et al. (2015b) and extend these models to describe the time-dependent behaviour of the system.

## 3.1 The pre-stellar phase

The cloud initially evolves magneto-hydrodynamically under self-gravity. The cloud fragments over roughly one free-fall time. We plot the density profile at the start of the simulation and after one free-fall time in Fig. 1. The dense clumps formed after one free-fall time cover only 2 per cent or less of the solid angle around the source. For this reason, the median density profile is significantly lower than the mean density profile.

As the turbulence dissipates, the cloud begins to contract, with dense clumps falling towards the centre of the cloud. The dissipation time is of the order of the crossing time for the velocity dispersion of the cloud, a few times the free-fall time of the cloud (Lee & Hennebelle 2016b). This effect was also seen in Geen et al. (2015b).

The densest clumps fall towards the centre of the cloud. They may be prevented from doing so by radiative feedback, though the precise behaviour of the clumps is complex and difficult to quantify using a simple model. Since the densest gas in the cloud is found in these clumps, if they are able to reach the position of the source, they will quench the H II region and cause it to ‘flicker’ (Peters et al. 2010).

## 3.2 Expansion of the H II region

Once the stars begin to radiate, they ionize the gas around them. This photoionized gas has an equilibrium temperature around $10^4$ K, set by the radiative cooling function. The ionization front reaches a hydrostatic limit called the Strömgren radius, proportional to $S_*^{1/3}$, where all of the photons are used to keep the gas photoionized. Due to a pressure difference with the neutral gas, the ionization front expands while maintaining photoionization equilibrium. In a uniform medium with no source of external pressure, the ionization front radius $r_i$ expands such that $dr_i/dt$ is proportional to $S_*^{1/7}$ (Matzner 2002).

Eventually, the front ‘stalls’ (i.e. is unable to expand further) due to ram pressure from external turbulence. In Geen et al. (2015b), we

present a model that estimates the radius at which this occurs, $r_{\text{stall}}$, for a virialized cloud with radius $r_{\text{cloud}}$ (see Appendix B). We find three regimes governed by the ratio $r_{\text{stall}}/r_{\text{cloud}}$. If this ratio is much larger than 1, the cloud is dispersed. If it is much smaller than 1, the ionization front is trapped by the cloud. If the ratio is close to 1, the ionization front escapes as a blister region but does not completely destroy the cloud.

The value of $r_{\text{stall}}$ depends on the density profile of the cloud. In Geen et al. (2015b), we fit a single spherically averaged radial power law to the density field in the simulated cloud, in addition to sampling the density of the simulated cloud as a function of radius and time. The single power-law fit provided an adequate match to the simulation results. However, in this paper, where we use a more massive cloud, we find that we must adopt the two-phase fit to the density profile of the cloud given by Franco et al. (1990), with a flat cloud core surrounded by a power-law density field. This is given by

$$n_{\text{ext}}(r \le r_c) = n_c \\ n_{\text{ext}}(r > r_c) = n_c(r/r_c)^{-w}, \qquad (1)$$

where $r_c = 3.6$ pc is the scaling radius given in the initial conditions and $n_c$ is the density at $r = r_c$. These values are found by fitting the spherically averaged mean density field outside $r_c$ to a power law with free parameters in $n_c$ and $w$ at the time the source is turned on. The fit gives $n_c = 1612$ atoms cm$^{-3}$ and $w = 2.93$, i.e. with a very steep transition from the cloud core to the diffuse medium outside.

There is some flexibility in deciding the optimal fit for the full 3D density profile. The fit to equation (1) is degenerate depending on the value of $r_c$ chosen. In addition, in cases where a large quantity of mass is in small clumps sufficiently far from the position of the source, the median density profile offers a better fit to the effective density field experienced by the H II region.

Since the power-law slope at $r > r_c$ is so steep, we argue that if $r_{\text{stall}}$ exceeds $r_c$, the ionization front is able to escape the cloud. This is equivalent to our limit comparing $r_{\text{stall}}$ and $r_{\text{cloud}}$ in Geen et al. (2015b).

## 3.3 Comparison to the simulations

We calculate the ratio $r_{\text{stall}}/r_c$ for each of our simulations. The values for N49-NSN, N50-NSN and N51-NSN are 0.82, 1.1 and 1.6. Each simulation is thus in (or close to) a regime in which we predict the front should stall, almost escape the cloud or expand more or less freely. Again, there is some error in these estimates depending on the quality of the fit of a spherical density profile to the simulations.

In Fig. 2, we repeat the calculation performed in fig. 8 of Geen et al. (2015b), solving equation (B1) numerically using the density field given in equation (1) to estimate the expansion of the ionization front. This model includes ram pressure from turbulent motions in the cloud, which we assume to be virialized. In Geen et al. (2015b), we set $v_{\text{ext}}$ to the escape velocity whereas here we use $v_{\text{ext}} = -v_{\text{vir}}(r) = -\sqrt{\frac{6GM(<r)}{5r}}$, which is 77 per cent of the escape velocity. Using the escape velocity instead results in a smaller value for $r_{\text{stall}}$. This model is labelled ‘Analytic’.

We also sample the spherically averaged density and radial velocity field at each time and radial position in simulation N00-NSN and solve equation (B1) using these inputs. We label this solution ‘Sampled’. We plot both these models against the median radius of the ionization front in our simulations at each timestep. Gas in the simulations is assumed to be ionized if its hydrogen ionization fraction is above 0.1.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparison of our simulation results (solid lines) and model predictions (dashed lines) for the median radial expansion of the ionization front with various photon emission rates. The left-hand plot labelled 'Analytic' uses the density profile fit described in equation (1). The plot on the right labelled 'Sampled' samples the density and radial velocity field in simulation N00-NSN to predict the behaviour of the ionization front. See Section 3.3 for more details.</figcaption></figure>

The behaviour of the ionization fronts in each of the simulations agrees well with the ‘Sampled’ models, with the exception of the $10^{50} \text{ s}^{-1}$ simulation, whose ionization front does not collapse as the model predicts. This is due to the fact that the dense clumps in this simulation are prevented from reaching the centre of the cloud, and hence the photons are able to escape over the majority of the lines of sight around the source.

The expansion in the ‘Analytic’ model follows the broad features of the ionization fronts with $10^{50}$ and $10^{51} \text{ s}^{-1}$ sources as they expand. The expansion of the Analytic model is too fast in the $10^{49}$ and $10^{50} \text{ s}^{-1}$ sources, suggesting either that our choice of density is too low or that the velocity of the gas is higher than our estimate. In addition, since we assume a static density field we do not capture the collapse of the ionization front in the $10^{49} \text{ s}^{-1}$ case. We present a term correcting for this collapse in Appendix D.

The equations governing the expansion of the ionization front in a turbulent medium are non-linear. For this reason, in Geen et al. (2015b), we solve these equations numerically. An analytic limit at which the front stalls is also given. It is useful to know over what time this stall radius is reached. In Appendices C and D, we show that this time-scale over which the front either stalls (or collapses, if the cloud is strongly accreting) is roughly equal to the free-fall time in the cloud core. This suggests that star formation in the cloud and the cloud destruction occur over similar time-scales, making it difficult to find an accurate analytic model for the point at which star formation is frozen out by the destruction of the cloud via radiative feedback.

In general, the models obey the broad behaviour defined by our limit in $r_{\text{stall}}$, with some uncertainty due to the complex structure of the cloud as the ratio $r_{\text{stall}}/r_c$ approaches 1 and the stalled front transitions to a freely expanding front. For the cloud simulated here, this means that the ionizing photons are only able to significantly disperse the cloud with a $10^{51} \text{ s}^{-1}$ source, though the $10^{50} \text{ s}^{-1}$ source is able to drive a large ionized bubble that escapes the cloud in certain directions but does not escape the simulation volume.

### 3.4 After the source is extinguished

After a few million years, the star no longer produces a significant amount of ionizing photons. This causes the H II region to cool on a time-scale governed by the density of the gas inside the H II region.

We discuss a simple model for this in Geen et al. (2015a), where the ionization front is no longer in ionization equilibrium. Once the source of pressure from photoheated gas is reduced, the shell around the ionization front can continue to expand due to momentum conservation, as derived by Hosokawa & Inutsuka (2006). If the front is stalled or collapses, the residual momentum will be close to zero. Indeed, once the source is turned off we observe a gradual re-collapse of the cloud.

The lowest density reached in the ionized gas in any simulation is $10 \text{ atoms cm}^{-3}$. The recombination time in fully ionized gas at a temperature of $10^4 \text{ K}$ with this density is 12 kyr, while the radiative cooling time is 2 kyr. Both time-scales are inversely proportional to the gas density. These time-scales are typically much shorter than the lifetime of the system. As such, when the most massive star in the cloud is extinguished, the ionized hydrogen will rapidly recombine.

Our results agree with the conclusions of Geen et al. (2015b). If the source is strong enough to resist the infall of clumps, the ionization front can expand into the external medium. Otherwise, it stalls or contracts. The boundary between these two cases is found at around $r_{\text{stall}} = r_c$.

## 4 AFTER THE SUPERNOVA

So far, we have discussed the evolution of the cloud and H II region up to the point where the source is extinguished, equivalent to the point at which the most massive star in the cluster reaches the end of its life. In this section, we discuss the phase of our simulations after this star goes supernova.

### 4.1 The structure of the cloud before and after the supernova

We first review the properties of the cloud at the point where the supernova is injected for the different photon emission rates where $r_{\text{stall}}/r_c$ is less than, roughly equal to and greater than 1.

In Fig. 3, we plot the density profiles around the source just prior to the supernova. Due to turbulent dissipation in the cloud and the infall of dense clumps, the profile becomes more peaked over time, particularly with the $10^{49}$ photons s$^{-1}$ source. Only the $10^{51} \text{ s}^{-1}$ source is able to significantly flatten the cloud profile. Similarly, as described in Section 3.3, only the $10^{51} \text{ s}^{-1}$ source is able to drive
<figure><img src="image_5.png" /><figcaption>Figure 3. Density profiles just before the supernova is launched for each of the photon emission rates. Top left: no source. Top right: $10^{49}$ photons s$^{-1}$. Bottom left: $10^{50}$ photons s$^{-1}$. Bottom right: $10^{51}$ photons s$^{-1}$. The plots are constructed as in Fig. 1.</figcaption></figure>

ionized flows out of the simulation volume, while the ionized gas around the weaker sources remains bounded by neutral gas. This means that the supernova shock must break through this material if it is to move out into the external medium.

Fig. 4 contains hydrogen column density maps before and after the supernova. As with the ionization front, the supernova remnant expands preferentially through lower density channels in the cloud. Since the $10^{51}$ s$^{-1}$ source was the only source capable of expelling a significant quantity of the cloud material – 80 per cent as a function of solid angle around the source – this is the only simulation in which the supernova successfully escapes the simulation volume and into the surrounding medium. In Fig. 5, we plot a slice through the temperature field around the source. Since the supernova explodes in dense gas, the remnant has cooled to around $10^4$ K after a few hundred kyr. Even in the $10^{51}$ s$^{-1}$ simulation, the temperature is reduced by reflection shocks of denser gas passing through the hot, diffuse bubble due to the asymmetric shape of the remnant and the recollapse of the cloud under gravity.

## 4.2 Analytic overview

The classical picture of the expansion of supernova remnants in a uniform medium is described by Chevalier (1977). Once the supernova blastwave leaves the surface of the star, there is an initially

kinetic phase as the ejecta moves outwards ballistically. Once the mass displaced by the supernova is roughly equal to the supernova ejecta mass, the supernova ejecta shocks against the surrounding medium and enters the adiabatic, ‘Sedov’ phase (Sedov 1946). As the remnant expands, it sweeps up the unshocked matter, forming a shell around the shock-heated gas. It also begins to cool radiatively, entering the pressure-driven snowplough phase. As radiative cooling becomes significant, the thermal pressure inside the shell drops to the point where the shell expands only due to momentum conservation in the momentum-conserving snowplough phase.

An equation for the time at which radiative cooling is complete and shell radius at which this occurs is given in section c (ii) of Cox (1972) as:

$$ 
\begin{align} 
t_{\Lambda} &= 5.0 (\epsilon_0)^{4/17} n_0^{-9/17} \times 10^4 \text{ yr}, \\ 
r_{\Lambda} &= 22.1 (\epsilon_0)^{5/17} n_0^{-7/17} \text{ pc}, 
\end{align} 
\qquad (2) 
$$

where $\epsilon_0$ is a unitless quantity given by the energy of the supernova divided by $10^{51}$ erg and $n_0$ is a unitless quantity given by the density in the surrounding medium divided by 1 atoms cm$^{-3}$ (assumed to be uniform). We find that for values typical in our simulation, the supernova should exit the Sedov phase before it leaves the core of the cloud, with $r_{\Lambda} = 3.3$ pc at $n_0 = 100$ and $r_{\Lambda} = 0.50$ pc at $n_0 = 10000$.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Hydrogen column density maps for each of the simulations containing supernovae at the time the supernova is injected (5.52 Myr after the start of the simulation, and 3 Myr after the source of ionizing radiation is turned on) and 1 Myr afterwards. The cyan line shows the maximum extent of ionized gas along the line of sight. This gas can be either photoionized in the case of the H II region or collisionally ionized in the case of the hot bubble inside the supernova remnant. The red dot shows the position of the source of ionizing radiation where one is included (note that this source is turned off at the times shown in these images). The total length of image is 86.3 pc in all dimensions.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** As in Fig. 4 but showing gas temperature as a slice through the centre of the simulation volume.</figcaption></figure>

Iffrig & Hennebelle (2015) find that an accurate prediction of the final momentum of the supernova remnant $p_Δ$ can be achieved by calculating the momentum of the remnant in the Sedov phase at $t_Δ$, where

$$p_Δ \propto n_0^{-2/17} \epsilon_0^{16/17}. \qquad (3)$$

In other words, the final momentum is weakly linked to density (around 1 to $3 \times 10^{43}$ g cm s$^{-1}$ in our cloud), and roughly proportional to the initial supernova energy. Iffrig & Hennebelle (2015) find that this simple argument compares well to simulations in both

uniform and turbulent molecular clouds. In the following section, we discuss how this compares to the simulations in this paper.

## 4.3 Momentum from the supernova blastwave

The supernova adds around $10^{43}$ g cm s$^{-1}$ of radial momentum to the system, roughly an order of magnitude lower than the total momentum in flows in the cloud (see Fig. 6). This estimate is made by subtracting the momentum in each simulation with a supernova from the momentum in an identical simulation without a
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Momentum in radial flows in simulations with and without supernova feedback. On the left is a plot showing the total radial momentum with and without a supernova included, as well as the hypernova ('HN') and multiple supernova ('MSN') phase. The right-hand plot shows estimates of the momentum added by the supernova blaswaves in each simulation. These are made by taking the momentum at a given time in a simulation containing one or more supernovae or hypernova, and subtracting the momentum at the same time in an analogous run without a supernova. The time in the right-hand plot uses a log scaling to emphasize the early evolution of the blastwave. The $10^{51}$ photons s$^{-1}$ results are taken from a simulation with twice the box length of the other simulations in order to prevent losing mass and momentum from the supernova shocks.</figcaption></figure>

supernova. Nonlinearities will add some error to the precise value found.

Most of the momentum in the system is found in gas above $100$ atoms cm$^{-3}$, with the exception of the $10^{51}$ photons s$^{-1}$ simulation, which has succeeded in destroying the cloud. When ionizing radiation is added, the supernova is able to add some momentum in the phase between $10$ and $100$ atoms cm$^{-3}$. Since the supernova leaves the adiabatic phase long before it escapes the cloud, its role is largely to accelerate clumps of gas away from the cloud rather than driving hot, diffuse winds out of the cloud.

We compare our results to estimates given in previous papers assuming the density at the position of the supernova, though there are strong density gradients around the position of the cloud (see Fig. 3). Estimates by Iffrig & Hennebelle (2015) and Kim & Ostriker (2015) give somewhere between $1.1$ and $1.7 \times 10^{43}$ g cm s$^{-1}$ of momentum for a cloud of density $10^3-10^4$ atoms cm$^{-3}$, the densities found in the cloud around the supernova in all runs except when a $10^{51}$ photons s$^{-1}$ source is included. Here, the density at the position of the source is closer to $100$ atoms cm$^{-3}$, although even here we find energy losses due to interactions with dense cloud material as the supernova remnant expands. By contrast, Iffrig & Hennebelle (2015), who use a similar setup with a $10^4$ M$_\odot$ cloud, find that their results agree well with simulations performed in a uniform box with negligible external pressure forces.

Our momentum is lower than the values found by Walch & Naab (2015), who use a similar cloud mass, although our cloud is more fragmented and turbulent. We posit that this is because we pre-evolve our cloud with an initial turbulent field, whereas Walch & Naab (2015) impose a fractal density field with zero initial velocity field. Both our simulations and Walch & Naab (2015) include gravity. Tests performed in a uniform medium reproduce other works more closely, although a small difference is found due to the use of non-equilibrium hydrogen and helium cooling (Rosdahl et al., in preparation).

We suggest that the lower momentum found in our simulations is due to the density and velocity structure of the cloud. Our cloud fragments into dense, turbulent clumps as it evolves. These clumps have a large amount of ram pressure that opposes shocks interacting with

them. For example, a parcel of fluid with density $10^4$ atoms cm$^{-3}$ moving at $5$ km s$^{-1}$ provides $3 \times 10^{-9}$ erg cm$^{-3}$ in ram pressure. The thermal pressure in the hot bubble of the supernova remnant is likely to be lower than this: $2 \times 10^{-9}$ erg cm$^{-3}$ for a monatomic gas at $10^7$ K and $1$ atoms cm$^{-3}$, which is a very high estimate for the density inside the hot bubble. In addition to this, McKee & Ostriker (1977) argue that the evaporation of dense clumps inside the hot bubble further reduces the energy of the supernova remnant.

Cioffi et al. (1988) argue that at the end of the supernova remnant’s life, the shock will merge with the surrounding material. Equation (4.5) in that paper gives the time that this happens as

$$t_{\text{merge}} = 7.6n_{0}^{-18/49}\epsilon_{0}^{31/98} \times 10^{5} \text{yr}, \quad (4)$$

assuming the characteristic velocity of the turbulence is around $10$ km s$^{-1}$, which is to within an order of magnitude the value in our cloud. For $n_0 = 100$ and $10\,000$ (in units of atoms cm$^{-3}$ as in equation 2), $t_{\text{merge}} = 0.14$ and $0.026$ Myr, respectively. This is consistent with the transition in our simulations where the added radial momentum from the supernova remnant begins to fluctuate as radial flows are transferred to turbulent motions in dense clumps. The only simulation where this does not happen is in the $10^{51}$ photons s$^{-1}$ simulation, where the cloud prior to the blast is highly porous (see Fig. 4).

We discount the possibility that the shell around the supernova remnant decelerates due to gravity from material still embedded inside the cloud after the blastwave has expanded. In Appendix E, we calculate that this should only become important for supernova remnants moving through densities above $10^3$ atoms cm$^{-3}$. Since the majority of the cloud is at a lower density, we expect this to only be a secondary effect for this system.

## 5 BEYOND $10^{51}$ ERG

The previously discussed results were for the case where a single supernova with energy $10^{51}$ erg was injected. In this section, we discuss two simple extensions to this model that include both a rare but powerful hypernova explosion, and a simple model of multiple supernovae exploding in the same cluster.
## 5.1 Hypernovae

According to Heger et al. (2003), above around $25 M_\odot$ stars become black holes rather than neutron stars, leading to weak supernovae. However, Nomoto et al. (2005) find that a small fraction of very massive stars with high rotation rates will explode as hypernovae, which inject 10 times or more the energy of typical supernovae. In this section, we discuss a simulation with a $10^{50}$ s$^{-1}$ source where we inject a hypernova of energy $10^{52}$ erg into the cloud. In Fig. 6, we find 10 times the momentum from the hypernova as the standard supernova model, agreeing with equation (3). In this case, the blast adds sufficient momentum that the signal is not lost in the turbulent motions of the cloud.

The total momentum in flows balancing the gravitational forces in the cloud is approximately the cloud mass multiplied by the virial velocity at the cloud edge. For a cloud of $10^5 M_\odot$ with a total radius of 20 pc (see Fig. 3), this is roughly $10^{44}$ g cm s$^{-1}$, approximately the amount of momentum deposited by the hypernova. A hypernova thus deposits sufficient momentum to counteract the gravitational forces binding the cloud. In practice, due to the clumpy nature of the medium and imperfect coupling of the hypernova blastwave to the clumps, we find that some cloud material is able to remain embedded in the cloud.

## 5.2 Multiple supernovae

In the cloud of $10^5 M_\odot$, there will typically be multiple massive stars in the embedded cluster. There will hence be multiple supernovae that explode in the same cloud. In this section, we analyse a simulation in which we inject nine further $10^{51}$ erg supernovae every 0.1 Myr after the first supernova. This is similar to the experiment performed by Kim & Ostriker (2015) in a turbulent environment. Each time a supernova is injected, it occurs in an environment swept out by the previous supernovae, gradually inflating the supernova remnant in short bursts.

In Fig. 6, we find that the momentum injected by these supernovae gives roughly the same end result as the hypernova. Again, this is consistent with equation (3), in which the final momentum is roughly proportional to the energy of the blast. We do not find the large increase in momentum with subsequent supernovae as Gentry et al. (2016) do in their 1D study, despite reaching a comparable spatial resolution (0.084 pc compared to 0.06 pc). Part of this is due to the fact that dense material remains embedded in the cloud even after the initial supernova. However, the claim in Gentry et al. (2016) that Eulerian codes suffer from overmixing near the shell boundary should be carefully studied, particularly in combination with 3D turbulence and full non-equilibrium radiative cooling.

For a cluster with mass $10^4 M_\odot$ (10 per cent of the total mass of the cloud), using an initial mass function (IMF) slope of -2.35 in the high-mass end, we estimate 160 stars above $8 M_\odot$. In the range of 4–20 Myr, this gives a mean time delay of 0.1 Myr, although due to the shape of the IMF, more stars will explode at later times. We thus expect this model to be an underestimate for the total injection of energy by supernovae into the cloud by a factor of 16. However, if stars form over around $t_{ff}$ in the cloud (2.53 Myr), these supernovae will occur too late to prevent star formation. Rather, they will be more effective at dispersing any remaining cloud material.

# 6 DISCUSSION

In this work, we have focused on the combined role of ionizing radiation and supernovae in destroying molecular clouds, and the

interaction between the two processes. However, there are a number of important aspects still to be explored in both reproducing and explaining star formation and feedback in molecular clouds.

A number of physical processes have not been included in this work. Radiation pressure is omitted, though there is some debate in the literature regarding its effectiveness. Agertz et al. (2013) and Hopkins et al. (2014) assume strong coupling between infrared photons and dust grains, while Krumholz & Thompson (2012), Sales et al. (2014), Rosdahl et al. (2015) and Haworth et al. (2015) argue that radiation pressure has a minimal effect in the regimes studied.

Stellar winds too are omitted. Dale et al. (2014) argue that winds are not likely to be as effective at driving outflows as ionizing radiation, although these authors neglect late stage stellar evolution where winds become stronger. This source of late-stage feedback is particularly important because the amount of UV photons drops as the star expands and cools. Corrections by Kudritzki & Puls (2000) for the wind velocity boost wind energies by up to a factor of 10 over values assuming a wind travelling at the escape velocity of the star. By combining recent models of stellar evolution by Ekström et al. (2012) with the boost from Kudritzki & Puls (2000), cumulative energies from stellar winds can exceed that from a supernova, though this energy will be spread out across the lifetime of the star rather than as a single burst of energy as in a supernova.

Mason et al. (2009) and de Mink et al. (2014) find that the majority of observed massive stars are in binaries. As such, many of them are rotating and winds may become more efficient in this regime. This means that as these stars lose their envelopes, they stay hot in the final stage of their evolution, producing more UV photons (Köhler et al. 2014).

The model for supernovae becomes more complex for very massive stars. Heger et al. (2003) and Nomoto et al. (2005) discuss the fate of stars above $25 M_\odot$, which can become either weak supernovae or very energetic hypernovae. Podsiadlowski et al. (2004) argue that the latter events are rare, though since they can deposit significantly more energy into the ISM than $10^{51}$ erg, they may be significant events, as we find in Section 5.1.

Finally, processes such as protostellar jets (see review by Frank et al. 2014) and X-ray emission from stellar-mass black holes (Mapelli & Zampieri 2014) can also increase the energy budget for stellar feedback in the cloud. For a review of stellar feedback processes in star-forming regions, see Dale (2015).

The time over which the star emits energy in various forms depends strongly on the stellar evolution models used. In this paper, we assume a lifetime of 3 Myr followed by a supernova for the most massive stars. We omit a period of around 1 Myr during which stellar winds, which we do not include in this paper, become important. Binary evolution, which affects a majority of observed very massive stars, also affects the lifetime of stellar winds and ionizing radiation from the stars. There are, thus, many open questions regarding how much energy is available from feedback in various forms, and over which time-scale.

Star formation in the cloud and the response of the cloud to radiation both occur on a scale of a free-fall time. There is hence a competition between the two processes that will, in part, set the star formation efficiency of the cloud. Since, in e.g. Matzner (2002), the expansion of the ionization front is proportional to $S_*^{1/7}$, the cloud is not highly sensitive to the precise number of photons produced by the cluster. Additionally, stars form in very dense regions that are underresolved by our simulations, rather than at the centre of the cloud as in this work. We will begin to address these question in future work.
In our simulations, we find that a $10^{51}$ erg supernova will add $10^{43}$ g cm s$^{-1}$ of momentum to a $10^5$ M$_\odot$ cloud, mostly in the dense phase. Adding more energy, either in a single hypernova or as multiple supernovae, adds proportionally more momentum. Eventually, enough momentum is added to unbind the cloud.

The role of supernovae in setting the star formation efficiency of the cloud is unclear. Since the first supernovae occur a few Myr after the first star is formed, they cannot immediately regulate star formation in the cloud. In addition, one supernova will not be enough to unbind a massive cloud. By contrast, ionizing radiation is capable of disrupting star-forming clouds, particularly, in the case studied here where a cluster producing $10^{51}$ photons s$^{-1}$ unbinds the entire cloud prior to the first supernovae. In this scenario, the supernovae will be injected directly into the diffuse medium.

There are two advantages to invoking supernova feedback in a cloud environment as opposed to ionizing radiation. The first is that they couple directly to the gas, and as such are more efficient at transferring their energy to the gas than ionizing photons. Walch et al. (2012) estimate an efficiency of energy from ionizing photons to kinetic energy in the gas as approximately 0.1 per cent, compared to the few per cent found by e.g. Chevalier (1974) for supernovae. The second is that they are capable of driving outflows at high velocities, whereas photoionized gas can only expand at around $10$ km s$^{-1}$.

Different environments will have different behaviours. For example, high-redshift H II regions will have higher temperatures since metal cooling is absent, but they will be embedded in denser, higher pressure environments. Very massive clouds will also be difficult to disperse. For example, Krause et al. (2016) find that only a large number of hypernovae are able to prevent the majority of a massive extragalactic cloud from turning into stars. More work must be done to extend the findings of this paper in local environments to more universal conditions.

Even if supernovae are less important in dispersing molecular clouds than ionizing radiation, provided that star formation continues after 3 Myr, they will begin to metal-enrich future generations of stars in the system. This has been proposed for our own Solar system by e.g. Cameron & Truran (1977) and Looney, Tobin & Fields (2006). A controlled study of the triggering of star formation by supernova blastwaves remains beyond the scope of this paper, since the system is highly nonlinear and we lack self-consistent star formation. Dale, Ercolano & Bonnell (2013) find that it is difficult to identify protostar formation as being 'triggered', as opposed to forming via gravitational collapse. A cursory analysis of our cloud would suggest that the same dense clumps form stars regardless of their interaction with the H II region, provided that they are not photoevaporated first. However, Boss et al. (2008), Gounelle et al. (2009) and Gritschneder et al. (2011) argue for a scenario where only supernova shocks are able to cause collapse and chemical mixing of protostars at a rate that can explain observed chemical abundances in the Solar system.

## 7 CONCLUSIONS

We perform a series of simulations in order to study how ionizing UV feedback from young star clusters is able to pre-process molecular cloud environments prior to the first supernova explosions. We simulate a $10^5$ M$_\odot$ cloud with embedded sources of ionizing radiation and a supernova. We compare our results to the analytic arguments in Paper I (Geen et al. 2015b) for the expansion of ionization fronts, and existing analytic models for supernova

blastwaves. The goal of this comparison is to provide insights into the behaviour of these systems with both simple, illustrative models and full radiative magnetohydrodynamic numerical simulations.

We find good agreement with the model from Paper I provided a good fit is found for the density field of the cloud. We find that the limit at which the ionization front stalls (i.e. stops expanding) is reached over approximately one free-fall time in the cloud. This introduces a competition between the two processes of star formation and stellar feedback. The expansion of H II regions in clouds is strongly dependent on the velocity field in the cloud, and it is thus crucial to capture the turbulent and infalling velocity profiles of the cloud gas if we are to understand how feedback regulates the cloud structure. Once the ionizing radiation is extinguished, the H II rapidly recombines and the remaining dense cloud gas begins to recollapse.

The environment that the supernova blastwave expands into depends strongly on the emission rate of ionizing photons from the cluster beforehand. An emission rate of $10^{51}$ photons s$^{-1}$, roughly equivalent to a cluster of $10^4$ M$_\odot$ or 10 per cent of the total cloud mass, is capable of disrupting the cloud, though dense clumps remain. If the number of stars formed is much lower, the ionizing photons will not be able to destroy the cloud and the supernova will transfer its momentum to the dense cloud gas rather than the diffuse interstellar medium. This suggests that for this cloud, a star formation efficiency of approximately 10 per cent is expected if the main source of feedback is from ionizing photons. The position of the source of ionizing photons and the supernova is highly important due to the effect of gas density on photon recombination and radiative cooling in general.

We inject supernovae as a single thermal pulse of $10^{51}$ erg. We also perform simulations with ten supernovae of the same energy 0.1 Myr apart or one hypernova of $10^{52}$ erg. The resulting total momentum from our supernovae is roughly $10^{43}$ g cm s$^{-1}$ per $10^{51}$ erg of injected energy. This is at the low end of the values given in other works, but it is not inconsistent provided the early evolution of the blastwave occurs in gas at around $10^4$ atoms cm$^{-3}$ or higher. Most of the momentum from a single supernova is deposited into the dense gas rather than as fast, hot, diffuse flows, except in cases where the ionizing photons have first swept away most of the cloud. Using analytic argument from previous authors, we posit that the strong turbulence in the cloud causes the radial momentum from the supernova in the post-Sedov phase to be absorbed by turbulent flows.

Supernovae will occur too late to prevent the bulk of star formation in the cloud. The emission rate of ionizing photons remains the most important quantity for determining whether stellar feedback can destroy the host cloud or not. However, over longer time-scales, multiple supernovae from a cluster remain an important source of momentum for flows in the ISM.

## ACKNOWLEDGEMENTS

The authors thank the referee for their useful comments that helped improve the clarity and accuracy of the manuscript. We would like to thank Olivier Iffrig, Romain Teissier, Andreas Bleuler, Alex Richings, Suzanne Madden, Yueh-Ning Lee, for useful discussions during the preparation of this paper. The simulations presented here were run on the machine irfucoast at CEA Saclay. This work has been funded by the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013). SG and PH are funded by Grant Agreement no. 306483 of this programme. JR is funded
by Grant Agreement 278594-GasAroundGalaxies of the same programme and the Marie Curie Training Network CosmoComp (PITN-GA-2009-238356).

## REFERENCES

Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

Audit E., Hennebelle P., 2005, A&A., 433, 1

Boss A. P., Ipatov S. I., Keiser S. A., Myhill E. A., Vanhala H. A. T., 2008, ApJ, 686, L119

Cameron A., Truran J., 1977, Icarus, 30, 447

Chabrier G., 2003, PASP, 115, 763

Chevalier R. A., 1974, ApJ, 188, 501

Chevalier R. A., 1977, ARA&A, 15, 175

Cioffi D. F., McKee C. F., Bertschinger E., 1988, ApJ, 334, 252

Cox D. P., 1972, ApJ, 178, 159

Dale J. E., 2015, Phys. Rev. B, 68, 1

Dale J. E., Ercolano B., Bonnell I. A., 2012, MNRAS, 424, 377

Dale J. E., Ercolano B., Bonnell I. A., 2013, MNRAS, 431, 1062

Dale J. E., Ngoumou J., Ercolano B., Bonnell I. A., 2014, MNRAS, 442, 694

de Mink S. E., Sana H., Langer N., Izzard R. G., Schneider F. R. N., 2014, ApJ, 782, 7

Didelon P. et al., 2015, A&A, 584, A4

Draine B. T., Woods D. T., 1991, ApJ, 383, 621

Dwarkadas V. V., 2007, ApJ, 667, 226

Ekström S. et al., 2012, A&A, 537, A146

Ferland G. J., 2003, ARA&A, 41, 517

Fierlinger K. M., Burkert A., Ntormousi E., Fierlinger P., Schartmann M., Ballone A., Krause M. G. H., Diehl R., 2015, MNRAS, 456, 710

Franco J., Tenorio-Tagle G., Bodenheimer P., 1990, ApJ, 349, 126

Frank A. et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. Univ. Arizona Press, Tucson

Fromang S., Hennebelle P., Teyssier R., 2006, A&A, 457, 371

Garcia-Segura G., Franco J., 1996, ApJ, 469, 171

Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015a, MNRAS, 448, 3248

Geen S., Hennebelle P., Tremblin P., Rosdahl J., 2015b, MNRAS, 454, 4484

Gentry E. S., Krumholz M. R., Dekel A., Madau P., 2016, MNRAS, preprint (arXiv:1606.01242)

Gounelle M., Meibom A., Hennebelle P., Inutsuka S. I., 2009, ApJ, 694, L1

Gritschneder M., Lin D. N. C., Murray S. D., Yin Q. Z., Gong M. N., 2011, ApJ, 745, 22

Haid S., Walch S., Naab T., Seifried D., Mackey J., Gatto A., 2016, MNRAS, 460, 2962

Haworth T. J., Harries T. J., Acreman D. M., Bisbas T. G., 2015, MNRAS, 453, 2278

Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ, 591, 288

Hennebelle P., Falgarone E., 2012, A&AR, 20, 55

Hopkins P. F., Keres D., Onorbe J., Faucher-Giguere C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Hosokawa T., Inutsuka S., 2006, ApJ, 646, 240

Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27

Iffrig O., Hennebelle P., 2015, A&A, 576, A95

Kahn F. D., 1954, Bull. Astron. Inst. Neth., 12, 187

Keto E., 2002, ApJ, 580, 980

Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99

Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900

Köhler K. et al., 2014, A&A, 573, A71

Körtgen B., Seifried D., Banerjee R., Vázquez-Semadeni E., Zamora-Avilés M., 2016, MNRAS, 459, 3460

Krause M. G. H., Charbonnel C., Bastian N., Diehl R., 2016, A&A, 587, A53

Krumholz M. R., Thompson T. A., 2012, ApJ, 760, 155

Kudritzki R.-P., Puls J., 2000, ARA&A, 38, 613

Larson R. B., 1969, MNRAS, 145, 271

Lee Y.-N., Hennebelle P., 2016a, A&A, 591, A30

Lee Y.-N., Hennebelle P., 2016b, A&A, 591, A31

Leitherer C., Ekström S., Meynet G., Schaerer D., Agienko K. B., Levesque E. M., 2014, ApJS, 212, 14

Looney L. W., Tobin J. J., Fields B. D., 2006, ApJ, 652, 1755

McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148

Mapelli M., Zampieri L., 2014, ApJ, 794, 7

Martizzi D., Faucher-Giguere C.-A., Quataert E., 2015, MNRAS, 450, 504

Mason B. D., Hartkopf W. I., Gies D. R., Henry T. J., Helsel J. W., 2009, AJ, 137, 3358

Matzner C. D., 2002, ApJ, 566, 302

Nomoto K., Maeda K., Tominaga N., Ohkubo T., Deng J., Mazzali P. A., 2005, Ap&SS, 298, 81

Ntormousi E., Hennebelle P., 2015, A&A, 574, A130

Osterbrock D. E., Ferland G. J., 2006, in Osterbrock D. E., Ferland G. J., eds, Astrophysics of gaseous nebulae and active galactic nuclei, 2nd. ed. University Science Books, Sausalito, CA

Peters T., Banerjee R., Klessen R. S., Low M.-M. M., Galván-Madrid R., Keto E. R., 2010, ApJ, 711, 1017

Podsiadlowski P., Mazzali P. A., Nomoto K., Lazzati D., Cappellaro E., 2004, ApJ, 607, L17

Raga A. C., Canto J., Rodriguez L. F., 2012, MNRAS, 419, L39

Rogers H., Pittard J. M., 2013, MNRAS, 431, 1337

Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380

Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188

Rosdahl J., Schaye J., Teyssier R., Agertz O., 2015, MNRAS, 451, 34

Sales L. V., Marinacci F., Springel V., Petkova M., 2014, MNRAS, 439, 2990

Sedov L. I., 1946, J. Appl. Math. Mech., 10, 241

Spitzer L., 1978, Physical Processes in the Interstellar Medium. Wiley-Interscience, New York

Sternberg A., Hoffmann T. L., Pauldrach A. W. A., 2003, ApJ, 599, 1333

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Teyssier R., 2002, A&A, 385, 337

Thornton K., Gaudlitz M., Janka H.-T., Steinmetz M., 1998, ApJ, 500, 95

Tremblin P. et al., 2014, A&A, 568, A4

Vacca W. D., Garmany C. D., Shull J. M., 1996, ApJ, 460, 914

Verner D. A., Ferland G. J., Korista K. T., Yakovlev D. G., 1996, ApJ, 465, 487

Walch S., Naab T., 2015, MNRAS, 451, 2757

Walch S. K., Whitworth A. P., Bisbas T., Wünsch R., Hubber D., 2012, MNRAS, 427, 625

Weidner C., Kroupa P., Bonnell I., 2009, MNRAS, 401, 275

Whitworth A., 1979, MNRAS, 186, 59

Williams J. P., McKee C. F., 1997, ApJ, 476, 166

## APPENDIX A: EMISSION FROM CLUSTERS

In this section, we estimate the approximate hydrogen-ionizing photon emission rate from star clusters of various masses. We sample stellar masses for a set of clusters of various total masses using a random Monte Carlo sampling. We do not a priori force the maximum stellar mass to be below a maximum value given by e.g. Weidner, Kroupa & Bonnell (2009). Ionizing hydrogen fluxes for each stellar mass are found by interpolating the results of Sternberg, Hoffmann & Pauldrach (2003), though the results do not differ strongly from values found using earlier work by Vacca et al. (1996). We plot these results in Fig. A1, along with a fit assuming that the cluster is a perfectly sampled stellar population. Below a few thousand M☉, the IMF is incompletely sampled, and as such statistical noise begins to cause a large spread in the results. However, the linear fit is still reasonable given the large spread of photon emission rates.
<figure><img src="image_9.png" /><figcaption>**Figure A1.** Hydrogen-ionizing photon emission rate for clusters of various masses. The scatter plot shows values calculated by sampling an IMF (e.g. Chabrier 2003) and assigning photon emission rates to each star using the values given by Sternberg et al. (2003). The size of the points is proportional to the most massive star in the cluster. The solid line is given by $S_* = kM_*$, where $M_*$ is the cluster mass, and $k$ is found by calculating $\sum S_*/\sum M_*$ for each cluster.</figcaption></figure>

## APPENDIX B: EXPANSION SOLUTION FOR THE IONIZATION FRONT

The spherically averaged expansion rate of the ionization front with radius $r_i$ at time $t$ in a turbulent cloud, as derived in Geen et al. (2015b) and based on Raga et al. (2012), is given by

$$ \frac{1}{c_i} \frac{dr_i}{dt} = F(r, t) - \frac{c_{\text{ext}}^2}{c_i^2} \frac{1}{F(r, t)} + \frac{v_{\text{ext}}(r, t)}{c_i}, \quad (\text{B1}) $$

where

$$ F(r, t) = \sqrt{\frac{n_i}{n_{\text{ext}}}} = \left(\frac{r_s}{r_i}\right)^{3/4} \left(\frac{n_c}{n_{\text{ext}}(r, t)}\right)^{1/2}. \quad (\text{B2}) $$

$c_{\text{ext}}$ is a term including the sound speed and turbulent motions in the gas just outside the shock radius and $v_{\text{ext}}$ is the velocity of the gas just outside the shock radius normal to the shock surface (assumed in 1D to be radial from the source position). $c_i$ is the sound speed in the ionized gas and $n_i$ is the density in the ionized gas. $r_s$ is the initial Strömgren radius, i.e. the radius at which the ionization front reaches equilibrium assuming a hydrostatic approximation. $n_c$ and $n_{\text{ext}}$ are as defined in equation (1).

## APPENDIX C: TIME-SCALE FOR IONIZATION FRONT EXPANSION

In this section, we provide a simplistic calculation for the typical time-scale for the expansion of the ionization front in equation (B1) as it approaches $r_{\text{stall}}$, where $\dot{r} = 0$. We make the simplifying assumption that $v_{\text{ext}} \gg c_{\text{ext}}$ (a solution assuming the reverse would be equally valid). We solve this equation for the flat core in the cloud, where the $n_{\text{ext}} = n_0$ is constant.

A full solution of these equations requires hypergeometric functions, which are difficult to interpret. Instead, we adopt the following simplistic estimate for the time-scale over which the ionization front stalls. This is taken to be the time that the solution to this equation assuming that $v_{\text{ext}} = 0$ (Spitzer 1978; Matzner 2002) reaches the same radius as $r_{\text{stall}}$. The reason we do this is that it provides a reasonable first-order estimate for the time at which $v_{\text{ext}}$ becomes a

<figure><img src="image_10.png" /><figcaption>**Figure C1.** Numerical solutions for the radial evolution of the ionization front in equation (B1) assuming cloud properties as given in Section 3. Overplotted are values for $r_{\text{stall}}$ and $t_{\text{stall}}$. Note that since equation (C1) makes the assumption that $r_s$ is reached over a negligible time, the match between the numerical solution and the intersection of $r_{\text{stall}}$ and $t_{\text{stall}}$ is not exact.</figcaption></figure>

limiting factor in the expansion of the ionization front. We compare this time to numerical solutions to equation (B1) in Fig. C1.

If $v_{\text{ext}} = 0$, and assuming the ionization front rapidly reaches $r_s$ (i.e. the recombination time is negligible), we can write

$$ r_i(t) = r_s \left( \frac{7 c_i t}{4 r_s} \right)^{4/7}. \quad (\text{C1}) $$

Alternatively, if $v_{\text{ext}}$ is non-negligible, $r_i$ tends towards a limit where $\dot{r} = 0$

$$ r_{\text{stall}} = r_s \left( \frac{c_i}{v_{\text{ext}}} \right)^{4/3}. \quad (\text{C2}) $$

Setting $r_i(t_{\text{stall}}) = r_{\text{stall}}$, and combining equations (C1) and (C2), we find

$$ \left( \frac{c_i}{v_{\text{ext}}} \right)^{4/3} = \left( \frac{7 c_i t_{\text{stall}}}{4 r_s} \right)^{4/7}. \quad (\text{C3}) $$

Invoking equation (C2) again to replace $r_s$ with $r_{\text{stall}}$, we find

$$ t_{\text{stall}} = \frac{4 r_{\text{stall}}}{7 v_{\text{ext}}}. \quad (\text{C4}) $$

We now calculate the value of $r_{\text{stall}}/v_{\text{ext}}$. We assume, as in Section 3, that $v_{\text{ext}}$ is the virial velocity at $r_{\text{stall}}$, i.e.

$$ v_{\text{ext}}^2 = \frac{6 GM}{5 r_{\text{stall}}}, \quad (\text{C5}) $$

where

$$ M = \frac{4}{3} \pi r_{\text{stall}}^3 \rho_0, \quad (\text{C6}) $$

where $\rho_0 = n_0 m_H / X$. The free-fall time in this cloud core is

$$ t_{\text{ff}} = \sqrt{\frac{3\pi}{32 G \rho_0}}, \quad (\text{C7}) $$

and hence we can write

$$ t_{\text{ff}} = \sqrt{\frac{6 \pi^2 r_{\text{stall}}}{5 8 v_{\text{ext}}} = 1.2 \frac{r_{\text{stall}}}{v_{\text{ext}}}}. \quad (\text{C8}) $$

Combining equations (C8) and (C4), we can write

$$ t_{\text{stall}} \approx 0.7 t_{\text{ff}}. \quad (\text{C9}) $$
Comparing this value as plotted on Fig. C1 (by eye), the ionization front will reach a value close to $r_{stall}$ over $\approx 2t_{stall}$. Hence, the time-scale over which the ionization front stalls, assuming it remains inside the core where $n_{ext} = n_0$, is roughly $1.4t_{ff}$, which is on the order of $t_{ff}$. Note that this is a very crude estimate, as it simplifies greatly the full equations that govern the ionization front. We discuss this result further in Section 3.3.

## APPENDIX D: COLLAPSING IONIZATION FRONTS

In our simulation N49-NSN, the ionization front collapses after expanding briefly. We introduce a simple spherically symmetric model invoking accretion on to the cloud core. Ntormousi & Hennebelle (2015), based on Larson (1969), give the time-dependent density of an accreting cloud as

$$n_{\text{ext}} = \frac{n_0}{\left(1 - \frac{r}{t_{\text{ff}}}\right)^2} \qquad (\text{D1})$$

Solving equation (B1) using this equation and assuming no external velocity field, we arrive at

$$r_i(t) = r_s \left\{ \frac{7 c_i}{4 r_s} t \left( 1 - \frac{t}{2t_{\text{ff}}} \right) \right\}^{4/7} \qquad (\text{D2})$$

Note that the density field has a singularity at $t_{ff}$. Beyond this point, the solution for $r_i$ becomes unphysical, and we instead keep $r_i = 0$.

We plot the solution to equation (D2) in Fig. D1 for the $10^{49}$ s$^{-1}$ source, whose ionization front collapses in our cloud. In addition, we solve equation (B1) numerically in the case where $v_{ext}$ is set to the virial velocity (as described in Section 3.3). As in Section 3.3, $t_{ff}$ is calculated for the cloud core rather than the cloud as a whole. We find that the numerical solution with $v_{ext} = v_{vir}$ agrees with the simulation results reasonably well, except for the flickering of the H II region due to the orbits of dense clumps passing through the source position.

<figure><img src="image_11.png" /><figcaption>**Figure D1.** Comparison between the radial expansion and collapse of the ionization front in the simulations and analytic models. The simulation is N49-NSN, i.e. using a $10^{49}$ s$^{-1}$ ionizing photon source in the cloud discussed in this paper. Analytic solutions are found to equation (B1) as described in Appendix D.</figcaption></figure>

Hence, in both the case where the front stalls and the case where it collapses, the time evolution of the ionization front is governed by the free-fall time in the cloud core.

## APPENDIX E: SHELL EXPANSION UNDER GRAVITY

In a spherically-symmetric solution, the supernova must entrain all the material in its path as it expands. However, in the 3D case, dense clumps of gas or (not included in these simulations) stars will remain embedded inside the supernova remnant. This provides an additional gravitational force on the shell. Similar models are derived for H II regions in Garcia-Segura & Franco (1996), Keto (2002) and Didelon et al. (2015).

In this model, we assume a shell moving outwards spherically, entraining all mass enclosed within it except for a fixed central mass. The mass of the shell is assumed to be the total mass $m(r)$ displaced by the shell at $r$. We assume a power-law density field with index $w$ and characteristic density $\rho_0$ and radius $r_0$, defined as $\rho(r) = \rho_0 r/r_0^{-w}$. Integrating, we find $m(r) = \frac{4}{3}\pi\rho_0 r_0^w r^{3-w}$. The reaction force from mass accretion on to the shell in a power-law density field can be written as

$$m\ddot{r} = -\dot{m}\dot{r} = -\frac{dm}{dr}\dot{r}^2 \qquad (\text{E1})$$

Dividing by $m$, we find

$$\ddot{r} = -\frac{(3-w)\dot{r}^2}{r} \qquad (\text{E2})$$

Note that this equation assumes the mass outside the supernova remnant is static. In our simulations this assumption is not too unreasonable, since the cloud is roughly virialized and the dissipation time-scale for the turbulence in the cloud that drives dynamic evolution in the cloud is longer than the time over which the supernova remnant evolves.

Including gravity, the equation becomes

$$\ddot{r} = -\frac{GM_c}{r^2} = -\frac{(3-w)\dot{r}^2}{r} \qquad (\text{E3})$$

where the central mass is $M_c$. Note that neither equation depends on the ambient density of the medium, although the density will set the initial velocity and radius of the shell as it enters the momentum-driven phase as given in equation (3).

Note that if $r$ shrinks, $M_c$ will also drop, whereas in reality, a contracting shell would retain its mass. Thus, the solution to equation (E3) after $r$ begins shrinking should be used with some caution. Rather, the value of this expression is determining at what point the shell stalls under gravity in the presence of a central gravitating mass.

In Fig. E1, we provide numerical solutions for this equation with and without the gravity for three values of $\rho_0$ with $w = 0$ and an initial shell velocity and radius calculated by assuming a momentum of $10^{43}$ g cm s$^{-1}$ and radius $r_\Lambda$ (equation 2). The mass of the embedded cluster $M_c$ is assumed to be a point mass of $10^4$ M$_\odot$, or 10 per cent of the total mass of the cloud. In our simulations, gravity should only become important above $10^3$ atoms cm$^{-3}$. This is a highly simplistic view of the properties of the central mass, so these results should be considered largely illustrative.
<figure><img src="image_12.png" /></figure>

**Figure E1.** Solutions to equation (E3) in Section E with and without gravity. We assume an initial ballistic momentum $10^{43}$ g cm s$^{-1}$, radius $r_c$ in equation (2) and $M_c = 10^4 M_{\odot}$.

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