# Feedback-regulated star formation and escape of LyC photons from mini-haloes during reionization

Taysun Kimm, Harley Katz, Martin Haehnelt, Joakim Rosdahl, Julien Devriendt, Adrianne Slyz

## ► To cite this version:

Taysun Kimm, Harley Katz, Martin Haehnelt, Joakim Rosdahl, Julien Devriendt, et al.. Feedback-regulated star formation and escape of LyC photons from mini-haloes during reionization. Monthly Notices of the Royal Astronomical Society, 2017, 466 (4), pp.4826-4846. 10.1093/mnras/stx052 hal-01582367

HAL Id: hal-01582367

https://hal.science/hal-01582367v1

Submitted on 28 Jun 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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-regulated star formation and escape of LyC photons from mini-haloes during reionization

Taysun Kimm,<fnref n="1" />★ Harley Katz,<fnref n="1" /> Martin Haehnelt,<fnref n="1" /> Joakim Rosdahl,<fnref n="2" /> Julien Devriendt<fnref n="3" />,<fnref n="4" /> and Adrianne Slyz<fnref n="3" />

Accepted 2017 January 9. Received 2016 December 20; in original form 2016 August 16

## ABSTRACT

Reionization in the early Universe is likely driven by dwarf galaxies. Using cosmological radiation-hydrodynamic simulations, we study star formation and the escape of Lyman continuum (LyC) photons from mini-haloes with $M_{\text{halo}} \lesssim 10^8 M_{\odot}$. Our simulations include a new thermo-turbulent star formation model, non-equilibrium chemistry and relevant stellar feedback processes (photoionization by young massive stars, radiation pressure and mechanical supernova explosions). We find that feedback reduces star formation very efficiently in mini-haloes, resulting in the stellar mass consistent with the slope and normalization reported in Kimm & Cen and the empirical stellar mass-to-halo mass relation derived in the local Universe. Because star formation is stochastic and dominated by a few gas clumps, the escape fraction in mini-haloes is generally determined by radiation feedback (heating due to photoionization), rather than supernova explosions. We also find that the photon number-weighted mean escape fraction in mini-haloes is higher (~20–40 per cent) than that in atomic-cooling haloes, although the instantaneous fraction in individual haloes varies significantly. The escape fraction from Pop III stars is found to be significant (≥10 per cent) only when the mass is greater than ~100 $M_{\odot}$. Based on simple analytic calculations, we show that LyC photons from mini-haloes are, despite their high escape fractions, of minor importance for reionization due to inefficient star formation. We confirm previous claims that stars in atomic-cooling haloes with masses $10^8 M_{\odot} \lesssim M_{\text{halo}} \lesssim 10^{11} M_{\odot}$ are likely to be the most important source of reionization.

**Key words:** galaxies: high-redshift – dark ages, reionization, first stars – early Universe.

## 1 INTRODUCTION

Observations of Lyman $\alpha$ opacities in the spectra of quasi-stellar objects (QSOs) at high redshift have shown unambiguously that the Universe becomes nearly transparent to Lyman continuum (LyC) photons ($\lambda \le 912$ Å) at $z \sim 6$ (Becker et al. 2001; Fan et al. 2001, 2006; McGreer, Mesinger & D’Odorico 2015). Several candidates are identified as a potential source of reionization, including dwarf galaxies (e.g. Couchman & Rees 1986; Madau, Haardt & Rees 1999), active galactic nuclei (e.g. Shapiro & Giroux 1987; Haiman & Loeb 1998), accretion shock (Dopita et al. 2011), globular clusters (Ricotti 2002; Katz & Ricotti 2013, 2014) and X-rays from accreting stellar-mass black holes (e.g. Madau et al. 2004; Ricotti & Ostriker 2004; Mirabel et al. 2011). Many studies agree

that the primary source of reionization is likely to be massive stars in dwarf galaxies (Haehnelt et al. 2001; Cowie, Barger & Trouille 2009; Fontanot et al. 2014; Madau & Fragos 2016, cf. Madau & Haardt 2015); however, the time-scale over which reionization occurred and the mass range of haloes that provided the majority of the ionizing photons are issues that remain unresolved (e.g. Bolton & Haehnelt 2007; Ahn et al. 2012; Kuhlen & Faucher-Giguère 2012; Wise et al. 2014).

The two critical ingredients for reionization are star formation and escaping LyC photons. The former describes how many LyC photons are available from massive stars, while the latter determines what fraction are actually used to ionize the intergalactic medium (IGM). Unsurprisingly, the prediction of both quantities is very challenging, as galaxy evolution involves highly non-linear processes, such as the interaction between the interstellar medium (ISM) and feedback from stars. For this reason, numerical studies often report discrepant results on the escape fraction. An early attempt by

<footnote n="1">¹Kavli Institute for Cosmology and Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK</footnote>

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

<footnote n="3">³Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="4">⁴Observatoire de Lyon, UMR 5574, 9 avenue Charles Andre, F-69561 Saint Genis Laval, France</footnote>

<footnote>* E-mail: tkimm@ast.cam.ac.uk</footnote>
Gnedin, Kravtsov & Chen (2008) suggested that the escape frac-
tion roughly increases with halo mass in the range 10$^{10}$–10$^{12}$ M$_{\odot}$,
because stellar discs in lower mass haloes tend to be embedded in
gaseous discs in their simulations. Wise & Cen (2009) also found
a positive correlation between the escape fraction and halo mass
in the range 10$^{6}$–10$^{10}$ M$_{\odot}$, but as Wise et al. (2014) pointed out,
this may have been caused by the initial strong starburst due to the
absence of cooling while constructing the initial conditions of the
simulations. On the contrary, by post-processing hydrodynamics
simulations with strong stellar feedback, Razoumov & Sommer-
Larsen (2010) concluded that low-mass haloes (~10$^{8}$ M$_{\odot}$) show a
higher escape fraction of ~100 per cent, while less than 10 per cent
of LyC photons escape from the massive haloes with mass 10$^{11}$ M$_{\odot}$.
A negative dependence on halo mass in the atomic-cooling regime
is also found in other large cosmological simulations based on
post-processing (Yajima, Choi & Nagamine 2011; Paardekooper,
Khochfar & Dalla Vecchia 2013). Large-scale simulations of-
ten predict very high escape fractions of $f_{\text{esc}} \gtrsim 50$ per cent in
atomic-cooling haloes, but these conclusions may be subject to
numerical resolution (Ma et al. 2015; Paardekooper, Khochfar &
Dalla Vecchia 2015) and how the escape fraction is measured (i.e.
whether or not it is photon number-weighted; Kimm & Cen 2014).
Paardekooper et al. (2015) pointed out that significant absorption
of LyC photons occurs on GMC scales (i.e. 10 pc). Recent theo-
retical work based on effective feedback with high numerical res-
olution (~10 pc) suggests that, on average, only ~10 per cent of
LyC photons escape from their host haloes with the mass range
$10^8$ M$_{\odot} \lesssim M_{\text{halo}} \lesssim 10^{11}$ M$_{\odot}$ (Kimm & Cen 2014; Ma et al. 2015;
Xu et al. 2016). The only exception to this relatively low escape
fraction found in numerical simulations is mini-haloes where 40–
60 per cent of the ionizing photons escape the galaxies and con-
tribute to the ionization of the IGM (Wise et al. 2014; Xu et al. 2016).

In observations, the leakage of LyC photons is measured via the
relative flux density ratio ($F_{UV}/F_{LyC}$) between the ionizing part
of the spectrum at 900 Å and the non-ionizing part at 1500 Å
(Steidel, Pettini & Adelberger 2001). Once absorption due to the
IGM is corrected (e.g. Inoue et al. 2014), one can estimate the
absolute escape fraction assuming a ratio of the intrinsic lumi-
nosity at 900 and 1500 Å appropriate for the observed multiband
photometric data. The detection of LyC photons in the local Uni-
verse is limited to starburst galaxies (Leitet et al. 2011, 2013;
Borthakur et al. 2014; Leitherer et al. 2016), where generally small
escape fractions (~5 per cent) are observed. Star-forming galax-
ies at $z \sim 1$ with LyC detection also show low escape fractions
of a few per cent (Siana et al. 2007, 2010; Bridge et al. 2010;
Rutkowski et al. 2016). Efficient LyC leakers ($f_{\text{esc}} \gtrsim 10$ per cent)
seem to be more common at higher redshift ($z \gtrsim 3$; e.g. Reddy
et al. 2016), but only a handful of cases are confirmed as ro-
bust detections that are not affected by contamination due to low-
redshift interlopers along the line of sight (Mostardi et al. 2015;
Leethochawalit et al. 2016; Shapley et al. 2016). The average, rel-
ative escape fraction in nearly all observations is found to be very
small, even at high redshift ($f_{\text{esc}}^{\text{rel}} \lesssim 2$ per cent; Vanzella et al. 2010;
Boutsia et al. 2011; Mostardi et al. 2015; Siana et al. 2015; Grazian
et al. 2016) and appears to be in tension with the $f_{\text{esc}} \sim 10$ per cent
needed to reconcile the observed luminosity of high-redshift galax-
ies with observational constraints on the evolution of the average
neutral hydrogen fraction. However, it should be noted that these es-
timates mostly focus on small galaxies of mass $M_{\text{star}} \le 10^8$ M$_{\odot}$ or
$M_{\text{UV}} \gtrsim -18$ at $z \ge 6$, whereas observed samples are biased towards
bright galaxies ($M_{\text{UV}} \lesssim -20$; e.g. Grazian et al. 2016). Because star
formation in small galaxies is more bursty than in bright galaxies

observed at lower redshift (e.g. Speagle et al. 2014), it is conceiv-
able that the star-forming clouds are disrupted more efficiently in
simulated galaxies, resulting in higher escape fractions (e.g. Kimm
& Cen 2014; Cen & Kimm 2015). Moreover, since the simulated
galaxies are more metal-poor than the observed bright galaxies, they
are likely less affected by dust compared to observed galaxies (e.g.
Izotov et al. 2016). Finally, as pointed out by Cen & Kimm (2015),
individual measurements of the escape fraction may underesti-
mate the 3D escape fraction, especially when the escape fraction is
small.

Unlike the observed LyC flux that conveys information about
the instantaneous escape fraction, the Thompson electron opti-
cal depth ($\tau_e$), derived from the polarization signal of cosmic
microwave background (CMB) photons, provides a useful mea-
sure of how extended reionization was in the early Universe. The
analysis of the nine-year Wilkinson Microwave Anisotropy Probe
(WMAP9) observations suggested a high electron optical depth of
$\tau_e = 0.089 \pm 0.014$ (Hinshaw et al. 2013), indicating that ionized
hydrogen (H II) bubbles are likely to have grown relatively early.
However, the observed number density of bright galaxies in the ul-
traviolet (UV, $M_{\text{UV}} \lesssim -17$) is unable to explain such a high $\tau_e$ (e.g.
Bunker et al. 2010; Finkelstein et al. 2010; Bouwens et al. 2012).
By taking a parametric form of the UV luminosity density, moti-
vated by observations of the Hubble Ultra Deep Field, Robertson
et al. (2013) showed that the inclusion of small dwarf galaxies with
$-17 \le M_{\text{UV}} \le -13$ can increase $\tau_e$ to a higher value of 0.07,
provided that 20 per cent of LyC photons escape from the dark
matter haloes. Wise et al. (2014) claim that mini-haloes of mass
$M_{\text{halo}} \le 10^8$ M$_{\odot}$, corresponding to $M_{\text{UV}} \gtrsim -13$, may be able to
provide a large number of LyC photons to the IGM as LyC photons
escape freely from their host halo. Because the mini-haloes emerge
first and they are abundant in the early Universe ($z \ge 15$), the au-
thors find that the resulting $\tau_e \approx 0.09$ can easily accommodate the
WMAP9 analysis, demonstrating the potential importance of mini-
haloes to reionization of the Universe (see also Ahn et al. 2012).

However, a more accurate modelling of dust emission in our
Galaxy (Planck Collaboration XV 2014) and the use of the low
frequency instrument on the Planck Satellite lead to a decrease in
the optical depth to $\tau_e = 0.066 \pm 0.016$ (Planck Collaboration
XIII 2016a). The latest results utilizing the high-frequency instru-
ment to measure the low-multipole polarization signal point to a
possibility of an even lower value of $\tau_e = 0.055 \pm 0.009$ (Planck
Collaboration XLVI 2016b). Furthermore, recent findings of signif-
icant Ly$\alpha$ opacity fluctuations on large scales in absorption spectra
of $5 \lesssim z \lesssim 6$ QSOs (Becker et al. 2015) and the observed rapid evo-
lution of Ly$\alpha$ emitters at $z > 6$ (Ono et al. 2012; Caruana et al. 2014;
Pentericci et al. 2014; Schenker et al. 2014; Tilvi et al. 2014; Matthee
et al. 2015) suggest that reionization may have ended later than pre-
viously thought (e.g. Chardin et al. 2015; Choudhury et al. 2015;
Davies & Furlanetto 2016, cf. Haardt & Madau 2012; Mesinger
et al. 2015). If the contribution from mini-haloes were important
for reionization, this may potentially be in tension with the reduced
$\tau_e$ measurement and the long Ly$\alpha$ troughs still observed at $z \sim 5.6$
(Becker et al. 2015). Therefore, in this study, we revisit the im-
portance of mini-haloes and assess their role in reionization using
state-of-the-art numerical simulations.

This paper is organized as follows. In Section 2, we describe
the physical ingredients used in our cosmological, radiation-
hydrodynamic simulations. The measurements of the escape frac-
tion and star formation in the simulations are presented in Section 3.
Section 4 discusses the mechanisms responsible for the escape
of LyC photons and whether or not the ionizing radiation from
mini-haloes is crucial to reionization of the Universe. We summarize our findings in Section 5.

## 2 SIMULATION

We use RAMSES-RT, a radiation hydrodynamics code with adaptive mesh refinement (Teyssier 2002; Rosdahl et al. 2013; Rosdahl & Teyssier 2015), to study reionization due to starlight in mini-haloes with $10^6 \lesssim M_{\text{halo}}/M_{\odot} \lesssim 10^8$. The cosmological initial conditions are generated using MUSIC (Hahn & Abel 2011), with the cosmological parameters ($\Omega_m = 0.288$, $\Omega_\Lambda = 0.712$, $\Omega_b = 0.045$, $H_0 = 69.33 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $n_s = 0.971$ and $\sigma_8 = 0.830$) consistent with the WMAP9 results (Hinshaw et al. 2013). We first run dark matter only simulations with volume $(2 \text{ Mpc}/h)^3$, and identify nine regions hosting a halo of mass $\approx 10^8 M_{\odot}$ at $7 \le z \le 11$. The initial conditions for the zoom-in regions are then generated with a higher dark matter resolution of 90 $M_{\odot}$ to resolve each halo with more than 10 000 dark matter particles. We ensure that the haloes are not contaminated by coarse dark matter particles.

We solve the Euler equations using an HLLC scheme (Toro, Spruce & Speares 1994), with the typical courant number of 0.7. The Poisson equation is solved using a multigrid method (Guillet & Teyssier 2011). For the transport of multiple photon groups, RAMSES-RT uses a moment-based method with M1 closure for the Eddington tensor (Rosdahl et al. 2013; Rosdahl & Teyssier 2015, see also Aubert & Teyssier 2008). We adopt a GLF scheme to solve the advection of the photon fluids. Because the hydrodynamics is fully coupled to the radiation, the computational time step is usually determined by the speed of light. Since we are interested in the escaping flux, which is a conserved quantity, we use a reduced speed of light approximation ($\tilde{c} = 3 \times 10^{-3}c$) to keep the computational cost low, where c is the full speed of light.

Each zoom-in simulation is covered with $128^3$ root cells, and we allow for further refinement of the computational grid to achieve a maximum physical resolution of 0.7 pc. To do so, we adopt two different refinement criteria. First, a cell is refined if the total baryonic plus dark matter inside each cell exceeds eight times the mass of a dark matter particle (i.e. 720 $M_{\odot}$). Secondly, we enforce that the thermal Jeans length is resolved by at least 32 cells until it reaches the maximum resolution. Although the use of the latter condition is computationally expensive, it makes our simulations more robust than previous simulations, as the turbulent properties of gas can be more accurately captured (Federrath et al. 2011; Turk et al. 2012; Meece, Smith & O'Shea 2014).

We identify dark matter haloes with the AdaptaHop algorithm (Aubert, Pichon & Colombi 2004; Tweed et al. 2009). The centre of a dark matter halo is chosen as the centre of mass of the star particles in the halo. If a halo is devoid of stars, we use the densest location of the halo. The virial mass and radius of a halo is computed such that the mean density within the virial sphere is equivalent to $\Delta_{\text{crit}}\rho_{\text{crit}}$, where $\rho_{\text{crit}}$ is the critical density of the universe ($3H(z)^2/8\pi G$), $\Delta_{\text{crit}} = 18\pi^2 + 82x - 39x^2$ is the virial overdensity (Bryan & Norman 1998), $x \equiv \Omega_m/(\Omega_m + a^3\Omega_\Lambda) - 1$, G is the gravitational constant and $H(z)$ is the Hubble constant at some redshift $z$.

### 2.1 Star formation

Star formation is modelled based on a Schmidt law (Schmidt 1959),

$$ \frac{d\rho_{\text{star}}}{dt} = \epsilon_{\text{ff}} \frac{\rho_{\text{gas}}}{t_{\text{ff}}}, \qquad (1) $$

where $\rho_{\text{gas}}$ is the density of gas and $t_{\text{ff}} = \sqrt{3\pi/32 G \rho_{\text{gas}}}$ is the free-fall time. The main parameter characterizing star formation is the star formation efficiency per free-fall time ($\epsilon_{\text{ff}}$). Local observations find that the efficiency is only a few per cent when averaged over galactic scales (e.g. Kennicutt 1998). However, recent findings from small-scale numerical simulations suggest that $\epsilon_{\text{ff}}$ depends on physical properties of the ISM (Padoan & Nordlund 2011; Federrath & Klessen 2012). Motivated by this, we adopt a thermo-turbulent star formation model in which $\epsilon_{\text{ff}}$ is determined on a cell-by-cell basis (Devriendt et al., in preparation). The details of the model will be presented elsewhere, and here we briefly describe the basic idea for completeness.

The most fundamental assumption in the thermo-turbulent model is that the probability distribution function (PDF) of the density of a star-forming cloud is well described by a lognormal distribution. By integrating the gas mass from some critical density above which gas can collapse ($\rho = \rho(s_{\text{crit}})$) to infinity ($\rho = \infty$) per individual free-fall time, $\epsilon_{\text{ff}}$ can be estimated as (e.g. Federrath & Klessen 2012)

$$ \epsilon_{\text{ff}} = \frac{\epsilon_{\text{ecc}}}{2\phi_t} \exp\left(\frac{3}{8}\sigma_s^2\right) \left[1 + \text{erf}\left(\frac{\sigma_s^2 - s_{\text{crit}}}{\sqrt{2\sigma_s^2}}\right)\right], \qquad (2) $$

where $\sigma_s^2 = \ln(1 + b^2 M^2)$ is the standard deviation of the logarithmic density contrast ($s \equiv \ln(\rho/\rho_0)$), $\rho_0$ is the mean density of gas, b is a parameter that depends on the mode of turbulence driving, $\epsilon_{\text{ecc}} \approx 0.5$ is the maximum fraction of gas that can accrete on to stars without being blown away by proto-stellar jets and outflows, $\phi_t \approx 0.57$ is a factor that accounts for the uncertainties in the estimation of a free-fall time of individual clouds and $M$ is the sonic Mach number. We assume a mixture of solenoidal and compressive modes for turbulence ($b \approx 0.4$). An important quantity in equation (2) is the critical density ($s_{\text{crit}}$) that may be regarded as the minimum density above which gas in the post-shock regions of a cloud is magnetically supercritical and thus can collapse (Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund 2011). Numerical simulations suggest that $s_{\text{crit}}$ may be approximated as (Padoan & Nordlund 2011; Federrath & Klessen 2012)

$$ s_{\text{crit}} = \ln(0.067\theta^{-2}\alpha_{\text{vir}}M^2), \qquad (3) $$

where $\theta$ is a numerical factor of order unity that encapsulates the uncertainty in the post-shock thickness with respect to the cloud size, and we adopt $\theta = 0.33$ that gives a best fit to the results of Federrath & Klessen (2012). Here $\alpha_{\text{vir}} \equiv 2E_{\text{kin}}/|E_{\text{grav}}|$ is the virial parameter of a cloud, which we take to be $\alpha_0 \approx 5(\sigma_{\text{gas}}^2 + c_s^2)/(\pi\rho_{\text{gas}}G\Delta x^2)$, where $\Delta x$ is the size of a computational cell, $c_s$ is the gas sound speed and $\sigma_{\text{gas}}$ is the turbulent gas velocity that is the tracer of the velocity gradient. Note that the resulting $\epsilon_{\text{ff}}$ can be larger than 1 (fig. 1 of Federrath & Klessen 2012) if the sonic Mach number is very high ($M \gtrsim 10$) in tightly gravitationally bound regions ($\alpha_{\text{vir}} \lesssim 0.1$). In practice, we find that clouds with such conditions are extremely rare in our simulations, and $\epsilon_{\text{ff}}$ typically ranges from 5 per cent to 20 per cent when star particles are created.

This thermo-turbulent model allows for star formation only if the thermal plus turbulent pressure is not strong enough to prevent the gravitational collapse of a gas cloud. This may be characterized by the turbulent Jeans length (Bonazzola et al. 1987; Federrath & Klessen 2012)

$$ \lambda_{\text{J,turb}} = \frac{\pi \sigma_{\text{gas}}^2 + \sqrt{36\pi c_s^2 G \Delta x^2 \rho_{\text{gas}} + \pi^2 \sigma_{\text{gas}}^4}}{6 G \rho_{\text{gas}} \Delta x}. \qquad (4) $$
In order for gas to be gravitationally unstable, the Jeans length needs to be smaller than the size of a computational cell ($\lambda_{J, turb} \le \Delta$). Note that star formation occurs only in the maximally refined cells, because our refinement strategy enforces the thermal Jeans length to be resolved by 32 cells until it reaches the maximum level of refinement.

Once a potential site for star formation is identified, we estimate $\epsilon_{ff}$ using equation (2) and determine the number of newly formed stars ($N_*$) based on a Poisson distribution

$$P(N_*) = \frac{\lambda^{N_*}}{N_*!} \exp(-\lambda), \qquad (5)$$

with a mean of

$$\lambda = \epsilon_{ff} \frac{\rho_{\text{gas}} \Delta x^3}{m_{\star, \text{min}}} \left( \frac{\Delta t_{\text{sim}}}{t_{\text{ff}}} \right), \qquad (6)$$

where $m_{\star, min}$ is the minimum mass of a star particle and $\Delta t_{sim}$ is the simulation integration time step. We adopt $m_{\star, min} = 91 M_{\odot}$ for Pop II stars,<fnref n="1" /> which would host a single supernova (SN) for a Kroupa initial mass function (IMF; Kroupa 2001). Of the stellar mass, 21 per cent is returned to the surrounding medium as a result of SN explosions, and 1 per cent is assumed to be newly synthesized to metals (i.e. a metal yield of 0.01).

The formation of Pop III stars is included following Wise et al. (2012b). We adopt a Salpeter-like IMF for masses above the characteristic mass ($M_{\text{char}}$), while the formation of low-mass Pop III stars is assumed to be inefficient,

$$\frac{dN}{d \log M} \propto M^{-1.3} \exp \left[ - \left( \frac{M_{\text{char}}}{M} \right)^{1.6} \right], \qquad (7)$$

where N is the number of Pop III stars per logarithmic mass bin. The precise determination of the characteristic mass is a matter of debate. Early studies suggested that the mass of protostellar clumps is ~100 $M_{\odot}$ (Abel, Bryan & Norman 2002; Bromm, Coppi & Larson 2002; Yoshida et al. 2006). Later, several groups point out that gas clumps may be fragmented further reducing the characteristic mass to ~40 $M_{\odot}$ (Turk, Abel & O'Shea 2009; Greif et al. 2012). However, recent radiation-hydrodynamics simulations report that several tens to a thousand solar masses of gas may collapse to form a Pop III star (Hirano et al. 2014; Hosokawa et al. 2016, cf. Lee & Yoon 2016). In this work, we adopt $M_{\text{char}} = 100 M_{\odot}$, consistent with the most recent simulations.

We assume that Pop III stars form only in a region where the gas metallicity is below $10^{-6} Z_{\odot}$. This means that at least one Pop III star will form in a dark matter halo during the initial gas collapse due to radiative cooling by molecular hydrogen. In principle, the external pollution by neighbouring haloes can suppress the formation of Pop III stars (e.g. Smith et al. 2015), but our simulated haloes are chosen to reside in an isolated environment and thus are not affected by neighbours. Note that more than one Pop III star can form in each halo if the first Pop III star does not explode and enrich the IGM/ISM or if a pristine gas cloud is accreted on to a dark matter halo through halo mergers.

## 2.2 Stellar feedback

Modelling the feedback from stars is essential to predict the escape of LyC photons in dwarf galaxies. In order for the LyC photons to

leave their host dark matter halo, feedback should clear away low-density channels or entirely blow out the birth clouds. Otherwise, the photons will simply be absorbed by neutral hydrogen inside of the halo. We include three different types of feedback (photoionization, radiation pressure from the absorption of UV and infrared, IR, photons, and Type II SN feedback) in our simulations.

### 2.2.1 Radiation feedback

Young, massive stars emit large amounts of ionizing photons that drive winds through various processes. Because the absorption cross-section of neutral hydrogen is so large ($\sigma_{\text{abs}} \sim 6 \times 10^{-18} \text{ cm}^2$), the presence of a small amount of hydrogen makes the ISM optically thick to photons with $E > 13.6 \text{ eV}$. When the ISM is fully ionized, dust becomes the next most efficient absorber, as its opacity in the UV wavelengths is large as well ($\kappa_{\text{abs}} \sim 1000 \text{ cm}^2 \text{ g}^{-1}$).

Of the several radiation feedback processes, photoionization is probably the most important mechanism that governs the dynamics of a giant molecular cloud (GMC; Dale et al. 2014; Lopez et al. 2014; Rosdahl & Teyssier 2015). LyC photons can ionize hydrogen, which heat the gas to $T \approx 2 \times 10^4 \text{ K}$. This creates an overpressurized H II bubble that lowers the density of the ambient medium and drives winds with velocities up to $10 \text{ km s}^{-1}$ (e.g. Krumholz, Stone & Gardiner 2007; Walch et al. 2012; Dale et al. 2014). To capture the dynamics of the H II region, the Stromgren sphere radius ($r_S$) should be resolved.

$$r_S = \left( \frac{3 \dot{N}_p}{4 \pi \alpha_B n_H^2} \right)^{1/3} \\ \approx 1.2 \text{ pc} \left( \frac{m_{\text{star}}}{10^3 M_{\odot}} \right)^{1/3} \left( \frac{n_H}{10^3 \text{ cm}^{-3}} \right)^{-2/3}, \qquad (8)$$

where $\dot{N}_p$ is the production rate of ionizing photons and $\alpha_B = 2.6 \times 10^{-13} \text{ cm}^3 \text{ s}^{-1}$ is the case B recombination rate coefficient at $T = 10^4 \text{ K}$. For the latter equality, we use $N_p = 5 \times 10^{46} \text{ s}^{-1}$ per $1 M_{\odot}$. Note that this scale describes the maximum distance within which recombination is balanced by ionization, assuming that photoheating does not affect the dynamics of the ISM.

When ionizing photons are absorbed by the neutral ISM, their momentum is transferred to the medium (e.g. Haehnelt 1995) at a rate

$$\dot{p}_v = \sum_{\text{groups}} \frac{F_i}{c} \left( \kappa_i + \sum_j \sigma_{ij} n_j \right), \qquad (9)$$

where $F_i$ is the photon flux (erg cm$^{-2}$s$^{-1}$) for the ith photon group, $\kappa$ is the dust opacity (cm$^2$ g$^{-1}$), $\sigma$ is the photoionization cross-section (cm$^2$) and $n_j$ is the number density (cm$^{-3}$) of the ion species j. For a Kroupa IMF, direct radiation pressure from ionizing radiation ($\lambda \le 912 \AA$) can impart momentum of up to ~40 km s$^{-1}$ M$_{\odot}$ if we integrate the number of ionizing photons from a simple stellar population of 1 M$_{\odot}$ until 50 Myr (Leitherer et al. 1999). The absorption of non-ionizing UV and optical photons by dust can further increase the total momentum input to ~190 km s$^{-1}$ M$_{\odot}$ within 50 Myr.

We also take into account radiation pressure by trapped IR photons. IR photons are generated when UV and optical photons are absorbed by dust or when molecular hydrogen is fluorescently excited by the absorption of Lyman-Werner photons and radiatively de-excited through forbidden rotational-vibrational transitions (see the Chemistry section). We assume that these IR photons are efficiently trapped only if the optical depth over the cell width by

<footnote n="1">¹ We discuss the possible impact of the IMF sampling, the neglect of runaway stars and the uncertainties in stellar evolutionary models on the determination of the escape fractions in Section 4.3.</footnote>
dust is high. The resulting trapped IR photon energy in each cell is modelled as (Rosdahl & Teyssier 2015)

$$E_{\text{IR,trapped}} = f_{\text{trapped}} E_{\text{IR}} = \exp\left(-\frac{2}{3\tau_{\text{d}}}\right) E_{\text{IR}}, \quad (10)$$

where $\tau_{\text{d}} = \kappa_{\text{sc}}\rho_{\text{gas}}\Delta x$ and $\kappa_{\text{sc}} \sim 5(Z_{\text{gas}}/Z_{\odot})\text{cm}^2\text{g}^{-1}$ is the scattering cross-section by dust (Semenov et al. 2003). This trapped IR radiation is then included as a non-thermal pressure term in the momentum equation, as

$$\frac{\partial \rho \mathbf{v}}{\partial t} + \nabla \cdot (\rho \mathbf{v} \otimes \mathbf{v} + (P + P_{\text{rad}})\mathbf{I}) = \dot{\mathbf{p}}_{\gamma} + \rho \nabla \Phi \quad (11)$$

where

$$P_{\text{rad}} = \frac{\tilde{c}}{c} \frac{E_{\text{IR,trapped}}}{3}. \quad (12)$$

The non-thermal pressure imparted by trapped IR radiation ($P_{\text{rad}}$) is also added to the energy equation. Note that these trapped IR photons are advected with the gas, while the remaining fraction of the IR energy density

$$E_{\text{IR,stream}} = \left[ 1 - \exp\left(-\frac{2}{3\tau_{\text{d}}}\right) \right] E_{\text{IR}} \quad (13)$$

is diffused out to the neighbouring cells and it is re-evaluated whether or not these photons are trapped by dust (Rosdahl & Teyssier 2015).

In this paper, we adopt the photon production rates of Pop II stars from Bruzual & Charlot (2003) assuming a Kroupa IMF. This is done by interpolating the spectral energy distributions for a given metallicity and age, and by counting the LyC photons from the spectrum of each star particle. The lifetime and the photon production rates for Pop III stars are taken by fitting the results of Schaerer (2002).

### 2.2.2 Type II SN explosions

We adopt the mechanical feedback scheme introduced by Kimm & Cen (2014) and Kimm et al. (2015) to model the explosion of massive ($M \ge 8 M_{\odot}$) Pop II stars. Based on Thornton et al. (1998, see also Blondin et al. 1998; Geen et al. 2015; Kim & Ostriker 2015; Martizzi, Faucher-Giguère & Quataert 2015), this model captures the correct radial momentum input from SN explosions at the snowplough phase ($p_{\text{SN, snow}}$),

$$p_{\text{SN,snow}} = 3 \times 10^5 \text{ km s}^{-1} \text{M}_{\odot} n_{\text{H}}^{-2/17} E_{51}^{16/17} Z'^{-0.14}, \quad (14)$$

by imparting momentum according to the stage of the Sedov–Taylor blast wave. Here $n_{\text{H}}$ is the hydrogen number density in units of cm$^{-3}$, $E_{51}$ is the explosion energy in units of $10^{51}$ erg and $Z' = \max[0.01, Z/0.02]$ is the metallicity of gas, normalized to the solar value.

Recently, Geen et al. (2015) have shown that the final radial momentum from an SN can be augmented by including photoionization from massive stars. Because the thermal energy that is liberated during the ionization process overpressurizes and decreases the density of the surroundings into which SNe explode, the final momentum from SNe in a medium pre-processed by ionizing radiation is found to be significantly larger than without it. Most notably, they find that the amount of momentum is nearly independent of the background density, indicating that more radial momentum should be imparted to the ISM than suggested by Thornton et al. (1998), especially when SNe explode in dense environments. In principle, this extra momentum should be generated by solving the full radiation hydrodynamics, but this requires the simulation to resolve the Stromgren radius. If young stars are embedded in a cloud denser than $10^3 \text{ cm}^{-3}$,

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Final radial momentum from an individual SN explosion in different environments ($p_{\text{SN}}$). Grey triangles represent the results from 1D hydrodynamic calculations (Thornton et al. 1998), while violet circles show the final momentum in the presence of ionizing radiation from 3D radiation-hydrodynamic simulations (Geen et al. 2015). The grey and violent dashed lines display a simple fit to these results ($p_{\text{SN}} \propto n_{\text{H}}^{-2/17}$, $p_{\text{SN}} \sim \text{const}$), respectively. The inclusion of photoheating enhances $p_{\text{SN}}$ by reducing the density of the ambient medium before SNe explode and also by imparting momentum from photoheating and direct radiation pressure. We take this effect into consideration explicitly only if the Stromgren sphere is under-resolved in our simulations. When the Stromgren sphere is resolved, the effect is captured self-consistently in our simulations. The solid lines with different colours illustrate three examples of the momentum that we would impart during the momentum-conserving phase at different resolutions in the presence of the radiation from a star cluster with $M_{\text{star}} = 10^3 \text{ M}_{\odot}$, and the dotted lines indicate the corresponding densities at which the Stromgren sphere becomes unresolved. Note that we inject momentum of the initial ejecta ($p_{\text{SN,ad}} = 4.5 \times 10^4 \text{ M}_{\odot} \text{ km s}^{-1}$) if the mass in the neighbouring cells is negligible compared to the ejecta mass (see equation 21).</figcaption></figure>

the effect of photoionization is likely to be underestimated with our parsec-scale resolution. In order to circumvent this issue, we adopt a simple fit to the results of Geen et al. (2015),

$$p_{\text{SN+PH}} = 4.2 \times 10^5 \text{ km s}^{-1} \text{M}_{\odot} E_{51}^{16/17} Z'^{-0.14}, \quad (15)$$

if the Stromgren sphere is underresolved ($\Delta x \gg r_{\text{S}}$). We adopt the dependence on the SN energy and the metallicity from Thornton et al. (1998). Note that the values taken from Geen et al. (2015) are lowered by a factor of $1.2^{16/17}$, as they use a 20 per cent larger SN explosion energy ($1.2 \times 10^{51}$ erg) compared to other studies (Thornton et al. 1998). The final momentum input during the snowplough phase is then taken from the combination of $p_{\text{SN, snow}}$ and $p_{\text{SN + PH}}$ by comparing the resolution of a computational cell ($\Delta x$) with the Stromgren radius ($r_{\text{S}}$), as

$$p_{\text{SN}} = p_{\text{SN,snow}} \exp\left(-\frac{\Delta x}{r_{\text{S}}}\right) + p_{\text{SN+PH}} \left(1 - \exp\left[-\frac{\Delta x}{r_{\text{S}}}\right]\right). \quad (16)$$

Fig. 1 illustrates three examples of the SN momentum that we would inject at different resolutions (1, 10 and 100 pc) as a function of
<table><caption>**Table 1.** Summary of simulation parameters.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td>L<sub>box</sub></td><td>2 Mpc h<sup>-1</sup></td><td>Simulation box size</td></tr><tr><td>Δx<sub>min</sub></td><td>0.7 pc</td><td>Physical size of the finest cell</td></tr><tr><td>m<sub>DM</sub></td><td>90 M<sub>☉</sub></td><td>Mass resolution of DM particles</td></tr><tr><td>m<sub>star</sub></td><td>91 M<sub>☉</sub></td><td>Mass resolution of Pop II star particles</td></tr><tr><td>λ<sub>J</sub>/Δx</td><td>32</td><td>Jeans length criterion for refinement</td></tr><tr><td>ε<sub>ff</sub></td><td>FK12</td><td>SF efficiency per free-fall time</td></tr><tr><td>N<sub>halo</sub></td><td>9</td><td>Number of zoom-in haloes</td></tr><tr><td>IMF</td><td>Kroupa (2001)</td><td>(i.e. 1 SN per 91 M<sub>☉</sub>)</td></tr></tbody></table>

density in the presence of the radiation from a star cluster with $10^3$ M☉.

Specifically, the model first calculates the mass ratio ($\chi$) between the swept-up mass ($M_{swept}$) and the ejecta mass ($M_{ej}$) along each $N_{nbor}$ neighbouring cell, as

$$ \chi \equiv dM_{\text{swept}}/dM_{\text{ej}}, \qquad (17) $$

where

$$ dM_{\text{ej}} = (1 - \beta_{\text{sn}})M_{\text{ej}}/N_{\text{nbor}}, \quad \text{and} \qquad (18) $$

$$ dM_{\text{swept}} = \rho_{\text{nbor}} \left(\frac{\Delta x}{2}\right)^3 + \frac{(1 - \beta_{\text{sn}})\rho_{\text{host}} \Delta x^3}{N_{\text{nbor}}} + dM_{\text{ej}}. \qquad (19) $$

Here $\beta_{sn}$ is a parameter that determines what fraction of the gas mass ($M_{ej} + \rho_{host} \Delta x^3$) is re-distributed to the host cell of an SN. In order to distribute the mass evenly to the host and neighbouring cells in the uniformly refined case, we take $\beta_{sn} = 4/52$. Note that since the maximum number of neighbouring cells is 48 if they are more refined than the host cell of an SN (see fig. 15 of Kimm & Cen 2014), we use $N_{nbor} = 48$. If the neighbours are not further refined, we simply take the physical properties ($\rho$, $v$, $Z$) of the neighbours assuming that they are refined.

We then use the mass ratio ($\chi$) to determine the phase of the Sedov-Blast wave. To do so, we define the transition mass ratio ($\chi_{tr}$) by equating $p_{SN}$ with the radial momentum one would expect during the adiabatic phase $p_{ad} = \sqrt{2\chi M_{ej} f_e E_{SN}}$, where $f_e \sim 2/3$ is the fraction of energy that is left in the beginning of the snowplough phase, as

$$ \chi_{tr} \approx \frac{900 n_{\text{H}}^{-4/17} E_{51}^{-2/17} Z'^{-0.28}}{f_e(M_{\text{ej}}/M_{\odot})}. \qquad (20) $$

If $\chi$ is greater than $\chi_{tr}$, we inject the momentum during the snowplough phase, whereas the momentum during the adiabatic phase is added to the neighbouring cell, as

$$ p_{\text{SN}} = \begin{cases} p_{\text{SN,ad}} = \sqrt{2\chi M_{\text{ej}} f_e(\chi) E_{\text{SN}}} & (\chi < \chi_{\text{tr}}) \\ p_{\text{SN}} & (\chi \ge \chi_{\text{tr}}) \end{cases}, \qquad (21) $$

where the fraction of energy left in the SN bubble ($f_e(\chi) \equiv 1 - \frac{\chi^{-1}}{3(\chi_{tr}-1)}$) is modified to smoothly connect the two regimes.

In order to account for the fact that the lifetime of SN progenitors varies from 3 to 40 Myr depending on their mass, we randomly draw the lifetime based on the integrated SN occurrence rate from STARBURST99 (Leitherer et al. 1999) using the inverse method, as in the MFBmp model from Kimm et al. (2015). Simulation parameters are summarized in Table 1.

### 2.2.3 Explosion of Pop III stars

The explosions of Pop III stars are modelled similarly as Pop II explosions, but with different energy and metal production rates.

<table><caption>**Table 2.** Properties of photon groups.</caption><thead><tr><th>Photon group</th><th>ε<sub>0</sub> (eV)</th><th>ε<sub>1</sub> (eV)</th><th>κ (cm<sup>2</sup> g<sup>-1</sup>)</th><th>Main function</th></tr></thead><tbody><tr><td>IR</td><td>0.1</td><td>1.0</td><td>5</td><td>Radiation pressure (RP)</td></tr><tr><td>Optical</td><td>1.0</td><td>5.6</td><td>10<sup>3</sup></td><td>Direct RP</td></tr><tr><td>FUV</td><td>5.6</td><td>11.2</td><td>10<sup>3</sup></td><td>Photoelectric heating</td></tr><tr><td>LW</td><td>11.2</td><td>13.6</td><td>10<sup>3</sup></td><td>H<sub>2</sub> dissociation</td></tr><tr><td>EUV<sub>H1,1</sub></td><td>13.6</td><td>15.2</td><td>10<sup>3</sup></td><td>H<sub>1</sub> ionization</td></tr><tr><td>EUV<sub>H1,2</sub></td><td>15.2</td><td>24.59</td><td>10<sup>3</sup></td><td>H<sub>1</sub> and H<sub>2</sub> ionization</td></tr><tr><td>EUV<sub>He1</sub></td><td>24.59</td><td>54.42</td><td>10<sup>3</sup></td><td>He<sub>1</sub> ionization</td></tr><tr><td>EUV<sub>He2</sub></td><td>54.42</td><td>∞</td><td>10<sup>3</sup></td><td>He<sub>2</sub> ionization</td></tr></tbody></table>

Stars with $40 M_{\odot} \le M_{\star} \le 120 M_{\odot}$ and $M_{\star} \ge 260 M_{\odot}$ are likely to implose without releasing energy and metals, while massive stars with $120 M_{\odot} < M_{\star} < 260 M_{\odot}$ may end up as a pair-instability SN (Heger et al. 2003). The explosions of stars less massive than $40 M_{\odot}$ are modelled as either normal Type II SN with $10^{51}$ erg if $11 M_{\odot} \le M_{\star} \le 20 M_{\odot}$ (Woosley & Weaver 1995) or hypernova if $20 M_{\odot} \le M_{\star} \le 40 M_{\odot}$ (Nomoto et al. 2006). For the explosion energy ($E_{SN, III}$) and returned metal mass ($M_{\star}$), we adopt the compilation of Wise et al. (2012a), which is based on Woosley & Weaver (1995), Heger & Woosley (2002) and Nomoto et al. (2006), as

$$ \frac{E_{\text{SN,III}}}{10^{51} \text{ erg}} = \begin{cases} 1 & [11 \le M_{\star} < 20] \\ (-13.714 + 1.086 M) & [20 \le M_{\star} \le 40] \\ (5.0 + 1.304 \times (M_{\text{He}} - 64)) & [140 \le M_{\star} \le 230] \\ 0 & \text{elsewhere} \end{cases} \qquad (22) $$

where $M_{\text{He}} = \frac{13}{24}(M_{\star} - 20)$ is the helium core mass, and

$$ \frac{M_{\text{z}}}{M_{\odot}} = \begin{cases} 0.1077 + 0.3383(M_{\star} - 11) & [11 \le M_{\star} < 20] \\ -2.7650 + 0.2794 M_{\star} & [20 \le M_{\star} \le 40] \\ (13/24)(M_{\star} - 20) & [140 \le M_{\star} \le 230] \\ 0 & \text{elsewhere} \end{cases} \qquad (23) $$

We neglect accretion and feedback from black holes formed by the implosion of massive Pop III stars.

## 2.3 Non-equilibrium photochemistry and radiative cooling

The public version of RAMSES-RT can solve non-equilibrium chemistry of hydrogen and helium species (H I, H II, He I, He II, He III and e⁻), involving collisional excitation, ionization and photoionization (Rosdahl et al. 2013). In order to take into account cooling by molecular hydrogen (H₂), which is essential to model gas collapse in mini-haloes, we have made modifications based on Glover et al. (2010) and Baczynski, Glover & Klessen (2015). Note that the photon number density and fluxes in the eight energy bins, which we describe in Table 2, are computed in a self-consistent way by tracing the photon fluxes from each star particle. The chemical reactions and radiative cooling are fully coupled with the eight photon groups. More details of the photochemistry will be presented in Katz et al. (2016) in terms of the prediction of molecular hydrogen in high-z galaxies.

Molecular hydrogen mainly forms on the surface of interstellar dust grains. However, in the early universe where there is little dust (e.g. Fisher et al. 2014), the formation is dominated by the reaction involving H⁻. These hydrogen molecules are dissociated
by Lyman-Werner photons with energy $11.2 \text{ eV} \le E \le 13.6 \text{ eV}$ or collisions with other species (H₁, H₂, He₁ and e⁻). Furthermore, H₂ can be ionized by photons with $E \ge 15.2 \text{ eV}$. We assume that all $H_2^+$ are immediately destroyed by dissociative recombination. To estimate the formation rate of molecular hydrogen through the H⁻ channel, we assume that the abundance of H⁻ is set by the equilibrium between the formation and destruction via associative detachment and mutual neutralization (e.g. Anninos et al. 1997), as

$$k_1 n_{\text{H}}; n_e = k_2 n_{\text{H}} - n_{\text{H}_1} + k_5 n_{\text{H}} - n_{\text{H}_2} + k_{13} n_{\text{H}} - n_e, \quad (24)$$

where $k_X$ is the reaction rate (reactions 1, 2, 5 and 13 in appendix B of Glover et al. 2010). This neglects the photodetachment of H⁻ by IR photons, and it may thus lead to the overestimation of the H⁻ abundance when Pop II stars are present (Cen 2017). However, we expect that cooling is dominated by metals once Pop III stars explode (see fig. 2 of Wise et al. 2014, for example), hence the gas collapse in mini-haloes is unlikely to be significantly affected.

Radiative cooling by hydrogen and helium species is directly computed from the chemical network. In particular, we include the cooling by molecular hydrogen following Halle & Combes (2013), which is largely based on Hollenbach & McKee (1979). In addition, Lyman-Werner photons also heat the gas when they photodissociate and photoionize H₂ or when they indirectly excite the vibrational levels of H₂ (see section 2.2.4 in Baczynski et al. 2015). Gas can cool further with the aid of metals, which we consider by interpolating look-up tables that are pre-computed with the CLOUDY code (Ferland et al. 1998) as a function of density, temperature and redshift (Smith et al. 2011). Finally, we also include photoelectric heating on dust by UV photons with $5.6 \text{ eV} \le h\nu \le 13.6 \text{ eV}$ (Bakes & Tielens 1994) following Baczynski et al. (2015, Section 2.2.5), as

$$\mathcal{H}_{\text{pe}} = 1.3 \times 10^{-24} \epsilon G_0 f_{\text{D/G}} n_{\text{H}} [\text{erg cm}^{-3} \text{ s}^{-1}], \quad (25)$$

where $G_0$ is the strength of the local intensity in each cell, normalized to the Habing field ($1.6 \times 10^{-3} \text{ erg s}^{-1} \text{ cm}^{-2}$; Habing 1968), and $f_{\text{D/G}} = 1$ is the dust-to-gas ratio, normalized to the local ISM value (Draine et al. 2007). The efficiency for the heating ($\epsilon$) is taken from Wolfire et al. (2003), as

$$\epsilon = \frac{4.9 \times 10^{-2}}{1 + 4.0 \times 10^{-3} \left( G_0 T^{1/2} / n_e \phi_{\text{PAH}} \right)^{0.73}} \frac{3.7 \times 10^{-2} \left( T / 10^4 \right)^{0.7}}{1 + 2.0 \times 10^{-4} \left( G_0 T^{1/2} / n_e \phi_{\text{PAH}} \right)}, \quad (26)$$

where $\phi_{\text{PAH}} = 0.5$ is a factor that controls the collision rates between electron and polycyclic aromatic hydrocarbon (PAH). We do not use any uniform UV or Lyman-Werner background radiation (e.g. Haardt & Madau 2012).

## 3 RESULTS

The main aim of this paper is to assess the contribution of mini-haloes to the reionization history of the universe. For this purpose, we investigate the evolution of nine dwarf galaxies in haloes of mass $10^6 \text{ M}_\odot \lesssim M_{\text{vir}} \lesssim 10^8 \text{ M}_\odot$ at $7 \le z \le 20$. In this section, we first describe the main features of the simulated galaxies, present the escape fraction of LyC photons at the virial radius and discuss the physical processes governing the evolution of the escape fraction.

### 3.1 Galactic properties of the dwarf population during reionization

Our simulated haloes begin forming Pop III stars when the halo mass approaches a few times $10^6 \text{ M}_\odot$. These Pop III stars disperse dense gas clouds and pollute the ISM and IGM with metals via energetic explosions. We find that the typical metallicity of the halo gas after the explosion of Pop III stars is $\sim 10^{-3}-10^{-2} \text{ Z}_\odot$, consistent with previous studies (Greif et al. 2010; Ritter et al. 2012). The enrichment of the dense medium ($n_H \ge 100 \text{ cm}^{-3}$) takes place more slowly than for the IGM ($\sim 10^{-4}-10^{-3} \text{ Z}_\odot$), as this gas mixes with the newly accreted, pristine material with primordial composition. Once the metal-enriched gas collapses, Pop II stars form in a very stochastic fashion. As the haloes become massive ($\sim 10^8 \text{ M}_\odot$) and a large amount of gas accumulates in the halo centres, the star formation histories become less bursty, compared to those in haloes with masses of a few times $10^7 \text{ M}_\odot$. Because stellar feedback violently disrupts star-forming clouds, the gas component of these mini-haloes show irregular morphologies rather than well-defined discs (Fig. 2). The resulting stellar metallicities in haloes with $\sim 10^8 \text{ M}_\odot$ range from $Z = 10^{-3}$ to $10^{-1.7} \text{ Z}_\odot$ with a median of $Z = 10^{-2.6} \text{ Z}_\odot$. We summarize several galactic properties in the nine simulated haloes in Table 3.

Overall, we find that star formation is very inefficient in the mini-haloes (Fig. 5). For example, haloes with masses $\sim 10^{7.5} \text{ M}_\odot$ form clusters of stars of $\sim 10^{3-4} \text{ M}_\odot$. At larger masses ($M_{\text{vir}} \sim 10^8 \text{ M}_\odot$), the efficiency becomes higher, but still $\lesssim 1$ per cent of baryons are converted into stars. We note that these results are consistent with the recent radiation-hydrodynamic calculations of Wise et al. (2014) and Xu et al. (2016). However, our results are slightly different from these studies in the sense that the least massive haloes ($M_{\text{vir}} \lesssim 10^7 \text{ M}_\odot$) appear to host progressively smaller amount of stars, whereas, in Wise et al. (2014), the stellar mass appears to saturate at a few times $10^3 \text{ M}_\odot$. We also find that the dispersion in the stellar mass-halo mass relation is smaller than that of Xu et al. (2016). This is likely to be due to the small number of samples used in this work. Indeed, Xu et al. (2016) find a larger dispersion than Wise et al. (2014) when they increase the number of galaxies from 32 to $\sim 2000$ simulated with the same assumptions.

The inefficient star formation in our simulations is due to strong stellar feedback. This can be inferred from Fig. 3 where the star formation histories are shown to be very bursty. The typical time-scale of star formation does not exceed $\sim 10 \text{ Myr}$ and is often smaller than 5 Myr when the halo mass is small. The recovery time after a burst of star formation in the mini-haloes is also large, ranging from $\sim 20$ to $\sim 200 \text{ Myr}$. Because radiation can drive winds by overpressurizing the ISM and early SNe can create cavities in the star-forming regions, we find that the majority of SNe explode in a low-density environment with $n_H \sim 2 \times 10^{-3} \text{ cm}^{-3}$, while stars form only in very dense environments with $n_H \ge 4 \times 10^4 \text{ cm}^{-3}$ (Fig. 4).

At present, there are no direct observational constraints on the evolution of galaxies in mini-haloes at high redshift, and thus it is not currently possible to validate our feedback model. Nevertheless, it is worth noting that the simulated galaxies follow the stellar mass-to-halo mass sequence obtained from Kimm & Cen (2014) (Fig. 5), which successfully reproduces the faint-end slope and normalization of the observed UV luminosity function at $z \sim 7$ without dust correction. It is also interesting to note that the predicted stellar masses are comparable to the local stellar mass-to-halo mass relation derived from the abundance matching technique (Behroozi, Wechsler & Conroy 2013) when extrapolated to the mini-halo mass
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Composite images of the density, H II fraction, and H I-ionizing photon density for seven simulated haloes out of the total nine samples. The central panel displays the location of the nine zoom-in haloes in the entire simulation box of length 2 Mpc h⁻¹ (comoving) at z = 7. We evolve each halo until its mass becomes large enough (M_halo ≈ 10⁸ M☉) to radiate energy away mainly by atomic transitions. Actively star-forming regions are shown as bright yellow colours, while highly ionized regions are displayed as light blue colours. Dark regions show mostly neutral gas. It can be seen that the central star-forming clump is disrupted by stellar feedback in many cases, leading to a high escape fraction of LyC photons. The bottom rightmost panel shows the projected gas density distribution of the H6 halo for comparison.</figcaption></figure>

<table><caption>**Table 3.** Summary of simulation results. All quantities are measured at the final redshift of each simulation. Column (1): ID of simulated halo. Column (2): virial mass of the dark matter halo. Column (3): total stellar mass of Pop II. Column (4): total gas mass inside a halo. Column (5): mass-weighted mean metallicity of star particles. Column (6): mass-weighted mean gas metallicity inside a halo. Column (7): total stellar mass of Pop III. Column (8): initial redshift to form Pop III stars. Column (9): initial redshift to form Pop II stars. Column (10): final redshift of each simulation.</caption><thead><tr><th>Halo ID</th><th>log M<sub>halo</sub><br>(M☉)</th><th>log M<sub>star, II</sub><br>(M☉)</th><th>log M<sub>gas</sub><br>(M☉)</th><th>⟨log Z<sub>star</sub>⟩<br>(Z☉)</th><th>log Z<sub>gas</sub><br>(Z☉)</th><th>log M<sub>star, III</sub><br>(M☉)</th><th>z<sub>PopIII</sub></th><th>z<sub>PopII</sub></th><th>z<sub>final</sub></th></tr></thead><tbody><tr><td>H0</td><td>8.05</td><td>4.53</td><td>6.93</td><td>-3.0</td><td>-2.7</td><td>3.20</td><td>12.4</td><td>10.6</td><td>7.0</td></tr><tr><td>H1</td><td>8.03</td><td>5.17</td><td>7.10</td><td>-2.8</td><td>-2.4</td><td>1.57</td><td>14.3</td><td>12.7</td><td>7.3</td></tr><tr><td>H2</td><td>8.08</td><td>5.26</td><td>7.06</td><td>-2.5</td><td>-2.0</td><td>3.35</td><td>18.1</td><td>15.0</td><td>10.6</td></tr><tr><td>H3</td><td>8.00</td><td>5.08</td><td>7.16</td><td>-2.7</td><td>-2.5</td><td>1.40</td><td>10.9</td><td>10.0</td><td>7.5</td></tr><tr><td>H4</td><td>8.03</td><td>4.72</td><td>6.32</td><td>-2.8</td><td>-2.6</td><td>3.50</td><td>14.7</td><td>13.4</td><td>7.0</td></tr><tr><td>H5</td><td>8.01</td><td>5.17</td><td>6.98</td><td>-2.7</td><td>-2.3</td><td>3.31</td><td>14.7</td><td>13.6</td><td>8.1</td></tr><tr><td>H6</td><td>7.96</td><td>5.70</td><td>7.17</td><td>-2.4</td><td>-2.0</td><td>2.16</td><td>15.6</td><td>13.8</td><td>7.3</td></tr><tr><td>H7</td><td>7.91</td><td>4.61</td><td>7.09</td><td>-3.0</td><td>-2.9</td><td>2.69</td><td>13.4</td><td>12.4</td><td>11.2</td></tr><tr><td>H8</td><td>7.85</td><td>4.17</td><td>7.00</td><td>-3.3</td><td>-3.2</td><td>1.59</td><td>12.5</td><td>10.7</td><td>8.1</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 3. Star formation histories of the dwarf galaxies in our simulated mini-haloes. The x-axis indicates the age of stars measured at the end of each simulation. Different colour coding denotes different galaxies. We split the sample for clarity. It can be seen that star formation is very stochastic. The recovery time from the stellar feedback ranges from ~20 to 200 Myr.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4. Local environments of star formation and Type II SN explosions in our mini-halo simulations. Even though stars form in very dense media, SN explosions occur in low-density regions. This is first because radiation lowers the ambient density before SN explosions, and secondly, because the lifetime of massive stars ranges from 3 to 40 Myr and late explosions take advantage of the early SN events.</figcaption></figure>

regime. The fact that the simulated galaxies are very metal-poor (~0.003 Z⊙) suggests that our results do not suffer from the overheating problem, as it usually leads to significantly higher metallicities of >0.1 Z⊙ (Wise et al. 2012a).

Our simulated galaxies are slightly more metal-poor (roughly a factor of 2) than the local dwarf population (Woo, Courteau & Dekel 2008). However, given the large scatter in the observed stellar metallicities, the difference is unlikely to be significant. Rather, we find that the simulated protogalaxies are more compact (by a factor of a few) than the local dwarf spheroids of comparable masses (e.g. Brodie et al. 2011). As can be observed in Fig. 2, the stellar

<figure><img src="image_7.png" /><figcaption>Figure 5. Predicted stellar mass as a function of halo mass in our simulations. Different colour codings and symbols correspond to nine different haloes. Note that we plot the results at various redshifts, and thus this may also be seen as an evolutionary sequence at $7 \le z < 18$. The empty symbols indicate the haloes hosting only Pop III stars, while the haloes with Pop II stars are shown as filled symbols. We also include the stellar mass in haloes outside the main halo if they are still within the zoom-in region and not contaminated by coarse dark matter particles. We find that our results follow the slope and normalization predicted by Kimm & Cen (2014) (black line), which reproduced the UV luminosity function at $z \sim 7$. The dashed line indicates an extrapolation of the Kimm & Cen (2014) results. Our results are also in fair agreement with simulations from Xu et al. (2016) (the grey squares). For comparison, we include the empirical sequence at $z \approx 0$ extrapolated to the mini-halo regime (Behroozi et al. 2013) (grey dotted line).</figcaption></figure>

components of high-z dwarves are dominated by a few clusters, and the resulting half-mass radii are found to be ~20-100 pc. The smaller size of the simulated galaxies may not be very surprising, given that the universe is denser and star formation per unit stellar mass is known to be more efficient at high redshifts (e.g. Speagle et al. 2014).

## 3.2 Escape fraction of LyC photons

We calculate the escape fraction by comparing the photon flux ($F_{ion}$) at the virial radius with the emergent flux from stars. Since photons travel with finite speed, we use the production rate at $t - R_{vir}/\tilde{c}$, where the time delay ($R_{vir}/\tilde{c}$) is roughly ~1 Myr. The escape fraction is then

$$f_{esc}(t) = \frac{\int d\Omega F_{ion}(t) \cdot \hat{r}}{\int dm_* \dot{N}_{ion}(t - R_{vir}/\tilde{c})}, \quad (27)$$

where $\Omega$ is the solid angle, $m_*$ is the mass of star particles, $\dot{N}_{ion}$ is the production rate of ionizing radiation per unit mass. Note that the choice of the radius (i.e. $R_{vir}$) at which the measurement is made is conventional, but it is desirable to measure the escape fraction at a radius large enough that it can be combined with estimates of the clumping factor from large-scale simulations (e.g. Pawlik, Schaye
& van Scherpenzeel 2009; Finlator et al. 2012; Shull et al. 2012; So et al. 2014) to study the reionization history of the universe.

### 3.2.1 A high average escape fraction in mini-haloes

Fig. 6 shows that LyC photons escape from mini-haloes quite efficiently after a burst of star formation. Not only Pop III stars, which are characterized by a constant photon production rate and an abrupt decrease, but also Pop II stars, characterized instead by an exponentially decaying rate, show a high escape fraction of $f_{\text{esc}} \sim 30–40$ per cent. In some cases, the escape fraction remains very low even after the formation of Pop II star clusters. Similarly, not all Pop III stars lead to a high escape fraction. Fig. 7 displays the photon number-weighted escape fraction of individual Pop III stars with different masses. It can be seen that only massive Pop III stars with $M_{\text{PopIII}} \gtrsim 100 M_{\odot}$ are able to provide LyC photons to the IGM (Whalen, Abel & Norman 2004), while almost all of the ionizing radiation from Pop III stars with $M_{\text{PopIII}} \lesssim 70 M_{\odot}$ is absorbed inside the virial radius.

The high escape fraction can be associated with the blow out of birth clouds due to radiation feedback (i.e. photoionization plus direct radiation pressure). This is especially evident for Pop III stars, as they tend to form in an isolated fashion and radiation is the only energy source while they are emitting LyC photons. Even for Pop II star clusters, we find that radiation feedback is the main culprit for creating the low-density, ionized channels through which LyC photons can escape. This is supported by the short time delay ($\lesssim 5$ Myr) between the peak of the photon production rate and the peak of the escape fraction. Even though the youngest SN occurs after 3.5 Myr, the stochasticity in our random sampling of the lifetime of SN progenitors is unlikely to explain the short delay. To substantiate this further, we show an example in Fig. 8 where the escape fraction increases from $f_{\text{esc}}(t_0) = 0$ per cent to $f_{\text{esc}}(t_1) \sim 20$ per cent within 3.7 Myr, during which no SN explosions occur. The dense, star-forming clouds are disrupted and LyC photons propagate to the virial radius, ionizing the neutral hydrogen in the halo (second row). Note that only the birth clouds are dispersed, while the average density of the halo gas is little affected by radiation.

We find that SN explosions enhance the escape of LyC photons from time to time by ejecting gas from the dark matter halo. As an illustration, in Fig. 8, we show the projected density distributions and the ionization fraction of hydrogen at several different epochs. After the birth clouds are dispersed and lifted by radiation feedback ($t = t_1$), the density of the gas beyond the galaxy actually increases, obscuring the LyC photons in the halo region ($t = t_2$). Once this gas is completely pushed out from the halo, the column density of neutral hydrogen along these solid angles becomes very small ($t = t_3$, bottom middle panel). As a result, even though the projected ionized fraction at $t = t_3$ appears to be smaller than at the previous stage ($t = t_2$), the actual escape fraction is larger. Nevertheless, the effect of SNe by creating the secondary peak does not play a significant role in increasing the total number of escaping photons, as the stellar populations become too old to generate a large amount of LyC photons. It should be noted, however, that the effects of SNe may be more substantial if several star clusters form simultaneously in more massive haloes ($M_{\text{halo}} \gtrsim 10^8 M_{\odot}$) and SNe in slightly older clusters generate strong winds that strip off gas in other star-forming clumps (e.g. Kimm & Cen 2014).

In Fig. 9, we show the average escape fraction as a function of the dark matter halo mass. As demonstrated in previous studies (Kimm & Cen 2014; Wise et al. 2014; Paardekooper et al. 2015), the

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Evolution of the escape fraction (black lines) and the photon production rates ($N_{\text{ph}}$) in units of $\text{#s}^{-1}$ (cyan shaded regions) in seven different haloes as a function of the age of the Universe. The two vertical dashed lines mark the time at which the mass of each halo becomes $10^7$ or $10^8 M_{\odot}$, respectively. The photons produced by Pop III and Pop II stars can be distinguished by the shape of the photon production rate; only Pop III stars exhibit a squarish evolution, because no LyC photons are generated once they explode or implode. Mini-haloes show a high escape fraction in general, although a significant variation can be seen. The escape fraction often exhibits a double peak for individual star formation events, which is explained in detail later in Fig. 8.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Escape fraction of individual Population III stars. We measure the mean escape fraction of Pop III stars by weighting the photon production rate during their lifetime. More massive Pop III stars show a higher escape fraction. We find that most LyC photons from low-mass Pop III stars ($M_{\odot} < 60 M_{\odot}$) are absorbed within the dark matter halo.</figcaption></figure>

instantaneous escape fraction varies significantly for a given halo mass. However, the photon number-weighted mean escape fraction in the mini-halo regime ($M_{\text{vir}} \le 5 \times 10^7$) is found to be generally very high (~20–40 per cent, Table 4). Note that this is a factor of a few larger than the average escape fraction predicted in atomic-cooling haloes ($f_{\text{esc}} \sim 10$ per cent; Kimm & Cen 2014), indicating that there is a dependence on the halo mass. It is also interesting to note that our results are in good agreement with Wise et al. (2014) and Xu et al. (2016), despite the significant differences in the modelling of star formation and feedback. Nevertheless, both sets of studies (this work and Wise et al. 2014 and Xu et al. 2016) allow for rapid star formation that leads to the disruption of gas clumps. Wise et al. (2014) adopted 7 per cent efficiency for star formation within a sphere of mean density $n_H \sim 10^3 \text{ cm}^{-3}$ per dynamical time, while our simulations employ an efficiency that varies (~5–20 per cent) according to the local turbulent and gravitational conditions. The common prediction of the high $f_{\text{esc}}$ suggests that there is a dominant feedback process, which is included in both studies, shaping the escape fraction in the mini-haloes. Indeed, as will be detailed later, we find that photoionization is the main culprit for the effective leakage of LyC photons (see Section 4).

### 3.2.2 Low escape fraction due to slow destruction of birth clouds

Even though the leakage of LyC photons in the mini-haloes is very efficient on average, mini-haloes are sometimes optically thick for the LyC photons to escape (Fig. 6). In order to understand the origin of the occasionally low escape fraction, we compute the optical depth centred on each young (<10 Myr) star by casting 3072 rays per stellar particle using the HEALPIX algorithm (Górski et al. 2005). We adopt the absorption cross-section for neutral hydrogen from Osterbrock & Ferland (2006), and the Small Magellanic Cloud-type extinction curve for dust (Weingartner & Draine 2001) assuming that the dust mass constitutes 40 per cent of the metal mass (e.g.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Origin of the double-peaked escape fraction for individual star formation events. The double-peaked feature is shown in the bottom right panel for clarity. The first column of the mosaic shows the evolution of the projected density distributions in the central region of the halo at four different times, as marked in the bottom right panel. The second and the third columns display the projected density and the fraction of ionized hydrogen distributions within a virial radius. The escape fraction at each time is also indicated in the top right corner. The bottom left panels show the density slice at two different epochs. One can see that the escape fraction becomes very high once the central star-forming clump is destroyed, and then drops as massive stars end their life. The escape fraction increases again when SNe blow out gas from the halo.</figcaption></figure>

Dwek 1998; Draine & Li 2007). We neglect the contribution from metals residing in hot gas with $T > 10^6$ K, as they are likely to be thermally sputtered (e.g. Draine & Salpeter 1979).

In Fig. 10, we examine the minimum distance within which 50 per cent and 99.9 per cent of the LyC photons would be absorbed. Also included in the bottom panel is the local density at which each young star is located. Since we are interested in the origin of the low escape fraction even after a burst of star formation, we restrict our analysis to $11 < z < 13$ in H1, $11.2 < z < 12$ in H2, $9.5 < z < 10$, $7 < z < 7.5$ in H3 and $7.5 < z < 8$ in H6. The plot demonstrates that the escape fraction is low in two circumstances. First, a high density of the local environment ($n_H \gtrsim 10^3 M_{\odot}$) indicates that stars are either just born out of dense clouds or radiation
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Escape fraction of LyC photons as a function of halo mass. Each point corresponds to the escape fraction in different haloes at various redshifts. Larger photon production rates are shown as redder colours. The photon number-weighted average is displayed as the star symbols with error bars representing the interquartile range (see Table 4). We find that the photon-weighted escape fraction is high (~20–40 per cent) in the mini-haloes, in good agreement with the previous studies (Wise et al. 2014; Xu et al. 2016), even though there is a considerable scatter in instantaneous measurements.</figcaption></figure>

<table><caption>**Table 4.** Photon number-weighted $f_{esc}$ at $7 \le z \le 15$.</caption><thead><tr><th>log M<sub>vir</sub></th><th>⟨f<sub>esc</sub>⟩</th></tr></thead><tbody><tr><td>6.40</td><td>0.373<sup>+0.206</sup><sub>-0.350</sub></td></tr><tr><td>6.83</td><td>0.391<sup>+0.211</sup><sub>-0.339</sub></td></tr><tr><td>7.07</td><td>0.360<sup>+0.228</sup><sub>-0.219</sub></td></tr><tr><td>7.49</td><td>0.230<sup>+0.086</sup><sub>-0.074</sub></td></tr><tr><td>7.95</td><td>0.050<sup>+0.147</sup><sub>-0.026</sub></td></tr></tbody></table>

is not strong enough to blow away the surrounding medium. In this case, the rays from these young star particles typically propagate to a small distance ($r \lesssim 10$ pc), implying that they get absorbed within the birth clouds. Such stars account for ~70 per cent of the low escape fraction. Secondly, stars with older ages ($t \gtrsim 2.5$ Myr) show progressively larger minimum distances (the middle panel), meaning that the birth clouds are dispersed at later stages. The trend that older stars are located in lower density environments supports this interpretation (the bottom panel). Once the dense structures are destroyed, photons propagate out to a larger radius ($r \gtrsim 10$ pc), but only a small fraction of the total ionizing radiation appears to be able to do so. The similarity of the density PDFs for 50 per cent and 99.9 per cent absorption probability in the top panel suggests that the photon luminosity is sufficiently weak for the ionization front to stall inside the halo. The composite image of H5 in Fig. 2 illustrates an example of such a case. Thus, the escape fraction can sometimes be low even after a burst of star formation mainly because the birth clouds are cleared away too slowly compared to the lifetime of massive stars (~5 Myr).

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Top panel: PDFs of the distance within which photons from stars younger than 10 Myr are absorbed when the escape fraction is low ($f_{esc} < 10^{-4}$) even after starbursts. These, for example, correspond to $11 < z < 11.5$ in H2 or $7 < z < 8$ in H6 (see the text). Also included as a dark red line is the PDF of the distance within which 99.9 per cent of the photons from stars of age $5 < t < 10$ Myr is absorbed. Bottom panel: PDF of the local density at which the young stars are located.</figcaption></figure>

# 4 DISCUSSION

## 4.1 Physical mechanism for driving the escape of LyC photons in mini-haloes

In the previous section, we showed that the high escape fraction is best explained by strong radiation feedback. This raises the question, which mechanism within this category is primarily responsible for regulating the escape fraction (i.e. photoheating, direction radiation pressure from UV photons or multiply scattered IR photons)? Whether or not multiply scattered IR photons govern the dynamics of the ISM in our primordial galaxies can be tested by examining the trapping factor in the dense regions of the ISM. The typical density of the star-forming clouds of size ~10 pc is ~$10^5$ cm$^{-3}$ in our simulations, and their metallicity ranges from $Z \sim 10^{-4} Z_{⊙}$ to $10^{-2} Z_{⊙}$. Assuming a dust opacity of $\kappa_{sc} \sim 5 (Z/Z_{⊙})$ cm$^2$ g$^{-1}$ (Semenov et al. 2003), IR photons in the star-forming clouds have a maximum optical depth of $\tau_d = \kappa_{sc} \Sigma_{gas} \lesssim 0.3$. Note that $\tau_d$ will be reduced further if the smaller scattering cross-section is considered at temperatures lower than ~100 K (Semenov et al. 2003). This indicates that the IR photons are not efficiently trapped inside the ISM, and that reprocessing gives a negligible momentum boost in the metal-poor, mini-haloes (cf. the min$\rho_{Q}$ case in Bieri et al. 2017).

The effects of photoheating and direct radiation pressure are difficult to disentangle, because they operate simultaneously. Nevertheless, it is possible to estimate the maximum extent within which each process can balance the external pressure in the ISM. Rosdahl & Teyssier (2015) show that if the ambient gas temperature ($T_0$) is significantly lower than that of the H II bubble ($T_{ion} \approx 2 \times 10^4$ K),
<figure><img src="image_13.png" /><figcaption>Figure 11. Comparison between the maximum radius within which photoionization can counter-balance the external pressure ($r_{PH}$) and the radius within which direct radiation pressure can overcome the external pressure ($r_{DP}$). Only regions cooler than $T \le 2 \times 10^4$ K are shown. Darker colours indicate the regions with stronger radiation field. The plot suggests that photoionization is more important mechanism than direct radiation pressure in most regions, although direct pressure can be slightly more effective than photoheating at emptying the gas in very dense regions.</figcaption></figure>

the gas density inside the bubble can be lowered due to overpressure, and the extent to which ionizing radiation can balance recombination becomes larger by a factor of $(T_{ion}/T_0)^{2/3}$. Thus, the gas can be pressure-supported within $r_{PH}$,

$$r_{PH} \approx 26 \text{ pc} \left( \frac{m_{\text{star}}}{10^3 \text{ M}_{\odot}} \right)^{1/3} \left( \frac{n_{\text{H},0}}{10^3 \text{ cm}^{-3}} \right)^{-2/3} \left( \frac{T_{\text{ion}}}{10^4 \text{ K}} \right)^{2/3} \left( \frac{T_0}{10^2 \text{ K}} \right)^{-2/3}, \quad (28)$$

where $n_{H,0}$ is the density of the ambient medium. On the other hand, the absorption of ionizing radiation directly imparts momentum that offsets the external pressure within $r_{DP}$,

$$r_{DP} = \sqrt{\frac{L}{4\pi n_{H,0} c k_B T_0}} \\ \approx 6.4 \text{ pc} \left( \frac{m_{\text{star}}}{10^3 \text{ M}_{\odot}} \right)^{1/2} \left( \frac{n_{\text{H},0}}{10^3 \text{ cm}^{-3}} \right)^{-1/2} \left( \frac{T_0}{10^2 \text{ K}} \right)^{-1/2}. \quad (29)$$

Here we use $L = 2 \times 10^{36} \text{ erg s}^{-1} \text{ per 1 M}_{\odot}$, adequate for a simple stellar population with age younger than 5 Myr. As pointed out Rosdahl & Teyssier (2015), the direct radiation pressure from a massive star cluster dominates over the pressure from the warm ionized gas only in a very dense medium, as $r_{PH}/r_{DP} \propto (L n_{H,0})^{-1/6}$.

In Fig. 11, we compare the two radii given by equations (28) and (29). We take the external pressure as the mass-weighted pressure in the neighbouring 26 cells. This approach is, no doubt, a simplification, given that the immediate neighbours may already be ionized and cannot represent the cold surroundings. Nevertheless, this is justifiable because we are interested in the mostly neutral region

<figure><img src="image_14.png" /><figcaption>Figure 12. Effects of photoheating (PH) and radiation pressure from UV photons (DP). The figure shows a slice of density (first and third rows) and ionized hydrogen fraction distributions (second row) from the simulations with and without PH and/or DP right after a Pop III star with $\approx 150 \text{ M}_{\odot}$ is formed. The virial radius of the mini-halo is $R_{vir} \sim 350 \text{ pc}$, which is approximately the size of the panels in the first and second rows. The panels in the third row are the zoomed-in images of the panels in the first row. It is evident that photoheating is responsible for the efficient escape of LyC photons.</figcaption></figure>

where the impact of direct radiation pressure can be strong. The intensity of ionizing radiation is calculated by combining the luminosity of each star residing in each cell. Note that we only consider cells with temperatures low enough to contain neutral hydrogen ($T \le 2 \times 10^4$ K). Fig. 11 suggests that photoheating is primarily responsible for blowing gas away in the ISM at densities lower than $n_H \sim 10^4 \text{ cm}^{-3}$, whereas direct radiation pressure dominates over photoheating at very high densities. It is worth noting that at such high densities, pressure equilibrium due to direct radiation pressure is achieved on very short time-scales (i.e. less than an Myr; see fig. 4 in Rosdahl & Teyssier 2015). Therefore, photoheating is more likely to be responsible for the decrease in gas density than direct radiation pressure during the lifetime of massive stars (cf. Haehnelt 1995).

To isolate the effects of photoionization heating and radiation pressure, we identify a Pop III forming halo showing a high escape fraction ($f_{esc} \sim 10$ per cent) and run three simulations varying the included feedback processes. In the simulation without photoheating (DP only), we compute the non-equilibrium chemistry without increasing the temperature by ionization of hydrogen and helium. Fig. 12 shows slices of density and ionized hydrogen fraction immediately preceding the explosion of a Pop III star with $M \approx 145 \text{ M}_{\odot}$. It is evident from this figure that the efficient escape of LyC photons can be attributed to heating by photoionization. The central star-forming cloud is effectively destroyed only in the runs with photoheating.

Interestingly, we find that the run without photoheating results in $f_{esc} = 0$, even though direct radiation pressure creates low-density holes in the central region of the dark matter halo. The density of the central hole is even lower than the run with photoheating,
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Origin of the high escape fraction. Top panel shows the photon number-weighted escape fraction of each star formation event as a function of the peak photon production rate ($\dot{N}_{LyC}$), which is essentially proportional to star formation rate. Each event involving Pop II and III stars is shown with different symbols, as indicated in the legend. Redder colours denote a larger $\dot{N}_{LyC}$. The escape fraction is weakly correlated with $\dot{N}_{LyC}$. Bottom panel: same as the top panel, but with respect to the peak photon production rate per lowest 10 per cent dense gas mass from each solid angle ($\frac{dN_{LyC}}{d\Omega} / \frac{dN_{dense}}{d\Omega}$). Using the lowest 10 per cent is motivated by the fact that LyC photons usually escape through a small opening angle with little neutral hydrogen. Here we define the dense gas as $n_H \ge n_{H, strom}$, where $n_{H, strom}$ is the minimum density below which gas can be fully ionized for a given $N_{LyC}$. The escape fraction is higher if there are more photons available per dense gas mass. We interpret this as being mainly due to photoionization feedback being more effective when there is less dense gas.</figcaption></figure>

because the UV photons accelerate the central gas more effectively due to lower temperature in the early stage of the expansion. Note that in order to receive a momentum boost from direct radiation pressure, atoms must be neutral and thus able to absorb ionizing radiation. Because UV photons from the Pop III star propagate and push the gas only in the radial direction, a dense shell structure that continually absorbs the ionizing radiation is formed. This is an important difference from the run that includes photoheating where the gas becomes overpressurized and develops low-density channels naturally.

What then determines the escape fraction of LyC photons in the mini-haloes? Given the importance of photoheating, the escape fraction is likely to depend on the strength of the radiation field. To investigate this issue, we measure the photon number-weighted escape fraction and the peak production rate of LyC photons for individual star formation episodes. Fig. 13 (top panel) shows that the escape fraction is not strongly correlated with the peak production

rate of LyC photons, indicating that it is not a simple function of star formation rate. However, taking the sample of Pop III stars alone, there is a trend that more massive stars show higher escape fractions, as seen in Fig. 7. Given that the masses of their host dark matter haloes are more or less the same (6.25 $\lesssim$ log $M_{halo}$ $\lesssim$ 7), we hypothesize that the effect of photoheating is likely to be mitigated if stars are enshrouded by a large amount of dense gas. To test this idea, we measure the gas mass for each solid angle assuming that a star cluster is a point source within the virial sphere. We exclude the gas less dense than $n_{strom}$, which is the maximum density below which recombination can be offset by the photons from Pop III or Pop II clusters. In addition, since photons in the outskirts of dense clouds are likely to escape more easily and LyC photons usually escape through a small opening angle with a low neutral fraction (e.g. Kim et al. 2013; Cen & Kimm 2015), we use the lowest 10 per cent of the projected mass for each solid angle. Note that this is essentially the same as probing when the escape fraction is greater than 10 per cent. The bottom panel of Fig. 13 reveals that mini-haloes exhibit a high escape fraction when they form a lot of stars while there is a small amount of enshrouding gas. This is not particularly surprising and is also physically sensible because the more photons the halo has per dense gas mass, the more likely the ISM will be ionized and overpressurized. One might wonder whether or not the correlation arises simply due to the slim geometry of enshrouding gas. We confirm that taking the mean of the dense gas mass for each solid angle rather than the lowest 10 per cent mass does not wash out the trend, but only weakens it slightly, suggesting that geometry is a secondary effect.

## 4.2 Importance of mini-haloes for reionization

Our simulations represent only a small volume of the observable Universe, and thus it is not possible to calculate $\tau_e$ directly. Instead, we take a simple analytic approach to infer the importance of the mini-haloes to the reionization history by computing the evolution of the mass fraction of ionized hydrogen ($Q_{HII}$), as

$$ \frac{dQ_{HII}}{dt} = \frac{\dot{n}_{ion}}{\langle n_{H} \rangle} - \frac{Q_{HII}}{t_{rec}(C_{HII})}, \qquad (30) $$

where $\langle n_H \rangle$ is the mean density of the universe in comoving scale. The recombination time-scale ($t_{rec}$) is a function of a clumping factor ($C_{HII}$) and temperature, as

$$ t_{rec}(C_{HII}) = [C_{HII} \alpha_B(T) f_e(z) \langle n_H \rangle (1+z)^3]^{-1}, \qquad (31) $$

where a correction factor ($f_e(z)$) is included to account for the additional contribution from 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 ‘effective’ clumping factor of $C_{HII} = 1 + \exp(-0.28 z + 3.59)$ at $z \ge 10$ or $C_{HII} = 3.2$ at $z < 10$ (Pawlik et al. 2009). Once $Q_{HII}$ is determined, the Thomson optical depth is calculated as

$$ \tau_e(z) = c \langle n_H \rangle \sigma_T \int_0^z f_e(z) Q_{HII}(z') \frac{(1+z')^2 dz'}{H(z')}. \qquad (32) $$

The most crucial term in equation (30) is the photon production rate per unit volume $\dot{n}_{ion}$ (#Mpc$^{-3}$ s$^{-1}$). To compute this quantity, we first generate dark matter halo mass functions as a function of redshift (Sheth & Tormen 2002), and assign stars to each dark matter halo by taking into account the following factors:

(1) Star formation efficiency: Kimm & Cen (2014) performed cosmological radiation-hydrodynamic simulations, and showed that
the UV luminosity function at $z \sim 7$ can be reproduced with efficient SN feedback. They find that the stellar mass to the dark matter halo mass may be approximated by the fit,

$$ \log M_{\text{star}} = -8.08 + 1.55 \log M_{\text{halo}} \quad (33) $$

(see Fig. 5, the black line). Our simulations suggest that the relation can be extended even to the mini-halo regime at $7 \lesssim z \lesssim 15$. Thus, we use this fit to compute the total stellar mass formed in each dark matter halo assuming no redshift dependence.

(2) Occupation fraction: not all of the mini-haloes hosts Pop II stars. This is essentially because early stars pre-heat the IGM and suppress the accretion on to the halo (Gnedin 2000). Moreover, the formation of molecular hydrogen, which is a pre-requisite for gas collapse in mini-haloes in the absence of metals, is suppressed by the Lyman-Werner radiation from the neighbouring haloes. Furthermore, gas in haloes can be blown out by efficient feedback processes from Pop III stars. To account for the occupation fraction in haloes with $M_{\text{halo}} \ge 10^7 M_{\odot}$, we use a redshift-dependent fit to the fraction of haloes that host stars in the Wise et al. (2014) simulations (see their fig. 2). This corresponds to $f_{\text{host}} \approx 1$ for $M_{\text{halo}} \ge 10^7 M_{\odot}$ at $z \gtrsim 17$, while $f_{\text{host}}$ drops to less than 0.01 in haloes with $10^7 M_{\odot}$ at $z \sim 7$. For haloes with $M_{\text{halo}} \le 10^{6.25} M_{\odot}$ and $M_{\text{halo}} \ge 10^8 M_{\odot}$, we adopt a fixed value of $f_{\text{host}} = 2.4 \times 10^{-4}$ and unity, respectively (Wise et al. 2014). The occupation fraction is considered by simply reducing the number density of dark matter haloes that can form stars in estimating $\dot{n}_{\text{ion}}$ (equation 30).

(3) Escape fraction: our estimates of the escape fraction in mini-haloes are in good agreement with the results performed with ENZO (Wise et al. 2014; Xu et al. 2016). However, the escape fractions in atomic-cooling haloes appears to be less certain, and thus we compare three different models: (i) an escape fraction that slowly decreases with the halo mass (Kimm & Cen 2014; Xu et al. 2016), (ii) an escape fraction showing a more rapid decline (e.g. Ma et al. 2015; Paardekooper et al. 2015) and (iii) an escape fraction with a lower limit of 20 per cent. For the fiducial case, we employ the fit,

$$ \log f_{\text{esc}} = 1.0 - 0.2 \log M_{\text{halo}}, \quad (34) $$

whereas a broken law is adopted for the ‘low’ case

$$ \log f_{\text{esc}} = \begin{cases} 0.17 - 0.2 \log M_{\text{halo}} & (M_{\text{halo}} > 10^{8.5}) \\ 9.52 - 1.3 \log M_{\text{halo}} & (10^{7.75} < M_{\text{halo}} \le 10^{8.5}) \\ 1.0 - 0.2 \log M_{\text{halo}} & (M_{\text{halo}} \le 10^{7.75}) \end{cases}. \quad (35) $$

The last model has the same escape fraction as the fiducial case, but with the minimum of $f_{\text{esc}} = 20$ per cent, as is often assumed in reionization studies. The three models are shown in Fig. 14 (the top left panel). Also included as a dot-dashed line is a model with an escape fraction that would be necessary to yield a later end of reionization as perhaps suggested by the Ly$\alpha$ opacity data presented by Becker et al. (2015).

$$ \log f_{\text{esc}} = 1.292 - 0.245 \log M_{\text{halo}}. \quad (36) $$

(4) Metallicity-dependent $\dot{N}_{\text{ion}}$: metal-rich main-sequence stars are cooler than metal-poor stars, and because of the larger opacity in the stellar atmosphere, an SSP of solar metallicity produces approximately a factor of 2.5 fewer LyC photons than an SSP with $0.02 Z_{\odot}$. In order to not over-estimate the photon budget from metal-rich, massive galaxies, we take into consideration the metallicity-dependent $\dot{N}_{\text{ion}}$ by using the local mass-metallicity relation for dwarf galaxies (Woo et al. 2008), as $log Z =$

$-3.7 + 0.4 \log (M_{\text{star}}/10^6 M_{\odot})$. Although this is certainly a simplification that needs to be tested against upcoming high-$z$ observations, we note that our simulated galaxies are metal-poor by about a factor of 2 compared to the relation, well within the scatter of the measurements. Once the metallicity of stars is determined, we use the metallicity-dependent number of LyC photons per $M_{\odot}$, as

$$ \log Q_{\text{LyC}} = 60.31 - 0.237 \log Z, \quad (37) $$

which is appropriate for the Padova asymptotic giant branch (AGB) models with a Kroupa IMF (Leitherer et al. 1999). Here we adopt a maximum mass cut-off of 120 $M_{\odot}$ for the metallicity range $0.0004 \le Z \le 0.05$. We do not extrapolate $Q_{\text{LyC}}$ for metallicities lower than $0.02 Z_{\odot}$. It is worth noting that for $Z = 0.004$, the Padova AGB model predicts $Q_{\text{LyC}} = 8.86 \times 10^{60}$, while the binary evolutionary model with a 300 $M_{\odot}$ cut-off gives a slightly higher estimate of $Q_{\text{LyC}} = 9.26 \times 10^{60}$ (Stanway, Eldridge & Becker 2016). Stellar models with strong rotation seem to generate even more photons $Q_{\text{LyC}} = 11.2 \times 10^{60}$ at $Z = 0.002$ (Topping & Shull 2015), but this leads to only a minor increase ($\sim 20$ per cent) in the LyC production rate, which does not change our conclusions significantly.

(5) Pop III stars: the largest uncertainty is perhaps the contribution of Pop III stars to reionization. Not only is the IMF of the Pop III stars not well constrained, but also only a handful of Pop III stars are simulated in this work to estimate their escape fraction. Thus, we examine two extreme cases: one in which Pop III stars do not contribute to reionization at all; and the other in which 40 per cent of the photons from Pop III stars with $M_{\text{PopIII}} \gtrsim 100 M_{\odot}$ are assumed to contribute to reionization. The former can be regarded as the case where only low-mass Pop III stars with $M_{\text{PopIII}} \lesssim 100 M_{\odot}$ form, whereas the latter corresponds to the case where one allows for the wide range of masses from 10 to 1000 $M_{\odot}$. Note that the assumption of the mass-dependent escape fraction is motivated by Fig. 7. The total number of LyC photons that Pop III stars can generate is computed by convolving the IMF (equation 7) and the escape fraction. This yields $\langle N_{\text{LyC,PopIII}} \rangle = 1.8 \times 10^{62} M_{\odot}^{-1}$. We then adopt the formation rate of Pop III from the large volume simulations of Xu et al. (2016) (the ‘Normal’ run) by using the fit, $\log \dot{M}_{\text{PopIII}}/(M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}) = -6.428 + 0.275 z - 0.011 z^2$. This rate is roughly an order of magnitude smaller than the rate employed in Wise et al. (2014), which is based on the simulation of an overdense region (‘Rare Peak’; Xu, private communication).

Once the number of photons escaping from a dark matter halo ($N_{\text{LyC}}^{\text{esc}}$) as a function of $M_{\text{halo}}$ and $z$ is computed, the total number of escaped photons until $z$ can be obtained by multiplying the dark matter halo mass function with $N_{\text{LyC}}^{\text{esc}}$. We use the Sheth & Tormen (2002) mass function for the mass range $M_{\text{halo}} = 10^{6.25} - 10^{12} M_{\odot}$. The time derivative of the integrated quantity corresponds to $\dot{n}_{\text{ion}}$ in equation (30).

The top right panel in Fig. 14 shows the evolution of $\dot{n}_{\text{ion}}$ for the different models. In the cases with the fiducial escape fraction, the number of LyC photons that escape from their host halo peaks at $z \sim 5$ and decreases at lower redshift, in remarkably good agreement with the observational findings at $3 \le z \le 4.7$ (Becker & Bolton 2013). Note that the peak is earlier than the peak of the star formation rate density ($z \sim 3$; e.g. Hopkins & Beacom 2006), because the photon budget at $z \lesssim 6$ is determined by more massive, metal-enriched galaxies where fewer ionizing photons are generated per mass of stars formed than in dwarf galaxies. The peak would shift to even earlier redshifts if the escape fraction decreases at lower
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Reionization models based on simple analytic calculations. We examine three different escape fractions (models I, II and III), the presence of the mini-haloes (model V) or Pop III stars (model VII) and the occupation fraction of the mini-haloes (model VI). Also included as model IV is the escape fraction favoured by the patchy reionization picture (Becker et al. 2015). Top left: different assumptions about the escape fraction used to calculate the reionization history. Different symbols indicate the escape fraction from different theoretical studies, as indicated in the legend. The solid, dashed and dotted line represent the fiducial-, high- and low-escape fraction model, respectively. Top right: the photon production rate per unit volume ($dn_{LYC}/dt$) as a function of redshift. We use the Sheth & Tormen (2002) mass functions and our theoretical results (the dashed line in Fig. 4, equation 33) to compute $dn_{LYC}/dt$. Bottom left: the electron optical depth ($\tau_e$) in different models. The optical depths inferred from the Planck experiments are shown as shaded regions (light grey: Planck Collaboration XIII 2016a; dark grey: Planck Collaboration XLVI 2016b). Bottom right: the fraction of ionized hydrogen. Our fiducial model predicts the reionization redshift of 6.7 and $\tau_e = 0.067$. We find that the mini-haloes produce a rather small amount of LyC photons, and thus are of minor importance for the reionisation of the Universe.</figcaption></figure>

redshifts. Although we assume a redshift-independent escape fraction, which is motivated by our simulations and simulations of Wise and collaborators at $z \gtrsim 7$, there is a hint that the escape fraction decreases with decreasing redshift. For example, from cosmological radiation-hydrodynamic simulations, Cen & Kimm (2015) measured that the escape fraction is decreased by a factor of $\sim 2$ at $z \sim 4$, compared to galaxies at $z \gtrsim 7$ (Kimm & Cen 2014). They attribute this to the fact that birth clouds are more slowly disrupted, as star formation becomes less stochastic at lower redshift. Thus, we caution the readers that $\dot{n}_{ion}$ at $z \le 6$ is probably still uncertain by a

factor of a few. Our redshift-independent assumption should thus be kept in mind in our discussion of the reionization history of the Universe at $z \gtrsim 6$.

Fig. 14 depicts the evolution of the electron optical depth ($\tau_e$) and the fraction of ionized gas ($Q_{HII}$). Several important conclusions can be drawn from this plot. First, the mini-haloes turn out to be of minor importance to reionization. Although their escape fraction is high, the photon production rate is intrinsically low, first because stellar feedback efficiently suppress star formation, and secondly because not every mini-halo hosts stars (Gnedin 2000; Wise et al. 2014).
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Left: contribution of photons from different halo masses to reionization as a function of redshift. We adopt the fiducial model for escape fractions (the black solid line in Fig. 14) and the stellar mass-to-halo mass relation from this work (the dashed line in Fig. 5). The photon production rate density from Pop III stars is shown as the black dashed line. The dotted line displays the maximum contribution from mini-haloes assuming that all of them host stars ($f_{\text{host}} = 1$). Right: the electron optical depth (top) and the fraction of ionized hydrogen (bottom) for models with different minimum halo masses, as indicated in the legend. The plots demonstrate that it is necessary to resolve haloes of mass $\approx 10^8 M_{\odot}$ to model reionization in simulations.</figcaption></figure>

We find that even if all the mini-haloes are assumed to host stars, $\tau_e$ will only increase slightly (model VI). Conversely, if we neglect the mini-haloes, there is no noticeable difference to the prediction of the reionization history. This is in stark contrast with the claim that mini-haloes provide enough photons in the early universe to achieve a high $\tau_e \gtrsim 0.1$ (Wise et al. 2014). The difference can be attributed to a large part to the reduced stellar masses in the smallest mini-haloes in this work, compared to Wise et al. (2014).

Similarly, we find that the contribution from Pop III stars is not substantial (see also Paardekooper et al. 2013, cf. Ahn et al. 2012; Maio et al. 2016). The model without Pop III stars (model VI) predicts a minor decrease (35 per cent and 10 per cent) in the number of ionizing photons at $z \sim 20$ and 10, respectively (see also Fig. 15). The resulting difference in $\tau_e$ is found to be much smaller than the $1\sigma$ error of the *Planck* measurements. This strongly suggests that even if the characteristic mass of Pop III stars is lower ($\sim 40 M_{\odot}$), which prohibits them from ionizing the IGM, the optical depth measurement may not be able to put strong constraints on the properties of Pop III stars (cf. Visbal, Haiman & Bryan 2015).

Instead, we find that the prediction of $\tau_e$ is primarily sensitive to the escape fraction of the atomic-cooling haloes ($M_{\text{halo}} \gtrsim 10^8 M_{\odot}$). Our fiducial model with the relatively high escape fraction and Pop III stars yields $\tau_e \approx 0.067$, which is in reasonable agreement with the *Planck* measurements (Planck Collaboration XIII 2016a; Planck Collaboration XLVI 2016b). We also find that the evolution of $Q_{H\bar{\imath}}$ lies between that of the ‘Late’ and the ‘Very Late’ reionization model of Choudhury et al. (2015) at $z \ge 8$ (see also Kulkarni et al. 2016), both of which are shown to be able to reproduce the rapid drop in the fraction of $Ly\alpha$ emitters at $z \gtrsim 7$ (Choudhury et al. 2015). On the contrary, the model with the Low’ escape fraction in the atomic-cooling regime fails to ionize the universe by $z \sim 3$ and underpredicts $\tau_e$, even in the presence of the mini-haloes and Pop III stars. This demonstrates that massive stars in atomic-cooling

haloes are likely to be the main sources responsible for reionizing the Universe. In order for this “Low’ escape fraction model to fully ionize the Universe by $z \sim 6$, a large number of AGN would be necessary at high redshifts to provide enough LyC photons (e.g. Madau & Haardt 2015). This scenario is certainly not yet ruled out, but it depends critically on the assumption about the evolution of the emissivity as a function of redshift, which is still a matter of debate (Chardin et al. 2015; Giallongo et al. 2015; Kashikawa et al. 2015; Kim et al. 2015).

When the escape fraction is assumed to always be greater than 20 per cent, the optical depth is predicted by our models to be rather high ($\tau_e = 0.075$) and reionization ends early ($z \sim 8$), which is less favoured by the latest *Planck* measurements ($\tau_e = 0.055 \pm 0.009$). Such a model also struggles to explain the evolution of the $Ly\alpha$ opacity (fluctuations) in QSO spectra at $z > 5$ as well as the rapid evolution of $Ly\alpha$ emitters (Becker et al. 2015; Chardin et al. 2015; Choudhury et al. 2015). In fact, even our fiducial model predicts reionization at $z \sim 6.7$, slightly earlier than in the fiducial model of Chardin et al. (2015) that fits well the photo-ionization rate inferred from $Ly\alpha$ forest data and perhaps suggesting that the escape of LyC photons in the atomic-cooling haloes needs to be suppressed further.

It should be noted, however, that there is a possibility that the simple analytic calculations are likely to predict too rapid propagation of $H\bar{\imath}$ bubbles in the late stage of reionization ($z \lesssim 8$) due to the assumption of infinite speed of light used in equation (30). In practice, photons from massive haloes will have to travel to encounter neutral hydrogen in regions devoid of sources, which inevitably delays the propagation of $H\bar{\imath}$ bubbles. If this is the case, reionization in the fiducial model is likely to end later than $z = 6.7$, bringing the optical depth in better agreement with the latest *Planck* results. However, the precise determination of the delay requires large-scale reionization simulations that are calibrated to match the photon production rates from haloes of different masses in this study, hence it is not clear yet how significant the effect will be. If the propagation of
H II bubbles is not significantly overestimated in our simple calculations, the escape fraction needs to be reduced further in more massive haloes in order for dwarf galaxies to reionize the universe only by $z \le 6$, as in our 'Late z_{reion}' model shown by the dot-dashed line in Fig. 14. Such escape fractions in atomic-cooling haloes are certainly within the uncertainty of current numerical results, and it would be worthwhile to re-examine escape fractions in more massive haloes ($M_{halo} \ge 10^{10} M_{\odot}$) with the physically well-motivated thermo-turbulent star formation model employed here, where stars form preferentially in gravitationally bound regions, and see how it compares with the results based on a simple density criterion for star formation (Kimm & Cen 2014).

In the left-hand panel of Fig. 15, we show the relative contribution of LyC photons from different components for our fiducial reionization model. It can be seen that more massive haloes dominate the photon budget at lower redshift in this model, and that reionization is driven by intermediate-mass atomic-cooling haloes ($10^8 M_{\odot} \lesssim M_{halo} \lesssim 10^{11} M_{\odot}$). The current theoretical estimate of Pop III star formation rates (Xu et al. 2016) suggests that Pop III stars are important only at $z \gtrsim 15$. The plot also demonstrates that the photon production rate density decreases at $z \le 5$ in our model, essentially because the number density of the intermediate-mass haloes with $10^8 M_{\odot} \lesssim M_{halo} \lesssim 10^{11} M_{\odot}$ does not increase anymore at lower redshift. One may wonder here whether the inclusion of more massive haloes with masses $M_{halo} \gtrsim 10^{12} M_{\odot}$ would change our conclusions, but we find that the photon production rate density from these haloes is significantly smaller than that from $M_{halo} \ge 10^{11} M_{\odot}$ at $z \ge 3$ if the local $M_{star}-M_{halo}$ relation (Behroozi et al. 2013) is assumed, as star formation becomes less efficient in haloes with $M_{halo} \ge 10^{12} M_{\odot}$.

Our results indicate that it will be necessary to resolve haloes of mass $\approx 10^8 M_{\odot}$ in order to correctly capture the expansion of ionized H II bubbles in large-scale simulations. The right-hand panels of Fig. 15 show that neglecting the contribution of LyC photons from haloes less massive than $10^9 M_{\odot}$ will only marginally delay the expansion, reasonably reproducing the electron optical depth. However, simulations that cannot resolve haloes with $M_{halo} \lesssim 10^{10} M_{\odot}$ would need to adopt higher escape fractions than the values derived in this work to provide photons early enough (e.g. Gnedin & Kaurov 2014; Pawlik, Schaye & Dalla Vecchia 2015; Ocvirk et al. 2016).

## 4.3 Caveats

Even though we include the most important physical ingredients in our simulations, it should be mentioned that several potential issues that can affect the predictions of star formation and the escape fraction are neglected.

First, we do not take into account the fact that the mass of the most massive star within a cluster is correlated with the total cluster mass (Weidner, Kroupa & Bonnell 2010; Kirk & Myers 2011). Given that the most massive stars are the most efficient at producing LyC photons, the number of photons per unit stellar mass (i.e. specific number) is likely to be overestimated in the simulated small Pop II clusters of mass $M \le 10^4 M_{\odot}$ where the IMF would not be fully sampled (see Kroupa et al. 2013, for a review). For example, at $Z = 0.0004$, the specific number of LyC photons from a $10^3 M_{\odot}$ cluster is smaller by ~60 per cent than that from a $10^4 M_{\odot}$ cluster. Neglecting this IMF sampling issue may have suppressed star formation more effectively than it should have. However, we do not expect that the escape fraction would be affected significantly by lowering $N_{LyC}$ by a factor of 2, given that the number of Lyα photons per dense gas mass spans three orders of magnitude already (Fig. 13). Furthermore, it should be noted that the grids of the stellar tracks do not extend to values lower than 0.02 $Z_{\odot}$. For comparison, star particles with $M < 10^4 M_{\odot}$ in our simulations often have the metallicity of $10^{-4}-10^{-3} Z_{\odot}$. We use the lowest grid point in order to not extrapolate the number of photons produced, but the use of more appropriate stellar grids is likely to enhance the Lyα photon production, possibly compensating the overestimation due to the incomplete IMF sampling.

Secondly, we neglect the contribution from runaway stars, which are thought to form through three-body interactions within clusters (e.g. Gies & Bolton 1986; Leonard & Duncan 1988; Fujii & Portegies Zwart 2011), and/or the ejection in a binary system after the explosion of the primary star (e.g. Blaauw 1961; Portegies Zwart 2000). Using an analytic model of disc galaxies, Conroy & Kratter (2012) claim that the inclusion of the runaway stars can, in principle, enhance the escape fraction by up to an order of magnitude if $f_{esc}$ without them is very small. However, Kimm & Cen (2014) compare the radiation-hydrodynamic simulations with and without the runaways, and show that the increase in $f_{esc}$ due to runaways is small (~20 per cent level) mainly because $f_{esc}$ is already high ($f_{esc} \sim 10$ per cent). In the same context, we do not expect that the runaways significantly affect the prediction of $f_{esc}$ in mini-haloes where $f_{esc}$ is even higher (~20–40 per cent).

Thirdly, our predictions are subject to the stellar evolutionary models. We take the production rate of LyC photons from the Padova model (Leitherer et al. 1999; Girardi et al. 2000) that is based on the evolution of single stars without rotation. However, it is known that many, if not all, massive stars live in binaries (e.g. Sana et al. 2012), and that their spin is non-negligible. The interactions of binary stars cannot only increase the total number of LyC photons, but also make the decline of the production rate slower, as the primary star transfers gas to the secondary, removes the hydrogen envelope and the stars merge (Stanway et al. 2016). Rotation can also increase the total number of ionizing photons and the lifetime of massive stars by regulating the mixing and fuelling of gas into the stellar core (e.g. Ekström et al. 2012). As aforementioned, the adoption of the recent stellar models can increase the production rate of LyC photons by ~10–20 per cent for an SSP (Topping & Shull 2015; Stanway et al. 2016), and is therefore not very important for the total photon production. However, these effects can potentially increase the escape fraction to the level of ~20 per cent in galaxies hosted by massive haloes with $10^{10}-10^{11} M_{\odot}$ (Ma et al. 2016), as stars older than $t \gtrsim 3$ Myr produce a non-negligible fraction of LyC photons. This could potentially bring the results of the simulations in tension with the rather late end of reionization suggested by the CMB data as well as the Lyα forest and Lyα emitter data, as dwarf galaxies would produce too many photons (see Fig. 14). Future radiation-hydrodynamic simulations with the self-consistent modelling of binary populations are necessary to better understand the role of (massive) binary stars and to see if theory and observations can still be reconciled.

Finally, due to the limited computational resources we had available, the work presented here lacks the systematic investigation of the effects of numerical resolution on the prediction of the escape of LyC photons. However, we note that our resolution and refinement strategy allows for star formation in sufficiently dense and gravitationally bound regions, which should be essential to not overestimate the escape fraction (e.g. Ma et al. 2015). Nevertheless, in order to make sure that this is really the case, we perform a resolution test by re-running the H1 halo with one fewer level of refinement (i.e. 1.4 pc resolution). We find that the escape fraction
in the mini-halo is still high ($f_{esc} \sim 40$ per cent), although the stellar mass is increased by 30 per cent. We also note that the recent work by Xu et al. (2016) shows that their 19 pc (comoving) resolution simulations ($\sim 2$ pc, physical) yield the results consistent with their previous simulations with a higher resolution of 1 pc (comoving; Wise et al. 2014). This suggests that the prediction of the escape fraction and star formation rates is likely to be converged in (sub-)parsec-scale simulations, as the physics of the collapse and subsequent disruption of star-forming clouds is reasonably well captured.

## 5 CONCLUSIONS

Using the radiation-hydrodynamic simulation code, RAMSES-RT, we investigate the importance of mini-haloes with $10^{6.25} \lesssim M_{halo} \lesssim 10^8 M_{\odot}$ for the reionization history of the Universe. For this purpose, we run nine zoom-in, cosmological simulations with non-equilibrium photochemistry involving molecular hydrogen, star formation based on the local conditions of the ISM (i.e. thermoturbulent model), radiation pressure from UV and IR photons, heating by photoionization and mechanical SN feedback. We measure the stellar mass and escape fractions from the simulations at $7 \le z < 18$ and examine the relative importance of mini-haloes for reionization by computing the evolution of the fraction of ionized hydrogen based on a simple analytic model. Our conclusions can be summarized as follows.

(i) We find that even though the instantaneous escape fraction is highly variable, the photon number-weighted average escape fraction is generally high ($\sim 20–40$ per cent) in mini-haloes (Fig. 9), in agreement with previous studies (Wise et al. 2014; Xu et al. 2016). The escape fraction is nevertheless occasionally low even after a burst of star formation, and we attribute this to the fact that in this case, the disruption of the birth clouds is too slow for LyC photons to escape efficiently (Fig. 10).

(ii) The process primarily responsible for the efficient escape of LyC photons in mini-haloes is heating due to photoionization (Fig. 11). Direct radiation pressure alone seems unlikely to create low-density channels through which ionizing radiation can escape. Rather, it simply compresses gas radially, resulting in dense neutral gas shells (Fig. 12). Because star formation is very stochastic with our thermo-turbulent model, strong radiation from young stars plays a significant role before SN explosions come into play. This leads to a short time-delay of $\lesssim 5$ Myr between the peak of the production rate of LyC photons and the escape fraction when ionizing radiation can efficiently escape into the IGM (Figs 6 and 8).

(iii) Only massive Pop III stars ($M_{PopIII} \gtrsim 100 M_{\odot}$) can play a role in ionizing the neutral gas beyond the virial radius. LyC photons from Pop III stars with lower masses ($M_{PopIII} \lesssim 70 M_{\odot}$) are mostly absorbed in the central region of their host halo (Fig. 7).

(iv) The escape fraction for individual star formation events is not strongly correlated with the star formation rate (Fig. 13). Instead, we find the production rate of ionizing photons per dense gas mass is the key to determining the escape fraction. Mini-haloes producing a large number of LyC photons per dense gas mass show a higher escape fraction. This also explains why the escape of LyC photons is significant in Pop III stars with $M_{PopIII} \gtrsim 100 M_{\odot}$.

(v) Star formation is efficiently regulated by feedback processes in the mini-haloes. The typical stellar mass in the mini-halo ranges from $10^3$ – a few times $10^4 M_{\odot}$, meaning that only less than a few per cent of baryons is converted into stars (Fig. 5). Our simulated galaxies follow the normalization and slope of the $M_{star} - M_{halo}$

sequence obtained from the radiation-hydrodynamic simulations by Kimm & Cen (2014), which is interestingly similar to the $z \sim 0$ empirical sequence by Behroozi et al. (2013) when extrapolated to the mini-halo regime.

(vi) Using a simple analytic approach, we find that the main driver of reionization is the LyC photons from atomic-cooling haloes ($10^8 M_{\odot} \lesssim M_{halo} \lesssim 10^{11} M_{\odot}$). Although mini-haloes are more abundant and their escape fraction is higher than that of the more massive haloes, their contribution to reionization is of minor importance, essentially because star formation is inefficient (Fig. 14). Even if 100 per cent of the mini-haloes are assumed to host stars, this does not alter our conclusions. Our fiducial model with a reasonably high escape fraction in the atomic-cooling regime predicts that the universe is fully ionized by $z \sim 6.7$ and $\tau_e = 0.067$. Although the optical depth estimation is consistent with the latest Planck measurement within the $2\sigma$ uncertainty, reionization is predicted to end relatively early, which may be in tension with the rather late end of reionization suggested by the most recent CMB data as well as the Ly$\alpha$ forest and Ly$\alpha$ emitter data. A larger escape of $\sim 20$ per cent in massive haloes with $10^9–10^{11} M_{\odot}$ would further raise this tension. Further investigation of the escape fraction in atomic-cooling haloes should shed light on this potential overproduction of LyC photons.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for helpful comments that improved the manuscript. We are grateful to John Wise and Renyue Cen for helpful conversations and Romain Teyssier for making his code RAMSES publicly available. We also thank John Wise and Hao Xu for sharing the occupation fraction and the Pop III star formation rate data with us. This work was supported by the ERC Advanced Grant 320596 ‘The Emergence of Structure during the Epoch of Reionization’ and the Spin(e) grants ANR-13-BS05-0005 of the French Agence Nationale de la Recherche (http://cosmicorigin.org). HK is supported by Foundation Boustany, the Isaac Newton Studentship and the Cambridge Overseas Trust. JR was funded by the European Research Council under the European Unions Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement 278594-GasAroundGalaxies, and the Marie Curie Training Network CosmoComp (PITN- GA-2009-238356). JD and AS’s research is supported by the funding from Adrian Beecroft, the Oxford Martin Schooland the STFC. This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure.

## REFERENCES

* Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93

* Ahn K., Iliev I. T., Shapiro P. R., Mellema G., Koda J., Mao Y., 2012, ApJ, 756, L16

* Anninos P., Zhang Y., Abel T., Norman M. L., 1997, New Astron., 2, 209

* Aubert D., Teyssier R., 2008, MNRAS, 387, 295

* Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376

* Baczynski C., Glover S. C. O., Klessen R. S., 2015, MNRAS, 454, 380

* Bakes E. L. O., Tielens A. G. G. M., 1994, ApJ, 427, 822

* Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023

* Becker R. H. et al., 2001, AJ, 122, 2850
Becker G. D., Bolton J. S., Madau P., Pettini M., Ryan-Weber E. V., Vene-
mans B. P., 2015, MNRAS, 447, 3402

Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

Bieri R., Dubois Y., Rosdahl J., Wagner A. Y., Silk J., Mamon G. A., 2017,
MNRAS, 464, 1854

Blaauw A., 1961, Bull. Astron. Inst. Neth., 15, 265

Blondin J. M., Wright E. B., Borkowski K. J., Reynolds S. P., 1998, ApJ,
500, 342

Bolton J. S., Haehnelt M. G., 2007, MNRAS, 382, 325

Bonazzola S., Heyvaerts J., Falgarone E., Perault M., Puget J. L., 1987,
A&A, 172, 293

Borthakur S., Heckman T. M., Leitherer C., Overzier R. A., 2014, Science,
346, 216

Boutsia K. et al., 2011, ApJ, 736, 41

Bouwens R. J. et al., 2012, ApJ, 752, L5

Bridge C. R. et al., 2010, ApJ, 720, 465

Brodie J. P., Romanowsky A. J., Strader J., Forbes D. A., 2011, AJ, 142,
199

Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Bryan G. L., Norman M. L., 1998, ApJ, 495, 80

Bunker A. J. et al., 2010, MNRAS, 409, 855

Caruana J., Bunker A. J., Wilkins S. M., Stanway E. R., Lorenzoni S., Jarvis
M. J., Ebert H., 2014, MNRAS, 443, 2831

Cen R., 2017, MNRAS, 465, L69

Cen R., Kimm T., 2015, ApJ, 801, L25

Chardin J., Haehnelt M. G., Aubert D., Puchwein E., 2015, MNRAS, 453,
2943

Choudhury T. R., Puchwein E., Haehnelt M. G., Bolton J. S., 2015, MNRAS,
452, 261

Conroy C., Kratter K. M., 2012, ApJ, 755, 123

Couchman H. M. P., Rees M. J., 1986, MNRAS, 221, 53

Cowie L. L., Barger A. J., Trouille L., 2009, ApJ, 692, 1476

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

Davies F. B., Furlanetto S. R., 2016, MNRAS, 460, 1328

Dopita M. A., Krauss L. M., Sutherland R. S., Kobayashi C., Lineweaver
C. H., 2011, Ap&SS, 335, 345

Draine B. T., Li A., 2007, ApJ, 657, 810

Draine B. T., Salpeter E. E., 1979, ApJ, 231, 77

Draine B. T. et al., 2007, ApJ, 663, 866

Dwek E., 1998, ApJ, 501, 643

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

Fan X. et al., 2001, AJ, 122, 2833

Fan X. et al., 2006, AJ, 132, 117

Federrath C., Klessen R. S., 2012, ApJ, 761, 156

Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011,
ApJ, 731, 62

Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B.,
Verner E. M., 1998, PASP, 110, 761

Finkelstein S. L., Papovich C., Giavalisco M., Reddy N. A., Ferguson H. C.,
Koekemoer A. M., Dickinson M., 2010, ApJ, 719, 1250

Finlator K., Oh S. P., Özel F., Davé R., 2012, MNRAS, 427, 2464

Fisher D. B. et al., 2014, Nature, 505, 186

Fontanot F., Cristiani S., Pfrommer C., Cupani G., Vanzella E., 2014,
MNRAS, 438, 2097

Fujii M. S., Portegies Zwart S., 2011, Science, 334, 1380

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

Giallongo E. et al., 2015, A&A, 578, A83

Gies D. R., Bolton C. T., 1986, ApJS, 61, 419

Girardi L., Bressan A., Bertelli G., Chiosi C., 2000, A&AS, 141, 371

Glover S. C. O., Federrath C., Mac Low M.-M., Klessen R. S., 2010,
MNRAS, 404, 2

Gnedin N. Y., 2000, ApJ, 535, 530

Gnedin N. Y., Kaurov A. A., 2014, ApJ, 793, 30

Gnedin N. Y., Kravtsov A. V., Chen H.-W., 2008, ApJ, 672, 765

Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke
M., Bartelmann M., 2005, ApJ, 622, 759

Grazian A. et al., 2016, A&A, 585, A48

Greif T. H., Glover S. C. O., Bromm V., Klessen R. S., 2010, ApJ, 716, 510

Greif T. H., Bromm V., Clark P. C., Glover S. C. O., Smith R. J., Klessen
R. S., Yoshida N., Springel V., 2012, MNRAS, 424, 399

Guillet T., Teyssier R., 2011, J. Comput. Phys., 230, 4756

Haardt F., Madau P., 2012, ApJ, 746, 125

Habing H. J., 1968, Bull. Astron. Inst. Neth., 19, 421

Haehnelt M. G., 1995, MNRAS, 273, 249

Haehnelt M. G., Madau P., Kudritzki R., Haardt F., 2001, ApJ, 549, L151

Hahn O., Abel T., 2011, MNRAS, 415, 2101

Haiman Z., Loeb A., 1998, ApJ, 503, 505

Halle A., Combes F., 2013, A&A, 559, A55

Heger A., Woosley S. E., 2002, ApJ, 567, 532

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

Hennebelle P., Chabrier G., 2011, ApJ, 743, L29

Hinshaw G. et al., 2013, ApJS, 208, 19

Hirano S., Hosokawa T., Yoshida N., Umeda H., Omukai K., Chiaki G.,
Yorke H. W., 2014, ApJ, 781, 60

Hollenbach D., McKee C. F., 1979, ApJS, 41, 555

Hopkins A. M., Beacom J. F., 2006, ApJ, 651, 142

Hosokawa T., Hirano S., Kuiper R., Yorke H. W., Omukai K., Yoshida N.,
2016, ApJ, 824, 119

Inoue A. K., Shimizu I., Iwata I., Tanaka M., 2014, MNRAS, 442, 1805

Izotov Y. I., Orlitová I., Schaerer D., Thuan T. X., Verhamme A., Guseva
N. G., Worseck G., 2016, Nature, 529, 178

Kashikawa N. et al., 2015, ApJ, 798, 28

Katz H., Ricotti M., 2013, MNRAS, 432, 3250

Katz H., Ricotti M., 2014, MNRAS, 444, 2377

Katz H., Kimm T., Sijacki D., Haehnelt M., 2016, preprint
(arXiv:1612.01786)

Kennicutt R. C., Jr, 1998, ApJ, 498, 541

Kim Y. et al., 2015, ApJ, 813, L35

Kim J.-h., Krumholz M. R., Wise J. H., Turk M. J., Goldbaum N. J., Abel
T., 2013, ApJ, 775, 109

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

Kimm T., Cen R., 2014, ApJ, 788, 121

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

Kirk H., Myers P. C., 2011, ApJ, 727, 64

Kroupa P., 2001, MNRAS, 322, 231

Kroupa P., Weidner C., Pf gramm-Altenburg J., Thies I., Dabringhausen J.,
Marks M., Maschberger T., 2013, in Oswalt T. D., Gilmore G., eds,
Planets, Stars and Stellar Systems, Vol. 5, Galactic Structure and Stellar
Populations. Springer Science+Business Media, Dordrecht, p. 115

Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250

Krumholz M. R., Stone J. M., Gardiner T. A., 2007, ApJ, 671, 518

Kuhlen M., Faucher-Giguère C.-A., 2012, MNRAS, 423, 862

Kulkarni G., Choudhury T. R., Puchwein E., Haehnelt M. G., 2016, MNRAS,
463, 2583

Lee H., Yoon S.-C., 2016, ApJ, 820, 135

Leethochawalit N., Jones T. A., Ellis R. S., Stark D. P., Zitrin A., 2016, ApJ,
831, 152

Leitet E., Bergvall N., Piskunov N., Andersson B.-G., 2011, A&A, 532,
A107

Leitet E., Bergvall N., Hayes M., Linné S., Zackrisson E., 2013, A&A, 553,
A106

Leitherer C. et al., 1999, ApJS, 123, 3

Leitherer C., Hernandez S., Lee J. C., Oey M. S., 2016, ApJ, 823, 64

Leonard P. J. T., Duncan M. J., 1988, AJ, 96, 222

Lopez L. A., Krumholz M. R., Bolatto A. D., Prochaska J. X., Ramirez-Ruiz
E., Castro D., 2014, ApJ, 795, 121

McGree I. D., Mesinger A., D'Odorico V., 2015, MNRAS, 447, 499

Ma X., Kasen D., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Kereš
D., Murray N., 2015, MNRAS, 453, 960
Ma X., Hopkins P. F., Kasen D., Quataert E., Faucher-Giguère C.-A., Kereš
D., Murray N., Strom A., 2016, MNRAS, 459, 3614

Madau P., Fragos T., 2016, preprint (arXiv:1606.07887)

Madau P., Haardt F., 2015, ApJ, 813, L8

Madau P., Haardt F., Rees M. J., 1999, ApJ, 514, 648

Madau P., Rees M. J., Volonteri M., Haardt F., Oh S. P., 2004, ApJ, 604, 484

Maio U., Petkova M., De Lucia G., Borgani S., 2016, MNRAS, 460, 3733

Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504

Matthee J., Sobral D., Santos S., Röttgering H., Darvish B., Mobasher B.,
2015, MNRAS, 451, 400

Meece G. R., Smith B. D., O’Shea B. W., 2014, ApJ, 783, 75

Mesinger A., Aykutalp A., Vanzella E., Pentericci L., Ferrara A., Dijkstra
M., 2015, MNRAS, 446, 566

Mirabel I. F., Dijkstra M., Laurent P., Loeb A., Pritchard J. R., 2011, A&A,
528, A149

Mostardi R. E., Shapley A. E., Steidel C. C., Trainor R. F., Reddy N. A.,
Siana B., 2015, ApJ, 810, 107

Nomoto K., Tominaga N., Umeda H., Kobayashi C., Maeda K., 2006, Nucl.
Phys. A, 777, 424

Ocvirk P. et al., 2016, MNRAS, 463, 1462

Ono Y. et al., 2012, ApJ, 744, 83

Osterbrock D. E., Ferland G. J., 2006, Astrophysics of Gaseous Nebulae
and Active Galactic Nuclei

Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2013, MNRAS, 429, L94

Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2015, MNRAS, 451,
2544

Padoan P., Nordlund Å., 2011, ApJ, 730, 40

Pawlik A. H., Schaye J., van Scherpenzeel E., 2009, MNRAS, 394, 1812

Pawlik A. H., Schaye J., Dalla Vecchia C., 2015, MNRAS, 451, 1586

Pentericci L. et al., 2014, ApJ, 793, 113

Planck Collaboration XV, 2014, A&A, 571, A15

Planck Collaboration XIII, 2016a, A&A, 594, A13

Planck Collaboration XLVI, 2016b, A&A, 596, A107

Portegies Zwart S. F., 2000, ApJ, 544, 437

Razoumov A. O., Sommer-Larsen J., 2010, ApJ, 710, 1239

Reddy N. A., Steidel C. C., Pettini M., Bogosavljevic M., Shapley A., 2016,
ApJ, 828, 108

Ricotti M., 2002, MNRAS, 336, L33

Ricotti M., Ostriker J. P., 2004, MNRAS, 352, 547

Ritter J. S., Safranek-Shrader C., Gnat O., Milosavljević M., Bromm V.,
2012, ApJ, 761, 56

Robertson B. E. et al., 2013, ApJ, 768, 71

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

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

Rutkowski M. J. et al., 2016, ApJ, 819, 81

Sana H. et al., 2012, Science, 337, 444

Schaerer D., 2002, A&A, 382, 28

Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P., 2014, ApJ, 795, 20

Schmidt M., 1959, ApJ, 129, 243

Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A,
410, 611

Shapiro P. R., Giroux M. L., 1987, ApJ, 321, L107

Shapley A. E., Steidel C. C., Strom A. L., Bogosavljević M., Reddy N. A.,
Siana B., Mostardi R. E., Rudie G. C., 2016, ApJ, 826, L24

Sheth R. K., Tormen G., 2002, MNRAS, 329, 61

Shull J. M., Harness A., Trenti M., Smith B. D., 2012, ApJ, 747, 100

Siana B. et al., 2007, ApJ, 668, 62

Siana B. et al., 2010, ApJ, 723, 241

Siana B. et al., 2015, ApJ, 804, 17

Smith B. D., Hallman E. J., Shull J. M., O’Shea B. W., 2011, ApJ, 731, 6

Smith B. D., Wise J. H., O’Shea B. W., Norman M. L., Khochfar S., 2015,
MNRAS, 452, 2822

So G. C., Norman M. L., Reynolds D. R., Wise J. H., 2014, ApJ, 789, 149

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS,
214, 15

Stanway E. R., Eldridge J. J., Becker G. D., 2016, MNRAS, 456, 485

Steidel C. C., Pettini M., Adelberger K. L., 2001, ApJ, 546, 665

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

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

Tilvi V. et al., 2014, ApJ, 794, 5

Topping M. W., Shull J. M., 2015, ApJ, 800, 97

Toro E. F., Spruce M., Speares W., 1994, Shock Waves, 4, 25

Turk M. J., Abel T., O’Shea B., 2009, Science, 325, 601

Turk M. J., Oishi J. S., Abel T., Bryan G. L., 2012, ApJ, 745, 154

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506,
647

Vanzella E. et al., 2010, ApJ, 725, 1011

Visbal E., Haiman Z., Bryan G. L., 2015, MNRAS, 453, 4456

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

Weidner C., Kroupa P., Bonnell I. A. D., 2010, MNRAS, 401, 275

Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296

Whalen D., Abel T., Norman M. L., 2004, ApJ, 610, 14

Wise J. H., Cen R., 2009, ApJ, 693, 984

Wise J. H., Abel T., Turk M. J., Norman M. L., Smith B. D., 2012a, MNRAS,
427, 311

Wise J. H., Turk M. J., Norman M. L., Abel T., 2012b, ApJ, 745, 50

Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J.,
Abel T., Smith B. D., 2014, MNRAS, 442, 2560

Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003,
ApJ, 587, 278

Woo J., Courteau S., Dekel A., 2008, MNRAS, 390, 1453

Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181

Xu H., Wise J. H., Norman M. L., Ahn K., O’Shea B. W., 2016, ApJ, 833,
84

Yajima H., Choi J.-H., Nagamine K., 2011, MNRAS, 412, 411

Yoshida N., Omukai K., Hernquist L., Abel T., 2006, ApJ, 652, 6

This paper has been typeset from a *TEX/LATex* file prepared by the author.