The visibility of Lyman α emitters: constraining reionization, ionizing photon escape fractions and dust

Abstract

1 INTRODUCTION

The epoch of reionization (EoR) begins when the first stars start producing neutral hydrogen (H i) ionizing photons, carving out an ionized Strömgren region in the neutral intergalactic medium (IGM) around themselves, and ends when all the H i in the IGM is ionized. Understanding the process of cosmic reionization is extremely important because in addition to marking the last major change in the ionization state of the Universe, it affected all subsequent structure formation through a number of radiative feedback effects (see e.g. Barkana & Loeb 2001; Ciardi & Ferrara 2005; Maio et al. 2011; Sobacchi & Mesinger 2013; Wyithe & Loeb 2013, and references therein). A broad picture that has emerged (Choudhury & Ferrara 2007) is one wherein hydrogen reionization is an extended process that starts at z ≈ 15 and is about 90 per cent complete by z ≃ 8; while it is initially driven by metal-free Population III stars in low-mass haloes ( ≤ 10⁸M_⊙), the conditions for star formation (SF) in these haloes are soon erased by a combination of chemical and radiative feedback by z ≃ 10. However, any attempt at modelling reionization necessarily requires a number of important assumptions regarding the number of H i ionizing photons produced by the source galaxy population depending on their physical properties, the fraction of these photons (f_esc) that can escape out of the galactic environment and contribute to reionization, and the clumping factor of the IGM at high redshifts, to name a few (Salvaterra, Ferrara & Dayal 2011). Given these assumptions, reionization scenarios need to be constantly updated as new data sets are acquired.

By virtue of their continually growing numbers, a valuable new data set is provided by high-redshift Lyman α emitters (LAEs); this is a class of galaxies identified by means of their Lyman α (Lyα) emission line at 1216 Å in the galaxy rest frame. Indeed, hundreds of LAEs have now been confirmed in the EoR: z ≃ 5.7 (Malhotra et al. 2005; Shimasaku et al. 2006; Hu et al. 2010; Curtis-Lake et al. 2012), z ≃ 6.6 (Taniguchi et al. 2005; Kashikawa et al. 2006, 2011; Hu et al. 2010; Ouchi et al. 2010) and z ≃ 7 (Iye et al. 2006; Ouchi et al. 2009; Stark et al. 2010; Pentericci et al. 2011). In addition to their number statistics, LAEs have been rapidly gaining popularity as probes of reionization and high-redshift galaxy evolution for two reasons: (a) the strength, width and continuum break bluewards of the Lyα line make their detection unambiguous, and (b) Lyα photons are extremely sensitive to attenuation by H i; the observed Lyα luminosity can then be used to infer the ionization state of the IGM at redshifts close to those of the emitter and hence to reconstruct, at least piecewise, the cosmic reionization history.

However, interpreting LAE data is complicated by a number of physical effects. First, the intrinsic Lyα luminosity depends on the total number of H i ionizing photons that are produced by a galaxy, depending on the star formation rate (SFR), age and metallicity of its stellar population (e.g. Santos 2004). Secondly, depending on the H i and dust contents in the interstellar medium (ISM), only a fraction (1 − f_esc) of these H i ionizing photons are able to contribute to the intrinsic Lyα luminosity by ionizing the ISM H i; the rest contribute to building the ionized H ii region around the galaxy. We briefly digress to note that the dependence of f_esc on galaxy properties has been the subject of much recent debate: while some authors find f_esc to decrease with an increase in the halo mass (Razoumov & Sommer-Larsen 2010; Yajima, Choi & Nagamine 2011; Ferrara & Loeb 2013), other works find the opposite trend (Gnedin, Kravtsov & Chen 2008; Wise & Cen 2009). The value of f_esc also remains poorly constrained with findings ranging from a few per cent (e.g. Gnedin et al. 2008) up to 20–30 per cent (e.g. Mitra, Ferrara & Choudhury 2013) or even higher (e.g. Wise & Cen 2009). Thirdly, only a fraction (f_α) of the Lyα photons produced inside a galaxy can escape out of it unattenuated by dust (Dayal, Ferrara & Gallerani 2008; Finkelstein et al. 2009; Nagamine et al. 2010; Forero-Romero et al. 2011; Dayal & Ferrara 2012). Fourthly, only a fraction (T_α) of the Lyα photons that emerge out of a galaxy are transmitted through the IGM and reach the observer. As expected, this transmission sensitively depends on the IGM H i ionization state and it has been shown that only galaxies residing in overlapping H ii regions would be observed as LAEs in the initial stages of reionization, i.e. reionization increases the observed clustering of LAEs (McQuinn et al. 2007a; Dayal et al. 2009). Fifthly, the Lyα IGM transmission calculation is complicated by the presence of peculiar velocities. Inflows (outflows) of gas into (from) the emitter blueshift (redshift) the Lyα line, leading to a decrease (increase) in the value of T_α along different lines of sight (LOS), strongly affecting the visibility of galaxies as LAEs (Verhamme, Schaerer & Maselli 2006; Iliev et al. 2008; Zheng et al. 2010; Dijkstra, Mesinger & Wyithe 2011). Indeed, Santos (2004) has shown that the IGM ionization state cannot be constrained in the presence of peculiar velocities. However, Dayal et al. (2011) have shown that the interpretation of LAE data is much more involved; a decrease in the Lyα transmission due to peculiar velocities and/or a highly neutral IGM can be compensated by an increase in f_α due to dust being clumped in the ISM of LAEs, as per the ‘Neufeld model’ (Neufeld 1991).

A number of past theoretical works have used one or more of the above ingredients to use LAEs as tracers of reionization. With semi-analytic modelling, a number of authors (e.g. Dijkstra, Lidz & Wyithe 2007; Dayal et al. 2008; Samui, Srianand & Subramanian 2009) have shown that the LAE Lyα luminosity functions (LFs) are consistent with a fully ionized IGM and can be explained solely by an evolution of the underlying dark matter (DM) halo mass functions. However, when considered in combination with the non-evolving observed ultraviolet (UV) LFs of LAEs between z ≃ 5.7 and 6.6 (Kashikawa et al. 2006), some of these works (Dijkstra et al. 2007; Ouchi et al. 2010) argue for an additional dimming of the Lyα line by about 30 per cent. A number of studies have used large pure DM simulations to study the clustering of LAEs (McQuinn et al. 2007b; Iliev et al. 2008; Orsi et al. 2008). Finally, a number of studies have been undertaken using cosmological hydrodynamic simulations with dust (Dayal et al. 2008, 2009; Nagamine et al. 2010; Forero-Romero et al. 2011), with radiative transfer (RT) without dust (Zheng et al. 2010) and RT with dust (Dayal et al. 2011; Forero-Romero et al. 2011; Duval et al. 2014).

This work is quite close in spirit to the calculations presented in Dayal et al. (2011) wherein the authors used (a) cosmological hydrodynamic simulations run with gadget-2 to obtain the physical properties of z ≃ 5.7 galaxies, (b) a dust model that took into account the entire SF history of each galaxy to calculate its dust enrichment and (c) an RT code (crash) to obtain the ionization fields to calculate the IGM Lyα transmission for each galaxy, in order to identify the simulated galaxies that would be observationally classified as LAEs. Using this model, the authors showed that the effects of dust and IGM are degenerate in affecting the visibility of LAEs such that a large f_α can be compensated by a small T_α and vice versa to yield a given value of the observed Lyα luminosity. However, the main caveat in that work was that the authors used a constant value of f_esc = 0.02 (Gnedin et al. 2008) for all galaxies in their calculations, and started their RT runs assuming the IGM gas to be in photoionization equilibrium with a uniform ultraviolet background (UVB) produced by unresolved sources corresponding to an average neutral hydrogen fraction, |$\langle \chi _{\rm H\,\small{I}} \rangle = 0.3$|⁠.

In this work, we substantially enhance the model presented in Dayal et al. (2011) to build a self-consistent model that couples cosmological smoothed particle hydrodynamic (SPH) simulations, dust modelling and a fast RT code (pcrash) to identify LAEs without making any prior assumptions on the IGM ionization state, f_esc and the ISM dust distribution. Starting from a uniform neutral IGM, consistent with the Haardt & Madau (1996) background imposed in simulations, we use five different values of f_esc = 0.05, 0.25, 0.5, 0.75, 0.95 in order to study how varying f_esc affects the visibility of LAEs through (a) the direct impact on the intrinsic Lyα luminosity, and (b) the IGM Lyα transmission that depends on the topology and extent of H ii regions as determined by f_esc. Comparing the statistics of the simulated LAEs to observations, our aim is to jointly constrain three fundamental parameters that determine the visibility of LAEs – f_esc, |$\chi _{\rm H\,\small{I}}$| (or T_α) and the relative escape of Lyα photons with respect to continuum photons from the dusty ISM of galaxies (f_α/f_c).

We begin by describing the hydrodynamical simulation, the dust model and the identification of galaxies as Lyman break galaxies (LBGs) through the UV continuum in Section 2. We follow this approach since UV photons are only attenuated by dust but unaffected by the IGM ionization state, simplifying their identification. In order to validate the simulations, we compare the simulated LBG UV LFs to observations in Section 3; a comparison of the simulated LBG stellar mass functions, stellar mass densities (SMDs) and specific star formation rates (sSFRs) with the observations is shown in Appendix A for completeness. Once the physical properties of the simulated galaxies have been validated, we choose a simulation snapshot at the edge of the EoR (z ≃ 6.7) and post-process it with a fast RT code (pcrash) with the five values of f_esc mentioned above, as explained in Section 4. We then compare the simulated LAE UV and Lyα LFs with observations (Kashikawa et al. 2011) in Section 5 in order to jointly constrain f_esc, |$\chi _{\rm H\,\small{I}}$| and f_α/f_c. Due to the model-dependent relation between |$\chi _{\rm H\,\small{I}}$| and the Lyα IGM transmission (T_α), we show our constraints also in terms of f_esc, T_α and f_α/f_c, before concluding in Section 6.

2 COSMOLOGICAL HYDRODYNAMIC SIMULATIONS

In this section, we describe the cosmological simulation used to obtain the physical properties of high-redshift galaxies and the semi-analytic model used to obtain their dust enrichment, in order to calculate their observed visibility in the UV.

2.1 The simulation

The hydrodynamical simulation analysed in this work has been carried out using the TreePM-SPH code gadget-2 (Springel 2005). The adopted cosmological model corresponds to the Λ cold dark matter universe with dark matter (DM), dark energy and baryonic density parameter values of (Ω_Λ, Ω_m, Ω_b) = (0.73, 0.27, 0.047), a Hubble constant H₀ = 100 h = 70 km s⁻¹ Mpc⁻¹ and a normalization σ₈ = 0.82, consistent with the results from Wilkinson Microwave Anisotropy Probe 5 (Komatsu et al. 2009). The simulation box has a size of 80 h⁻¹ comoving Mpc (cMpc), and contains 1024³ DM particles, and initially the same number of gas particles; the mass of a DM and gas particle is 3.6 × 10⁷ h^{− 1}M_⊙ and 6.3 × 10⁶ h^{− 1}M_⊙, respectively. The softening length for the gravitational force is taken to be 3 h⁻¹ comoving kpc and the value of the smoothing length for the SPH kernel for the computation of hydrodynamic forces is allowed to drop at most to the gravitational softening.

The runs include the SF prescriptions of Springel & Hernquist (2003) such that the ISM is described as an ambient hot gas containing cold clouds, which provide the reservoir for SF, with the two phases being in pressure equilibrium; the relative number of stars of different masses is computed using the Salpeter initial mass function (IMF; Salpeter 1955) between 0.1 and 100 M_⊙. The density of the cold and of the hot phase represents an average over small regions of the ISM, within which individual molecular clouds cannot be resolved by simulations sampling cosmological volumes. The runs also include the feedback model described in Springel & Hernquist (2003) which includes (a) thermal feedback: supernovae (SN) inject entropy into the ISM, heat up nearby particles and destroy molecules, (b) chemical feedback: metals produced by SF and SN are carried by winds and pollute the ISM, and (c) mechanical feedback: galactic winds powered by SN. In the case of mechanical feedback, the momentum and energy carried away by the winds are calculated assuming that the galactic winds have a fixed velocity of 500 km s⁻¹ with a mass upload rate equal to twice the local SFR, and carry away a fixed fraction (50 per cent) of the SN energy (for which the canonical value of 10⁵¹ erg is adopted). Finally, the run assumes a metallicity-dependent radiative cooling (Sutherland & Dopita 1993) and a uniform redshift-dependent UVB produced by quasars and galaxies as given by Haardt & Madau (1996).

Galaxies are recognized as gravitationally bound groups of at least 20 total (DM+gas+star) particles using the Amiga Halo Finder (ahf; Knollmann & Knebe 2009). When compared to the standard Sheth–Tormen mass function (Sheth & Tormen 1999), the simulated mass function is complete for halo masses M_h ≥ 10^9.2 M_⊙ for z ≃ 6–8; galaxies above this mass cut-off are referred to as the ‘complete sample’. Of this complete sample, we identify all the ‘resolved’ galaxies that contain a minimum of 4N gas particles, where N = 40 is the number of nearest neighbours used in the SPH computations; this is twice the standard value of 2N gas particles needed to obtain reasonable and converging results (e.g. Bate & Burkert 1997). We impose an additional constraint and only use those resolved galaxies that contain at least 10 star particles so as to get a statistical estimate of the composite spectral energy distribution (SED). For each resolved galaxy used in our calculations (with M_h ≥ 10^9.2 M_⊙, more than 4N gas particles and a minimum of 10 star particles), we obtain the properties of all its star particles, including the redshift and mass/metallicity at formation; we also compute global properties including the total stellar mass (M_*), gas mass (M_g), DM mass (M_h), mass-weighted stellar metallicity (Z_*) and the mass weighted stellar age (t_*).

2.2 Dust model

The evidence for dust at high redshifts comes from observations of damped Lyα systems (Pettini et al. 1994; Ledoux, Bergeron & Petitjean 2002) and from the thermal dust emission from SDSS QSOs (Omont et al. 2001; Bertoldi & Cox 2002). Although dust is produced both by SN and evolved stars, several works (Todini & Ferrara 2001; Dwek, Galliano & Jones 2007) have shown that the contribution of asymptotic giant branch (AGB) stars to the total dust content is negligible at z ≳ 6, since the age of the Universe is shorter than the typical evolutionary time-scales of AGB stars above this redshift. We therefore make the hypothesis that Type II SN (SNII) are the primary dust factories and compute the total dust mass, M_d, in each of our simulated galaxies as described in Dayal, Ferrara & Saro (2010).

2.3 Identifying LBGs

To identify the simulated galaxies that could be observed as LBGs at z ≃ 6–8, we start by computing their UV luminosities. We consider each star particle to form in a burst, after which it evolves passively. The total SED, including both the stellar and nebular continuum, is computed for each star particle via the population synthesis code STARBURST99 (Leitherer et al. 1999), using its mass, stellar metallicity and age. The total intrinsic UV luminosity, |$L_{\rm c}^{\rm int}$| (at 1500 Å in the galaxy rest frame) is then calculated for each progenitor by summing the SEDs of all its star particles.

Continuum photons can be absorbed by dust within the ISM and only a fraction, f_c, escape out of any galaxy unattenuated by dust. However, these photons are unaffected by the ionization state of the IGM, and all the continuum photons that escape out of a galaxy can then reach the observer, so that the observed continuum luminosity can be expressed as |$L_{\rm c}^{\rm obs} = L_{\rm c}^{\rm int} \times f_{\rm c}$|⁠; this can be translated into an absolute magnitude, M_UV. In accordance with current observational criterion, at each redshift z ≃ 6–8, resolved simulated galaxies with M_UV ≤ −17 are identified as LBGs.

3 COMPARING THE SIMULATIONS WITH LBG OBSERVATIONS

Once we have identified the simulated galaxies that would be detected as LBGs in the snapshots at z ≃ 6–8, we compare their UV LFs, stellar mass functions, SMDs and sSFRs to the observed values in order to validate the simulations used in this work. In this section, we show a comparison between the intrinsic and dust-attenuated (observed) theoretical UV LFs and the data. The theoretical stellar mass functions, SMD and sSFR, and their comparison with observed values are shown in Appendix A.

We calculate the intrinsic UV LFs by binning simulated LBGs on the basis of their intrinsic (i.e. dust-unattenuated) UV magnitudes and dividing this by the width of the UV bin (0.5 dex) and the volume of the box. We find that the intrinsic UV LF shifts towards brighter luminosities and higher number densities with decreasing redshift from z ≃ 8 to 6, as expected from the hierarchical structure formation scenario where successively larger systems build up with time from the merger of smaller systems. As seen in Fig. 1, we find that the simulated intrinsic LBG UV LFs are overestimated with respect to the data, with the overestimation increasing with decreasing redshift, hinting at the increasing dust enrichment of these galaxies; we note that the intrinsic UV LF at z ≃ 8 is already in agreement with the observations, requiring no dust correction at this redshift.

Figure 1.

UV LFs at z ≃ 6, 7 and 8, from left to right, as marked in each panel. In all panels, the solid black (dashed red) histograms show the dust-corrected (intrinsic) simulated UV LFs with error bars showing the Poissonian error and symbols showing the observed data. The observed UV LFs have been taken from (a) z ≃ 6: Bouwens et al. (2007, empty squares), McLure et al. (2009, filled circles); (b) z ≃ 7: Bouwens et al. (2011, filled squares), McLure et al. (2010, empty circles), McLure et al. (2013, filled circles), Oesch et al. (2010, filled triangles) and (c) z ≃ 8: Bradley et al. (2012, filled triangles), Bouwens et al. (2011, filled squares), McLure et al. (2010, empty circles) and McLure et al. (2013, filled circles).

Open in new tabDownload slide

We then use the observed (i.e. dust-attenuated) UV luminosity obtained for each galaxy at z ≃ 6–8 (see Sections 2.2 and 2.3) to build the observed UV LF. It is encouraging to note that both the slope and the amplitude of the dust-attenuated simulated LFs are in agreement with the observations for z ≃ 6 and 7 as seen from Fig. 1; as mentioned before, the simulated UV LF at z ≃ 8 requires no dust to match the observations. As can be seen from the same figure, the effects of dust on continuum photons at z ≃ 6, 7 are most severe for the most massive/luminous galaxies; indeed while f_c ≃ 0.8 for galaxies in halo masses of M_h ≤ 10¹⁰ M_⊙, the value drops steadily thereafter such that f_c ≃ 0.01 for the largest galaxies at both z ≃ 6, 7. Further, as a result of galaxies typically being less massive, younger, and hence less dust enriched with increasing redshift, averaged over all LBGs f_c drops from ≃ 1 at z ≃ 8 to 0.6 at z ≃ 7 and 0.5 at z ≃ 6. We briefly digress to note that these average f_c values are about a factor of 2 higher than the values inferred in Dayal et al. (2010). This is due to the different feedback models implemented in these two simulations: while only 25 per cent of the SN energy was used to power outflows in the simulation used in Dayal et al. (2010), 50 per cent of the SN energy has been used to power outflows in the simulation used in this work. As a result of the much larger energy inputs that power SN winds in driving out gas from the galaxy, the typical stellar masses and SFRs obtained from the simulation used in this work are about a factor of 2 lower than those presented in Dayal et al. (2010).

Returning to our discussion regarding the observed LBG UV LFs, we notice that these also shift to progressively higher luminosities and/or higher number densities with decreasing redshift. Dayal et al. (2013) have shown that this evolution depends on the luminosity range probed: the steady brightening of the bright end of the LF is driven by genuine physical luminosity evolution due to a fairly steady increase in the UV luminosity (and hence SFRs) in the most massive LBGs; the evolution at the faint end arises due to a mixture of both positive and negative luminosity and density evolution as these putative systems brighten and fade, and continuously form and merge into larger systems.

Finally, we find that the best-fitting Schechter parameters to the dust-attenuated simulated UV LFs are faint-end slope α = (−1.9 ± 0.2, −2.0 ± 0.2, −1.8 ± 0.2) and the knee of the LF M_{UV, *} = ( − 19.8 ± 0.3, −19.7 ± 0.2, −19.9 ± 0.3) at z ≃ (6, 7, 8). These values are in agreement with the values α = (−1.71 ± 0.11, −1.90 ± 0.14, −2.02 ± 0.22), M_{UV, *} = ( − 20.04 ± 0.12, −19.9 ± 0.2, −20.12 ± 0.37) and α = (−1.74 ± 0.16, −2.01 ± 0.21, −1.91 ± 0.32), M_{UV, *} = ( − 20.24 ± 0.19, −20.14 ± 0.26, −20.1 ± 0.52) inferred observationally by McLure et al. (2009, 2013) and Bouwens et al. (2011) at z ≃ (6, 7, 8), respectively.

4 SIMULATING LAEs

In the previous section (and Appendix A), we have validated that the simulated galaxy population at z ≃ 6–8 is in agreement with a number of high-z LBG observations. We now use these galaxy populations to identify the simulated galaxies that could be detected as LAEs. We start by describing the RT code (pcrash) used to obtain reionization topologies for f_esc = 0.05, 0.25, 0.5, 0.75, 0.95, which are utilized to calculate the IGM Lyα transmission. We then describe how these transmission values, and the effects of dust on Lyα photons, are taken into account, so as to obtain the observed Lyα luminosity from the intrinsic value for each galaxy. We note that all the RT and Lyα calculations are carried out using a single snapshot of the hydrodynamical simulation at the edge of the reionization epoch, at z ≃ 6.7.

4.1 Reionizing the Universe with pcrash

As mentioned in Section 1, the progress of reionization critically depends on the total number of H i ionizing photons that can escape a galaxy and ionize the IGM around it. As expected, this depends both on the total number of H i ionizing photons produced by a galaxy, as well as the fraction that can escape the galactic environment (f_esc). While the simulated SFRs are in reasonable agreement with observations (see Appendix A), giving us a handle on the H i ionizing photon production rate, the value of f_esc remains only poorly understood: using a variety of theoretical models, the value of f_esc has been found to range between 0.01 and 0.8 (see e.g. Ricotti & Shull 2000; Fujita et al. 2003; Razoumov & Sommer-Larsen 2006; Wise & Cen 2009).

Due to the poor theoretical constraints available on the value of f_esc, we do not make any prior assumptions on its value. We instead post-process the hydrodynamical simulation snapshot at z ≃ 6.7 with the RT code pcrash using five different values of f_esc = 0.05, 0.25, 0.50, 0.75 and 0.95 to study how varying f_esc affects the progress of reionization.

We now briefly describe pcrash, and interested readers are referred to Ciardi et al. (2001), Maselli, Ferrara & Ciardi (2003), Maselli, Ciardi & Kanekar (2009) and Partl et al. (2011) for more details. pcrash is a 3D grid-based MPI-parallelized RT code that utilizes a combination of ray tracing and Monte Carlo schemes. It follows the time evolution of the IGM gas properties including ionization fractions (H i and He i, He ii) and the temperature. A large number of point sources with different spectra emit photon packages into the medium. The relevant matter–radiation interactions (photoionization, recombination, collisional ionization) as well as the heating/cooling of the IGM (collisional ionization cooling, recombination cooling, collisional excitation cooling, bremsstrahlung, Compton cooling/heating, adiabatic cooling, photoheating) are calculated. In addition, the code includes diffuse radiation following recombinations of ionized atoms. For the RT runs, we use the density and temperature fields obtained from the z ≃ 6.7 snapshot of the hydrodynamical simulation as the inputs for pcrash. For each resolved source in the SPH simulation, the luminosity and spectrum were calculated by summing up the spectra of all their star particles, as described in Section 2.3; we processed 31 855 resolved sources from the snapshot at z ≃ 6.7. The spectrum of each source was binned into 125 frequency intervals ranging from 3.21 × 10¹⁵ to 3.30 × 10¹⁶ Hz. Each source emitted 10⁷ rays and we used 10⁶ timesteps per 500 Myr on a 128³ grid; we used the mean local clumping factor of each cell obtained from the underlying density distribution from the SPH simulation to account for the clumpiness of the IGM (Raičević & Theuns 2011). Depending on the escape fraction, the runs were carried out on 64 or 128 processors on AIP clusters. Each simulation was run until the IGM global hydrogen fraction dropped to |$\langle \chi _{\rm H\,\small{I}} \rangle \lesssim 10^{-4}$|⁠, so that the entire progress of reionization could be mapped.

These two effects can be clearly seen in Figs 2 and 3: the former and latter show the ionization fields obtained for different f_esc by running pcrash for a given time (50 Myr) and until the IGM is ionized to a level of |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.1$|⁠, respectively. From Fig. 2, it can be clearly seen that reionization proceeds faster for increasing f_esc values: after 50 Myr of running pcrash, while the largest ionized region built has a size of about 15 cMpc for f_esc = 0.05, the entire box is almost reionized for f_esc = 0.95. Indeed, the average value of |$\chi _{\rm H\,\small{I}}$| drops steadily with increasing f_esc such that |$\langle \chi _{\rm H\,\small{I}} \rangle$| =(0.97, 0.82, 0.62, 0.37, 0.20) for f_esc = (0.05, 0.25, 0.5, 0.75, 0.95), respectively. Increasing f_esc by a factor of 19, from 0.05 to 0.95, has the effect of decreasing the neutral hydrogen fraction; while the IGM is essentially neutral for f_esc = 0.05, it is almost ionized for f_esc = 0.95.

Figure 2.

Maps (80 h⁻¹ Mpc on a side) of the spatial distribution of H i in a 2D cut through the centre of the simulated box at z ≃ 6.7 obtained by running pcrash for 50 Myr using the f_esc value marked above the panel. The colour bar shows the values (in log scale) of the H i fraction. As can be seen, reionization proceeds faster with increasing f_esc values.

Open in new tabDownload slide

Figure 3.

Maps (80 h⁻¹ Mpc on a side) of the spatial distribution of H i in a 2D cut through the simulated box at z ≃ 6.7 obtained by running pcrash using the f_esc value marked above the panel, until the average neutral hydrogen fraction drops to |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.1$|⁠. The colour bar shows the values (in log scale) of the H i fraction. As seen, while the topology of reionization looks very similar for the same average H i values, the degree of ionization close to any source increases with increasing f_esc.

Open in new tabDownload slide

From Fig. 3, we find that the topology of reionization looks similar for all the five f_esc values by the time the average H i fraction has dropped to |$\langle \chi _{\rm H\,\small{I}} \rangle =0.1$|⁠. However, the time taken to reach this average ionization state is very different for the five runs; while the run with f_esc = 0.95 took 80 Myr to reach this ionization level, the runs with f_esc = 0.05 took about 10 times longer (≃800 Myr) to build up these ionized regions. We also note that although the spatial distribution of the ionization fields looks very similar for the varying f_esc values, the level of ionization in any given cell increases with increasing f_esc as seen from the same figure (see also equation 3).

To summarize, we find that increasing f_esc both accelerates the progress of reionization and leads to higher ionization fractions in the ionized regions. A combination of both these factors leads to a larger IGM Lyα transmission T_α with increasing f_esc, as explained in the following.

4.2 Identifying LAEs

In this section, we start by explaining how we calculate the intrinsic Lyα luminosity produced by each simulated galaxy at z ≃ 6.7. We then show how we calculate the Lyα attenuation by ISM dust and IGM H i to determine the fraction of this intrinsic luminosity that is finally observed.

4.2.1 Intrinsic Lyα luminosity

4.2.2 Effects of dust

The intrinsic Lyα luminosity produced by a galaxy has to penetrate through dust in the ISM, with only a fraction f_α emerging out of the galactic environment. In Section 2.2, we have shown calculations for the total dust enrichment of each simulated galaxy and the resulting value of f_c. The relation between f_α and f_c depends on the adopted extinction curve if dust is homogeneously distributed; it also depends on the differential effects of RT on Lyα and UV photons if dust is inhomogeneously distributed/clumped (Neufeld 1991; Hansen & Oh 2006). While some pieces of evidence exist that the SN extinction curve (Bianchi & Schneider 2007) can successfully be used to interpret the observed properties of the most distant quasars (Maiolino et al. 2006) and gamma-ray bursts (Stratta et al. 2007), the effect of dust inhomogeneities in enhancing the escape fraction of Lyα photons with respect to continuum photons through the Neufeld effect (Neufeld 1991) remains controversial. On the one hand, Finkelstein et al. (2009) have found tentative observational evidence of f_α/f_c > 1 for z ≃ 4.5 LAEs, and Dayal et al. (2009, 2011) have shown that models require f_α/f_c > 1 to reproduce LAE data at z ≲ 6. On the other hand, Laursen, Duval & Östlin (2013) have shown that a value of f_α > f_c requires very special circumstances (no bulk outflows, very high metallicity, very high density of the warm neutral medium and a low-density and highly ionized medium) that are unlikely to exist in any realistic ISM; moreover, they show that a value of f_α/f_c > 1 results in very narrow Lyα lines that are sensitive to infalling gas, resulting in low T_α values. Due to its poor understanding, the relative escape fraction f_α/f_c is left as a free parameter in our model and its value is fixed by matching the theoretical LAE Lyα and UV LFs to the observations, as shown in Section 5 that follows.

4.2.3 Effects of the IGM

We also remind the reader that T_α is fixed from the ionization state of the IGM as explained above. Then, the only free parameter that we are left with to identify LAEs is the relative escape fraction of Lyα and continuum photons, i.e. f_α/f_c. This final free parameter is fixed by matching the simulated LFs to the observed ones, as explained in Section 5.

Before proceeding to a comparison between the theoretical results and observational data, we briefly digress to show the relation between T_α and |$\langle \chi _{\rm H\,\small{I}} \rangle$| for our five f_esc values. As can be seen from Fig. 4 (see also Section 4.1), the mean IGM transmission (〈T_α〉) averaged over all galaxies increases with decreasing |$\langle \chi _{\rm H\,\small{I}} \rangle$| values: as ionized regions grow larger, they allow all galaxies to transmit more of their Lyα flux through the IGM. The inhomogeneity of the ionized regions also leads to a huge variation in T_α as seen from the 1σ dispersions in the same figure. Finally, our mean IGM transmission values are consistent with the results obtained by Jensen et al. (2013), Iliev et al. (2008) and Dijkstra et al. (2007) for f_esc = 0.05; we note that Lyα transmission calculations have not been performed for escape fractions higher than this value.

Figure 4.

Mean transmission of all galaxies (averaged over 48 LOS) in the simulation snapshot at z ≃ 6.7 as a function of the mean IGM hydrogen neutral fraction. Red (dash–dotted), magenta (long dashed), dark blue (short dashed), light blue (dotted) and green (long dash–dotted) lines show the relation for f_esc = 0.05, 0.25, 0.5, 0.75 and 0.95, and error bars indicate the 1σ dispersion.

Open in new tabDownload slide

5 JOINT CONSTRAINTS ON REIONIZATION TOPOLOGY AND DUST

From Sections 4.1 and 4.2, it is clear that f_esc affects the observed Lyα luminosity both through its effect on the intrinsic Lyα luminosity and through its effect on the ionization state of the IGM. Further, for a given f_esc and ionization state, the only free parameter left in our model is the escape fraction of Lyα photons compared to UV-continuum photons, f_α/f_c. In this section, we use the LAE data accumulated by Kashikawa et al. (2011) to simultaneously constrain and find best fits to our three parameters, f_esc, |$\langle \chi _{\rm H\,\small{I}} \rangle$| and f_α/f_c. We use pcrash outputs for |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.9,0.75,0.5,0.25,0.1,0.01$| and 10⁻⁴ for each of the five f_esc = 0.05, 0.25, 0.5, 0.75, 0.95 at z ≃ 6.7. For each such combination, we start from the special case of assuming f_α/f_c = 0.68 which corresponds to dust being homogeneously distributed in the ISM in Section 5.1, and progress to the more general scenario of varying f_α/f_c to best fit the observations, in Section 5.2.

5.1 A special case: homogeneous dust

As shown in equation (4), the intrinsic Lyα luminosity produced by a galaxy depends on the fraction of H i ionizing photons that ionize the H i in the ISM. It is thus naturally expected that |$L_\alpha ^{\rm int}$| decreases with increasing f_esc since this leads to fewer photons being available for ionization in the ISM. This behaviour can be seen from panel (a) of Fig. 5, where we find that the amplitude of the Lyα LF drops by a factor of about 40 as f_esc increases from 0.05 to 0.95. This drop seems surprising at the first glance: as f_esc increases by a factor of about 19, the intrinsic Lyα LF would be expected to drop by the same amount. However, this result can be easily explained by the fact that we identify galaxies as LAEs based on the Lyα luminosity being larger than 10⁴² erg s⁻¹ and an EW ≥ 20 Å. As f_esc increases, fewer and fewer galaxies are able to fulfil these selection criteria, leading to an enhanced drop in the amplitude of the Lyα LF.

Figure 5.

The cumulative Lyα LFs. Panel (a) shows the intrinsic Lyα LFs for escape fractions f_esc = 0.05, 0.25, 0.50, 0.75, as marked. The lines in the other panels show the results for homogeneously distributed SNII dust (f_α/f_c = 0.68) for different |$\langle \chi _{\rm H\,\small{I}} \rangle$| states (⁠|$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.90$|⁠, 0.75, 0.50, 0.25, 0.10, 0.01, 10⁻⁴) obtained using the fixed f_esc value marked above the panel. In each panel, black open circles show the Lyα LFs at z ≃ 6.5 inferred observationally by Kashikawa et al. (2011). Due to the EW selection criterion, no LAEs are identified for f_esc = 0.95.

Open in new tabDownload slide

The ratio f_α/f_c depends on many physical effects (see e.g. Laursen et al. 2013), including the distribution of the Lyα sources in the dusty ISM (Scarlata et al. 2009), the ISM H i column densities and velocities. However, a natural consequence of simulating cosmological volumes is that we are unable to resolve the ISM of individual galaxies. We therefore start with the special scenario wherein dust is homogeneously distributed in the ISM and Lyα photons are not scattered by dust and H i within the ISM. Since SNII are assumed to be the main sources of dust in our model, we deduce the ratio f_α/f_c = 0.68 purely from the corresponding SN extinction curve (Bianchi & Schneider 2007). This fixed ratio is then applied uniformly to all galaxies in the simulation snapshot at z ≃ 6.7. Once f_α has been fixed using this ratio, the only two parameters that can affect the slope and amplitude of the Lyα LFs are f_esc and |$\chi _{\rm H\,\small{I}}$|⁠.

We start by discussing the Lyα LF evolution as a function of |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠. As expected from equation (6), for a given f_esc, τ_α decreases with decreasing |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠, leading to an increase in the IGM Lyα transmission. This naturally results in an increase in the amplitude of the Lyα LF by making galaxies more luminous in the observed Lyα luminosity, |$L_\alpha ^{\rm obs}$|⁠, as can be seen from panels (b)–(e) of Fig. 5. However, the effect of decreasing |$\chi _{\rm H\,\small{I}}$| is the strongest on the visibility of the faintest galaxies and decreases with increasing luminosity. This is because the spatial scale imposed by the Gunn–Peterson damping wing on the size of the H ii region corresponds to a separation of about 280 kpc (physical) at z ≃ 7 (Miralda-Escude 1998). The most massive galaxies are easily able to carve out H ii regions of this size as a result of their relatively large H i ionizing photon output. However, in the initial stages of reionization, only those faint galaxies that are clustered are able to build H ii regions of this size and become visible as LAEs (McQuinn et al. 2007b; Dayal et al. 2009). Hence, for a given f_esc, the faint-end slope of the Lyα LF becomes steeper with decreasing |$\langle \chi _{\rm H\,\small{I}} \rangle$| values as seen from panels (b)–(e) of Fig. 5. We also see that the slope and amplitude of the Lyα LFs start converging for |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 10^{-2}$|⁠; this corresponds to the average H i fraction at which galaxies typically inhabit H ii regions such that the red part of the Lyα line escapes unattenuated by H i.

We now discuss the evolution of the Lyα LFs with f_esc. As we have already noted in Fig. 3, although the spatial distribution of the ionization fields looks very similar for the varying f_esc values for a given |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠, the |$\chi _{\rm H\,\small{I}}$| value in any given cell decreases with increasing f_esc. This leads to an increase in T_α with f_esc such that for |$\langle \chi _{\rm H\,\small{I}} \rangle$| ≃ 0.1, averaged over all LAEs, <T_α> = (0.43, 0.45, 0.48, 0.59) for f_esc = (0.05, 0.25, 0.50, 0.75). However, this slight increase in T_α is compensated by the decrease in the intrinsic Lyα luminosity with increasing f_esc as shown in panel (a) of Fig. 5. As a result, the amplitude of the resulting Lyα LFs decreases slightly as f_esc increases from 0.05 to 0.5 for a given |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠; however, the slopes are very similar since f_esc affects the intrinsic luminosity of all simulated galaxies by the same factor. The Lyα LF drops steeply between f_esc = 0.5 and 0.75, and there are no simulated galaxies that would be identified as LAEs for f_esc = 0.95. This steep drop for f_esc > 0.5 arises due to our imposed criterion of the observed Lyα EW being larger than 20 Å. For a continually star-forming galaxy with an age of about 200 Myr, a stellar metallicity of Z_* = 0.02 Z_⊙ and f_esc = 0, the intrinsic Lyα EW has a value of about 114 Å; we note that this intrinsic EW would decrease for increasing values of Z_* and f_esc. For a value of f_esc = (0.75, 0.95), this intrinsic EW then reduces to about (28.5, 5.7) Å. Then, taking into account the effect of the IGM transmission, very few (none) of our simulated galaxies fulfil the EW selection criterion for f_esc = 0.75 (0.95) leading to a steep drop in the number density of LAEs; the total number density of LAEs drops to zero for f_esc = 0.95.

Finally, we note that both the slope and amplitude of the simulated Lyα LFs are in agreement with observations for f_esc = 0.05 and |$\langle \chi _{\rm H\,\small{I}} \rangle \le 0.1$|⁠; for f_esc ≥ 0.25, the simulated Lyα LFs are underestimated with respect to the data, even for a fully ionized IGM.

The UV LFs are not directly dependent on the IGM ionization state. However, both their slope and amplitude are affected by |$\langle \chi _{\rm H\,\small{I}} \rangle$| since these are the UV LFs for the simulated galaxies that would be identified as LAEs on the basis of |$L_\alpha ^{\rm obs}$| and the EW. It is interesting to note that with our deeper Lyα luminosity limit of |$L_\alpha ^{\rm obs}\ge 10^{42}$| erg s⁻¹ compared to the observational limit of 10^42.4 erg s⁻¹, we find that the cumulative UV LFs keep rising until M_UV ≤ −19 for f_esc ≤ 0.5 as seen from panels (a)–(c) of Fig. 6. From the same panels, we find that the bright-end slope of the UV LF (M_UV ≤ −21.5) matches the observed UV LFs for all values of |$\langle \chi _{\rm H\,\small{I}} \rangle$| and f_esc ≤ 0.5. Further, the UV LFs start converging at much higher values of |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.25$|⁠, compared to the value of 0.1 required for the Lyα LFs to settle. Similar to the behaviour seen for the Lyα LF, for a given f_esc ≤ 0.5, the faint end of the UV LF steepens with decreasing |$\langle \chi _{\rm H\,\small{I}} \rangle$| as more and more of these galaxies are able to transmit enough of their Lyα flux to be visible as LAEs. Again, for a given |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠, while the amplitude of the UV LFs decreases slightly between f_esc = 0.05 and 0.5, the slope remains the same. The steep drop in the number of LAEs for f_esc ≃ 0.75 naturally leads to a drop in the UV LF shown in panel (d) of the same figure; as for the Lyα LF, the UV LF amplitude is zero for f_esc = 0.95 due to no simulated galaxies being classified as LAEs.

Figure 6.

The cumulative UV LFs for simulated LAEs using a homogeneous dust case (f_α/f_c = 0.68). The lines in the panels show the UV LFs for different |$\langle \chi _{\rm H\,\small{I}} \rangle$| states (⁠|$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.90$|⁠, 0.75, 0.50, 0.25, 0.10, 0.01, 10⁻⁴) obtained using the fixed f_esc value marked above the panel. In each panel, black open circles show the LAE UV LFs at z ≃ 6.5 inferred observationally by Kashikawa et al. (2011). The vertical grey dotted lines indicate the limiting 2σ magnitude in the z′ band (see Kashikawa et al. 2011); measurements for fainter magnitudes may be uncertain because the corresponding z′-band magnitudes are no longer reliable. Due to the EW selection criterion, no LAEs are identified for f_esc = 0.95.

Open in new tabDownload slide

Finally, we note that by jointly reproducing the observed Lyα and UV LFs, we can constrain for f_esc = 0.05 one of the most important parameters for reionization: |$\langle \chi _{\rm H\,\small{I}} \rangle \le 0.1$| or 〈T_α〉_LAE ≥ 0.43.

However, we caution the reader that this result is only valid in the case where the ISM dust is homogeneously distributed. Since using only the UV LFs allows a much larger parameter range of |$\langle \chi _{\rm H\,\small{I}} \rangle \le 0.25$|⁠, a question that arises is whether this larger parameter space could also be reconciled with the Lyα LFs given clumped dust that could increase |$L_\alpha ^{\rm obs}$| by increasing f_α. We carry out these calculations in Section 5.2 that follows.

5.2 Clumped dust

A number of works have shown that inhomogeneously distributed (clumped) dust in the ISM can enhance the escape fraction of Lyα photons compared to UV photons which are not resonantly scattered by H i (Neufeld 1991; Hansen & Oh 2006). However, the amount by which the ratio f_α/f_c can be enhanced by clumped dust sensitively depends on the amount of dust that is bound in cold neutral gas clouds and the number of these encountered by Lyα photons within the galaxy. Since such sub-grid parameters cannot be resolved by our cosmological simulation, we use f_α/f_c as a free parameter and explore the range required to reproduce the data for the different f_esc and |$\langle \chi _{\rm H\,\small{I}} \rangle$| combinations discussed in Section 5.1.

As shown in Section 5.1, an increasing f_esc reduces the intrinsic Lyα luminosity for all galaxies, leading to a drop in the amplitude of the Lyα LF. However, since |$L_{\alpha }^{\rm obs}\propto (1-f_{\rm esc})f_{\alpha } T_\alpha$|⁠, a decrease in the intrinsic Lyα luminosity (due to an increasing f_esc) can be compensated by a larger fraction of Lyα photons escaping out of the galaxy (i.e. a higher f_α) and/or a larger amount being transmitted through the IGM (i.e. a larger T_α). Hence, for a given |$\langle \chi _{\rm H\,\small{I}} \rangle$| (i.e. T_α) by allowing values of f_α/f_c > 0.68, the amplitude of the Lyα LF can be boosted. However, it is important to note that since we use a constant value of f_α/f_c for all galaxies, it only changes the amplitude of the Lyα LF at the bright end; raising this ratio makes the faint-end slope slightly steeper since a number of galaxies which are at our detection limit (L_α ≥ 10⁴² erg s⁻¹ and EW ≥ 20 Å) become visible as LAEs with increasing f_α/f_c.

It is interesting to note that even using f_α/f_c ratios such that f_α saturates to 1, the simulated Lyα LFs are well below the observations for |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.75$|⁠, for all f_esc = 0.05 to 0.95, as seen from panels (a)–(e) of Fig. 7. This is because the H ii regions built by faint galaxies are too small to allow a large IGM Lyα transmission, even for f_esc = 0.95. We now discuss the results for varying f_esc: as mentioned in Section 5.1, the Lyα and UV LFs can be well reproduced using the homogeneous dust case for f_esc = 0.05 for |$\langle \chi _{\rm H\,\small{I}} \rangle \le 0.1$|⁠. In the case of clumped dust, the simulated LFs match the observations for f_esc = 0.05 for |$\langle \chi _{\rm H\,\small{I}} \rangle$| as high as 0.5, provided that the decrease in T_α is compensated by f_α/f_c increasing to 1.2 (see Table 1). For f_esc = 0.25, 0.5, within a 1σ error the simulated LFs can be reconciled with the observations for |$\langle \chi _{\rm H\,\small{I}} \rangle \le 0.5$|⁠, given that f_α/f_c increases to 1.6, 1.8, respectively, as shown in Table 1; however, due to a decrease in the intrinsic Lyα luminosity and EW, the faint-end slope of the simulated LFs is always underestimated with respect to the data for |$\langle \chi _{\rm H\,\small{I}} \rangle \ge 0.5$|⁠. As f_esc increases to 0.75 and above, the amplitude of the Lyα LFs cannot be boosted up to the observed value, even for f_α = 1; although the H ii regions built by these galaxies are very large as a result of the large H i ionizing photon escape fraction leading to large T_α values, the EW of the Lyα line prevents most of these galaxies from being classified as LAEs (see also Section 5.1 above).

Figure 7.

The cumulative Lyα LFs. The panels show the best-fitting simulated results for the f_α/f_c values shown in Table 1. The panels show the results for different |$\langle \chi _{\rm H\,\small{I}} \rangle$| states (⁠|$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.90$|⁠, 0.75, 0.50, 0.25, 0.10, 0.01, 10⁻⁴) obtained using the fixed f_esc value marked above the panel. In each panel, black open circles show the Lyα LFs at z ≃ 6.5 inferred observationally by Kashikawa et al. (2011).

Open in new tabDownload slide

We therefore find a three-dimensional degeneracy between the IGM ionization state |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠, the escape fraction of H i ionizing photons (f_esc) and the clumped/inhomogeneous distribution of dust in the ISM of galaxies (governing f_α/f_c) such that a decrease in |$L_\alpha ^{\rm int}$| with increasing f_esc can be compensated by a larger T_α and f_α/f_c. However, we find that the |$\langle \chi _{\rm H\,\small{I}} \rangle$| values which reproduce the observations decrease with increasing f_esc because an increasing f_α cannot compensate for the decrease in the intrinsic Lyα luminosity. This degeneracy can be seen from Fig. 8, where we show contours of 1σ-5σ deviations from the observations (Kashikawa et al. 2011), obtained by computing the χ² values for all the combinations of f_esc = (0.05, …, 0.95), |$\langle \chi _{\rm H\,\small{I}} \rangle =(10^{-4},\ldots , 0.90)$| and f_α/f_c = (0.60, …, 1.8), as shown in Table 1.

Figure 8.

1σ-5σ probability contours (black to light grey, respectively) for the combinations of |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠, f_esc and f_α/f_c shown in Table 1 that best fit the observed Lyα LF data at z ≃ 6.5 (Kashikawa et al. 2011). Within a 1σ error, we cannot distinguish between |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 10^{-4}$| to 0.5, f_α/f_c from 0.6 to 1.8 and f_esc ranging between 0.05 and 0.5. There exists a degeneracy between the IGM ionization state, escape fraction of H i ionizing photons and the dust clumping inside high-redshift galaxies such that a decrease in the IGM Lyα transmission can be compensated by an increase in the intrinsic Lyα luminosity produced, and the fraction of this luminosity that can escape the galaxy.

Open in new tabDownload slide

The relation between |$\langle \chi _{\rm H\,\small{I}} \rangle$| and 〈T_α〉 is complex due to the dependency upon the reionization topology as shown in Dijkstra et al. (2011; see also Fig. 4 and Jensen et al. 2013) and the assumed line shape (e.g. Gaussian, double-peaked; Laursen, Sommer-Larsen & Razoumov 2011; Jensen et al. 2013; Duval et al. 2014). For this reason, we also show the three-dimensional degeneracy in terms of the mean transmission across all LAEs (〈T_α〉_LAE) for a given combination of f_esc and f_α/f_c in Fig. 9.

Figure 9.

1σ-5σ probability contours (black to light grey, respectively) for the combinations of |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠, f_esc and f_α/f_c shown in Table 1 that best fit the observed Lyα LF data at z ≃ 6.5 (Kashikawa et al. 2011). Within a 1σ error, we cannot distinguish between 〈T_α〉_LAE ≃ 0.38 to 0.5, f_α/f_c from 0.6 to 1.8 and f_esc ranging between 0.05 and 0.5. There exists a degeneracy between the IGM Lyα transmission, escape fraction of H i ionizing photons and the dust clumping inside high-redshift galaxies such that a decrease in the IGM Lyα transmission can be compensated by an increase in the intrinsic Lyα luminosity produced, and the fraction of this luminosity that can escape the galaxy.

Open in new tabDownload slide

We note that for a given |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠, the mean transmission of all LAEs 〈T_α〉_LAE decreases for increasing f_α/f_c: with rising f_α/f_c more galaxies, with lower values of T_α, are identified as LAEs, leading to a decline in 〈T_α〉_LAE. This effect can be seen for the maximum 〈T_α〉_LAE values, within a 4σ error, in each panel of Fig. 9, and it demonstrates the dependency of 〈T_α〉_LAE on f_α/f_c.

From Figs 8 and 9, we see that within a 1σ error, the observed Lyα and UV LFs can be reproduced for f_esc values ranging between 0.05and0.5 and |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 10^{-4}{\rm -}0.5$| or 〈T_α〉 ≃ 0.38-0.5, provided that the decreasing intrinsic Lyα luminosity (with increasing f_esc) and decreasing IGM transmission (with increasing |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠) are compensated by an increasing f_α/f_c = 0.6 to 1.8. This implies that if dust is indeed clumped in the ISM of high-redshift galaxies such that the relative Lyα escape fraction can be as high as 1.8 times the UV escape fraction, we cannot differentiate between a universe which is either completely ionized or half neutral (or where 38–50 per cent of LAE Lyα radiation is transmitted through the IGM) or an f_esc ranging between 5 and 50 per cent. This parameter space is much larger than that allowed for a homogeneous dust distribution, where matching to the observations requires f_esc = 0.05 and |$\langle \chi _{\rm H\,\small{I}} \rangle \le 0.1$| (〈T_α〉_LAE ≥ 0.43).

6 CONCLUSIONS

In this work, we couple state-of-the-art cosmological simulations run using gadget-2 with a dust model and an RT code (pcrash) to build a physical model of z ≃ 6.7 LAEs to simultaneously constrain the IGM reionization state, the escape fraction of H i ionizing photons and ISM dust distribution (homogeneous/clumped).

We start by validating the properties of the simulated galaxy population against observations of high-redshift (z ≃ 6–8) LBGs. Identifying bound structures using the ahf, simulated galaxies which are complete in the halo mass function (M_h ≥ 10^9.2 M_⊙) contain a minimum of 4N gas particles (N = 40) and a minimum of 10 star particles are identified as the ‘resolved’ population. For each such resolved galaxy, we use the mass, age and metallicity of each of its star particles to build the composite spectra. We use a dust model (Dayal et al. 2010) that calculates the SNII dust enrichment depending on the entire SF history of a given galaxy, to obtain the associated escape fraction of UV photons. Once these calculations are carried out, galaxies with a dust-attenuated UV magnitude M_UV ≤ −17 are identified as LBGs at z ≃ 6–8. Comparing the predictions by the simulation against LBG observations, we find that the slope and amplitude of the dust-attenuated UV LFs are in good agreement with the data at z ≃ 6, 7 while the UV LF at z ≃ 8 requires no dust to match the observations, such that f_c ≃ (0.5, 0.6, 1.0) at z ≃ (6, 7, 8), averaged over all LBGs. The theoretical Schechter parameters are found to be α = (−1.9 ± 0.2, −2.0 ± 0.2, −1.8 ± 0.2) and M_{UV, *} = ( − 19.8 ± 0.3, −19.7 ± 0.2, −19.9 ± 0.3) at z ≃ (6, 7, 8), and are consistent with those inferred observationally (McLure et al. 2009, 2013; Bouwens et al. 2011). We also find that the theoretical LBG stellar mass functions, SMDs and sSFRs are in excellent agreement with the data at z ≃ 6–8 (see Appendix A).

Once that the physical properties and dust enrichment of the simulated galaxy population have been validated against the mentioned observational data sets, we proceed to modelling LAEs. This requires further information on (a) the fraction of the H i ionizing photons produced by a galaxy that are absorbed by ISM H i (1 − f_esc) and contribute to the Lyα emission line; the rest (f_esc) emerge out and contribute to reionization, (b) the relative effect of ISM dust on Lyα and UV photons (f_α/f_c), and (c) the level to which the IGM is ionized by the galaxy population (⁠|$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠), in order to calculate the IGM Lyα transmission. Since f_esc remains poorly constrained with proposed values ranging between 1 and 100 per cent (e.g. Ricotti & Shull 2000; Gnedin et al. 2008; Ferrara & Loeb 2013), we use five different values of f_esc = 0.05, 0.25, 0.5, 0.75, 0.95 to run the RT code (pcrash) in order to obtain various ‘reionization’ states of the simulated volume as it evolves from being fully neutral (⁠|$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 1$|⁠) to completely ionized (⁠|$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 10^{-4}$|⁠); pcrash follows the time evolution of the ionization fractions (H i and He ii, He iii ) and temperature given the spectra of the source galaxy population. Once these calculations are carried out, galaxies with |$L_{\alpha }^{\rm obs}\ge 10^{42}$| erg s⁻¹ and EW ≥ 20 Å are identified as LAEs; the only free parameter of the model (f_α/f_c) is then fixed by matching the simulated LFs to the observations.

We start from the simplest case of considering the SNII produced dust to be homogeneously distributed in the ISM. The SNII extinction curve yields f_α/f_c = 0.68. Using this model, we find that for a given f_esc, the IGM Lyα transmission increases with decreasing |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠, with this increase being most important for the smallest galaxies; as |$\langle \chi _{\rm H\,\small{I}} \rangle$| decreases, even small galaxies are able to ionize a large enough H ii region around themselves to be visible as LAEs. Hence, the Lyα LF becomes steeper with decreasing |$\langle \chi _{\rm H\,\small{I}} \rangle$| for any given f_esc value. Further, for a given |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠, although the ionization fields look very similar for different f_esc values, the degree of ionization (⁠|$\chi _{\rm H\,\small{I}}$|⁠) in any cell increases with f_esc due to an increase in the H i ionizing photon emissivity. This leads to a slight increase in T_α with f_esc for a given |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠. However, this increase is more than compensated by the decrease in the intrinsic Lyα luminosity with increasing f_esc (as less photons are available for ionizations in the ISM), as a result of which the Lyα LF decreases slightly in amplitude with f_esc for a given |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠. Finally, comparing simulated Lyα and UV LFs with observations, we find that for a homogeneous ISM dust distribution and f_esc = 0.05, we can constrain |$\langle \chi _{\rm H\,\small{I}} \rangle \le 0.1$| or 〈T_α〉_LAE ≥ 0.43.

We also explore the scenario of dust being inhomogeneously distributed/clumped in the ISM of high-redshift galaxies, boosting up the ratio of f_α/f_c. For a given |$\langle \chi _{\rm H\,\small{I}} \rangle$| value, an increase in f_esc leads to a drop in the amplitude of the Lyα LF due to a decrease in |$L_\alpha ^{\rm int}$|⁠, as mentioned above. In the clumped dust scenario, this decrease can be compensated by an increase in f_α. For |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.5$|⁠, f_esc values as high as 0.5 can be reconciled with the data, given that f_α/f_c increases from ≃ 0.6 to 1.8 as f_esc increases from 0.05 to 0.5. We thus find a three-dimensional degeneracy between |$\langle \chi _{\rm H\,\small{I}} \rangle$|⁠, f_esc and f_α/f_c such that a decrease in |$L_\alpha ^{\rm int}$| with increasing f_esc can be compensated by a larger T_α and f_α/f_c. Due to the model dependence (e.g. Lyα line profile, reionization topology) of the |$\langle \chi _{\rm H\,\small{I}} \rangle$|–〈T_α〉 relation, we also rephrase our result by converting the |$\langle \chi _{\rm H\,\small{I}} \rangle$| values to the corresponding mean transmission values for LAEs (〈T_α〉_LAE). Within a 1σ error, we find that allowing for clumped dust, the data can be reproduced for a wide parameter space of 〈T_α〉_LAE ≃ 0.38-0.50 or |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.5{\rm -}10^{-4}$|⁠, f_esc ≃ 0.05-0.5 and f_α/f_c ≃ 0.6-1.8, i.e. if dust is indeed clumped in the ISM of high-redshift galaxies such that the relative Lyα escape fraction can be as high as 1.8 times the UV escape fraction, we cannot differentiate between a universe which is either completely ionized or half neutral or where on average 38–50 per cent of the Lyα radiation is transmitted, or has an f_esc ranging between 5 and 50 per cent. We caution the reader that although the constraints in terms of 〈T_α〉_LAE account for the model-dependent relation between |$\langle \chi _{\rm H\,\small{I}} \rangle$| and 〈T_α〉, the value of 〈T_α〉_LAE sensitively depends on the galaxies that are classified as LAEs for a given combination of f_esc, |$\langle \chi _{\rm H\,\small{I}} \rangle$| and f_α/f_c.

We now discuss some of the caveats involved in this work. First, an enhancement of the escape fraction of Lyα photons with respect to continuum photons via the Neufeld effect (Neufeld 1991) in a clumpy ISM remains controversial: on the one hand, Hansen & Oh (2006) showed that such enhancement was possible if the Lyα photons encountered a certain number of ‘clumps’. On the other hand, Laursen et al. (2013) have shown that enhancing the ISM Lyα escape fraction requires very specific conditions (no bulk outflows, very high metallicity, very high density of the warm neutral medium and a low-density and highly ionized medium) that are unlikely to exist in any real ISM.

Secondly, as a consequence of simulating cosmological volumes, we are unable to resolve the ISM of the individual galaxies to be able to model the Lyα line profile that emerges out of any galaxy. We have therefore assumed the emergent Lyα line to have a Gaussian profile with the width being set by the rotation velocity of the galaxy. However, several observational (Tapken et al. 2007; Yamada et al. 2012) and theoretical works (Verhamme et al. 2008; Jensen et al. 2013) have found the emerging Lyα line profile to be double peaked; this profile results in higher IGM Lyα transmission values compared to a Gaussian profile (Jensen et al. 2013). An increase in T_α due to such a double-peaked profile allows our theoretical LFs to fit the observations with slightly lower f_α/f_c values.

Thirdly, we assume a constant escape fraction f_esc for all galaxies. Although the value of f_esc and its dependence on galaxy properties (such as the stellar mass) are only poorly known, recent work hints at an f_esc that decreases with increasing mass (Paardekooper et al. 2011; Ferrara & Loeb 2013). If f_esc is indeed higher for low-mass galaxies, the ionization topology would become more homogeneous with the H ii regions being built by low-mass (high-mass) galaxies being larger (smaller) than the constant f_esc case. This would have two effects: on the one hand, the IGM Lyα transmission would be enhanced (decreased) for low-mass (high-mass) galaxies, leading to a steeper Lyα LF. On the other hand, due to its direct dependence on f_esc, the intrinsic Lyα luminosity would decrease (increase) for low-mass (high-mass) galaxies, leading to a flattening of the Lyα LF. The relative importance of these two effects is the subject of a work currently in preparation.

Fourthly, we assume that dust in high-redshift galaxies is solely produced by SNII. However, as shown by Valiante et al. (2009), AGB stars can start contributing to significantly, and even dominate, the dust budget in regions of very high SF (e.g. QSO/quasars) after 150-500 Myr. However, a significant part of this extra dust could be destroyed by forward SN shocks (Bianchi & Schneider 2007) that have also not been considered in this work.

In ongoing calculations, our aim is to see whether the parameter space allowed by ‘clumped’ dust can be narrowed with a combination of LAE clustering and through 21 cm brightness maps that will be directly testable with upcoming observations with the Subaru Hyper-Suprime Camera and the LOw Frequency ARray (LOFAR).

The authors thank the referee M. Dijkstra for his insightful and constructive comments, Gustavo Yepes for providing the hydrodynamical simulation and Benedetta Ciardi for a collaboration in developing pcrash. This work has been performed under the HPC-EUROPA2 project (project number: 228398) with the support of the European Commission – Capacities Area – Research Infrastructures. AH thanks the European transnational HPC access programme for an eight-week visit to the Institute of Astronomy (Edinburgh) where a part of this work was carried out. PD acknowledges the support of the European Research Council. The hydrodynamical simulation has been performed on the Juropa supercomputer of the Jülich Supercomputing Centre (JSC).

REFERENCES

APPENDIX A: COMPARING THE SIMULATIONS WITH LBG OBSERVATIONS

A1 Stellar mass functions and SMD

We build simulated LBG mass functions by summing up LBGs on the basis of their stellar mass, and dividing this by the width of the bin and the simulated volume; for consistency with observations, this calculation is carried out for all LBGs with M_UV ≤ −18. As expected from the hierarchical model, small-mass systems are the most numerous, with the number density decreasing with increasing mass. Further, the mass function shifts to progressively lower masses with increasing redshifts, as fewer massive systems have had the time to assemble. Both trends can be seen clearly from Fig. A1: at z ≃ 6, galaxies with M_* ≃ 10^8.5 M_⊙ are two orders of magnitude more numerous as compared to more massive galaxies with M_* ≃ 10^10.5 M_⊙. Secondly, while there are systems as massive as 10^10.7 M_⊙ in the simulated volume at z ≃ 6, such systems had not the time to evolve by z ≃ 8 where the most massive systems have M_* ≃ 10^9.4 M_⊙. From the same figure, we see that the slope and amplitude of the simulated stellar mass functions are in agreement with the observations for M_* ≥ 10^7.5 M_⊙; this number marks the limit of our simulation resolution for galaxies containing at least 4N gas particles and 10 star particles, as a result of which, the number density drops below it.

Figure A1.

Stellar mass function for LBGs with M_UV ≤ −18 at z ≃ 6–8. The dotted (red), dashed (blue) and solid (green) lines represent the simulated stellar mass functions for z ≃ 6, 7 and 8, respectively. The shaded areas denote the respective Poisson errors. Squares (circles) represent the observational stellar mass functions at z ≃ 6 (7) inferred by González et al. (2011), corrected for completeness but not for dust.

Open in new tabDownload slide

As a result of the stellar mass functions shifting to progressively lower masses and number densities with increasing redshift (see Fig. A1), the total SMD in the simulated volume decreases with increasing redshift, as shown in Fig. A2, falling by a factor of about 40 from 10⁷ M_⊙ Mpc^{− 3} at z ≃ 6 to 10^5.4 M_⊙ Mpc^{− 3} at z ≃ 9.5. As expected from the agreement between the simulated stellar mass functions and those inferred observationally by González et al. (2011), the simulated SMD values are also in agreement with the observations. However, as pointed out by Schaerer & de Barros (2010) and Stark et al. (2013), including the effects of nebular emission can lead to a decrease in the observationally inferred stellar mass values; this effect causes a slight discrepancy between the SMD values obtained from our model and those inferred by Stark et al. (2013) at z ≃ 6, as shown in Fig. A2.

Figure A2.

Total SMD for LBGs with M_UV ≤ −18. Filled black points represent the simulated SMDs; Poissonian errors are too small to be visible on the plot. The other symbols show the observed SMD values inferred by González et al. (2011, empty red squares), Labbé et al. (2010b, empty green circles), Labbé et al. (2010a, empty green circles) and Stark et al. (2013, empty blue triangles).

Open in new tabDownload slide

Figure A3.

sSFR for LBGs with M_UV ≤ −18. Filled black points represent the simulated sSFRs; Poissonian errors are too small to be visible on the plot. Observed sSFR values plotted have been taken from González et al. (2010, empty red squares), Stark et al. (2009, empty green triangles), Yabe et al. (2009, empty magenta triangles), Schaerer & de Barros (2010, empty orange circles) and Stark et al. (2013, empty blue triangles).

Open in new tabDownload slide

A2 Specific star formation rates

The specific star formation rate (sSFR) is an excellent indicator of the current SFR compared to the entire past SF history of the galaxy. This quantity also has the advantage that it is relatively independent of assumptions regarding the IMF, metallicity, age and dust content since both the SFR and stellar mass are affected similarly by these parameters. Recently, a number of observational groups have suggested that the sSFR settles to a value consistent with 2-3 Gyr^{− 1} at z ≃ 3-8 (Stark et al. 2009; González et al. 2010; McLure et al. 2011; Labbé et al. 2013; Stark et al. 2013). This result is quite surprising given that theoretically, the sSFR is expected to trace the baryonic infall rate that scales as (1 + z)^2.25 (Neistein & Dekel 2008). However, Labbé et al. (2013) and Stark et al. (2013) have shown that at least a part of this discrepancy can be accounted for by the addition of nebular emission lines; indeed, the sSFR can rise by a factor of about 5 between z ≃ 2 and z ≃ 7 when nebular emission is taken into account.

Physically, the least massive galaxies tend to have the highest sSFR: as a result of their low M_* values, even a small amount of SF can boost their sSFR compared to that of more massive systems (see also Dayal & Ferrara 2012; Dayal et al. 2013). Since the mass function progressively shifts to lower masses with increasing redshifts, and smaller mass galaxies have larger sSFR values, the sSFR rises with increasing redshift (see Fig. A3), from sSFR ≃ 5 Gyr^{− 1} at z ≃ 6 to sSFR ≃ 12 Gyr^{− 1} at z ≃ 9.5, in accordance with other theoretical results (Biffi & Maio 2013; Dayal et al. 2013; Salvaterra et al. 2013). We note that in the redshift range of overlap (z ≃ 6-7), the theoretically inferred sSFR values are in agreement with observations.

APPENDIX B: COSMIC VARIANCE

In Fig. A4, we show the cosmic variance associated with the simulated Lyα LFs. The cosmic variance is estimated by computing the Lyα LFs for eight sub-volumes of our simulation volume, such that each sub-volume is 1.87 × 10⁵ Mpc³ and thus comparable to the volume observed by Kashikawa et al. (2011). As |$\langle \chi _{\rm H\,\small{I}} \rangle$| evolves from an IGM which is predominantly neutral at |$\langle \chi _{\rm H\,\small{I}} \rangle$| ≃ 0.75 to one wherein reionization is complete (⁠|$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 10^{-4}$|⁠), the value of T_α increases in all sub-volumes, leading to a decrease in the cosmic variance at the faint end of the LF; the variance at the bright end arises due to few objects being massive/luminous enough to occupy these bins in any given sub-volume.

Figure A4.

Simulated cumulative Lyα LFs and the associated cosmic variance for homogeneous dust (f_α/f_c = 0.68) and f_esc = 0.05. The panels show the results for different neutral hydrogen fractions, |$\langle \chi _{\rm H\,\small{I}} \rangle \simeq 0.90$|⁠, 0.75, 0.50, 0.25, 0.10, 0.01, 10⁻⁴, as marked. In each panel, the solid red line shows the Lyα LF for the whole simulation box of volume (80 h⁻¹ Mpc)³ with the grey shaded areas showing the associated cosmic variance obtained by splitting the entire simulation box into eight equal sub-volumes. In each panel, the open circles show the Lyα LFs at z ≃ 6.5 observed by Kashikawa et al. (2011). As the IGM becomes more ionized, the value of T_α increases in all sub-volumes, leading to a decrease in the cosmic variance at the faint end of the LF; the variance at the bright end arises due to the low numbers of massive galaxies found in any given sub-volume.

Open in new tabDownload slide