Fluctuating escape fraction of high-z galaxies

Fluctuating feedback-regulated escape fraction of ionizing radiation in low-mass, high-redshift galaxies


Low mass galaxies are thought to provide the bulk of the ionizing radiation necessary to reionize the Universe. The amount of photons escaping the galaxies is poorly constrained theoretically, and difficult to measure observationally. Yet it is an essential parameter of reionization models. We study in detail how ionizing radiation can leak from high redshift galaxies. For this purpose, we use a series of high resolution radiation hydrodynamics simulations, zooming on three dwarf galaxies in a cosmological context. We find that the energy and momentum input from the supernova explosions has a pivotal role in regulating the escape fraction, by disrupting dense star forming clumps, and clearing sight lines in the halo. In the absence of supernovae, photons are absorbed very locally, within the birth clouds of massive stars. We follow the time evolution of the escape fraction, and find that it can vary by more than six orders of magnitude. This explains the large scatter in the value of the escape fraction found by previous studies. This fast variability also impacts the observability of the sources of reionization: a survey even as deep as would miss about half of the underlying population of Lyman-continuum emitters.

radiative transfer – methods: numerical – galaxies: dwarfs – galaxies: formation – galaxies: high redshift – dark ages, reionization, first stars.

1 Introduction

The Epoch of Reionization (EoR) is a transition era in the history of the Universe during which the first structures formed. Until the formation of the first stars and galaxies at redshift , all the gas in the Universe was neutral. The appearance of these first objects marks the end of the so-called ‘Dark Ages’, and as the first luminous sources formed, they started to emit ionizing radiation, creating ionized H ii regions around them. As the first galaxies assembled, these bubbles gradually grew and merged, until the Universe was fully ionized at . The presence of large volumes of neutral hydrogen in the high redshift Universe has initially emerged from the observations of the Gunn-Peterson trough (Gunn & Peterson, 1965) in the spectra of high redshift quasars (Becker et al., 2001; Fan et al., 2001, 2006). It has later been confirmed through the measurement of the Thomson scattering optical depth of Cosmic Microwave Background (CMB) photons on free electrons in the intergalactic medium (IGM) with satellites like the Wilkinson Microwave Anisotropy Probe (WMAP, Hinshaw et al., 2013) and Planck. The latest results of the Planck mission report a Thomson optical depth of , equivalent to a redshift of instantaneous reionization (Planck Collaboration et al., 2016).

The nature of the sources responsible for the production of ionizing ultra-violet (UV) photons necessary to reionize of the Universe is still subject to debate. While quasars are extremely bright, they are likely too rare at to account for all the ionizing budget of reionization (e.g. Willott et al., 2010; Fontanot et al., 2012; Grissom et al., 2014; Haardt & Salvaterra, 2015). The simplest alternative is the stellar reionization scenario (see for instance Kuhlen & Faucher-Giguère 2012, or the review by Barkana & Loeb 2001), which postulates that high-redshift star-forming galaxies are accountable for the production of the bulk of the ionizing photons. With the advent of extremely deep surveys such as the Hubble Ultra Deep Field (UDF) survey (Koekemoer et al., 2013), the total production of ionizing photons from star-forming galaxies can be constrained by observations. Based on such campaigns, Robertson et al. (2013) showed that it is necessary to extrapolate the UV luminosity function (LF) of galaxies down to very faint magnitudes of typically (see also Kuhlen & Faucher-Giguère, 2012; Finkelstein et al., 2012) to simultaneously match constraints on the reionization history and the Thomson optical depth. While the recent results of the Planck mission reduced the tension between these two probes, and thus the need for early reionization, the contribution of a yet undetected population of star-forming galaxies to the ionizing budget of the Universe at might still be major if not dominant, unless the fraction of ionizing photons escaping the galaxies is very high. Similarly, using semi-analytical modelling, Mutch et al. (2016) found that galaxies residing in haloes of mass are dominant contributors of the ionizing budget of the Universe before reionization is complete. This point is strengthened by the recent observational constraints on the faint end of the UV LF at high redshift (see for example Atek et al., 2015), showing no cut-off down to magnitude , consistent with what had been found previously at lower redshift (Alavi et al., 2014). Several studies showed that these galaxies have bluer UV-continuum slopes than their more massive counterparts (Bouwens et al., 2009, 2012, 2014), with no significant redshift evolution. This indicates that there is a large reservoir of faint galaxies that must be accounted for when computing the cosmic ionizing budget (see also Bouwens et al., 2016).

There is a wealth of reionization models, which either use an analytical framework driven by observations of high- galaxies (Madau et al., 1999; Haardt & Madau, 2012; Robertson et al., 2013; Robertson et al., 2015; Duncan & Conselice, 2015), are semi-analytical (e.g Benson et al., 2006; Wyithe & Cen, 2007; Dayal et al., 2008; Mitra et al., 2015; Mutch et al., 2016) or are based on large scale simulations (e.g. Gnedin, 2014; Ricotti et al., 2002; Iliev et al., 2006; Iliev et al., 2014; Ocvirk et al., 2016), and all of them use directly or indirectly, the fraction of ionizing photons that escape the galaxy to ionize the IGM, as a key parameter. In order to quantify the contribution of low mass galaxies, one needs to assess and understand the value of , which is one of the main sources of uncertainty in the modelling of the reionization history.

Direct observations of Lyman continuum (LyC) photons leaking from galaxies have turned out to be extremely challenging. Assessing the escape fraction from individual galaxies is even more difficult, and most studies approach the problem via stacking. In the Local Universe, extremely few LyC leaking galaxies have been found, all exhibiting low escape fractions (Bergvall et al., 2006; Leitet et al., 2013; Izotov et al., 2016, see also Heckman et al. (2011)), and for some other candidates, only upper limits could be secured (Heckman et al., 2001). Borthakur et al. (2014) recently reported the finding of a discovery with , but taking dust effects into account reduced this to , further indicating that all low- LyC leakers have low . At intermediate redshifts, most efforts have given null results or relatively loose upper limits (e.g. Malkan et al., 2003; Siana et al., 2007, 2010; Cowie et al., 2009), but all report a low upper limit on . At higher , analysis of individual detections or stacked observational samples shows in general a higher (Shapley et al., 2006; Iwata et al., 2009; Vanzella et al., 2010, 2015; Nestor et al., 2013; Mostardi et al., 2013, 2015), sometimes above 50% (Vanzella et al., 2016), even though some studies found a much lower value for the typical (Boutsia et al., 2011; Grazian et al., 2016). The discrepancy might be due to the high contamination rate of the high- samples (Siana et al., 2015), and to the very uncertain properties of dust at these redshifts. Bergvall et al. (2013) suggested that selection effects might bias the low- quest for LyC leaking galaxies.

In parallel to the observational attempts to determine the escape fractions of star-forming galaxies, a strong theoretical effort has been made over the past two decades. Analytical estimates or calculations using very idealized geometries have found a fairly low value for in typical galaxies (Wood & Loeb, 2000; Clarke & Oey, 2002; Fernandez & Shull, 2011), quite sensitive to the gas distribution in the ISM (Ciardi et al., 2002). More recently, several simulations have been undertaken, either in an idealized galaxy setup or in a full cosmological context (e.g. Wise & Cen, 2009; Wise et al., 2014; Razoumov & Sommer-Larsen, 2010; Kim et al., 2013; Paardekooper et al., 2013; Kimm & Cen, 2014; Yajima et al., 2014; Ma et al., 2015; Kimm et al., 2017). These numerical experiments can give important insights into the escape of radiation in a more realistic context. However, these simulations yield vastly different results: for instance, while Gnedin et al. (2008) and Wise & Cen (2009) found that increases with the mass of the halo hosting the galaxy, other studies such as Yajima et al. (2011), Kimm & Cen (2014), Wise et al. (2014) or Paardekooper et al. (2015) found the opposite trend, and Ma et al. (2015) found no trend at all. In practice, it is very difficult to compare directly the results from these studies. The halo masses span over six orders of magnitude, and the numerical methods employed are tremendously different: Smooth Particle Hydrodynamics vs. grid codes, coupled radiative transfer vs. post-processing, various degrees of numerical resolutions, etc. Maybe even more decisively, all simulations modelling galaxies in a cosmological context must rely on sub-resolution recipes to describe star formation and the various associated feedback processes. It is extremely hard to compare directly the results of these different simulations, since even the detailed modelling of a single physical process will strongly affect the outcome of a simulation (see for instance Kimm et al. (2015) for a comparison of various implementations of supernova feedback). Furthermore, for a given simulation, there is an intrinsic scatter in the value of : Kimm & Cen (2014) found that there is a large spread in the vs. halo mass relationship, and Paardekooper et al. (2015) showed that the escape fraction seems to be affected by a large variety of physical parameters.

Using a sample of three galaxies simulated with very high resolution and state of the art subgrid models, we propose to go a step forward in studying the detailed mechanism leading to the escape of ionizing radiation, with the ultimate goal of understanding the scatter in . The combination of our subgrid models for star formation and feedback has so far only been employed in simulations of mini-haloes, which have been found to be of minor importance for reionization (Kimm et al., 2017). Our study extends the mass ranged probed by this model to haloes of masses around a few . Our new model for star formation, which causes a very bursty and clustered star formation, also extends the work of Kimm & Cen (2014), which focuses on a mass range similar to the one presented in this study. This paper is organized as follows: in Sect. 2, we describe our simulation methodology and our sub-resolution recipes. In Sect. 3, we present general properties of our simulated sample. We study in Sect. 4 the evolution of the ionizing escape fraction at the halo scale, and in Sect. 5, we discuss the mechanisms that controls the escape from the galaxy. We discuss possible improvements of this work in Sect. 7, and summarize our results in Sect. 8.

2 Description of the suite of simulations

We introduce here the suite of simulations that analyse use through this work. We focus on a selection of three haloes of masses up to that we follow down to , when the Universe was 1 Gyr old.

2.1 RHD simulations with Ramses-RT

We run cosmological simulations with Ramses-RT (Rosdahl et al., 2013; Rosdahl & Teyssier, 2015), a public multi-group radiative transfer (RT) extension of the adaptive mesh refinement (AMR) code Ramses3 (Teyssier, 2002). The evolution of the collisionless dark matter (DM) and stellar particles is followed using a particle-mesh solver with cloud-in-cell interpolation. The evolution of the gas is followed by solving the Euler equations using a second-order Godunov scheme. We use the HLLC Riemann solver (Toro et al., 1994), with a MinMod total variation diminishing scheme to reconstruct the intercell fluxes. For all the simulations, a Courant factor of 0.8 has been used. We use a standard quasi-Lagrangian refinement strategy, in which a cell is refined if , where , and are respectively the DM, gas and stellar densities in the cell, is the cell size, and is the mass of the highest resolution DM particle. In a DM-only run, this would allow refinement as soon as there are at least 8 high-resolution DM particles in a cell.

The RT module propagates the radiation emitted by stars in three groups (describing the average H i , He i and He ii ionizing photons) on the AMR grid using a first-order Godunov method with the M1 closure for the Eddington tensor. We highlight that this moment-based method can handle an arbitrary number of ionizing sources. However, because the radiation is evolved with the same timestep as the gas, we use the reduced speed of light approximation, with a reduced speed of light of , where is the real speed of light. The ionizing photons are produced by star particles at each timestep using the Galaxev model of Bruzual & Charlot (2003), assuming a Chabrier (2003) initial mass function (IMF). Note that while most of the radiation is emitted by stars younger than 5 Myr, we continue to take into account the photons emitted by older star particles. For this set of simulations, we use the on-the-spot approximation, assuming that any ionizing photon emitted by recombination will be absorbed locally. We hence ignore straight-to-ground level recombination and the associated emission of ionizing radiation from the gas. The coupling with the hydrodynamical evolution of the gas is done by incorporating the local radiation while computing the non-equilibrium thermochemistry for the hydrogen and helium.

2.2 Initial conditions

We assume a flat CDM cosmology consistent with the Planck results ( and , Planck Collaboration et al., 2014). We select three haloes from a comoving Mpc DM-only simulation with particles () initialized at . The haloes are selected in the final output at with HaloMaker (Tweed et al., 2009). They all are relatively isolated (there is no halo with more than half its mass within more than 10 virial radii), with quiet merger histories and masses of approximately and at .

The haloes are then resimulated for 1 Gyr at much higher resolution using the ‘zoom-in’ technique. The initial conditions for the resimulations are produced using the Music4 code (Hahn & Abel, 2011). For each resimulation, we refine the mass distribution such that the DM particle mass is in the zoomed-in region, equivalent to a particle simulation. This corresponds to an AMR coarse grid level of in the zoomed-in region, and we allow for 9 more levels of refinement, giving a most refined cell size of pc.

Since we do not track the formation of population III (pop III) stars and early metal enrichment, we assume an initial gas phase metallicity everywhere, consistent with metal enrichment from primordial mini-haloes (Whalen et al., 2008).

2.3 Physical modelling

Star formation

Because of the huge dynamical range involved in galaxy and star formation, we cannot afford to resolve all the stages of the gravitational collapse from the cosmological scale down to the formation of individual stars. Following what is usually done in galaxy formation studies, we use a subgrid recipe to model the collapse of a gas cloud below the resolution limit. In this section, we briefly describe the main ideas behind our star formation model. The reader interested in the detailed implementation is referred to Devriendt et al. (in prep, see also Kimm et al. 2017) for a complete discussion of the model.

The main physical ingredient of this model is the idea that at the star-forming cloud scale, turbulence in the ISM can act as an effective extra pressure support. We thus define a ‘turbulent Jeans length’ as


with the velocity dispersion of the gas computed using the velocity gradients around the cell, is the gas density and is the local sound velocity. We identify star forming sites as cells with , i.e. cells which are unstable even with this additional turbulent support.

For these star forming cells, we use an approach similar to that of Rasera & Teyssier (2006) to convert of gas into star particles during one timestep , where is the local star formation efficiency and is the gas free fall time. The primary difference with the recipe of Rasera & Teyssier (2006) is that the star formation efficiency is not a constant, but rather derived from the local properties of the gas (Hennebelle & Chabrier, 2011; Federrath & Klessen, 2012):


where characterizes the turbulent density fluctuations, is the critical density above which the gas will be accreted onto stars (where is the mean density of the cloud), and . The exact formulae are given in Federrath & Klessen (2012), using their multi-scale model based on Padoan & Nordlund (2011), coined ‘multi-ff PN’ model. As a result, instead of being smoothly distributed in the ISM, the star forming gas will be gathered in clumps, leading to a very clustered star formation process. We adopt a minimum stellar particle mass of , leading to the explosion of at least one supernova per particle assuming a Chabrier (2003) IMF. On average, for the most massive galaxy in our sample, the mean mass for the star particles is around .

Feedback from massive stars

Our simulations implement both radiative feedback from stars and type II supernova feedback. For the radiative feedback, we only take photoionization heating into account, since it is believed to be the dominant channel of radiative feedback (Rosdahl et al., 2015). More specifically, star particles emit ionizing radiation as described in 2.1, which will inject energy and heat the surrounding medium.

We use the recent mechanical supernova (SN) feedback implementation described in (Kimm & Cen, 2014; Kimm et al., 2015, see also Rosdahl et al. 2017), with one instantaneous supernova event per star particle after a 10 Myr delay, consistent with the typical lifetime of a SNii progenitor of . The main motivation for this model is to describe correctly the momentum transfer from the SN at all stages of the Sedov blast wave, from the early adiabatic expansion to the late snowplough phase. In practice, the amount of momentum deposited in the surrounding cells depends on the properties (density and metallicity) of the gas in these cells. Note that because the gas distribution can be anisotropic around a star particle, the SN explosion can in principle result in an anisotropic energy deposition as well. A feature of our new star formation model is that the stars will form in a very clustered fashion, and as a result the total amount of energy deposited by the supernovae from a group of star particles is very similar to what would happen with more massive star particles.

Gas cooling and heating

Ramses-RT features non-equilibrium cooling by explicitly tracking the ionization state of 5 species (H, H, He, He, He). The primordial cooling rate is directly computed from these abundances. We include an extra cooling term for metals, using tabulated rates computed with Cloudy (last described in Ferland et al., 2013) above  K. We also account for energy losses via metal line cooling below  K following the prescription of Rosen & Bregman (1995). We approximate the effect of the metallicity by scaling linearly the metal cooling enhancement. Note that we do not take into account the impact of the local ionizing flux on the metal cooling. We use an homogeneous metallicity floor of in the whole box to account for the lack of molecular cooling and to allow the gas to cool down below  K. We ensured after the fact that our choice of initial metallicity leads to a redshift of first star formation in the progenitors of our main haloes roughly similar to that found by studies of minihaloes including H cooling (e.g. Kimm et al., 2017). Since we do not include any meta-galactic UV background in our simulation either, we only take into account the local photoheating rate. While this will not be valid at the end of the Epoch of Reionization, when the ionized bubbles are overlapping, this holds for our isolated galaxies that are likely to reionize themselves.

2.4 Computing the escape fraction

A last methodological point to discuss is the way we measure the instantaneous escape fraction. With Ramses-RT, we have access to the local photon flux anywhere at any time. We can then simply compute the escape fraction as the total flux crossing a given boundary divided by the intrinsic luminosity of the sources. Unless stated otherwise, we always compute the escape fraction at the virial radius.

Figure 1: Mass-weighted projection of the gas density (top row) and zoomed stellar surface density (bottom row), as indicated by the black boxes, for the three galaxies at the end of each simulation, ordered by decreasing mass from left to right: and . The black circle shows the virial radius of each halo, while the physical scale and the redshift of the snapshot are indicated in the bottom left and right corners, respectively.

However, since the speed of light is finite (and even considerably reduced in our simulations), it takes time for photons to propagate from their source to the point at which we compute the escape fraction. For a steady source of photons, this delay would be unimportant, but as we will see in Sec. 3, the star formation is very bursty in our simulations. This means that the galaxy will flicker. A photon crosses 10 kpc in approximately 30 kyr at the full speed of light. At our reduced speed of light, this timescale is increased by a factor 100, meaning that in our simulation, the radiation emitted in the centre of a halo with kpc would only reach the virial radius after 3 Myr, which is comparable to the lifetime of massive stars, and of the same order of magnitude as the duration of the bursts. It is therefore necessary to take this delay into account when we compute the escape fraction.

We compute the angle averaged escape fraction at a distance from the centre of the halo at any time as


where is the outgoing flux of photons, is the local radial direction from the centre of the halo, is the mass of a star particle in unit mass, and is the ionizing photon production rate per unit mass for a simple stellar population of age . We take if the local flux of photons is indeed going outward and 0 otherwise, which prevents us from computing negative contributions to the escape fraction coming e.g. from a satellite or a star-forming clump infalling onto the main galaxy. This is making the assumption that all the photon sources are located at the centre of the halo, leading to a common time delay for all stars.

While this is not always valid in theory (e.g. a star-forming clump closer than the others to the virial radius), we will check in Sect. 4.2 whether this is a good approximation. For this purpose, we introduce a second way of estimating the escape fraction, through ray-tracing. We generate directions using the HEALPix (Hierarchical Equal Area isoLatitude Pixelation, Górski et al., 2005) decomposition of the sphere into equal-area pixels. We typically use , or directions, unless stated otherwise. For each direction , we can compute the optical depth at a given distance for each star particle , where is the average cross-section for the H i -ionizing radiation, and is the H i column density seen from the star in the direction . As , any line of sight with would be optically thick to ionizing photons. The luminosity-weighted angle-averaged escape fraction can then be expressed as


with the angle-averaged transmission for the star particle, and its ionizing luminosity. While this does not take the halo light crossing time into account, it is in principle more robust since it does not rely on the assumption that all sources are at the centre of the halo. This can also be employed to estimate the individual escape fraction of each star particle. We will discuss in Sect. 4.2 the differences between these two methods on the estimated value of .

3 Galaxy properties

Figure 2: Left: Stellar mass - halo mass relationship for the three haloes. Each halo is indicated by a different symbol: filled circles ( halo), squares ( halo) and triangles ( halo), the colour of the which indicates the redshift. The three diagonal dotted lines indicate constant baryonic fractions 10%, 1% and 0.% from top to bottom. At low mass, episodes of star formation are followed by long periods of quiet accretion. Roughly between 1% and 10% of the baryons in each halo is converted into stars. Right: Relationship between stellar mass and gas mass within . After each star formation episode, most of the gas is ejected out of the galaxy.

In Fig. 1, we illustrate from left to right the gas density and stellar surface density of the three galaxies at the end of each simulation, in decreasing order of halo mass. The stellar distribution is very clumpy and irregular, as expected for such low mass galaxies. While all three haloes seem to be embedded in large scale filaments, we see at least for the first two haloes a disturbed gas morphology, shaped by the strong SN feedback affecting these galaxies.

In Fig. 2, we show the mass assembly history of each of the three simulated haloes by comparing the stellar mass to the total halo mass (on the left) or to the gas mass within 10% of the virial radius (on the right). The coloured symbols represent the positions of the haloes on the - plane at a given redshift. The diagonal dotted lines indicate, from top to bottom, 10%, 1% and 0.1% of the baryonic mass fraction. Overall, in our simulations, roughly between 1% and 10% of the baryons are locked in stars at any given time. This is slightly higher than what has been found by previous studies, but given the lack of observational constraints at for low-mass galaxies, the large scatter in the - relation allows for haloes such as the one we present here. The recent work of Miller et al. (2014, Fig. 6) shows that in the Local Volume, low-mass galaxies have higher stellar mass to halo mass ratios than predicted by abundance matching. While these galaxies are not exact analogues of higher ones, this is consistent with our findings, and despite the large uncertainties, this mitigates the importance of differences between our simulations and other results e.g. from abundance matching techniques.

We note two striking features of Fig. 2. First, for all three galaxies, the first episodes of star formation occur at (or similarly at ), independent of the redshift, which is the typical ‘atomic cooling’ limit. This is most likely a resolution issue: we re-ran the lowest mass halo at a higher resolution, and found that the first stars formed earlier. Second, we see that the stellar mass increases only episodically, especially at low masses. As an example, let us focus on the ‘large’ halo. At , a first dramatic star formation event happens, pushing the stellar mass of the galaxy to . This episode is followed by a long plateau during which star formation is shut off as the halo still grows in mass, until (approximately 120 Myr later). This discontinuous stellar mass assembly directly results from the supernova feedback, which strongly reduces (or even shuts off) star formation by ejecting most of the gas out of the galaxy, so the galactic gas reservoir has to replenish before fuelling any new star formation (right panel). At the same time, the growth of the halo mass (and also the total gas mass) is virtually unaffected.

Figure 3: Star formation (blue, solid line), outflow at (orange, dotted line) and infall at (green, dash-dotted line) history of the most massive halo as a function of time. Each star formation episode is followed by a massive outflow event.

The bursty behaviour of the star formation is further illustrated as a blue filled surface in Fig. 3 for the most massive halo, where we show its evolution as a function of time. Here, we define the star formation rate as the mass of stars newly formed averaged over 10 Myr. We also display the outflow rate at the virial radius with the orange dash-dotted line. Typically, during a star formation episode, the galaxy can reach a SFR as high , and more typically , each episode lasting for 20 to 50 Myr. This burstiness has been seen in other numerical studies (e.g. Wise & Cen, 2009; Kimm & Cen, 2014; Hopkins et al., 2014; Shen et al., 2014), and is regulated by supernovae explosions that heat the gas or even eject it, and delay the formation of new dense, star-forming clumps. Indeed, we see that approximately 10 Myr after the beginning of each star formation episode, the outflow rate rises up to 2 to 4 , denoting the presence of strong winds, with mass loading factor above unity. The evolution of the outflow rate with time follows the same highly variable pattern as the SFR: because star formation is shut off, after all the massive stars have ended their lives as supernovae, the galactic winds gradually weaken.

These winds are responsible for the presence of the plateaus in in Fig. 2: the strong outflows quench star formation for some time. For the ‘medium’ halo, this lasts for about 200 Myr, while for the ‘large’ halo, the plateau only lasts for 100 Myr. One reason for this is that the infall rate at the virial radius (the green dotted line in Fig. 3) is typically higher by a factor of 2 to 3 for the large halo compared to the medium one. In the large halo, the time needed to refill the gas reservoir after the massive outflow is therefore significantly shorter, and star formation therefore resumes faster.

4 Halo escape fraction

We now discuss the ionizing properties of the galaxies in our sample with the goal of understanding the large galaxy to galaxy variation in the escape fraction that has been found in other studies.

4.1 Time evolution of the escape fraction

Figure 4: Evolution of the escape fraction of ionizing radiation with time for all three haloes. The upper panel is in linear scale, while the lower panel is in log scale. All haloes display a very bursty behaviour, varying from 0 to sometimes more than 60% in a very short time. The smallest halo only reaches during a short burst at Myr and therefore remains indistinguishable from 0 in the upper panel.

In Fig. 4, we show the evolution of the angle averaged escape fraction of H i -ionizing radiation for the three haloes, which has a very bursty behaviour. The most massive halo of our study (in red) experiences six episodes during which , each of them lasting for 10 to 50 Myr. These episodes are followed by periods of much lower () escape fraction. The medium halo (in blue) experiences a similar alternation of high and low escape fraction, with higher than 20% for long periods. For the smallest halo (in green), there is only one peak of escape fraction, around , but the galaxy only reaches a relatively low escape fraction of . This is because the galaxy undergoes only one star formation event. Except for that small halo, our values are broadly consistent with the findings of Kimm & Cen (2014) and Paardekooper et al. (2015) who find that for haloes of , is typically between and , but can be as high as . We must however note that these studies present quantities averaged over an ensemble of galaxies (from a few hundreds for Kimm & Cen 2014 to a few ten thousands for Paardekooper et al. 2015), and the intrinsically stochastic nature of the very strong variations in prevents us from comparing a time-averaged value to the ensemble average of the aforementioned studies.

Figure 5: Evolution of the star formation rate (blue), outflow rate at (dash-dotted, orange) and escape fraction (red) for the most massive halo. The escape fraction starts to rise at the same time as the outflow rate, and typically 10 Myr after the beginning of a star formation event.

The quick variations over several orders of magnitude of the escape fraction reflect the burstiness of the star formation histories of the three galaxies that we discussed in Sect. 3. We illustrate the correlations between these two quantities by plotting in Fig. 5 the escape fraction (red), the star formation rate (blue) and the outflow rate (dash-dotted orange) for the most massive halo in our sample. Associated with the bursts of star formation, we see that the escape fraction systematically jumps to its highest value right at the onset on the wind. This sharp transition indicates that as soon as there is a hole in the ISM, radiation is able to escape. This typically happens with a 10 Myr delay with respect to the beginning of a new event of star formation, corresponding to the age at which star particles experience the supernova explosion with our modelling (see Sect. 3). While this time delay of Myr has been used in other studies (e.g. Kimm & Cen, 2014; Wise et al., 2014; Ma et al., 2015), we must tread carefully in that direction, and not take this number at face value: this is the direct result of our subgrid recipe for supernovae, which explode 10 Myr after the birth of the star particle. We also note that the peak of is usually reached much more rapidly than the maximum outflow rate: this is because while we compute both quantities at the virial radius, radiation travels this distance much faster than the gas carried in outflows. In addition to limiting the star formation, the supernova explosions completely alter the morphology of the gas in and around the galaxies. The combination of outflows and shock-heating of the gas clears lines of sight around the galaxy and allow ionizing radiation to freely flow into the IGM.

Fig. 6 presents the instantaneous, intrinsic ionizing luminosity of the most massive galaxy () in red, and the remaining luminosity that escapes the halo () in blue. The general trend is the same as for the evolution of the escape fraction: both the emitted and the escaped luminosity vary quickly over several orders of magnitude. However, we note that the phases during which (high escape fraction) do not necessarily correspond to peaks of the emissivity of the galaxy. For instance, right before Gyr, reaches 15% and the photon production rate is roughly ionizing photons per second. The galaxy produces a similar amounts of photons around 920 Myr, but there barely reaches 0.1%. The maximum injection of ionizing photons into the IGM does not necessarily correspond to the peaks of star formation (and thus of intrinsic luminosity). Therefore, a galaxy experiencing a strong star formation episode will not always contribute significantly to the ionizing budget of reionization.

Figure 6: Photons emitted (blue) and escaping into the IGM (red) per second from the most massive halo. broadly follows , but can drop significantly lower. Not all peaks in have corresponding peaks in , meaning that not all episodes of star formation contribute to reionization.

Overall, if we assume that our simulated sample is reasonably representative of the low mass, high- galaxy population, we can expect that the fast paced variations of during the galaxies history result in a large scatter for the escape fraction for a homogeneous population of galaxies with e.g. the same mass. Indeed, the exact value of is highly dependent on the phase (pre-starburst, starburst or post-starburst) the galaxy is going through. Even at fixed , a galaxy population will likely be out of phase, leading to a very different from galaxy to galaxy.

4.2 Escape fraction estimator

We introduced in Sect. 2.4 a second estimator for the escape fraction using ray-tracing in order to test the robustness of our measurements of the escape fraction. We present in Fig. 7 the angle averaged, luminosity-weighted escape fraction computed using eq. 4 for the most massive halo. The red filled area represents the ray-tracing estimator of , while the solid black line is the flux estimator used in the previous figures. We plot both in linear (upper panel) and logarithmic scale (lower panel) to better illustrate the amplitude of the variations.

We could in principle expect some differences : the flux-based estimator assumes that all stars are at the centre of the halo, and may also give imperfect estimates as the moments method we use for the radiative hydrodynamics is known for being too diffusive (Iliev et al., 2006, 2009; Rosdahl et al., 2013), leading to spurious leakage of radiation around opaque obstacles that should in principle cast a shadow. As we can see, this is not the case, and it is very reassuring that the two estimators are in very good agreement: the variations are seen at the same times and have very similar amplitudes.

We explain this by the fact that the stars are concentrated at the centre of the halo. Seen from the virial radius, the behaviour of the galaxy is very close to that of a central point source. The effect of the ray-crossing issue of M1 methods (i.e. when two opposed collimated beams would collide instead of cross each others) seems to be marginal, again because of the light sources are numerous and all very central.

Figure 7: Comparison of the two estimators of : based on ray casting in red, and based on the local ionizing flux in black. The two estimators agree very well.

4.3 Directionality

Figure 8: Full-sky Mollweide projection of the gas flows (upper panel) and of the escape fraction (lower panel) at from the centre of the most massive halo at , when . In the upper panel, the gas mass flux is positive for outflowing gas and negative for infall. Photons mostly escape through channels carved by SN feedback.

In the upper panel of Fig. 8, we show the angular distribution of gas flows reaching the virial radius for a snapshot at , when . The wind seems to develop in mostly three directions (because of the Mollweide projection, the rightmost and leftmost patches look like they are disconnected), with a large part of the sky blocked by infalling gas. On the lower panel, the angular distribution of is displayed using the same projection. It shows clearly that radiation escapes the halo through channels created by the strong winds: the high- patches follow the same morphology as the outflows.

The importance of those high-outflow channels has been already studied e.g. by Fujita et al. (2003), who found that the formation of outflows is a necessary condition for the creation of low column density direction through which ionizing radiation can escape the galaxy. This is in line with the more recent findings of Gnedin et al. (2008), Wise & Cen (2009), Kim et al. (2013) and Paardekooper et al. (2015) who showed that the ionizing radiation escapes anisotropically, favouring low column density regions resulting from such outflows.

4.4 Contribution to the ionizing budget

As we have discussed in Sect. 4.1, at a given stellar mass, the escape fraction varies a lot and is not in phase with the evolution of the star formation rate. This results in situations in which a galaxy can have a high and still contribute only very little to the ionization of the IGM. In order to make a step forward in quantifying the actual contribution of small galaxies to the ionizing budget of the Universe, we compute the ionizing duty cycle of our galaxies. We define the duty cycle as the fraction of the time spent by all galaxies in a mass bin with .


where is the time spent by the galaxy in the mass bin , and is the time spent in the mass bin with . We define similarly .

Figure 9: Ionizing duty cycle of the two most massive haloes, as a function of the galaxy (upper panel) and halo (lower panel) mass. The blue (red) dashed line shows the fraction of the time a galaxy of a given mass spends with higher than () photons per second. On average, a more massive galaxy spends more time leaking Lyman continuum photons.

We show these quantities in the upper panel of Fig. 9, for four regularly log-spaced mass bins, centred on and , with in red, in blue. There is a clear trend that more massive galaxies tend to spend more time in a Lyman-leaking phase. At the very low-mass end of our plot, the duty cycle drops to zero: this corresponds to the first stellar population, formed in a single short episode, lasting around 20 Myr (e.g. see Fig. 3). When the first massive, Lyman-continuum producing stars end their lives, the production of ionizing photons drop rapidly, and is therefore quite low when the supernovae start to create clear channels. At these early times, the amount of photons reaching the IGM before the supernova phase is low as well. Similarly, the lower panel of Fig. 9 shows the ionizing duty cycle as a function of halo mass, and , for four regularly log-spaced mass bins centred on and . This is in strong disagreement with the recent results of Kimm et al. (2017), who simulated the formation of mini-haloes and found that can be very high even before the first supernova. However, their study focuses on much lower mass haloes, reaching at most. They also include a description of pop III stars and molecular hydrogen, which can lead to significant differences for haloes below the atomic cooling limit of .

On the other side of the mass spectrum, when , the duty cycle saturates, with and . This happens when the galaxy starts to reach a more stable regime, with feedback-regulated star formation, but where the star formation is rarely completely shut off. This means that the gas in the galaxy will usually not have the time to settle down in between two starburst events, leading on average to a more steady production of ionizing photons.

5 Local escape of ionizing radiation

Figure 10: Top row: Temperature maps of the largest galaxy of our sample before (left), at the beginning (middle) and a few million years after (right) a massive feedback event. Bottom row: Luminosity-weighted distribution of the Strömgren radius of each stellar particle (see text), as a function of the age of the particle. The thick horizontal black line indicates the size of the most refined cell. Areas below the line correspond to radiation absorbed inside the emission cell. Each pixel is colour-coded by the total ionizing luminosity.

In Sect. 4, we have found that the escape fraction is mainly regulated by the interplay between supernovae, gas accretion, and clustered star formation, and we have limited our analysis to the total escape fraction of each galaxy, computed at the virial radius. Because it is strongly linked with ISM-scale phenomena, we will now concentrate on the intra-halo processes leading to the escape of ionizing radiation. For this purpose, we will mainly focus on the most massive galaxy in our sample.

5.1 Does radiation escape from the emission clouds?

Before reaching the IGM, the ionizing photons must travel through the complex distribution of gas inside and around the galaxy. Kimm & Cen (2014); Ma et al. (2015); Paardekooper et al. (2015) showed that when is low, most of the photons are absorbed very locally, within 100 pc of their emission point, and that it is indeed crucial to resolve properly the ISM in order to study properly the escape of ionizing radiation. We find very similar results with our simulations, and we suggest that this is the reason behind the strong coupling between the supernova feedback and evolution of .

While the idea that clear channels for ionizing radiation from young stars is created by ionizing radiation itself is appealing, Geen et al. (2015b) showed that depending on the structure of the cloud and the strength of the source, the H ii region may or may not expand beyond the boundaries of the cloud. We note that even with our high resolution, we cannot properly resolve the internal structure of these H ii regions. Nevertheless, we can still make a first step in the analysis of the local escape of photons from the star forming clumps. If we define the Strömgren radius5 of an ionizing source as the radius of the sphere within which the rate of recombination is balanced by the ionizing luminosity, we can express as


where is the rate of ionizing photons emitted by the source, is the density, and is the case B recombination rate, given by Hui & Gnedin (1997). A young star particle (before 3-4 Myr) of 135 (our mass resolution) will yield approximately ionizing photons per second. For a typical diffuse ISM density of cm, this gives pc, much larger than our typical cell size of 7 pc. However, inside a star forming cloud, the density can reach cm in our simulations, resulting in pc, meaning that the radiation should all be absorbed inside the emission cell.

While this depends on the inner structure of the cloud, this should in principle prevent the apparition of high escape fraction episodes for most star forming clouds. We argue that the leakage of ionizing radiation from such dense clouds is the result of supernovae in the neighbourhood of the young star cluster. The explosion removes gas from the cluster, effectively lowering the local gas density and thus making it possible for radiation to escape before the end of the few Myr lifetime of the massive stars. In Fig. 10, we compare three consecutive snapshots: right before a massive feedback event, right after, and when the outflow reaches the virial radius. The upper panel shows the temperature maps for the three snapshots, where the outflow can be traced by the expansion of the hot (dark red) region. In the lower panel, we present the corresponding luminosity-weighted distribution of Strömgren radii for the star particles as a function of their age. Note that this assumes that each star particle is alone in its cell. Because stars form in clump, this will underestimate the effective Strömgren radius of each star. Before the onset of supernovae, we see that most of the ionizing radiation is produced by young stars embedded in dense cells (the red region in the bottom left corner of the plot), resulting in much smaller than the cell size. In the middle column, right after the onset of supernovae, the ionizing emissivity is dominated by young stars for which is several thousands time larger than the cell size (the dark red spot at 5 Myr and ). While the value of is approximate (e.g. because the ISM is not homogeneous), this means that a lot of radiation will escape from the galaxy. The rightmost column shows the expansion of the outflow: the stellar population is ageing, and while radiation continues to leak (meaning a high ), the net ionizing emissivity is much lower. We can match that sequence to Fig. 5 and 6: at 850 Myr, there is noticeable burst of star formation, followed by the development of a large outflow. At the beginning, the intrinsic ionizing emissivity is high, but is very low. Then, just as the outflow rate starts to increase, the escaped emissivity rises, and . In the final stages of the succession of events, the star formation rate drops and there are no more young massive stars left. Both the intrinsic and escaped emissivity drop, so remains high, of the order of 50%.

Figure 11: Evolution of for the medium halo with SN feedback and a time delay of 10 Myr (3 Myr) between star formation and the supernova in blue (red), and without feedback in black.

As an additional line of evidence that SN feedback is a crucial element of the journey of ionizing photons from the stars to the IGM, we present in Fig. 11 a comparison of for three simulations of the intermediate mass halo with variations on the feedback. The blue curve shows the fiducial simulation presented before, the red curve corresponds to a run where the delay between star formation and the supernova phase has been reduced to 3 Myr, and the black curve corresponds to a case without supernovae. Interestingly, the total mass of stars formed in the simulations with feedback is very close, with less than 15% difference after 1 Gyr. For the two simulations with feedback, the evolution of is qualitatively similar: rapidly alternating, and reaching values as high as 30% to 50%. Even if it is difficult to quantitatively assess the influence of the delay because of the randomness of the starburst events, these two runs present behaviours strongly contrasting with the no-feedback run, for which almost no radiation escape. In that last case, even though much more stars are formed over the course of the simulation, is always below 1%. This indicates that, at least at the pc resolution of our simulations, radiation feedback from the young stars in not strong enough to ionize the star-forming cloud, and thus that it is SN feedback that permits radiation to escape.

Figure 12: Luminosity contributed by stars in a clear environment (filled green) compared to the total emitted luminosity (blue). The escaped luminosity is indicated by a solid black line and is closer to green curve than to the total emitted luminosity.

This analysis indicates that the escape of ionizing radiation is mostly a local phenomenon, controlled by whether or not the photons can escape the cloud in which they have been emitted. Since the typical size for molecular clouds is of the order of a few tens of parsecs, we can consider a star particle to be in a clear environment if pc. In Fig. 12, we compare the total intrinsic ionizing luminosity (in blue) to that of stars in such a clear environment (in filled green). We show the escaped luminosity as a red line. Overall, high escape fraction episodes happen in phases where the ionizing luminosity is dominated by stars in a clear environment. We can note that in most cases, the evolution of the luminosity contributed by stars in a clear environment follows closely that of the escaped luminosity, strengthening our scenario in which the escape of ionizing radiation is mainly regulated on the cloud scale.

Figure 13: Distribution of the luminosity-weighted individual escape fraction from young stars for a snapshot where the global is low on the upper panel, and high on the lower panel. The thick lines denote the individual escape fraction at the virial radius, and the thin lines show the local escape fraction at 100 pc.

We continue this analysis in Fig. 13, where we compare the escape fraction on a local scale (filled areas) and at the virial radius (thick solid lines), for a snapshot in which the galaxy is leaking ionizing radiation (lower panel) and for another where no radiation escapes the halo (upper panel). For both panels, has been estimated using the ray-tracing method described in Sect. 2.4, and we plot the luminosity-weighted distribution of the transmissivity per star particle averaged over all directions , defined as , where the optical depth is integrated over both 100 pc and . As depicted on the upper panel, it is clear that when no radiation escapes the halo, it is because all the photons are absorbed locally, within 100 pc of their emission site. When the cloud has been destroyed and all the photons escape locally (lower panel), most of them will be able to reach the IGM. This is very consistent with previous results of e.g. Kimm & Cen (2014); Ma et al. (2015); Paardekooper et al. (2015).

5.2 Absorption within the halo and in the CGM

Simulations such as the ones we present in this work can be used to calibrate large scale models of reionization, or analytical estimates of the history of reionization. The key figure in these models is usually the angle averaged halo escape fraction, which represent the amount of ionizing photons that will escape the DM halo they originate from. This is typically the view that is adopted in semi-analytical models. In the previous sections, we have argued that the escape of ionizing radiation is mostly regulated on a local scale by various stellar feedback processes. A fraction, however, is still absorbed within the halo. This can be either because of intervening clouds or just because the halo is not fully transparent. For example, in the case of the snapshot corresponding to the the lower panel of Fig. 13, about of the hydrogen in the halo is neutral, which be responsible for the absorption in the halo. The LyC leaking episodes typically correspond to the development of strong galactic winds, that significantly alter the gas distribution inside and around the galaxy, especially in low mass galaxies where SN feedback is thought to be the most relevant (see e.g. Teyssier et al., 2013). These winds can reach several virial radii, somewhat challenging the implicit assumption that the halo is a black box which radiative efficiency can be described with only the parameter. We attempted to use radial profiles to determine whether or not the virial radius is characteristic of the escape of radiation (e.g. does reach a plateau at the virial radius?), and it appears that is probably not a relevant distance for photon propagation. This is not surprising: we have shown that the escape of radiation is a local process, and gaseous haloes around galaxies extend beyond anyway.

Figure 14: Upper panel: Comparison of the typical absorption distance (in red) to (in blue). When , most of the photons can escape the halo and reach the IGM. Lower panel: Distribution of the angle-averaged around individual star particles for the two snapshots of Fig. 13.

From the viewpoint of the ionizing sources, it makes more sense to elucidate at which distance photons are absorbed. With this in mind, we employ the same ray-tracing technique we described in Sect. 2.4 with 192 rays per star particle, but instead of computing the optical depth at a given distance, we evaluate the distance at which the optical depth reaches . We compare in the upper panel of Fig. 14 the luminosity-weighted average (in red) to the virial radius of the halo (in blue). The lower panel presents a comparison of the distribution of this characteristic scale for the two snapshots already illustrated in Fig. 13, employing a higher resolution (we used 3072 rays per particle, corresponding to a HEALPix level of 4). We note that the outputs where the average is high correspond to the LyC leaking episodes, and ultimately reaching the peak of when . During these phases, there is very little absorption within the halo. The blue histogram corresponds to a snapshot at Myr, where reaches its highest value, around 75%. The photons emitted by almost every star particle travel more than ten times farther than the virial radius before being absorbed. On the contrary, the red histogram corresponds to an episode of very low (at Myr), and on average, the radiation is absorbed within less than 0.01 pc, well inside the emission cell. Looking again at the upper panel, we see that for most of the time, the typical absorption distance reaches 10 – 1000 pc. This corroborates our findings that channels created by SN explosions favour the escape of radiation. Unless the galaxy undergoes a massive, coordinated feedback event, there will still be absorption within the halo, even outside of the emission cloud.

6 Observability of small galaxies

Due to the opacity of the high redshift IGM to ionizing radiation, very deep surveys even with the next generation of telescopes such as the James Webb Space Telescope (JWST) will not be able to directly measure the ionizing flux coming from the sources of reionization. It is therefore necessary to assess the non-ionizing properties of such galaxies at frequencies that will be observed by these instruments. For this purpose, we compute the non-ionizing, rest-frame UV magnitude around for the three galaxies of our sample. The luminosity of each star particle is given by the models of Bruzual & Charlot (2003) and is a function of its age and metallicity, but we do not include any dust absorption. To compute the luminosity of the galaxy, we sum the luminosities of each star particle within the virial radius. We present in Fig. 15 the evolution of the absolute UV magnitude expressed in the AB magnitude system (Oke & Gunn, 1983).

Figure 15: Evolution of the UV magnitude at for the three simulated galaxies. The variability follows that of the star formation rate, albeit on longer timescales.

Overall, the peaks of UV luminosity broadly follow the episodes of star formation, with a much slower decline. We find that the most massive galaxy can reach absolute UV magnitudes as high as , but spends most of its time at . This lower limit is comparable to the deepest surveys using gravitational lensing to probe the very faint end of the luminosity function at high redshift (e.g. Atek et al., 2015), and will be within the reach of the next generation of surveys using JWST as proposed e.g. by Finkelstein et al. (2015). Recently, Huang et al. (2016) found a Lyman- emitting galaxy at with a detectable UV continuum around at . They inferred a mass of and a star formation rate of . This is very similar to our most massive galaxy in its bursting phase, hinting that we are already starting to observe galaxies just like the ones presented in this work.

Figure 16: Comparison of the UV luminosity and the (escaped) ionizing luminosity for each snapshot of the three simulations. The redshift of the snapshot is indicated by the colour coding. The apparent lack of correlation is explained the much faster decrease of ionizing luminosity compared to non-ionizing UV.

However, Fig. 15 shows that there is a large variability in the UV magnitude, and we have seen in Sect. 4 that the escape of ionizing radiation is also highly variable following the star formation episodes. The correlation between the escaped ionizing luminosity and the UV magnitude is relatively poor, as is illustrated in Fig. 16. The figure compares the UV magnitude and for each snapshot of all three simulations, the colour coding indicates the redshift of the snapshot. We see that at any time when a galaxy shines at in the UV, its ionizing luminosity spans over more than six orders of magnitude, ranging between and ionizing photons per second. This lack of correlation can be explained by two factors: first, there is a delay between the peak of ionizing emissivity and escape of radiation due to the necessity to pierce the ISM, and second, the ionizing luminosity decreases much faster than the UV luminosity as a function of stellar age. This results in episodes with star formation (meaning high and intrinsic ionizing emissivity) but very small , leading to a low .

Figure 17: Estimation of the completeness of a sample limited by UV magnitude. The lines show the fraction of time that a galaxy emitting a given ionizing luminosity spends at a magnitude above -12, -14 or -16 (respectively in blue, green or red).

Conversely, this has consequences for the selection of Lyman-leakers: a galaxy releasing at least ionizing photons per second in the IGM can be as bright as or fainter than . Fig. 17 shows this selection effect more quantitatively using a representation similar to Fig. 9 : it shows the fraction of the time that a galaxy contributing spends shining brighter than a magnitude of -12 (in blue), -14 (in green) or -16 (in red). Assuming that the full set of snapshots from our simulations represents correctly a population of faint galaxies, this can be read as some sort of completeness function. For instance, let us focus on galaxies emitting between and ionizing photons per second, which represent the bulk of the galaxies in the mass range considered in this work, according to Kimm & Cen (2014, see their fig. 5). A survey as deep as would select all of them, while a survey limited at would only see 10% of them. This is of course a very rough estimate and should not be taken as a prediction for future surveys. Nevertheless, this resonates with the recent study of Hartley & Ricotti (2016), who used semi-analytical modelling and halo matching techniques to estimate the impact of the burstiness of the sources of reionization. They found that models with bursty star formation predict brighter galaxies at fixed halo mass.

We have so far only discussed quantities integrated in all directions, while we have shown in Fig. 8 that the escape of ionizing radiation can be very anisotropic. This will naturally lead to a larger scatter at fixed UV magnitude. The existence of preferred directions of escape will also result in anisotropic H ii regions, and this could affect the visibility of Lyman- photons at high redshift. However, since the Lyman- line is resonant, a proper treatment with dedicated radiative transfer is necessary. We plan to address this in a further study.

7 Discussion

7.1 On the importance of resolving the ISM

We have shown in this study that a large fraction of the photons emitted by young stars is absorbed in their immediate neighbourhood, corroborating the finding of Kimm & Cen (2014) and Paardekooper et al. (2015), who argued that the escape of radiation is mostly governed at the cloud scale. Similarly, Ma et al. (2015) found that the time-dependent, multiphase structure of the ISM in the vicinity of stars is an important element that regulates the escape of ionizing radiation. Earlier, Ciardi et al. (2002), Clarke & Oey (2002), and then Fernandez & Shull (2011) discussed how the porosity or clumpiness of the ISM affects . We wish now to discuss several aspects of the ISM scale physics, and how they would impact our results.

Inner structure of molecular clouds

Before going any further, we would like to discuss the question of the resolution of our series of simulation. Star formation being a process that occurs on scales well below the parsec, no cosmological simulation will be able to resolve the formation of individual stars in the near future, and even though the available computing power continues to increase, current simulations are limited. But what resolution is high enough? Wise & Cen (2009) suggest that a resolution of 0.1 pc is needed to resolve the early phases of the expansion of the ionization front, but this is unattainable for the current generation of cosmological simulations. Rosdahl et al. (2015) reached the rather depressing conclusion that, with a standard star formation model, H ii regions can never be properly resolved, just due to the fact that there is a limited amount of gas available in a cell for star formation. In light of this, we limit ourselves to a maximum resolution of pc. This allows us to start to resolve the internal structure of our galaxies, which is mandatory to describe the sites where radiation gets absorbed.

Recently, Kimm et al. (2017) presented a study of the escape of radiation from mini-haloes, reaching the extremely high resolution of 0.7 pc. While their setup is similar to ours, they find that on some occasions, photons manage to escape even before the star particles reach their supernova phase. They interpret this as a sign that radiative feedback from the young stars can disrupt the star-forming clouds very early. We do not see this happening in our simulation, even though we include most of the relevant processes. Indeed, Kimm et al. (2017) find that the main channel of radiative feedback is through photoionization, which we also include in this work. This is in part due to our resolution, as we do not resolve the inner structure of the clouds. Most importantly, they focus on much lower mass haloes (reaching a maximum mass of ) for which they include the feedback from massive pop III stars. In this study, we only considered more massive haloes above the atomic cooling limit.

Using high resolution simulations of molecular clouds with radiative feedback, Dale et al. (2012) found that only a fraction of the photons emitted during the lifetime of a cloud will escape to the more diffuse ISM, and that fraction varies from to depending on the properties of the cloud. The work of Geen et al. (2015b) extends this analysis, focusing on the expansion of an ionization front in a turbulent, magnetized, self-gravitating molecular cloud. They focus on a series of very high resolution ( pc) simulations of a cloud embedded in a 27 pc box. They find that depending on the structure of the cloud and on the gas inflow rate, radiation can be completely trapped inside the cloud. Obviously, our cosmological scale simulation does not capture the inner structure of the cloud: the whole simulation volume of Geen et al. (2015b) is described in our simulation by a few star-forming cells. Nevertheless, they provide an analytical framework to estimate the maximum expected size of an H ii region in a molecular cloud (see their eq. 15). This formalism relies on only a few parameters describing the density profile of molecular clouds, and once these parameters have been calibrated against high resolution, galactic scale simulations of molecular clouds, we could in principle apply this to our star-forming cells. This would yield a more realistic description of the escape of ionizing radiation from the stellar birth clouds, which are currently assumed to be mostly homogeneous. For instance, the simulation of a massive cloud presented in Howard et al. (2017) showed that environment of stellar clusters even inside the cloud plays an important role in the regulation of . Because in the simulations of Geen et al. (2015b), radiation leaks out of the clouds before the end of the main sequence evolution of the most massive stars, we would expect that the instantaneous increases, and that the time delay between the peaks of SFR and of would be reduced. However, developing this kind of framework would require a significant amount of work, since the study of Geen et al. (2015b) only focuses on low-mass clouds.

Stellar modelling

While other sources of pre-SN feedback such as the inclusion of stellar winds or radiation pressure from infrared radiation have been suggested in the past, it is unclear that they might strongly affect the escape of ionizing photons. For instance, Geen et al. (2015a) suggested that the size of the gaseous shell formed by the main sequence stellar feedback is largely controlled by the radiative photoheating from the star, dominating over stellar winds. We note however that the study of Paardekooper et al. (2011) indicates that the pre-processing of the ISM in the vicinity of the ionizing sources by stellar feedback before SNe increases the ability of photons to escape the star forming cloud.

A more crucial assumption we made is the fact that all supernovae explode after 10 Myr. While this is the lifetime of a star, which would be a typical SNii progenitor, we should in principle sample the stellar lifetime properly, or at least use a variable time delay between the formation of a star particle and the release of energy due to the SN. Indeed, depending on the mass of the progenitor, SNe can occur as soon as 3 Myr, or as late as 40 Myr. However, we checked that changing the time delay in the simulation to 3 Myr has only a small impact on the escape of radiation: we compare in Fig. 11 the evolution of for the run with a delay of 10 Myr (in blue) and 3 Myr (in red), and the behaviour is qualitatively the same.

Finally, a note of caution regarding the stellar model we use is in order. We would like to draw the attention of the reader to the fact that we use a very small mass for our star particles, of only . At these scales, the IMF cannot be fully sampled: Crowther et al. (2010) found, for instance, several stars with masses above 150 in the R136 star cluster in the Milky Way. We can obviously not capture the effect of these extremely massive stars with our setup. It would be necessary to stochastically sample the IMF, and use dedicated stellar evolution and SED models. However, as a result of our star formation and feedback schemes, the stars are very clustered spatially in our simulation, and there are groups of almost coeval star particles, somehow mitigating this issue.

On top of the various physical ingredients resulting in stellar feedback, the stellar population model itself might be of importance for a proper estimation of the ionizing budget. For instance, Topping & Shull (2015) showed that stellar rotation increases the production of ionizing photons, over a longer period. Stanway et al. (2016) proposed that binary interactions might result in a similar enhancement of the ionizing luminosity. Recently, Ma et al. (2016) implemented this in their simulations, and find significantly higher escape fractions at all times. Runaway stars have also been advanced as a way to increase the photon production rate in high- galaxies Conroy & Kratter (2012), but the work of Kimm & Cen (2014) and Ma et al. (2015) ruled out any significant impact of such stars.

7.2 Other limitations

Apart from the small scale processes discussed in the previous section, our simulations make two supplementary simplifying assumptions: we include no external UV background, and we assume that in our galaxies, there is no additional feedback from a central, massive black hole.

Because we are mostly interested in the escape of radiation from the complex gas structure around low mass galaxies, we do not include any contribution from an external UV background. Recently, Simpson et al. (2013) studied the impact of various physical processes on the formation of a dwarf galaxies (), and found that the inclusion of a meta-galactic UV background has mainly the effect of removing diffuse gas around the halo, and that it is SN feedback that is mostly responsible of removing the dense gas that lies in the way of ionizing radiation. Paardekooper et al. (2015) showed that an external ionizing background results in the photoevaporation of gas in low mass haloes, helping to increase . At the same time, the star formation is strongly suppressed, therefore keeping the total contribution of very small haloes quite low. Overall, the main effect of the the external background is to remove the supply of fresh neutral gas, which would delay or block later star formation. Even so, the simulations of Rosdahl & Blaizot (2012) constrain the global effect of the UV background to the gas at densities below , as denser gas is self-shielded. As a safety check, we reran the intermediate mass halo with an extra ionizing background, and found very little difference. This is because the haloes under study are too massive for star formation to be completely shut down.

Our simulations do not model the formation of massive black holes in the centre of our dwarf galaxies. This is aligned with most of the numerical studies focusing on the process of reionization by stellar sources. Indeed, while ubiquitously found in massive galaxies, it is only recently that Reines et al. (2013) discovered that massive black holes are not infrequent in low mass galaxies. Because high redshift galaxies are much more gas rich than their local counterparts, these black holes are likely to be more active. While we do not expect a strong direct contribution of these faint active galactic nuclei to the ionizing budget, the additional feedback due to the presence of a central massive black hole can in principle strongly affect the gas morphology, therefore affecting the escape of ionizing radiation.

8 Conclusions

Current constraints on the observed population of high redshift galaxies seem to indicate that low-mass galaxies play a significant role in the supply of ionizing photons needed to reionize the Universe by . Yet, their radiative efficiency, and in particular their escape fraction , is poorly constrained. To address this, we have performed a set of cosmological zoom simulations of high redshift dwarf galaxies, with very high spatial resolution ( pc), and including a physically motivated subgrid description of star formation and supernova feedback. We have used the multifrequency, hydro-coupled radiative transfer module of Ramses to follow on the fly the emission of ionizing photons and their journey through the galactic halo to reach the intergalactic medium. While we do not have the statistical significance to assess directly the importance of the contribution of dwarf galaxies to the ionizing budget of reionization, we studied in great detail the processes affecting the evolution of the escape fraction, and our results are the following:

  • The combination of spatially clustered star formation and efficient supernova feedback regulates the stellar mass of our simulated galaxies.

  • For galaxies with stellar masses of , the value of the instantaneous can vary by more than 6 orders of magnitude in Myr at a fixed (Fig. 4). This variability can explain the large scatter in vs. relation observed in simulations of cosmological volumes.

  • There is a delay between the peak of star formation and the peak of corresponding to the lifetime of massive stars (Fig. 5). This illustrate that the escape of ionizing radiation is mostly controlled by supernova feedback.

  • The escape of radiation from the halo is enabled by supernova feedback, resulting in a very anisotropic distribution of (Fig. 8).

  • Without supernova feedback, all the radiation is absorbed locally within the stellar birth cloud (Fig. 10 and 14). When a massive star ends as a supernova, the mechanical energy removes gas from the neighbouring star forming clouds, thus clearing sight lines for radiation to escape.

  • While the possibility for a galaxy to leak ionizing radiation is mainly determined by cloud-scale properties, the gas distribution in the circumgalactic medium plays an important role in the determination of the exact value of at the virial radius, modulating the amount of ionizing radiation escaping the ISM.

  • The UV luminosity of a galaxy is only a very poor proxy for its ionizing emission: a UV-bright galaxy can be in a non-leaking phase (Fig. 16). Other diagnostics are needed to select promising Lyman-continuum emitters candidates.

  • Reciprocally, the mismatch between ionizing and non-ionizing luminosities means that even relatively deep UV limited surveys will systematically miss an important fraction of the sources of reionization.

This work depicts a scenario in which the escape fraction is a galaxy-scale quantity strongly tied to processes happening at the molecular cloud level. Our study shows how stellar feedback is necessary to rearrange the gas in the halo to allow the release of radiation in the IGM. While photoionization on its own is not the dominant process in creating channels through which ionizing radiation can escape, it is still important to model it correctly. Indeed, previous work (e.g. Rosdahl et al., 2015, especially their Fig. 3) have shown that photoionization can contribute to smooth the ISM, thus altering the possible escape routes for ionizing photons.

The burstiness of the star formation history of such low mass galaxies might strongly impact their observability and inferred properties (see e.g. Domínguez et al., 2015), the study of which is critical for the design of large surveys with the upcoming James Webb Space Telescope. We plan to study this in further details in a future work.


The authors wish to thank the anonymous referee for comments that helped improve this manuscript. We would like to thank Thibault Garel, Sam Geen, Taysun Kimm, Anne Verhamme and all the participants of the ORAGE project for useful discussion. We are very grateful to Leindert Boogaard for his careful bug-hunting. MT and JB acknowledge support from the Lyon Institute of Origins under grant ANR-10-LABX-66, and from the ORAGE project from the Agence Nationale de la Recherche under grand ANR-14-CE33-0016-03. MT acknowledges funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 614199, project ‘BLACK’), and the hospitality of the University of Oxford and New College through the award of a Balzan Fellowship. JR was funded by the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 278594, project ‘GasAroundGalaxies’) and the Marie Curie Training Network CosmoComp (PITN-GA-2009-238356). This work has made extensive use of the PyMSES6 analysis package (Guillet et al., 2013).


  1. pubyear: 2017
  2. pagerange: Fluctuating feedback-regulated escape fraction of ionizing radiation in low-mass, high-redshift galaxiesFluctuating feedback-regulated escape fraction of ionizing radiation in low-mass, high-redshift galaxies
  3. https://bitbucket.org/rteyssie/ramses/
  4. https://bitbucket.org/ohahn/music/
  5. While the Strömgren analysis only holds for an homogeneous medium with a static source, and assumes no backreaction of the radiation on the gas, we take it as a lower limit for the radius of the ionization front. Indeed, if the front expands further, it will do so quicker than the lifetime of the massive stars. Since we do not resolve properly the internal structure of star forming clouds, this is a first order approximation that allows analysis without introducing more free parameters.
  6. https://bitbucket.org/dchapon/pymses/


  1. Alavi A., et al., 2014, \apj, 780, 143
  2. Atek H., et al., 2015, \apj, 814, 69
  3. Barkana R., Loeb A., 2001, \physrep, 349, 125
  4. Becker R. H., et al., 2001, \aj, 122, 2850
  5. Benson A. J., Sugiyama N., Nusser A., Lacey C. G., 2006, \mnras, 369, 1055
  6. Bergvall N., Zackrisson E., Andersson B.-G., Arnberg D., Masegosa J., Östlin G., 2006, \aap, 448, 513
  7. Bergvall N., Leitet E., Zackrisson E., Marquart T., 2013, \aap, 554, A38
  8. Borthakur S., Heckman T. M., Leitherer C., Overzier R. A., 2014, Science, 346, 216
  9. Boutsia K., et al., 2011, \apj, 736, 41
  10. Bouwens R. J., et al., 2009, \apj, 705, 936
  11. Bouwens R. J., et al., 2012, \apjl, 752, L5
  12. Bouwens R. J., et al., 2014, \apj, 793, 115
  13. Bouwens R. J., Smit R., Labbé I., Franx M., Caruana J., Oesch P., Stefanon M., Rasappu N., 2016, \apj, 831, 176
  14. Bruzual G., Charlot S., 2003, \mnras, 344, 1000
  15. Chabrier G., 2003, \pasp, 115, 763
  16. Ciardi B., Bianchi S., Ferrara A., 2002, \mnras, 331, 463
  17. Clarke C., Oey M. S., 2002, \mnras, 337, 1299
  18. Conroy C., Kratter K. M., 2012, \apj, 755, 123
  19. Cowie L. L., Barger A. J., Trouille L., 2009, \apj, 692, 1476
  20. Crowther P. A., Schnurr O., Hirschi R., Yusof N., Parker R. J., Goodwin S. P., Kassim H. A., 2010, \mnras, 408, 731
  21. Dale J. E., Ercolano B., Bonnell I. A., 2012, \mnras, 424, 377
  22. Dayal P., Ferrara A., Gallerani S., 2008, \mnras, 389, 1683
  23. Domínguez A., Siana B., Brooks A. M., Christensen C. R., Bruzual G., Stark D. P., Alavi A., 2015, \mnras, 451, 839
  24. Duncan K., Conselice C. J., 2015, \mnras, 451, 2030
  25. Fan X., et al., 2001, \aj, 122, 2833
  26. Fan X., et al., 2006, \aj, 132, 117
  27. Federrath C., Klessen R. S., 2012, \apj, 761, 156
  28. Ferland G. J., et al., 2013, \rmxaa, 49, 137
  29. Fernandez E. R., Shull J. M., 2011, \apj, 731, 20
  30. Finkelstein S. L., et al., 2012, \apj, 758, 93
  31. Finkelstein S. L., Dunlop J., Le Fevre O., Wilkins S., 2015, preprint, (arXiv:1512.04530)
  32. Fontanot F., Cristiani S., Vanzella E., 2012, \mnras, 425, 1413
  33. Fujita A., Martin C. L., Mac Low M.-M., Abel T., 2003, \apj, 599, 50
  34. Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015a, \mnras, 448, 3248
  35. Geen S., Hennebelle P., Tremblin P., Rosdahl J., 2015b, \mnras, 454, 4484
  36. Gnedin N. Y., 2014, \apj, 793, 29
  37. Gnedin N. Y., Kravtsov A. V., Chen H.-W., 2008, \apj, 672, 765
  38. Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, \apj, 622, 759
  39. Grazian A., et al., 2016, \aap, 585, A48
  40. Grissom R. L., Ballantyne D. R., Wise J. H., 2014, \aap, 561, A90
  41. Guillet T., Chapon D., Labadens M., 2013, PyMSES: Python modules for RAMSES, Astrophysics Source Code Library (ascl:1310.002)
  42. Gunn J. E., Peterson B. A., 1965, \apj, 142, 1633
  43. Haardt F., Madau P., 2012, \apj, 746, 125
  44. Haardt F., Salvaterra R., 2015, \aap, 575, L16
  45. Hahn O., Abel T., 2011, \mnras, 415, 2101
  46. Hartley B., Ricotti M., 2016, \mnras, 462, 1164
  47. Heckman T. M., Sembach K. R., Meurer G. R., Leitherer C., Calzetti D., Martin C. L., 2001, \apj, 558, 56
  48. Heckman T. M., et al., 2011, \apj, 730, 5
  49. Hennebelle P., Chabrier G., 2011, \apjl, 743, L29
  50. Hinshaw G., et al., 2013, \apjs, 208, 19
  51. Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, \mnras, 445, 581
  52. Howard C., Pudritz R., Klessen R., 2017, \apj, 834, 40
  53. Huang K.-H., et al., 2016, \apjl, 823, L14
  54. Hui L., Gnedin N. Y., 1997, \mnras, 292, 27
  55. Iliev I. T., et al., 2006, \mnras, 371, 1057
  56. Iliev I. T., et al., 2009, \mnras, 400, 1283
  57. Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, \mnras, 439, 725
  58. Iwata I., et al., 2009, \apj, 692, 1287
  59. Izotov Y. I., Orlitova I., Schaerer D., Thuan T. X., Verhamme A., Guseva N. G., Worseck G., 2016, \nat, 529, 178
  60. Kim J.-h., Krumholz M. R., Wise J. H., Turk M. J., Goldbaum N. J., Abel T., 2013, \apj, 775, 109
  61. Kimm T., Cen R., 2014, \apj, 788, 121
  62. Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, \mnras, 451, 2900
  63. Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017, \mnras,
  64. Koekemoer A. M., et al., 2013, \apjs, 209, 3
  65. Kuhlen M., Faucher-Giguère C.-A., 2012, \mnras, 423, 862
  66. Leitet E., Bergvall N., Hayes M., Linné S., Zackrisson E., 2013, \aap, 553, A106
  67. Ma X., Kasen D., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Kereš D., Murray N., 2015, \mnras, 453, 960
  68. Ma X., Hopkins P. F., Kasen D., Quataert E., Faucher-Giguère C.-A., Kereš D., Murray N., Strom A., 2016, \mnras, 459, 3614
  69. Madau P., Haardt F., Rees M. J., 1999, \apj, 514, 648
  70. Malkan M., Webb W., Konopacky Q., 2003, \apj, 598, 878
  71. Miller S. H., Ellis R. S., Newman A. B., Benson A., 2014, \apj, 782, 115
  72. Mitra S., Choudhury T. R., Ferrara A., 2015, \mnras, 454, L76
  73. Mostardi R. E., Shapley A. E., Nestor D. B., Steidel C. C., Reddy N. A., Trainor R. F., 2013, \apj, 779, 65
  74. Mostardi R. E., Shapley A. E., Steidel C. C., Trainor R. F., Reddy N. A., Siana B., 2015, \apj, 810, 107
  75. Mutch S. J., Geil P. M., Poole G. B., Angel P. W., Duffy A. R., Mesinger A., Wyithe J. S. B., 2016, \mnras, 462, 250
  76. Nestor D. B., Shapley A. E., Kornei K. A., Steidel C. C., Siana B., 2013, \apj, 765, 47
  77. Ocvirk P., et al., 2016, \mnras,
  78. Oke J. B., Gunn J. E., 1983, \apj, 266, 713
  79. Paardekooper J.-P., Pelupessy F. I., Altay G., Kruip C. J. H., 2011, \aap, 530, A87
  80. Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2013, \mnras, 429, L94
  81. Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2015, \mnras, 451, 2544
  82. Padoan P., Nordlund Å., 2011, \apj, 730, 40
  83. Planck Collaboration et al., 2014, \aap, 571, A16
  84. Planck Collaboration et al., 2016, \aap, 594, A13
  85. Rasera Y., Teyssier R., 2006, \aap, 445, 1
  86. Razoumov A. O., Sommer-Larsen J., 2010, \apj, 710, 1239
  87. Reines A. E., Greene J. E., Geha M., 2013, \apj, 775, 116
  88. Ricotti M., Gnedin N. Y., Shull J. M., 2002, \apj, 575, 33
  89. Robertson B. E., et al., 2013, \apj, 768, 71
  90. Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, \apjl, 802, L19
  91. Rosdahl J., Blaizot J., 2012, \mnras, 423, 344
  92. Rosdahl J., Teyssier R., 2015, \mnras, 449, 4380
  93. Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, \mnras, 436, 2188
  94. Rosdahl J., Schaye J., Teyssier R., Agertz O., 2015, \mnras, 451, 34
  95. Rosdahl J., Schaye J., Dubois Y., Kimm T., Teyssier R., 2017, \mnras, 466, 11
  96. Rosen A., Bregman J. N., 1995, \apj, 440, 634
  97. Shapley A. E., Steidel C. C., Pettini M., Adelberger K. L., Erb D. K., 2006, \apj, 651, 688
  98. Shen S., Madau P., Conroy C., Governato F., Mayer L., 2014, \apj, 792, 99
  99. Siana B., et al., 2007, \apj, 668, 62
  100. Siana B., et al., 2010, \apj, 723, 241
  101. Siana B., et al., 2015, \apj, 804, 17
  102. Simpson C. M., Bryan G. L., Johnston K. V., Smith B. D., Mac Low M.-M., Sharma S., Tumlinson J., 2013, \mnras, 432, 1989
  103. Stanway E. R., Eldridge J. J., Becker G. D., 2016, \mnras, 456, 485
  104. Teyssier R., 2002, \aap, 385, 337
  105. Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, \mnras, 429, 3068
  106. Topping M. W., Shull J. M., 2015, \apj, 800, 97
  107. Toro E. F., Spruce M., Speares W., 1994, Shock Waves, 4, 25
  108. Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, \aap, 506, 647
  109. Vanzella E., et al., 2010, \apj, 725, 1011
  110. Vanzella E., et al., 2015, \aap, 576, A116
  111. Vanzella E., et al., 2016, \apj, 825, 41
  112. Whalen D., O’Shea B. W., Smidt J., Norman M. L., 2008, \apj, 679, 925
  113. Willott C. J., et al., 2010, \aj, 139, 906
  114. Wise J. H., Cen R., 2009, \apj, 693, 984
  115. Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, \mnras, 442, 2560
  116. Wood K., Loeb A., 2000, \apj, 545, 86
  117. Wyithe J. S. B., Cen R., 2007, \apj, 659, 890
  118. Yajima H., Choi J.-H., Nagamine K., 2011, \mnras, 412, 411
  119. Yajima H., Li Y., Zhu Q., Abel T., Gronwall C., Ciardullo R., 2014, \mnras, 440, 776
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description