THE CONNECTION BETWEEN REDDENING, GAS COVERING FRACTION, AND THE ESCAPE OF IONIZING RADIATION AT HIGH REDSHIFT^∗

Models were fit to the composite spectra in order to deduce the reddening, average neutral gas column density, and covering fraction for galaxies in our sample. To accomplish this, we constructed an intrinsic template from a combination of the Rix et al. (2004) models with the 2010 version of Starburst99 far-UV spectral synthesis models (Leitherer et al. 1999, 2010). We adopted a stellar metallicity of $0.28{Z}_{\odot }$, based on current solar abundances (Asplund et al. 2009), as this metallicity was found to best reproduce the UV stellar wind lines, including C iii λ1176 and the N v λ1240 and C iv λ1549 P-Cygni profiles. Furthermore, a constant star formation history and an age of 100 Myr were assumed based on results from stellar population modeling of $z\sim 3$ LBGs presented in Reddy et al. (2012). Nebular continuum emission was added to the stellar template to produce the final "intrinsic template," which we denote by the symbol m.

To account for absorption by neutral hydrogen, we added H i Lyman series absorption to the intrinsic template assuming a Voigt profile with a Doppler parameter b = 125 km s⁻¹ (the wings of the absorption lines are insensitive to the particular choice of b) and varying column densities $N({\rm{H}}\,{\rm{I}})$.⁸ The ${\rm{H}}\,{\rm{I}}$ absorption was blueshifted by 300 km s⁻¹ to account for the observed positions of these lines in the composites. The observed blueshift indicates that most of the neutral hydrogen seen in absorption is outflowing at several hundred km s⁻¹. Models that include ${\rm{H}}\,{\rm{I}}$ absorption are indicated by the symbol $m({\rm{H}}\,{\rm{I}})$. Lastly, we included Lyman–Werner absorption with varying column densities and covering fractions of ${{\rm{H}}}_{2}$ as described in Paper I. This absorption, which is unresolved, decreases the level of the continuum at $\lambda \simeq 912\mbox{--}1120\,\mathring{\rm{A}}$. Models that include both ${\rm{H}}\,{\rm{I}}$ and ${{\rm{H}}}_{2}$ are referred to by the symbol $m({\rm{H}}\,{\rm{I}}+{{\rm{H}}}_{2})$.

We considered two different cases when fitting the stacked spectra. In the first case, we calculated the average reddening of the composite assuming a 100% covering fraction of dust. This is analogous to deriving the reddening of a galaxy assuming a starburst attenuation curve (e.g., Calzetti et al. 2000; Reddy et al. 2015). In the subsequent discussion, the "updated" versions of the Calzetti et al. (2000) and Reddy et al. (2015) curves refer to those that are updated for our new measurements of the far-UV shape of the dust curve as presented in Paper I. The intrinsic template was reddened by various amounts ($E(B-V)=0.00\mbox{--}0.60$) assuming the updated attenuation curves. The best-fit $E(B-V)$ was obtained through ${\chi }^{2}$ minimization by comparing the reddened templates—i.e.,

$$\begin{eqnarray}&&{m}_{{\rm{final}}}={10}^{-0.4E(B-V)k}\times m,\end{eqnarray} \tag{ 1 }$$

where k is the attenuation curve—to the composite spectrum in the "Composite $E(B-V)$" wavelength windows listed in Table 1. These windows were chosen to avoid prominent interstellar and stellar absorption/emission features. Note that the windows used to measure $E(B-V)$ are unaffected by ${\rm{H}}\,{\rm{I}}$ or ${{\rm{H}}}_{2}$ absorption, thus obviating the need to fit for this absorption when determining the reddening of the composite.

Table 1.  Wavelength Windows for Fitting

Type of Fitting	λ (Å)
Composite $E(B-V)$ a	1268–1293
	1311–1330
	1337–1385
	1406–1521
	1535–1543
	1556–1602
	1617–1633
$N({\rm{H}}\,{\rm{I}})$ b	967–974
	1021–1029
	1200–1203
	1222–1248
$N({{\rm{H}}}_{2})$ c	941–944
	956–960
	979–986
	992–1020

Notes.

^aWindows used to determine the best-fit $E(B-V)$ corresponding to a composite spectrum. ^bWindows used to determine the best-fit column density $N({\rm{H}}\,{\rm{I}})$ and covering fraction ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ of neutral hydrogen. ^cWindows used to determine the best-fit column density $N({{\rm{H}}}_{2})$ and covering fraction ${f}_{{\rm{cov}}}({{\rm{H}}}_{2})$ of molecular hydrogen.

Download table as:  ASCIITypeset image

In Section 4, we present evidence from the composite spectra that indicates a non-unity covering fraction of gas for galaxies in our sample. In this second case, we fit models to the composite spectra assuming non-unity covering fractions of dust, and neutral and molecular hydrogen. The models were computed as:

$$\begin{eqnarray}{m}_{{\rm{final}}} & = & {10}^{-0.4E{(B-V)}_{{\rm{los}}}k}\\ & & \times [({f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})-{f}_{{\rm{cov}}}({{\rm{H}}}_{2}))m({\rm{H}}\,{\rm{I}})\\ & & +{f}_{{\rm{cov}}}({{\rm{H}}}_{2})m({\rm{H}}\,{\rm{I}}+{{\rm{H}}}_{2})]\\ & & +(1-{f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}}))m.\end{eqnarray} \tag{ 2 }$$

Here, ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ is the covering fraction of dust and neutral gas, and ${f}_{{\rm{cov}}}({{\rm{H}}}_{2})$ is the covering fraction of ${{\rm{H}}}_{2}$. This equation assumes that the ${{\rm{H}}}_{2}$ is shielded by ${\rm{H}}\,{\rm{I}}$ and that only the portion of the spectrum that emerges through the ${\rm{H}}\,{\rm{I}}$ and ${{\rm{H}}}_{2}$ gas is reddened. As shown in Sections 4.3 and 5, this particular formalism is motivated by the fact that the ISM is expected to be porous and that the covering fraction of ${{\rm{H}}}_{2}$ is much smaller than that of ${\rm{H}}\,{\rm{I}}$ (see Paper I and below). The last term in Equation (2) corresponds to the portion of the light that exits the galaxy unabsorbed, which is simply the intrinsic template multiplied by $1-{f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$. Finally, k in this equation now refers to a line-of-sight extinction curve in lieu of an attenuation curve, and $E{(B-V)}_{{\rm{los}}}$ refers to the corresponding line-of-sight reddening. These adjustments assume that a foreground uniform dust screen describes the dust-covered portions of the galaxy.

The models were computed assuming a range of line-of-sight reddening, and covering fractions and column densities of ${\rm{H}}\,{\rm{I}}$ and ${{\rm{H}}}_{2}$. In all cases, the models were smoothed to the spectral resolution of the composites, and modified by the CGM and IGM opacity appropriate for the average redshift of each composite. The average opacity was determined by randomly drawing line-of-sight absorber column densities from the distribution measured in the Keck Baryonic Structure Survey (Rudie et al. 2013), and taking an average over many lines-of-sight.⁹ The best-fitting $E{(B-V)}_{{\rm{los}}}$, ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, and $N({\rm{H}}\,{\rm{I}})$ were found by comparing the models to the stacked spectra in the "Composite $E(B-V)$" and "$N({\rm{H}}\,{\rm{I}})$" windows listed in Table 1. As discussed in Paper I, the "$N({\rm{H}}\,{\rm{I}})$" windows were chosen to include the core and the wings of Lyβ and Lyγ, and the wings of Lyα, while excluding regions of metal line absorption (e.g., C iii λ977 and C ii λ1036). The region between Lyα and N v λ1238 is particularly sensitive to the gas column density and is included in the fitting windows. The higher order Lyman transitions (e.g., Lyδ, Ly, etc.) were not used in the fitting due to line overcrowding. Furthermore, as noted above, the fitting windows do not include any interstellar metal absorption lines. Similarly, the best-fit values of ${f}_{{\rm{cov}}}({{\rm{H}}}_{2})$ and $N({{\rm{H}}}_{2})$ were determined by comparing the models to the stacked spectra in the "$N({{\rm{H}}}_{2})$" windows listed in Table 1. As discussed in Paper I, the "$N({{\rm{H}}}_{2})$" windows were chosen to avoid regions of strong metal line and ${\rm{H}}\,{\rm{I}}$ absorption.

We determined the extent to which uncertainties in the redshifts and the resolution of the individual galaxy spectra lead to non-zero fluxes in the ${\rm{H}}\,{\rm{I}}$ line cores in the stacked spectrum. To accomplish this, we added ${\rm{H}}\,{\rm{I}}$ absorption of varying column densities $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]=20.0\mbox{--}22.0$ to the intrinsic template, assuming (as above) a Voigt profile with a Doppler parameter b = 125 km s⁻¹, and a 100% covering of gas, resulting in an "${\rm{H}}\,{\rm{I}}$-absorbed template." Thus, this ${\rm{H}}\,{\rm{I}}$-absorbed template has saturated black ${\rm{H}}\,{\rm{I}}$ line cores. The ${\rm{H}}\,{\rm{I}}$-absorbed template was then smoothed to have spectral resolutions varying from 1.5 Å (the average rest-frame resolution inferred for the LyC spectra) up to a factor of two times worse (3.0 Å). For each value of the spectral resolution, we generated 121 versions of the smoothed ${\rm{H}}\,{\rm{I}}$-absorbed templates—one for each of the 121 LyC spectra contributing to the composite at $\lambda \lt 1150$ Å—where each version was shifted in wavelength according to a normal distribution with $\sigma =125$ km s⁻¹. This velocity spread is based on the uncertainties in the systemic redshifts when they are derived using interstellar absorption lines (see Steidel et al. 2010), and the vast majority of the LyC spectra have systemic redshifts based on the latter. These 121 versions were then averaged together to produce a composite model. Based on these simulations, we would have had to systematically underestimate the average resolution of the LyC spectra by $\gtrsim 50 \%$ in order to match the ${\rm{H}}\,{\rm{I}}$ line depths for $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]\gt 20.0$, a highly unlikely situation given that this 50% systematic difference is much larger than the rms uncertainty ($\lt 0.1$ Å) in determining the spectral resolution from the widths of the night sky lines. For the actual value of the resolution of our spectra and our typical redshift uncertainty, the cores of the ${\rm{H}}\,{\rm{I}}$ lines remain black, implying that none of the observed residual flux density can be attributed to spectral resolution or redshift uncertainties. We conclude that systematic uncertainties in spectral resolution cannot account for the residual line flux density at the ${\rm{H}}\,{\rm{I}}$ line cores.

We ruled out foreground contamination—i.e., low-redshift objects that are unresolved from the LBGs in the ground-based photometry and spectroscopy—as a possible origin for the residual line flux observed at the ${\rm{H}}\,{\rm{I}}$ line cores based on multiple lines of reasoning. First, we examined the contamination statistics for the sample of LBGs with spatially resolved Hubble Space Telescope (HST) multi-filter (UVJH) and ground-based NB3420 imaging (to target LyC emission at $z\sim 3$) analyzed in Mostardi et al. (2013, 2015). Their sample contains 16 LBGs with HST imaging, four of which have NB3420 detections. Of these four, two were identified as having foreground contaminants, such that the offending low-redshift interlopers would not have been resolved separately from the background LBGs in the ground-based photometry. Of the 12 objects that were not detected in the NB3420 imaging, only one suffers from foreground contamination. Thus, the contamination rates for NB3420 detections and non-detections are 2/4 and 1/12, respectively. Extrapolating to the full LBG sample of Mostardi et al. (2013), where five of 49 objects were detected in the narrow-band imaging, implies that the fraction of LBGs with contaminated photometry is $2/4\times 5/49+1/12\times 44/49\approx 0.13$. Note that what we have just computed is the foreground contamination rate. However, the quantity of relevance in the present context is not the contamination rate, but the actual fraction of flux contributed by foreground contaminants.

To calculate this flux contamination, we examined the photometry for the 16 LBGs with HST imaging from Mostardi et al. (2015). We summed the UVJH flux densities for the 16 objects, both including and excluding foreground contaminants that would have been unresolved from the LBGs in the ground-based imaging. The ratio of the total flux in each band including foreground contamination to that excluding this contamination is 2.38, 1.02, 1.03, and 1.02, respectively, for the UVJH bands. Accounting for the higher foreground contamination rate among the 16 LBGs ($3/16\approx 0.19$) relative to the rate among the full LBG sample ($\approx 0.13$, see above), and adopting the reasonable assumption that the flux contamination per object is roughly constant, the aforementioned ratios become 1.92, 1.01, 1.02, and 1.01, respectively, for the UVJH bands. In other words, the flux contamination fraction for the rest-frame non-ionizing UV/optical flux is $\lt 2 \%$. The flux contamination fraction in the U-band is substantially larger, a point that we revisit in Section 7. Based on these statistics, the contaminating flux density across the wavelength range spanning the Lyβ, Lyγ, and Lyδ lines is a factor of $\approx 5$ lower than the residual flux densities measured at the aforementioned line centers.

Note that the previously discussed contamination rates were calculated for LBGs that were spectroscopically confirmed using relatively shallow spectroscopic data (Reddy et al. 2008; Steidel et al. 2011; Mostardi et al. 2013). Spectroscopic blends may not be as readily identifiable in these shallow spectra compared to the substantially deeper LyC spectra analyzed here. As indicated in Section 2, the $\simeq 7\mbox{--}10\,{\rm{hr}}$ spectroscopic integrations obtained with LRIS-B allowed us to remove a handful ($\approx 5 \%$) of LBGs that were confused with foreground interlopers. In fact, this contaminating rate is consistent with the expected number of chance superpositions of foreground objects, as discussed in C. C. Steidel et al. (2016, in preparation), implying that the present sample has been expunged of most contaminants. Thus, the contamination rates calculated above should be regarded as conservative upper limits.

There is additional evidence that suggests that the residual line flux at the ${\rm{H}}\,{\rm{I}}$ line centers is intrinsic to the LBGs. Namely, in the subsequent analysis, we show that this residual flux varies systematically with $E(B-V)$ such that redder galaxies exhibit lower levels of residual emission (Section 5), a trend that would not be expected if foreground interlopers were dominating the observed flux. For example, of the three LBGs identified with foreground contamination from the sample of Mostardi et al. (2015), two have interlopers that result in slightly redder colors compared to the corrected photometry, and one has an interloper that results in a color that is indistinguishable from the corrected color. Notwithstanding the small statistics, it follows that an increase in contamination fraction for galaxies with redder $E(B-V)$ should lead to a corresponding increase in residual flux, a trend that is opposite to the one actually observed (Section 5).

To summarize our arguments, we find that systematic uncertainties in the resolution of our spectra, combined with random uncertainties in systemic redshifts, cannot explain the level of residual flux seen in the stacked spectra. Furthermore, calculations of the flux contamination fraction for the sample of LBGs examined in Mostardi et al. (2013, 2015) imply that $\gt 80 \%$ of the residual flux must be intrinsic to the LBGs and not a result of foreground interlopers. If the statistics of how the colors of LBGs are altered due to foreground interlopers extend to larger samples, then significant contamination would be expected to lead to a higher residual flux for redder galaxies, which is contrary to the trend that we observe, as discussed in Section 5.

We conclude with some words of caution regarding the covering fractions derived here. First, while the composite may indicate ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})\lt 1$, low resolution spectroscopy may mask the presence of very narrow absorption components that are opaque to LyC radiation.¹⁰ However, if the presence of such narrow (and black) components were the rule rather than the exception for galaxies in our sample, we would expect little correlation between the covering fraction derived from low resolution spectra and the Lyα escape fraction, counter to our observations (Sections 5 and 6).

Second, the method for computing ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ assumes that the different velocity components of the ${\rm{H}}\,{\rm{I}}$ gas are spatially coincident when viewed in projection. If this is not the case—i.e., if optically thick gas that absorbs far from line center covers sightlines that are disjoint from those where the absorption at line center dominates—then ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ deduced based on the overall depth of the ${\rm{H}}\,{\rm{I}}$ absorption will underestimate the true ${\rm{H}}\,{\rm{I}}$ covering fraction. As we show in Section 7, the two effects just mentioned (i.e., unresolved narrow absorption and spatially disjoint gas at different velocities) must be limited in scope. Even so, in light of this discussion, one must consider that the covering fractions derived in our subsequent analysis may be lower limits to the true fractions.

To demonstrate the requirement of a non-unity covering fraction, we modeled the composite spectrum of all galaxies in our sample (Figure 1) using the same procedure as above in fitting $E(B-V)$ and ${\rm{H}}\,{\rm{I}}$ (we ignored ${{\rm{H}}}_{2}$ for this exercise), but where we now assumed a unity covering fraction of both dust and gas. The depth of the Lyβ line can be reproduced with $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]=15.9$, but such a low column density does not match the combined red damping wing of Lyα and P-Cygni N v absorption, nor the depths of the higher Lyman series lines (Figure 2). Alternatively, the assumption of $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]=21.0$ gas reproduces the combined red damping wing of Lyα and P-Cygni N v absorption, but implies optically thick gas and hence saturated ${\rm{H}}\,{\rm{I}}$ lines that should reach zero intensity at the line cores for a unity covering fraction.

Zoom In Zoom Out Reset image size

Thus, it is clear that neither low or high column density gas with complete covering can simultaneously reproduce the residual flux at the ${\rm{H}}\,{\rm{I}}$ line centers and the damping wings of Lyα. Rather, the roughly equal depths of the Lyman series lines combined with the damping wings of Lyα suggest optically thick gas with $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]\gtrsim 21.0$ and ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})\lt 1$. Consequently, the depths of the ${\rm{H}}\,{\rm{I}}$ lines and the red damping wing of Lyα are most sensitive to changes in the covering fraction of optically thick ${\rm{H}}\,{\rm{I}}$.

From a physical standpoint, a partial covering of dust and gas can result from supernovae chimneys or galactic fountains, whereby supernovae explosions heat the ISM and drive hot gas into the galactic halo (e.g., Putman et al. 2012 and references therein). Aside from direct observations of such hot halos in the Milky Way and other nearby galaxies (e.g., Lynds & Sandage 1963; Shapiro & Field 1976; Armus et al. 1995; Lehnert & Heckman 1995; Strickland et al. 2004), the higher SFR surface densities of—and the ubiquitous presence of outflows in—high-redshift galaxies may increase the porosity of the ISM. Furthermore, simulations of high-redshift galaxies show that the LyC escape fraction is essentially independent of dust obscuration, with the ionizing photons escaping the galaxies through channels that are essentially free of dust (Gnedin et al. 2008; Razoumov & Sommer-Larsen 2010; Ma et al. 2016), as depicted in Figure 3. In this illustration, the column of gas is sufficient to completely bound the ionized region by neutral gas, where the bulk of the neutral and ionized gas along the line of sight are entrained in an outflow (see also Figure 1 in Zackrisson et al. 2013).

Zoom In Zoom Out Reset image size

Using the procedure laid out in Section 3, we fit a spectral model to the composite of all galaxies (full composite), assuming a non-unity covering fraction of dust and gas, and an SMC extinction curve. The model was obtained by simultaneously varying $E{(B-V)}_{{\rm{los/SMC}}}$, ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, $N({\rm{H}}\,{\rm{I}})$, ${f}_{{\rm{cov}}}({{\rm{H}}}_{2})$, and $N({{\rm{H}}}_{2})$, until a best-fit was reached. The best-fit parameters are $E{(B-V)}_{{\rm{los/SMC}}}=0.096$, $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]=21.0$, ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})=0.96$, $\mathrm{log}[N({{\rm{H}}}_{2})/{{\rm{cm}}}^{-2}]=20.9$, and ${f}_{{\rm{cov}}}({{\rm{H}}}_{2})=0.03$. By perturbing the composite spectrum by its errors and refitting these realizations, we estimated the errors in $E{(B-V)}_{{\rm{los}}}$, ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, and $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]$ to be $\sigma (E{(B-V)}_{{\rm{los}}})\approx 0.01$, $\sigma ({f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}}))\approx 0.01$, and $\sigma (\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}])\approx 0.20$ dex, respectively. Nominally, these errors include the covariance between the parameters. However, as ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ is determined primarily by the depth of the ${\rm{H}}\,{\rm{I}}$ lines and $N({\rm{H}}\,{\rm{I}})$ is constrained primarily by the wings, there is negligible covariance between the covering fraction and column density of ${\rm{H}}\,{\rm{I}}$. Furthermore, as noted in Paper I, the ${{\rm{H}}}_{2}$ lines are unresolved and therefore we cannot accurately determine the errors in the ${{\rm{H}}}_{2}$ covering fraction and column density. In other words, for the models to match the continuum level of the stacked spectrum in the "$N({{\rm{H}}}_{2})$" windows, the ${{\rm{H}}}_{2}$ column density may be increased while the covering fraction is decreased, or vice versa. As noted above, ${f}_{{\rm{cov}}}({{\rm{H}}}_{2})$ may also suffer a systematic bias depending on which IGM opacity prescription we adopt.

As shown in Figure 2 and in Paper I, a Voigt profile provides a poor description for the shape of the ${\rm{H}}\,{\rm{I}}$ lines (i.e., the Voigt profile underestimates the width of the lines when simultaneously matching their depths). To obtain a better fit for these lines, we derived the shape of the ${\rm{H}}\,{\rm{I}}$ line profile from the composite spectrum itself. To recover the ${\rm{H}}\,{\rm{I}}$ line profile, the intrinsic spectrum was multiplied by $1-{f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, modified for the opacity of the IGM, and subtracted from both the composite and its model fit, with the resulting spectra shown in Figure 4. This step allows us to remove the contribution of unabsorbed continuum prior to modeling the ${\rm{H}}\,{\rm{I}}$ lines. We then used the model fit subtracted by the intrinsic spectrum (i.e., red curve in Figure 4) to define "continuum" points by which the Lyβ, Lyγ, and Lyδ line profiles were normalized. There is a single continuum point each for Lyβ and Lyγ and we treated the continuum as constant across these line profiles. For Lyδ, the continuum was estimated as a linear function passing through the two continuum points bracketing the line. The normalized line profiles, transformed into velocity space (where the minimum of each absorption trough is forced to zero velocity), and cast in terms of the optical depth, are shown in Figure 4.

Zoom In Zoom Out Reset image size

The optical depth profile shown in Figure 4 was assumed for all of the Lyman series lines, except Lyα. Because the damping wings are weaker and more difficult to measure for the higher Lyman series lines due to line crowding, we modeled the damping wings of Lyα using a range of ${\rm{H}}\,{\rm{I}}$ column densities with a Doppler parameter b = 125 km s⁻¹—as noted earlier, the wings of the absorption are insensitive to the particular choice of b.

In order to investigate trends between continuum color excess and covering fraction, we divided the galaxies in our sample into five bins of $E(B-V)$. Specifically, the intrinsic template was reddened with different values of $E(B-V)$ assuming the Reddy et al. (2015) attenuation curve and a 100% covering fraction of dust, shifted to the observed frame for the redshift of each galaxy, multiplied by the IGM opacity relevant for that redshift, and then convolved with the G and ${ \mathcal R }$ filters. The best-fit $E(B-V)$ was taken to be the value that resulted in a model-generated $G-{ \mathcal R }$ that best matched the observed $G-{ \mathcal R }$ color, where the latter was corrected for Lyα emission or absorption as measured from the spectra. Composite spectra for objects in each bin were constructed according to the procedures laid out in Paper I and Section 2.

In fitting the stacked spectra assuming a porous ISM, we considered both the SMC and Milky Way extinction curves in the modeling (i.e., for the line-of-sight reddening) to span the possible shapes of the extinction curve in both low and high metallicity environments. The fitting was accomplished by simultaneously varying $E{(B-V)}_{{\rm{los}}}$, $N({\rm{H}}\,{\rm{I}})$ (which affects the degree of damping of the Lyα line), ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ (which affects the depths of all the Lyman series lines), and ${f}_{{\rm{cov}}}({{\rm{H}}}_{2})$, until the differences between the composite spectra and the models were minimized in the wavelength regions indicated by the "Composite $E(B-V)$," "$N({\rm{H}}\,{\rm{I}})$," and "$N({{\rm{H}}}_{2})$" entries in Table 1. As the ${{\rm{H}}}_{2}$ lines are unresolved, we fixed $\mathrm{log}[N({{\rm{H}}}_{2})/{{\rm{cm}}}^{-2}]=20.9$ in the fitting. We verified that the covering fractions of neutral gas obtained with and without including ${{\rm{H}}}_{2}$ in the model fits were similar within $\approx 1 \%$. The covering fractions of ${{\rm{H}}}_{2}$ are significantly smaller, with a typical value of ${f}_{{\rm{cov}}}({{\rm{H}}}_{2})\lesssim 0.05$. Additionally, each composite spectrum was fit assuming a 100% covering fraction of dust in order to derive its $E(B-V)$. The $E(B-V)$ computed from the composite spectra agree well with the average $E(B-V)$ of objects contributing to those spectra. For simplicity, in the following discussion the color excess derived assuming a unity covering fraction of dust is referred to as the "continuum reddening," $E(B-V)$, while the reddening deduced assuming a partial covering fraction of dust is referred to as the "line-of-sight reddening," $E{(B-V)}_{{\rm{los}}}$. The continuum reddenings derived assuming the updated attenuation curves of Reddy et al. (2015) and Calzetti et al. (2000) are denoted by "$E{(B-V)}_{{\rm{MOS}}}$" and "$E{(B-V)}_{{\rm{Calz}}}$," respectively.

The fits assuming the SMC extinction curve and a non-unity covering fraction of dust and neutral gas are shown in Figure 5, and the best-fit parameters and their errors are summarized in Table 2. The errors were calculated as the 1σ dispersion in the values obtained by perturbing the composite spectra according to their error spectra and re-fitting these realizations. Not surprisingly, the new ${\rm{H}}\,{\rm{I}}$ optical depth profile adopted for Lyβ and the higher Lyman series transitions (Section 5.2) provides a significantly better fit for these lines—as well as the continuum level at λ < 950 Å—than a single Voigt profile, though the inferred covering fractions remain identical (see Figure 1 for the model fit to the full composite).

Zoom In Zoom Out Reset image size

Table 2.  Best-fit $E(B-V)$ and ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$

$E{(B-V)}_{{\rm{MOS}}}$ a	$E{(B-V)}_{{\rm{Calz}}}$ a	$E{(B-V)}_{{\rm{los/SMC}}}$ b	${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ b	$\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]$ b	$E{(B-V)}_{{\rm{los/MW}}}$ c	${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ c	$\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]$ c
0.098 ± 0.005	0.104 ± 0.006	0.052 ± 0.002	0.92	20.3	0.100 ± 0.002	0.93	20.4
0.179 ± 0.005	0.184 ± 0.006	0.069 ± 0.002	0.95	20.3	0.138 ± 0.002	0.95	20.5
0.197 ± 0.005	0.203 ± 0.008	0.091 ± 0.001	0.96	20.5	0.174 ± 0.002	0.96	20.8
0.276 ± 0.006	0.285 ± 0.005	0.101 ± 0.001	0.97	21.0	0.187 ± 0.004	0.97	21.2
0.292 ± 0.006	0.301 ± 0.007	0.112 ± 0.001	0.97	21.0	0.245 ± 0.004	0.97	21.4

Notes.

^aContinuum reddening assuming the updated Reddy et al. (2015) and Calzetti et al. (2000) attenuation curves. ^bBest-fitting $E{(B-V)}_{{\rm{los}}}$, ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, and $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]$, assuming a non-unity covering fraction of neutral gas and dust and the SMC extinction curve. The uncertainties in the ${\rm{H}}\,{\rm{I}}$ covering fractions are $\sigma ({f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}}))=0.01$. The uncertainties in $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]$ are $\sigma (\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}])=0.2$ dex. ^cBest-fitting $E{(B-V)}_{{\rm{los}}}$, ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, and $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]$, assuming a non-unity covering fraction of neutral gas and dust and the Milky Way extinction curve. The uncertainties in the ${\rm{H}}\,{\rm{I}}$ covering fractions are $\sigma ({f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}}))=0.01$. The uncertainties in $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]$ are $\sigma (\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}])=0.2$ dex.

Download table as:  ASCIITypeset image

Several notable results from the fitting are summarized in Figure 6 and Table 3. Significant correlations are found between the continuum reddening and the line-of-sight reddening, the column density of ${\rm{H}}\,{\rm{I}}$, and the covering fraction of ${\rm{H}}\,{\rm{I}}$. For simplicity in the subsequent analysis, and unless indicated otherwise, we assume the updated Reddy et al. (2015) curve for the continuum reddening and the SMC extinction curve for deriving the line-of-sight extinction and the neutral gas covering fraction and column density. Other combinations of the extinction/attenuation curves yield similar results to those obtained with the default choices specified above.

Zoom In Zoom Out Reset image size

Table 3.  Summary of Correlation Analysis and Fits^a

x variable	y variable	a	b	ρb	pc
$E(B-V)$	${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$	−5.214 ± 0.736d	0.312 ± 0.074d	0.90	0.04
$\mathrm{log}E(B-V)$	$\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]$	1.60 ± 0.53e	21.76 ± 0.39e	0.85	0.07
$\mathrm{log}E(B-V)$	$\mathrm{log}E{(B-V)}_{{\rm{los}}}$	0.614 ± 0.025f	−0.632 ± 0.015f	0.96	0.01

Notes.

^aParameters are presented for the correlations obtained when assuming the updated Reddy et al. (2015) stellar attenuation curve and the SMC extinction curve. Other combinations of the attenuation/extinction curves result in correlations with similar fits, with similar Pearson correlation coefficients and similar probabilities for null correlations between the considered quantities. ^bPearson correlation coefficient. ^cProbability of a null correlation. ^dFree parameters a and b, as in Equation (15). ^eSlope and intercept, denoted by a and b, respectively, in the linear fit to $\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]$ versus $\mathrm{log}E(B-V)$. ^fSlope and intercept, denoted by a and b, respectively, in the linear fit to $\mathrm{log}E{(B-V)}_{{\rm{los}}}$ versus $\mathrm{log}E(B-V)$.

Download table as:  ASCIITypeset image

To compute the significance of the correlations between $E(B-V)$ and $E{(B-V)}_{{\rm{los}}}$, $N({\rm{H}}\,{\rm{I}})$, and ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, we calculated the Pearson correlation coefficient and probability of a null correlation between the various quantities while perturbing the measured values many times according to their measurement uncertainties. The values reported in Table 3 are the median Pearson correlation coefficients and median probabilities of a null correlation between $E(B-V)$ and ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, $N({\rm{H}}\,{\rm{I}})$, and $E{(B-V)}_{{\rm{los}}}$. The inferred ${\rm{H}}\,{\rm{I}}$ covering fractions assuming a partial covering fraction of both dust and gas are $\gtrsim 92 \%$, substantially larger than those deduced assuming a unity covering fraction of dust (${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})\simeq 0.70\mbox{--}0.80;$ see Paper I). In the latter case, as the entire spectrum is reddened including the uncovered portion, the line depths can be reproduced with smaller covering fractions of ${\rm{H}}\,{\rm{I}}$. With a non-unity covering fraction of dust, larger ${\rm{H}}\,{\rm{I}}$ covering fractions are required to reproduce the depths when unreddened continuum is contributed from the uncovered portion of the galaxy spectrum. At any rate, our results imply that the reddening of the UV continuum slope reflects an increasing covering fraction of ${\rm{H}}\,{\rm{I}}$ and dust, where the covering gas is also characterized by a higher line-of-sight reddening.

Furthermore, the column density of ${\rm{H}}\,{\rm{I}}$ increases with continuum reddening. The ratio of the neutral gas column density and the line-of-sight reddening is only marginally correlated with continuum reddening with a Pearson correlation coefficient of $\rho =0.56$ and a probability of null correlation of p = 0.28 (Figure 7). The mean ratio of the neutral gas column density and the line-of-sight reddening for our sample is $\langle N({\rm{H}}\,{\rm{I}})/E{(B-V)}_{{\rm{los}}}\rangle =(5.8\pm 3.3)\times {10}^{21}$ and $(5.4\pm 3.7)\times {10}^{21}$ cm⁻² mag⁻¹ for the SMC and Milky Way extinction curves, respectively, where the errors represent the standard deviation in the five values contributing to the mean. For context, the average values found for SMC and MW sightlines are $\sim 2.1\times {10}^{22}$ and $\sim 4.5\times {10}^{21}$ cm⁻² mag⁻¹, respectively (Welty et al. 2012). A note of caution is that the $\langle N{({\rm{H}}{\rm{I}})/E(B-V)}_{{\rm{los}}}\rangle$ values calculated here are analogous to, but not the same as, the gas-to-dust ratio. Specifically, we do not make any accounting of either the molecular or ionized gas in our column density estimates, and much of the line-of-sight reddening can be attributed to dust in the outflowing ionized gas in high-redshift galaxies (see Section 6.2).

Zoom In Zoom Out Reset image size

Given the correspondence between Lyα strength and reddening/covering fraction noted in previous studies (e.g., Shapley et al. 2003; Kornei et al. 2010; Jones et al. 2013; Wofford et al. 2013; Rivera-Thorsen et al. 2015; Trainor et al. 2015), we calculated the absolute Lyα escape fraction,

$$\begin{eqnarray}&&f{({\rm{Ly}}\alpha )}_{{\rm{abs}}}=\displaystyle \frac{L{({\rm{Ly}}\alpha )}_{{\rm{obs}}}}{L{({\rm{Ly}}\alpha )}_{{\rm{int}}}},\end{eqnarray} \tag{ 3 }$$

where $L{({\rm{Ly}}\alpha )}_{{\rm{obs}}}$ and $L{({\rm{Ly}}\alpha )}_{{\rm{int}}}$ are the average observed and intrinsic Lyα line luminosities, respectively, for the galaxies contributing to each composite shown in Figure 5. The observed line luminosity for each spectral stack was calculated as follows. We first redshifted the composite to the mean redshift of objects contributing to that composite. The spectrum was then converted to ${f}_{\lambda }$ units. The observed Lyα line luminosity was calculated by directly integrating the flux over just the Lyα emission profile (i.e., from the points where the emission line intersects with the underlying Lyα absorption profile) using the IRAF task splot. Naturally, this calculation will not account for the potentially large fraction of Lyα (up to 80%; Steidel et al. 2011) that may be resonantly scattered into a diffuse halo (see discussion below).

The intrinsic Lyα line luminosity was estimated from the SFR. The SFR was calculated using the UV luminosity density at 1500 Å, corrected for dust based on the average $E(B-V)$ of objects contributing to each composite and assuming the Reddy et al. (2015) attenuation curve¹¹ , and converted to SFR assuming the Kennicutt (1998) relation for a Salpeter (1955) IMF. The intrinsic Lyα luminosity was then estimated from the SFR using a modified form of the relation given in Kennicutt (1998), where we assumed that the intrinsic ratio of the Lyα-to-Hα line intensities is $f({\rm{Ly}}\alpha )/f({\rm{H}}\alpha )\simeq 8.7$ for Case B recombination and ${T}_{{\rm{e}}}={10}^{4}$ K. This yields the following relationship between SFR and Lyα line luminosity, again assuming a Salpeter (1955) IMF:

$$\begin{eqnarray}&&L{({\rm{Ly}}\alpha )}_{{\rm{int}}}\,[{\rm{erg}}\,{{\rm{s}}}^{-1}]=\displaystyle \frac{{\rm{SFR}}\,[{M}_{\odot }\,{{\rm{yr}}}^{-1}]}{9.1\times {10}^{-43}}.\end{eqnarray} \tag{ 4 }$$

The measurement error in the mean escape fraction was estimated by perturbing the composite many times according to its error spectrum, and calculating the observed and intrinsic Lyα luminosities for each of these realizations. Note that the observed Lyα line luminosity and SFR are calculated using the spectrum itself; we made no corrections for slit loss, and thus the escape fractions derived here assume that the Lyα and UV continuum fluxes arise from spatially coincident regions. The resulting $f{({\rm{Ly}}\alpha )}_{{\rm{abs}}}$ are indicated for each composite in Figure 5.

Not surprisingly, $f{({\rm{Ly}}\alpha )}_{{\rm{abs}}}$ correlates with continuum reddening (and hence with ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}});$ Figure 8; see also Steidel et al. 2010) with a Pearson correlation coefficient of $\rho =-0.97$ and a probability of a null correlation of p = 0.005. For four of the five $E(B-V)$ bins, the absolute escape fraction of Lyα photons agrees well with our expectation if such photons only escaped through clear sightlines in the ISM; namely, $f{({\rm{Ly}}\alpha )}_{{\rm{abs}}}\approx 1-{f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ (Figure 8). The bluest bin of $E(B-V)$ exhibits a larger $f{({\rm{Ly}}\alpha )}_{{\rm{abs}}}$ than this expectation. This may suggest that a significant fraction of the Lyα flux may be resonantly scattered into the spectroscopic slit aperture and/or may arise from Lyα photons scattering out of resonance, though we note that the kinematic profile of the Lyα line for the galaxies in the bluest $E(B-V)$ bin is not substantially different than that of the other bins. Moreover, it is unclear why the fraction of Lyα photons scattered into the spectroscopic slit would be much larger for galaxies with bluer $E(B-V)$. A more detailed assessment of how Lyα photons escape these galaxies (as judged through the kinematic profile of the line) must await more precise systemic redshifts derived from nebular emission lines.

Zoom In Zoom Out Reset image size

In the context of Figure 3, one might expect the fraction of Lyα emission exiting through low column density sightlines to correlate more directly with ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ than the diffuse Lyα component; the latter is not accounted for in our calculation. However, based on the results of the previous section, both $N({\rm{H}}\,{\rm{I}})$ and $E{(B-V)}_{{\rm{los}}}$ increase in tandem with ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$. Hence, Lyα photons that diffuse far from the sites of star formation face an increasing probability of absorption by dust with increasing ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, implying that the fraction of Lyα scattered into a diffuse halo and ultimately exiting the ISM/CGM of the galaxy is also likely to correlate with ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$. This inference could have also been deduced from the results of Steidel et al. (2011), who find that—when considering resonantly scattered Lyα—the attenuation of Lyα photons is consistently larger than that of the UV continuum, and the fraction of the total Lyα produced from star formation that diffuses into a halo correlates inversely with the spectroscopic Lyα flux. At any rate, the computed $f{({\rm{Ly}}\alpha )}_{{\rm{abs}}}$ imply that the selection of galaxies based on strong Lyα emission will also tend to isolate those with lower gas covering and higher LyC escape fractions. In accordance with expectations, the conditions that give rise to strong Lyα emission also facilitate the escape of LyC radiation (e.g., Dijkstra et al. 2016).

Given the prevalent use of low-ionization interstellar absorption lines as proxies for the neutral gas covering fraction in high-redshift galaxies (e.g., Shapley et al. 2003; Jones et al. 2013; Alexandroff et al. 2015; Trainor et al. 2015), it is instructive to directly compare the covering fractions deduced from the low-ionization lines with those computed from the depth of the ${\rm{H}}\,{\rm{I}}$ lines. To accomplish this, the five composites shown in Figure 5 were divided by their respective best-fitting models to produce normalized composites. This method of normalization has the advantage of removing any stellar absorption components (e.g., C ii $\lambda 1334$) that may affect equivalent width measurements of the interstellar absorption lines. The normalized composites were transformed into velocity space around the Si ii $\lambda 1260$, Si ii $\lambda 1527$, C ii $\lambda 1334$, and Al ii $\lambda 1670$ IS lines. These include the strongest and most commonly used IS lines observed in the rest-frame UV spectra of high-redshift galaxies.

These features are shown in Figure 9, with the different colors indicating the normalized spectrum in each bin of $E(B-V)$. If the two Si ii transitions at 1260 and 1527 Å are optically thin, then the ratio of their equivalent widths will be ${W}_{{\rm{1260}}}/{W}_{{\rm{1527}}}\gtrsim 6$. The observed ratio is ${W}_{{\rm{1260}}}/{W}_{{\rm{1527}}}\simeq 1$, implying that the lines are saturated and hence their depths are sensitive to the covering fraction of Si ii-enriched material. Other strong low-ionization IS absorption features are shown in Figure 9, including C ii $\lambda 1334$, and Al ii $\lambda 1670$ lines. The similarity in the velocity widths and depths of these additional lines to those of the saturated Si ii transitions discussed above suggests that the former may be also arise from the same gas as the latter. If this is the case, then collectively the depths of the low-ionization lines indicate a covering fraction that rises from ${f}_{{\rm{cov}}}(\mathrm{ion})\approx 0.30$ to 0.65 proceeding from the bluest to reddest $E(B-V)$ bin (see also Shapley et al. 2003). The errors in ${f}_{{\rm{cov}}}(\mathrm{ion})$ were calculated by perturbing many times the normalized composite spectra by the normalized error spectra, and measuring the depth of the absorption lines for each of these realizations. The dispersion in the absorption line depths yielded the measurement uncertainty in ${f}_{{\rm{cov}}}(\mathrm{ion})$, but of course does not include systematic errors associated with the specific model for the continuum level.

Zoom In Zoom Out Reset image size

The covering fractions of the low-ionization species are correlated with, but are systematically lower than, those derived from the ${\rm{H}}\,{\rm{I}}$ lines (lower right panel of Figure 9). There are several possible explanations for why the low-ionization lines suggest lower covering fractions. First, if most of the metals inhabit narrow and unresolved absorption systems, then the depth of the metal lines will underpredict the covering fraction. Likewise, if the metal-enriched gas has a substantial velocity gradient across the face of the galaxy, then the derived covering fractions may be lower limits (e.g., see discussion in Section 4.3). While we cannot definitively rule out either of these possibilities, the similarity in the blueshift ($\approx 200$ km s⁻¹) and velocity spread (FWHM $\approx 800$ km s⁻¹; e.g., compare Figure 9 to the right panel of Figure 4) of the low-ionization and ${\rm{H}}\,{\rm{I}}$ lines suggests that the two components are kinematically coupled, where the former covers a smaller portion of the galaxy's continuum. Such a configuration may arise if dense and metal-enriched pockets of gas are embedded within the outflowing ${\rm{H}}\,{\rm{I}}$. Consequently, some substantial fraction of the outflowing gas may be metal-poor.

Additionally, we note that ${f}_{{\rm{cov}}}({\rm{ion}})$ increases by a factor of $\approx 1.5$ in response to a small ($\approx 5 \%$) increase in ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$. This behavior suggests that as ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ increases, a larger fraction of the neutral gas is enriched with Si ii, implying a corresponding increase in the average dust-to-gas ratio of the outflowing gas. Notwithstanding the fact that the line-of-sight reddening will be sensitive to the total column of dust (including that of the ionized gas), an increase in the dust-to-gas ratio of the outflowing gas for objects with redder continuum would result in a strong correlation between $E{(B-V)}_{{\rm{los}}}$ and $E(B-V)$, consistent with what we actually observe (left panel of Figure 6).

Figure 10 shows the profiles of the two Si iv $\lambda \lambda 1393,1402$ transitions for each bin of reddening. As the ionization potential for Si iv (q = 33.5 eV) is much higher than that required for H-ionization, the associated transitions arise from H-ionized gas. The lines are blueshifted by $\approx 160$ km s⁻¹, slightly smaller than the blueshifts derived from the low-ionization and ${\rm{H}}\,{\rm{I}}$ lines, though a more detailed comparison of the kinematics of the neutral, low-, and high-ionization gas will rely on future precise measurements of the systemic redshifts of galaxies in our sample. The ratio of the equivalent widths of the Si iv lines is ${W}_{{\rm{1393}}}/{W}_{{\rm{1402}}}\approx 2$, indicating that the absorption is consistent within the errors of being on the linear part of the curve of growth. As the depths of the lines do not correlate obviously with the reddening of the continuum, we surmise that the outflowing ionized gas has an optical depth $\tau \simeq 1$. Note that the errors in ${W}_{{\rm{1393}}}/{W}_{{\rm{1402}}}$ do not rule the possibility of the lines being saturated, in which case the covering fraction of the ionized material is ${f}_{{\rm{cov}}}$(Si iv) ≳ 30%. More generally, the presence of these high-ionization IS absorption lines suggests an ionized medium mixed with metals and dust.

Zoom In Zoom Out Reset image size

To summarize, we find that the absolute escape fraction of Lyα photons, $f{({\rm{Ly}}\alpha )}_{{\rm{abs}}}$, inversely correlates with ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, suggesting that the neutral gas covering fraction plays a significant role in modulating the escape of Lyα photons. The low-ionization interstellar doublet Si ii λλ1260, 1527 is saturated, exhibits similar absorption depths and velocity profiles to those of other low-ionization IS absorption lines, and suggests a covering fraction of low-ionization material, ${f}_{{\rm{cov}}}(\mathrm{ion})$, that increases with continuum reddening (and ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$) from $\approx 0.30$ for the bluest galaxies in our sample to $\approx 0.65$ for the reddest ones. However, ${f}_{{\rm{cov}}}(\mathrm{ion})$ is systematically lower than ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, while the low-ionization lines have similar blueshifts and velocity profiles as the ${\rm{H}}\,{\rm{I}}$ lines, suggesting that galaxy outflows consist in part of a neutral medium embedded with clumps of metal-enriched gas. Moreover, the rate of increase of ${f}_{{\rm{cov}}}(\mathrm{ion})$ relative to ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ implies that the average column of dust increases with covering fraction, consistent with our observation that the line-of-sight reddening correlates with continuum reddening. These results imply that a significant fraction of the neutral medium may be chemically pristine, and that low-ionization interstellar absorption lines may severely underpredict the neutral gas covering fraction. A similar analysis for the high-ionization absorption lines, namely Si iv λλ1393, 1402, suggests that the high-ionization gas is not optically thick and has ${f}_{{\rm{cov}}}$(Si iv) ≳ 30%.

LyC photons may be extinguished by dust absorption or through hydrogen photoelectric absorption. In the case of dust absorption, our new measurement of the far-UV shape of the dust curve at high redshift, presented in Paper I, implies a lower attenuation of LyC photons than those inferred from simple extrapolations of the Calzetti et al. (2000) and Reddy et al. (2015) attenuation curves. Alternatively, the high column densities ($\mathrm{log}[N({\rm{H}}\,{\rm{I}})/{{\rm{cm}}}^{-2}]\gtrsim 20.3$) inferred from the fitting presented in Section 5 suggest that the covered portions of a galaxy are optically thick to LyC radiation, with implied optical depths at 900 Å of ${\tau }_{900}\simeq 1150$. Thus, for the average galaxy, photoelectric absorption rather than dust dominates the depletion of LyC photons. Similar results have been suggested by simulations of escaping ionizing radiation from high-redshift galaxies (Gnedin et al. 2008; Razoumov & Sommer-Larsen 2010; Ma et al. 2016). Accordingly, in the geometry illustrated in Figure 3, LyC photons will only escape from dust-free holes in the ISM. The very high optical depth of the covering neutral gas suggests that Lyα photons will also predominantly escape from regions cleared of dust and gas, an expectation that is verified by the strong correlation between $f{({\rm{Ly}}\alpha )}_{{\rm{abs}}}$ and reddening/covering fraction (Section 6.1).

LyC extinction may have a more significant effect in galaxies where the radiation field is sufficient to completely ionize the hydrogen (i.e., corresponding to a "density-bounded" case; see Zackrisson et al. 2013), but such galaxies must be relatively rare in our sample given the apparent frequency of high column density neutral gas as reflected in the composite spectra. We can connect the gas covering fraction to the absolute escape fraction of LyC photons, where the latter is commonly written as:

$$\begin{eqnarray}f{({\rm{LyC}})}_{{\rm{abs}}} & = & \exp [-{\tau }_{{\rm{ISM}}}({\rm{LyC}})]\\ & & \times {10}^{-0.4k({\rm{LyC}})E(B-V)}.\end{eqnarray} \tag{ 5 }$$

The first term in this expression accounts for the opacity of the ISM, where ${\tau }_{{\rm{ISM}}}({\rm{LyC}})$ is the optical depth at LyC wavelengths. The second term accounts for the obscuration of LyC photons by dust, where $E(B-V)$ is the reddening and $k({\rm{LyC}})$ is the dust attenuation curve evaluated at LyC wavelengths. In the case where the dust and neutral gas completely cover the continuum, the opacity of the ISM can be combined with the dust attenuation curves calculated in Paper I to compute $f{({\rm{LyC}})}_{{\rm{abs}}}$.

In the case of incomplete covering, the absolute escape fraction of LyC photons can be written as:

$$\begin{eqnarray}f{({\rm{LyC}})}_{{\rm{abs}}} & = & {f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})\exp [-{\tau }_{{\rm{ISM}}}^{{\rm{cov}}}({\rm{LyC}})]\\ & & \times {10}^{-0.4k({\rm{LyC}})E(B-V{)}_{{\rm{los}}}^{{\rm{cov}}}}\\ & & +(1-{f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}}))\exp [-{\tau }_{{\rm{ISM}}}^{{\rm{uncov}}}({\rm{LyC}})]\\ & & \times {10}^{-0.4k({\rm{LyC}})E(B-V{)}_{{\rm{los}}}^{{\rm{uncov}}}},\end{eqnarray} \tag{ 6 }$$

where the first and second terms correspond to the covered and uncovered portions contributing to the galaxy spectrum, respectively. In this expression, $k({\rm{LyC}})$ is a line-of-sight extinction curve evaluated at LyC wavelengths, and $E(B-V{)}_{{\rm{los}}}^{{\rm{cov}}}$ and $E(B-V{)}_{{\rm{los}}}^{{\rm{uncov}}}$ indicate the line-of-sight reddening appropriate to the covered and uncovered continuum, respectively.

Similarly, the absolute escape fraction of UV continuum photons (e.g., typically measured between 1500 and 1700 Å) is

$$\begin{eqnarray}&&f{({\rm{UV}})}_{{\rm{abs}}}={10}^{-0.4k({\rm{UV}})E(B-V)},\end{eqnarray} \tag{ 7 }$$

where $k({\rm{UV}})$ refers to the attenuation curve relevant for the stellar continuum, assuming a unity covering fraction of dust.¹²

These equations are cast most conveniently in terms of the "relative" escape fraction defined as

$$\begin{eqnarray}&&{f}_{{\rm{rel}}}=\displaystyle \frac{f{({\rm{LyC}})}_{{\rm{abs}}}}{f{({\rm{UV}})}_{{\rm{abs}}}}\end{eqnarray} \tag{ 10 }$$

(Steidel et al. 2001). The relationship between the relative escape fraction and the observed ratio of LyC-to-UV flux density is

$$\begin{eqnarray}{f}_{{\rm{rel}}} & = & \displaystyle \frac{S{({\rm{LyC}})}_{{\rm{obs}}}/S{({\rm{LyC}})}_{{\rm{int}}}}{S{({\rm{UV}})}_{{\rm{obs}}}/S{({\rm{UV}})}_{{\rm{int}}}}\\ & & \times \exp [{\tau }_{{\rm{IGM}}}({\rm{LyC}})]\times \exp [{\tau }_{{\rm{CGM}}}({\rm{CGM}})],\end{eqnarray} \tag{ 11 }$$

where $S{({\rm{LyC}})}_{{\rm{obs}}}$ and $S{({\rm{UV}})}_{{\rm{obs}}}$ are the observed flux densities of LyC and UV continuum photons, respectively; $S{({\rm{LyC}})}_{{\rm{int}}}$ and $S{({\rm{UV}})}_{{\rm{int}}}$ are the intrinsic flux densities of LyC and UV continuum photons, respectively, as typically estimated from stellar population synthesis models; and the last two terms indicate the further depletion of LyC photons due to the opacity of the IGM and CGM.

Combining Equations (7), (10), and (11) yields the following expression for the observed LyC-to-UV flux density ratio:

$$\begin{eqnarray}{\left[\displaystyle \frac{S({\rm{LyC}})}{S({\rm{UV}})}\right]}_{{\rm{obs}}} & = & {\left[\displaystyle \frac{S({\rm{LyC}})}{S({\rm{UV}})}\right]}_{{\rm{int}}}\\ & & \times \exp [-{\tau }_{{\rm{IGM}}}({\rm{LyC}})]\\ & & \times \exp [-{\tau }_{{\rm{CGM}}}({\rm{LyC}})]\\ & & \times \left[\displaystyle \frac{f{({\rm{LyC}})}_{{\rm{abs}}}}{{10}^{-0.4k({\rm{UV}})E(B-V)}}\right],\end{eqnarray} \tag{ 12 }$$

where $f{({\rm{LyC}})}_{{\rm{abs}}}$ is given by Equation (6).

For the average $z\sim 3$ galaxy in our sample, the first term of Equation (6) is negligible owing to the large neutral gas column density and hence the large value of ${\tau }_{{\rm{ISM}}}^{{\rm{cov}}}$, as noted above. Hence, if the uncovered portions of the galaxy continuum are characterized by very low gas and dust column densities, as would be expected in the simple geometry of Figure 3, then the absolute escape fraction of LyC photons is simply

$$\begin{eqnarray}&&f{({\rm{LyC}})}_{{\rm{abs}}}\approx 1-{f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}}).\end{eqnarray} \tag{ 13 }$$

We noted in Section 4.3 that $f{({\rm{LyC}})}_{{\rm{abs}}}$ may be lower than that indicated by Equation (13). We discuss this further in Section 7.4.3. In the present context, Equation (12) becomes

$$\begin{eqnarray}{\left[\displaystyle \frac{S({\rm{LyC}})}{S({\rm{UV}})}\right]}_{{\rm{obs}}} & = & {\left[\displaystyle \frac{S({\rm{LyC}})}{S({\rm{UV}})}\right]}_{{\rm{int}}}\\ & & \times \exp [-{\tau }_{{\rm{IGM}}}({\rm{LyC}})]\\ & & \times \exp [-{\tau }_{{\rm{CGM}}}({\rm{LyC}})]\\ & & \times \left[\displaystyle \frac{1-{f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})}{{10}^{-0.4k({\rm{UV}})E(B-V)}}\right].\end{eqnarray} \tag{ 14 }$$

In practice, ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ is difficult to constrain in the absence of direct spectroscopy of the ${\rm{H}}\,{\rm{I}}$ absorption lines—e.g., as noted in Section 6.2, the neutral gas covering fractions deduced from low-ionization IS absorption lines may severely underpredict those derived from the ${\rm{H}}\,{\rm{I}}$ lines. Thus, the utility of an empirical relation between ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ and $E(B-V)$ becomes obvious: the continuum reddening is relatively easy to measure, even with just a few broadband photometric filters covering the rest-frame UV. As discussed in Section 5.4 and shown in Figure 6 and Table 3, we find a significant correlation between ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ and $E(B-V)$. Thus, a suitable functional relationship between these two variables can be used to cast Equations (6) and (14) in terms of $E(B-V)$. The functional form relating ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ and $E(B-V)$ is the subject of the next subsection.

The middle and right panels of Figure 6 show the relationship between ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ and $E(B-V)$. The Pearson correlation coefficients when we assume the SMC and Milky Way extinction curves in computing ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ are $\rho =0.90$ and 0.87, respectively, with probabilities for no correlation between ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ and $E(B-V)$ of p = 0.04 and 0.05, respectively. To uncover the functional form of ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ versus $E(B-V)$, we reddened the intrinsic template assuming different values of $E(B-V)$ and the updated and original forms of the Reddy et al. (2015) and Calzetti et al. (2000) attenuation curves. For each of these reddened templates, we determined which combination of the covering fraction and line-of-sight reddening—assuming either the SMC or Milky Way extinction curves—provided the best match in the "Composite $E(B-V)$" wavelength windows listed in Table 1, and one additional window covering the LyC region: $\lambda =800\mbox{--}900\,\mathring{\rm{A}}$. The model expectations for various assumptions of the reddening and line-of-sight dust curves are shown in the middle and right panels of Figure 6. For clarity, not all of the combinations of attenuation/extinction curves are shown, as they all result in a similar functional behavior between ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ and $E(B-V)$.

As expected, the covering fraction asymptotes to zero and unity for small and large values of the continuum reddening, respectively. Motivated by this behavior, we adopted the following functional form for ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ versus $E(B-V)$:

$$\begin{eqnarray}&&{f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})=1-\exp [a\times E{(B-V)}^{b}],\end{eqnarray} \tag{ 15 }$$

where a and b are free parameters. Table 3 summarizes the free parameters obtained by fitting the measured values of ${f}_{{\rm{cov}}}$ and $E(B-V)$ for the updated Reddy et al. (2015) attenuation curve and the SMC extinction curve applicable to the continuum and line-of-sight reddening, respectively.

To compute these fits, we randomly perturbed the color excess and covering fraction measurements according to their respective errors and fit each of these realizations with Equation (15). The average and standard deviation in the best-fit parameters among these realizations are the values reported in Table 3. The 68% confidence intervals in the fits were computed as the 1σ dispersion in the fits to the realizations at each value of $E(B-V)$. These confidence intervals do not vary significantly over the range of $E(B-V)$ probed in our study for different assumptions of the attenuation curve, as the fits are forced to go through the measured data.

Finally, we can combine Equations (14) and (15) to obtain an expression for the observed LyC-to-UV flux density ratio in terms of $E(B-V)$ for the geometry illustrated in Figure 3:

$$\begin{eqnarray}{\left\langle \displaystyle \frac{S({\rm{LyC}})}{S({\rm{UV}})}\right\rangle }_{{\rm{obs}}} & = & {\left\langle \displaystyle \frac{S({\rm{LyC}})}{S({\rm{UV}})}\right\rangle }_{{\rm{int}}}\\ & & \times \langle \exp [-{\tau }_{{\rm{IGM}}}({\rm{LyC}})]\\ & & \times \exp [-{\tau }_{{\rm{CGM}}}({\rm{LyC}})]\rangle \\ & & \times \left[\displaystyle \frac{\exp [a\times \langle E(B-V){\rangle }^{b}]}{{10}^{-0.4k({\rm{UV}})\langle E(B-V)\rangle }}\right].\end{eqnarray} \tag{ 16 }$$

We have written this equation in terms of averages to emphasize that the relation applies for an ensemble of galaxies. Adopting values for the IGM and CGM opacities appropriate for the redshifts of interest allows us to use Equations (6), (12), and (16) in several useful ways, as we discuss in Section 7.4.3.

To uncover the extent to which the relationship between ${[S({\rm{LyC}})/S({\rm{UV}})]}_{{\rm{obs}}}$ and $E(B-V)$ may change over a larger range of $E(B-V)$ than accessible in the current sample, we examined how several of the other factors relevant in the escape of LyC photons, namely $N({\rm{H}}\,{\rm{I}})$ and $E{(B-V)}_{{\rm{los}}}$, vary with $E(B-V)$. As noted in Section 5.4, and shown in Figure 6 and Table 3, several significant correlations are found between these factors and $E(B-V)$. We fit $\mathrm{log}N({\rm{H}}\,{\rm{I}})$ versus $\mathrm{log}E(B-V)$ and $\mathrm{log}E{(B-V)}_{{\rm{los}}}$ versus $\mathrm{log}E(B-V)$ with linear functions. The best-fitting intercepts and slopes, and their uncertainties, are summarized in Table 3. We used these linear fits to estimate $N({\rm{H}}\,{\rm{I}})$ and $E{(B-V)}_{{\rm{los}}}$ for galaxies with $E(B-V)\lesssim 0.09$, and used Equation (6) to calculate $f{({\rm{LyC}})}_{{\rm{abs}}}$, and Equation (12) to calculate ${[S({\rm{LyC}})/S({\rm{UV}})]}_{{\rm{obs}}}$. In the next section, we examine how ${[S({\rm{LyC}})/S({\rm{UV}})]}_{{\rm{obs}}}$ varies with continuum reddening.

7.4.1. Calculating the Models

To calculate the predicted $\langle S({\rm{LyC}})/S({\rm{UV}}){\rangle }_{{\rm{obs}}}$ from our model, we first calculated the expected opacities of the IGM and CGM by considering the neutral gas column density distribution in the Lyα forest at $z\sim 2.0\mbox{--}2.8$, as well as the excess ${\rm{H}}\,{\rm{I}}$ absorption around $\sim {L}^{* }$ galaxies at the same redshifts, as quantified in Rudie et al. (2013). Using 1000 random realizations of these column density distributions, we calculate the attenuation of 900 Å photons to be $\langle \exp [-{\tau }_{{\rm{IGM}}}(900)]\times \exp [-{\tau }_{{\rm{CGM}}}(900)]\rangle =0.37$ for galaxies at z = 3.05, corresponding to the mean redshift of the 121 galaxies that contribute to the composite at $\lambda \lt 1150\,\mathring{\rm{A}}$. The dispersion in the transmissivity among the 1000 realizations at this wavelength is $\approx 0.33$, implying that the error in the mean transmissivity for the 121 galaxies is $\approx 8.1 \%$.

Next, we use the observed relations listed in Table 3 to estimate ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$, $E{(B-V)}_{{\rm{los}}}$ and $N({\rm{H}}\,{\rm{I}})$, for different values of $E(B-V)$. Thus, Equations (6), (12) and (16) can be written purely in terms of $E(B-V)$, the line-of-sight extinction and attenuation curves evaluated at the wavelengths of LyC and UV photons, respectively, and the intrinsic LyC-to-UV flux density ratio. To examine the possible range of models given the constraints imposed by our data, we calculated an "extrapolated" model which assumes that the neutral gas is optically thick to ionizing photons at all values of $E(B-V)$, and that the uncovered portions of the galaxies are effectively free of dust and gas. In this case, the model is given simply by Equation (16). In the "fiducial" model, we assumed that the neutral gas column density varies with $E(B-V)$ according to the relations listed in Table 3. Finally, in the "fixed" model, we assumed that the ratio $N{({\rm{H}}{\rm{I}})/E(B-V)}_{{\rm{los}}}$ is fixed, and that $E{(B-V)}_{{\rm{los}}}$ scales with $E(B-V)$ according to the relations listed in Table 3. All three models were computed assuming the combination of the updated Reddy et al. (2015) attenuation curve and either the SMC or Milky Way extinction curves, and are shown in Figure 11. Other combinations of the stellar attenuation and line-of-sight extinction curves yield quantitatively similar results (see discussion below). All models shown are computed by taking the UV luminosity density at 1500 Å and are cast in terms of the ratio of the observed and intrinsic 900-to-1500 Å flux density ratios.

Zoom In Zoom Out Reset image size

7.4.2. Behavior of the Models

7.4.3. Constraints on the Intrinsic LyC-to-UV Flux Density Ratio

We preface our discussion by noting that deep far-UV spectra of the kind used for our study would be useful to obtain over a broader dynamic range of $E(B-V)$, in particular focusing on the bluer galaxies that are more analogous to typical galaxies at $z\gtrsim 5$ and which may have larger ionizing escape fractions. Far-UV spectra may be used to directly measure the covering fraction of neutral gas at $E(B-V)\lesssim 0.09$ using the technique described in Section 5, and will enable more accurate estimates of the column densities and line-of-sight reddening for such galaxies. Nonetheless, in the following, we illustrate how such a model may be used to discern the mechanisms by which LyC photons escape from high-redshift galaxies.

The most direct application of our model is that it can be used to predict the average relative escape fraction of ionizing photons for a given stellar population model. For example, if we assume $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}=1/6$, a value typical of that adopted in previous studies (e.g., Siana et al. 2007; Nestor et al. 2011), then a sample of galaxies at z = 3 with $\langle E(B-V)\rangle =0.15$ would be predicted to have an observed relative escape fraction of $\langle S(900)/S(1500){\rangle }_{{\rm{obs}}}\approx 0.011$. The relative escape fraction corrected for attenuation by the IGM and CGM is 0.031, and the absolute escape fraction is $\approx 5 \%$. We emphasize that the model should only be applied to an ensemble of galaxies where the stellar populations and different configurations of the gas and stars from galaxy-to-galaxy, and the stochasticity in the Lyα forest along different sightlines, can be averaged out.

A second and potentially powerful application of our model is immediately obvious. Namely, with direct measurements of the average observed LyC-to-UV flux density ratio for an ensemble of galaxies with some average continuum reddening, $\langle E(B-V)\rangle$, one may invert Equation (16) to calculate $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$. The importance of direct empirical constraints on the intrinsic LyC-to-UV flux density ratio is underscored by the fact that the different flavors of stellar population models (e.g., arising from different ages and metallicities for a given star formation history, and/or the treatment of stellar rotation/binarity; Eldridge & Stanway 2009; Brott et al. 2011; Levesque et al. 2012; Leitherer et al. 2014) predict a factor of $\approx 5$ variation in the intrinsic LyC-to-UV flux density ratios of $\simeq 1.3\mbox{--}6.4$ (e.g., Nestor et al. 2013). The actual range in the predictions of the stellar population models are larger than a factor of $\approx 5$ if we also consider variations in the stellar initial mass function (IMF) and star formation history. Clearly, independent information from the rest-frame UV continuum spectra, rest-frame optical nebular emission line spectra, and/or stellar population modeling may be used to narrow the range of possible $\langle S({\rm{LyC}})/S({\rm{UV}}){\rangle }_{{\rm{int}}}$ (Steidel et al. 2016). Nevertheless, there is substantial utility in directly constraining the intrinsic LyC-to-UV flux density ratio without having to extrapolate the spectra of massive stars to wavelengths blueward of what is typically accessible from observations.

To demonstrate how our model can be used to constrain the intrinsic LyC-to-UV flux density ratio, we considered the recent measurement of the relative escape fraction of galaxies in our sample by C. C. Steidel et al. (2016, in preparation). For the subsample of galaxies with deep LyC spectra—where the sample has been cleaned of foreground contamination based on the ground-based spectroscopy (see Section 2)—these authors measured the ratio of the flux densities at 900 and 1500 Å, $\langle S(900)/S(1500){\rangle }_{{\rm{obs}}}=0.019\pm 0.002$. With this measurement, we can invert Equation (16) to compute $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$. To account fully for the measurement and random errors that enter into Equation (16), we performed a Monte Carlo simulation by perturbing each variable in the equation by its associated measurement or random uncertainty. As noted above, the mean combined IGM and CGM attenuation at the mean redshift of the 121 galaxies of $\langle z\rangle =3.05$ is $\langle \exp [-{\tau }_{{\rm{IGM}}}({\rm{LyC}})]$ × $\exp [-{\tau }_{{\rm{CGM}}}({\rm{LyC}})]\rangle$ $=\,0.370\pm 0.030$. The mean reddening of these galaxies is $\langle E(B-V)\rangle =0.212\pm 0.003$, assuming the Reddy et al. (2015) dust attenuation curve. For this attenuation curve, $k(1500\,\mathring{\rm{A}} )=8.687\pm 0.200$ as calculated in Paper I.

We performed 10,000 simulations where, each time, we randomly perturbed the IGM+CGM attenuation, the mean reddening, the dust attenuation curve, and the observed relative escape fraction according to normal distributions centered on the mean values with σ equal to the associated uncertainties. In addition, we randomly perturbed ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ and $E(B-V)$ for each composite according to their respective uncertainties and fit Equation (15) to these perturbed values. We then calculated $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$ for each of the 10,000 trials. The resulting distributions of $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$ are shown in Figure 12. We find $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}=0.199\pm 0.038$ and 0.240 ±0.047 assuming the Calzetti et al. (2000) and Reddy et al. (2015) attenuation curves, respectively (with an SMC line of sight extinction in both cases). The simulations indicate that we can use our methodology to constrain the intrinsic 900-to-1500 Å flux density ratio to $\lesssim 20 \%$ random error for the sample of galaxies considered here.

Zoom In Zoom Out Reset image size

In addition to this random error, there are several sources of systematic error in our estimate of $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$, including the choice of the line-of-sight extinction curve (e.g., SMC or MW) and the choice of the stellar reddening curve (e.g., Reddy et al. 2015 or Calzetti et al. 2000). The magnitude of the bias in $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$ that stems from these effects is $\lesssim 2 \%$ and $\approx 25 \%$, respectively. Additional systematic uncertainty may arise from the choice of IGM prescription. For example, Inoue et al. (2014) find a mean IGM transmission of $\approx 0.5$ for 900 Å photons emitted from a source at z = 3, a value at the high end of published measurements. On the other hand, the combined IGM+CGM opacity predicted by Rudie et al. (2013) is 0.37. This suggests a systematic uncertainty associated with the IGM prescription of $\approx 30 \%$.

The largest source of systematic uncertainty relates to our inability to fully resolve the ${\rm{H}}\,{\rm{I}}$ absorption. Specifically, our estimate of ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ may be a lower limit as a result of several effects, some of which are discussed in Section 4.3. First, the relatively low spectral resolution of our data may mask the presence of narrow absorption line systems with a higher covering fraction than that deduced from the spectra. Second, if the different velocity components of the optically thick ${\rm{H}}\,{\rm{I}}$ are not covering the same lines-of-sight, then the integrated covering fraction may be larger than the value found for the bulk of the absorbing gas. Third, we are unable to discern whether non-negligible columns of ${\rm{H}}\,{\rm{I}}$ and dust exist in regions that we would otherwise refer to as "clear" sightlines. Fourth, the ionized gas may not be completely bounded by ${\rm{H}}\,{\rm{I}}$ (i.e., the "average" configuration of the ${\rm{H}}\,{\rm{I}}$ and H ii may be intermediate between the radiation- and density-bounded scenarios). In this latter case, attenuation by dust in the ionized gas may deplete LyC photons before they escape the galaxy. For these reasons, the effective covering fraction, ${f}_{{\rm{cov}}}^{{\rm{eff}}}({\rm{H}}\,{\rm{I}})$, may be larger than the value derived from the composite spectrum and hence the calculated value of $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$ may be a lower limit.

To further evaluate the range of possible $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$ allowed by the data, we compared our estimate with the predictions of several stellar population synthesis models (Figure 13). We considered Bruzual & Charlot (2003) 0.28 ${Z}_{\odot }$ and ${Z}_{\odot }$ models with a Salpeter (1955) IMF (power law index of $\alpha =2.35$) over the mass range $0.1\leqslant {M}_{* }\leqslant 100$ ${M}_{\odot };$ and the most recent version of the Starburst99 (S99) models which adopt weaker stellar winds relative to the the previous version (Leitherer et al. 2014). We generated S99 models using the Kroupa (2001) IMF with different power law slopes of 2.30, 2.00, and 1.70, over the mass range $0.5\leqslant {M}_{* }\leqslant 100$ ${M}_{\odot }$, and with a metallicity of 0.14 ${Z}_{\odot }$. Lastly, we consider the "Binary Population and Spectral Synthesis" (BPASS; Stanway et al. 2016) models that include binary evolution, with a power law slope of the IMF of 2.35 over the mass range $0.5\leqslant {M}_{* }\leqslant 300$ ${M}_{\odot }$, and metallicities of 0.07 and 0.14 ${Z}_{\odot }$. All metallicities assume the solar abundances measured in Asplund et al. (2009).

Zoom In Zoom Out Reset image size

The effect of binary evolution (as in the BPASS models) is to allow efficient mass transfer from a massive star to its (secondary) companion, causing the latter to spin up quickly and, for a low stellar metallicity where the stellar wind is sufficiently weak, retain its angular momentum (Cantiello et al. 2007). Consequently, the secondary becomes fully mixed and undergoes quasi-homogeneous evolution (Yoon & Langer 2005; Brott et al. 2011; Eldridge & Stanway 2012) whereby the main sequence lifetime and effective temperature are increased. In such models, the ionizing photon flux per non-ionizing UV photon flux may increase by up to 60% relative to single star evolution (e.g., Stanway et al. 2016; see also Topping & Shull 2015). A similar effect can be achieved by decreasing the power law index of the IMF, which results in a larger mass fraction in more massive, hotter stars.

Figure 13 shows the limit on $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$ based on ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})\approx 0.96$ derived from the composite spectrum of all 933 galaxies in our sample. Our inference of the gas covering fraction favors an intrinsic LyC-to-UV flux density ratio that is larger than the value predicted from the Bruzual & Charlot (2003) stellar population models, and one that is consistent with those obtained with the newer versions of the S99 models that adopt weaker massive star stellar winds, a low stellar metallicity, and/or a flatter slope of the IMF, or those models that include the effects of binary evolution. For the reasons stated above, the covering fractions derived here are likely lower limits, and therefore it is instructive to gauge the degree to which the covering fraction may be higher and still be consistent with the predictions of stellar population synthesis models.

In particular, if the effective covering fraction is ${f}_{{\rm{cov}}}^{{\rm{eff}}}({\rm{H}}\,{\rm{I}})\approx 0.97$ or 0.98 instead of 0.96, then the inferred $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$ increases from 0.24 to 0.31 and 0.51, respectively. Note that the latter is $\approx 50 \%$ larger than the value predicted by the low-metallicity BPASS models and/or the low-metallicity and shallow-α IMF S99 models. Thus it is apparent, though ironic, that it is those galaxies with the largest covering fractions (as is the case for the average galaxy at $z\sim 3$) which provide the most sensitive constraints on the intrinsic LyC-to-UV flux density ratio, since even small increases in ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ will lead to very large changes in the inferred $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}$ (Figure 13). Steidel et al. (2016) suggest that stellar population models that include the effects of binary evolution are required to simultaneously reproduce the rest-frame far-UV continuum, stellar, and nebular features, and the rest-frame optical nebular emission line strengths measured for galaxies at $z\sim 2$. If we adopt such models as the most appropriate, then it implies that our estimate of ${f}_{{\rm{cov}}}({\rm{H}}\,{\rm{I}})$ for typical galaxies at $z\sim 3$ is not biased low by more than $\approx 1 \%$. Stated another way, the combined effects of additional ${\rm{H}}\,{\rm{I}}$ absorption and dust attenuation cannot deplete more than 25% of the photons exiting through what we would otherwise refer to as "clear" sightlines based on the depths of the ${\rm{H}}\,{\rm{I}}$ lines in the composite spectra.

This analysis underscores the fact that if the simple geometry of Figure 3 is valid on average for an ensemble of galaxies, then deep UV spectroscopy for large samples of typical star-forming galaxies at high redshift should enable tight constraints on the intrinsic LyC-to-UV flux density ratio. Hence, our framework provides an alternative and powerful independent method of constraining the intrinsic LyC-to-UV flux density ratio, which is currently one of the most uncertain parameters entering into the calculation of the absolute escape fraction of ionizing photons in high-redshift galaxies.

7.4.4. Deviations from the Assumed Geometry

The geometry depicted in Figure 3 may not apply for all galaxies, and, in reality, more complicated configurations of the dust, gas, and stars may be expected. In this context, the model described here provides a useful benchmark to investigate changes in geometry and hence alternate pathways for the escape of ionizing radiation. To illustrate this, we show in Figure 14 $\langle S(900)/S(1500){\rangle }_{{\rm{obs}}}$ as a function of $E(B-V)$ assuming the fiducial model for the SMC extinction curve. We have assumed $\langle S(900)/S(1500){\rangle }_{{\rm{int}}}=0.35$, a value typical of the BPASS models and the S99 models with a flat IMF slope, and IGM and CGM opacities relevant for galaxies at z = 3.05. Direct measurements of $\langle S(900)/S(1500){\rangle }_{{\rm{obs}}}$ and $\langle E(B-V)\rangle$ for an ensemble of galaxies may be used to ascertain the mechanisms by which LyC photons escape galaxies. For example, if these average values lie in the green region, then it suggests that, at least on average, those galaxies have either a smaller intrinsic 900-to-1500 Å flux density ratio than the one we have adopted here, or that the "uncovered" portions of the continuum are not completely free of dust and/or gas. Similarly, objects whose average values place them in the purple region are objects that may have a larger intrinsic 900-to-1500 Å flux density ratio or optically thin neutral gas that allows the escape of additional LyC photons. Finally, galaxies whose average values lie in the cyan region must have a larger intrinsic 900-to-1500 Å flux density ratio on average than the one assumed here. Further discrimination of the mechanisms by which LyC photons escape galaxies, or of the value of the intrinsic LyC-to-UV flux density ratio, can be accomplished with deep spectroscopy, such as that presented here. From this spectroscopy, we have found that typical $z\sim 3$ galaxies can be characterized by optically thick outflowing neutral ${\rm{H}}\,{\rm{I}}$, a finding that already narrows the range of possible avenues by which ionizing radiation can escape these galaxies.

Zoom In Zoom Out Reset image size