© ESO 2018

1 Introduction

Star-forming galaxies are ideal laboratories to understand how the early universe became reionized. Galaxies likely reionized the universe because quasars are too rare at high redshifts (Fontanot et al. 2012, 2014). Compact galaxies with intense star formation rates produce large amounts of ionizing photons which, under certain circumstances, escape the interstellar medium (ISM) and ionize the intergalactic medium (IGM). To reionize the universe, studies suggest that 10−20% of the ionizing photons produced by star-forming galaxies must escape galaxies (Ouchi et al. 2009; Robertson et al. 2013; Dressler et al. 2015). However, it has been challenging to detect Lyman continuum (LyC) radiation from individual galaxies.

Three unambiguous observations of ionizing photons have been reported at z ~ 3 (Vanzella et al. 2015; de Barros et al. 2016; Shapley et al. 2016; Bian et al. 2017). The observed escape fractions (near 50 %) are greater than those required to reionize the universe. Additionally, there are nine low-redshift (z < 0.3) galaxies with 1% $\le f_{\text{esc}}^{\text{LyC}}\le$ 13% (Bergvall et al. 2006; Leitet et al. 2013; Borthakur et al. 2014; Leitherer et al. 2016; Izotov et al. 2016a,b; Puschnig et al. 2017), and one recent detection at z = 0.37 with $f_{\text{esc}}^{\text{LyC}}\approx 46\%$ (Izotov et al. 2018). The low number of detections emphasizes the difficulty of detecting Lyman continuum emitters (LCEs).

Zackrisson et al. (2013) proposed two theoretical models to explain how ionizing photons escape galaxies. In the first scenario, low H I column densities allow Lyman continuum photons to pass through without being completely absorbed; this is called the density-bounded scenario (Jaskot & Oey 2013; Nakajima & Ouchi 2014). This scenario manifests itself as low H I column densities (< 10¹⁸ cm⁻²). In the second scenario, ionizing photons leak into the IGM through holes in the neutral gas (Heckman et al. 2011). This scenario is called the picket-fence model and relies on a patchy neutral gas. A patchy neutral ISM manifests as H I absorption lines with a covering fraction less than one. It is unclear which of these scenarios describes how ionizing photons leak from galaxies. Constraining neutral gas properties, especially the H I covering fraction and column density, is an effective way to disentangle how ionizing photons escape galaxies.

H I absorption lines in the rest-frame far ultraviolet (the Lyman series: 912–1026 Å) directly probe the H I covering fraction and column density. However, the Lyman series is challenging to observe for several reasons. First, this requires deep rest-frame far-ultraviolet observations blueward of Lyα, which is notoriously difficult to observe at low redshifts. Second, the Lyman series is unavailable at redshift z > 3 because the Lyα forest completely absorbs this region. Third, several observational obstacles need to be accounted for to measure the H I properties. In particular, foreground contamination (Vanzella et al. 2010), intervening absorbers, stellar continuum from the galaxy itself, ISM absorption lines, and, at low redshifts, Milky Way and geocoronal emission need to be identified.

Interstellar medium metal absorption lines (i.e., Si II 1260 Å, C II 1334 Å) are easier to observe than the Lyman series (Heckman et al. 2011; Alexandroff et al. 2015). However, using metal absorption lines to trace the H I assumes that metals directly probe the neutral gas. Recent studies indicate that ISM metal lines may have a factor of 2 times smaller covering fractions than H I absorption lines (Reddy et al. 2016b). As a result, metal absorption lines may not trace the H I.

In this article, we directly observe the Lyman series of individual high and low-redshift star-forming galaxies to determine their neutral gas properties, and to compare H I measurements to ISM metal properties. For the first time, we measure the H I properties of spectroscopically confirmed LyC emitters to determine which physical process enables ionizing photons to escape galaxies. A companion paper uses these observed H I properties to predict the escape fractions of ionizing photons (Chisholm et al. 2018, hereafter Paper II).

This paper is organized as follows: Sect. 2 describes the observational data. Section 3 defines the various models and equations used to fit the stellar continua and UV absorption lines. In Sect. 4 we use synthetic spectra to illustrate how accurately we recover H I column densities and covering fractions from observations. Section 5 discusses the measured H I covering fractions of the LyC emitters, the relation between the H I and Si II covering fractions, the effects of the assumed dust geometry, and comparisons to previous studies. We summarize our results in Sect. 6.

2 Observed data

We studied the neutral gas properties of a sample of 18 star-forming galaxies listed in Table 1. Our selection was driven by the need to observe the Lyman series, i.e., available rest-frame UV spectroscopy between Lyman-β and the Lyman limit. We selected the low-redshift galaxies observed with the Cosmic Origins Spectrograph (COS) on the Hubble Space Telescope (HST; Green et al. 2012) with this wavelength coverage. Given the sensitivity and wavelength range of the G130M grating on COS, the Lyman series is observable with a spectral resolution R≳1500 for galaxies at z > 0.18. Therefore, we selected the 15 galaxies at low redshifts (z < 0.36) with Lyman series observations from the COS archive. Each galaxy is at a redshift such that at least the Lyβ line is observable with the COS G140L or G130M gratings.

Many of our low-redshift galaxies were originally targeted to observe possible LyC emission, but only nine of them are confirmed LyC emitting galaxies. Four galaxies are from the Green Pea sample of Henry et al. (2015) and two are Lyman break analogs from Heckman et al. (2011). The other nine galaxies are known LyC emitters (J1503+3644, J0925+1409, J1152+3400, J1333+6246, J1442–0209, J0921+4509, Tol1247–232, Tol0440–381, and Mrk 54 from Izotov et al. 2016a,b; Borthakur et al. 2014; Leitherer et al. 2016). The nine leakers with COS/HST data were reduced using CALCOS v2.21 and a custom method for faint COS spectra (Worseck et al. 2016). The other COS/HST data were reduced with CALCOS v2.20.1 and the methods from Wakker et al. (2015). The five Izotov et al. (2016a,b) spectra were smoothed with a 5 pixel boxcar.

Finally, we also included three gravitationally lensed galaxies at z ≈ 3 (Stark et al. 2008; Koester et al. 2010). These galaxies are part of the Magellan Evolution of Galaxies Spectroscopic and Ultraviolet Reference Atlas (MEGaSaURA; Rigby et al. 2018) and were selected because they are the only galaxies in the MEGaSaURA sample with a signal-to-noise ratio (S/N) greater than 2 near the Lyman series. These are moderate resolution (R ~3000) spectra observed with the MagE spectrograph on the Magellan Telescopes (Marshall et al. 2008). Instead of the full galaxy names, we used the following short names for the two sources: SGAS J1226 = SGAS J122651.3+215220 and SGAS J1527 = SGAS J152745.1+065219.

The complete sample is summarized in Table 1, where we list a few host properties and (nominal) spectral resolutions of the observations. An upper limit on the [N II]/Hα ratio constrains the metallicity for SGAS J1527 (12 + log(O/H) < 8.5; Wuyts et al. 2012); for SGAS J1226 these lines are not accessible from the ground. Meanwhile, Stark et al. (2008) measured the metallicity of the Cosmic Eye using the R₂₃ index.

3 Ultraviolet spectral fitting methods and results

We now describe the theory and methods we used to fit the UV spectra including stellar continua, ISM absorption lines, and Milky Way absorption lines. The method is then applied to simulated data (in particular to determine systematic errors) and to the observed spectra.

3.1 Ultraviolet continuum and interstellar absorption line modeling

3.1.1 Adopted geometries and basic formulas

To describe the radiation transfer in the host galaxy, we adopted the classical picket-fence model with different assumptionson the geometric distribution of gas and dust. In practice we considered two cases: (a) the picket-fence model with a uniform foreground dust screen, and (b) a clumpy picket-fence model where dust is only in the neutral gas clumps. Two parameters describe the two models: the dust attenuation (here parametrized by E_B-V) and the geometric covering fraction of neutral gas (C_f), which is defined as the fraction of the total lines of sight of the emitted UV radiation that are intercepted by neutral gas in the direction toward the distant observer.

In (a) both the radiation emerging from the gas clumps (with geometric coverage C_f) and radiation directly escaping (1 − C_f) is attenuated by a uniform foreground dust screen. In (b), only a fraction, C_f, of radiation is processed through the gas clumps, imprinting interstellar absorption lines and attenuating the stellar continuum. The rest escapes unaltered. Gas within these clumps is assumed to be homogeneous and the interclump medium has a negligible column density of neutral gas, i.e., is assumed to be completely transparent.

These simple models have already been examined and assumed by other authors (e.g., Zackrisson et al. 2013; Borthakur et al. 2014; Vasei et al. 2016; Reddy et al. 2016b). We, therefore, only briefly list the main equations used in our spectral modeling. For a picket-fence model with a uniform foreground dust screen (a) the emergent flux F_λ is $$F_{\lambda }=F_{\lambda }^{\star }\times {10}^{-0.4k_{\lambda }{E}_{\text{B-V}}}\times \left(C_f\exp (-{\tau }_{\lambda })+(1-C_f)\right),$$(1)

where $F_{\lambda }^{\star }$ is the intrinsic stellar emission prior to alteration by the ISM, k_λ describes the attenuation law, and τ_λ is the optical depth of the interstellar absorption lines. For a picket-fence model, with a clumpy gas distribution (b), the emergent flux becomes $$F_{\lambda }=F_{\lambda }^{\star }\times {10}^{-0.4k_{\lambda }{E}_{\text{B-V}}}\times C_f\exp (-{\tau }_{\lambda })+F_{\lambda }^{\star }\times (1-C_f),$$(2)

where the second term describes the unattenuated, directly escaping radiation. This light is unattenuated because the holes are assumed to be free of gas and dust. For high covering fractions (C_f → 1) or low attenuations(E_B-V → 0), Eq. (1) and Eq. (2) are identical.

We define the residual flux, R, as the ratio of the flux density at the observed wavelength of the line to the continuum flux density. The R gives the fraction of light unabsorbed by the neutral gas. For saturated lines (τ_λ ≫ 1) the residual flux becomes $$R=1-C_f,$$(3)

for a uniform dust screen (a), and $$R=\frac{(1-C_f)}{{10}^{-0.4k_{\lambda }{E}_{\text{B-V}}}{C}_f+(1-C_f)},$$(4)

for case (b) (Vasei et al. 2016). We note that R of the uniform screen model is always greater than, or equal to, the R of the clumpy model (see the discussion in Sect 5.3).

Finally, the absolute escape fraction of radiation for monochromatic radiation, $f_{\text{esc}}=F_{\lambda }/F_{\lambda }^{\star }$, becomes $$f_{\text{esc}}=10-0.4k_{\lambda }{E}_{\text{B-V}}\times (C_f\exp (-\tau \lambda )+(1-C_f)),$$(5)

for the uniform dust screen (a), and $$f_{\text{esc}}=10-0.4k_{\lambda }{E}_{\text{B-V}}\times C_f\exp (-\tau \lambda )+(1-C_f),$$(6)

in clumpy geometry (b). It is important to note that the values of C_f and E_B-V differ a priori between the two model geometries, except for the case of high covering fractions (C_f → 1) or low attenuation (E_B-V → 0). They must be determined consistently from fits to the observed data, adopting the equations corresponding to the assumed geometry.

3.1.2 Fitting method

The main spectral region modeled here is the Lyman series, H I absorption lines from Lyβ to the Lyman break (≈912−1050 Å). In practice, owing to the lower S/N close to the Lyman break, the bluest Lyman line that we include is Ly6 (930 Å). Figure 1 emphasizes the complicated nature of the reddest Lyman series line, Lyβ: strong stellar continuum features (red dashed line) blend with the broad H I interstellar absorption (blue line) and the weaker O I interstellar absorption line (orange line). Consequently, bluer Lyman series lines (especially Lyγ) have fewer complications because of their simpler stellar continua.

We modeled ISM metal absorption lines from O I, O VI, Si II, C II, and C III (listed in Table 2). O I has a similar ionization structure as H I, such that O I absorption lines blend with all H I lines, except for O I 989 and 1039 Å. These two lines constrain the O I profile.

First, we fit the stellar continuum. We started with an initial stellar model using a linear combination of ten single-age stellar continuum models with ages of 1, 2, 3, 4, 5, 8, 10, 15, 20, and 40 Myr. We also used stellar continuum metallicities of 0.05, 0.2, 0.4, 1, or 2 Z_⊙. These spectra were drawn from the fully theoretical STARBURST99 library (S99; Leitherer et al. 1999), computed with the WM-Basic method (Leitherer et al. 2010), and have a spectral resolution R(S99) ≈ 2500. We chose the stellar continuum metallicity closest to the ISM metallicity (Table 1). The STARBUST99 models use a Kroupa initial mass function with a high (low) mass exponent of 2.3 (1.3), a high-mass cutoff at 100 M_⊙, and the stellar evolution tracks with high mass loss from Meynet et al. (1994). We fit for a linear combination of the stellar continuum flux, F_S99, as $$F^{\star }={\Sigma }_{i=1}^{10}{X}_i{F}_i^{99},$$(7)

where X_i are the linear coefficients for a given age (i) and $F_i^{99}$ are the STARBURST99 theoretical stellar continuum models for a given age.

Absorption lines of different ions are added using Voigt profiles defined by the velocity shift v, b-parameter, column density N, and C_f. The metal covering fraction is initially fixed to 1 to reduce the number of free parameters. We included the O I absorption lines that are directly blended with the Lyman series. Each ion is considered independent, as are its parameters. For each galaxy, we tested whether including the remaining metal lines listed in Table 2 improves the Lyman series fits. If, by eye, they did not improve the fits, we did not include the lines. Finally, the STARBURST99 models are convolved with the nominal spectral resolution of the observations.

We accounted for dust attenuation using the attenuation law from Reddy et al. (2016a), a uniform dust screen model, and fitting for the dust attenuation parameter (E_B-V; similar to Chisholm et al. 2015). The linear combination of stellar continuum models, interstellar absorption lines, and dust attenuation produces the final fitted spectrum.

We fit the data using an IDL routine based on the nonlinear least squares method, MPFIT (Markwardt 2009). The MPFIT method returns the best fit and errors for E_B-V, b, v, N, and C_f of each ion. For observed spectra, the first step consists of masking the ISM absorption lines and the contaminating Milky Way absorption lines, geocoronal emission, and intervening absorbers. We applied these masks on the data and fit for the linear combinationof dust-attenuated STARBURST99 models. In a second step, we fixed the stellar continuum and fit for the ISM absorption lines and Milky Way absorption lines adjacent to the ISM lines. We fit for all of the observed Lyman series lines up to Ly6, provided that they are not near geocoronal emission or intervening absorbers, and do not have a S/N below one. Since the simulated data do not contain the extra complications of adjacent Milky Way and geocoronal lines, all parameters are simultaneously fit in one step.

3.2 Fitting simulated data

We tested our fitting method with both noise-free and noisy synthetic data to determine how well the estimated parameter errors represent the actual errors. Additionally, we tested how these errors depend on the S/N and spectral resolution. In Appendix A we fully describe the generation of synthetic spectra, but here we summarize the steps. The synthetic spectra were produced for two scenarios: one assuming the picket-fence model (log (N_H I[cm⁻²]) = 20, C_f = 0.9) and one describing a uniform ISM in the density-bounded regime (log (N_H I[cm⁻²]) = 17.57, C_f = 1). Both scenarios correspond to an escape fraction f_esc = 0.1, regardless of the chosen dust distribution, since we set E_B-V = 0 (Eqs. (5) and (6)).

We created synthetic spectra using the parameters of the picket-fence and density-bounded regime (see Table A.1) for 7 different spectral resolutions between R = 600−15 000. For each of these 14 set-ups (7 spectral resolutions and two scenarios) we created 50 different realizations by adding various sets of random Gaussian noise. The level of Gaussian noise was chosen to produce a final S/N per pixel between 2− 50, in 7 total S/N steps. In total we created 98 different configurations (2 different scenarios, 7 spectral resolutions, and 7 S/Ns), and each configuration has 50 individual spectra, for a total of 4900 synthetic spectra. Figures B.1 and B.2 show synthetic spectra and the best-fit models for high and low spectral resolution (R = 15 000 and R = 600) for three S/N configurations (noise-free, 10 and 2).

We fit the synthetic spectra with the methods outlined in Sect. 3.1.2. From this fitting we define two types of error:

The statistical error (Err_stat). This is the individual errors returned by MPFIT.

The systematic error (Err_syst). This is the deviation of the 50 parameter estimates of each scenario from the actual parameter value that created the original line profile. The Err_syst is therefore a function of S/N and resolution. Err_syst is calculatedas $$Er{r}_{\text{syst}}=\sqrt{\frac{1}{50}{\sum }_{i=1}^{50}{(x_i-x_T)}^2},$$(8)

where x_i is the estimated parameter and x_T is the true parameter value. The systematic error accounts for differences between the estimated parameters and the actual parameters that are unaccounted for by MPFIT. We express Err_syst as a percent error, defined as the deviation of the measured parameter from the true parameter divided by the true parameter value. This fractional error is more broadly applicable to a variety of measurement values.

By definition, Err_syst is not included in the statistical uncertainties reported by MPFIT, but it critically describes the ability of MPFIT to recover the actual parameter values. We account for Err_syst by defining the total error of each parameter (Err_tot) as the quadratic sum of Err_syst and Err_stat as follows: $$Er{r}_{\text{tot}}=\sqrt{Er{r}_{\text{syst}}^2+Er{r}_{\text{stat}}^2}.$$(9)

The systematic errors derived here are included in the errors of the H I covering fraction and column density (Cols. 3 and 6 of Table 3).

3.3 Fitting observed data

Using the fitting method described in Sect. 3.1.2, we fit the observed spectra of the 18 galaxies in our sample. The resulting fit parameters E_B-V, C_f, v, b, and N are listed in Table 3, and the corresponding best-fit spectra are shown in Appendix B. The reported uncertainties on the H I column densities and covering fractions account for the systematic errors derived from simulations.

We note that Tol0440–381 and Mrk54 have much lower redshifts than the rest of the sample. Consequently there are more Milky Way absorption lines near the Lyman series. This high density of Milky Way Lyman series absorption lines means that we cannot accurately fit the Lyman series with Voigt profiles. Rather, we used nonparametric fits to describe the H I properties (see Sect. 3.3.1).

The H I column densities are affected by large uncertainties because the low S/N and insufficient resolution of the spectra do not allow a reliable estimate of N_H I given the saturation of the H I absorption lines (see Sect. 4.1). Therefore direct N_H I measurements from the Lyman series are not reliable and an indirect method is used (Sect. 3.3.2). To examine whether saturation and resulting degeneraciesaffect the C_f or metal column densities that we derived, we refit the spectra fixing the H I column density to N_H I = 10¹⁸ cm⁻². This N_H I corresponds to an optically thick portion of the curve of growth for the observed Lyman series lines (from Lyβ to Ly6). We find that all the fit parameters (except for b(H I)) are consistent within 1σ with the results listed in 3. The fact that b changes is consistent with b being degenerate with N_H I at these column densities. This indicates that saturated H I absorption does not affect the measured velocities and covering fractions.

3.3.1 Nonparametric measurements of the H I covering fraction

Fitting the Lyman series using the radiative transfer equation (Eq. (1)) is one way to measure the H I covering fraction. However, this assumes that the absorption profiles follow a single Lorentzian velocity distribution. Observed absorption profiles arising from galactic outflows are highly non-Lorentzian (Heckman et al. 2000; Pettini et al. 2002; Shapley et al. 2003; Weiner et al. 2009; Chisholm et al. 2017b). Consequently, we also measured C_f from the residual flux of the Lyman series lines after removing the stellar continuum (see Eq. (3)). This assumes a uniform dust geometry and that the Lyman series is fully saturated. The Lyman series is saturated from N_H I≳10¹⁶ cm⁻² (for Lyβ to Ly6). Importantly, this nonparametric approach does not assume a velocity distribution of the H I gas and accounts for arbitrary line profiles. In other words, the nonparametric approach does not assume that the gas has the same covering fraction at every velocity (cf. discussion in Vasei et al. 2016).

The covering fraction is derived as the maximum of (1 − F_Gau) in a velocity range chosen by eye near the deepest part of each Lyman series line. The F_Gau is the stellar continuum removed flux, modified by a Gaussian kernel centered on zero with standard deviation corresponding to the error array. We measured (1 − F_Gau) 1000 times, where each time the flux value of each pixel is determined from a different noise distribution. We then took the median and standard deviation of this distribution as the C_f value and uncertainty for each Lyman series transition. Table C.3 lists the C_f derived from the residual flux of each Lyman series transition in each galaxy. We then defined the C_f (H I) from the residual flux as the error weighted mean (M_W) of the i observed Lyman series transitions and the C_f error as the error on M_W as $$\begin{array}{l}{M}_W=\frac{{\sum }_{i=1}^n{C}_f{}_i\times {\omega }_i}{{\sum }_{i=1}^n{\omega }_i}\text{\hspace{0.17em}with\hspace{0.17em}}{\omega }_i=\frac{1}{{\sigma }_i^2},\hfill \\ {\sigma }_{M_W}=\sqrt{\frac{1}{{\sum }_{i=1}^n{\omega }_i}}.\hfill \end{array}$$(10)

The corresponding values, denoted as C_f depth, are reported in Table 3. The C_f depth values are consistent, within the errors, with the C_f values derived from the fits in Sect. 3.3. Except for Tol0440–381 and Mrk54, where we did not fit the absorption lines with the method in Sect. 3.1.2, we computed the final C_f (H I) as the weighted mean of the C_f depth and the fitted C_f values (column 8 in Table 3).

3.3.2 Indirect measurements of the H I column density

Since the Lyman series lines are saturated, but not damped, and the spectra have insufficient spectral resolution and S/N direct H I column density measurements are largely unconstrained, as discussed further in Sect. 4.1. Therefore, indirect methods of measuring N_H I are needed. We used the O I absorption lines, constrained when possible by the unsaturated O I 1039 Å absorption line, to derive the O I column density. Using the known metallicity (12 + log(O/H)) of the galaxy, we then indirectly inferred the hydrogen column density. This approach assumes that the emission-line based oxygen abundance, tracing the chemical composition of the ionized gas, is identical to that of the neutral gas. If the metallicity was lower in the neutral gas, for example, because of the presence of some pristine gas (see, e.g., Lebouteiller et al. 2013, and references therein) or incomplete mixing of the metals (see Sect. 5.2), the resulting N_H I would be higher (and more saturated) than inferred here.

The O I fit parameters b, v, N_O I, and the H I column density derived using the metallicities from Table 1, are listed in Table 4. A curve of growth analysis indicates that O I 1039 Å saturates at N_O I > 10^16.5 cm⁻², whereas all but one of our N_O I are less than 10^16.2 cm⁻². If the O I 1039 line is saturated, then the N_H I values in Table 4 would underestimate the actual N_H I (i.e., N_H I is not lower than the quoted values). Figure B.3 demonstrates that O I 1039 Å is resolved from the C II 1036 Å and the O VI 1038 Å absorption lines. However, O I is undetected for J0925+1409 (see Fig. B.5), and we used O I 1302 Å for Tol0440–381 and Mrk54. We could not determine N_H I for the three MEGaSaURA sources, since the O I line was either saturated, or there was not a literature metallicity. The resulting H I column densities are found to be log(N_{H I}) ~ 18.6−20 cm⁻².

We fixed the O I covering fractions to the H I values (i.e., C_f (O I) = C_f(H I)). This is plausible because O I and H I have similar ionization potentials and their ionization fractions are locked together by charge exchange. We examined how tying the covering fractions impact the derived O I column densities. As an extreme case, we fixed the O I covering fraction to 1 for J1152+3400, which has one of the lowest C_f (H I), and J0921+4509, which has the lowest N_H I derived from O I. We found that log(N_OI) is reduced by 0.30 and 0.15 dex, respectively, which is comparable to the N_O I errors. Therefore, we conclude that tying C_f(O I) to C_f(H I) does not drastically change the measured O I column density.

3.3.3 The Si II covering fraction

We measured the Si II covering fraction (C_f(Si II)) using the Si II 1190 Å doublet and the Si II 1260 Å singlet. For the Si II 1260 Å singlet, C_f(Si II) was derived using the same procedure as for the residual flux measurements of the Lyman series (Sect. 3.3.1). A different approach, detailed in Chisholm et al. (2017b), was adapted for the doublet. As the two transitions share the same C_f and the ratio of the two optical depths is given by the ratio of their oscillator strengths, f, the velocity-resolved covering fraction was measured from a system of equations (Hamann et al. 1997) as $$C_f(v)=\frac{F_W{(v)}^2-2F_W(v)+1}{F_S(v)-2F_W(v)+1},$$(11)

where F_W is the continuum subtracted flux of the weaker doublet line (Si II 1190 Å) and F_S is the continuum subtracted flux of the stronger doublet line (Si II 1193 Å). This method accounts for the possibility that the lines are not saturated. We measured the covering fraction 1000 times by varying the flux in a similar way to the residual flux measurements of the Lyman series. The median and standard deviation of this distribution is taken as the C_f(Si II) value and error (Table 5). Using the Si II 1260 Å singlet assumes that the line is saturated, whereas the doublet method accounts for unsaturated absorption. We derived statistically consistent values using both methods (all within 1σ) for all the galaxies for which we can measure C_f(Si II) with both methods. This indicates that Si II 1260 Å is indeed saturated.

3.4 Effect of ultraviolet attenuation law on covering fraction

A priori themeasurements made here also depend on the dust attenuation law used. In this study, we used the attenuation law from Reddy et al. (2016a) because it is defined blueward of Lyα. We refit J1503+3644 and GP0911+1831, which are two galaxies with high E_B-V and low C_f, using a Small Magellanic Cloud (SMC) attenuation law¹ still assuming a uniform dust foreground. The SMC law is significantly steeper than the Reddy et al. (2016a) law. With the SMC dust law, we measured C_f(H I) of 0.653 ± 0.109 and 0.744 ± 0.156, respectively. These are consistent, within 1 σ, with the C_f (H I) estimated using the Reddy et al. (2016a) attenuation law (0.634 ± 0.109 and 0.731 ± 0.150). The E_B-V values change based upon the attenuation law used to match the observed continuum, but these changes do not affect the measured C_f. We therefore conclude that the adopted attenuation law does not significantly impact the measured covering fractions.

4 Recovering H I properties from simulated spectra

We simulated synthetic spectra and fit these mock H I lines with the method in Sect. 3.2. Comparing the fitted results with the parameters that created the spectra characterizes how accurately the method returns the H I parameters. We discuss these results in the context of the H I column density (Sect. 4.1) and the covering fraction (Sect. 4.2). This discussion illustrates that C_f is accurately measured for most resolutions and S/Ns, while the H I column density has large uncertainties. These simulations are especially helpful for planning future observations by determining the S/Ns and resolutions required to accurately measure the covering fractions of LyC emitters.

4.1 H I column densities

The synthetic spectra allow us to quantify how accurately our method reproduces the H I properties. For S/Ns (per pixel) less than 10 and resolutions less than 3000, the simulations have N_H I percent errors greater than 300% (Table C.1). Even at higher resolutions (R = 15 000), the percent error is greater than 200%, unless the S/N is greater than 5. For the lower S/Ns typical of our observations, we measured order of magnitude systematic uncertainties on N_H I. These large uncertainties are inherent because the Lyman series transitions saturate for H I column densities between N = 10¹⁶ cm⁻² to N = 10²² cm⁻². For these so-called Lyman limit systems, high quality and very high-resolution spectra (R ~ 30 000) are needed to constrain N_H I with Voigt fitting methods (see, e.g., O’Meara et al. 2007). Therefore, we conclude that N_H I cannot be directly fitted from the Lyman series absorption lines. However, the O I absorption lines included in our fits remain unsaturated for N_O I < 10^16.5 cm⁻² and do not suffer from these large uncertainties. Consequently, the neutral column density is most accurately inferred by converting N_O I into N_H I using the gas-phase metallicity, as in Sect. 3.3.2.

4.2 H I covering fractions

Conversely, the simulations show that C_f(H I) has a low percent error, under typical observing conditions. At R > 1500 and S∕N> 5, the C_f systematic percent errors are less than 6% of the measured value (Fig. 2, Table C.2). Therefore, the neutral gas covering fractions are accurately recovered from our observational conditions.

The James Webb Space Telescope (JWST) will accurately measure C_f of metal absorption lines from high-redshift galaxies. The JWST’s Near-Infrared Spectrograph (NIRSpec) is expected to have R ~ 3000 (1000) in the high- (medium-) resolution configurations. This means that the systematic errors will be 3% (6%) of the measured C_f for S∕N = 5 observations, illustrating the feasibility of measuring C_f from high-redshift galaxies with JWST.

Fig. 1 Total fit (red solid line) of the Lyβ absorption line for the galaxy GP 1244+0216 (black data). The contributions from the stellar continuum (dashed red line), Lyβ absorption line (blue solid line), and O I absorption line (orange solid line) blend together. O I absorption lines at 989 and 1039 Å are unblended and robustly constrain the O I profile. All three components need to be accounted for when fitting the Lyman series.

Fig. 2 Color map of the systematic percent error of the covering fraction as a function of the resolution and S/N. The synthetic spectra are created with R < 120 km s−1 combine a theoretical stellar continuum spectra with R(S99) = 120 km s−1 and absorption lines of spectral resolution R. The covering fraction is recovered to within 5% of the estimated parameter for all observations within the dark blue region.

5 Covering fraction of LyC emitters and comparison sources

We now examine the derived covering fractions from the Lyman series and the Si II absorption lines. Then we discuss different geometrical model assumptions and compare our results to earlier work.

5.1 Leakers have low neutral gas covering fractions

The H I covering fraction describes the porosity of the neutral gas and demonstrates whether the neutral gas is clumpy. A smaller C_f (H I) means that there are more low-density channels for ionizing photons to escape through.

Our sample has H I covering fractions ranging from 0.50 to unity (Fig. 3). Only 5 of the 18 galaxies (28%) have an H I covering fraction consistent with unity at 1σ. The low C_f (H I) values are likely because the sample is biased: 15 of these 18 galaxies were targeted as potential LyC leaker or for being particularly strong line emitters. Since a nonunity C_f is a possible LyC escape mechanism, it is not too surprising that many of these galaxies have low C_f (H I). The three galaxies that were not targeted as LyC leakers (the MEGaSaURA galaxies) have C_f (H I) consistent with unity. In contrast, the galaxies with the highest confirmed escape fractions of ionizing photons have the lowest C_f (H I) values (Fig. 3). The leakers have a median C_f (H I) = 0.62 ± 0.10, while the unknown leakers have a median C_f (H I) = 0.95 ± 0.10. Lyman continuum emitters have low H I covering fractions, which allows LyC photons to escape through low-density channels.

The high N_H I values, estimated from the O I column densities (Table 4), further emphasize that ionizing photons must escape through holes in the H I. N_H I, calculated from N_O I, is greater than 10^18.6 cm⁻² for the entire sample. Even without converting into N_H I, the high N_O I values require unphysically high metallicities (12 + log(O/H) > 10) for ionizing photons to escape through low-density regions. At these column densities the Lyman series and Lyman limit are saturated. Ionizing photons cannot pass through the neutral gas unabsorbed. In other words, even LyC galaxies have optically thick H I; ionizing photons must escape through low-density channels.

Low C_f(H I) values indicate that the escape of LyC photons is dominated by a patchy ISM or the picket-fence model. However, we find that a low covering fraction is not the only parameter leading to a high f_esc. For example, Mrk 54 has a C_f (H I) ~ 0.5 and a f_esc < 1%. Dust crucially impacts the LyC escape fraction by removing ionizing photons (see Eq. (5), which assumes a uniform dust screen). Consequently, the escape fractions cannot simply be inferred from the measured covering fractions, but requires a joint determination of C_f and E_B-V for a given set of geometrical assumptions (e.g., cases (a) or (b) from in Sect. 3.1.1). This is discussed in depth in Sect. 5.3 and in Paper II, in which we show that the LyC escape fractions derived using the Lyman series absorption features are consistent with the directly observed escape fractions.

Using Eq. (5), we can predict which of the nine galaxies in the sample without measured LyC emission should emit ionizing photons. To emit LyC photons, both a low E_B-V and a low C_f are required. J0926+4427, J1429+0643, and GP1054+5238 are the best candidates in our sample to leak ionizing photons (in order of most likely to emit ionizing photons). Follow-up observations should expect to find f_esc values between 0.02−0.06.

5.2 H I and Si II covering fractions

Low-ionization metal absorption lines are often used as proxies for the neutral gas covering fraction (Shapley et al. 2003; Jones et al. 2013; Alexandroff et al. 2015; Trainor et al. 2015). Metal lines are in redder portions of the spectra. Consequently, instruments can more efficiently observe metals lines at low redshifts and the Lyman forest does not obscure metal lines at high redshifts. Therefore, metal absorption lines are ideal probesof neutral gas properties when the Lyman series is not observed. However, recent observations of z ~ 3 stacked spectra have suggested that metal absorption lines have covering fractions a factor of two smaller than H I absorption lines (Reddy et al. 2016b). Here, we test whetherC_f(Si II) traces C_f(H I).

The C_f(Si II) is systematically lower than C_f(H I) and has a mean offset of 0.25 (Fig. 4). C_f(H I) and C_f(Si II) are linearly related at the 3σ significance level (p-value < 0.001; Pearson’s correlation coefficient of 0.79) as $$C_f(\text{H}\text{I})=(0.63\pm 0.12)\times C_f(\text{Si}\text{II})+(0.44\pm 0.07).$$(12)

While C_f(Si II) is not equal to C_f(H I), Eq. (12) estimates C_f(H I) from C_f (Si II). This empirically derived relation estimates C_f(H I) from measurements of the metal line covering fraction (here from Si II lines), even if the Lyman series lines are not accessible. This emphasizes a clear practical advantage of the uniform dust screen geometry: the stellar attenuation can be estimated without directly measuring C_f (H I). The stellar attenuation then defines the continuum level from which C_f (metal) can be estimated and, in turn, C_f(H I) estimated from Eq. (12). Consequently, C_f(H I) can be estimated using C_f(Si II) and stellar continuum fits of high-redshift galaxies, where the Lyman series is not directly observable. This is not possible in the clumpy geometry because it requires C_f(H I) to fit the stellar continuum. In Paper II we used the Lyman series absorption lines and dust attenuation to reproduce the observed Lyman continuum escape fraction.

There are numerous reasons why C_f(Si II) could be related to, but not equal to, C_f(H I) (see Reddy et al. 2016b, for an in-depth discussion). First, narrow, high covering fraction, metal lines may be unresolved in low-resolution spectra (~1500). While this is possible for low-resolution data, it is less likely for the high-resolution COS data. The Si II lines could possibly be unsaturated, but the doublet method accounts for this possibility and does not remove the systematic offset. Alternatively, the Si II ionization potential overlaps with, but is not equal to, H I. Therefore, the Si II and H I gas may trace similar, but not equal, gas. However, Reddy et al. (2016b) have found that the Si II and H I line profiles and central velocities are similar (also see Chisholm et al. 2016). This indicates that the two transitions are comoving and trace similar gas. Finally, the Si II lines are also potentially affected by scattering and fluorescent emission in-filling, although it is difficult to predict whether the amount the emission increases the covering fraction (cf. Prochaska et al. 2011; Scarlata & Panagia 2015).

Another explanation is that the neutral gas and metal-enriched gas are not fully mixed. In this situation, Si II only probes metal-rich regions, while the H I gas probes both high and low metallicity gas. Therefore, the neutral gas covers more of the background source than the metals alone. At lower metallicities this is exaggerated because there is less metal-rich gas. This leads to fewer high-density metal regions to absorb the background continuum, and systematically reduces the metal covering fraction. To test the effect of metallicity on the relation between the neutral and metal covering fractions, we fit a multiple linear relationship between C_f (H I), C_f (Si II), and 12 + log(O/H). We find a significant (3σ, p-value < 0.001, Pearson’s correlation coefficient of 0.79) trend between these three parameters that scales as $$\begin{array}{l}{C}_f(\text{H}\text{I})=\left(1.8\pm 0.8\right)-\left(0.18\pm 0.10\right)\times \left[12+\log (\text{O/H})\right]\\ +\left(0.75\pm 0.16\right)\times C_f(\text{Si}\text{II}).\end{array}$$(13)

While Eq. (13) is only marginally more significant than the simpler relation between C_f (Si II) and C_f (H I), it physicallyexplains the relation between the neutral and metal covering fractions.

The difference between Eqs. (12) and (13) is the metallicity dependence of the correction added to C_f (Si II). As the metallicity increases, the required correction to convert C_f(Si II) into C_f(H I) decreases. While at lower metallicities, the factor added to C_f(Si II) must be larger. Physically, this means that at lower metallicities the Si II traces a smaller fraction of the total H I. Therefore, the linear relationship between C_f(H I) and C_f(Si II) may arise because metal-rich clumps cover the background source differently than the H I gas.

The empirical relationship, Eq. (12), is recovered from Eq. (13) when the median metallicity of the sample (8.06) is used. However, Eq. (12) would change if the median 12 + log(O/H) of the sample changes. This suggests that the empirical relationship between C_f(H I) and C_f(Si II) is only constrained near the metallicities of the sample. For divergent 12 + log(O/H) values, the metallicity should be accounted for when estimating C_f(H I).

Fig. 3 Histogram of the H I covering fraction (Cf(H I)) from the 18 galaxies in our sample. We split the sample into galaxies known to leak LyC photons (blue) and galaxies without measured LyC emission (black). The leakers have the lowest Cf (H I) and the unknown emitters have the highest Cf(H I).

Fig. 4 Comparison of the Si II and H I covering fractions. The black solid line shows a one-to-one relation and the dotted line is the fitted relation (Eq. (12)). The average offset between Cf (H I) and Cf (Si II) is 0.25.

5.3 Is the dust geometry clumpy or uniform?

In Sect. 5.1, we emphasized that dust heavily contributes to the escape of ionizing photons. However, since the impact of dust depends on the assumed geometry, we explore the two cases described above (Sect. 3.1.1), namely (a) a uniform dust screen and (b) a clumpy geometry dust model. To test the impact of these assumed geometries on C_f and E_B-V, we modify our fitting routine from (a) a uniform dust screen model to (b) a clumpy model. We refit two galaxies: J0921+4509 and J1152+3400, which have an average C_f(H I) and low C_f (H I), respectively. Adopting a clumpy geometry (b), we derive E_B-V = 0.236 and 0.239 (versus E_B-V = 0.224 and 0.134 in model (a)) and C_f (H I) = 0.976 and 0.912 (versus C_f = 0.754 and0.506 in model (a)), respectively. These values and Eq. (6) predict a LyC escape fraction of f_esc = 0.024 and 0.088 at λ = 912 Å, which is consistent, within 1σ, with the values derived using the uniform dust screen model (f_esc = 0.017 ± 0.007% and 0.092 ± 0.027%).

As just shown, using the uniform dust screen leads to a lower E_B-V and C_f (H I) compared to the clumpy model. In the clumpy model, all of the light that escapes through holes in the neutral gas is unattenuated. Therefore, the holes in the neutral gas must be smaller, while any continuum passing through the neutral gas must be more heavily attenuated to match the observed flux. C_f (H I) and E_B-V must take values that fit both the geometry and the observed flux (extending to Lyman continuum). Both quantities are therefore model dependent. However, disentangling the actual geometric model is not mandatory to predict f_esc with C_f (H I) and E_B-V, as long as the geometry is consistently modeled (see Paper II).

On the other hand, the observed residual flux of the H I absorption lines provide clues to the geometric model. Equation (4) shows that the residual flux of the H I absorption lines are wavelength dependent in the clumpy geometry. Therefore, if the residual flux of the individual Lyman series absorption lines changes with wavelength (for example from Lyβ to Ly6), then the data would suggest a clumpy dust geometry. Figure 5 shows the relationship between the residual flux of a saturated absorption line (R) and the neutral covering fraction for various values of A_λ = k_λE_B-V in a clumpy geometry. As A_λ increases (larger E_B-V or a bluer transition), the residual flux decreases as dust removes more continuum photons. In clumpy dust geometry, R is less than 1–C_f if there is dust. For a covering fraction of 0.8, the measured R is 0.2 when A_λ is 0, 0.4 when A_λ is 0.1, and 0.7 when A_λ is 0.2 (Fig. 5). The C_f estimated from 1–R is underestimated in the clumpy dust geometry case.

Since the residual intensity increases with A_λ and k_λ, this effect, if present, should be detectable in the Lyman series absorption lines. Nevertheless, the differences between theresidual flux of the Lyman series should be small and most easily detected in galaxies with large E_B-V. Table 6 shows the expected residual flux of saturated Lyman series lines, assuming an f_esc = 0.10 (C_f = 0.9 in the clumpy geometry). When E_B-V is 0.1, the difference between the Lyβ and Ly6 residual flux is less than 0.02. For E_B-V = 0.30, the offset rises to 0.05. While we do not measure a trend in our observations (Table C.3), the measured residual intensities have large uncertainties. Our simulations imply that S∕N = 30 and R = 15 000 are required to distinguish the two geometries. We cannot determine the dust geometry with the current data and both a uniform dust screen or clumpy geometry are allowed by the observations.

Fig. 5 Residual flux R of a saturated absorption line from Eq. (4) for various values of Aλ = kλ EB-V in a clumpy ISM geometry. When dust is present, the residual flux is higher than (1 − Cf). However, in a uniform dust screen geometry, the residual flux is always R = (1 − Cf). The covering fraction derived from the observed residual flux and a clumpy dust geometry is underestimated when dust is present.

5.4 Comparison with other studies

Reddy et al. (2016b) recently determined the H I covering fraction of z ~ 3 Lyman break galaxies from Lyman series fits. Analyzing stacked spectra, they found high C_f(H I) ~ 0.92−0.97. They also examined the relation between C_f(H I) and the C_f of metal lines (Si II, C II, and Al III), finding a much lower covering fraction for the metals than for H I. The differences between their C_f(metal) and their C_f(H I) are nearly twice as large as we find for C_f(Si II) in Fig. 4 (their C_f(metal) is 0.4–0.6 smaller than C_f(H I), versus our average of 0.25).

This discrepancy arises from an inconsistency in determining C_f (metal). To fit their C_f(H I), Reddy et al. (2016b) adopted a clumpy geometry with dust only in the H I. The continuum level used to obtain the normalized composite spectra in their study is determined as $$F_{\lambda }^{\text{cont}}=F_{\lambda }^{\star }{10}^{-0.4E_{\text{B-V}}{k}_{\lambda }}\times C_f(\text{H}\text{I})+F_{\lambda }^{\star }[1-C_f(\text{H}\text{I})].$$(14)

In this geometry, the residual flux of a saturated metal absorption line is $$R(\text{metal})=\frac{F_{\lambda }}{F_{\lambda }^{\text{cont}}}=\frac{1-C_f(\text{metal})}{{10}^{-0.4E_{\text{B-V}}{k}_{\lambda }}\times C_f(\text{H}\text{I})+1-C_f(\text{H}\text{I})}.$$(15)

We note that Eq. (15) reduces to Eq. (4) when C_f (metal) = C_f(H I). We stress that the continuum normalization impacts the determination of C_f from R because R depends on the assumed dust geometry. In the uniform dust screen geometry (a; the geometry that we predominately use, i.e., in Fig 3), the covering fraction is simply 1–R. Conversely, in the clumpy model (b), the covering fraction of a saturated absorption lines is not simply 1–R. The covering fraction cannot be read off from the profile; a measured E_B-V is required to determine C_f (see the large differences in Fig. 5; Vasei et al. 2016). However, Reddy et al. (2016b) measured C_f (metal) as 1–R from the continuum normalized spectra. Meanwhile, they determined C_f (H I) from the full fits of their spectra, which properly accounts for their assumed geometry. This means that C_f (metal) and C_f (H I) are calculated with different geometric assumptions.

The inconsistency in the derivation of C_f(metal) and C_f(H I) produces the large mismatch between our observations and Reddy et al. (2016b). We correct the observed residual flux from Reddy et al. (2016b) to recompute consistently the C_f (metal) values using Eq. (15) and their measured E_B-V (using their SMC extinction law) and C_f(H I) values (Fig. 6). The large offset between the Si II 1260 Å and H I is halved (compare the dashed blue line to the dashed black line), although C_f(metal) is still systematically smaller by ~0.2. This offset is similar to the offset that we measure in Fig. 4 (green dashed line in Fig. 6). The discrepancy between C_f(metal) and C_f(H I) emphasizes that C_f can only be read off from the residual flux of absorption lines in the uniform dust screen geometry because this continuum normalization accounts for the dust attenuation. In other geometries, the dust attenuation must be accounted for when measuring C_f.

The small remaining offset (~0.2) between the C_f(H I) and C_f (metal) is comparable to the C_f (Si II)–C_f(H I) relation found in Eq. (12). To demonstrate this, we convert our C_f (H I) and C_f(Si II) values, measured in a uniform geometry, to a clumpy geometry, and recalculate how C_f (Si II) scales with C_f(H I) (green line in Fig. 6). This recalculated relation is within the errors of the Reddy et al. (2016b) Si II 1260 Å relation (blue dashed line). The consistent scaling between C_f (H I) and C_f(Si II) from our study and Reddy et al. (2016b) demonstrates first that C_f (Si II)≠C_f(H I), and second that empirical relations can convert C_f(Si II) into C_f(H I) for a variety of galaxy types. Both findings are crucial for indirect measurements of the Lyman continuum escape fraction, as discussed in Paper II.

There is still a C_f spread among the different transitions: strong Si II 1260 Å and C II lines have similar C_f, but weaker Si II 1527 and Al III lines have smaller C_f. Puzzlingly, the Si II 1260 Å and Si II 1527 Å lines arise from the same transition, and should have the same C_f. Different C_f (Si II 1260 Å) and C_f (Si II 1527 Å) would be expected if Si II 1527 Å was not saturated, but Reddy et al. (2016b) showed that the Si II 1527 Å/Si II 1260 Å equivalent width ratio is unity. Consequently, both transitions are saturated. Alternatively, since Si II 1260 Å is a stronger line than Si II 1527 Å (f value of 1.22 versus 0.133), there could be more optically thick Si II 1260 Å regions along the line of sight than optically thick Si II 1527 Å regions. Intersecting more regions along the line of sight covers more of the background continuum, such that C_f (Si II 1260 Å) is larger than C_f(Si II 1527 Å). This differential covering fraction is similar to the proposed physical mechanism creating the difference between C_f (Si II) and C_f (H I) in Sect. 5.3. In the Si II 1260 Å and Si II 1527 Å case, the number of absorbers is different owing to differences in line strength, not metallicity.

Fig. 6 Metal covering fractions from Si II, C II, and Al III (as given by the line types in the legend) as a function of the H I covering fraction (Cf(H I)) from Reddy et al. (2016b). Cf(H I) is always measured in a clumpy geometry. The black lines show 1–R measured from the dereddened stellar continuum normalized stacks. This assumes a clumpy geometry with a uniform dust screen. The blue lines show when the Cf (metal) is derived in the same geometry as Cf(H I). When calculated with consistent geometries, Cf(Si II 1260 Å) (blue dashedline) is offset from Cf(H I) by a similaramount as our observations in Fig. 4. The green dashed line shows our fitted relationship between Cf(Si II) and Cf (H I) from Fig. 4, after converting it to clumpy geometry. The Reddy et al. (2016b) Cf (H I)–Cf(Si II) relationship is consistent, within the errors, with ours. The red dashed line shows a one-to-one relation.

6 Conclusions

We have analyzed the H I, O I, and Si II low ionization interstellar absorption lines from a sample of 18 star-forming galaxies with rest-frame ultraviolet spectroscopy of the Lyman series. The majority of the sources have COS spectra taken with the HST and are at low redshift (z < 0.4). Our sample includes nine Lyman continuum leaking galaxies. We fit the stellar continuum, dust attenuation, H I Lyman series absorption lines, and several low ionization absorption lines (Si II λλ1190, 1193 and λ1260 Å, O I λ 1039 Å in particular). These fits determine the UV attenuation, as well as column densities and covering fractions of neutral hydrogen and metals (Sect. 3). Additionally, we applied our fitting method to synthetic Lyman series ISM absorption lines to investigate the systematic errors of the covering fractions and column densities. Then, we studied the observed Lyman series lines to constrain the H I properties. The direct H I properties were compared to indirect estimates of the neutral gas properties using O I and Si II absorption lines.

Our main results are summarized as follows:

• The H I covering fraction is accurately recovered from the Lyman series. Synthetic spectra recover the covering fractions with low systematics for a wide range of S/N and resolutions (S∕N ≥ 2 and R > 3000, or S∕N ≥ 5 and R ≥ 600; Fig. 2). Future observatories, such as JWST, will accurately measure covering fractions of high-redshift galaxies.

• The observed H I lines are found to be saturated in all galaxies. Assuming the same O/H abundance in the neutral and ionized gas, we derive H I column densities of log(N_H I) ~ 18.6−20 cm⁻² from theO I absorption lines (Sect. 3.3.2).

• The H I column densities derived for the known LyC leakers are too high to allow ionizing photons to escape. Rather, we find that the LyC emitting galaxies have H I covering fractions below unity. Ionizing photons escape through opticallythin holes in a clumpy ISM. The median H I covering fraction of confirmed LyC emitting galaxies is 0.62, as compared to 0.95 for galaxies that do not have LyC detections.

• The Si II covering fraction is systematically lower than the H I covering fraction (Fig. 4). However, the Si II covering fraction is found to scale linearly with the H I covering fraction (Eq. (12)). We show that this relation is compatible with the relationship from Reddy et al. (2016b) of stacked z ~ 3 spectra. Thus, with an empirical correction, the Si II absorption lines can be used to determine the H I covering fraction when the Lyman series is not observable. This is especially powerful at high redshift when the Lyman series cannot be observed.

• The assumed dust geometry (here a uniform screen) impacts the measured covering fractions and dust attenuations (Fig. 5), but it does not impact the inferred escape fractions of ionizing photons (Sect. 5.3). Crucially, the geometric covering fraction cannot be read off from the residual line flux in all geometries. A consistent fitting of the dust attenuation, continuum, and absorption lines is required to determine the covering fraction properly.

By relating the H I covering fraction to metal absorption lines, we have provided the framework to measure the neutral gas covering fractions in star-forming galaxies. This will be of particular interest at high redshift, where the Lyman series and Lyman continuum are not observable and indirect measurements are required to measure the escape fraction of ionizing photons. In a companion paper (Chisholm et al. 2018), we use the absorption lines to accurately predict the escape fraction of the known LyC leakers. These methods also yield consistent results for a sample of z ~ 2 − 2.4 galaxies. Our analysis emphasizes that UV spectra of sufficient quality (S∕N > 5) from JWST and extremely large telescopes may constrain how the Universe became reionized.

Acknowledgements

We would like to thank the anonymous referee for carefully reading our manuscript and for giving constructive comments which improved the quality of the paper.

Appendix A Simulated spectra

Here we detail the procedure we used to create the simulated spectra in Sect. 3.2. We simulated synthetic spectra to mimic the observed Lyman series. Therefore, we created synthetic spectra using a single STARBURST99 stellar continuum model of R(S99) ≈ 2500 (Leitherer et al. 1999, 2010) with a stellar age of 3 Myr, a solar metallicity, and without dust attenuation. This stellar age has a strong O IV P-Cygni profile (broad blueshifted absorption and redshifted emission) that blends with the Lyβ absorption line (see Fig. 1). The synthetic spectra include the wavelength regime between 910–1050 Å with a pixel separation of 0.02 Å (6 km s⁻¹ at 1035 Å). These wavelengths include the Lyman series, from Lyβ to the Lyman break (≈911.8 Å). We imprinted on top of the stellar continuum ISM absorption lines from H I, O I, O VI, and C II (see 2) with Voigt profiles. The Doppler broadening parameters (b), column densities (N), and covering fractions (C_f) are created with independent values for each element (see Table A.1). We fixed the velocities of each line at 0 km s⁻¹. The Einstein-A coefficients and the oscillator strengths of each transition were taken from the NIST Atomic Database².

We defined two categories of synthetic spectra corresponding to an ionizing escape fraction of 10%. This value is drawn from the largest known publicly available low-redshift LyC escape fractions (Izotov et al. 2016a,b). For the first category, we assumed holes in an otherwise optically thick ISM allow ionizing photons to escape. We called this the picket-fence model. For the picket-fence model, the H I parameters are N_H I = 10²⁰ cm⁻² (e^−τ_H → 0) and C_f (H I) = 0.90. For the second scenario, we assumed a low H I column density allows ionizing photons to escape (the density-bounded regime). We fixed N_H I = 10^17.57 cm⁻² (e^−τ_H ≈ 0.10) and C_f (H I) = 1.0. Since we fixed E_B-V = 0 in all our simulations, the geometry does not change our simulations (Sect. 3.1.1).

To reproduce the observations, we created synthetic spectra at various instrumental resolutions and S/Ns. We simulated seven resolutions, R = λ∕δλ (where δλ is the FWHM of the spectra), of 15 000, 6000, 3000, 1500, 1000, 750, and 600. For the synthetic spectra having R > R(S99), this is performed by combining the theoretical stellar continuum with R(S99) = 2500 withabsorption lines of spectral resolution R. The range of resolutions chosen correspond to spectrographs such as the COS on HST, the Low and High Resolution Imaging Spectrometer (LRIS/HRIS) on Keck, MagE on Magellan, MUSE, or the upcoming NIRSpec on JWST. For each R, we generated seven S/Ns (2, 5, 10, 20, 30, 40, and 50) in a sufficiently large range to study typically delivered noise levels. At each S/N, we generated 50 different sets of random Gaussian noise and added the noise to the same synthetic spectra. This produces a sample of 50 synthetic spectra where only the random noise varies between the spectra. Figures B.1 and B.2 show simulated spectra with resolutions of 15 000 and 600 for three S/Ns (noise-free, 10 and 2). These results are discussed in Sect. 4.

Appendix B Fits

We provide fits of the synthetic and observed spectra detailed in this article. Fig. B.1 and Fig. B.2 show the synthetic spectra for three S/Ns (∞, 10 and 2) with R = 15 000 and R = 600, respectively. In Figs. B.3–B.20 we present the fits for the 18 galaxies in our sample. The top panels of each figure show the full wavelength coverage; the observed flux is indicated in black and the total fit in red. Gray flux indicates that these regions were masked during the fit. The green lines show the error on the observed flux. The lower panels zoom in on individual Lyman series lines. The upper portions of each plot indicate the fitted Lyman series lines (blue), ISM metal absorption lines (solid black line), and Milky Way absorption lines (black dashed lines). The line labels are gray if they are not fit. Geocoronal regions are denoted by shaded gray regions.

Fig. B.1 Fits(in red) of the synthetic spectra at high resolution (R = 15 000) forthe picket-fence model (Cf(H I) = 0.9, NH I = 1020cm−2) with S/N equal to ∞ (top panel), 10 (middle panel), and 2 (bottom panel). We exclude the gray portion when fitting. In the upper portion of the top panel, we indicate and label the fitted absorption lines.

Fig. B.2 Same as Fig. B.1 but for a spectral resolution of R = 600.

Fig. B.3 Fit (in red solid line) of the COS G130M spectrum of the galaxy J0921+4509. Black is the observed flux included in the routine either to fit the stellar continuum or the ISM absorption lines. Gray portions are masked out for both steps; gray shaded regions indicate those masked because of geocoronal emission. The flux error array appears in green. We display the ISM and Milky Way absorption lines as solid and dotted lines in the upper portion of each panel. Black or blue labels indicate that the lines are fit, whereas gray labels indicate that they are not fit. When present, red labels indicate lines that are not detected. Reference for this observation: Borthakur et al. (2014).

Fig. B.4 Same as Fig. B.3 but for the COS G140L spectrum of J1503+3644. Reference for this observation: Izotov et al. (2016b).

Fig. B.5 Same as Fig. B.3 but for the COS G140L spectrum of J0925+1409. Reference for this observation: Izotov et al. (2016a).

Fig. B.6 Same as Fig. B.3 but for the COS G140L spectrum of J1152+3400. Reference for this observation: Izotov et al. (2016b).

Fig. B.7 Same as Fig. B.3 but for the COS G140L spectrum of J1333+6246. Reference for this observation: Izotov et al. (2016b).

Fig. B.8 Same as Fig. B.3 but for the COS G140L spectrum of J1442–0209. Reference for this observation: Izotov et al. (2016b).

Fig. B.9 Same as Fig. B.3 but for the COS G140L spectrum of Tol1247–232. Reference for this observation: Leitherer et al. (2016).

Fig. B.10 Best fit for the stellar continuum for the COS G140L spectrum of Tol0440–381. This galaxy has a low redshift and the Milky Way absorption lines contaminate the Lyman series. Consequently, we do not fit for the ISM absorption lines. Reference for this observation: Leitherer et al. (2016).

Fig. B.11 Same as Fig. B.10 but for the COS G140L spectrum of Mrk54. Reference for this observation: Leitherer et al. (2016).

Fig. B.12 Same as Fig. B.3 but for the COS G130M spectrum of J0926+4427. Reference for this observation: Heckman et al. (2011).

Fig. B.13 Same as Fig. B.3 but for the COS G130M spectrum of J1429+0643. Reference for this observation: Heckman et al. (2011).

Fig. B.14 Same as Fig. B.3 but for the COS G130M spectrum of GP0303–0759. Reference for this observation: Henry et al. (2015).

Fig. B.15 Same as Fig. B.3 but for the COS G130M spectrum of GP1244+0216. Reference for this observation: Henry et al. (2015).

Fig. B.16 Same as Fig. B.3 but for the COS G130M spectrum of GP1054+5238. Reference for this observation: Henry et al. (2015).

Fig. B.17 Same as Fig. B.3 but for the COS G130M spectrum of GP0911+1831. Reference for this observation: Henry et al. (2015).

Fig. B.18 Same as Fig. B.3 but for the MagE spectrum of SGAS J122651.3+215220. Reference for this observation: Rigby et al. (2018).

Fig. B.19 Same as Fig. B.3 but for the MagE spectrum of SGAS J152745.1+065219. Reference for this observation: Rigby et al. (2018).

Fig. B.20 Same as Fig. B.3 but for the MagE spectrum of the Cosmic Eye. Reference for this observation: Rigby et al. (2018).

Appendix C Additional tables

C.1 Measured systematic errors from single synthetic spectra

We tabulate the systematic errors measured from the synthetic spectra (Sect. 4). Table C.1 gives the N_H I systematic errors, and Table C.2 gives the C_f systematic errors. The percent errors plotted in Fig. 2 are computed by subtracting 0.9 from the errors and dividing the difference by 0.9 (the C_f used to generate the spectra).

C.2 Lyman series residual flux

Table C.3 lists the residual flux of the individual Lyman series lines (Sect. 3.3.1).

References

All Tables

All Figures

Fig. 1 Total fit (red solid line) of the Lyβ absorption line for the galaxy GP 1244+0216 (black data). The contributions from the stellar continuum (dashed red line), Lyβ absorption line (blue solid line), and O I absorption line (orange solid line) blend together. O I absorption lines at 989 and 1039 Å are unblended and robustly constrain the O I profile. All three components need to be accounted for when fitting the Lyman series.

In the text

Fig. 2 Color map of the systematic percent error of the covering fraction as a function of the resolution and S/N. The synthetic spectra are created with R < 120 km s−1 combine a theoretical stellar continuum spectra with R(S99) = 120 km s−1 and absorption lines of spectral resolution R. The covering fraction is recovered to within 5% of the estimated parameter for all observations within the dark blue region.

In the text

	Fig. 3 Histogram of the H I covering fraction (Cf(H I)) from the 18 galaxies in our sample. We split the sample into galaxies known to leak LyC photons (blue) and galaxies without measured LyC emission (black). The leakers have the lowest Cf (H I) and the unknown emitters have the highest Cf(H I).
In the text

	Fig. 4 Comparison of the Si II and H I covering fractions. The black solid line shows a one-to-one relation and the dotted line is the fitted relation (Eq. (12)). The average offset between Cf (H I) and Cf (Si II) is 0.25.
In the text

Fig. 5 Residual flux R of a saturated absorption line from Eq. (4) for various values of Aλ = kλ EB-V in a clumpy ISM geometry. When dust is present, the residual flux is higher than (1 − Cf). However, in a uniform dust screen geometry, the residual flux is always R = (1 − Cf). The covering fraction derived from the observed residual flux and a clumpy dust geometry is underestimated when dust is present.

In the text

Fig. 6 Metal covering fractions from Si II, C II, and Al III (as given by the line types in the legend) as a function of the H I covering fraction (Cf(H I)) from Reddy et al. (2016b). Cf(H I) is always measured in a clumpy geometry. The black lines show 1–R measured from the dereddened stellar continuum normalized stacks. This assumes a clumpy geometry with a uniform dust screen. The blue lines show when the Cf (metal) is derived in the same geometry as Cf(H I). When calculated with consistent geometries, Cf(Si II 1260 Å) (blue dashedline) is offset from Cf(H I) by a similaramount as our observations in Fig. 4. The green dashed line shows our fitted relationship between Cf(Si II) and Cf (H I) from Fig. 4, after converting it to clumpy geometry. The Reddy et al. (2016b) Cf (H I)–Cf(Si II) relationship is consistent, within the errors, with ours. The red dashed line shows a one-to-one relation.

In the text

	Fig. B.1 Fits(in red) of the synthetic spectra at high resolution (R = 15 000) forthe picket-fence model (Cf(H I) = 0.9, NH I = 1020cm−2) with S/N equal to ∞ (top panel), 10 (middle panel), and 2 (bottom panel). We exclude the gray portion when fitting. In the upper portion of the top panel, we indicate and label the fitted absorption lines.
In the text

	Fig. B.2 Same as Fig. B.1 but for a spectral resolution of R = 600.
In the text

Fig. B.3 Fit (in red solid line) of the COS G130M spectrum of the galaxy J0921+4509. Black is the observed flux included in the routine either to fit the stellar continuum or the ISM absorption lines. Gray portions are masked out for both steps; gray shaded regions indicate those masked because of geocoronal emission. The flux error array appears in green. We display the ISM and Milky Way absorption lines as solid and dotted lines in the upper portion of each panel. Black or blue labels indicate that the lines are fit, whereas gray labels indicate that they are not fit. When present, red labels indicate lines that are not detected. Reference for this observation: Borthakur et al. (2014).

In the text

	Fig. B.4 Same as Fig. B.3 but for the COS G140L spectrum of J1503+3644. Reference for this observation: Izotov et al. (2016b).
In the text

	Fig. B.5 Same as Fig. B.3 but for the COS G140L spectrum of J0925+1409. Reference for this observation: Izotov et al. (2016a).
In the text

	Fig. B.6 Same as Fig. B.3 but for the COS G140L spectrum of J1152+3400. Reference for this observation: Izotov et al. (2016b).
In the text

	Fig. B.7 Same as Fig. B.3 but for the COS G140L spectrum of J1333+6246. Reference for this observation: Izotov et al. (2016b).
In the text

	Fig. B.8 Same as Fig. B.3 but for the COS G140L spectrum of J1442–0209. Reference for this observation: Izotov et al. (2016b).
In the text

	Fig. B.9 Same as Fig. B.3 but for the COS G140L spectrum of Tol1247–232. Reference for this observation: Leitherer et al. (2016).
In the text

	Fig. B.10 Best fit for the stellar continuum for the COS G140L spectrum of Tol0440–381. This galaxy has a low redshift and the Milky Way absorption lines contaminate the Lyman series. Consequently, we do not fit for the ISM absorption lines. Reference for this observation: Leitherer et al. (2016).
In the text

	Fig. B.11 Same as Fig. B.10 but for the COS G140L spectrum of Mrk54. Reference for this observation: Leitherer et al. (2016).
In the text

	Fig. B.12 Same as Fig. B.3 but for the COS G130M spectrum of J0926+4427. Reference for this observation: Heckman et al. (2011).
In the text

	Fig. B.13 Same as Fig. B.3 but for the COS G130M spectrum of J1429+0643. Reference for this observation: Heckman et al. (2011).
In the text

	Fig. B.14 Same as Fig. B.3 but for the COS G130M spectrum of GP0303–0759. Reference for this observation: Henry et al. (2015).
In the text

	Fig. B.15 Same as Fig. B.3 but for the COS G130M spectrum of GP1244+0216. Reference for this observation: Henry et al. (2015).
In the text

	Fig. B.16 Same as Fig. B.3 but for the COS G130M spectrum of GP1054+5238. Reference for this observation: Henry et al. (2015).
In the text

	Fig. B.17 Same as Fig. B.3 but for the COS G130M spectrum of GP0911+1831. Reference for this observation: Henry et al. (2015).
In the text

	Fig. B.18 Same as Fig. B.3 but for the MagE spectrum of SGAS J122651.3+215220. Reference for this observation: Rigby et al. (2018).
In the text

	Fig. B.19 Same as Fig. B.3 but for the MagE spectrum of SGAS J152745.1+065219. Reference for this observation: Rigby et al. (2018).
In the text

	Fig. B.20 Same as Fig. B.3 but for the MagE spectrum of the Cosmic Eye. Reference for this observation: Rigby et al. (2018).
In the text