ESCAPE FRACTION OF IONIZING PHOTONS DURING REIONIZATION: EFFECTS DUE TO SUPERNOVA FEEDBACK AND RUNAWAY OB STARS

Taysun Kimm; Renyue Cen

doi:10.1088/0004-637X/788/2/121

1. INTRODUCTION

Gunn & Peterson (1965) predicted that Lyα absorption would give rise to a sudden drop of continuum flux at wavelengths shorter than 1216 Å if a tiny amount of neutral hydrogen is present along the line of sight. The dramatic clearing of the Gunn–Peterson trough from the observation of quasars at z ∼ 6 demonstrates that hydrogen in the universe is highly ionized at z ≲ 6 (Becker et al. 2001; Fan et al. 2001, 2006). Polarization signals from the comic microwave background (CMB) also suggest that a large fraction of hydrogen may already be ionized by z ∼ 10–12 (Komatsu et al. 2011; Planck Collaboration et al. 2013). Yet, the detailed processes on how reionization has occurred remain unclear.

In the standard ΛCDM universe, dwarf galaxies form early (e.g., Somerville et al. 2003) and could dominate the budget of hydrogen ionizing photons at the epoch of reionization. Photons that escape from the porous interstellar medium (ISM, Clarke & Oey 2002), driven by supernova (SN) explosions (McKee & Ostriker 1977), to the intergalactic medium (IGM) create H ii bubbles, which expand as more stars form. The eventual percolation of H ii bubbles would mark the end of the cosmological reionization (e.g., Gnedin 2000a; McQuinn et al. 2007; Shin et al. 2008). This stellar reionization scenario has been studied extensively, both (semi-)analytically (e.g., Madau et al. 1999; Miralda-Escudé et al. 2000; Barkana & Loeb 2001; Bianchi et al. 2001; Cen 2003; Wyithe & Loeb 2003; Somerville et al. 2003; Bolton & Haehnelt 2007; Wyithe & Cen 2007; Kuhlen & Faucher-Giguère 2012; Robertson et al. 2013) and numerically (e.g., Gnedin 2000a; Razoumov et al. 2002; Ciardi et al. 2003; Fujita et al. 2003; Trac & Cen 2007; Gnedin et al. 2008; Wise & Cen 2009; Razoumov & Sommer-Larsen 2010; Yajima et al. 2011; Paardekooper et al. 2013). It appears that dwarf galaxies are the most plausible source of the ionizing photons, provided that the escape fraction is significant ( ${f_{\rm esc}}>10 \%$ ). Active galactic nuclei also contribute to ionizing photons in both the ultraviolet (UV) and X-ray bands but are generally believed to be subdominant to stellar sources (Haehnelt et al. 2001; Wyithe & Loeb 2003; Schirber & Bullock 2003; Faucher-Giguère et al. 2008; Cowie et al. 2009; Willott et al. 2010; Fontanot et al. 2014). The strong accretion shock present in massive halos ( ${{M}_{\rm vir}}\gtrsim 10^{10.5}\, {{M}_\odot }$ ) may also produce a nonnegligible amount of hydrogen ionizing photons in the vicinity of the galactic gaseous disk (Dopita et al. 2011).

The major uncertainty in the dwarf galaxy-driven reionization picture is the escape fraction of ionizing photons. Observationally, this is difficult to probe, because the hydrogen ionizing photons escaping from dwarf galaxies will get easily absorbed by the IGM during reionization (z ≳ 7). Besides, it requires a large sample of galaxies to obtain a statistically significant estimate of the escape fraction (f_esc). Nevertheless, it is worth noting that galaxies at higher redshift often exhibit a larger relative escape fraction ( $f_{\rm esc}^{\rm rel}$ ), which is defined as the ratio of the escape fraction at 900 Åand 1500 Å, than their low-z counterparts (Siana et al. 2010). Observations of star-forming galaxies at z ≲ 1 indicate that the relative escape fraction is only a few percent (Leitherer et al. 1995; Deharveng et al. 2001; Malkan et al. 2003; Siana et al. 2007, 2010; Cowie et al. 2009; Bridge et al. 2010). The only exception reported so far is Haro 11, which shows f_esc ∼ 4%–10% (Bergvall et al. 2006). On the other hand, a non-negligible fraction (∼10%) of star-forming galaxies at z ∼ 3 reveals a high escape of $f_{\rm esc}^{\rm rel} \ge 0.5$ (Shapley et al. 2006; Iwata et al. 2009; Nestor et al. 2011, 2013; Cooke et al. 2014). For typical Lyman break galaxies at z ∼ 3 in which 20%–25% of UV photons are escaping (Reddy et al. 2008), the relative fraction corresponds to a high escape fraction of ${f_{\rm esc}}\sim 0.1$ . Given that galaxies are more actively star forming at high redshift (e.g., Bouwens et al. 2012; Dunlop et al. 2013), it has been suggested that there may be a correlation between star formation rate and f_esc, and possibly evolving f_esc with redshift (Kuhlen & Faucher-Giguère 2012).

Predicting the escape fraction in theory is also a very challenging task. This is essentially because there is little understanding on the structure of the ISM at high-z dwarf galaxies. Numerical simulations are perhaps the most suited to investigate this subject, but different subgrid prescriptions and/or finite resolution often lead to different conclusions. Using an adaptive mesh refinement (AMR) code, ART (Kravtsov et al. 1997), with SN-driven energy feedback, Gnedin et al. (2008) claim that the angle-averaged escape fraction increases with galaxy mass from 10⁻⁵ to a few percents in the range 10¹⁰ ≲ M_gal ⩽ 4 × 10¹¹. They attributed this trend to the fact that more massive galaxies have smaller gas-to-stellar scale-height than lower mass galaxies in their simulations. On the other hand, Razoumov & Sommer-Larsen (2010) argue based on cosmological TreeSPH simulations (Sommer-Larsen et al. 2003) that more than 60% of the hydrogen ionizing photons escape from dwarf galaxies in dark matter halos of M_halo = 10⁸–10⁹ M_☉. More massive halos of 10¹¹ M_☉ are predicted to have a considerably smaller f_esc (≲ 10%). A similar conclusion is reached by Yajima et al. (2011). It should be noted, however, that resolution could potentially be an issue in these two studies in the sense that their resolution of a few hundreds to thousands of parsec is unable to resolve most star-forming regions and hence capture obscuring column densities and a porous ISM. Wise & Cen (2009) performed cosmological radiation hydrodynamic simulations employing very high resolution (0.1 pc), and found that the neutral hydrogen column density varies over the solid angles from $N_{{\rm H\,\scriptsize{I}}}\sim 10^{16}\, {\rm cm^{-2}}$ to 10²² cm⁻² with the aid of SN explosions and photo-ionization. Because of the porous ISM, a high f_esc of ∼40% is achieved in small halos of M_halo = 10⁷–10^9.5 M_☉. Wise et al. (2014) show that an even higher fraction (∼50%) of hydrogen ionizing photons escapes from minihalos of M_halo = 10^6.25–10⁷ M_☉.

Another potentially important source of ionizing radiation is runaway OB stars that are dynamically displaced from their birthplace. The runaway OB stars are normally defined by their peculiar motion (v_pec ⩾ 30 km s⁻¹; Blaauw 1961), and roughly 30% of OB stars are classified as runaways in the Milky Way (Stone 1991; Hoogerwerf et al. 2001; Tetzlaff et al. 2011). Although the fraction is still uncertain, their peculiar speed of 〈v_pec〉 ∼ 40 km s⁻¹ means that the runaway OB stars can, in principle, travel away from the birthplace by ∼200 pc in 5 Myr, making them an attractive source for the ionizing photons. The runaway OB stars are thought to originate from a three-body interaction with other stars in a young cluster (Leonard & Duncan 1988), and/or from a SN explosion of a companion in a binary system (Blaauw 1961). Conroy & Kratter (2012) evaluated the impact of the inclusion of runaway OB stars on f_esc using a simple analytic argument, and concluded that the runaway OB stars may enhance f_esc by a factor of up to ∼4.5 in halos with M_halo = 10⁸–10⁹ M_☉.

The aim of this study is to investigate the importance of the aforementioned two processes by measuring the escape fraction from high-resolution cosmological radiation hydrodynamics simulations. First, given that modeling the SN explosion as thermal energy is well known to have the artificial radiative cooling problem (e.g., Katz 1992; Slyz et al. 2005), we expect that the role of the SN is likely to be underestimated in some cosmological simulations (e.g., Gnedin et al. 2008). With a new physically based SN feedback model that captures all stages of the Sedov explosion from the free expansion to the snowplow phase, we study the connection between the escape of ionizing photons and feedback processes in dwarf galaxies. Second, we extend the idea by Conroy & Kratter (2012) and quantify the impact from the runaway OB stars on reionization in a more realistic environment.

We first describe the details of our cosmological radiation hydrodynamics simulations including the implementation of runaway OB stars in Section 2. We present the feedback-regulated evolution of the escape fraction and the impact of the inclusion of runaway OB stars in Section 3. We summarize and discuss our findings in Section 4. Our new mechanical feedback from SN explosions is detailed in Appendices A–C.

2. METHOD

2.1. Hydrodynamics Code

We make use of the Eulerian AMR code, ramses (Teyssier 2002; ver. 3.07), to investigate the escape of ionizing radiation from high-z galaxies. ramses is based on the fully threaded oct-tree structure (Khokhlov 1998), and uses the second-order Godunov scheme to solve Euler equations. The hydrodynamic states reconstructed at the cell interface are limited using the MinMod method, and then advanced using the Harten–Lax–van-Leer contact wave Riemann solver (HLLC; Toro et al. 1994). We adopt a typical Courant number of 0.8. The Poisson equation is solved using the adaptive particle-mesh method. Gas can effectively cool down to 10⁴ K by atomic and metal cooling (Sutherland & Dopita 1993). Below 10⁴ K, metal fine-structure transitions, such as [C ii] 158 μm, can further lower the temperature down to 10 K, as in Rosen & Bregman (1995). We set the initial metallicity to 2 × 10⁻⁵, as primordial SNe can quickly enrich metals in mini-halos of mass 10⁷ M_☉ (e.g., Whalen et al. 2008), which our simulations cannot resolve properly.

We use the multi-group radiative transfer (RT) module developed by Rosdahl et al. (2013) to compute the photoionization by stars. The module solves the moment equations for three photon packets (H ii, He ii, and He iii ionizing photons) using a first-order Godunov method with M1 closure for the Eddington tensor. We adopt the Harten–Lax–van-Leer (HLL; Harten et al. 1983) intercell flux function. Ionizing photons from each star are taken into consideration in every fine step. Note that an advantage of the moment-based RT is that it is not limited by the number of sources. The production rate of the ionizing photon varies with time for a given initial mass function (IMF; Leitherer et al. 1999, see also Rosdahl et al. 2013). The majority of the ionizing photons are released in ∼5 Myr of stellar age. We adopt the production rate equivalent to that of Kroupa IMF (Kroupa 2001) from the Starburst99 library (Leitherer et al. 1999).¹ The radiation is coupled with gas via photo-ionization and photo-heating, and a set of nonequilibrium chemistry equations for H ii, He ii, and He iii are solved similarly as in Anninos et al. (1997). We assume that photons emitted by recombination are immediately absorbed by nearby atoms (case B). The speed of light is reduced for the speed-up of the simulations by 0.01 (e.g., Gnedin & Abel 2001). This is justifiable because we are mainly interested in the flux of escaping photons at the virial sphere.

2.2. Cosmological Simulations

We carry out cosmological simulations to investigate the escape fraction in realistic environments. For this purpose, we generate the initial condition using the music software (Hahn & Abel 2011), with the WMAP7 cosmological parameters (Komatsu et al. 2011): (Ω_m, Ω_Λ, Ω_b, h, σ₈, n_s = 0.272, 0.728, 0.045, 0.702, 0.82, 0.96). A large volume of (25 Mpc h⁻¹)³ is employed to include the effect of the large-scale tidal field. To achieve high mass resolution, we first run dark matter-only simulations with 256³ particles, and identify a rectangular region of 3.8 × 4.8 × 9.6 Mpc (comoving) that encloses two dark matter halos of ≃ 1.5 × 10¹¹ M_☉ at z = 3. Then, we further refine the mass distribution of the zoomed-in region, such that the mass of a dark matter particle is m_dm = 1.6 × 10⁵ M_☉, which corresponds to 2048³ particles in effect. Despite that we purposely select the region in which two massive dark matter halos are present at z = 3, a comparison with the number of dark matter halos per volume predicted by Jenkins et al. (2001) shows that our simulated box represents an average region of the universe at z = 7 (Figure 1).

**Figure 1.** Dark matter halo mass function from the zoomed-in region of the ${\mathsf {FR}}$ run at z = 7. Comparison with Jenkins et al. (2001) mass function at the same epoch indicates that our simulated volume represents the average region of the universe.
Download figure:
Standard image High-resolution image

The level of the root grid in the zoomed-in region is 11, consistent with the dark matter resolution. Further 12 levels of refinement are triggered if the dark matter plus baryon mass in a cell exceeds 8 times the mass of a dark matter particle. We keep the minimum physical size of a cell to Δx_min = 25 Mpc h⁻¹/2²³ = 4.2 pc over the entire redshift. However, this refinement criterion is not optimized to resolve the structure of the ISM, unless extremely high mass resolution is adopted. For example, for a gas cell of n_H = 10 cm⁻³, the criterion will come into play only if the size of the cell is larger than ∼160 pc. In order to better resolve the structure of the ISM, we enforce a cell with n_H ⩾ 1 cm⁻³ to be resolved on 8Δx_min = 34 pc. In a similar context, we apply more aggressive refinement criterion for the star-forming gas in such a way that gas with n_H = 100 cm⁻³ (800 cm⁻³) is always resolved on a 8.5 pc (4.2 pc) cell. We adopt very high stellar mass resolution of ≈49 M_☉. This means that a star particle with the minimum mass will produce a single SN event for the Chabrier IMF.

We run two sets of cosmological simulations, ${\mathsf {FR}}$ and ${\mathsf {FRU}}$ , with the identical initial condition down to z = 7 (Table 1). Both runs include star formation, metallicity-dependent radiative cooling (Sutherland & Dopita 1993; Rosen & Bregman 1995), thermal stellar winds, mechanical feedback from SN explosions, and photoionization by stellar radiation. The runaway OB stars are included only in the ${\mathsf {FRU}}$ run. In Figure 2, we show an example of the growth of H ii bubbles in the ${\mathsf {FR}}$ run. Our simulated region is nearly ionized at z = 7.

**Figure 2.** Expansion of the H ii bubble in a cosmological simulation ( ${\mathsf {FR}}$ ). The three panels show the evolution of the density-weighted fraction of ionized hydrogen of the zoomed-in region. The horizontal size of the figure is 9.5 Mpc (comoving).
Download figure:
Standard image High-resolution image

Table 1. Summary of Cosmological Simulations

Model	SNII	RT	Runaways	Δx_min	m_{star, min}	m_dm
Model	SNII	RT	Runaways	(pc)	(M_☉)	(10⁵ M_☉)
FR	$\checkmark$	$\checkmark$	⋅⋅⋅	4.2	49	1.6
FRU	$\checkmark$	$\checkmark$	$\checkmark$	4.2	49	1.6

Download table as: ASCII Typeset image

Dark matter (sub) halos are identified using the Amiga halo finder (Ahf; Gill et al. 2004; Knollmann & Knebe 2009). Ahf first constructs the adaptive meshes based on the particle distribution, finds the density minima, and determines physical quantities based on a virial overdensity (Δ_vir). Gravitationally unbound particles are removed iteratively if they move faster than the local escape velocity during this procedure. The virial radius is defined such that the mass enclosed within the virial sphere is the virial overdensity times the critical density of the universe times the volume, i.e., ${{M}_{\rm vir}}(z) = \Delta _{\rm vir}(z) \rho _{\rm crit}(z) 4 \pi r_{\rm vir}^3 / 3$ . We take Δ_vir = 177, appropriate for a Λ-dominated universe at z > 6 (Bryan & Norman 1998). This results in 796, 443, and 183 dark matter halos of mass ${{M}_{\rm vir}}\ge 10^{8}\,{{M}_\odot }$ immune to the contamination by coarse dark matter particles (m_dm > 1.6 × 10⁵ M_☉) at z = 7, 9, and 11, respectively.

2.3. Star Formation and Feedback

Stars form in a very dense, compact molecular core. Infrared extinction maps of nearby interstellar cores indicate that their size ranges from 0.01 to 0.4 pc (e.g., Alves et al. 2007; Könyves et al. 2010), which is difficult to resolve in current cosmological simulations. Nevertheless, studies of gravitational collapse in converging flows (Gong & Ostriker 2011) seem to suggest that a gravitationally bound cloud is likely to experience runaway collapse no matter how the collapse is initiated. In a similar spirit, we assume that stars would form in a cell if the following conditions are met simultaneously (e.g., Cen & Ostriker 1992).

1.
The flow is convergent ( $\boldsymbol {\nabla }\cdot \rho {\boldsymbol {v}} <0$ ).
2.
The cooling time is shorter than the dynamical time.
3.
The gas is Jeans unstable.
4.
The number density of hydrogen exceeds the threshold density n_th = 100 cm⁻³.

The last condition is motivated by the density of a Larson–Penston profile (Larson 1969; Penston 1969) at 0.5Δx, $\rho _{\rm LP}\approx 8.86 c_s^2 / \pi \,G\,\Delta x^2$ , where c_s is the sound speed and Δx is the size of the most refined cell. Star particles are created based on the Schmidt law (Schmidt 1959), $\dot{\rho }_{\star } = \epsilon _{\rm ff} \, \rho _{\rm gas} / t_{\rm ff}$ , assuming that 2% of the star-forming gas ( epsilon _ff) is converted into stars per its free-fall time (t_ff) (Krumholz & Tan 2007; Kennicutt 1998). The mass of each star particle is determined as $m_\star =\alpha \, N_p \rho _{\rm th} \, \Delta x_{\rm min}^3$ , where ρ_th is the threshold density for star formation, Δx_min is the size of the most refined cell, and α is a parameter that controls the minimum mass of a star particle. N_p is the number of star particles to be formed in a cell, which is drawn from a Poisson random distribution, $P(N_p) = (\lambda ^{N_p} / N_p!) \exp \left(-\lambda \right)$ . Here the Poissonian mean (λ) is computed as λ ≡ epsilon _ff(ρΔx³/m_{⋆, min})(Δt_sim/t_ff), where Δt_sim is the simulation time step, and m_{⋆, min} is the minimum stellar mass (i.e., N_p = 1).

We describe the SN feedback using a new physical model which captures the SN explosion at all stages from the early free expansion to the final momentum-conserving snowplow phase. Briefly, we deposit radial momentum to the cells affected by SN feedback, conserving energy appropriately. The amount of input momentum is determined by the stage the blast wave is in, which in turn is dependent upon the physical condition (density and metallicity) of the gas being swept up and simulation resolution. The virtue of our scheme is that an approximately (within 20%) correct amount of momentum is imparted to the surrounding gas regardless of the resolution. Thus, this prescription should be useful to cosmological simulations, especially those with finite resolution that potentially suffer from the artificial radiative cooling. The details of our implementation and a simple test are included in Appendices A–C.

The frequency of a SN per solar mass is estimated assuming the Chabrier IMF (Chabrier 2003). For the simple stellar population with a low-(high-)mass cutoff of 0.1 (100) M_☉, the total mass fraction between 8 and 100 M_☉ is 0.317, and the mean SN progenitor mass is 15.2 M_☉ on the zero-age main sequence. At the time of the explosion, we also deposit newly processed metals into the surrounding. The mass fraction of newly synthesized metals in stellar ejecta is taken to be 0.05, following Arnett (1996). A star particle is assumed to undergo the SN phase after the main sequence lifetime of the mean SN progenitor (10 Myr; Schaller et al. 1992). As discussed in Slyz et al. (2005), allowing for the delay between the star formation and explosion (i.e., stellar lifetimes) is crucial to the formation of hot bubble in the ISM. We find that the physically based SN feedback employed in this study drives stronger galactic winds than the runs with thermal feedback or kinetic feedback that are valid only under certain conditions (Dubois & Teyssier 2008, see below). Stellar winds from massive stars are modeled as thermal input, based on Leitherer et al. (1999).

2.4. Runaway OB Stars

Our implementation of runaway OB stars is largely motivated by Tetzlaff et al. (2011), who compiled candidates of runaway stars younger than 50 Myr for the 7663 Hipparcos sample. By correcting the solar motion and Galactic rotation, they found that the peculiar space velocity of the stars may be decomposed into two Maxwellian distributions intersecting at 28 km s⁻¹. Assuming that each Maxwellian distribution represents a kinematically distinctive population, they estimated the fraction of the runaways to be ∼27.7% ± 1.9 for the sample with full kinematic information. The dispersion of the Maxwellian distribution is measured as 24.4 km s⁻¹ for the high-velocity group.

Since either runaway OB stars formed through the explosion of a SN in a binary or those dynamically ejected in a cluster are not resolved in our simulations, we crudely approximate this by splitting a star particle into a normal (70% in mass) and a runaway particle (30%) at the time of star formation. While the initial velocity of the normal star is chosen as the velocity of the birth cloud, we add a velocity drawn from the Maxwellian distribution on top of the motion of the birth cloud for runaway particles. To do so, we generate the distribution following the Maxwellian with the dispersion of σ_v = 24.4 km s⁻¹ and the minimum space velocity of v_3D = 28 km s⁻¹ using the rejection method (Press et al. 1992). The direction of the runaway motion is chosen randomly for simplicity. A similar approach is taken by Ceverino & Klypin (2009) to study the formation of disk galaxies in a cosmological context.

2.5. Estimation of Escape Fraction

The fraction of escaping ionizing photons (f_esc) is measured by comparing the photon flux at the virial radius and the photon production rate from young massive stars. Since the speed of light is finite, there is a small delay in time between the photons produced by the stars and the photons escaping at the virial sphere. In order to take this into account, we use the photon production rate at earlier time (t − r_vir/c'), where c' is the reduced speed of light used in the simulations. The escape fraction is then computed as

$\begin{equation} f_{\rm esc}(t) \equiv \frac{\int d\Omega \, \boldsymbol {F}_{\rm ion}(t) \cdot \hat{r} \Theta (\boldsymbol {F}_{\rm ion}\cdot \hat{r})}{\int dm_* \, \dot{N}_{\rm ion} (t-r_{\rm vir}/c^{\prime })}, \end{equation} \tag{ 1 }$

where $\boldsymbol {F}_{\rm ion}$ is the ionizing photon flux, dΩ is the solid angle, m_* is the mass of each star particle, $\dot{N}_{\rm ion}(t)$ is the photon production rate of a simple stellar population of age t per solar mass, and Θ is the Heaviside step function. Here, we approximate the delay time to be a constant, r_vir/c', for each halo assuming that the central source is point-like. Since only outflowing photons are considered in Equation (1), we find that a minor fraction (∼5%) of galaxies exhibit f_esc greater than 1. This happens mostly when there is little absorbers left in the halo after disruptive SN explosions. In this case, we randomly assign f_esc between 0.9 and 1.0. We confirm that the photon production rate-averaged escape fraction, which is the most important quantity in this study, is little affected by this choice even if the net flux is used, and thus we decide to take a simpler method.

Dust can also affect the determination of the escape of the hydrogen ionizing photons. However, given that our simulated galaxies are very metal poor (0.002–0.05 Z_☉) and galaxies with lower metallicity have a progressively lower amount of dust (Lisenfeld & Ferrara 1998; Engelbracht et al. 2008; Galametz et al. 2011; Fisher et al. 2014), it is unlikely that dust decreases the escape fraction substantially. Thus, we neglect the absorption of hydrogen ionizing photons by dust in this study.

3. RESULTS

3.1. Feedback-regulated Escape of Ionizing Photons

Cosmological hydrodynamics simulations often suffer from the artificial over-cooling problem in forming disk galaxies (e.g., Kimm et al. 2011a; Hummels & Bryan 2012), mainly because the energy from SN explosions is radiated away before it is properly transferred to momentum due to inadequate resolution of the multi-phase ISM. This directly affects the escape of ionizing photons. Motivated by this challenge, we have implemented a SN feedback scheme that reasonably approximates the Sedov blast waves from the free expansion to snowplow stages. In Figure 3, we present the baryon-to-star conversion efficiency (f_⋆ ≡ M_star/(Ω_bM_vir/Ω_m) of the central galaxies in dark matter halos at z = 7 from the ${\mathsf {FR}}$ run. It shows that our new physically motivated SN feedback is very effective at suppressing star formation. For example, the most massive halo with M_vir ∼ 3 × 10¹⁰ M_☉ at z = 7 shows f_⋆ ≈ 0.08. Although the direct comparison may be difficult due to a different initial condition used, it is worth noting that the conversion efficiency is about a factor of 7 smaller than that found in the NutFB run (Kimm et al. 2011a; see Figure 13), shown as a star in Figure 3. We note that the momentum input from SN explosions used in the NutFB run is a factor of 3–4 smaller compared with that at the end of the cooling phase (see Appendices A–C; Blondin et al. 1998). For lower mass halos, the conversion efficiency is found to be even lower, reaching M_star/M_vir ≲ 0.01 Ω_b/Ω_m at M_vir ∼ 10⁹ M_☉. It is also interesting to note that the conversion efficiency at M_vir ⩾ 10¹⁰ M_☉ also agrees reasonably well within error bars with the semi-analytic results obtained to reproduce the observed stellar mass function, star formation rate, and cosmic star formation rate density (e.g., Behroozi et al. 2013, Figure 7). As the feedback becomes more effective and fewer stars are formed, the stellar metallicity of these high-z galaxies would be lower. We find that the most massive galaxy in our z = 7 sample (M_star = 4 × 10⁸ M_☉) has a stellar metallicity of 0.05 Z_☉. This is at least factor of 2–3 smaller than the prediction by Finlator et al. (2011) at the same epoch. Kimm & Cen (2013) also investigated UV properties of z = 7 galaxies of stellar mass 5 × 10⁸–3 × 10¹⁰ M_☉ using a SN energy-based feedback scheme, and found that stellar metallicities are generally higher than those found in the ${\mathsf {FR}}$ run. Kimm & Cen (2013) found that the stellar metallicity for galaxies of mass 4 × 10⁸ M_☉ falls in the range of 0.1–0.5 Z_☉. The gas metallicities (Z_gas) are also different in the two simulations. The gas metallicity of the ISM within 2.56 kpc for the 4 × 10⁸ M_☉ galaxies is 0.083 Z_☉ in the FR run, which is about a factor of 3three lower, on average, than that of Kimm & Cen (2013) (Z_gas = 0.1–0.7 Z_☉). These comparisons lead us to conclude that our physically based feedback scheme is effective in alleviating the overcooling problem.

**Figure 3.** Baryon-to-star conversion efficiency at z = 7 from the ${\mathsf {FR}}$ (blue) and the ${\mathsf {FRU}}$ (orange) runs. Only central galaxies are shown. The cosmic mean (Ω_b/Ω_m = 0.165) is shown as a black solid line. Also included as a star is the stellar fraction measured from the NutFB simulation (Kimm et al. 2011a). Our mechanical feedback from SN explosions is more effective at regulating star formation, compared with previous studies injecting thermal or kinetic energy (see the text).
Download figure:
Standard image High-resolution image

One may wonder whether stars form inefficiently in these small halos (10⁸ ≲ M_vir ≲ 10⁹ M_☉) because gas accretion is suppressed due to the ionizing background radiation (Shapiro et al. 1994; Thoul & Weinberg 1996; Gnedin 2000b; Dijkstra et al. 2004; Sobacchi & Mesinger 2013; Noh & McQuinn 2014). However, this is unlikely the case, given that galaxies in the atomic cooling halos are fed mainly by dense filaments and satellites at high redshift (e.g., Powell et al. 2011), which are self-shielded from the background radiation (Faucher-Giguère et al. 2010; Rosdahl & Blaizot 2012). Even in the absence of the self-shielding, Geen et al. (2013) find no clear sign that reionization suppresses star formation in such halos at z > 6. Wise et al. (2014) also show that the fraction of baryons in a 10⁸–10⁹ M_☉ halo is reduced only by less than a factor of two compared with the cosmic mean in their cosmological radiation hydrodynamics simulations with thermal SN feedback and reionization. Indeed, we confirm that our mechanical SN feedback is primarily responsible for the low conversion efficiency by directly comparing the stellar mass of the dwarf galaxies between the simulations with and without ionizing radiation (see Appendices A–C).

We now present the time evolution of star formation rate and ionizing photon escape fraction of two randomly chosen relatively massive galaxies in Figure 4. The plot corroborates that the feedback from stars governs the evolution of galaxies. The top and bottom panels show the evolution of specific star formation rate (sSFR $\equiv \dot{M}_{\rm star}/M_{\rm star}$ ) and instantaneous f_esc of the central galaxy in dark matter halos of mass 3 × 10¹⁰ and 10¹⁰ M_☉, respectively. The SFR is computed by averaging the mass of newly formed stars over 3 Myr. It is evident that star formation is episodic on a timescale of 10–30 Myr with both the frequency and oscillation amplitude decreasing with increasing stellar mass. This means that SN explosions effectively control the growth and disruption of star-forming clouds. When the galaxies are small (t_H ≲ 0.5 Gyr), the explosions even completely shut down the star formation across the galaxies, as stars form only in a few dense clouds. During these quiet periods, f_esc is kept high ( ${f_{\rm esc}}\gtrsim 0.2$ ). On the other hand, massive galaxies contain many star-forming clumps, as can be seen in the projected density plot (middle row). The fact that the episodic star formation history becomes more smooth at late times indicates that these clumps are not entirely susceptible, but somewhat resilient to the SN explosions arising from neighboring star clusters.

**Figure 4.** Evolution of the escape fraction (f_esc) and specific star formation rate (sSFR) in two massive halos from the ${\mathsf {FR}}$ run. Black solid lines in the top and bottom panels indicate the escape fraction measured at the virial radius at each snapshot as a function of the age of the universe. We denote the logarithmic stellar mass at different times by orange text. Black dashed lines correspond to the photon number-weighted average of f_esc by that time (<f_esc >). Blue shaded regions display the sSFR in Gyr⁻¹. One can see that there is a delay between the peak in f_esc and sSFR due to the delay in the onset of the strong outflow. The middle panels show an example of this delay identified in the top panel (a,b). The projected density of gas and the fraction of ionized hydrogen are shown in both cases, as indicated in each panel. Interestingly, the volume filling fraction of the neutral hydrogen within 0.2 R_vir is found to be 25% large in the snapshot (b), indicating that f_esc depends not only by the volume-filling, circumgalactic neutral gas, but also dense star-forming gas. We do not display the physical quantities if M_vir ⩽ 10⁸ M_☉.
Download figure:
Standard image High-resolution image

More importantly, we find that there is a time delay between the peak of f_esc and sSFR. This is essentially because massive stars with M ≈ 15 M_☉ explode ∼10 Myr after their birth in our simulation. Let us suppose a dense cloud that just begins to form stars. Since the gas flow is usually convergent in these regions, the density of the gas will rise with time, and so does the SFR. This means that more and more massive stars will explode as time goes on. Once enough SNe that can significantly redistribute the birth cloud go off, SFR will begin to drop, and f_esc will increase. Note that the increase in the number of SNe continues even after the peak of SFR, as massive stars live ∼10 Myr. Once the massive stars formed at the peak of SFR evolve off, star formation will be further suppressed as a result of the destruction of the star-forming clouds, and strong outflows are likely to be produced, thus maximizing f_esc. Therefore, the time delay stems from the interplay between the build-up of a noncoeval star cluster and subsequent SN explosions after the lifetime of the massive stars (∼10 Myr). The projected density distributions of gas at two snapshots, one of which displays the peak in sSFR (a) and the other shows the peak in f_esc (b), substantiates that it is indeed the strong outflow that elevates f_esc (middle row). When sSFR is at the peak value, the central galaxy appears relatively quiet (panel a), whereas strong outflows are seen when f_esc is highest and sSFR drops rapidly (panel b). As one can read from the figure, this mis-match of SFR and f_esc means that a large amount of ionizing photons at the peak of SF are absorbed by their birth clouds. Although f_esc is high in the early time (t_H ≲ 0.5 Gyr), the photon number-weighted mean f_esc (dashed lines) stays at around 10% level in these two examples.

We present statistical results of the escape fraction in Figure 5. Since there are a limited number of galaxies in our simulated volume and f_esc varies significantly on ∼10 Myr, we compute the median and interquartile range of f_esc by combining the results from seven consecutive snapshots spanning 21 Myr. Several features can be gleaned from this figure. First, although there is a considerable scatter, high-z galaxies exhibit a high f_esc on the order of 10%, which is normally required by semi-analytic calculations of reionization to ionize the universe by z ∼ 6 (Wyithe & Cen 2007; Shull et al. 2012; Robertson et al. 2013). Second, there is a hint that photons can escape more easily in the galaxies hosted by lower mass halos. We attribute this to the fact that feedback from stars efficiently destroys a few star-forming clouds that are responsible for the total SF in smaller halos, as opposed to larger ones in which young massive stars are buried in many star-forming clouds that are relatively resilient to the SN feedback arising from neighboring star clusters. As shown in the top and bottom panels of Figure 4, when galaxies are small, the entire star formation can be suppressed due to the energetic outflows driven by SN explosions. Third, we find that f_esc is slightly higher at lower redshift for a given halo mass, consistent with Paardekooper et al. (2013). This is essentially because the mean density of the gas is smaller at lower redshift, and the impact from SNe becomes more effective.

Note that high f_esc does not necessarily mean that more photons would leave their host halo. Star clusters older than ∼5 Myr would not contribute significantly to the total ionizing photon budget even if their f_esc is 1. The more relevant quantity for reionization should take into account the photon production rate, and we find that the (weak) redshift dependence of f_esc disappears when the photon escape rate is plotted (right panel in Figure 5). Since the instantaneous measurement of f_esc could be misleading, we also present the photon production rate-weighted, time-averaged escape fraction, ${\left<f_{\rm esc}\right>}({\le} t_{\rm H}) \equiv \int _0^{t_{\rm H}} \dot{N}_{\rm ion}(t) f_{\rm esc}(t) dt / \int _0^{t_{\rm H}} \dot{N}_{\rm ion}(t) dt,$ in Figure 6 (left panel). This is a better quantity to be used for the semi-analytic calculations of reionization than f_esc from Figure 5. Overall, we find that the time-averaged escape fraction at z = 7 is around ∼10%, regardless of the halo mass in the range considered. Also included as the black dotted line in Figure 6 is the photon production rate-weighted average of f_esc of all the samples at different times ( ${\left<f_{\rm esc}\right>}(t)$ ). Again, the value is found to fluctuate around 10%, but no clear sign of redshift dependence is detected.

**Figure 6.** Left: photon production rate-weighted escape fraction, ${\left<f_{\rm esc}\right>}$ , averaged over the age of the universe (t_H) in the ${\mathsf {FR}}$ run. The effective escape fraction in different halo mass bins is shown as different color codings, as indicated in the legend. We also display the photon rate-averaged escape fraction of the whole sample at each snapshot ( ${\left<f_{\rm esc}\right>}(t)$ ) (black dotted line), as opposed to the time-averaged quantities (solid and dashed lines). We find the effective escape fraction to be ∼10%, regardless of the halo mass and redshift. Altogether, 11.4% of the photons produced until z = 7 have escaped from halos of ${{M}_{\rm vir}}\ge 10^8\,{{M}_\odot }$ . Right: relative contribution of halos of different mass ranges to the total number of ionizing photons measured at the virial radius. The contribution is computed by taking into account the cumulative number of photons produced and the cumulative number of photons escaped from halos of relevant mass range until t ⩽ t_H.
Download figure:
Standard image High-resolution image

The relative contributions from halos of different masses to the total escaping ionizing photons are compared in Figure 6 (right panel). As the small structures form first in the ΛCDM universe, the small halos of mass ${{M}_{\rm vir}}\le 10^{8.5}\,{{M}_\odot }$ dominate down to z ∼ 9. More massive halos and galaxies emerge later, and their cumulative contribution becomes comparable with that of the smallest halos ( ${{M}_{\rm vir}}\le 10^{8.5}\,{{M}_\odot }$ ) by z = 7. In our simulations, 14 of the most massive halos supply more ionizing photons than the 556 smallest halos with ${{M}_{\rm vir}}\le 10^{8.5}\,{{M}_\odot }$ at z = 7. This is mainly because f_⋆ is much higher in the more massive halos than in the small halos, while the effective escape fraction is similar. The typical number of escaping photons per second in halos with ${{M}_{\rm vir}}\sim 10^{8.5}\,{{M}_\odot }$ is $f_{\rm esc}\,\dot{N}_{\rm ion}\sim 10^{49}\,{\rm s^{-1}}$ , whereas the number can increase up to $f_{\rm esc}\,\dot{N}_{\rm ion}\sim 10^{52}\,{\rm s^{-1}}$ in the most massive halos ( ${{M}_{\rm vir}}> 10^{10}\,{{M}_\odot }$ ; Figure 5, right panel). Notice, however, that this does not necessarily translate to their relative role to the reionization of the universe. Small halos at high redshift may make a more significant contribution to the Thompson optical depth (Wyithe & Cen 2007; Shull et al. 2012; Kuhlen & Faucher-Giguère 2012; Robertson et al. 2013).

It is noted that the recombination timescale corresponding to the mean density of the universe at z ∼ 10 (n_H ∼ 10⁻³ cm⁻³) is relatively long (∼50–100 Myr),² and thus the halo gas around a galaxy may be kept partially ionized even though it is irradiated by the galaxy intermittently. Figure 4 (the second panel in the middle row) indeed shows that a large fraction of the IGM in the vicinity of the central galaxy is largely ionized despite the fact that instantaneous f_esc is low. Although we do not include the whole distribution of the ionized hydrogen inside the halo, we confirm that the halo gas between 2 kpc and 12 kpc (virial radius) is fully ionized apart from the small region taken by cold filamentary gas. In fact, the volume filling fraction of the neutral hydrogen (f_v) inside $0.2\,{{R}_{\rm vir}}$ (∼2.3 kpc) is found to be ∼25% larger in the snapshot (b) (f_v ≈ 0.04) than that in the snapshot (a), suggesting that dense star-forming gas plays a more important role in determining the escape fraction than volume-filling diffuse neutral gas.

Figure 7 demonstrates the importance of resolving the ISM in predicting the escape of ionizing photons. In order to estimate the optical depth by neutral hydrogen in the vicinity of each star particle (<100 pc), we spawn 768 rays per particle using the Healpix algorithm (Górski et al. 2005). Each ray carries the spectral energy distribution determined by the age and mass of the star particle (Leitherer et al. 1999). As the ray propagates, we compute the absorption of the Lyman continuum by neutral hydrogen as $F_{\rm abs} (\nu) = F_{\rm int} (\nu) \exp {\left[-\tau _{{\rm H\,\scriptsize{I}}} (\nu)\right]}$ , where $\tau _{{\rm H\,\scriptsize{I}}}$ ( $=N_{{\rm H\,\scriptsize{I}}} \sigma _{{\rm H\,\scriptsize{I}}}$ ) is the optical depth and $\sigma _{{\rm H\,\scriptsize{I}}}$ is the hydrogen ionization cross section (Osterbrock & Ferland 2006). We then combine the attenuated spectral energy distributions propagated out to 100 pc from each star particle, and measure the remaining number of ionizing photons ( $N_{\rm ion,tot}^{\rm final}$ ) per galaxy. This is compared with the initial number of ionizing photons ( $N_{\rm ion,tot}^{\rm int}$ ) to obtain the effective optical depth as $\tau _{\rm eff, 100\,pc} \equiv \ln \left(N_{\rm ion,tot}^{\rm int} / N_{\rm ion,tot}^{\rm final} \right)$ . Figure 7 shows the distribution of the effective optical depth by the nearby gas for the galaxies with a low escape fraction ( ${f_{\rm esc}}< 0.1$ ) at z ∼ 8. We find that τ_{eff, 100 pc} shows a wide distribution ranging from 0.01 to ∼100, with the photon production rate-weighted averages of τ_{eff, 100 pc} = 3.8 and 1.9 for less ( $10^8 < {{M}_{\rm vir}}\le 10^9\,{{M}_\odot }$ ) and more massive ( $10^9 < {{M}_{\rm vir}}\le 10^{10.5}\,{{M}_\odot }$ ) halo groups, respectively. This indicates that the number of escaping photons is reduced by a factor of 7–45 due to the gas near young stars in galaxies with the small f_esc. In this regard, one may find it reconcilable that results from cosmological simulations with limited resolutions (e.g., Fujita et al. 2003; Razoumov & Sommer-Larsen 2010; Yajima et al. 2011) often give discrepant results.

To summarize, we find that there is a time delay between the peak of star formation activity and the escape fraction due to the delay in the onset of effective feedback processes that can blow birth clouds away. Because of the delay, only 11.4% of the ionizing photons could escape from their host halos when photon production rate-averaged over all halos at different redshifts, despite the fact that the instantaneous f_esc could reach a very high value temporarily. Halos of different masses ( $8\le \log {{M}_{\rm vir}}\le 10.5$ ) contribute comparably per logarithmic mass interval to reionization, and a photon production rate-averaged escape fraction ( ${\left<f_{\rm esc}\right>}(t)$ ) shows a weak dependence on redshift in the range examined (see Kuhlen & Faucher-Giguère 2012).

3.2. Escape Fraction Enhanced by Runaway OB Stars

Ionizing photons cannot only escape from their birth clouds by destroying them through feedback processes, but also emerge from runaway OB stars displaced from the birth clouds. If we take the typical velocity of the runaway OB stars ∼40 km s⁻¹ (Stone 1991; Hoogerwerf et al. 2001; Tetzlaff et al. 2011), they could travel a distance of ∼200 pc in 5 Myr. Conroy & Kratter (2012) examined the possible ramification of the inclusion of the runaway OB stars using a simple analytic formulation, and concluded that f_esc can be enhanced by a factor of up to 4.5 from ${f_{\rm esc}}\approx 0.02\hbox{--}0.04$ to ${f_{\rm esc}}\approx 0.06\hbox{--}0.18$ in halos of mass $10^8 \lesssim {{M}_{\rm vir}}\lesssim 10^{9}\,{{M}_\odot }$ . Given the complexity of the ISM dynamics (e.g., McKee & Ostriker 2007), it would seem prudent to examine this issue in greater details in realistic environments. To do so, we have performed a twin cosmological simulation of the ${\mathsf {FR}}$ run by designating 30% of mass in each stellar particle as a separate runaway particle and dynamically follow their motion.

Figure 8 shows an example of the difference in environment between runaway and nonrunaway particles in a galaxy in a 3 × 10¹⁰ M_☉ halo at z = 7. At this redshift, the central galaxy shows ${f_{\rm esc}}=0.14$ . The average hydrogen number density for runaways younger than 5 Myr (n_H ∼ 130 cm⁻³) is found to be roughly 20 times smaller than that of nonrunaways (n_H ∼ 3000 cm⁻³). Given that these stars will explode in the next 5–10 Myr, the fact that the local density of some runaway OB stars is smaller than nonrunaways suggests that the impact from SN explosions will be enhanced. Indeed, we find that the stellar mass of the galaxies in halos of mass ${{M}_{\rm vir}}\gtrsim 10^9\,{{M}_\odot }$ is smaller by a factor of 1.7 on average, compared with that from the ${\mathsf {FR}}$ run (see Figure 3). For galaxies in smaller halos, there is no clear hint that the runaway OB stars help suppress the star formation. This is partly because runaway OB stars cannot only provide energy but also distribute metals more efficiently, which can increase the cooling rate in halos. Comparison of the temperature distribution between the two runs further substantiates the claim that runaway OB stars help regulate the star formation (Figure 9). The volume of T ⩾ 10⁵ K gas inside the zoomed-in region in the ${\mathsf {FRU}}$ run (≈7 kpc³, physical) is 30% larger than that in the ${\mathsf {FR}}$ run.

**Figure 9.** Comparison of the temperature distribution in the run without (top, ${\mathsf {FR}}$ ) and with (bottom, ${\mathsf {FRU}}$ ) runaway OB stars at z = 10.2. The white bar measures 100 kpc (proper). The ${\mathsf {FRU}}$ run shows bigger hot bubbles (30%) with T ⩾ 10⁵ K than the ${\mathsf {FR}}$ run, suggesting that runway OB stars affect the regulation of star formation.
Download figure:
Standard image High-resolution image

The left panel in Figure 10 shows the instantaneous f_esc measured at three different redshifts from the ${\mathsf {FRU}}$ run. Again, less massive galaxies tend to exhibit a higher f_esc, which can be attributed to the fact that star formation in smaller halos is more easily affected by the energetic explosions. As expected, the inclusion of the runaway OB stars increases the instantaneous escape fraction on average. The photon production rate-weighted average of f_esc (right panel in Figure 10) shows this more clearly. In our fiducial run ( ${\mathsf {FR}}$ ), 11.4% of the ionizing photons produced escaped from the halos of mass ${{M}_{\rm vir}}\ge 10^8\,{{M}_\odot }$ at z ⩾ 7. On the other hand, the ${\mathsf {FRU}}$ run yields higher ${\left<f_{\rm esc}\right>}$ of 13.8%, which is enhanced by 22% compared with that of the ${\mathsf {FR}}$ run. Although this increase is not as large as claimed in Conroy & Kratter (2012), the contribution from the runaway OB stars is certainly significant. Similarly, as in the ${\mathsf {FR}}$ run, no clear dependence of ${\left<f_{\rm esc}\right>}$ on halo mass is found.

It is interesting to discuss possible origins of the significantly different enhancement in the escape fraction due to runaway OB stars found in our simulations compared with the estimate by Conroy & Kratter (2012). First, while their model predicts f_esc of nonrunaways to be about 2%–4% in halos of mass $10^8 \le {{M}_{\rm vir}}\le 10^9 \, {{M}_\odot }$ , we find that the self-regulation of star formation via SN explosions leads to a high escape of ∼10% in our fiducial model ( ${\mathsf {FR}}$ ). Second, while their model finds that runaway OB stars are found to have high ${f_{\rm esc}}$ (=30%–80%), our results imply that the mean escape fraction of ionizing photons from runaway OB stars is about 20% ( $11.4\%\times 70\% + {\it 20\%}\times 30\%\approx 13.8\%$ ). We also make a more elaborate estimate as follows. We measure the optical depth in the Lyman continuum for the gas inside each halo along 768 sightlines per star particle, and combine the attenuated spectral energy distributions. These are used to count the number of hydrogen ionizing photons for runaways and nonrunaways separately. We find that the relative contribution from the runaways to the total number of escaping photons is comparable with that of the non-runaways. Considering that the runaway particle is assumed to explain only 30% of all the OB stars, the net ${f_{\rm esc}}$ for the runaways can be estimated to be roughly 23% (=13.8%/2/0.3). This is twice higher chance of escaping than the nonrunaways, but much smaller than computed in the analytic model. If the escape fraction of nonrunaway OB stars were 2% in our simulations, the total escape fraction would become 2% × 70% + 23% × 30% = 8.3%, corresponding to an increase of a factor of 4.2. It is thus clear that most of the discrepancies arise in a large part due to different escape fraction values for nonrunaway OB stars and also due to different escape fraction values for runaway OB stars.

Although ${\left<f_{\rm esc}\right>}$ is 22% larger in the ${\mathsf {FRU}}$ run than ${\mathsf {FR}}$ , the cumulative number of photons escaped in halos with ${{M}_{\rm vir}}\ge 10^8\,{{M}_\odot }$ by z = 7 (N_ion ≈ 1.3 × 10⁶⁹) is found to be similar to that of the ${\mathsf {FR}}$ run (N_ion ≈ 1.6 × 10⁶⁹). This is because star formation is suppressed in relatively massive halos ( ${{M}_{\rm vir}}\ge 10^9\,{{M}_{\rm vir}}$ ).

One question is whether or not enough photons escape to keep the universe at z ∼ 7 ionized. The critical photon rate density that can balance the recombination of ionized hydrogen is

$\begin{equation} \dot{n}_{\rm ion}^{\rm crit} = \alpha _{\rm B} \, n_e \, n_{{\rm H\,\scriptsize{II}}} \simeq 10^{47.2} C_{{\rm H\,\scriptsize{II}}} (1+z)^3 \, {\rm [s^{-1}\,Mpc^{-3}]}, \end{equation} \tag{ 2 }$

where α_B is the case B recombination coefficient, n_e is the number density of electron, $n_{{\rm H\,\scriptsize{II}}}$ is the number density of ionized hydrogen, and $C_{{\rm H\,\scriptsize{II}}} \equiv \left<n_{{\rm H\,\scriptsize{II}}}^2\right>/\left<n_{{\rm H\,\scriptsize{II}}}\right>^2$ is the clumping factor of ionized gas. For a choice of the clumping factor $C_{{\rm H\,\scriptsize{II}}}\sim 3$ (Pawlik et al. 2009; Raičević et al. 2011) and the temperature T = 20000 K, $\dot{n}_{\rm ion}^{\rm crit} = 10^{50.4}\,[(1\,{+}\,z)/8]^3 \, {\rm s^{-1}\,Mpc^{-3}}$ . Figure 11 shows that the escaped photons in ${\mathsf {FRU}}$ can balance the recombination at z ⩽ 9. We find that the photon rate density at z ∼ 7 is $\dot{n}_{\rm ion}=10^{50.7-50.9} \, {\rm s^{-1}\,Mpc^{-3}}$ , consistent with observational findings. Ouchi et al. (2009) estimated the ionizing photon density to be $\log \dot{n}_{\rm ion} \simeq 49.8\hbox{--}50.3$ by integrating the UV luminosity function (UVLF) down to M_UV = −18 (lower) or L = 0 (upper estimate) with a slope of α = −1.72 at z ∼ 7 with ${f_{\rm esc}}=20\%$ . If the slope found in the more recent literature (McLure et al. 2013), α = −1.90, is used, the maximum photon rate density derived would increase to $\log \dot{n}_{\rm ion} \simeq 50.8$ , which is in agreement with our estimation. Note that the photons escaping from halos of mass ${{M}_{\rm vir}}\ge 10^8\,{{M}_\odot }$ account for more than 90% of the total escaping photons if the baryon-to-star conversion efficiency derived in our simulation is extrapolated to smaller halos ( ${{M}_{\rm vir}}<10^{8.5}\,{{M}_\odot }$ , see below), and hence our results should be compared with the maximum photon rate density. Given that their chosen <f_esc > is closed to what our simulation yields (13.8%), the agreement implies that SFRs of the galaxies are well reproduced in our simulation. Indeed, we find that our simulated UVLF measured at 1500 Å (rest-frame) shows excellent agreement with the LF with the slope of α = −1.90 (McLure et al. 2013) down to M₁₅₀₀ = −13 (Figure 12). Here, we neglect the effect of dust extinction, as the galaxies in our sample are very metal-poor (Z_star ≲ 10⁻³).

**Figure 11.** Balance between the ionizing photons escaping from the dark matter halo and the recombination rate in the ${\mathsf {FRU}}$ run. The thick gray line shows the balance condition when the clumping of $C_{{\rm H\,\scriptsize{II}}}=3$ is used. Enough photons to keep the universe ionized escape from the halo after z ∼ 8.
Download figure:
Standard image High-resolution image

**Figure 12.** Rest-frame ultraviolet luminosity function from the ${\mathsf {FRU}}$ run at z = 7. Error bars denote the Poissonian error. Observational data from Bouwens et al. (2011) and McLure et al. (2013) are shown as the shaded region and empty squares, respectively. Also included as solid and dashed lines are the Schechter fits to the data provided in these studies.
Download figure:
Standard image High-resolution image

In Figure 13, we plot the product of photon number-weighted escape fraction ( ${\left<f_{\rm esc}\right>}$ ) and baryon-to-star conversion efficiency (f_⋆ ≡ Ω_mM_star/Ω_bM_vir) at z = 7. Notice that we include all stars within the virial radius of a dark matter halo in this measurement. Since there is little evolution in ${\left<f_{\rm esc}\right>}$ with redshift (Figure 10, right panel), we combine ${\left<f_{\rm esc}\right>}$ of the halos in the same mass range at 7 ⩽ z < 20 to obtain the mean escape fraction as a function of halo mass (Table 2). We then use a simple fit to the mean, as

$\begin{equation} \log {\left<f_{\rm esc}\right>}({{M}_{\rm vir}}) \approx -0.510 - 0.039 \log {{M}_{\rm vir}}. \end{equation} \tag{ 3 }$

We limit our fit to the sample with ${{M}_{\rm vir}}\ge 10^{8.5}\,{{M}_\odot }$ , where each halo is resolved with ∼2000 dark matter particles and more. There is a trend that more massive halos contribute more to the total number of ionizing photons per mass, which essentially reflects the fact that low-mass halos are inefficient in forming stars (see also Figure 3). The average ${\left<f_{\rm esc}\right>}f_{\star }$ of different halo masses can be fitted with

$\begin{equation} \log {\left<f_{\rm esc}\right>}f_{\star } \approx -7.342 + 0.474\, \log {{M}_{\rm vir}}, \end{equation} \tag{ 4 }$

shown as the red dashed line in Figure 13. We note that ${\left<f_{\rm esc}\right>}f_{\star }$ becomes as low as ∼5 × 10⁻⁴ in small halos ( ${{M}_{\rm vir}}\sim 10^{8.5}\,{{M}_\odot }$ ), which is roughly 40 times smaller than the results from Wise & Cen (2009; ${\left<f_{\rm esc}\right>}f_\star \approx 0.02$ ). The difference can be attributed to two factors. First, our <f_esc > is smaller by a factor of ∼three to four than that of Wise & Cen (2009). This is probably due to the fact that their cosmological runs start from the initial condition extracted from adiabatic simulations in which no prior star formation is included. Since radiative cooling and star formation are suddenly turned on at some redshift, the gas in the halo rapidly collapses and forms too many stars in their cosmological runs. This is likely to have resulted in stronger starbursts in the galaxies, leading to a higher escape probability. Second, because of the same reason, f_⋆ is considerably higher in the Wise & Cen (2009) halos than in our halos. For halos of masses with ${{M}_{\rm vir}}\sim 10^{8.5}\,{{M}_\odot }$ , we find that f_⋆ ≈ 0.003, which is smaller by a factor of ∼10 than those in Wise & Cen (2009). Indeed, we find fairly good agreement with the latest determination of ${\left<f_{\rm esc}\right>}f_{\star }$ in halos of ${{M}_{\rm vir}}\sim 10^{8.5}$ by Wise et al. (2014), who model star formation self-consistently in their cosmological radiation hydrodynamics simulations.

**Figure 13.** Product of the stellar mass fraction within the virial radius of a dark matter halo ( $f_\star =\Omega _{\rm m} {{M}_{\rm star}}/ \Omega _{\rm b} {{M}_{\rm vir}}$ ) at z = 7 and halo mass-dependent photon production rate-averaged escape fraction from the cosmological simulation with runaway OB stars ( ${\mathsf {FRU}}$ ). Averages are shown as red empty squares, with the simple regression (dashed line). A smaller number of photons is escaped per unit mass in smaller halos, reflecting the results that star formation is inefficient in the low-mass halos.
Download figure:
Standard image High-resolution image

Table 2. Photon Number-weighted f_esc at 7 ⩽ z ≲ 15 from the FRU Run

log M_vir	<f_esc >
8.25	0.144 ± 0.038
8.75	0.146 ± 0.064
9.25	0.148 ± 0.077
9.75	0.128 ± 0.069
10.25	0.113 ± 0.079

Download table as: ASCII Typeset image

It is worth mentioning that adopting high spatial resolution (or gravitational softening length) is important to accurately predict the escape fraction. If the resolution is not high enough to capture the rapid collapse of gas clouds, the resulting star formation histories would become less episodic, leading to a longer time delay between the peak of star formation and escape fraction. This in turn would reduce the fraction of escaping photons. To examine this issue, we run two additional simulations with the identical initial condition and other parameters, but with one less or more level of refinement, corresponding to 8.5 pc or 2.1 pc (physical) resolution, respectively. We find that the run with the lower resolution yields a factor of two smaller mean escape fraction at z = 9 ( ${\left<f_{\rm esc}\right>}=7.6\%$ ; see Appendices A–C). On the contrary, higher resolution run exhibits a comparable mean escape fraction of ${\left<f_{\rm esc}\right>}=13.9\%$ at z = 10, suggesting that the results are reasonably converged for the parameters used in the ${\mathsf {FRU}}$ run.

4. DISCUSSION

Recent studies show that the escape fraction should be larger than 20% to re-ionize the universe by z = 6 matching the Thomson optical depth inferred from the CMB (Kuhlen & Faucher-Giguère 2012; Shull et al. 2012; Robertson et al. 2013). This can be obtained by numerically solving the simple differential equation for the H ii bubble

$\begin{equation} \frac{d Q_{{\rm H\,\scriptsize{II}}} }{dt} = \frac{\dot{n}_{\rm ion}}{\left<n_{\rm H}\right>} - \frac{Q_{{\rm H\,\scriptsize{II}}}}{t_{\rm rec}(C_{{\rm H\,\scriptsize{II}}})}, \end{equation} \tag{ 5 }$

where $Q_{{\rm H\,\scriptsize{II}}}$ is the volume filling fraction of the bubble, <n_H > is the comoving mean density of the universe, and $t_{\rm rec} (C_{{\rm H\,\scriptsize{II}}})= [ C_{{\rm H\,\scriptsize{II}}}\,\alpha _{\rm B}(T)\, f_e\, \left<n_{\rm H}\right> \,(1+z)^3]^{-1}$ is the recombination timescale for a given clumping factor and temperature. Here f_e is a correction factor that accounts for the additional contribution of singly (z > 4) or doubly (z < 4) ionized helium to the number density of electron (e.g., Kuhlen & Faucher-Giguère 2012). We adopt a redshift-dependent clumping factor of $C_{{\rm H\,\scriptsize{II}}} = 1 + \exp (-0.28\, z \,+\,3.59)$ at z ⩾ 10 or $C_{{\rm H\,\scriptsize{II}}} = 3.2$ at z < 10 following Pawlik et al. (2009). Once $Q_{{\rm H\,\scriptsize{II}}}$ is determined, the Thomson optical depth as a function of redshift can be calculated as

$\begin{equation} \tau _e (z)= \int _0^z c \left<n_{\rm H}\right>\,\sigma _T\,f_e\,Q_{{\rm H\,\scriptsize{II}}}(z^{\prime }) \frac{(1+z^{\prime })^2 dz^{\prime }}{H(z^{\prime })}, \end{equation} \tag{ 6 }$

where σ_T is the Thomson electron cross section, and H(z) is the Hubble parameter. We follow the exercise by using the ionizing photon density from Figure 11 to examine whether our models provide a reasonable explanation for the reionization history. For $\dot{n}_{\rm ion}$ at z < 7, we extrapolate based on the simple fit to the results in Figure 11. This simple experiment indicates that the universe can be re-ionized by z = 7.25. However, the evolution of the photon density from the ${\mathsf {FRU}}$ run predicts a smaller volume filling fraction of the H ii bubble at z = 10 ( $Q_{{\rm H\,\scriptsize{II}}}=12\%$ ), compared with other analytic models ( $Q_{{\rm H\,\scriptsize{II}}}\gtrsim 20\%$ ; e.g., Shull et al. 2012) that could reproduce the CMB measurement (τ_e ∼ 0.09; Komatsu et al. 2011). Consequently, the FRU run yields the Thomson optical depth of τ_e = 0.065, which is consistent only within 2σ with the CMB measurement. This implies that more ionizing photons are required to escape from halos at high redshift to explain the reionization history of the universe.

The deficiency of ionizing photons may in part be attributed to the fact that our simulations cannot resolve the collapse of small-mass halos ( ${{M}_{\rm vir}}\lesssim 10^8\,{{M}_\odot }$ ) due to finite mass resolution. Paardekooper et al. (2013) argue that reionization is driven by dwarf-sized halos of masses ${{M}_{\rm vir}}=10^7\hbox{--}10^8\,{{M}_\odot }$ with high ${\left<f_{\rm esc}\right>}$ of ≈0.4–0.9. Similarly, Wise et al. (2014) find that the ionizing photons from the minihalos with ${{M}_{\rm vir}}=10^{6.25}\hbox{--}10^{8.25}\,{{M}_\odot }$ is crucial at reproducing the Thompson optical depth from the CMB measurements. In order to examine the importance of the minihalos in light of our new results, we estimate the optical depth as a function of the minimum halo mass that can contribute to reionization. To do so, we use the theoretical halo mass functions at different redshifts (Jenkins et al. 2001), convolved with the baryon-to-star conversion efficiency measured at z = 7 from the ${\mathsf {FRU}}$ run for ${{M}_{\rm vir}}\ge 10^{7.5}\,{{M}_\odot }$ and Wise et al. (2014) for ${{M}_{\rm vir}}<10^{7.5}\,{{M}_\odot }$ (orange line in the middle panel of Figure 14), to derive the increase in the stellar mass density with redshift. The number of escaping ionizing photons is then calculated by multiplying the number of photons produced with the halo mass-dependent escape fraction based on our results and Wise et al. (2014), as (Figure 14, top panel)

$\begin{equation} \log {\left<f_{\rm esc}\right>}=\left\lbrace \begin{array}{@{}ll@{}}-0.51 - 0.039\,\log {{M}_{\rm vir}}&\! (\log {{M}_{\rm vir}}\ge 8.5) \\ 2.669 - 0.413\,\log {{M}_{\rm vir}}&\! (7 \le \log {{M}_{\rm vir}}< 8.5) \\ -0.222 &\! (\log {{M}_{\rm vir}}< 7)\\ \end{array} \right.. \end{equation} \tag{ 7 }$

We neglect the contribution from rare massive halos with ${{M}_{\rm vir}}>10^{12}\, {{M}_\odot }$ . Figure 14 (orange line, bottom panel) shows that the minihalos of ${{M}_{\rm vir}}<10^7\,{{M}_\odot }$ can indeed provide enough photons to match τ_e inferred from the CMB measurement. While the ionizing photons from ${{M}_{\rm vir}}>10^7\,{{M}_\odot }$ only gives τ_e = 0.072, the additional photons arising from the minihalos augment the optical depth to 0.122. However, we note that this sensitively depends on the assumption on the baryon-to-star conversion efficiency in the minihalos. For example, when the stellar-mass–halo-mass relation found in the ${\mathsf {FRU}}$ is extrapolated to the minihalos (blue line in the bottom panel), the optical depth for the entire halos is only τ_e = 0.073. Given that these minihalos would host a handful of star particles with m_star ∼ 10²–10³ M_☉ in current numerical simulations, it is unclear how the mass resolution affects the conversion efficiency, and further investigations on star formation in the minihalos will be useful to better understand their relative role to the total ionizing budget.

In our simulation, we approximate that massive stars (M > 8 M_☉) evolve off and explode after 10 Myr. We note that this is roughly the timescale of the delay between the peak of star formation and escape fraction. In reality, the SN can emerge as early as ∼3 Myr for a simple population (Schaller et al. 1992). Stellar winds, photo-ionization, and radiation pressure acting on electron and dust can come into play even earlier. Walch et al. (2012) claim that a 10⁴ M_☉ molecular cloud of the radius 6.4 pc can be dispersed on a 1–2 Myr timescale by the overpressure of H ii regions. Moreover, it is also plausible that the ionization front instabilities may lead to the higher escape probability of ionizing photons (Whalen & Norman 2008). If these mechanisms played a role in shaping the evolution of individual molecular clouds, the escape fraction measured in our simulations would have been higher than 14%. In this regard, our photon number-weighted mean is likely to represent the minimum escape of ionizing photons. When a higher <f_esc > of 30% is assumed for the star formation history in the ${\mathsf {FRU}}$ run, dark matter halos of ${{M}_{\rm vir}}>10^8\,{{M}_\odot }$ alone can achieve τ_e = 0.076, suggesting that a more precise determination of the escape fraction is as equally important as resolving ultra-faint galaxies with M₁₅₀₀ > −13. Future studies focusing on the interplay between the feedback processes will shed more light on the reionization history of the Universe.

5. CONCLUSIONS

The escape fraction of hydrogen ionizing photons is a critical ingredient in the theory of reionization. Despite its importance, only a handful of studies examined the escape fraction ( ${f_{\rm esc}}$ ) of high-z galaxies in a cosmological context (Wise & Cen 2009; Razoumov & Sommer-Larsen 2010; Yajima et al. 2011; Paardekooper et al. 2013; Wise et al. 2014). To better understand the physics behind the escape of ionizing photons and quantify f_esc, we have carried out two zoomed-in cosmological radiation hydrodynamics simulations of 3.8 × 4.8 × 9.6 Mpc³ box (comoving) with the Ramses code (Teyssier 2002; Rosdahl et al. 2013) with high spatial (∼4 pc, physical) and stellar mass resolution of 49 M_☉. Because energy-based feedback from SN explosions suffers from the artificial radiative cooling if the cooling length is under-resolved, we have implemented a new mechanical feedback scheme that can approximate all stages of a SN explosion from the free expansion to snowplow phase. With the physically based feedback model, we have investigated the connection between the regulation of star formation and corresponding evolution of the escape of ionizing photons. We have also explored the relative importance of runaway OB stars to the escape fraction by comparing the twin simulations with ( ${\mathsf {FRU}}$ ) and without ( ${\mathsf {FR}}$ ) runaways. Our findings can be summarized as follows.

1.
When a dense cloud begins to form a cluster of stars, the escape fraction is negligible. As energetic explosions by massive stars follow after ∼10 Myr, it blows the star-forming gas away, increasing the instantaneous escape fraction (f_esc) to ≳10%. Although f_esc is kept high in this phase, subsequent star formation is markedly suppressed, and only a small number of photons escapes from their host dark matter halo (Figure 4). This time delay between the peak of star formation and the escape fraction is crucial in predicting the actual escape probability of ionizing photons. While the instantaneous f_esc can easily attain ≳ 30% in halos of mass ${{M}_{\rm vir}}\ge 10^8\,{{M}_\odot }$ on average (Figure 5), the photon number-weighted mean of the escape fraction (<f_esc >) is found to be 11.4% (Figure 6).
2.
f_esc tends to be higher in less massive halos and at lower redshift for a give halo mass (Figure 5). This is essentially because less dense and smaller galaxies are more susceptible to SN explosions. However, the photon production rate-averaged escape fractions show no clear dependence on halo mass and redshift, again implying that the interplay between star formation and the delay in the onset of negative feedback is more important in determining the actual escape probability.
3.
Absorption of ionizing photons by neutral hydrogen in the ISM is significant (Figure 7). For galaxies with a low escape fraction ( ${f_{\rm esc}}<10\%$ ), the effective optical depth by the gas within 100 pc from each young star particles is found to be τ_{eff, 100 pc} ∼ 1.9–3.8 at z ∼ 8. The nearby neutral gas alone can reduce the number of ionizing photons by 7–45 in this case, demonstrating the importance of properly resolving the ISM to predict a more accurate escape fraction.
4.
Our physically based SN feedback effectively regulates star formation. Only 0.1% to 10% of the baryons are converted into stars in galaxies at z = 7 (Figure 3). The energetic explosions sometimes completely shut down star formation when galaxies are small. The baryon-to-star conversion ratio is smaller in less massive halos. Consequently, halos of different masses contribute comparably to the total number of ionizing photons escaped by z = 7 (Figure 6).
5.
Inclusion of runaway OB stars increases the escape fraction to ${\left<f_{\rm esc}\right>}=13.8\%$ from 11.4% (Figure 10). Since the runaway OB stars tend to move to lower density regions, photons from them have a higher chance of escaping. Moreover, as the runaway OB stars explode in a less dense medium, feedback from SNe becomes more effective, resulting in reduced star formation in halos ${{M}_{\rm vir}}\ge 10^9\,{{M}_\odot }$ , compared with the ${\mathsf {FR}}$ run. Because of the balance between the increase in <f_esc > and the decrease in star formation, the total number of ionizing photons escaped by z = 7 is found to be comparable in the two runs.
6.
A sufficient amount of photons escape from the dark matter halos with ${{M}_{\rm vir}}\ge 10^8\,{{M}_\odot }$ to keep the universe ionized at z ⩽ 9. The simulated UV luminosity function with a faint end slope of −1.9 is consistent with observations.

We thank an anonymous referee for constructive suggestions that improved this paper. We are grateful to Julien Devriendt, Sam Geen, Chang-Goo Kim, Eve Ostriker, Adrianne Slyz, and John Wise for insightful discussions. Special thanks go to Romain Teyssier and Joakim Rosdahl for sharing their radiation hydrodynamics code with us. Computing resources were provided in part by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and in part by Horizon-UK program through DiRAC-2 facilities. The research is supported by NSF grant AST-1108700 and NASA grant NNX12AF91G.

APPENDIX A: A NEW PHYSICAL SCHEME FOR SUPERNOVA FEEDBACK

It is well established that feedback from SN explosions is crucial to understanding many aspects of galaxy evolution (e.g., Dekel & Silk 1986). Detailed implementation of SN feedback has progressed over time. Early works have included it by depositing internal energy into the parent cell of the SN particle, but the energy is found to be radiated away quickly, making the feedback ineffective (e.g., Katz 1992; Abadi et al. 2003; Slyz et al. 2005; Hummels & Bryan 2012). To alleviate the problem, Cen & Ostriker (2006) distributed the thermal energy over 27 neighboring cells weighted by specific volume, which has the virtue of being able to mimic propagation of an explosion that is preferentially channelled into more diffuse regions. Some groups have adopted unconventional schemes to regulate star formation and gas dynamics. These include increasing the total energy from the SN (Thacker & Couchman 2000), making the SN event more episodic and stronger (Scannapieco et al. 2006; Dalla Vecchia & Schaye 2012), decoupling hydrodynamic interactions of super-wind particles with the ambient medium (Oppenheimer & Davé 2006), turning off gas cooling after deposition of SN energy for a significant period of time (Mori et al. 1997; Stinson et al. 2006; Governato et al. 2007), or disabling cooling when the turbulence measured locally is significant (Teyssier et al. 2013). Alternatively, the explosion can be provided in kinetic form (Navarro & White 1993; Dubois & Teyssier 2008; Dalla Vecchia & Schaye 2008). This is not immune to the artificial cooling problem, because the kinetic energy is converted into heat through shocks immediately after launching given the typical high Mach number, ambient density, and limited numerical resolution.

So far, there is no satisfactory scheme for modeling the SN feedback that is robust and does not strongly depend on simulation resolution. We have devised a new physical scheme that is reasonably accurate. Before going into details, it is useful to gain a basic physical understanding of the Sedov explosion. The Sedov explosion consists of four stages. In the first stage, the SN ejecta sweep up an insignificant amount of mass compared with the initial ejecta mass, and both energy and momentum are conserved. In the second stage, the swept-up ISM mass is now comparable to or exceeds the initial ejecta mass so the explosion now enters a self-similar phase. Cooling has not set in so energy is conserved up to this stage. As a result, the radial outward momentum increases as the square root of the total shell mass. The third stage is the cooling phase, which is normally very brief in an isobaric gas. The last stage is called the snow-plow phase, when the total linear radial outward momentum is conserved. In essence, the problem of excessive cooling in some simulations is a result of not giving the surrounding gas an adequate amount of momentum that is commensurate with the stage it is supposed to be in, given the physical conditions. For example, in some cosmological simulations with limited resolution, a single cell is already larger than the expected radius at which the snow-plow phase commences. In this case, a physically correct way would be to deposit the maximum expected momentum, as at the end of the cooling phase (i.e., the end of the third stage), to the surrounding cells. If one instead deposits thermal energy with no momentum or a combination of thermal energy and an inadequate amount of momentum, the final momentum would fall short of the true expected momentum. We now describe our new physical scheme in detail.

The momentum input from SN explosions is modeled as

$\begin{equation} dp=\left\lbrace \begin{array}{@{}cc@{}}\sqrt{2 \chi (\Omega) M_{\rm ej}\,f_e\,E_{\rm SN}}\times \frac{d\Omega }{4\pi } & \left[ \chi (\Omega) \le \chi _{\rm tr}(\Omega)\right]\\ p_{\rm SN}(E_{\rm SN},\,n_{\rm H}, Z) \times \frac{d\Omega }{4\pi } & \left[\chi (\Omega) > \chi _{\rm tr}(\Omega)\right] \\ \end{array} \right., \end{equation} \tag{ A1 }$

where dΩ is the solid angle subtended by a neighboring cell, E_SN(= N_SN 10⁵¹ erg) is the total explosion energy, M_ej is the total ejecta mass, n_H and Z are the hydrogen number density and gas metallicity of the cell into which the blast wave propagates, and $f_e = 1-\frac{\chi -1}{3(\chi _{\rm crit}-1)}$ is introduced to smoothly connect the two regimes. The most important quantity in this equation is

$\begin{equation} \chi \equiv dM_{\rm shell}(\Omega) / dM_{\rm ej}(\Omega), \end{equation} \tag{ A2 }$

which is the ratio of the shell (ejecta plus mass swept up) to ejecta mass along some solid angle (Ω, Ω + dΩ; see Equation (A6)). The first part of Equation (A1) is valid up to the cooling phase and the second after that. A new and critical element is that we differentiate the energy-conserving and momentum-conserving phase by introducing a mass ratio at the transition (χ_tr). In actual implementations, the swept-up mass in direction Ω is taken as the sum of the gas mass in the adjacent cell in direction Ω and some fraction of the gas mass in the SN cell (see below).

A simple estimate of momentum budget per SN at the free expansion phase is $\sqrt{2 m_{\rm SNII} e_{\rm SN}} \approx 3.9\times 10^4 \, {\rm km\, s^{-1}} \,{{M}_\odot }$ if the typical SN progenitor mass of $\bar{m}_{\rm SNII}=15.2\, {{M}_\odot }$ on the zero-age main sequence, appropriate for the Chabrier IMF (Chabrier 2003) with the lower (upper) mass of 8 (100) M_☉, is used. On the other hand, the momentum of the SN bubble at the end of the adiabatic phase is much higher (Chevalier 1974; Cioffi et al. 1988; Blondin et al. 1998),

$\begin{equation} p_{\rm SN} (E,\,n_{\rm H}) \approx 3\times 10^5\, {\rm km\, s^{-1}} \,{{M}_\odot }\, E_{51}^{16/17} n_{\rm H}^{-2/17}, \end{equation} \tag{ A3 }$

where E₅₁ is the energy in the unit of 10⁵¹ erg. In addition, Thornton et al. (1998) showed that the momentum input is a decreasing function of the metallicity of ambient medium (p_SN∝f(Z) ≡ max [Z/Z_☉, 0.01]^−0.14). Below 0.01 Z_☉ atomic cooling is primarily responsible for radiative energy loss, and no dependence with metallicity is seen. Motivated by this, we take the momentum

$\begin{equation} p_{\rm SN} (E,\,n_{\rm H},\,Z) \approx 3\times 10^5\, {\rm km\, s^{-1}} \,{{M}_\odot }\, E_{51}^{16/17} n_{\rm H}^{-2/17} f(Z), \end{equation} \tag{ A4 }$

during the momentum-conserving phase (i.e., χ > χ_tr). It is worth noting that the momentum transfer from SN explosions at the snowplow phase is essential to understanding the self-regulation of star formation in the local ISM (e.g., Shetty & Ostriker 2012; Kim et al. 2013).

The mass ratio at the transition (χ_tr) is estimated by equating Equation (A4) to $\sqrt{2\,\chi _{\rm tr} \, M_{\rm ej} \, N_{\rm SN} e_{\rm SN,tr}}$ , where e_{SN, tr} ≈ 6.76 × 10⁵⁰ erg is the kinetic energy at the transition (Blondin et al. 1998), and N_SN is the total number of SN events in a cell. This yields

$\begin{eqnarray} \chi _{\rm tr} & = \frac{900}{0.676 \left(\bar{m}_{\rm SNII }/{{M}_\odot }\right)} \, E_{\rm 51}^{-2/17} n_0^{-4/17} \, f^2(Z) \nonumber \\ & \simeq 87.54 \, E_{\rm 51}^{-2/17} n_0^{-4/17} \, Z^{\prime -0.28}, \end{eqnarray} \tag{ A5 }$

where Z' = max(Z/Z_☉, 0.01). Note that Equation (A1) is correct when energy is conserved, and provides a good approximation until χ = χ_tr. On the other hand, all the momentum available (p_SN) will be transferred to the surrounding gas when SNe explode in a cell in which a large amount of gas exists (χ > χ_tr).

Along with the momentum, mass (ejecta plus mass swept up) and energy are added to the neighboring cells. To illustrate this, let us suppose that a cell in which a SN sits is surrounded by further refined cells (Figure 15). In this case, the total number of adjacent cells except for the ones near 8 vertices is 48 in three-dimensional space. Note that even when the neighbors are not further refined, the following scheme can be applied pretending that the cells are composed of refined cells with the same physical properties. We assume that the gas mass entrained from the SN cell and ejecta are evenly distributed (by the volume) to the cells that are directly affected by the SN (shaded region in Figure 15). The total shell mass entrained in a neighboring cell is then

$\begin{equation} dM_{\rm shell} = \rho _i \Delta x_i^3 + (1-\beta _{\rm sn}) \left(\rho _{\rm sn} \Delta x_{\rm sn}^3 + M_{\rm ej} \right) \frac{d\Omega }{4\pi }, \end{equation} \tag{ A6 }$

where ρ_sn and Δx_sn are the density and size of the SN cell, respectively, M_ej = α ∑m_⋆ is the ejecta mass from SN explosions in the SN cell, and α ≃ 0.317 is the mass fraction of 8–100 M_☉ stars for a simple stellar population with the Chabrier IMF. The corresponding stellar ejecta in the neighboring cell is dM_ejecta = (1 − β_sn)M_ejdΩ/4π. We take a simple approximation that dΩ/4π = 1/48. In Equation (A6), β_sn determines how much fraction of the gas entrained is left behind in the SN cell. In the case of unigrid or smoothed particle hydrodynamics (SPH) simulations, this number can simply be chosen as the ratio of the volume taken by the SN cell and total volume that are directly affected, i.e., β_sn = V_sn/(V_sn + V_neighbor). We take β_sn = 4/56 so that mass is evenly distributed when the neighboring cell has the same level of refinement as the SN cell.

**Figure 15.** Schematic plot of a grid structure around a star particle that undergoes an SN explosion. Shaded regions display cells that are directly affected by the SN explosion. The cell with the light gray color represents a low-density gas with χ < χ_tr into which a smaller amount of momentum is deposited than other dense cells.
Download figure:
Standard image High-resolution image

Note that this scheme can easily be implemented in SPH simulations by computing the local mass loading through some solid angle for the N nearest particles. While this paper was being written, we came to know that a similar approach for SN feedback is implemented by Hopkins et al. (2013), based on the cooling radii of individual blast waves, $r_{\rm cool}\approx {\rm 28\,pc}\, E_{51}^{0.29} n_0^{-0.43}\,f(Z)$ . It is useful to point out that there are two significant differences between the two models. First, the momentum input at solar metallicity is about twice larger in Hopkins et al. (2013) than the measurements from a set of high-resolution hydrodynamics simulations of SN explosions that we use (Thornton et al. 1998). Second, their input momentum has a steeper dependence on metallicity (p_SN∝Z^−0.27) than that our case (p_SN∝Z^−0.18; Thornton et al. 1998). We note that the combination may result in the overestimation of the impact of SN explosions by a factor of ∼2.4 for the gas with Z ≈ 0.01–0.1 Z_☉ in their case.

Finally, we ensure that the total energy is conserved before and after the momentum injection; the residual (total minus kinetic) surplus energy is added as thermal energy to the affected cells in our case.

In order to examine how the new model compares with the standard energy feedback scheme, we perform idealized simulations by placing a SN in a uniform medium of number density n_H = 10 cm⁻³ with solar metallicity. The radiative cooling is included in all calculations. The size of the simulated box is set to 128 pc (or 256 pc for the coarsest resolution run) while increasing the number of cells from 16³ to 512³. Then we measure the radial momentum from the explosion at the momentum-conserving phase (0.1 Myr; Thornton et al. 1998). Note that the mass swept up by SN ejecta is different depending on the resolution in our new feedback model, which in turn change the input momentum to the adjacent cells at the time of explosion. For example, χ varies from 1.0006 to 164.3 for the runs with 0.25–16 pc, respectively. The corresponding momentum input based on Equation (A1) is 1.0003–7.136 times $\sqrt{2m_{\rm SNII} e_{\rm SN}}$ , and thus the different resolution runs correspond to different stages of the explosion from the adiabatic to momentum-conserving phase in practice. Figure 16 shows that approximately the same amount of momentum is transferred to the surrounding medium with our new mechanical feedback scheme. For comparison, we also run the same set of simulations with thermal feedback which distributes the SN energy into the 27 surrounding cells. We find that although the amount of momentum at the snowplow phase is roughly 20% smaller than the prediction by 1D hydrodynamics simulations (Thornton et al. 1998), the results from the thermal and our mechanical model are very similar when high resolution is employed (Δx ≲ 4 pc). On the other hands, the thermal feedback exhibits the well-known overcooling problem in lower resolution runs (Δx > 8 pc). We also perform the same experiment in gas with n_H = 100 cm⁻³, and find that at least 1pc resolution is required to properly model the momentum transfer with thermal feedback, while the momentum transfer is again not sensitive to the resolution with our mechanical feedback scheme.

**Figure 16.** Momentum transfer from a single SN event in a uniform medium of density n_H = 10 cm⁻³ and solar metallicity as a function of different spatial resolution (Δx_min). Note that the radiative cooling is included in this calculation. The radial momentum is measured at the momentum-conserving phase (0.1 Myr). The dashed line shows the momentum in the shell predicted from the one dimensional hydrodynamic calculation by Thornton et al. (1998). Orange stars and green circles show the momentum transfer from the explosion with our new model and thermal feedback often used in the literature, respectively.
Download figure:
Standard image High-resolution image

APPENDIX B: SUPPRESSION OF STAR FORMATION IN LOW-MASS HALOS

In the main text, we show that stars form inefficiently in a low-mass halo with $10^8\lesssim {{M}_{\rm vir}}\lesssim 10^9\,{{M}_\odot }$ , and attribute this to the fact that SN feedback is very effective. However, photoionizing background radiation may also prevent gas from collapsing into the small halos (e.g., Shapiro et al. 1994; Thoul & Weinberg 1996). In order to substantiate our claim, we carry out a hydrodynamics+N-body version of the ${\mathsf {FRU}}$ run down to z = 9, with a radiative transfer module turned off. Other physical and numerical parameters including the initial condition are kept fixed. Note that we expect a stronger impact from reionization in the ${\mathsf {FRU}}$ run, as the escape of ionizing photons is more significant than in the ${\mathsf {FR}}$ run. Figure 17 shows that only a small fraction of baryons (<0.01Ω_b/Ω_m) is converted into stars in the small halos even in the absence of ionizing radiation, demonstrating that SN feedback is primarily responsible for the low conversion efficiency found in our simulations.

**Figure 17.** Comparison of the galaxy stellar mass from the simulations with (orange) and without (green) ionizing radiation at z = 9. This plot shows that the low baryon-to-star conversion efficiency of the central galaxy in our small mass halos ( ${{M}_{\rm vir}}\lesssim 10^9\,{{M}_\odot }$ ) is not mainly due to the photoionizing background radiation, but due to supernova feedback.
Download figure:
Standard image High-resolution image

APPENDIX C: EFFECT OF SPATIAL RESOLUTION

We examine the effect of varying the spatial resolution on the mean escape fraction in Figure 18. All other parameters are fixed as in the ${\mathsf {FRU}}$ run. We run the lower (higher) resolution simulation down to z = 9 (z = 10), and compare the photon number-weighted time average of the escape fraction ( ${\left<f_{\rm esc}\right>}(<t_{\rm H})$ ) with that from the ${\mathsf {FRU}}$ run. While the lower resolution run (8.5 pc, dotted) shows a lower mean escape of ${\left<f_{\rm esc}\right>}=7.6\%$ than the ${\mathsf {FRU}}$ , a similar fraction (13.9%) of ionizing photons is escaped from halos in the higher resolution run (2.1 pc, dot-dashed), indicating that the escape fraction in our fiducial run is reasonably converged for the parameters used. We carefully inspect the cause of the lower escape fraction in the 8.5 pc run, and find that this is due to a longer time delay between the peak of star formation and escape fraction as gas clouds collapse slowly compared with the higher resolution runs (2.1 or 4.2 pc).