ON THE REDSHIFT EVOLUTION OF THE Lyα ESCAPE FRACTION AND THE DUST CONTENT OF GALAXIES

We now proceed to compile various estimates of f^Lyα_esc as a function of redshift, but first we present the formalism. We continue with the Hayes et al. (2010a) definition of f^Lyα_esc: the sample-averaged, "volumetric" escape fraction. This quantity is defined as the ratio of observed to intrinsic Lyα luminosity densities (ρ_L), derived by integration over luminosity functions (LFs), as in Equation (1):

$$\begin{equation} f_{\rm {esc}}^{\rm {Ly}\alpha } = \frac{\rho _{\rm {L,Ly}\alpha }^{\rm {Obs}}}{\rho _{\rm {L,Ly}\alpha }^{\rm {Int}}} = \frac{\int _{L_{\rm {lo}}}^\infty \Phi (L)_{\rm {Ly}\alpha }^{\rm {Obs}} \times L \times dL}{\int _{L_{\rm {lo}}}^\infty \Phi (L)_{\rm {Ly}\alpha }^{\rm {Int}} \times L \times dL}, \end{equation} \tag{ 1 }$$

where Φ(L) are the standard LFs.⁷ Thus, f^Lyα_esc is not simply a re-scaling of the LF by L (constantly scaling the escape fraction of all galaxies) or by Φ (the duty cycle scenario; see Nagamine et al. 2010 for examples of both of these methods). Instead, since f^Lyα_esc is simply defined as the ratio of luminosity densities, it can be thought of as the fraction of Lyα photons that escape from the survey volume, regardless of whether all galaxies show low f^Lyα_esc, or whether only a fraction of galaxies are in the Lyα-emitting phase with high individual f^Lyα_esc (see arguments in Tilvi et al. 2009). By this definition, f^Lyα_esc also includes any possible effect that the IGM may have on the measured Lyα emission from galaxies. However, it is clear that the bulk of the evolution of f^Lyα_esc with redshift found in this paper can clearly not be attributed to variations of the IGM transmission.

Where possible (i.e., z < 2.3) we make a direct comparison between Lyα and Hα. We apply the most appropriate dust correction to Hα and multiply by the case B recombination ratio of Lyα/Hα = 8.7 (Brocklehurst 1971) in order to obtain the intrinsic Lyα. That is,

$$\begin{equation} f_{\rm {esc}}^{\rm {Ly}\alpha } (z<2.3) = \frac{\rho _{\rm {L,Ly}\alpha }^{\rm {Obs}}}{8.7 \times \rho _{\rm {L,H}\alpha }^{\rm {Int}}} = \frac{\rho _{\rm {L,Ly}\alpha }^{\rm {Obs}}}{8.7 \times 10^{0.4 E_{B-V} k_{6563} } \cdot \rho _{\rm {L,H}\alpha }^{\rm {Obs}}}, \end{equation} \tag{ 2 }$$

where E_B−V must be the dust attenuation computed for the Hα-emitting sample, and k₆₅₆₃ is the extinction coefficient at the wavelength of Hα. Superscripts "Int" and "Obs" refer to the intrinsic and observed quantities.

At z ≳ 2.3, we are unable to obtain Hα LFs in order to use line ratios to estimate f^Lyα_esc, and instead the estimate is derived from the UV continuum. This is a less elegant method since the conversion between UV and Lyα requires the assumption of a metallicity, initial mass function (IMF), and evolutionary stage. However, in light of the fact that higher redshift Hα studies will remain impossible until the arrival of the James Webb Space Telescope, this is the only way to proceed. It is fortunate that there is no evidence that IMFs should differ between Lyα- and UV-selected populations, although metallicities have been shown to be around 0.2 dex lower (e.g., Cowie et al. 2010), which translates into a difference of ≲20% in the intrinsic Lyα/UV ratio (Leitherer et al. 1999). For "normal" metallicities and IMFs, and assuming that on average star formation is ongoing at equilibrium, this method is the same as taking the ratio of Lyα/UV star formation rate densities (SFRDs),$\dot{\rho }_{\star }$:

$$\begin{equation} f_{\rm {esc}}^{\rm {Ly}\alpha }(z>2.3) = \frac{\dot{\rho }_{\star, \rm {Ly}\alpha }^{\rm {Obs}}}{\dot{\rho }_{\star }^{\rm {Int}}} = \frac{\dot{\rho }_{\star,\rm {Ly}\alpha }^{\rm {Obs}}}{10^{0.4 E_{B-V} k_{\rm {UV}} } \times \dot{\rho }_{\star,\rm {UV}}^{\rm {Obs}}}, \end{equation} \tag{ 3 }$$

where now E_B−V must be the extinction seen by the UV-selected population and k_UV is the extinction coefficient in the UV.

The UV is of course not the only wavelength we can use for this experiment, but we choose to work exclusively with UV LFs since they (1) are so abundant in the literature, (2) have reasonably well understood selection functions, and (3) span an appropriately large range in redshift. We adopt UV measurements at redshifts most appropriate to our compiled Lyα data and dust attenuations derived from these samples themselves. We further adopt the dust attenuation law of Calzetti et al. (2000), and the star formation rate (SFR) calibrations of Kennicutt (1998). These calibrations assume a stabilized star formation episode at a constant rate for longer than around 100 Myr, with a Salpeter IMF (mass limits between 0.1 and 100 M_☉), and a complete ionization efficiency (no leaking and no destruction of ionizing photons by dust). In general, we assume that "UV" refers to the rest-frame wavelength of 1500 Å, where the extinction coefficient computed from the relationship of Calzetti et al. (2000) is 10.3. We want to emphasize that the definition of f^Lyα_esc we are using for high-redshift galaxies includes any effect that would decrease the number of observed Lyα photons with respect to the number expected from the SFR derived from the UV continuum level. The leaking of ionizing photons, as we will discuss later, would therefore imply an f^Lyα_esc value below unity, even if 100% of the Lyα photons effectively produced in the galaxy are able to escape without being affected by resonant trapping or destruction by dust.

The goal of this study is to determine the total volumetric escape fraction of a given volume, and ideally would include the very faintest systems. In practice, this would require the integration of the LFs down to zero, which, depending on the observational limits of a given survey and the redshift-dependent values of both L_⋆ and α, may include large extrapolations (or may even be divergent). It is therefore vital that our study employs lower integration limits that are (1) self-consistent between the populations, (2) include a sufficiently meaningful fraction of ρ_L, and (3) are not dominated by overextrapolation and uncertainties in the faint-end slope.

At z = 2, 3, and >4, several studies of the ρ_L,UV have been published, and here we adopt those of Reddy et al. (2008) and Bouwens et al. (2009), respectively. Both perform integrations down to 0.04L^z=3_⋆,UV and integrate to the same numerical lower limit at all redshifts. The lower limit is, of course, somewhat arbitrary, but is designed to find a reasonable medium between including a large fraction of the total luminosity/SFR density, and preventing (possible over–) extrapolation by integrating to zero. In this sense, it reflects the observational limits of the UV surveys.

Admitting that this number is somewhat arbitrary, we adopt the same approach and use 0.04L^z=3_⋆ − ∞ as the range for all of the integrations of the UV LF. For M^z=3_⋆ = −21.0 (AB), the corresponding lower luminosity limit is 4.36 × 10²⁷ erg s⁻¹ Hz⁻¹ (unobscured SFR = 0.6 M_☉ yr⁻¹). By adopting this limit, our results can easily be cross-checked against the available literature. At redshift 3 for the UV LF of Reddy et al. (2008), this range incorporates 70% of the luminosity, compared to an infinite integration under the LF.

Deciding upon a lower limit for the Hα LF is trickier, since it is difficult to know if we are extracting comparable samples of galaxies. There is no available z = 3 Hα LF, but if we adopt that compiled at z = 2.2 in Hayes et al. (2010b), and set the lower limit to 0.04 L_⋆, we obtain 4.6 × 10⁴¹ erg s⁻¹. This corresponds to much higher unobscured SFR than the lower UV limit at 3.5 M_☉ yr⁻¹. However, the UV and Hα-selection functions naturally recover galaxies of different dust contents; if we translate these limits to "true" SFRs for the respective samples, we obtain limits of 2.6 and 6.0 M_☉ yr⁻¹ for the UV and Hα, respectively. These limits differ by a factor of over two in SFR, but still are not able to account for the differing populations of galaxies that survive the respective selection functions; were the dustier galaxies that are selected by Hα able to enter the UV-selected catalogs, the increased average dust content would bring these values even closer together. We also argue that to some extent, the overall shape of the UV and Hα LFs must be governed by the same physical processes and, regardless of the exact dust content, selecting galaxies brighter than a certain fraction of the characteristic luminosity should recover objects with similar underlying SFRs. Ultimately this argument is backed up in Section 2.3 when we find very similar UV- and Hα-derived SFRs in the local universe, and by the very similar SFRDs derived by the two tracers in Reddy et al. (2008) and Hayes et al. (2010b). Naturally, by cutting both LFs at the same fraction of L_⋆, we recover similar fractions of the luminosity density compared with integration from zero (70%).

For Lyα, the situation is more complicated still: cutting at the same intrinsic SFR would mean that we do not include Lyα emission at lower luminosities. This is now not simply a matter of dust attenuation, but also includes radiation transport effects. Since we expect the line to be systematically weakened, applying a cut at the corresponding SFR to that of Hα or the UV would cause us to miss much of this light. The best way to proceed, therefore, is to adopt the same philosophy as above, and adopt 0.04 L^z=3_⋆,Lyα. By selecting the LF of Gronwall et al. (2007), we obtain a lower limit of 1.75 × 10⁴¹ erg s⁻¹. Should f^Lyα_esc = 1, this would correspond to an SFR of just 0.15 M_☉ yr⁻¹. However, in Hayes et al. (2010a) we determined a volumetric f^Lyα_esc of just 5%, and scaling this SFR up by a factor of 20 brings it to 2.9 M_☉ yr⁻¹, almost perfectly into line with the UV-derived 2.6 M_☉ yr⁻¹ discussed above. Naturally, this integration from 0.04 L_⋆ again includes ≈70% of the total luminosity density (compared with integrating from zero).

In summary, selecting the optimal integration limits is a non-trivial process; yet we argue that by adopting these limits we should be selecting very similar samples of galaxies, at least with respect to their unobscured SFR. The lower limits are 4.36 × 10²⁷ erg s⁻¹ Hz⁻¹ (UV), 4.6 × 10⁴¹ erg s⁻¹ (Hα), and 1.75 × 10⁴¹ erg s⁻¹ (Lyα). We have ensured that these limits include the bulk of the luminosity density but are not dominated in uncertainty by extrapolation in the faint end, although we have also confirmed that integration to zero in fact has only very minor effects on the final measurements of f^Lyα_esc.

All of the assembled data and the derived f^Lyα_esc measurements are summarized in Table 1 and Figure 1. The measurements of E_B−V relevant to each of the Hα or UV measurements are derived from data in the same publication as the Hα or UV LF data themselves (with one exception, which is discussed in the following paragraph). In this subsection, we provide the necessary motivation for our choices and comments on the various samples.

Zoom In Zoom Out Reset image size

Table 1. Lyα Escape Fractions with Redshift

Lyα Quantities			Intrinsic Quantities				Derived Results
z	Reference	$\dot{\rho }_{\star }$	z	Reference	EB−V	$\dot{\rho }_{\star }$	fLyαesc[% ]	Comment
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
Estimates based upon Lyα and Hα LFs.
0.2–0.35	De 08	(3.79 ± 1.69) × 10−4	0.2–0.35	TM 98	0.33	(0.0303 ± 0.017)	(1.25 ± 0.90)	1 mag at Hα
0.2–0.4	Co 10	(8.33 ± 2.60) × 10−5	0.2–0.35	TM 98	0.33	(0.0303 ± 0.017)	(0.275 ± 0.18)	1 mag at Hα
2.2	Ha 10	⋅⋅⋅	2.2	Ha 10	0.22	⋅⋅⋅	(5.3 ± 3.8)	Multi dimensional M.C.
Estimates based upon Lyα and UV LFs
2.5	Ca 11	(7.08 ± 0.81) × 10−3	〈2.3〉	Re 08	0.15	(0.201 ± 0.022)	(3.51 ± 0.56)
3.1	Gr 07	(8.50 ± 5.32) × 10−3	〈3.05〉	Re 08	0.14	(0.116 ± 0.017)	(7.33 ± 4.71)
3.1	Ou 08	(5.54 ± 2.91) × 10−3	〈3.05〉	Re 08	0.14	(0.116 ± 0.017)	(4.78 ± 2.61)
3.7	Ou 08	(4.78 ± 1.14) × 10−3	〈3.8〉	Bo 09	0.14	(0.089 ± 0.011)	(5.36 ± 1.43)
3.8	Ca 11	(8.71 ± 1.00) × 10−3	〈3.8〉	Bo 09	0.14	(0.089 ± 0.011)	(9.77 ± 1.64)
4.5	Da 07	(3.22 ± 1.25) × 10−3	〈4.7〉	Ou 04	0.075	(0.025 ± 0.011)	(12.6 ± 7.17)
4.86	Sh 09	(2.35 ± 3.17) × 10−3	〈4.7〉	Ou 04	0.075	(0.025 ± 0.011)	(9.24 ± 13.0)
5.65	Ca 11	(8.53 ± 3.44) × 10−3	〈5.9〉	Bo 09	0.029	(0.022 ± 0.005)	(38.1 ± 17.2)
5.7	Ou 08	(6.76 ± 4.77) × 10−3	〈5.9〉	Bo 09	0.029	(0.022 ± 0.005)	(30.2 ± 22.2)
6.6	Ou 10	(4.73 ± 1.24) × 10−3	6.5	Bo 07	0.012	(0.016 ± 0.008)	(30.0 ± 17.8)	UV Interpolated
7.0	Iy 06	(1.07 ± 1.16) × 10−3	7.0	Bo 09	0.010	(0.012 ± 0.008)	(8.96 ± 11.5)	UV Interpolated
7.7	Hi 10	(0+88.5−0) × 10−3	7.7	Bo 10	0.0	(0.005 ± 0.002)	(0+50.6−0)	UV Interpolated

Notes. For the Hα-based estimates, we use the integrated luminosity densities directly; SFRD measurements are presented just for homogeneity with the UV estimates. $\dot{\rho }_{\star }$ units are M_☉ yr⁻¹ Mpc⁻³ and E_B−V is in magnitudes. The references are expanded as: Bo 09 = Bouwens et al. (2009); Ca 11 = Cassata et al. (2011); Co 10 = Cowie et al. (2010); Da 07 = Dawson et al. (2007); De 08 = Deharveng et al. (2008); Gr 07 = Gronwall et al. (2007); Ha 10 = Hayes et al. (2010a); Hi 09 = Hibon et al. (2010); Iy 08 = Iye et al. (2006); Ou 04 = Ouchi et al. (2004); Ou 08 = Ouchi et al. (2008); Ou 10 = Ouchi et al. (2010); Sh 09 = Shioya et al. (2009); Re 08 = Reddy et al. (2008); TM 98 = Tresse & Maddox (1998). References for E_B−V measurements are the same as for the intrinsic SFRD (i.e., that listed in the fifth column) with the exception of the 〈z〉 = 0.3 points in which E_B−V is adopted from Kennicutt (1992).

Download table as:  ASCIITypeset image

No instrumentation can perform a Lyα-selected survey in the very nearby universe; so we begin at z ≈ 0.2–0.4 with the Lyα LFs presented in both Deharveng et al. (2008) and Cowie et al. (2010). At these redshifts Hα LFs are available, and therefore we proceed using Equation (2). We adopt the Hα LF of Tresse & Maddox (1998), and correct it for dust attenuation by applying the one magnitude of extinction that is representative of local Hα-selected galaxies (Kennicutt 1992). For security and consistency with higher redshift measurements, we also examine the z = 0.3 UV LFs of Arnouts et al. (2005), which we correct for dust using the method of Meurer et al. (1999) and the β slope measured by Schiminovich et al. (2005) in the same sample as Arnouts et al. (2005), finding extremely consistent numbers.

Beyond the very nearby universe, no further Lyα information is available before z = 2, where we adopt our own measurement of f^Lyα_esc=5.3% ± 3.8% (Hayes et al. 2010a), based upon Hα and individually estimated E_B−V.

It is already at this juncture in redshift that we lose the possibility of using Hα, and therefore we proceed using published UV LFs and Equation (3). Our next step is to take the Lyα LF of Cassata et al. (2011; 〈z〉 = 2.5) which we contrast against the dust-corrected ρ_L,UV of Reddy et al. (2008, 〈z〉 = 2.3). For this, and all subsequent points from Cassata et al. (2011), we adopt the values of L_⋆ that are uncorrected for IGM attenuation. It is re-assuring that the measurements at z = 2.2 and z ≈ 2.5 (which are based upon Hα and UV, respectively) give very consistent numbers. Furthermore, in a very recent submission (Blanc et al. 2010), an additional Lyα LF has been presented at 1.9 < z < 2.8, the integrated Lyα luminosity density which differs from our own result by ≈25%.

We then continue with the Reddy et al. (2008) UV data at 〈z〉 = 3.05, which we use to compute f^Lyα_esc for the z = 3.1 Lyα samples of Gronwall et al. (2007) and Ouchi et al. (2008).

At z ∼ 4, we have available Lyα LFs from Ouchi et al. (2008, z = 3.7) and Cassata et al. (2011, z = 3.9), and UV LFs from Bouwens et al. (2007, 〈z〉 = 3.8). We also use the z = 4.5 and 4.86 Lyα LF points from Dawson et al. (2007) and Shioya et al. (2009), which we normalize by the dust-corrected UV point at z = 4.7 from Ouchi et al. (2004).

The next redshift to examine is the popular z ≈ 5.7 Lyα window. Here, we adopt the UV data point from the i-dropout sample of Bouwens et al. (2007, 〈z〉 = 5.9), and the Lyα LF Ouchi et al. (2008, 〈z〉 = 5.7), which is in good agreement with those of Shimasaku et al. (2006), Ajiki et al. (2006), and Tapken et al. (2006). We also add the highest redshift LF from Cassata et al. (2011) at 〈z〉 = 5.65.

Finally, we assemble a few z>6 samples. We adopt the z = 6.5 point from Ouchi et al. (2010, which includes the sample of Kashikawa et al. 2006), and the measurement of Iye et al. (2006) at z = 7.0, which has also been compiled in Ota et al. (2008). Here, we adopt Bouwens et al. (2010b) UV measurement at 〈z〉 = 6.8 for comparison. It should be noted that at this redshift the dust-corrected and uncorrected measurements of Bouwens et al. (2010b) converge. We adopt the most optimistic estimate at z = 7.7 from Hibon et al. (2010), for which we interpolate between the Bouwens et al. (2010b) at z = 6.8 and 8.2 UV data points. Hibon et al. (2010) present Schechter parameters for four Lyα LFs, based upon various assumptions about the rate of contamination by lower redshift galaxies. By assuming all of their candidates are real (their sample a), we find a Lyα escape fraction of (33.5^+50.6_−33.5)%. We also briefly examine their subsample b, in which only four of the seven objects are real. For all of their subsamples, the numbers are insufficient to provide meaningful errors on the luminosity density and by our standard error procedure we derive f^Lyα_esc = (22.2⁺¹⁷⁰⁷_−22.2)%. Furthermore, it should be noted that in the Hibon et al. (2010) sample, the lower limits obtained on the Lyα equivalent width are in the range 6–15 Å, with their continuum-detected object showing W_Lyα = 13 Å. Thus, at an acceptable confidence limit, none of their seven objects would actually survive the canonical W_Lyα cut of 20 Å that is typically employed in narrowband surveys. Including these data is therefore not straightforward, but in order to treat them as consistently as possible with the lower redshift points, we have to set the z ≈ 7.7 Lyα escape fraction to zero, but adopted a characteristic error of 50.6% as derived from their most optimistic sample. We note that this limit is likely extremely high.

All of our measurements of f^Lyα_esc are listed in Table 1 and shown graphically in Figure 1, which is the main result of this paper.

It should always be borne in mind that we are compiling results from different survey teams, who may adopt different techniques for data reduction and photometry, derivation of the LFs, and incompleteness corrections. For example, Malhotra & Rhoads (2004) find reasonable agreement at z ≈ 5.7 between the narrowband-selected Lyα LFs of Rhoads & Malhotra (2001) and Ajiki et al. (2004), and the lensing-based survey of Santos et al. (2004). However, the z = 5.7 LF of Shimasaku et al. (2006), on which the study of Kashikawa et al. (2006) is based (see Section 5.2), find a strong disagreement at the faint end between their own LF and the compilation of Malhotra & Rhoads (2004). As commented by Shimasaku et al. (2006), the likely cause for this discrepancy lies in (1) the lack of incompleteness corrections, which are unmentioned in any of the 2004 articles, (2) the differences in equivalent-width-based selection criteria, and (3) the large cosmic variance which Ouchi et al. (2008) noted can be of factors of ≈2 in fields as large as 1 deg². Our results are sensitive to all of these considerations.

It is only now that sufficiently large samples of Lyα-emitting galaxies are presented in the literature for this study to be undertaken, and we are fortunate that a good fraction of our data must contain internal self-consistency. For example, three of our data points (at z = 3.1, 3.7, and 5.7) are drawn from a single paper (Ouchi et al. 2008) in which the methodologies must be internally consistent, and the basic trend can be seen in these data alone. A fourth point at z = 6.6 comes from Ouchi et al. (2010) where similar self-consistency is to be expected. In the same fashion, three further points are taken from Cassata et al. (2011) where internally the same methodology must have been adopted at each redshift. It is certainly encouraging that, for example at z = 5.7 the measurements of f^Lyα_esc based upon Cassata et al. (2011) and Ouchi et al. (2008) are practically indistinguishable, despite the fact that they are based upon completely different methods: blind spectroscopy and narrowband imaging, respectively. The z = 2.2 and 2.5 points of Cassata et al. (2011) and Hayes et al. (2010a) are similarly indistinguishable (and also robust against the same fundamental methodological difference of blind spectroscopy versus narrowband imaging), as are the z = 3.1 points of Ouchi et al. (2008) and Gronwall et al. (2007), both based on narrowband imaging.

Any study of the galaxy population benefits from targeting spatially disconnected, independent pointings in order to beat down cosmic variance. By adopting the studies of various authors pointed all over the extragalactic sky, this study is able to benefit from the inclusion of a large number of independent fields.

Figure 1 reveals a general and significant trend for f^Lyα_esc to increase with increasing redshift. Beginning in the very local universe we see $f_{\rm {esc}}^{\rm {Ly}\alpha } \sim 0.01$ or lower for nearby star-forming objects. This increases to around ≈5%–10% by redshift of ≈3–4, and further to ≈30%–40% by redshift 6. In order to quantify this trend we fit an analytical function to these data points, choosing a power law of the form f^Lyα_esc(z) = C × (1 + z)^ξ—we obtain coefficients of C = (1.67^+0.53_−0.24) × 10⁻³; ξ = (2.57^+0.19_−0.12). Note that we do not include any z>6 points in our fit since previous studies suggest that it is around this redshift that an appreciable fraction of the intergalactic hydrogen becomes neutral, and may in principle affect the Lyα LF. For more discussion on this see Section 5.2. To insure that the fit is not biased by the presence of two z ≈ 0.3 points that lie around 8 Gyr from z ≈ 2, we repeat the fit after excluding these points, finding C = (4.79^+5.68_−0.69) × 10⁻⁴; ξ = (3.38^+0.10_−0.37). Clearly, the fit is affected by these points, but their exclusion actually results in a more rapid evolution with redshift.

The apparent trend begins to break beyond redshift 6, but it is initially very slow. Over the redshift interval from 5.7 to 6.5, f^Lyα_esc stabilizes, but decreases again to just ≈10% at z = 7. The redshift 7 point from Iye et al. (2006) is confirmed, whereas none of the sample of redshift 7.7 candidates from Hibon et al. (2010) have confirmations by spectroscopy, and this upper error bar must be regarded as an optimistic upper limit.

Finally, we perform a simple experiment with the best-fit relationship to the f^Lyα_esc–z trend, and extrapolate to estimate the redshift at which f^Lyα_esc reaches unity. This would carry the implication that the ISM of the average galaxy has become effectively devoid of dust, and since dust is a byproduct of the star formation process, must also correspond to a time of approximately primeval star formation. It is interesting, therefore, that we find f^Lyα_esc = 1 at z = 11.1^+0.8_−0.6, which is consistent with the redshift of the instantaneous reionization of the universe based upon Wilkinson Microwave Anisotropy Probe data (z = 11 ± 1.4; Dunkley et al. 2009).

Naturally, this is not the first time that f^Lyα_esc has been estimated and several other studies based on a wide array of methods have attempted to pin down the same quantity at different redshifts. For example, at redshifts of 5.7 and 6.5, we compute f^Lyα_esc of around 40% and 30%, respectively. Based upon the fitting of spectral energy distributions (SEDs) to stacked broadband fluxes, Ono et al. (2010) estimate f^Lyα_esc=(36⁺⁶⁸₋₃₅)% and (4⁺¹⁸⁰_−3.8)% at the same redshifts. Although derived from an interesting approach, the uncertainties are still too large to provide a useful comparison.

Like us, Nagamine et al. (2010) compared observed Lyα LFs (Ouchi et al. 2008, in this case, which we also use) with intrinsic estimates, having derived this intrinsic LF from smoothed particle hydrodynamical (SPH) models of galaxy formation. They adopt two methods of scaling the intrinsic to the observed LFs, the first of which they call "escape fraction," which is a scaling to the data points along the luminosity axis, and assumes all galaxies have the same f^Lyα_esc. This method finds f^Lyα_esc = 10% at z = 3, which is certainly consistent with our estimates based on the z = 3.1 LF of Gronwall et al. (2007) and similar to but slightly higher than our estimate based on Ouchi et al. (2008). At z = 6, however, Nagamine et al. (2010) require an escape fraction of just 15% which is lower than our estimates of 30%–40%, and discrepant with our estimates at around the 2σ level. Nagamine et al. (2010) also test a "duty-cycle" scenario (an LF scaling along the Φ axis) in which only a fraction of the SPH galaxies are "on" as Lyα emitters, but emit 100% of their Lyα photons. Note that in these two extreme scenarios, there is no requirement for the integral over the scaled LF to be equivalent. Nagamine et al. (2010) present duty cycles of 0.07 and 0.2 at z = 3 and 6, respectively. However, before they compute these scalings the observed LFs are shifted along the luminosity axis by IGM attenuation factors of 0.82 (z = 3) and 0.52 (z = 6). This re-scaling needs to be removed before making a comparison with our estimate, and in the duty-cycle scenario the volumetric escape fractions that one would infer from the study of Nagamine et al. (2010) are 6% at z = 3 and 10% at z = 6. Again this agrees very well with our measurement at z ≈ 3, but compared with our estimates at z = 6 is an underestimate of around the same magnitude as their escape fraction method.

In contrast, using similar SPH galaxy formation models but modified prescriptions for Lyα production and transmission, as well as a different reionization history, Dayal et al. (2009) find Lyα escape fractions of 30% at both z = 5.7 and 6.5, which corresponds exactly with our measurements. Similar values of f^Lyα_esc ∼ 23%–33% have also been obtained in the followup work of Dayal et al. (2010), although they also include an IGM transmission of T_α = 0.48. Throughout this paper we have made sure not to apply any IGM correction, since the value of T_α remains poorly constrained, even theoretically, and from an observational perspective there is no strong evidence for exactly how close the IGM comes to a narrow Lyα line. As with the Madau (1995) prescription, it is likely that this IGM transmission is too low when considering lines that are systematically redshifted by the kinematics of the ISM, which would drive up these theoretical estimates of the Lyα escape fraction.

Adopting a similar method of LF scaling by luminosity, Le Delliou et al. (2005) found that an escape fraction of 2% was sufficient to match observed Lyα LFs with their predictions based upon semi-analytical models between z = 2 and 6, with the same machinery able to predict the clustering properties of Lyα emitters (Orsi et al. 2008). This is at the lower end of being consistent with our z = 3 measurements, and should the same escape fraction hold at z = 0.3, would also be consistent with our estimates in the nearby universe. However, the Le Delliou et al. (2005) escape fraction is highly inconsistent with our estimates at higher redshift. These semi-analytical models, using the prescription of Baugh et al. (2005), categorized star formation as occurring in two discrete modes, with a normal Salpeter IMF (α = −1.35) assigned to quiescent star formation and a flat IMF (α = 0) for bursting systems. This flat IMF increases the ionizing photon production at a given SFR by a factor of 10 and was implemented as a requirement in order to reproduce the population of submillimeter-selected galaxies at z>2. However, as noted by Le Delliou et al. (2006), the fraction of total star formation that occurs in bursts increases from 5% at z = 0 to over 80% at z = 6, and thus their model implies that by the z = 5.7 points; effectively all stars are formed in environments where ionizing photons are greatly overproduced compared with the present day. Should this requirement of the flat IMF be removed and the Salpeter IMF applied throughout, the intrinsic rate of production of ionizing photons would be decreased by a factor of three at z = 3.1 where the star formation is shared evenly between bursting and quiescent systems. This would bring the f^Lyα_esc estimate to 11% at this redshift. At z = 6, f^Lyα_esc = 16% would be found by replacing the flat IMF with the Salpeter IMF. These numbers are indeed very similar to the SPH models of Nagamine et al. (2010), but inconsistent with those of Dayal et al. (2009) and our own estimates based upon observation. It is interesting to point out, however, that the IMF assumption has little effect on the z ≈ 0.3 points where, in their model, the quiescent mode of star formation dominates.

The evolution in measured f^Lyα_esc is substantial, covering approximately two orders of magnitude, and no doubt holds vital information about the physical nature of galaxies at various cosmic epochs. As we will show in Section 5, the most likely explanation for this evolution is the decrease of the average dust content of galaxies. However, from a physical perspective many effects may enter. For example, galaxies may also contain less neutral hydrogen to scatter photons, show faster outflows, or become more clumpy. The inferred increase may alternatively be mimicked by galaxies becoming younger on average, having low and decreasing metallicities, or forming stars with IMFs that become more biased in favor of massive, ionizing stars. On the other hand, the scattering of Lyα photons by a neutral IGM and the general leakage of ionizing photons (Lyman continuum; LyC) are expected to increase with increasing redshift, and would both serve to lower the perceived Lyα escape fraction (although the "true" f^Lyα_esc of galaxies, i.e., before the IGM, would not be affected).

Regrettably, we are not able to measure any of these quantities directly from this compilation of data. We have, however, assembled data that show a number of trends with redshift: the Lyα and UV luminosity densities and the dust contents. These we have combined to show how f^Lyα_esc evolves; yet in order to extract the maximum of information from these, we need to examine another possible trend: how f^Lyα_esc correlates with the dust content. Thus, we delay a detailed discussion of what drives the f^Lyα_esc–z trend until Section 5 and now proceed to discuss the effects of radiation transport and dust absorption.

5.1.1. The Evolving Dust Content of Galaxies

We showed in the previous section that f^Lyα_esc of individual galaxies is anti-correlated with the measured E_B−V (Figure 2). Given that the typical E_B−V evolves with redshift (see Table 1), we may indeed expect a positive correlation between f^Lyα_esc and redshift. This is exactly what Figure 1 shows, where it is clear that the Lyα escape fraction increases smoothly and monotonically out to z ∼ 6. Thus, it appears that this increase in f^Lyα_esc is the result of the dust content of the star-forming galaxy population decreasing with redshift. We now take the measured values of E_B−V from the various samples (listed in Table 1), and use them to compute the f^Lyα_esc that would be expected from the three conversions between E_B−V and f^Lyα_esc discussed in the previous section (Calzetti et al. 2000; an empirical fit with one free parameter from Hayes et al. 2010a; and an empirical fit with two free parameters from this study). We show the measured escape fractions together with these predictions in Figure 3.

Zoom In Zoom Out Reset image size

We first discuss the predictions based upon the Calzetti et al. (2000, red dotted line), which is clearly discrepant with the observations at around the 3σ level at every redshift. Obviously, this is to be expected since Lyα photons resonantly scatter and it is unlikely that the dust is distributed in a uniform screen. The one-dimensional fit from Hayes et al. (2010a) offers substantial improvement and is able to describe the observations between redshifts 0 and 4. This reasoning is circular for the redshift 2 points where the f^Lyα_esc–E_B−V relationship was derived, but we stress the tautology is present only at this redshift. This relationship is not able to explain any of the data points at redshift above 4, where it systematically overpredicts the Lyα escape fraction.

As redshift increases, the dust content of galaxies is clearly shown to change, and could we plot Figure 2 at redshifts higher than 3, we would expect galaxies to cluster successively further toward the upper left corner of the plot. Since the Hayes et al. (2010a) f^Lyα_esc–E_B−V fit is forced through the (E_B−V,f^Lyα_esc) = (0,1) coordinate and a high value of k_Lyα is found, the predicted escape fraction evolves very quickly with redshift. Indeed, these predictions evolve much faster than the data, as f^Lyα_esc is forced for unphysical reasons toward unity.

When we introduce the new f^Lyα_esc–E_B−V fit with two free parameters and allow C_Lyα≠1, the agreement between the measured and observed Lyα escape fractions is striking: it agrees with essentially every data point, within the error bars, between redshifts 0 and 6.6. We should point out that it is not clear that the use of the average E_B−V for a sample should by necessity reproduce the volumetric escape fraction. Due to variations of the dust contents and ISM of individual galaxies, and the associated impact upon the transfer of Lyα and the selection of galaxies, it is plausible that the average f^Lyα_esc could have been skewed substantially from the data points. Indeed, close examination of the f^Lyα_esc–E_B−V relationship (C_Lyα = 1; k_Lyα=17.8) from Hayes et al. (2010a) reveals that it does not perfectly intersect the center of the f^Lyα_esc data point (z = 2.2 point in Figure 3) from the same survey, despite f^Lyα_esc, average E_B−V, and the coefficients of the f^Lyα_esc–E_B−V relationship all having been derived entirely from this one data set. This most likely results from the weighting across the population from which the average E_B−V is computed (the representative E_B−V is not an average weighted by the intrinsic Lyα luminosity), exactly the effect under discussion. However, the fact that such tight agreement is seen between the observational estimates and those derived from our fit suggests that such a bias in the selection of the populations is not at play here.

Again we stress that the relationship we derived between f^Lyα_esc and E_B−V in Section 4 includes the effects of resonance scattering, and thus in some manner the neutral gas content, its kinematics, and relative geometry all enter the relationship, which holds even when the measured optical color excess on the stellar continuum is zero. There is no reason to assume that these quantities are constant with redshift and we could, for example, envisage situations where the gas content, feedback properties, or clumpiness evolve and thereby change k_Lyα or C_Lyα. However, the tight agreement between our observed f^Lyα_esc values and those computed from the f^Lyα_esc-E_B−V relationship provides no evidence for the evolution of these properties (at least if the gas content does change it does not take part in the Lyα- scattering process). The evolution of f^Lyα_esc across almost the entire observable universe can be explained cleanly within the confines of this simple model, as mainly due to a dust content that evolves with redshift.

5.1.2. Other Effects

Beyond a redshift of around 5.7, the measured value of f^Lyα_esc begins to decline, although initially this decline is weak and the deviation from our best-fit relationship at redshift 6.5 is not significant. Adding the z = 6.5 point of Ouchi et al. (2010) and Kashikawa et al. (2006) to our fit does not change the result. However, the z = 7 point lies at just 8%, and is around 2σ below both the best-fit f^Lyα_esc–z relationship (Figure 1) and the predictions at this redshift based upon the f^Lyα_esc–E_B−V relationship (Figure 3). The z = 7.7 point formally takes the value of f^Lyα_esc = 0, and is presented with an extremely conservative error that is likely to be grossly overestimated (see Section 2.3). In comparison to z = 5.7, f^Lyα_esc has declined by a factor of at least two by z = 7. We have so far attributed the increase in f^Lyα_esc to an evolution in the dust content of galaxies, and it would be an extravagant departure from this evolutionary trend were ISM evolution to suddenly cause a sharp drop in f^Lyα_esc at z>6. Several other mechanisms are, however, naturally able to explain this break in the trend.

5.2.1. Leaking Ionizing Radiation

5.2.2. Neutralizing the Intergalactic Medium

So far, we have been taking advantage of the fact that we have measurements of the dust extinction in our samples of Hα- and UV-selected galaxies. We have used this to infer the intrinsic SFRD of the populations, and from there calculated f^Lyα_esc using Equations (2) and (3). These equations connect the quantity of f^Lyα_esc with E_B−V, via the ratio of the Lyα- and UV-derived measurements of the SFRD, $\dot{\rho }_{\star }$. However, in Section 4 we defined an alternative relationship between f^Lyα_esc and E_B−V, based upon analyzing individual galaxies, where we provide an empirical relationship between these two quantities (Equation (4)) and relate them simply through coefficients. Thus, we have four quantities ($\dot{\rho }_{\star }^{\rm Ly\alpha }$, $\dot{\rho }_{\star }^{\rm {UV}}$, E_B−V, and f^Lyα_esc) that are related by the various coefficients discussed in the previous sections.

In all the previous sections, we have made use of the measured values of E_B−V, but instead we could ignore this measurement and invert the problem: use the observed Lyα and uncorrected UV SFRDs at a given redshift to estimate E_B−V, using Equation (4) as a closure relation. Thus, substituting Equation (4) into Equation (2), we can write

$$\begin{equation} E_{B-V} = \frac{1}{0.4 (k_\lambda - k_{\rm {Ly}\alpha })} \times \log _{10} \left(\frac{\dot{\rho }_{\star,\rm {Ly}\alpha }^{\rm {Obs}}}{\dot{\rho }_{\star,\rm {UV}}^{\rm {Obs}} \times C_{\rm {Ly}\alpha }} \right). \end{equation} \tag{ 6 }$$

Out to z ≈ 6 we take the data compiled in Table 1, and compute the observed SFRD from either Hα or the UV, depending on the redshift. We then use Equation (6) to estimate the sample-averaged E_B−V at each redshift, independently of the attenuation measurements themselves. In short, we ignore the fact that these E_B−V measurements have been made, and see if we can recreate them. We show the result as black data points in the upper panel of Figure 4, with the actual measurements shown by the small gray symbols. We then hypothesize that the dust content of the universe may decrease exponentially, and adopting a function of the form E_B−V(z) = C_EBV ·  exp(z/z_EBV), we fit the coefficients C_EBV = 0.386 and z_EBV = 3.42. Or, the e-folding redshift scale for the E_B−V evolution is ≈3.4. We show this relationship in Figure 4 with the thick red line.

Zoom In Zoom Out Reset image size

By performing this experiment we are throwing away observational information and the plot becomes somewhat noisier, but nevertheless it resembles an inverted version of Figure 1. Fundamentally, the plot shows a decrease in the E_B−V of galaxies as redshift increases, which is consistent with the measurements (Table 1 and gray points). This decrease in the dust content of galaxies with redshift is already much discussed in the literature for LBGs at z ∼ 2–7 based upon a gradual blueing of the UV slopes (e.g., Hathi et al. 2008; Bouwens et al. 2009). At higher redshift, there is however a tendency for our new method to estimate higher E_B−V compared to the measurements obtained directly from the UV stellar continuum. In the upper panel of Figure 4, we also show the best-fitting relationships derived in the previous sections (black solid line), where we take the redshift evolution of f^Lyα_esc and use Equation (4) to convert to E_B−V using our best-fit coefficients—naturally, this line almost perfectly reproduces the gray points.

The red line (fit to these data) runs slightly flatter than the black one (combined fits from the previous sections) and suggests a slightly higher E_B−V, and therefore dust content, than measured at the highest redshifts in Bouwens et al. (2009). At z ∼ 3–6 however, it runs lower than the measurements of Hathi et al. (2008), who obtain slightly higher dust attenuations from LBG samples.

It is interesting to further investigate how the dust obscurations we derive compare with other estimates. The SPH modeling of Nagamine et al. (2010) and Dayal et al. (2009) already discussed in Section 3.2 both predict higher dust attenuations than measured in the z = 6 dropout populations at E_B−V = 0.15. Detailed SED modeling of z∼ 6–8 galaxies by Schaerer & de Barros (2010) also suggest the presence of dust in some high-z LBGs. Here, we estimate E_B−V ≈ 0.08 based upon the new methodology. Similarly the semi-analytical approach developed in Baugh et al. (2005) find E_B−V ∼ 0.1 at z>3 when examining the LBG population, which is certainly compatible with our estimates in the redshifts 3–5 domain.

Since the empirical relationship derived between f^Lyα_esc and E_B−V relates the two quantities directly, for completeness we convert our E_B−V redshift estimates to f^Lyα_esc through Equation (4). This enables us to approximately re-create the main observational result of this article, Figure 1, which we show in the lower panel of Figure 4. Here, we show the f^Lyα_esc estimates derived in this section with black shapes, with the original points from Figure 1 shown in gray. In short the difference between the two sets of points is that in the gray ones, the measured dust attenuation has been applied to the UV SFRD in the computation of f^Lyα_esc whereas in the black points, this quantity has been estimated directly from the observed SFRDs, using Equation (6). As with Figure 1 this shows f^Lyα_esc increasing with redshift, but the actual estimates of the dust attenuation in the individual samples have not been used in the derivation of this figure. The overall trend of Figure 1 is maintained, although significant scatter has been added to the plot. It shows that even were no E_B−V measurements available, our main result would have taken the same form and the overall trend would have been the same.