ABSORPTION-LINE SPECTROSCOPY OF GRAVITATIONALLY LENSED GALAXIES: FURTHER CONSTRAINTS ON THE ESCAPE FRACTION OF IONIZING PHOTONS AT HIGH REDSHIFT

The systemic redshift is an important measure since it enables us to examine the spectra for the presence of outflowing gas. Accurate systemic redshifts are usually obtained from stellar absorption lines or the nebular emission lines arising from H ii regions. For our sample, a redshift from nebular [O ii] λ3727 emission is only available for MS 1358 (Swinbank et al. 2009). However, since our spectra primarily target low-ionization absorption lines and Lyα emission, it is more practical to use these to estimate the systemic redshift (which is in reasonable agreement with that derived from nebular emission for MS 1358). Jones et al. (2012) determined an offset of Δv_IS = 190 km s⁻¹ between the systemic redshift (z_sys) determined from the C iii λ1176 stellar photospheric line and that measured using the low-ionization absorption lines (z_IS) using a composite spectrum of 81 LBGs at z ∼ 4. Following Jones et al. (2013), we apply these velocity offsets to the average weighted centroid of low-ionization absorption lines to determine the systemic redshift of each galaxy. We used all prominent low-ionization absorption lines available, namely, Si ii λ1260, O i λ1302, Si ii λ1304, C ii λ1334, and Si ii λ1527. The uncertainty is a square sum of the standard deviation of redshifts measured from different absorption lines, and the rms scatter is ∼125 km s⁻¹ between redshifts derived from absorption lines and those from nebular emission (Steidel et al. 2010). Together with the redshift of Lyα, the results are listed in Table 1.

We measure the SFRs of our lensed targets using the extinction-corrected 1600 Å flux from the observed DEIMOS spectra. To account for dust extinction, we measure the UV continuum slope(β)⁹ of the spectra via a linear fit to magnitude as a function of logarithmic wavelength over 1250–1650 Å in rest frame. It is not possible to estimate the UV slope using photometric data for three of our recently observed objects due to the absence of one or more of the key photometric bands. However, for MS 1358 we can compare our spectroscopically measured slope with that derived from available HST/ACS F850LP and HST/WFC3 F110W data. We find reasonable agreement given the uncertainties arising from the faint and extended nature of this object. Our UV slopes range from β ∼ −1.4 to −2.8, in good agreement with those determined for large samples of LBGs at similar redshifts (e.g., Bouwens et al. 2009; Castellano et al. 2012). The Calzetti et al. (2000) extinction law is applied to translate the UV slope to the extinction at 1600 Å. The SFR is then calculated from the extinction-corrected 1600 Å flux according to the relation given by Kennicutt (1998), in which a Salpeter (1955) initial mass function was assumed. Finally, we used the magnification from the lens models described in Section 3.3 to derive the intrinsic SFR. For the three lensed galaxies from Jones et al. (2013), we adopt a similar procedure, except the UV slopes were determined from the appropriate rest-frame UV HST photometry. For Abell 2390_H3a, the UV slope is calculated from HST/WFPC2 F814W and Keck/NIRC H-band photometry reported in Bunker et al. (2000). For Abell 2390_H5b and J1621+0647, the UV slopes are calculated from HST/WFPC2 F814W, HST/WFC3 F125W, and F110W images.

For the arc in the M1358+62 cluster, we can compare our derived SFR to that estimated from other methods in the literature. For the northwest arc (Arc C in Figure 1), we measure SFR_UV = 594 ± 202 M_⊙ yr⁻¹ (uncorrected for lensing). The result is consistent with that derived from an [O ii] emission line flux with no reddening correction at SFR_O[ii] = 525 ± 55 M_⊙ yr⁻¹ (Swinbank et al. 2009). Similarily, the unlensed SFRs for Abell 2390_H3a and Abell 2390_H5b of 95 ± 52 M_⊙ yr⁻¹ and 13 ± 5 M_⊙ yr⁻¹ are consistent with estimates reported in Pelló et al. (1999) using SED model fittings of ∼60 and ∼10 M_⊙ yr⁻¹, respectively.

Gravitational lensing enables us to directly measure covering fractions for galaxies at high redshift without the need for lens models since the magnification involved has no wavelength dependence. However, to derive the intrinsic properties such as the SFR, an accurate mass model for the lensing cluster is required. For each galaxy cluster we constructed a mass model using the light-traces-mass (LTM) method (full details in Zitrin et al. 2009, 2015). The model consists of a galaxy component and a smooth dark matter (DM) component. The galaxy component is a superposition of all cluster member galaxies, where the model assumes that the mass of each galaxy scales with its luminosity and can be described via an elliptical-shaped power-law profile. A core is often also introduced, especially for the brightest central galaxies (BCGs). The mass-to-light ratio is fixed for all cluster members except the BCG. The DM component is a smoothed version of the galaxy component, and the two components are added with a relative weight optimized in the minimization procedure. In addition, a two-component external shear is added to allow for further flexibility. We use a Markov Chain Monte Carlo to construct the mass model by minimizing differences between the predicted and actual position of multiple images. The multiple images and their spectroscopic redshifts were obtained from the literature, i.e., Swinbank et al. (2009), Zitrin et al. (2011), Smith et al. (2005), Gladders et al. (2002), and Rzepecki et al. (2007).

Using the mass model, we transfer the rest-frame UV image into the source plane and determine the magnification μ from the associated flux ratio. Uncertainties arise from both the precision of the model, which can be estimated from the predicted versus observed position of multiple images, and systematic effects associated with the LTM model compared, for example, with a model based on analytic profiles (see Zitrin et al. 2015). We adopt an uncertainty of 20% when the magnification μ < 10 and 50% if μ > 10. The magnification μ = 12.1 ± 6.1 derived for MS 1358 is consistent with the value of μ ≃ 12.5 derived by Swinbank et al. (2009), offering a valuable confirmation of our methods.

The goal of this paper is to interpret the covering fraction of low-ionization gas in terms of new constraints on the escape fraction of ionizing photons in an enlarged sample of 4 < z < 5 galaxies. We will introduce two methods for making the connection between the covering and escape fraction. However, first it is important to introduce some definitions.

The most important measure of the escape fraction is that used in calculating whether galaxies are capable of driving cosmic reionization. This is the absolute escape fraction, defined as

$$\begin{eqnarray*}&&{f}_{\mathrm{esc},\mathrm{abs}}\equiv \displaystyle \frac{{F}_{\mathrm{LyC},\mathrm{out}}}{{F}_{\mathrm{LyC},\mathrm{int}}},\end{eqnarray*}$$

that is, the ratio between the LyC flux density escaping the galaxy, usually measured at 900 Å, and the internal stellar flux density at the same wavelength. Most of the escape fractions reported in numerical simulations and theoretical calculations refer to this absolute value. However, observationally it is very difficult to estimate the internal stellar flux density at 900 Å since it requires modeling the SED beyond the range of direct observation. Instead, a more accessible estimate is the relative escape fraction, defined as

$$\begin{eqnarray*}&&{f}_{\mathrm{esc},\mathrm{rel}}\equiv \displaystyle \frac{{({f}_{\mathrm{LyC}}/{f}_{1500})}_{\mathrm{out}}}{{({f}_{\mathrm{LyC}}/{f}_{1500})}_{\mathrm{int}}}=\displaystyle \frac{{f}_{\mathrm{esc},\mathrm{abs}}}{{10}^{-0.4{A}_{1500}}}.\end{eqnarray*}$$

Since f_LyC,out can be additionally absorbed by neutral hydrogen in the IGM along the line of sight, we can write ${f}_{\mathrm{esc},\mathrm{rel}}$ as a function of the observed LyC flux(${f}_{\mathrm{LyC},\mathrm{obs}}$) as

$$\begin{eqnarray*}&&{f}_{\mathrm{esc},\mathrm{rel}}=\displaystyle \frac{{({f}_{\mathrm{LyC}}/{f}_{1500})}_{\mathrm{obs}}}{{({f}_{\mathrm{LyC}}/{f}_{1500})}_{\mathrm{int}}}{e}^{{\tau }_{{}_{{\rm{H}}{\rm{I}},\mathrm{IGM}}}}.\end{eqnarray*}$$

Notwithstanding the fact that the relative escape fraction can be obtained from good observations, modeling of the SEDs is still required to determine the intrinsic ratio of the 900 and 1500 Å luminosities, which results in additional uncertainties.

We can estimate the covering fraction of neutral hydrogen in a galaxy by measuring the amount of nonionizing UV radiation absorbed by metals following the approach described in Heckman et al. (2011) and Jones et al. (2013). Realistically, we must assume that galaxies are covered partially by H i regions with different velocities relative to the systemic value. Neutral clouds of a particular velocity will have a covering fraction of ${f}_{\mathrm{cov}}$(v). Photons of each low-ionization absorption line, when shifted to match the cloud velocity, will be able to escape if they travel through holes in the neutral gas or will otherwise be attenuated. The underlying principle is that low-ionization gas represents a faithful tracer of neutral hydrogen, an assumption discussed in detail by Jones et al. (2013), Henry et al. (2015), and Vasei et al. (2016). Trainor et al. (2015) find that low-ionization absorption-line widths correlate with the inferred escape of Lyα photons, and Rivera-Thorsen et al. (2015) find a strong correlation with 21 cm line widths that directly trace the neutral hydrogen. However, Reddy et al. (2016) claim, from their recent z ≃ 3 sample, that the outflowing gas may be metal-poor, in which case the true H i covering fraction could be underestimated somewhat by our technique.

The relevant radiative transfer equation is

$$\begin{eqnarray}&&\displaystyle \frac{I(v)}{{I}_{0}}=(1-{f}_{\mathrm{cov}}(v))+{f}_{\mathrm{cov}}(v){e}^{-\tau (v)},\end{eqnarray} \tag{ 1 }$$

where I₀ is the continuum level and τ is the optical depth of each absorption line,

$$\begin{eqnarray}&&\tau (v)=\displaystyle \frac{\pi {e}^{2}}{{m}_{e}c}N(v){f}_{{lu}}{\lambda }_{{lu}}=\displaystyle \frac{N(v)}{3.768\times {10}^{14}}{f}_{{lu}}{\lambda }_{{lu}},\end{eqnarray} \tag{ 2 }$$

where f_lu is an oscillator strength of each transition from level l to u (available from the NIST Atomic Spectra Databased), λ is the transition wavelength in Å, and N(v) is the ion column density in cm⁻² (km s⁻¹)⁻¹. Equation (2) assumes that none of the excited states are populated, i.e., N_u/N_l ≈ 0. Combining Equations (1) and (2), the attenuation in the continuum-normalized spectra at each velocity of an absorption line is then a function of two variables: ${f}_{\mathrm{cov}}$(v) and column density of absorbers N(v) with velocity v. If we have more than two absorption lines originating from the same ion, we can solve for the two variables.

This picture of a "picket fence" for the interstellar medium (ISM), whereby each individual absorbing cloud forms a "picket" for each velocity bin, is consistent with that determined to be an appropriate model for local analogs of high-redshift star-forming galaxies. A good example follows observations of the Ca ii H and K absorption lines in clouds around hot OB stars that resolve into several radial velocity components (Adams 1949). Moreover, Heckman et al. (2011) found that all four local Lyman break analogs with high residual intensities in the low-ionization absorption lines are a better match to the picket fence model compared to a uniform shell model.

Our DEIMOS spectra cover three absorption-line transitions of Si ii at 1260, 1304, and 1526 Å, enabling us to solve for both ${f}_{\mathrm{cov}}$(v) and N(v). We bin the absorption lines into several velocity components, with the width of each velocity bin chosen to be equal to the spectral resolution at that wavelength. We limit the usage of Si ii λ1304 to only those systems whose velocity v ≳ −200 km s⁻¹ in order to avoid contamination by O i λ1302. We use a brute-force (grid search) technique to find the best-fit parameters for each velocity bin. The likelihood of each pair of parameters is calculated from the least-squares residual ${\chi }^{2}=\sum {({I}_{\mathrm{obs}}-{I}_{N,{f}_{c}})}^{2}$. We adopt a prior range for ${f}_{\mathrm{cov}}$(v) of [0, 1], while the prior for N(v) is adjusted so that the range at each velocity covers the whole posterior probability distribution. An exception is made when the lines are optically thick, in which case only a lower limit on the column density can be obtained. Examples of posterior distributions for both optically thin and optically thick regimes are shown in Figure 3.

Zoom In Zoom Out Reset image size

In our second method, we simplify the covering fraction model given in Section 4.2 by assuming that all clouds are optically thick τ ≫ 1. Equation (1) then becomes

$$\begin{eqnarray}&&{f}_{\mathrm{cov}}(v)=1-\displaystyle \frac{I(v)}{{I}_{0}}.\end{eqnarray} \tag{ 3 }$$

This enables us to include absorption lines from other species present in our spectra, including O i λ1302 and C iiλ1334. We will see that the optically thick case is a reasonable assumption, as results from method 1 show that most of the Si ii absorption lines are optically thick.¹⁰ At each cloud velocity, we then calculate the covering fraction as an inverse-variance weighted mean of the line residuals.

The covering fractions estimated from low-ionization absorption profiles at each velocity resolution element are shown in Figure 4. In general, we find that the absorbing gas is spatially inhomogeneous: the absorption-line profiles show only partial covering of the stellar radiation with a maximum absorption depth of ∼50%–100%. Formally this translates to geometric escape fractions of up to 50%. Only two of the seven lensed galaxies, A2219 and J1621, are consistent with a complete (100%) covering fraction at any velocity and hence a zero escape fraction. Without accounting for possible spatial variations in the covering fraction as a function of velocity, which will be examined in Section 4.5, we find that a median value of covering fraction in our enlarged sample is f_cov,obs = 66%. If this represents the true H i covering fraction, it corresponds to ${f}_{\mathrm{esc},\mathrm{rel}}=34 \%$. However, to translate these measures to the desired absolute escape fractions, it is necessary to account for dust extinction since these galaxies are not dust-free, as indicated by the fact that their UV slopes β > −2 (see Table 2). Accounting for reddening based on the dust-in-cloud model (see Vasei et al. 2016), the corresponding median ${f}_{\mathrm{esc},\mathrm{abs}}$ and weighted standard deviation of our sample become 19% ± 6%.

Zoom In Zoom Out Reset image size

At face value such an absolute escape fraction, if typical of sources at higher redshift, would imply that early star-forming galaxies are readily capable of maintaining reionization (e.g., Robertson et al. 2015). Our median value derived from low-ionization absorption profiles is consistent with though somewhat larger than that inferred via recombination rate studies of the IGM using QSO absorption lines at a similar redshift, which imply f_esc ∼ 12% (Becker et al. 2015). If f_esc is indeed 12%, then we would expect an average ${f}_{\mathrm{cov},\mathrm{obs}}\sim 78 \%$ for our sample, accounting for reddening. It is likely that our estimate is higher than the true value due to several violations of the assumed simple "picket fence" model of absorption, as discussed in Jones et al. (2013) and Vasei et al. (2016). We now explore these possibilities quantitatively in order to estimate the degree by which our ${f}_{\mathrm{cov},\mathrm{obs}}$ measures may overestimate the true escape fraction.

First, it is possible that the nonionizing stellar radiation could be associated with a lower covering fraction of gas than that for the ionizing starlight. This would arise if, for example, the H i and low-energy ions preferentially cover younger stars. One way to estimate the potential importance of this effect is to estimate the fraction of detected UV continuum radiation (at ∼1500 Å), which arises from starlight whose photons are incapable of ionizing hydrogen. For a galaxy with a history of continuous star formation over 100 Myr, only around 10% of the UV radiation arises from spectral types later than B (Parravano et al. 2003), so this cannot be a significant effect. Furthermore, the Si ii transitions show no evidence of a higher covering fraction at shorter wavelengths. Such an effect would result in weaker absorption of Si ii λ1526 compared to λ1304 in an optically thick gas, which is not observed in the majority of our combined sample where both lines are free from sky contamination. We therefore conclude that a differential covering fraction of ionizing stars introduces <10% systematic error on our derived values.

A second uncertainty arises from the fact that scattered photons are capable of "filling in" the absorption profiles, mimicking a reduced covering fraction. In the extreme case where line filling is induced from optically thin gas, the Si ii λλ1304 and 1526 absorption profiles would be filled by 50% of the flux seen in the Si ii* λλ1309 and 1533 lines, respectively. The observed fine-structure emission-line fluxes therefore allow us to directly estimate the magnitude of this effect. In Abell 2219, where these emission lines are strongest, we find that scattering could only reduce the inferred covering fractions by at most 13%. On average across the full sample, we find that line filling can change our derived ${f}_{\mathrm{cov},\mathrm{obs}}$ by at most 4% (i.e., to ${f}_{\mathrm{cov},\mathrm{obs}}=70 \%$), which is small compared to the quoted uncertainties. The actual effect of line filling must be even less than this limit, as most cases indicate optically thick absorption.

In summary, accounting for the possibilities of wavelength-dependent covering fraction and line filling discussed above, we can revise the original covering fraction to at most ${f}_{\mathrm{cov},\mathrm{obs}}=66+10+4=80 \%$, implying ${f}_{\mathrm{esc},\mathrm{abs}}\sim 11 \% .$ However, we emphasize that these corrections represent the maximum conservative limits allowed by our data, and there is no conclusive evidence that the effects discussed actually reduce f_esc at all. Therefore, the difference between our covering fraction measurement and the ${f}_{\mathrm{esc},\mathrm{abs}}\sim 12 \%$ implied by IGM studies, if real, can only be reconciled if these are the only important effects.

The most challenging and remaining uncertainty arises from the possibility that absorbing gas at different velocities may cover different spatial lines of sight (e.g., Rivera-Thorsen et al. 2015). To examine this possibility, we would need to examine the spatial variation of the absorption-line profiles. Normally this would currently be impractical in any source at z > 4, but, by good fortune, this is possible in our most distant lensed source MS 1358. In the next section, we show that it is likely that such spatial variations may be the largest cause of difference between our derived covering fractions and the true ionizing escape fraction.

In the specific case of the highly magnified target MS 1358 at z = 4.93, we can examine the covering fraction of low-ionization gas in discrete subcomponents of the galaxy. When the region being examined spectroscopically is physically small, the angular distribution of clouds of different velocities should be less of a concern. In addition, if a particular region has a high SFR, the possibility that its stellar continuum is affected by nonionizing radiation should also be reduced. By examining the spectroscopic variation among the star-forming clumps in MS 1358, we can therefore address several of the uncertainties in the inferred escape fraction from our methods discussed above.

Adopting the lensing magnification from our LTM model (Section 3.3) and physical parameters for MS 1358 calculated in Swinbank et al. (2009), its total stellar mass is ∼(7 ± 2) × 10⁸ M_⊙ and its integrated SFR is ∼75 ± 15 M_⊙ yr⁻¹. The source itself is only a few kpc² in size. Since the individual clumps are conveniently linearly distributed along the arc, we have extracted spectra of each of three regions for specific clumps identified by Swinbank et al. (2009). Due to the limited spatial resolution of our ground-based spectroscopy and the S/N of our data, we group three of the smaller clumps (C1-3) together in extracting the spectra. The resulting arrangement targeting three regions (C0, C1–3, and C4) is shown in Figure 5 alongside the [O ii] image derived from the Swinbank et al. (2009) integral field unit observations.

Zoom In Zoom Out Reset image size

Importantly, we find that there is a significant variation in the low-ionization absorption profiles across the source (Figure 6). The inferred covering fraction ranges from near complete (100%) in C0 to only ≃40% in C4. The variation within MS 1358 is comparable to that across the integrated analysis of our full sample shown in Figure 4. As expected, it suggests that the escape fraction is governed by small-scale physics associated with star-forming regions. The variation between individual clumps indicates that covering fractions derived from integrated spectra may not be representative of the sum of all regions, and we can examine the degree of this effect. The SFR-weighted total of individual clumps is f_cov ∼ 0.9, whereas we find f_cov ∼ 0.8 based on the integrated light (Figure 7). It is possible that even higher covering fractions would be found at finer spatial resolution, although we note that Borthakur et al. (2014) found ${f}_{\mathrm{esc},\mathrm{rel}}\approx 1-{f}_{\mathrm{cov}}$ in a compact source of similar size to the clumps studied here (roughly a few hundred parsecs). We conclude that spatial variation in the covering fraction therefore reduces the limit on f_esc by a factor of ≃2 in this source. To help determine the extent to which this source is representative of the larger sample, in the following sections we examine trends in our sample as a function of both local and integrated physical properties.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since both Lyα radiation and LyC radiation are scattered by clouds of neutral hydrogen, the low column density region from which LyC can escape should also provide an escape route for Lyα photons. A suite of radiative transfer simulations in clumpy ISM models by Dijkstra & Gronke (2016) predicts a strong correlation between the escape fractions of Lyα and LyC. Indeed, observations indicate an anticorrelation between the low-ionization absorption depth and Lyα emission line equivalent width (W_Lyα) (e.g., Shapley et al. 2003; Jones et al. 2012, 2013; Erb et al. 2014). If this reflects the correlation between the Lyα and LyC escape fractions, not only does it imply that low-ionization absorption depths are a valuable proxy for the LyC escape fraction, but it also provides observational support for the clumpy ISM models.

We plot the relation between W_Lyα and the low-ionization absorption depth in Figure 7, combining both the integrated measures and those for the clumps identified in MS 1358 (Section 4.5). Most of the new data points support the anticorrelation found in Jones et al. (2013), as shown by the dashed line. However, two clumps within MS 1358 (C0 and C1+C2+C3) appear as outliers, although they are consistent with the large scatter found in the Dijkstra & Gronke (2016) simulations. The origin of the large scatter in the simulations is the cloud-covering factor (the average number of clouds along each sightline around the region). There is a higher Poisson probability that a line of sight would be free from obscuration when the cloud-covering factor is low. If we interpret Figure 7 according to the simulation model, we can infer that the two clumps in MS 1358 have a higher number of clouds given their increased SFRs.

Reddy et al. (2016) find that galaxies in their z ≃ 3 sample with redder UV continua tend to have larger covering fractions of neutral hydrogen. They propose that this correlation, modeled according to contributions from both photoelectric absorption and dust extinction, could form the basis of a new method to infer the escape fraction of high-redshift galaxies, thereby bypassing the assumption of our method that the opacity of low-ionization metal lines represents a reasonable proxy for that of neutral hydrogen. Although a correlation between the line-of-sight extinction, as inferred from the color excess $E(B-V)$, and f_cov seems reasonable, for the integrated spectra for our full sample, we nonetheless found no strong dependence between our absolute covering fraction and $E(B-V)$. However, the three data points from the resolved data on MS 1358 suggest that regions with higher covering fractions within the galaxy have redder UV continua and also stronger SFRs, as shown in Table 2. We will explore the variation of absorption-line strength and SFRs in the next section.

Our spatially resolved observations of MS 1358 presented in Section 4.5 suggested that the implied escape fraction of ionizing photons depends sensitively on the small-scale structure of its star-forming regions. Here we consider the origin of these variations in terms of both the local and integrated SFRs.

The SFR may have two opposing effects on the escape fraction. First, the SFR is directly correlated with the column density of gas, including LyC-absorbing H i as well as molecular H₂ and dust (Bigiel et al. 2008, 2010; Krumholz et al. 2012) However, strong feedback and radiation pressure in high-SFR regions may clear the absorbing clouds, resulting in a higher escape fraction. Understanding the balance between these two effects in early galaxies is critical in resolving the question of whether galaxies are capable of maintaining an ionized IGM.

According to our data shown in Figure 8, we see tentative evidence of a modest anticorrelation between SFR and the inferred escape fraction. Formally this trend is nonzero at the 2.8σ level, with a correlation coefficient of −0.53 ± 0.19 for the galaxy-integrated measurements. The trend is apparently stronger (though statistically consistent) for the individual clumps in MS 1358 as compared to that seen for the integrated values across the galaxies in the sample. If real, this difference might be expected since the escape fraction inferred for an entire galaxy is necessarily a spectroscopic measure averaged over clumps with a range of both low and high SFRs. The observed anticorrelation can be interpreted in two independent ways. First, regions of low SFR are associated with lower column densities of both interstellar and outflowing material (including H i and dust). Therefore, the probability of holes permitting leaking LyC photons is higher. Alternatively, regions with lower SFR may represent those whose rates were higher in the past and are now in decline. The associated delay time increases the chance that feedback effects can clear any covering clouds. Although both explanations may be valid, we consider the latter more likely because only a small column density of H i gas is needed to extinguish any leaking LyC radiation. Regardless of the surface density of low-ionization gas, it seems that feedback will be the more effective factor in governing both a low SFR and a smaller H i covering fraction.

Zoom In Zoom Out Reset image size

Most radiative transfer hydrodynamical simulations of the resolved ISM support our conjecture that a time delay will increase the escape fraction. Star formation may become increasingly "burst-like" at high redshift (Stark et al. 2009), and the escape fraction will then fluctuate on short timescales, out of phase with the SFR. Kimm & Cen (2014) find a delay of ≃10 Myr for a covering cloud to be dispersed by feedback from supernovae, while Ma et al. (2015) found a slightly shorter delay of ≃3–10 Myr. A larger delay will increase the strength of the anticorrelation between SFR and escape fraction.

One cautionary remark is relevant in considering how to interpret Figure 8 given that the SFR values are inferred from gravitationally lensed galaxies using our lensing model (Section 3.3). While the intrinsic SFR can be reliably estimated by dividing the observed value by the appropriate magnification, the covering fraction is measured from the absorption-line depths in the image plane. Hence, there is a potential bias such that the covering fraction so measured will refer to that of the regions with the highest magnifications. It would only be possible to appropriately transfer the absorption-line spectra into the source plane with resolved spectroscopy using an integral field unit spectrograph. Such a bias would, however, only be significant if there is a significant magnification gradient along an arc, e.g., when it is close to a critical curve. In our sample, this issue is not important. Both A1689 and A2219 have low and uniform magnifications, and the individual MS 1358 clumps are sufficiently small to have uniform magnifications within each region. When MS 1358 is considered as an entirety, we excluded Arc B (Figure 1), which has a high variation in magnification. MACS 0940 is the only system where the variation in magnification is significant. However, its SFR is relatively low, so such a bias is unlikely to affect its location on the relation in Figure 8.

It is clear that, despite the significant uncertainties, our moderate covering fractions imply absolute escape fractions that are higher than the averages observed via direct measures of the LyC radiation at lower redshift (e.g., Mostardi et al. 2015). On the other hand, if star formation is increasingly burst-like at high redshift, where galaxies are physically more compact, it is quite likely from the arguments discussed above that there is an increase in ${f}_{\mathrm{esc}}$ with redshift. Theoretical studies present conflicting predictions on the question of evolution in ${f}_{\mathrm{esc}}$. Ferrara & Loeb (2013) present an analytical calculation that illustrates why the average escape fraction should decrease with time during the reionization era. With an assumed value for the efficiency of star formation, they estimate that ${f}_{\mathrm{esc}}$ decreased by a factor of 3 from z = 10 to z = 6. However, sample-limited radiative hydrodynamical simulations show no significant variation over the same redshift interval (Ma et al. 2015).

In compiling the data for Figure 9, we include observations at redshifts z < 4 from Jones et al. (2013) (not shown in Table 1), Reddy et al. (2016), and Patrício et al. (2016). Across this increased redshift range, the figure shows that we cannot yet detect any significant time dependence. This is because the scatter in covering fraction is dominant and apparently constant over time, in agreement with our deduction that the escape fraction depends largely on the small-scale structure of clouds, their physical association with star formation regions, and the short timescale of activity, and is consistent with the behavior of simulated galaxies in Ma et al. (2015).

Zoom In Zoom Out Reset image size

As shown in Figure 9, absorption-line measurements of comparable quality to those presented here are still scarce at z = 2–3. Recent studies of the "cosmic horseshoe" at z = 2.4 further highlight the challenge in comparing direct and indirect methods. The covering fraction was found to be ∼60% using methods equivalent to ours (Quider et al. 2009), yet imaging below the Lyman limit indicates ${f}_{\mathrm{esc},\mathrm{rel}}\lt 0.08$ (Vasei et al. 2016). This might be reconciled if f_cov overestimates f_esc by a factor of 5, or by low IGM transmission (∼20% probability; Vasei et al. 2016), or a combination of both effects. At present the cause is not clear, and it will be of great interest to examine low-ion covering fractions in more low-redshift sources with direct f_esc measurements. This is becoming feasible following recent LyC detections in nearby galaxies using HST/COS (e.g., Izotov et al. 2016; Leitherer et al. 2016). Reddy et al. (2016) present new measures of the covering fraction of both the low-ionization gas and neutral hydrogen for stacked spectra derived from their Keck survey at z ≃ 3, providing a valuable complementary approach to that discussed here. Although we have highlighted the possible difficulties of inferring a mean covering fraction from a stack of many spectra with diverse kinematic absorption-line profiles (such as the variety shown in Figure 4), Reddy et al. (2016) draw from a much larger sample and have questioned some of the assumptions inherent in the use of low-ionization gas alone. However, within the observational uncertainties at this stage, both approaches give reasonably consistent results for the inferred escape fraction at z > 3.

Despite significant observational and theoretical effort, direct evidence for a rising escape fraction at high redshift is still lacking due to the scarcity of LyC detections. Our data set suggests a median ${f}_{\mathrm{esc},\mathrm{abs}}=19 \% \pm 6 \%$ at z = 4–5 if we interpret the covering fraction as ${f}_{\mathrm{cov}}=1-{f}_{\mathrm{esc},\mathrm{rel}}$. However, spatial variation may lower this value to ${f}_{\mathrm{esc},\mathrm{abs}}\approx 10 \%$ if the resolved properties of MS 1358 are typical of our sample, and other effects discussed in Section 4.4 may reduce the value slightly further (to ∼7%). These estimates are certainly compatible with independent determinations from the IGM, as well as an increase over the limits at z ∼ 2. Furthermore, the demographic trends in Figures 7 and 8 indicate that we may expect escape fractions to increase with redshift, commensurate with increasing W_Lyα (e.g., Law et al. 2012). We speculate that a key difference may be that star formation is more burst-like (e.g., higher specific SFR; Schreiber et al. 2015) at higher z. Intense short-duration starburst episodes may be capable of clearing sightlines through the local ISM, resulting in increased ionizing escape during periods of lower activity, giving rise to the trend in Figure 8.