© ESO 2022

1. Introduction

Since the discovery of the Gunn–Peterson effect (Gunn & Peterson 1965; Becker et al. 2001), understanding how the Universe was reionized around z ≈ 6 − 9 (Planck Collaboration XLVII 2016) requires not only identification of the sources that contributed to that change, but also knowledge of the physical mechanisms under the last major cosmic phase transition. There are two principal instigators of cosmic reionization. The first is massive stars that emit sufficient ionizing photons to keep hydrogen in an ionized state (Robertson et al. 2015), and the second is active galactic nuclei (AGN, Madau & Haardt 2015). Depending on the principal source for reionization (galaxies or AGN, or both; see Dayal et al. 2020), the thermal history of the late Universe changes accordingly (e.g., Kulkarni et al. 2017). For example, models primarily driven by star-forming (SF) galaxies generally have delayed reionization compared to AGN-dominated reionization, where the Universe is reionized at a higher redshift (e.g., see Finkelstein et al. 2019).

To establish the contribution of galaxies and AGN in terms of ionization rate, one needs to quantify three main properties: (1) the ionizing photon production efficiency (ξ_ion), (2) the ionizing photon escape fraction ($f_{\mathrm{esc}}^{\mathrm{abs}}$), and the number density of such galaxies or AGN at high z. The ionizing production efficiency (ξ_ion) is the intrinsic number of ionizing photons emitted by the stars per unit UV luminosity, while $f_{\mathrm{esc}}^{\mathrm{abs}}$ is the fraction of these photons escaping from them to the intergalactic medium (IGM).

Although always powerful in terms of ionizing emissivity ($f_{\text{esc}}^{\text{abs}}\sim 1$, Stevans et al. 2014; Lusso et al. 2015), AGNs appear to be very rare at high redshifts (Fontanot et al. 2014; Matsuoka et al. 2018) and with lower accretion rates and emissivities than in the late Universe (Dayal et al. 2020); although there is still the possibility of having an extending low-luminosity AGN population that is numerous enough to provide a substantial contribution to the total ionizing budget (Giallongo et al. 2015; Grazian et al. 2020). Star-forming galaxies are much more numerous at high-z, and the question now becomes whether more massive or compact low-mass SF galaxies contribute the most. On the one hand, the less numerous but massive SF galaxies with enormous gas and dust reservoirs at high-z can intrinsically generate a large amount of Lyman continuum (LyC) photons (Naidu et al. 2020), but these photons can be absorbed and attenuated within the gaseous reservoirs of their host galaxies at the same time. On the other hand, compact low-mass SF galaxies are more numerous at high z (Rosdahl et al. 2018; Finkelstein et al. 2019). Additionally, they take advantage of their weak gravitational potential and feedback effects to clean the interstellar medium (ISM) from gas and dust, increasing the fraction of LyC photons escaping the galaxy (Trebitsch et al. 2017). The true contribution to reionization from different galaxy types is still under debate, but models require a threshold of $f_{\text{esc}}^{\text{abs}}=5-20\%$ for M_UV ≲ −15 galaxies to drive reionization around z ∼ 7 (Robertson et al. 2013, 2015).

Over the last decade, several campaigns have searched for LyC emitters at low and intermediate redshifts. Using the Far Ultraviolet Spectroscopic Explorer (FUSE), Bergvall et al. (2006) and Leitet et al. (2013) reported the first LyC detections, although later re-observations revealed that some of them were affected by foreground scattered light. Using the Hubble Space Telescope (HST), about 30 new leakers were later discovered at z ∼ 0.3 (Leitherer et al. 2016; Borthakur et al. 2016; Izotov et al. 2016a,b, 2018a,b, 2021; Puschnig et al. 2017; Wang et al. 2019; Malkan & Malkan 2021), with $f_{\mathrm{esc}}^{\mathrm{abs}}$ values spanning from a few to 70%. The majority consist of compact SF galaxies with relatively high star formation rates (SFR), strong optical emission lines, and ionizing photon production efficiencies comparable to their high-z analogs (Schaerer et al. 2016). Nevertheless, it remains challenging to robustly detect faint LyC flux at the sensitivity limit of HST (see e.g., the re-analysis of Chisholm et al. 2017 versus Leitherer et al. 2016 and Puschnig et al. 2017).

At intermediate redshifts, investigating the LyC leakage is more difficult, as LyC photons are absorbed by Lyman limit systems in the IGM, which makes the detection extremely challenging beyond z ≥ 4. The presence of low-redshift interlopers also renders potential LyC detections at high z difficult (Vanzella et al. 2010; Inoue et al. 2014; Naidu et al. 2018). Rivera-Thorsen et al. (2019) recently reported LyC emission from the Sunburst Arc, a lensed compact dwarf galaxy at redshift z ∼ 2.2, with an escape fraction of $f_{\text{esc}}^{\text{abs}}=32\%$ averaged over 11 different gravitationally lensed HST images. Similar works are Bian et al. (2017) and the Cosmic Horseshoe in Vasei et al. (2016). The Keck Lyman Continuum Survey (KLCS, Steidel et al. 2018, see also Shapley et al. 2016) spectroscopically confirmed 13 individual LyC detections around z ∼ 3.1 in 9 independent fields (see the follow-up HST observations in Mostardi et al. 2015; Pahl et al. 2021). Using HST/F336W and narrow-band photometry, another 12 Lyman break galaxies (LBGs) and AGNs – enclosed in the SSA22 proto-cluster field (cf., Shapley et al. 2006; Micheva et al. 2017a,b) – were discovered to be leaking in the LymAn Continuum Escape Survey (LACES) by Fletcher et al. (2019). In Vanzella et al. (2015, 2016) and de Barros et al. (2016), the authors used the CANDELS U-band and HST/F336W imaging to select LyC leaker candidates in the GOODS-S field, resulting in the robust detection of the Ion2 galaxy at z = 3.2. With a similar methodology, Ion1 was discovered by Vanzella et al. (2012) and Ion3 found by Vanzella et al. (2018), which currently holds the record as the most distant SF galaxy observed to leak LyC photons (z ∼ 4). The recent re-analysis of Ion1 by Ji et al. (2020) suggests a 5%–10% escape, although this was originally measured at 60% (Vanzella et al. 2012). At z ∼ 3, three individual detections were performed by Grazian et al. (2016, 2017) in the CANDELS/ GOODS-N, EGS, and COSMOS fields thanks to deep ground-based observations with the Large Binocular Telescope (LBT). Finally, the recently launched AstroSAT telescope discovered a remarkable z = 1.42 SF galaxy that emits ionizing photons at a rest-frame of 650 Å (Saha et al. 2020). Some other tentative detections have been conducted by Siana et al. (2015) and Smith et al. (2020).

So far, low-z galaxy studies have shown that individual SF galaxies can have $f_{\mathrm{esc}}^{\mathrm{abs}}$ larger than the 10%–20% threshold for galaxies to reionize the Universe. However, at high z, the actual number of individual LyC detections is not only low but biased to the sightlines with the highest IGM transmissivity. As the number density and the properties of such high-z leakers are not yet known (and we cannot extrapolate them from the few clear detections), stacking analyses have become the most popular method in reionization studies, and are used to obtain the average $f_{\mathrm{esc}}^{\mathrm{abs}}$ of the SF galaxy population. Deep broad-band (Grazian et al. 2016, 2017; Rutkowski et al. 2017; Japelj et al. 2017; Naidu et al. 2018), narrow-band imaging (Guaita et al. 2016), and spectral stacking (Marchi et al. 2018; Meštrić et al. 2021) have yielded relatively low upper limits for the escape fraction of the SF population at z = 2 − 4 (< 5%–10% for M_UV ≲ −19 galaxies, cf., Naidu et al. 2018; Alavi et al. 2020).

In order to detect and quantify LyC emission at high z (z ≳ 4), including in the Epoch of Reionization, the development of indirect techniques that estimate the LyC escape are required. Low-z galaxies are the most suitable laboratories for testing indirect LyC indicators, because the LyC, the far-UV (FUV) stellar, and the interstellar absorption spectra, as well as the rest-frame optical nebular lines, can be observed at the same time and with better S/N than at high-z, allowing spatially resolved studies.

Starting from the closest to the Lyman limit (912 Å), the depth of the Lyman series lines directly links the neutral gas column density to the geometrical distribution of the same gas along the line of sight. The residual flux of the high-order Lyman absorption lines has been related with the neutral gas porosity (see Reddy et al. 2016b, 2022). Properly accounting for dust attenuation, these lines have proved to be an effective LyC tracer (Gazagnes et al. 2018, 2020; Chisholm et al. 2018; Steidel et al. 2018), although simulations show different results (Mauerhofer et al. 2021). Very prominent when present, the Lyαλ1216 line – accessible from the ground at z ≳ 2 – can also be used because Lyα and LyC photons interact with the same neutral gas. Jaskot & Oey (2014) and Verhamme et al. (2015) proposed that the Lyα peak separation can trace the presence of low-density gas: the physical mechanism driving the leakage is expected to be connected for both Lyα and LyC (Rivera-Thorsen et al. 2015; Verhamme et al. 2015; Dijkstra et al. 2016). Indeed, the separation and shape of double- and triple-peaked Lyα profiles have been associated with LyC leakage in observations (Rivera-Thorsen et al. 2017; Verhamme et al. 2017; Vanzella et al. 2020; Izotov et al. 2020, 2021). However, the Lyman series are still efficiently suppressed by IGM absorption at high-z and are therefore almost impossible to observe during the Epoch of Reionization. The use of indirect LyC indicators redwards Lyα is essential.

The most promising indirect tracers of LyC escape to date are the UV low-ionization state (LIS) absorption lines (Erb 2015; Chisholm et al. 2018) and the Mg IIλλ2796,2803 emission doublet (Henry et al. 2018; Chisholm et al. 2020). On the one hand, measurements of some UV metal absorption line depths such as Si IIλ1260 (Gazagnes et al. 2018) and C IIλ1334 (Grimes et al. 2009) have been shown to uniquely trace the neutral gas covering fraction (Reddy et al. 2016b, 2022). The observed depth of a given absorption transition depends on the optical depth and the covering fraction. If the optical depth is high enough that the lines are saturated, then the covering fraction can be inferred from the observed spectral line depths. On the other hand, the intrinsically bright Mg II nebular transitions have an ionization potential similar to that of warm H I. This makes the Mg II and LyC transitions optically thick at similar H I column densities (for typical Mg/H abundances), making Mg IIλλ2796,2803 a good tracer of LyC escape as well.

In the rest-frame optical range, some previous studies suggest that both low-z (Jaskot & Oey 2013; Izotov et al. 2018b) and high-z leakers (Nakajima et al. 2016) have high [O III]λλ4959,5007 versus [O II]λλ3726,3729 ratios (O₃₂, see also Faisst 2016). Subsequent efforts showed that this quantity is not a good predictor of the LyC escape (Bassett et al. 2019; Nakajima et al. 2020), mainly because O₃₂ alone cannot probe the outer extent of the ionization regions, the region that is most sensitive to the optical depth and thus chiefly affects the escape of LyC radiation. Different kinematic effects can also impact the LyC leakage, complicating the use of emission lines in the diagnostics (McKinney et al. 2019; Hogarth et al. 2020). Worth mentioning is also the [S II]λλ6717,6731 doublet (Wang et al. 2019; Ramambason et al. 2020), which shows a deficit in leakers and was recently investigated more thoroughly by Wang et al. (2021).

Finally, other indicators have been investigated in the literature, such as the SFR surface density (Heckman et al. 2001; Clarke & Oey 2002; Izotov et al. 2018b, 2021; Naidu et al. 2020; Cen 2020), and UV/optical lines such as He Iλ3889 (Izotov et al. 2017; Guseva et al. 2020) or O Iλ6300 (Ramambason et al. 2020), but these require further calibrations to be used as robust predictors of LyC escape. A detailed analysis of the various indirect LyC diagnostics is presented in Flury et al. (2022b).

All these indicators are sensitive to the ISM properties, and we can therefore expect correlations between different ISM-related quantities and LyC escape. Revealing the ISM properties of leakers versus nonleakers, one could illustrate the limitations and strengths of the previous tracers. In summary, two scenarios describing the ISM structure compete in the escape of LyC radiation (Zackrisson et al. 2013). The first, called density-bounded regions (Jaskot & Oey 2013; Nakajima & Ouchi 2014), envision a massive stellar population surrounded by a homogeneous and low-density neutral gas media, so that photons can leak from the galaxy at a rate that is inversely proportional to the hydrogen column density. Alternatively, in the “picket-fence” model (Heckman et al. 2001, 2011; Vasei et al. 2016; Reddy et al. 2016b, 2022; Gazagnes et al. 2018; Steidel et al. 2018), stars are surrounded by an optically thick ISM and the ionizing radiation can only escape through low-column-density channels: the ISM is not homogeneous but patchy. So far, local observations of leakers (Gazagnes et al. 2018, 2020) and also simulations (Trebitsch et al. 2017; Mauerhofer et al. 2021) favor a patchy ISM, and only the galaxies with the highest $f_{\mathrm{esc}}^{\mathrm{abs}}$ have a mostly ionized ISM consistent with a density-bounded scenario (Jaskot et al. 2019; Ramambason et al. 2020; Kakiichi & Gronke 2021). On top of the geometric distribution, other aspects like photo-ionized outflows and feedback effects (Chisholm et al. 2017; Hogarth et al. 2020; Carr et al. 2021) also change the ISM structure, affecting the profile of emission/absorption lines and the escape of LyC radiation.

In this paper, we statistically explore the impact of the neutral and low-ionized ISM absorption on LyC escape, taking a novel approach and using a large set of low-z LyC galaxy observations. We use spectra from the HST Low-z Lyman Continuum Survey (LzLCS, see Sect. 2), a set of 66 SF galaxies of which 35 have a secure LyC detection. The methodology for modelling the stellar continuum of the galaxies in the sample, the adopted absolute LyC photon escape fraction, and the absorption line measurements are presented in Sect. 3. The absorption line results are shown in the same section, and the possible inclusion of systematic errors is studied through absorption line simulations. In Sect. 4, we discuss the fitting results, secondary products, and limitations of the modelling. The main empirical relations between H I and LIS lines and different galaxy properties are then given in Sects. 5 and 6. A physical interpretation of the results and literature comparisons are presented in Sect. 7. We finally conclude in Sect. 8, also highlighting the relevance of these results for future high-z searches.

Throughout this paper, a standard ΛCDM cosmology is used, with a matter density parameter Ω_M = 0.3, a vacuum energy density parameter Ω_Λ = 0.7, and Hubble constant of H₀ = 70 km s⁻¹ Mpc⁻¹. All magnitudes are in the AB system, and we adopt a solar metallicity value of 12 + log(O/H)_⊙ = 8.69 (Z_⊙ = 1).

2. Data

In this study, we use UV spectra in the observed frame wavelength range 800–1950 Å from the Low-Redshift Lyman Continuum Survey (LzLCS, Flury et al. 2022a), the most comprehensive spectroscopic campaign so far to trace the LyC emission of galaxies in the nearby Universe (z ∼ 0.3). In addition, we include and reanalyze previous LyC emitters from the studies of Izotov et al. (2016a,b, 2018a,b, 2021) and Wang et al. (2019) as a comparison sample. For all the sources, we also use ancillary data derived in a homogeneous fashion by the LzLCS team (see Sect. 2.3).

2.1. The Low-z Lyman Continuum Survey

The Low-z Lyman Continuum Survey (LzLCS, see Flury et al. 2022a) is composed of 66 SF galaxies at z ∼ 0.3 whose UV spectra were observed by the Cosmic Origin Spectrograph (COS) on board the HST. At z ∼ 0.3, the COS/G140L grating not only probes the LyC but also traces an entire window of wavelengths redder than Lyα up to 1500 Å rest-frame.

The parent sample of LyC-emitter candidates was drawn from Sloan Digital Sky Survey sources (SDSS Data Release 15, Blanton et al. 2017) with GALEX counterparts in both the FUV (0.1530 μm) and near-UV (NUV; 0.2310 μm) photometric bands (Morrissey et al. 2007). Applying the methods described in Jaskot et al. (2019), the equivalent widths and fluxes of different optical emission lines were measured in the SDSS spectra, and the SF galaxies in the redshift range 0.22 ≤ z < 0.45 were selected by applying usual BPT (Baldwin et al. 1981) diagram diagnostics (HST/COS traces the LyC region at z > 0.22). Star-formation rate surface densities (Σ_SFR) were measured from the Hα and Hβ fluxes in the SDSS spectra and the SDSS u-band half-light radius, and UV β-indices for this subsample were computed using GALEX photometry. These measurements were later improved with the COS data in hand.

The final sample was built by selecting the galaxies that fulfill one or more of the following criteria (see Izotov et al. 2016a): high [O III]5007/[O II]3727 flux ratios (O₃₂ > 3), high SFR densities (Σ_SFR > 0.1 M_⊙ yr⁻¹ kpc⁻²), and blue UV continuum slopes (β < −2), resulting in 66 galaxies. Spectra of each galaxy were taken under the GO 15626 HST large program (P.I. Jaskot), with around two orbits per object on average. COS acquired each object via NUV imaging and centered its 2.5″ diameter spectroscopic aperture on the peak of the NUV emission. The G140L grating was then used to cover the 800–1950 Å observed wavelength range, with a nominal full width at half maximum (FWHM) resolution of R ≡ λ/Δλ ∼ 1000 at 1100 Å, increasing up to R ∼ 1500 at 1900 Å.

Science spectra were reduced following the methods described in Flury et al. (2022a) using the standard CALCOS pipeline (v3.3.9). The tackling of the systematic errors and the LyC signal-to-noise ratios (S/Ns) was improved thanks to the use of the FAINTCOS custom software (Makan et al. 2021), which properly estimates the dark current and scattered geo-coronal Lyα, co-adding spectra from multiple visits to improve the S/N. Finally, each spectrum was binned to 0.5621 Å to yield a sampling of ≳2 pixels per resolution element. The spectra were corrected for Milky Way extinction using the Galactic E_B−V estimates from the Green (2018) dust maps and the Fitzpatrick (1999) attenuation law.

The LyC flux (f^LyC) was measured as the mean over a 20 Å window as close as possible to but shortward of 900 Å in the rest-frame (typically 850–900 Å), while avoiding wavelengths redward of the observed frame of 1180 Å to reduce the impact of geo-coronal Lyα emission. Over the 66 galaxies in the LzLCS, 35 were significantly detected (LyC flux S/N ≥ 2) and 31 were not detected at LyC wavelengths. In the detected galaxies, LyC luminosities span from 10³⁸ to 10⁴⁰ erg s⁻¹. In the undetected sample, the average upper uncertainty on the LyC flux density was reported as the 1σ upper limit. See Flury et al. (2022a) for details.

2.2. Low-z LyC emitters in the literature

In addition to the new LzLCS sources, we add 20 SF galaxies with LyC observations from the sample of Izotov et al. (2016a,b, 2018a,b, 2021). Some of these sources are among the strongest leakers in the literature at z ∼ 0.3 (i.e., with the highest ionizing escape fractions).

Originally, the Izotov et al. (2016a,b, 2018a,b) galaxies were selected from SDSS in a similar way to some of the LzLCS in terms of compactness, FUV brightness, and ionizing features in the optical spectra (EW(Hβ) > 100 Å and O₃₂ ≥ 5). The selected galaxies exhibit narrow (≲130 km s⁻¹) double-peaked Lyα emission and high (up to 60%) Lyα escape fractions in most cases, low metallicities (1/10–1/3 Z_⊙), low stellar masses (0.15 − 6 × 10⁹ M_⊙), and high SFRs (14–80 M_⊙ yr⁻¹) and SFR densities (2–35 M_⊙ yr⁻¹ kpc⁻²). The resulting LyC escape in those galaxies extends from a few up to ∼70%. Here, we used these Izotov et al. SF galaxies as a reference sample for the strongest leakers in the new LzLCS.

Similarly, the Izotov et al. (2021) sample consists of 9 dust-poor and low-mass SF galaxies (≲10⁸ M_⊙), with very high specific SFRs (∼150–630 Gyr⁻¹), showing a wide range in Lyα peak velocity separations (∼200 − 500 km s⁻¹). Their published LyC escape fractions reach ∼35%.

In order to increase the dynamic range for of some the galaxy properties in our working sample, we finally also included three more SF galaxies observed by Wang et al. (2019), which were selected in terms of their [S II] deficiency (SDSS spectra), compactness ($\le 0\stackrel{″}{.}5$ in the SDSS u-band), and relatively high GALEX FUV fluxes. These galaxies are more massive (∼3 × 10¹⁰ M_⊙), more metal-rich (0.9–1 Z_⊙), dustier, and less extreme in terms of ionization properties than the Izotov et al. (2016a,b, 2018a,b) sample, but also have exceptionally high SFR densities (35–700 M_⊙ yr⁻¹ kpc⁻²). Their estimated escape fractions are significantly lower (up to 3%).

The LzLCS team re-processed and re-reduced these 23 COS/G140L spectra with the same methodology as that described above in order to maintain the homogeneity and consistency of the results (see Flury et al. 2022a). This results in 15 additional detections and 8 nondetections out of a total of 23 additional galaxies. The same total sample is described and analyzed in depth in Flury et al. (2022a,b).

2.3. Other LzLCS science products

The LzLCS team measured galaxy properties of every source in the LzLCS and in the selected literature samples, either from the SDSS spectra and the GALEX photometry, or directly from the COS spectra. In particular, for this work we use:

• Optical dust attenuation parameters: E_B − V-Balmer.

• Gas-phase metallicities: 12 + log(O/H).

• Lyα equivalent widths: W_Lyα.

Optical attenuations were computed using an iterative Balmer decrement methodology from the measured SDSS Hβ to Hϵ fluxes (following Izotov et al. 1994, see also Flury & Moran 2020). The gas-phase metallicities (relative oxygen abundances) were obtained using a direct-method approach using the extinction-corrected hydrogen and oxygen optical lines. Finally, Lyα equivalent widths were measured by integrating the COS spectral flux close to the Lyα line (continuum linearly fitted within a 100 Å window). Errors on every quantity come from 10⁴ Monte-Carlo realizations of the corresponding spectra. See Flury et al. (2022a) for more details. The LzLCS team is also measuring UV morphological parameters from the COS acquisition images (Ji et al., in prep.).

3. Methods

Measuring the properties of the ISM features in the UV requires accurate determination of the stellar continuum. In this section, we describe our methodology for modelling the UV stellar continuum of extragalactic spectra (Sect. 3.1, following Chisholm et al. 2019) and absolute LyC escape fractions in detail. The expressions used to determine the equivalent widths, residual fluxes, and corresponding escape fractions from the lines are also listed (Sect. 3.2).

3.1. Stellar continuum modelling

The stellar continuum modelling is achieved by fitting every COS spectrum with a linear combination of multiple bursts of single-age and single-metallicity stellar population models. We use the fully theoretical STARBURST99 single-star models without stellar rotation (S99, Leitherer et al. 2011a, 2014) using the Geneva evolution models (Meynet et al. 1994), and computed with the WM-BASIC method (Pauldrach et al. 2001; Leitherer et al. 2010). The S99 models use a Kroupa initial mass function (Kroupa 2001) with a high-(low-)mass exponent of 2.3 (1.3), and a high-mass cutoff at 100 M_⊙. The spectral resolution of the models is R ≡ λ/Δλ ∼ 2500, and remains constant at ≤ 1500 Å rest-frame.

A nebular continuum contribution was added to every single population by processing the stellar continuum models through the CLOUDY V17.0 code (Ferland et al. 2017), assuming similar gas-phase and stellar metallicities, an ionization parameter of log(U) = − 2.5, and a volume hydrogen density of n_H = 100 cm⁻³. The nebular continuum is negligible for the wavelengths that are typically fitted in this paper, but it can impact longer wavelength data (≳ 1200 Å rest-frame).

We chose four different metallicities (0.05, 0.2, 0.4 and 1 Z_⊙) and seven ages for each metallicity (1, 2, 3, 4, 5, 8, and 10 Myr) as a representative set of 28 models for our low-z UV spectra. In detail, ages are chosen to densely sample the stellar ages where the stellar continuum features appreciable change (Chisholm et al. 2019).

The observed spectra are first placed into the rest-frame by dividing by (1 + z). Both the observed spectra and the models are then normalized by the median flux in the 1070–1100 Å rest-frame wavelength interval (free of stellar and ISM features), and all the fits are performed in the same wavelength range: 975–1345 Å (rest-frame). Finally, the models are convolved by a Gaussian kernel¹ to the COS/G140L nominal spectral resolution (R ∼ 1000).

We first mask out the nonstellar features, the spectral regions that are affected by host-galaxy and Milky Way ISM absorptions, and geo-coronal emissions. We also mask bad pixels with zero flux, as well as pixels with S/N < 1. We then fit the UV stellar continuum (F^⋆(λ)) with a linear combination of multiple S99 models ($F_{i,j}^{S99}$) plus the nebular continuum. Adopting a uniform foreground dust-attenuation model for the galaxies (see Sect. 3.2), this yields:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill F^{\star }(\lambda )& ={10}^{-0.4k_{\lambda }{E}_{\mathrm{B}-\mathrm{V}}}{\displaystyle \sum _{i,j}}{X}_i{F}_{i,j}^{S99},\hfill \\ \hfill i& \equiv 1,2,3,4,5,8,10\mathrm{Myr}\hfill \\ \hfill j& \equiv 0.05,0.2,0.4,1Z_{\odot }\hfill \\ \hfill \end{array}\end{array}$$(1)

where 10^−0.4A_λ is the UV attenuation term, A_λ = k_λE_B-V, and k_λ is given by the Reddy et al. (2016a) attenuation law (R16 hereafter), k_λ = 2.191 + 0.974/λ(μm). We also tested the influence of different attenuation laws on the fitting results, in particular using the Small Magellanic Cloud (SMC) law instead (Prevot et al. 1984; Bouchet et al. 1985). We discuss these results in the next Sect. 4. The light-fraction coefficients (X_i) determine the weight of the single-stellar population within the fit, and the best fit is chosen through a nonlinear χ² minimization algorithm with respect to the observed data (LMFIT² package, Newville et al. 2016). Errors are derived in a Monte-Carlo way, varying observed pixel fluxes by a Gaussian distribution whose mean is zero and standard deviation is the 1σ-error of the flux in the same pixel, and re-fitting the continuum over 500 iterations per spectrum (enough realizations to sample the posterior continuum so that it approaches “Gaussianity” on each pixel). The adopted 1σ-errors on the spectra are approximated as the squared-root average of the original asymmetric spectral uncertainties (Flury et al. 2022a).

Figure 1 shows one of the outputs of our fitting routine, the observed and fitted spectrum for J091703, with the stellar continuum in red overlaying the observed spectrum in black. The reduced-χ² of the fit is ${\chi }_{\nu }^2=0.93$, with a median of 1.03 for the whole sample. This galaxy is dominated by a young and low-metallicity population whose average light-weighted age is 4.81 Myr and light-weighted metallicity results in 0.39 Z_⊙. J091703 is moderately attenuated, with a UV dust-attenuation parameter of E_B-V = 0.082 mag. We highlight with dark-blue labels the two main stellar features that help our algorithm to estimate the age and the metallicity of the stellar populations at those wavelengths: the O VIλ1032 and N Vλ1240 P Cygni stellar wind profiles (see also the C IIλ1175 photospheric line). As, in general, the FUV spectrum is dominated by young massive stars, we do not expect large changes between the STARBURST99 and other models including binaries (e.g., BPASS, Eldridge et al. 2017). Detailed results regarding the stellar populations will be discussed in a forthcoming paper (Chisholm et al., in prep.).

Fig. 1. Stellar continuum modelling for the galaxy J091703. Stellar continuum fit (red), observed spectrum (black), and error spectrum (yellow); spectral regions masked during the fit are shown in gray. The spectra, shown in Fλ units, have been normalized over 1070–1100 Å. Different ISM absorption lines and stellar features are indicated with black and dark-blue vertical lines in the top part of the figure, respectively. Geo-coronal emissions are labeled at the bottom. The best-fit dust attenuation color-excess (EB−V, in mag, R16 law), light-weighted stellar age (Myr) and metallicity (Z⊙) inferred for this source are given in the inset. Uncertainties come from a consistent Monte-Carlo error propagation along the fitting routine.

All the data sets are fitted using the same method. To conclude, the main derived quantities of interest for the present paper are the UV dust attenuation parameters (E_B−V) and the fitted stellar spectral energy distributions (SEDs). The intrinsic SEDs (before UV attenuation) are used below to derive important quantities such as the absolute LyC escape fraction ($f_{\mathrm{esc}}^{\mathrm{abs}}$) and intrinsic UV luminosities. The UV-fit $f_{\mathrm{esc}}^{\mathrm{abs}}$ values derived here are compared to other escape fraction estimates in Flury et al. (2022a,b).

Fiducial photon escape fractions. Due to the fact that the high-resolution S99 models are only well defined above the Lyman limit (912 Å), we use the low-resolution S99 models (R ∼ 100) to extend to the LyC region (Leitherer et al. 2011b). These low-resolution models are linearly combined with the same light-fractions and attenuation parameters as the previous fits.

The absolute LyC escape fraction at ∼ 900 Å (denoted by $f_{\mathrm{esc}}^{\mathrm{abs}}$^{, 900}≡ $f_{\mathrm{esc}}^{\mathrm{abs}}$ in the following) is then defined as the ratio between the modeled intrinsic LyC flux (i.e., the flux given by the low-resolution models) and the observed (measured, see Sect. 2) LyC flux for every galaxy (f^LyC):

$$\begin{array}{c}\hfill f_{\mathrm{esc}}^{\mathrm{abs}}={\displaystyle \frac{f^{\mathrm{LyC}}}{\left({\sum }_{i,j}{X}_i{F}_{i,j}^{S99}\right)(\mathrm{LyC})}}·\end{array}$$(2)

The modeled LyC flux is computed as the mean flux over the same 20 Å window as the observed LyC flux, which is typically at 850–900 Å, although varying from galaxy to galaxy (Flury et al. 2022a). For galaxies with nondetection in the LyC, we adopt the corresponding 1σ COS/G140L upper limit. Error bars on $f_{\mathrm{esc}}^{\mathrm{abs}}$ come from the simultaneous sampling of f^LyC and the modeled flux. In the case of LyC detections, these error bars are dominated by the uncertainty in the modeled flux. Figure 2 shows the low-resolution S99 modeled flux for J091703 (blue shaded region; see the figure caption) and emphasises this approach. The wavelength window over which the LyC flux is measured is shown by the gray region. Dividing the mean observed flux value in this window (black square) by the modeled intrinsic stellar flux (purple square) yields the absolute photon escape of this source ($f_{\text{esc}}^{\text{abs}}=0.{161}_{-0.055}^{+0.073}$).

Fig. 2. Stellar continuum modelling in the LyC region for the J091703 galaxy. The high-resolution stellar continuum fit (down to 925 Å) and the observed J091703 flux and error are plotted in red, black, and yellow lines, respectively. The stellar continuum is now reproduced with models of lower resolution (blue line and shaded area), allowing us to go below the Lyman limit (< 912 Å, dashed vertical line). The intrinsic (unreddened) flux is shown by the purple dashed line. In the LyC window (gray-shaded), the measured, modeled, and modeled-intrinsic mean LyC fluxes are indicated with filled squares and corresponding error bars. The ISM, stellar lines, and geo-coronal lines are displayed with the same color coding as in Fig. 1.

Throughout this paper, we use the definition of the photon escape given by Eq. (2). However, we note that similar trends and results are obtained if we use other methods to determine $f_{\mathrm{esc}}^{\mathrm{abs}}$ (see Flury et al. 2022b). The main quantities (E_B−V attenuation, stellar age, metallicity, and LyC escape fraction) derived for our sample are listed in Appendix A.

Together with the LyC flux significance (i.e., the S/N in the LyC flux), the inclusion of the escape fraction makes it possible to distinguish between strong, weak, and nonleakers. As in Flury et al. (2022b), we define strong leakers as those sources with secure LyC detection (significance > 5) and $f_{\text{esc}}^{\text{abs}}\ge 0.05$. Weak leakers are considered to be galaxies with fair or marginal LyC detection (significance ≥ 2), but not falling in the strong category. Finally, nonleakers are defined as undetected galaxies in the LyC (significance < 2). The total sample of 89 galaxies is then divided into 20, 30, and 39 strong, weak, and nonleakers, respectively (see also Appendix A).

3.2. Absorption lines

3.2.1. Equivalent widths and residual fluxes

In the rest-frame spectral range covered by our observations (∼800–1500 Å), several lines of interest are present. Here we choose the following: five Lyman series lines of hydrogen, from Lyβ to Ly6 (H Iλ1026, H Iλ973, H Iλ950, H Iλ938, H Iλ930), and six low-ionization state (LIS) lines (Si IIλ989, Si IIλ1020, Si IIλλ1190,1193, Si IIλ1260, O Iλ1302, C IIλ1334). Given the actual COS spectral resolution, some of the mentioned lines are blended with other lines at adjacent rest-frame wavelengths. In particular, the Lyman series lines are blended with O I transitions, and O Iλ1302 is also folded with the Si IIλ1304 line. Si IIλ989 is blended with O Iλ988 and N IIIλ989.

Our goal is to measure the equivalent widths and residual fluxes of these different lines (when present) for all the galaxies in the LzLCS and literature samples. We first find the wavelength of the minimum depth of the line (minimum flux). We then define a ±1250 km s⁻¹ integration window in velocity space (Δλ_i, chosen by visual inspection) centered on this wavelength, over which the equivalent width W_λi is computed as:

$$\begin{array}{c}\hfill W_{{\lambda }_i}={\displaystyle {\int }_{\mathrm{\Delta }{\lambda }_i}}(1-f_{\lambda }/F^{\star }(\lambda ))\mathrm{d}\lambda ,\end{array}$$(3)

where f_λ is the observed spectral flux density, and F^⋆(λ) is the modeled stellar continuum. We also compute the residual flux of the line, R(λ_i) = ⟨f_λ/F^⋆⟩, as the median over an interval of ±150 km s⁻¹ around the minimum flux (to avoid biasing towards minimum depth pixels affected by noise). Errors on the observed spectra plus the error on the stellar continuum fitting are incorporated within every absorption line uncertainty. More precisely, we take the median and standard deviation of the W_λi and R_λi distributions (over the 500 realizations of the perturbed spectra and subsequent continuum). We consider these as the final measurement and error values of each line.

The drop in the COS/G140L sensitivity redward of 1600 Å in the observed frame reduces the S/N of the measured flux, thus providing poorer constraints on the stellar continuum and lowering the S/N of the LIS lines at wavelengths longer than Lyαλ1216 (∼ 1600 Å at z = 0.3). Then, and only for the LIS measurements, we prefer to locally model the stellar continuum by a linear function. Equivalent widths and residual fluxes for Si IIλ1260, O Iλ1302, and C IIλ1334 lines are measured assuming this linear continuum, and using the same previous equations.

Ideally, we also want to determine the covering fraction of the lines. In general, the relation between the residual flux and the covering fraction depends on (1) the assumed gas and dust geometries, (2) the different velocity components of the gas clouds that potentially contribute to the line shape, and (3) the degree of saturation of the line, which we discuss in the forthcoming sections.

3.2.2. The covering fraction: assumptions about the ISM geometry

The escape of ionizing radiation and the UV absorption features of a galaxy can be explained by the so-called picket-fence model (e.g., Heckman et al. 2001, 2011; Vasei et al. 2016; Reddy et al. 2016b; Gazagnes et al. 2018; Steidel et al. 2018). This model describes a stellar ionizing radiation field surrounded by a nonuniform or patchy ISM, where the neutral and metallic gas are distributed in high-column-density clouds, with low-column-density channels in between. Two gas and dust relative geometries can be assumed. In the first, a uniform foreground dust-screen envelops the whole galaxy (A) while in the second, dust is hosted in the discrete gas clouds of the ISM (B).

A key parameter of the picket-fence model is the geometrical distribution of the gas parameterized by the covering fraction, C_f(λ_i). This variable represents the fraction of sight lines that are covered by optically thick gas at wavelength λ_i, in contrast to sight lines with low amounts of gas or no gas at all. Observationally, the covering fraction can be related to the residual flux of individual absorption lines (R(λ_i) = ⟨f_λ/F^⋆⟩). The functional form of this dependency changes with the gas and dust geometry (e.g., Vasei et al. 2016; Gazagnes et al. 2018), and always assumes absorption lines are saturated (i.e., the column density of the medium is high enough so that it is optically thick at the λ_i transition). Furthermore, given the limited spectral resolution, we assume that the lines are described by a single gas component or, in other words, that all velocity components of the gas have the same covering fraction.

For geometry A, the gas covering is simply the complement of the residual flux, and can be measured from the line depths as

$$\begin{array}{c}\hfill R({\lambda }_i)=1-C_f({\lambda }_i),\end{array}$$(4)

while for geometry B, the relation depends also on the dust attenuation (see detailed derivation in Vasei et al. 2016; Reddy et al. 2016b; Gazagnes et al. 2018):

$$\begin{array}{c}\hfill R({\lambda }_i)={\displaystyle \frac{1-C_f({\lambda }_i)}{C_f(\text{H}{\displaystyle \text{I}})\times {10}^{-0.4k_{{\lambda }_i}{E}_{\mathrm{B}-\mathrm{V}}^B}+(1-C_f(\text{H}\text{I}))}},\end{array}$$(5)

where $E_{\text{B-V}}^B$ is the dust attenuation measured according to geometry B, and C_f(H I) is the neutral gas covering fraction. The H I and the dust have the same covering fraction according to this model. Both the stellar continuum and the residual flux of the lines must be determined consistently with the same geometry (under the same gas and dust distribution), but the actual S/N for most of the spectra does not allow a stellar continuum fitting using geometry B (Gazagnes et al. 2018; Steidel et al. 2018). Therefore, in this paper we adopt geometry A. However, our absorption line calculations will provide insights into the underlying gas and dust geometry and other ISM properties of low-z LyC-emitting galaxies (see Sect. 5).

3.2.3. Saturation regime for the H I and LIS lines

Applying Eq. (4) (or Eq. (5)) to compute the gas covering fraction intrinsically assumes that absorption lines are saturated. The condition of saturation is tested using the equivalent width ratios for lines of the same ion (see Appendix B).

The saturation of H I is studied through the ratio of the Lyman series lines Lyβ (H Iλ1026) and Lyδ (H Iλ950). The Lyβ-to-Lyδ equivalent width comparisons reveal that most of the sources fall in the optically thick (saturation) region of the diagram (see Draine 2011), while some of them may be compatible with an optically thin regime, although many of them have large uncertainties at the same time. We caution that some of the discrepancies in this analysis can be due to slight underestimations of the continuum level, or because of the O I contribution, whose lines are blended with the Lyman series. For the metals, we use the Si II as a representative ion, and compare the ratio of Si IIλ1260 to Si IIλλ1190,1193. Nevertheless, this comparison is subject to large measurement errors and, although we can be statistically sure about the H I saturation, we can neither discard nor assume saturation for the metallic lines.

Henceforth, we assume that all the H I transitions and metal lines considered are saturated. A more detailed explanation of this reasoning is presented in Appendix B.

3.2.4. Estimation of systematic errors

We finally study the possible systematic errors on our measurements. When measuring absorption line properties from observed spectra, the instrumental resolution (R) tends to make the lines broader and so residual fluxes can be overestimated. In addition, the observed S/N of every pixel – due to instrumental sensitivity, noise, etc. – also impacts the line shape and can therefore affect our estimates.

To account for these effects, we simulate different absorption lines (Lyδ, Lyβ, Si IIλ1260, C IIλ1334) with a standard Voigt profile (Draine 2011) and typical saturated column densities and Doppler broadening in SF galaxies, a single velocity component, and a fixed covering fraction of C_f = 0.85 for all of them. A foreground dust-screen geometry of the picket-fence model is adopted. We then degrade to the COS/G140L instrumental resolution (R ∼ 1000) and finally add a set of different random Gaussian S/N per pixel. After that, the residual flux is measured 500 times for every S/N value, and converted to a covering fraction using Eq. (4). The median covering fraction is compared to the original input value to give the relative difference. We conclude that, for the typical S/N measured in our sample at the continuum level close to the lines in question (e.g., close to H Iλ950 the continuum is S/N ∼ 4) the systematic error on the covering fraction spans between 5% and 20%, with a bias towards lower coverings. However, this value does not dominate over original errors and is within the original error bars we are reporting for single covering fraction measurements (see Sect. 5). In addition, given the actual instrumental resolution, systematic deviations do not change significantly with respect to C_f = 0.85 when other input C_f is chosen in the simulations. Even in the case when the residual flux of single lines could be slightly overestimated, the general trends in this study will hold, because all the spectra in the working sample have the same resolution. We therefore leave our initial residual flux measurements and errors unaltered. A more exhaustive explanation of the estimation of systematic errors is presented in Appendix C.

3.3. Photon escape fraction from the Lyman series and LIS metallic lines

In the above-described picket-fence model, the absolute LyC photon escape fraction is related to the H I covering fraction of saturated lines and the dust attenuation (Heckman et al. 2001, 2011; Vasei et al. 2016; Reddy et al. 2016b; Gazagnes et al. 2018, 2020; Chisholm et al. 2018; Steidel et al. 2018). The expression of the predicted $f_{\mathrm{esc}}^{\mathrm{abs}}$ depends on the assumed gas and dust geometry. Nonetheless, as demonstrated by Gazagnes et al. (2018), Chisholm et al. (2018), although E_B−V and C_f(H I) are model dependent, regardless of the use of geometry (A) or (B), the attenuation and covering fraction terms balance each other so that the predicted $f_{\mathrm{esc}}^{\mathrm{abs}}$ remains almost unaltered (e.g., it does not strongly depend on the assumed geometry). This is because, for the same observed flux, model (B) would require larger dust obscuration but also higher gas covering fraction than model (A) in order to equally match the data.

Consequently, we decide to adopt a uniform dust screen model (A) consistently for the continuum and the $f_{\mathrm{esc}}^{\mathrm{abs}}$ determination from the absorption lines. Covering fractions are therefore simply given by Eq. (4), and finally the predicted photon escape from H I is (cf., Gazagnes et al. 2018; Chisholm et al. 2018)

$$\begin{array}{c}\hfill f_{\mathrm{esc}}^{\mathrm{abs},\mathrm{HI}}={10}^{-0.4k_{912}{E}_{\mathrm{B}-\mathrm{V}}}\times (1-C_f(\text{H}\text{I})),\end{array}$$(6)

where k₉₁₂ depends on the assumed attenuation law at 912 Å (e.g., k₉₁₂ ≈ 12.87 for R16). C_f(H I) is the average H I covering fraction. To compute this, we first use Eq. (4) to figure out individual covering fraction measurements for each line. Then, in order to improve the S/N in single measurements of C_f(λ_i), we use an inverse-variance-weighted mean definition for the H I covering fraction:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill C_f(\text{H}\text{I})& \pm \mathrm{\Delta }{C}_f(\text{H}\text{I})={\displaystyle \frac{{\sum }_{i=1}^N{C}_{\mathit{fi}}\times \frac{1}{{\sigma }_{C_i}^2}}{{\sum }_{i=1}^N\frac{1}{{\sigma }_{C_i}^2}}}\pm \sqrt{\frac{1}{\sum _{i=1}^N{\sigma }_{C_i}^2}},\hfill \\ \hfill i& \equiv \text{H}\text{I}\lambda 930,\text{H}\text{I}\lambda 938,\text{H}\text{I}\lambda 950,\text{H}\text{I}\lambda 973,\text{H}\text{I}\lambda 1026\hfill \end{array}\end{array}$$(7)

where N ≤ 5 is the number of Lyman lines that can be measured on a single spectra, and σ_Ci is the measurement error (from Monte-Carlo sampling of the spectra+continuum, see previous section). These methods are tested against the fiducial photon escape in Sect. 7.

4. UV attenuation and absolute photon escape of the LzLCS galaxies

From the direct output parameters of the spectral fits we obtain information on the light-weighted stellar age and metallicity of the stellar population, and the UV dust attenuation. The UV dust attenuation significantly impacts the derived absolute escape fractions, and so here we mostly focus on the relation between these two variables. The stellar population properties of the LzLCS sample will be discussed in a separate paper (Chisholm et al., in prep.).

The fitted E_B−V ranges between ∼0.02 and 0.45, with a mean value of ∼0.15 (0.07 and 0.24 as the 0.16 and 0.64 quantiles). Figure 3 shows the E_B−V derived from the UV–SED fits and a comparison to the ones determined from the Balmer decrement measurements using the optical SDSS spectra and the Cardelli et al. (1989) attenuation law (CCM). For the E_B−V from the UV fits, we examined two different attenuation laws: R16 and SMC, chosen to bracket the expected UV attenuation in SF galaxies³ (Shivaei et al. 2020). A direct comparison between E_B−V values obtained from SMC and CCM is straightforward, because both laws have identical k_λ at optical wavelengths. In the case of R16, we used the Reddy et al. (2015) parametrization at λ > 1500 Å converted to CCM law via the usual Balmer Decrement formulae⁴. When R16 is used, the overall E_B−V derived from the UV spectra is higher than that from the CCM law typically by a factor of ∼1.2, while it is a factor of ∼2.6 lower when SMC is chosen.

Fig. 3. Comparison between the optical attenuation parameter EB−V (color excess, in mag) obtained from the Balmer decrements (CCM attenuation law) and the values derived from the UV–SED fits to the COS spectra. Here we explore two different attenuation laws: R16 (red) and SMC (blue symbols). Median error bars are plotted in the upper left. Circles represent the LzLCS galaxies while diamonds are the Izotov et al. (2016a,b, 2018a,b, 2021) and Wang et al. (2019) sources. The theoretical ×0.44 shift between the stellar and nebular attenuation, which assumes a foreground screen of dust (Calzetti et al. 2000), is also shown.

A systematic shift between the nebular (obtained from emission line ratios) and the stellar attenuation (from the slope of the continuum) may be expected when a foreground dust-screen distribution is in play (see Calzetti et al. 2000). The stellar attenuation is typically a factor of 0.44 lower than the nebular attenuation (see Fig. 3). We only find this trend when the SMC law is used; R16 derived values agree better with the E_B−V from the Balmer decrements. Moreover, an SMC-like attenuation law has appeared more suitable for high-z SF galaxies (Salim et al. 2018) and low-metallicity starburst (Shivaei et al. 2020), requiring even steeper laws for leakers (Izotov et al. 2016b).

Distinguishing the attenuation law from the data used here is not possible, because the fit quality (reduced-χ²) is similar in both cases. Therefore, and as we mentioned in Sect. 3, we use the R16 attenuation law as a default in this study, motivated by the fact that this law is defined below 1000 Å by a significant number of galaxies. This translates to a factor two lower attenuation at 900 Å than estimates derived from extrapolations of dust attenuation curves such as the SMC, and will impact the determination of the escape of ionizing photons from SED fitting, as we explore below (see also Reddy et al. 2016b; Chisholm et al. 2019). Finally, to avoid a strong dependence on the assumed attenuation law, we refer to the UV attenuation defined as A_λ = k_λE_B-V (Eq. (1)) instead of E_B−V. Still, differences in A_λ are found for the two laws at all wavelengths, being systematically higher for R16. This reveals that the fitted stellar parameters (X_i) are slightly lower for the SMC law, meaning that both fits can match the observations similarly. In any case, both approaches yield the same light-averaged ages and metallicities (further investigations in Chisholm et al., in prep.).

In Fig. 4 we show the derived $f_{\mathrm{esc}}^{\mathrm{abs}}$ values (Eq. (2)), assuming both the R16 and the SMC law (see also Flury et al. 2022a). The derived LyC escape fractions vary over a wide range, with values up to 60% (85%) for the LzLCS (literature) sample. Typically, upper limits are below 0.01. Clearly, the choice of the attenuation law has some impact on the derived $f_{\mathrm{esc}}^{\mathrm{abs}}$: the relative difference between the two $f_{\mathrm{esc}}^{\mathrm{abs}}$ versions ranges from 2% to 40% ×$f_{\text{esc}}^{\text{abs}}$, with systematically higher values obtained for the SMC law. This systematic shift is smaller for high values of the absolute escape fraction, which are found in galaxies with relatively low UV attenuation (Sect. 6.4). The shift to higher $f_{\mathrm{esc}}^{\mathrm{abs}}$ for the SMC law is due to fact that this law is steeper than R16, which means that a smaller UV attenuation is needed to reproduce the overall, reddened observed spectral shape. Therefore, the intrinsic stellar flux (before attenuation by dust) is lower with the SMC law, and therefore $f_{\mathrm{esc}}^{\mathrm{abs}}$ is higher.

Fig. 4. 1:1 comparison between the absolute LyC photon escape fractions ($f_{\mathrm{esc}}^{\mathrm{abs}}$) from the UV–SED fits using R16 and SMC dust-attenuation laws. LzLCS leakers and nonleakers are indicated by blue points and black arrows (upper limits), respectively. Red squares are the Izotov et al. (2016a,b, 2018a,b, 2021) and Wang et al. (2019) sources. Median relative error bars are shown in the lower right corner.

Finally, our fitting algorithm is not devoid of degeneracies between stellar age, metallicity, and the E_B−V (see discussion in Gazagnes et al. 2020). The S/N of the original spectra impacts the ability of the fitting code at the time of solving this degeneracy, where only moderate (and high) S/N spectra will give a reliable estimation of the stellar population properties (Chisholm et al. 2019).

5. Empirical results and inferences from UV absorption lines

We now search for the empirical relations between the H I/LIS absorption lines, their equivalent widths (Sect. 5.1), residual fluxes, and some galaxy properties (Sect. 5.2). We also compare the observed trends with those in high-z galaxies (Sect. 5.3).

5.1. LyC leakers have weak absorption lines

The equivalent width (W_λ) distributions of the major hydrogen and metal absorption lines are shown in Fig. 5, where the strong leakers (dark gray) are separated from the weak and nonleaker distributions (light gray). The error bars on the histograms come from a proper estimation of the Poissonian confidence intervals for every W_λ bin, taking into account the uncertainty on single measurements.

Fig. 5. Equivalent width distributions (Wλ, in Å) for the measured absorption lines. Dark-gray histograms correspond to the strongest leakers (20) in the LzLCS plus published joint samples, while light-gray histograms represent the weak leakers (30) and nondetections in the LyC window (39), out of the total of 89 galaxies. Strong leakers are defined as galaxies with significant detection in the LyC and $f_{\text{esc}}^{\text{abs}}$ ≥ 0.05. Weak leakers also have significant LyC detection but $f_{\text{esc}}^{\text{abs}}$ < 0.05. Nonleakers are not detected in the LyC. Median equivalent-width errors are indicated in all the panels and for each source type.

Overall, the equivalent widths reach up to 4–6 Å for the strongest lines. The number of detections (defined here as value/error  > 2) out of the total of 89 galaxies in the sample are 62, 51, 54, 56, and 50 for the Lyman series, and 32, 35, 23, 20, 13, and 9 for the LIS lines, respectively. The main cause for nondetection is the contamination by geo-coronal emission in the Lyman series range and the decrease in S/N towards longer wavelengths. The presence of negative equivalent widths can be due to different reasons. For some H I lines, such as H Iλ1026, the continuum near the line is complex to model because of contributions from a stellar O VI P Cyngi profile as well as numerous adjacent ISM absorption lines (in particular O VI, C II, O I). Slight underestimations of the continuum level would lead to negative equivalent widths. In the case of the LIS lines, the low S/N and overlap of some fluorescence transitions within the integration window would also produce some negative W_λ values.

The equivalent width distributions of leakers peaks at lower values than the nonleaking distributions in general, suggesting that leakers have weaker H I and LIS absorption lines than nonleakers. Table 1 shows the median of the W_λ distributions for the entire sample (89 galaxies), strong leakers (20), and weak+nonleakers only (69). For example, for H Iλ1026 (Si IIλλ1190,1193), the median equivalent width value for leakers is 0.58 Å (0.17 Å), while nonleakers show median equivalent widths of 1.97 Å (1.21 Å), respectively. To test whether the two distributions differ or not, one can perform the Anderson-Darling test. For the previous lines, the null hypothesis that assumes leaker and nonleaker populations to come from the same parent distribution can be rejected at 97.5% level, with a significance of p_val ≤ 10⁻². Based on smaller samples Chisholm et al. (2017) already noted that LyC emitters show weaker H I and LIS lines (see also Mauerhofer et al. 2021, in simulations).

5.2. H I and LIS absorption line strengths and UV attenuation are driven by the gas covering fraction

The observed differences in the H I and LIS lines between LyC emitters and nonemitters suggests a connection between the mechanisms that favor the LyC leakage and those that lead to changes in the absorption lines profiles. Therefore, we explore several empirical correlations between H I and LIS equivalent widths and other galaxy properties, such as the intrinsic absolute magnitude at 1500 Å ($M_{1500}^{\mathrm{int}}$), the UV attenuation (A_UV), the measured gas-phase metallicity (12 + log(O/H)), and the average H I and LIS residual flux (R(H I), R(LIS)). $M_{1500}^{\mathrm{int}}$ is produced by taking the mean of the dust-free S99 modeled flux over a 100 Å window around 1500 Å, and converting this into an AB absolute magnitude given the redshift and the distance modulus. A_UV is taken from the fits, A_UV = k_UV × E_B-V – at 1500 Å –, and finally residual fluxes are computed using Eq. (4) and subsequently Eq. (7) for H I and similarly for the LIS lines (Si IIλλ1190,1193, Si IIλ1260, O Iλ1302, C IIλ1334).

For the equivalent widths, we take the inverse-variance-weighted average of the single H I and LIS line widths (W_H I, W_LIS):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill W_K& \pm \mathrm{\Delta }{W}_K={\displaystyle \frac{{\sum }_{i=1}^N{W}_{\lambda i}\times \frac{1}{{\sigma }_{W_i}^2}}{{\sum }_{i=1}^N\frac{1}{{\sigma }_{W_i}^2}}}\pm \sqrt{\frac{1}{\sum _{i=1}^N{\sigma }_{W_i}^2}},\hfill \\ \hfill K& \equiv \text{H}\text{I}\hfill \\ \hfill i& \equiv \text{H}\text{I}\lambda 930,\text{H}\text{I}\lambda 938,\text{H}\text{I}\lambda 950,\text{H}\text{I}\lambda 973,\text{H}\text{I}\lambda 1026\hfill \\ \hfill K& \equiv \mathrm{LIS}\hfill \\ \hfill i& \equiv \text{Si}\text{II}\lambda \lambda 1190,1193,\text{Si}\text{II}\lambda 1260,\text{O}\text{I}\lambda 1302,\text{C}\text{II}\lambda 1334,\hfill \end{array}\end{array}$$(8)

where N is the number of H I or LIS lines that can be measured on a single spectra and σ_Wi represents the error on the individual equivalent widths measurements (Monte-Carlo sampling). At this point it becomes necessary to take into account the following aspects. Some lines, such as H Iλ973 and Si IIλ989, are contaminated by geo-coronal emission in many cases (those values are removed from the average). We exclude Si IIλ1020 from the average because this line has a very low oscillator strength and is usually very weak and likely optically thin. We also exclude Si IIλ989 because of the geo-coronal contamination and the blending with the high-ionization-state N IIIλ989 line.

Averaging the H I and LIS lines of the same ion is reasonable owing to the fact that most of those lines have similar oscillator strengths or are likely saturated, and so their equivalent widths are similar. Even LIS transitions for very different ions (e.g., Si II versus O I) will only differ by a factor of two in SF galaxies (e.g., Chisholm et al. 2016). Averaging LIS line measurements has also frequently been reported for high-redshift galaxies (e.g., Shapley et al. 2003; Du et al. 2018, 2021; Trainor et al. 2019) to improve the S/N of single line measurements. As mentioned above, in our case the average is strictly necessary for the LIS lines redder than 1200 Å rest-frame (or below 1600 Å in the observed frame), where the S/N is usually very low due to the decrease in the COS instrument sensitivity.

In Figs. 6 and 7, the H I and LIS equivalent widths are plotted against the four quantities described above: $M_{1500}^{\mathrm{int}}$, A_UV, 12 + log(O/H), R(H I), and R(LIS). LzLCS sources are represented by circles while the published sample is shown using diamonds. The sources are color coded by our $f_{\mathrm{esc}}^{\mathrm{abs}}$ determined from the UV spectral fits. Overall, strong absorption lines are only found at high A_UV, high metallicity, and in UV-luminous sources (where $M_{1500}^{\mathrm{int}}$ is dust corrected). In other words, $W_{\text{H}\text{I}}$ and W_LIS both correlate with the UV attenuation and metallicity, but anti-correlate with the intrinsic UV absolute magnitude and residual fluxes. High $f_{\mathrm{esc}}^{\mathrm{abs}}$ sources are located at the lowest H I/LIS equivalent-width regime (see also Sect. 6.2). The variations in the absorption line properties with the LyC escape fraction are discussed further below.

Fig. 6. Empirical trends between the H I equivalent width (WH I) and different observational properties (from top to bottom and left to right): intrinsic absolute magnitude at 1500 Å ($M_{1500}^{\mathrm{int}}$), UV attenuation (AUV) at 1500 Å, measured gas-phase metallicity (12 + log(O/H)) and the average H I residual flux (R(H I)). Points are color-coded by their absolute LyC photon escape fraction, $f_{\mathrm{esc}}^{\mathrm{abs}}$ (and upper limits values on $f_{\mathrm{esc}}^{\mathrm{abs}}$ are indicated inside the color bar for reference). The Izotov et al. (2016a,b, 2018a,b, 2021) and Wang et al. (2019) sources are displayed with diamonds and also considered in the correlations (see text).

Fig. 7. Empirical trends between the LIS equivalent width (WLIS) and different observational properties: $M_{1500}^{\mathrm{int}}$, AUV, 12 + log(O/H) and R(LIS), with sources displayed and color-coded as in Fig. 6. Additionally, results from Shapley et al. (2003, z ∼ 3) and Harikane et al. (2020, z ∼ 6) composites are overplotted. The light-blue, blue, and dark-blue lines correspond to the different redshifts between z = 2 − 4 from Du et al. (2018).

In order to quantify the above-mentioned trends, we computed the Kendall correlation coefficient for every pair of variables described above (τ, see Table 2). For a sample of 89 objects, we consider correlations to be significant if p_val ≲ 2.275 × 10⁻² (i.e., 2σ) and strong if |τ| ≳ 0.2. For W_H I, the strongest and most significant correlation is found to be with the UV intrinsic absolute magnitude $M_{1500}^{\mathrm{int}}$. On the other hand, W_LIS shows the strongest trend with the LIS residual flux.

Generally, three different mechanisms can change the H I/LIS equivalent widths – if we assume the lines are formed in single gas clouds, with a single velocity and gas column density: (1) the H I/metal column densities and metallicity, (2) the Doppler broadening (velocity dispersion of the cloud), and (3) the H I/LIS covering fraction. Doppler broadening (2) is challenging to address, given our spectral resolution. This said, and based on visual inspection, we do not see any evidence for the highest equivalent widths being among the widest lines, but they are among the deepest. The effect of the covering fraction (3) is clearly visible in the bottom right panels of Figs. 6 and 7, where the residual flux is linked to the covering fraction – for a uniform dust geometry of the picket-fence model – and so the equivalent width is directly related to the line depth. It remains to be seen whether or not other physical effects can also play a role.

For example, in a study of nearby SF galaxies, Chisholm et al. (2017) suggested that the decrease in the W_LIS is set by a combination of small H I column densities and low metallicity. The relation between the H I column density (N_H I) and the column density of metals (N_Z) scales with metallicity (Z) as N_z ∼ Z × N_H I. A decrease in metallicity (Z) would produce a weaker W_LIS, leading to a natural W_LIS − A_UV dependency. However, a change in metallicity cannot alter the H I line widths and thus explain the observed $W_{\text{H}\text{I}}$– A_UV trend. Furthermore, as the H I lines must be saturated – based on the ratio of equivalent widths of different H I and Si II transitions (Appendix B) – column density must be secondary. The observations of the Lyman series lines and the fact that they behave in a very similar fashion to those of the LIS lines therefore exclude option (1) to explain the trends shown in Figs. 6 and 7.

The effect of metallicity is examined in more detail by restricting our analysis to the sources at a constant gas phase metallicity of 12 + log(O/H) = 8.2 ± 0.1 (25 galaxies). Doing so, the trends seen in Figs. 6 and 7 remain the same. As an example, in Fig. 8 we show the resulting W_LIS − R(LIS) behavior when the sample is limited in metallicity (colored points) against the whole LzLCS+literature sample (gray points in the backgorund). Quantitatively, a Kendall test performed over the W_LIS − A_UV relation for example (see inset in Fig. 8), still reveals a similar strong and significant correlation when the sample is restricted in metallicity (τ ∼ 0.25, p_val ∼ 10⁻³). On the other hand, when restricting to constant R(LIS), the statistical connection between W_LIS and A_UV disappears, and similarly for the other variables. This excludes the metallicity (1) as the parameter driving the LIS correlations. We therefore conclude that variations of the geometric covering fraction (3) are the most likely parameter driving the observed trends both for the H I and LIS lines, and also for the UV attenuation (with a possible cloud-to-cloud velocity dispersion, and the different column densities in front of the stars contributing to the scatter). These findings are similar to the conclusions reached in (Du et al. 2021) from the analysis of UV absorption lines and attenuation in LBGs at z ∼ 2 − 3. Du et al. (2021) find that metallicity only affects the LIS equivalent widths when the W_Lyα is fixed at a constant value or, in other words, when the H I covering fraction is also roughly constant.

Fig. 8. LIS equivalent widths (WLIS) versus residual fluxes (R(LIS), main panel) and UV attenuation (AUV, inset), when the working sample is restricted to a constant metallicity of 12 + log(O/H) = 8.2 ± 0.1 (color-coded as in Fig. 6). The global trend including all the sources is shown with gray background symbols.

5.3. Comparison with high-redshift galaxies

Few studies have explicitly shown high-order Lyman series trends with other observables (Henry et al. 2015; Reddy et al. 2016b; Gazagnes et al. 2018, 2020; Steidel et al. 2018; McKinney et al. 2019), because these lines are in the extreme FUV and hard to detect, often only accessible in the nearby Universe (e.g., see Lecavelier des Etangs et al. 2004; Lebouteiller et al. 2006; Grimes et al. 2009, using the FUSE satellite). In general, our H I and LIS equivalent width lines correlate similarly with the other galaxy properties, supporting a direct link between neutral gas and metallic gas distributions (see next Sect. 6.2), and the claim that the neutral gas content of galaxies is well traced by LIS species (Shapley et al. 2003). The similarity in the behavior of H I and LIS lines also suggests that metallicity cannot be the major contributor to these trends, whose effect would no longer be visible if the lines were saturated.

In Fig. 7, we compare the observed W_LIS and UV magnitudes with results from Jones et al. (2012) and Harikane et al. (2020), who used LBG composite spectra at z ∼ 4 and z ∼ 6, respectively. Here, Jones et al. (2012) and Harikane et al. (2020) reported values are not dust-corrected, i.e., they are observed AB absolute magnitudes. As seen in Fig. 7 (top left), the agreement is fair with the high-z literature results, suggesting that the ISM structure and properties of the low-z galaxies might be similar to those at high redshift.

Continuing, similar empirical relations between W_LIS and A_UV (or E_B−V) were found by Shapley et al. (2003), Du et al. (2018) and Pahl et al. (2020), although the relations are generally shifted to lower reddening values. In the same Fig. 7, we also show the results from Shapley et al. (2003)⁵. The continuum sample selection in Shapley et al. (2003) against our continuum plus emission line selected sample could be a possible reason for higher W_λ at fixed A_UV. Moreover, the use of a broad-band SED fitting approach to infer the reddening or the spectral stacking could lead to the current differences and, less importantly, the use of different attenuation laws. Leitherer et al. (2011b) and Grimes et al. (2009) show similar trends with the 12 + log(O/H), but shifted to higher metallicities (typical of local starburst galaxies). The trend between the W_LIS and LIS residual fluxes (bottom right panel) is compatible with the results of Harikane et al. (2020) who use a z ∼ 6 galaxy composite. Finally, another possible reason for the small systematic differences in the observed equivalent widths at low- and high−z could be the different spectral resolution and aperture effects, which would systematically lead to stronger observed lines for the low-z galaxies.

6. Relations between the UV absorption lines, Lyα, and LyC escape

Here, we present the empirical relations between the UV absorption line, the Lyα emission line, and the Lyman continuum escape (Sect. 6.1). We then investigate the connection between H I and LIS covering fractions (Sect. 6.2). Finally, we study the effect of the neutral and low-ionized gas covering fraction on the leakage of ionizing radiation (Sect. 6.3) and the influence of dust (Sect. 6.4).

6.1. Empirical correlations between the H I and LIS absorption lines, Lyα and LyC

Although Lyα and LyC photons interact with the same neutral gas, gas kinematics can influence the Lyα output, while it has almost no effect (directly) on LyC emission. The interplay between the Lyα emission and the LyC leakage has been widely studied (Rivera-Thorsen et al. 2015, 2017; Verhamme et al. 2015, 2017; Dijkstra et al. 2016; Marchi et al. 2018; Izotov et al. 2020, 2021; Gazagnes et al. 2020). Hence, it is a natural step to also study the behavior of the H I/LIS lines with the Lyα line. We use the Lyα equivalent widths (W_Lyα) measured by the LzLCS team from the COS spectra (Flury et al. 2022a).

In Fig. 9, the W_Lyα is compared to the different H I and LIS line properties. Very strong and significant correlations (first column) are found between W_Lyα and W_H I, and with W_LIS (see Table 2). W_Lyα and R(H I) are also moderately correlated (however they have larger uncertainties). In summary, larger Lyα equivalent widths are found in galaxies where the weaker H I and LIS absorption lines are, that is, W_Lyα is anti-correlated with $W_{\text{H}\text{I}}$ and W_LIS but correlates with R(H I) and R(LIS). The color-coding demonstrates the correlation between W_Lyα and $f_{\mathrm{esc}}^{\mathrm{abs}}$. Galaxies with higher $f_{\mathrm{esc}}^{\mathrm{abs}}$ are statistically the strongest Lyα emitters (LAEs), that is, they show large equivalent widths Flury et al. (2022b). Additionally, larger scatter at the high-W_Lyα end at fixed W_H I and W_LIS may result from boosted (or at least variations in) intrinsic Lyα production.

Fig. 9. Relation between the Lyα (WLyα), H I (top), and LIS line (bottom) equivalent widths and residual fluxes, respectively. Sources are displayed and color-coded as in Fig. 6. Solid and dashed lines represent the linear fits to the relations found in Steidel et al. (2018), Du et al. (2018), and Gazagnes et al. (2020) for comparison. Gray stars, the green polygon, and green crosses correspond to Shapley et al. (2003), Jones et al. (2012, 2013) results.

Steidel et al. (2018) and Gazagnes et al. (2020) studied the connection between the Lyα equivalent width and the H I residual flux. In Fig. 9 (top right) we show how their fitted linear relations compare to our data (gray and black lines). Both correlations are less steep than ours, but a more prominent difference exists when comparing with Steidel et al. (2018). As discussed in Gazagnes et al. (2020), sample selection is likely the key to understand these differences. Steidel et al. (2018) photometrically selected LBGs with rest-frame UV bands, meaning that they are generally brighter, and have higher masses and/or higher SFRs than the low-redshift galaxies studied by Gazagnes et al. (2020) and here. Another possible explanation is the scattering of the blue side of the Lyα line by the IGM, which would produce lower W_Lyα in the Steidel et al. (2018) sample. Additionally, flux aperture losses would also lead to a decrease in the W_Lyα of the former sample with respect to the LzLCS.

Finally, important to note is the bias introduced by the redshifts of the two samples: from nearby z ∼ 0.3 to distant z ∼ 3 galaxies, with intrinsic different properties. After all, the current sample may include sources with larger W_Lyα than the bulk of the LBGs at z ∼ 3, whereas in Steidel et al. (2018), sources were selected to have brighter continua.

The Lyα and LIS connection has been widely studied by stacking LBGs at z ∼ 3 − 5 (Shapley et al. 2003; Jones et al. 2012; Du et al. 2018; Pahl et al. 2020) and also with individual detections at slightly lower redshifts (Trainor et al. 2019; Du et al. 2021, at z ∼ 2 − 3). For comparison, we include the linear regression between the Lyα and LIS equivalent widths found by Du et al. (2018) (Fig. 9, bottom left), and data from stacked spectra in Shapley et al. (2003) and Jones et al. (2012). High-z LBG measurements seem to be more consistent with low-$f_{\mathrm{esc}}^{\mathrm{abs}}$ or nonleaking galaxies: strong leakers occupy a different portion of this plot, with the highest W_Lyα. Similar trends between the Lyα and metal equivalent widths were found in Jaskot et al. (2019) for individual Green Pea-like galaxies.

The W_Lyα correlation with the metal residual flux is more difficult to interpret – mainly because of the larger error bars on R(LIS) – but it is still moderately significant. This result may suggest that the trend between Lyα and LIS absorption equivalent width is predominantly due to variations in the covering fraction of neutral hydrogen (see below), which modulates the Lyα profile and is directly connected with the LIS covering fraction (see Sect. 6.2). The analysis of Gazagnes et al. (2020) indicates that the gas kinematics could also explain part of the observed scatter, as well as the distribution of different gas column densities in front of the stars. For comparison, we include the W_Lyα–R(LIS) correlation from the Jones et al. (2013) compilation of z ∼ 2 − 4 galaxies, although the latter sample lacks strong LAEs, as opposed to ours, which explains the differences.

As the LIS absorption lines are a probe of neutral gas, there is a direct link between the Lyα profile and LIS lines (Shapley et al. 2003). From the apparent nonevolution of the W_Lyα–W_LIS relation with redshift, Du et al. (2018) and Pahl et al. (2020) suggest that Lyα and the LIS absorption features are fundamentally related, indicating that the gas giving rise to LIS absorption modulates the radiative transfer of Lyα photons. We now know that, rather than metallicity, variations of the line covering fractions are the main cause for changes in W_λ. In order to examine the role of the gas covering fraction in W_Lyα from the low-redshift sample, we again restrict our sample to the sources at a constant gas phase metallicity of 12 + log(O/H) = 8.2 ± 0.1 (25 galaxies). With this restriction, we find the same trends as in Fig. 9, confirming that the observed W_Lyα–W_LIS relation is not primarily due to metallicity variations, but it is ultimately driven by variations in the H I covering fraction (with the different column densities in front of the ionizing sources contributing to the scatter). Given the locus of the strongest leakers in the W_Lyα–W_LIS plane, having large Lyα equivalent widths but weak LIS lines, this diagram can be used to select potential LyC emitters in high-redshift surveys.

6.2. H I and metals are not evenly distributed, but trace each other

Although with different ionization potentials, it is well-known that H I and LIS absorption lines in SF galaxies generally show similar kinematics (line profiles, central velocities, etc.; see Dessauges-Zavadsky et al. 2010; Trainor et al. 2015; Chisholm et al. 2016; Reddy et al. 2016b). This indicates that the LIS and neutral gas are comoving and trace similar gas content. Here, we proceed to study the relation between the observed neutral and metallic gas geometrical distributions through the covering fraction.

Covering fractions for every line were derived from the corresponding line residual fluxes using Eq. (4), and applying Eq. (7) for H I and similarly for the LIS lines. This assumes a uniform dust-screen geometry in the picket-fence model, a single velocity component, and saturated lines. In Fig. 10, the LIS covering C_f(LIS) is compared to the neutral gas covering fraction obtained from the set of Lyman series lines C_f(H I). Generally, the LIS covering fraction is lower, typically by a factor of two, compared to that of the H I lines. The median covering fraction is 0.82 for H I and 0.46 for the LIS. A linear relation between the two is clear from Fig. 10, such that the LIS covering fraction increases with increasing H I covering fraction.

Fig. 10. Empirical correlation between the LIS and H I covering fractions. Red, yellow, black, and dark blue lines correspond to the 1:1, Reddy et al. (2016b), and Gazagnes et al. (2018) relations and our Bayesian linear fit (LINMIX), respectively. The points are color-coded according to the gas-phase metallicity of each source, 12 + log(O/H), and typical error bars are shown in the bottom-right corner. The blue shaded area represents the 2σ confidence interval including the intrinsic scatter of the fit. The linear correlation coefficient (r2) and intrinsic scatter over the fit (σy) are given. The Izotov et al. (2016a,b, 2018a,b, 2021) and Wang et al. (2019) sources are displayed in diamonds, and also included in the fit.

To quantify this, we fit a linear regression model to the data using the Bayesian fitting code LINMIX⁶ (Kelly 2007), which accounts for errors on both variables, also including censored data (upper or lower limits). Further details are provided in Appendix D. We obtain

$$\begin{array}{c}\hfill C_f(\text{H}\text{I})=(0.63\pm 0.19)\times C_f(\mathrm{LIS})+(0.54\pm 0.09).\end{array}$$(9)

The resulting fit is quite significant, with a linear correlation coefficient of r² = 0.77 (+0.15, −0.22), and a negligible estimated intrinsic scatter of σ_y = 0.004.

Our result is compared to the Gazagnes et al. (2018) fit, where a similar approach was followed but for a small subset of the LyC emitters (dashed black line). Both linear relations have similar slopes, but slightly different normalizations. Possible differences could be due to the use of the mean of different LIS lines – in our case – instead of single Si IIλ1260 line in the previous work, or due to the effect of the resolution, which was corrected for by Gazagnes et al. (2018). The measured residual fluxes can be affected by the spectral resolution, and therefore the actual linear relation also depends on the resolution. However, as discussed in Sect. 3.2, the resolution effect is subdominant in this sample. Our results are also compatible with those of Reddy et al. (2016b), although spanning a wider range of C_f values. As sketched in the previous paragraph, the strong leakers are located not only with the lowest H I but also the lowest LIS covering fractions, having 0.77 (H I) and 0.41 (LIS) as median values.

Multiple explanations can lead to these covering fraction differences. As explained in Reddy et al. (2016b), gas kinematics can influence the measured residual fluxes of metal lines. If most of the metals reside in narrow and unresolved absorption regions, the limited spectral resolution could underpredict the gas covering fraction estimated from the residual intensity of the lines. On the other hand, if the metal-enriched gas has a large velocity gradient in the line of sight, the derived covering fractions may be lower limits (see discussion in Mauerhofer et al. 2021). Nevertheless, the similar blueshifted velocity and spread of H I and LIS lines suggests that both phases are not separated, but coupled in the ISM, although once again, the limited spectral resolution does not allow us to confirm this definitely. Orientation effects for lines created in the circumgalactic medium (CGM) instead, where H I and LIS are not particularly well-coupled, could also contribute to the scatter in the relation. Additionally, LIS lines could also be affected by scattering (Mauerhofer et al. 2021, in simulations) and fluorescence of LIS emission in-filling (Scarlata & Panagia 2015), although this effect is very difficult to estimate.

Gazagnes et al. (2018) suggested that a different distribution of the metals, which would also depend on metallicity, could explain the current relation. In an ISM where H I is not fully mixed with the metals, H I would probe high and low-metallicity regions, whereas metals would only probe the higher metallicity areas. Therefore, the H I covering fraction would be higher than that of the metals. To test this scenario, we color-coded the points in Fig. 10 according to the 0.25, 0.5, and 0.75 gas-phase metallicity (12 + log(O/H)) quantiles. Overall, sources with low H I and LIS covering fractions have lower metallicities, although the scatter and error bars blur this behavior. If we divide the sample into high- (12 + log(O/H) > 8.1) and low-metallicity (12 + log(O/H) ≤ 8.1) subsamples, we find average covering values of C_f(H I) = 0.85 (0.78) and C_f(LIS) = 0.52 (0.43) for high (low) metallicities, confirming the previous hypothesis. Nevertheless, high- and low-metallicity subsamples approximately follow the same linear relation predicted by Eq. (9) (similar slopes and intercepts). This excludes a significant dependence of the previous C_f(LIS) − C_f(H I) relation on metallicity, but as the measured metallicities are integrated over the whole galaxy, we cannot exclude the inhomogeneous mixing scenario.

Finally, there is still the possibility that the lower end of this C_f(H I) − C_f(LIS) correlation is mainly driven by unsaturated lines. There is no simple way to address this question rather than the analysis already made in Appendix B. In there, we conclude that the number of unsaturated lines sources must be minimal, but without certainty for the LIS lines. For example, for galaxies with 12 + log(O/H) = 7.8 (among our lower metallicities) and assuming that Si II roughly traces the H I fraction, one has log(Si/H)∼ − 5.7, adopting the nearly constant⁷ observational relation of log(Si/O)∼ − 1.5 (see Izotov & Thuan 1999). According to a curve-of-growth analysis (Draine 2011), the LIS lines are likely saturated at log(N_Si/cm⁻²) > 14, roughly corresponding to log(N_H I/cm⁻²) > 19.7, which can be accomplished for most of our sources but not for all. Therefore, at neutral column densities below log(N_H I/cm⁻²) ∼ 19.7, where H I is still optically thick, some of the LIS lines could be optically thin. For these unsaturated metal lines, Eq. (4) no longer applies, and the covering fraction cannot be derived from the residual flux of the lines.

In any case, the H I residual fluxes can still be roughly predicted from the LIS line residual fluxes using Eq. (9), a trend which is also found in simulations (Mauerhofer et al. 2021). This is important at high z, where the Lyman series is not easily observable due to the increase in IGM opacity, and so the H I residual fluxes are not measurable and the $f_{\mathrm{esc}}^{\mathrm{abs}}$ cannot be obtained from them. It should be finally considered that the relation between H I and LIS residual fluxes (Eq. (9)) may not be universal. However, so far the available data at high z are compatible with our findings.

6.3. The porosity of the ISM favors the escape of LyC photons

As discussed in Sect. 3.2, H I lines are likely saturated which means high neutral gas column densities of log(N_H I/cm⁻²) ≳ 15−16. Gazagnes et al. (2018) measured the H I column density of 6 out of 11 of our Izotov et al. (2016a,b, 2018a,b) leakers using the weak O Iλ1039 line and found column densities of log(N_H I/cm⁻²) ≳ 18, which means that the neutral gas medium is also opaque to the LyC photons. Therefore, a density-bounded leaking scenario as proposed by Zackrisson et al. (2013), Jaskot & Oey (2013), Nakajima & Ouchi (2014), for example, can be excluded in principle. In contrast, the existence of saturated but nonblack H I lines suggests a nonunity H I covering fraction, where the radiation can escape through holes in the ISM, although these channels may not be free of gas and dust, and therefore may not be fully optically thin (Gazagnes et al. 2018; Steidel et al. 2018). Such a configuration can also explain the high ionization parameters and other properties associated with Green Pea galaxies (Jaskot et al. 2019), and is needed to explain the observed emission line ratios of LyC emitters (Ramambason et al. 2020).

The relationship – already seen in previous figures – between the LyC absolute escape fraction ($f_{\mathrm{esc}}^{\mathrm{abs}}$) and the residual flux, R(H I) or R(LIS), is shown in Fig. 11. We performed a survival Kendall correlation analysis in order to take into account the upper limits on the escape fraction and the lower limits on the residual fluxes (see Akritas & Siebert 1996; Flury et al. 2022b). A strong (τ = −0.358) and significant correlation (p_val = 1.428 × 10⁻⁶) between $f_{\mathrm{esc}}^{\mathrm{abs}}$ and R(H I) is found (see Table 3). The strongest leakers tend to have the highest residual fluxes or, conversely, the lowest H I covering fractions. Such behavior is expected from the picket-fence model, where the neutral gas is distributed in clumps around the light sources.

Fig. 11. Correlation between the LyC absolute photon escape fraction ($f_{\mathrm{esc}}^{\mathrm{abs}}$) and the H I (left) and LIS line (right) residual fluxes. Purple circles represent the leakers, while the black triangles are the nondetections – and for which the $f_{\mathrm{esc}}^{\mathrm{abs}}$ is an upper limit. The size of the points corresponds to the relative uncertainty on the escape fraction ($f_{\text{esc}}^{\text{abs}}/f_{\mathrm{esc}}^{\mathrm{abs},\mathrm{err}}$), and the median 1σ uncertainty of the residual flux is placed in the bottom-right corner of each panel. Red squares correspond to the Izotov et al. (2016a,b, 2018a,b, 2021) and Wang et al. (2019) galaxies (see legend). Overplotted are the curves predicted by Eq. (6) (right) and Eq. (12) (left) for EB-V= 0, 0.1, 0.2, and 0.3 mag, respectively (dashed lines). The black solid line in the right panel corresponds to the best linear fit, and the gray-shaded area is the 2σ confidence interval including intrinsic scatter (Eq. (10), see text).

Clearly, H I residual flux overestimates the photon escape in most of the cases: $f_{\text{esc}}^{\text{abs}}\le 1-C_f(\mathrm{H}\mathrm{I})$. This is expected, because the presence of dust needs to be accounted for when inferring the $f_{\mathrm{esc}}^{\mathrm{abs}}$ (Gazagnes et al. 2018; Chisholm et al. 2018; Steidel et al. 2018). In the same figure, we plot the predictions of the Eq. (6) for E_B−V = 0, 0.1, 0.2, and 0.3 (R16 law). Visually, one can see how our data points follow the curves on top given by the picket-fence scenario (uniform dust-screen). At a given R(H I), variations in the LyC photon escape can be explained by different amounts of dust reddening.

The LyC photon escape fraction behaves similarly to R(LIS) ($f_{\text{esc}}^{\text{abs}}$ < 1 − C_f(LIS)), and shows a strong correlation too (Fig. 11, τ = 0.236, p_val. = 1.457 × 10⁻³), where leakers have higher residual LIS flux than nonleakers. The differences between the two plots is principally introduced by the C_f(H I)−C_f(LIS) relation (Sect. 6.2). There are, on the other hand, several possible scenarios in which the escape of ionizing photons could differ from the measured residual flux. First, as soon as the ISM conditions move from optically thick (saturated lines) to optically thin regimes (linear), the optical depths of the gas clouds are close to one, and so the covering fraction is not well-defined. The presence of residual gas and dust in the low column density channels can also contribute to the same effect (see discussion in Gazagnes et al. 2020). Secondly, the distribution of sources emitting ionizing photons and nonionizing UV continuum may not be the same. From detailed radiation transfer simulations, Mauerhofer et al. (2021) found that around 15%–30% of the light at for example 1020 Å is produced by stars that do not contribute to the budget of ionizing radiation. Hence, these continuum photons would be affected by screens of gas and/or dust that do not influence the LyC escape fraction. Thirdly, as already mentioned in the previous section, the gas kinematics and turbulence can spread the absorption line profiles, increasing the residual fluxes and causing the covering fractions to be underestimated (also Mauerhofer et al. 2021). All these reasons would contribute to the scatter and the decreased ability to infer the $f_{\mathrm{esc}}^{\mathrm{abs}}$ of galaxies (see Sect. 7.1) using the H I/LIS lines. More generally, a relation between the LyC escape and the residual fluxes of the lines, interpreted as the escape fractions of the lines in question, can always be expected (Mauerhofer et al. 2021), even when $f_{\mathrm{esc}}^{\mathrm{abs}}$ cannot be analytically predicted via the picket-fence model (see Sect. 7.4).

Given the strong correlation between $f_{\mathrm{esc}}^{\mathrm{abs}}$ and the residual line flux, we can also expect correlations with H I or LIS equivalent widths. These relations are displayed in Fig. 12. The $W_{\text{H}\text{I}}$ is significantly scattered (left panel), mainly driven by the nonleaker population (τ = −0.425, p_val = 1.082 × 10⁻⁸). On the LIS side (right panel), a tight correlation between the $f_{\mathrm{esc}}^{\mathrm{abs}}$ and the W_LIS exists, with τ = −0.331 and a p−value of p_val. = 8.273 × 10⁻⁶.

Fig. 12. The correlation between the LyC absolute photon escape fraction ($f_{\mathrm{esc}}^{\mathrm{abs}}$) and the H I (left) and LIS lines (right) equivalent widths. The sizes and colors of the points have the same meaning as in Fig. 11. The median 1σ uncertainty on the equivalent width is placed at the bottom of every panel. The black solid line in the right panel corresponds to the best linear fit (Eq. (11)), and the gray-shaded area represents the 2σ confidence interval, including intrinsic scatter. The upper limit prediction from Mauerhofer et al. (2021) using the C IIλ1334 line is plotted in orange.

It is of great interest to derive empirical relations using the LIS lines in order to predict the LyC escape fraction of galaxies. Nonetheless, caution is required, and all the previous limitations of the method and the current scatter in the correlations must be considered (see Figs. 11 and 12). For the LIS lines, a Bayesian linear fit to the data in log-linear space was performed using LINMIX (Kelly 2007), which provides a proper treatment of the $f_{\mathrm{esc}}^{\mathrm{abs}}$ upper limits. The resulting relation between the absolute LyC escape fraction and the LIS residual flux is

$$\begin{array}{c}\hfill log{f}_{\mathrm{esc}}^{\mathrm{abs}}=(-4.93\pm 0.65)+(5.57\pm 1.11)\times \mathrm{R}(\mathrm{LIS}),\end{array}$$(10)

where the linear correlation coefficient of the fit is r² = 0.69 (+0.18, −0.21), with an important intrinsic scatter of σ_y = 0.20. For the $f_{\mathrm{esc}}^{\mathrm{abs}}$ to W_LIS relation, and applying once again the LINMIX code to the log-linear data space, we get

$$\begin{array}{c}\hfill log{f}_{\mathrm{esc}}^{\mathrm{abs}}=(-0.35\pm 0.32)-(1.58\pm 0.34)\times W_{\mathrm{LIS}}/\mathrm{\AA },\end{array}$$(11)

with r² = 0.88 (+0.06, −0.09) and σ_y = 0.22. In Fig. 12 (right), our results are compared to the ones derived by Mauerhofer et al. (2021) using the C IIλ1334 line. These latter authors simulated a low-mass high-z galaxy and generated mock C IIλ1334 observations along different sight lines. The C IIλ1334 equivalent width (W_C II) was then computed together with the escape fraction for every line of sight. Mauerhofer et al. (2021) concluded that the equivalent width can be treated as an upper limit to the absolute photon escape as: $f_{\text{esc}}^{\text{abs}}<0.5-0.4\times W_{\mathrm{C}\mathrm{I}\mathrm{I}}$. Compared to our linear fit, we can see that both results are similar, although the Mauerhofer et al. (2021) relation follows a slightly different slope. In the latter work, the absorption line measurements were carried out in a similar way to that described in the present paper but only for a single line. Also, the shift could be due to the modeling of the lines, the dust location within the simulated galaxy (the density of dust is proportional to the metallicity located in the highest density H I clouds), or the different intrinsic properties (SFR, ionizing emissivity, metallicity, etc.) between our low-z sample and the z ∼ 3 simulated galaxy. Kendall survival statistics quantifying the correlations of $f_{\mathrm{esc}}^{\mathrm{abs}}$ versus H I/LIS equivalent widths and residual fluxes are computed and presented in Table 3.

6.4. Dust attenuation and LyC leakage

The impact of dust extinction on the escape of ionizing photons is not well understood. Based on simulations, Kimm et al. (2019) and Mauerhofer et al. (2021) suggested that the absorption due to dust is subdominant compared to that due to neutral hydrogen. Physically, LyC photons are destroyed in the same way as nonionizing photons at slightly longer wavelengths, but if the dust resides only in regions of high neutral gas density (usually assumed in the simulations), the high ionization cross-section and large abundance of the hydrogen atom is such that H I likely dominates the absorption of LyC photons in galaxies more than dust does. However, if the dust is still present in the ionized regions or channels (as in the uniform dust-screen scenario), LyC photons would significantly attenuate before escaping the galaxy.

Chisholm et al. (2018) observationally demonstrated that dust needs to be accounted for in order to infer the absolute escape from the Lyman series. It is essential to correct for the UV attenuation in order to derive the intrinsic LyC photon production, which yields the absolute escape fraction by comparison to the observed LyC flux. Therefore, the “effective” UV attenuation A_UV enters any derivation of the absolute escape fraction, whereas the relative escape fraction only accounts for the differential attenuation between 912 Å and nonionizing flux (see e.g., Siana et al. 2015; Steidel et al. 2018).

To investigate the role of the dust in the escape fraction in the LzLCS sample, in Fig. 13 we plot $f_{\mathrm{esc}}^{\mathrm{abs}}$ against the dust attenuation at 912 Å: A₉₁₂ = k₉₁₂ × E_B-V (R16). A strong correlation of τ = −0.367 is found between the two variables (p_val = 7.624 × 10⁻⁷). Most of the sources are below the 10^−0.4A_UV curve (i.e., $f_{\text{esc}}^{\text{abs}}\le {10}^{-0.4A_{\text{UV}}},C_f(\mathrm{H}\mathrm{I})=0$ following Eq. (6)) which indicates, as mentioned in the previous section, that the covering fraction and dust both overestimate the LyC photon escape of galaxies when considered independently (Heckman et al. 2001, 2011; Chisholm et al. 2018; Gazagnes et al. 2020).

Fig. 13. Correlation between the LyC absolute photon escape fraction ($f_{\mathrm{esc}}^{\mathrm{abs}}$) and the UV dust-attenuation: A912 = k912 × EB-V. Labeling is as in Fig. 11. The median 1σ uncertainty on AUV is plotted on the right of the plot. Overplotted are the three lines predicted by Eq. (6) for Cf(H I) = 0.5,0.75, and 0.95, respectively (dashed lines), using k912 = 12.87 (R16 law).

In Fig. 13, the results given by Eq. (6) for C_f(H I) = 0.5,0.75, and 0.95 are also shown. Clearly, the spread in the $f_{\mathrm{esc}}^{\mathrm{abs}}$-dust plane can be explained by joint changes of E_B−V and the H I covering fraction of galaxies. We use this expression below to predict the $f_{\mathrm{esc}}^{\mathrm{abs}}$ values of the LzLCS and high-z galaxies. Simultaneously considering the neutral gas covering (given by the depths of the Lyman lines) and the dust attenuation (from the UV spectral fits) provides the most accurate photon escape prediction of this work.

We note that a correlation between $f_{\mathrm{esc}}^{\mathrm{abs}}$ and E_B−V is expected because the $f_{\mathrm{esc}}^{\mathrm{abs}}$ values come from de-reddening the S99 models (see Eqs. (1) and (2)) and comparing the predicted LyC flux to the observed one. Ideally, an independent measurement of the UV attenuation, which is here obtained from fitting the UV spectra, would test if this correlation is indeed as strong as shown here. We also note that the exact value of the UV attenuation (or E_B−V) depends in principle on the assumed geometry and the attenuation law, for which we have only considered two simple cases (see Sect. 4). Probably the most direct verification of the robustness of these relations comes from the comparison of the absolute escape fractions derived using two independent methods, namely the one presented here, which uses the stellar UV spectrum, and the method used by Flury et al. (2022a) (see also Izotov et al. 2016b). In the latter study, $f_{\mathrm{esc}}^{\mathrm{abs}}$ is computed from the observed LyC flux and the extinction-corrected Hβ flux, together with Hβ equivalent width which is used as an indicator of the age of the underlying ionizing stellar population. As shown in Flury et al. (2022a), the overall agreement between the two methods is quite good, verifying the strong correlation that we see here between the escape and the dust-attenuation.

Finally, it is worth mentioning that the correlations between the absolute LyC escape fraction, residual fluxes, equivalent widths, and the UV attenuation shown in Figs. 11–13 also hold when one considers the observed flux ratio f_LyC/f₁₁₀₀ instead of $f_{\mathrm{esc}}^{\mathrm{abs}}$. In particular, the fact that both f_LyC/f₁₁₀₀ and $f_{\mathrm{esc}}^{\mathrm{abs}}$ correlate with the UV attenuation indicates that dust is primarily in the clouds, as shown in Sect. 5.2 (see also Sect. 7.4), and hence the amount of attenuation increases with the covering fraction.

For easier comparison with future studies, and to better illustrate the behavior of the LyC escape with the previous parameters, a summary plot of Figs. 11–13 can be seen in Appendix E (Fig. E.1). To this end, we computed the median value of R(H I), R(LIS), W_H I, W_LIS and E_B−V in three $f_{\mathrm{esc}}^{\mathrm{abs}}$ intervals, specifically at $f_{\text{esc}}^{\text{abs}}$ ∈ [0 − 5%), [5 − 10%), and [10−100%). Table E.1 emphasizes these results, and a quantitative analysis of these trends, comprising the sample median values, will be further investigated in future LzLCS publications using stacked spectra.

7. Discussion

Having established relations between the UV absorption lines and LyC escape, we now examine the power of these lines to predict the ionizing photon escape fraction of galaxies (Sect. 7.1). We also present a physical picture of the ISM of low-z LyC emitters based on our UV absorption lines and $f_{\mathrm{esc}}^{\mathrm{abs}}$ results as well as earlier results in the literature (Sect. 7.2). We finally discuss some implications and comparisons for the high-z Universe (Sect. 7.3), and the caveats of our modeling approach (Sect. 7.4).

7.1. Predicting the Lyman-continuum photon escape with the Lyman series and LIS lines

We want to test whether the H I and, more particularly, measurements of the LIS lines could be used to effectively predict the $f_{\mathrm{esc}}^{\mathrm{abs}}$ of galaxies. Under the assumption of a foreground screen of dust in the picket-fence model, Chisholm et al. (2018) determined that the parameters that contribute the most to predicting the photon escape are the dust attenuation and the H I covering fraction, implying that the column density negligibly affects the line depths (if the lines are saturated). Therefore, the accurate determination of both E_B−V and C_f(H I) is necessary. In this section, we carry out a similar exercise as in Chisholm et al. (2018) and Gazagnes et al. (2020), but taking advantage of our larger LzLCS data set, and thus improving the robustness of the method and providing more statistical meaning. First, we take the derived E_B−V from our SED fits as well as the weighted-averaged H I covering fractions, C_f(H I). The predicted photon escape from the H I lines is computed from Eq. (6), $f_{\mathrm{esc}}^{\mathrm{abs},\mathrm{H}\mathrm{I}}$.

Thus, the predicted $f_{\mathrm{esc}}^{\mathrm{abs},\mathrm{H}\mathrm{I}}$ value is compared to the fiducial value derived from the spectral UV fits (Sect. 3.3) in Fig. 14 (top). LINMIX determines a clear correlation between the two, with a correlation coefficient of r² = 0.97 (+0.02, −0.07). The intrinsic scatter is negligible (σ_y ∼ 10⁻⁴), which means that the errors on the variables can fully explain the scatter in the 1:1 relation. The median of the predicted $f_{\mathrm{esc}}^{\mathrm{abs}}$ for leakers does not change with respect to the fiducial $f_{\mathrm{esc}}^{\mathrm{abs}}$ distribution, although a small systematic shift to higher predicted escape values exists for the $f_{\text{esc}}^{\text{abs}}$ < 10% distribution, where a significant fraction of the sources are upper limits in $f_{\mathrm{esc}}^{\mathrm{abs}}$. In conclusion, $f_{\mathrm{esc}}^{\mathrm{abs}}$ can be well recovered from the Lyman series and dust attenuation.

Fig. 14. Comparison of predictive methods for the LyC absolute escape fraction. Top: Lyman series photon escape (Eq. (6)) against the fiducial S99-derived values. Bottom: same as top panel but for the photon escape inferred using the LIS line covering fraction (Eq. (9)). Purple points correspond to the LzLCS leakers while the red squares are the Izotov et al. (2016a,b, 2018a,b, 2021) and Wang et al. (2019) galaxies. Black triangles correspond to nonleakers. Median relative error bars are shown in the upper right corners of each panel. The subplots at the bottom illustrate the overall accuracy of the predicted $f_{\mathrm{esc}}^{\mathrm{abs}}$ (i.e., Δ$f_{\text{esc}}^{\text{abs}}$/(1 + $f_{\text{esc}}^{\text{abs}}$)) as a function of the fiducial $f_{\mathrm{esc}}^{\mathrm{abs}}$.

In any case, our original purpose is to decipher whether or not the $f_{\mathrm{esc}}^{\mathrm{abs}}$ can be predicted using the LIS lines. To approach this, we first take advantage of the intrinsic relation between the residual flux (or covering factor) of H I and metals, which is empirically given by our Eq. (9). The H I is then computed by substituting every measured C_f(LIS) in that expression, and then the predicted photon escape from LIS, $f_{\mathrm{esc}}^{\mathrm{abs},\mathrm{LIS}}$, is determined with a formula similar to that in Eq. (6):

$$\begin{array}{c}\hfill f_{\mathrm{esc}}^{\mathrm{abs},\mathrm{LIS}}={10}^{-0.4k_{912}{E}_{\mathrm{B}-\mathrm{V}}}\times (1-[a\times C_f(\mathrm{LIS})+b]),\end{array}$$(12)

where the coefficients [a, b] = [(0.63 ± 0.19),(0.54 ± 0.09)] are given by Eq. (9). The results are plotted against the fiducial $f_{\mathrm{esc}}^{\mathrm{abs}}$ in Fig. 14 (bottom). The gray arrows point towards $f_{\mathrm{esc}}^{\mathrm{abs},\mathrm{LIS}}\approx 0$, for those sources with C_f(LIS) ≈ 1 (zero residual flux) and/or very high E_B−V, falling below the axis limits. Still, a strong linear correlation is present when using the LIS lines (r² = 0.76 (+0.20, −0.32)), with no significant intrinsic scatter.

As discussed earlier (Sect. 6.3), several mechanisms can lead to disagreement between the real covering fraction and the covering derived using the residual flux of the lines: ISM regions close to the optically thin limit, the mismatch between the distribution of UV ionizing and nonionizing sources, the gas kinematics and turbulence, and the ISM geometry can all potentially contribute to the scatter in the relations shown in Fig. 14 (see Mauerhofer et al. 2021).

A more detailed look into Fig. 14 reveals that for the highest $f_{\text{esc}}^{\text{abs}}$ ≥ 0.1 sources, the Lyman series and low-ionization covering methods systematically underestimate the photon escape. The effect is more noticeable when one plots the overall accuracy of the predicted $f_{\mathrm{esc}}^{\mathrm{abs}}$ (Δ$f_{\text{esc}}^{\text{abs}}$/(1 + $f_{\text{esc}}^{\text{abs}}$)) as a function of the fiducial $f_{\mathrm{esc}}^{\mathrm{abs}}$ itself (Fig. 14, bottom panels). This result, which was already noted by Gazagnes et al. (2020) and is supported by hydrodynamical simulations by Kakiichi & Gronke (2021), suggests that the leakage could be partially due to a density-bounded ISM, where the LyC photons can escape through optically thin log(N_H I)cm⁻² ≲ 10^17.2 regions to the LyC.

Overall, we propose three different approaches to predict the LyC absolute photon escape from the LIS lines, and demonstrate the limitations on each: (1) using the residual flux of the absorption lines (Eq. (10)), (2) their equivalent widths (Eq. (11)), and (3) a combination of residual flux and dust attenuation (Eq. (12)). When compared directly to $f_{\mathrm{esc}}^{\mathrm{abs}}$, the third method is ideally recommended in all the cases, as it presents the most stringent correlation coefficient of the fit and the least intrinsic scatter.

A similar methodology to (3) allowed Chisholm et al. (2018) to predict the $f_{\mathrm{esc}}^{\mathrm{abs}}$ of five high-z galaxies detected in LyC, and without Lyman-series observations. In particular, for the Cosmic Horseshoe (Belokurov et al. 2007) at z = 2.38, these authors predicted a $f_{\mathrm{esc}}^{\mathrm{abs}}$ of around 1%, in agreement with the $f_{\text{esc}}^{\text{abs}}$ < 2% limit imposed by HST imaging (Vasei et al. 2016).

At high z, low-mass galaxies are characterized by hosting a metal- and dust-poor ISM similar to our low-z analogs, where E_B-V < 0.1 at z ∼ 3 − 6 (Schaerer et al. 2013). Taking this attenuation as representative, and using the median H I and metal LIS covering fraction for the leakers in our sample (0.70 and 0.38, respectively), this gives $f_{\text{esc}}^{\text{abs}}$ ≥ 0.07 for H I and LIS predicted values. This is below but similar to the 10% theoretical lower limit for faint SF galaxies to dominate reionization at z ∼ 6 − 9 (Robertson et al. 2013)⁸, and is in good agreement with the 9% (6%) average value derived in Steidel et al. (2018), Pahl et al. (2021) at the same redshifts; and also Rosdahl et al. (2018) in simulations. However, it should also be stressed that the predictions for $f_{\mathrm{esc}}^{\mathrm{abs}}$ given by Eq. (12) are only reliable in a statistical sense. Given the scatter in the above-mentioned relations, predictions for any individual object are not very reliable.

We finally point out that higher resolution and higher S/N UV spectra would help to obtain a better measurement of the covering fractions and ionic column densities (fundamental to disentangling the role of the ISM structure and gas kinematics in LyC leakage), which would improve the $f_{\mathrm{esc}}^{\mathrm{abs}}$ determination from the SED modeling (subject to quite large uncertainties in some cases), and to better constrain the stellar population and dust properties (E_B−V).

7.2. The absorbing ISM of LyC emitters

Low-metallicity (subsolar) galaxies typically have lower masses (≲10¹⁰ M_⊙), higher gas fractions, and less dust than the typical local SF galaxy (Trebitsch et al. 2017). Low-mass systems have lower gravitational potentials, which increases the efficiency of stellar feedback and leads to more bursty star-formation episodes (Muratov et al. 2015). The low-metallicity and burstiness of their star-formation history means that these sources also have high ionizing production efficiencies (intrinsic number of ionizing photons per UV luminosity or mass, equivalently) and moderate ionizing escape fractions (Ma et al. 2020). This, together with their higher number density at all epochs (with respect to more massive galaxies), makes them the most promising analogs of the sources of cosmic reionization (Finkelstein et al. 2019).

Different mechanisms can remove the gas and dust from such galaxies, creating holes in the ISM with appreciably lower gas and dust column densities. Our observations suggest that LyC photons primarily escape through these channels. These holes are likely created by stellar winds and supernova explosions after starburst episodes, or they are photoionized channels induced by turbulence during the starburst period (Kakiichi & Gronke 2021). They could also be embedded in an intrinsically low-density ISM, which results from catastrophic cooling effects around low-metallicity star clusters (Jaskot et al. 2019).

The picket-fence model describes an ionization-bounded galaxy (Zackrisson et al. 2013; Vasei et al. 2016; Reddy et al. 2016b) where the ionizing photons can partially escape through these holes in the ISM. The fraction of sight lines covered by optically thick neutral gas clouds is the covering fraction. We now summarize a scenario where the H I covering fraction, C_f(H I), can explain most of the observed empirical trends found in the present study.

The presence of (likely) saturated but nonblack H I and LIS absorption lines in our spectra (Sect. 5.2) implies the existence of low-column-density channels, and therefore a nonzero covering fraction for most of the galaxies. Measurements of the H I column density of Green-Pea galaxies lead to similar conclusions in Gazagnes et al. (2018) and McKinney et al. (2019), and in the recent, detailed study of the local galaxy Haro11 by Östlin et al. (2021).

The LzLCS sample reveals many strong correlations between the UV absorption lines and other secondary galaxy parameters. We find that the measured H I and LIS equivalent widths scale with (1) the intrinsic absolute magnitudes at 1500 Å, $M_{1500}^{\mathrm{int}}$, and with (2) the UV attenuation, A_UV (or E_B−V, equivalently), and scale inversely with (3) the H I and LIS residual fluxes (Figs. 6 and 7). We also demonstrated that metallicity plays a minor role in these relations.

The inverse connection between the equivalent width of the lines and the residual flux (3) naturally comes from the definition of a saturated absorption line profile (Draine 2011). In contrast, the W_λ − A_UV relation (2) can be explained by invoking a clumpy gas-to-dust geometry (Shapley et al. 2003; Steidel et al. 2018; Du et al. 2018, 2021; Pahl et al. 2020), where the total reddening is related to the fraction of the total sight lines with optically thick neutral gas (Vasei et al. 2016; Reddy et al. 2016b). A galaxy with a larger fraction of sight lines covered by optically thick neutral gas (a higher C_f(H I)) also encounters more dust and the spectrum is more attenuated (A_UV increases), because the dust resides in the same gas. In other words, the same gas which gives rise to W_LIS also reddens the continuum (Jones et al. 2012).

Aside from the W_λ − A_UV connection, there are other lines of evidences for a clumpy gas-to-dust geometry in our sample. For instance, in a clumpy geometry, transitions for the same ion at different wavelengths will have different residual fluxes (see Eq. (5) and Gazagnes et al. 2018). We tentatively see this trend when looking at the median residual flux of every Lyman line in the LzLCS. Values vary by as much as 0.71, 0.81, 0.83, and 0.89 from Ly5 to Lyβ, but typical errors of ∼0.13 obscure robust conclusions. The hypothesis of a clumpy distribution of gas and dust has been widely proposed in the literature, and it has special relevance for example when trying to explain the observed Lyα-to-Hα ratios in SF galaxies (see Scarlata et al. 2009; Jaskot et al. 2019, and references therein), and the Lyα peak-velocity separation (Gronke et al. 2016; Orlitová et al. 2018; Jaskot et al. 2019).

When Lyα varies from weak absorption to strong emission, the strengths of the H I and LIS absorption lines decrease (Fig. 9 in Sect. 6.1, see also Shapley et al. 2003; Jones et al. 2012; Rivera-Thorsen et al. 2015; Du et al. 2018, 2021; Trainor et al. 2019; Pahl et al. 2020). Invoking our picket-fence model, the higher the covering fraction, the higher the fraction of Lyα photons that are resonantly scattered out of the line of sight and absorbed by dust, resulting in a weaker Lyα emission with stronger H I and LIS profiles (Steidel et al. 2018; Gazagnes et al. 2020). Similarly, using z ∼ 3 − 5 composites, Jones et al. (2012) demonstrated that other host properties, such as the luminosity, have little effect on the W_Lyα– W_LIS variations. The relationships between W_Lyα–W_LIS and W_Lyα–A_UV are presumably invariant across cosmic time, as both are affected by variations in the same H I covering (Du et al. 2018). This also suggests that the Lyα production efficiency does not change drastically with redshift (Trainor et al. 2019; Pahl et al. 2020).

Moreover, the existence of double- or triple-peaked Lyα profiles in leakers (Verhamme et al. 2017; Izotov et al. 2020, 2021) can be associated with the characteristic bimodal distribution of H I column densities in the patchy ISM (Verhamme et al. 2015; Mauerhofer et al. 2021; Kakiichi & Gronke 2021), or the presence of completely ionized channels (Rivera-Thorsen et al. 2017). This may lead to the observed correlations between the escape of Lyα and LyC photons (e.g., Verhamme et al. 2017; Marchi et al. 2018; Steidel et al. 2018; Gazagnes et al. 2020; Pahl et al. 2021).

While the picket-fence model allows us to predict (at least on average) the LyC escape fraction relatively well, we find that the absorption lines tend to underestimate the escape fractions for the strongest leakers with $f_{\mathrm{esc}}^{\mathrm{abs}}$ ≳0.1 − 0.2 (see also Gazagnes et al. 2020). This suggests that there could be a secondary scenario where LyC emission is regulated by a density-bounded ISM in these galaxies, as suggested by several previous studies (Zackrisson et al. 2013; Jaskot & Oey 2013; Nakajima & Ouchi 2014; Jaskot et al. 2019). Indeed, from simulations, Kakiichi & Gronke (2021) show that a switch from an ionization-bounded regime to one that is density-bounded may naturally happen at evolved stages in the history of a galaxy. This transition occurs when turbulence ceases and the ionization radius then increases with time in a more static situation. Indications for the existence of the density-bounded regime for leakers with the highest $f_{\mathrm{esc}}^{\mathrm{abs}}$ was presented by Gazagnes et al. (2020) and Ramambason et al. (2020), but whether this also quantitatively agrees with the present observations will be examined in future publications.

7.3. The LyC escape fraction of high-z galaxies

With the empirically calibrated and theoretically motivated relations between the residual flux, UV attenuation, and the absolute LyC escape fraction (Eq. (12)), we can estimate $f_{\mathrm{esc}}^{\mathrm{abs}}$ from literature measurements of the LIS absorption lines for high-z galaxies, and compare them to local measurements from the LzLCS.

Figure 15 presents the result of this exercise, showing the predicted $f_{\mathrm{esc}}^{\mathrm{abs}}$ from individual and stacked galaxy spectra reporting residual flux measurements of the UV LIS lines as a function of the observed absolute UV magnitude (the latter without a correction for UV attenuation). To compute the LyC escape fraction, we use the reported LIS residual flux, and then Eq. (9) to estimate the C_f(H I), and the R16 attenuation law to estimate the E_B−V from the average dust corrections given in Bouwens et al. (2020), as a function of redshift⁹. We also show the $f_{\mathrm{esc}}^{\mathrm{abs}}$ measurements at low z from the present paper (gray, red circles and downward point triangles), and the latest update by Pahl et al. (2021) of the LyC escape fractions measured directly from stacks of z ∼ 3 galaxy spectra from the Keck Lyman Continuum Survey (KLCS, black circles; Steidel et al. 2018). Downward arrows means predicted $f_{\text{esc}}^{\text{abs}}$ ≈ 0, for those sources with C_f(LIS) ≈ 1 (zero residual flux) and/or very high E_B−V, falling below the axis limits.

Fig. 15. Estimated and measured LyC absolute escape fractions of high- and low-z galaxies as a function of the observed $M_{1500}^{\mathrm{obs}}$. Open circles and squares show the $f_{\mathrm{esc}}^{\mathrm{abs}}$ predictions from the UV LIS lines for individual galaxies and stacks, respectively, taken from Harikane et al. (2020) at z ∼ 6, in pink; Sugahara et al. (2019) at z ∼ 5, in dark red; Leethochawalit et al. (2016) at z ∼ 4 − 5, in blue; Jones et al. (2013) at z ∼ 2 − 4, in light green; and Jones et al. (2012) at z ∼ 4, in green. Direct $f_{\mathrm{esc}}^{\mathrm{abs}}$ measurements from this paper (LzLCS and literature) are plotted in light gray and red; z ∼ 3 measurements from Pahl et al. (2021) measurements are displayed with filled black symbols, for comparison. Downward-pointing arrows represent sources with estimated $f_{\text{esc}}^{\text{abs}}$ ≈ 0.

As seen in Fig. 15, all our estimates yield absolute LyC escape fractions $f_{\text{esc}}^{\text{abs}}$ < 0.1 at z ≳ 4 for relatively bright $M_{\text{1500}}^{\text{obs}}$ ≲ −21 galaxies. Measurements from individual sources show large scatter, as also found at low z. For the highest redshift available (z ∼ 6), the stack and the galaxy J0210–0523 of Harikane et al. (2020), with $M_{\text{1500}}^{\text{obs}}$ ∼ −23, and our method predict negligible escape fractions, $f_{\text{esc}}^{\text{abs}}$ ≲ 0.003. From the z ∼ 4 stack of Jones et al. (2012) with $M_{\text{1500}}^{\text{obs}}$ ∼ −21, we estimate $f_{\text{esc}}^{\text{abs}}$ ∼ 0.05, which is comparable to the value measured at z ∼ 3 for galaxies with similar UV magnitude (Pahl et al. 2021). We reiterate the fact that our $f_{\mathrm{esc}}^{\mathrm{abs}}$ values of the high-z galaxies are lower than those estimated in Jones et al. (2012), Harikane et al. (2020), because we account for both the residual flux and the UV attenuation (Gazagnes et al. 2018; Chisholm et al. 2018). Escape fractions of ∼10% are found at $M_{\text{1500}}^{\text{obs}}$ ∼ −20 from the KLCS stacks, and Pahl et al. (2021) find decreasing $f_{\mathrm{esc}}^{\mathrm{abs}}$ towards UV-brighter galaxies, which overall appears consistent with our estimates. However, the available data are still scarce, and whether $f_{\mathrm{esc}}^{\mathrm{abs}}$ increases further in fainter galaxies must be confirmed with future observations (see Bian & Fan 2020).

The overlap between the high- and low-z data in observed UV magnitude is quite limited (Fig. 15) and meaningful sample averages remain to be obtained at low redshift. It is therefore difficult to draw conclusions on the possible redshift evolution of $f_{\mathrm{esc}}^{\mathrm{abs}}$ from this figure, because possible trends may be a consequence of a general luminosity evolution or due to different detection limits within the high-z samples. Furthermore, the adopted UV attenuations for the high-redshift samples are unconstrained, making the comparison of $f_{\mathrm{esc}}^{\mathrm{abs}}$ as a function of the intrinsic UV magnitude also difficult to interpret. In any case, we note that both individual galaxies and stacks with $f_{\text{esc}}^{\text{abs}}$ ∼ 5%−10% are found over a range of four magnitudes (from $M_{\text{1500}}^{\text{obs}}$ ∼ −19 to −22). Objects with $f_{\text{esc}}^{\text{abs}}$ > 0.1 are (so far) found only at magnitudes fainter than $M_{\text{1500}}^{\text{obs}}$ < −20 at low z. However, some individual galaxies at z > 3 are bright ($M_{\text{1500}}^{\text{obs}}$ < −21) and have high LyC escape fractions (e.g., Ion2 in Vanzella et al. 2016). Thus, bright high-redshift galaxies may have markedly different LyC escape properties. Overall, the low-z results concur, albeit with the low LyC escape fractions found in bright galaxies ($M_{1500}^{\mathrm{obs}}$≲ − 20) by many studies using stacks at high z (see e.g., Grazian et al. 2016; Guaita et al. 2016; Marchi et al. 2017; Japelj et al. 2017).

Another quantity to be compared with high-z studies is the observed “Lyman decrement” between the LyC and the nonionizing UV flux, that is, the commonly used (f _LyC/f₁₅₀₀) _out corrected for the IGM attenuation. This quantity can directly lead to the ionizing emissivity of galaxies when integrating over the UV luminosity function (see e.g., Pahl et al. 2021). In Fig. 16 we compile the available data from high-z stacks, for which we plot the IGM-corrected measurements of (f _LyC/f₁₅₀₀) _out. We use stacked spectra to compensate for differences in the mean IGM transmission. The data comprise results from a variety of surveys mostly carried out at redshift z ∼ 2.5 − 4, plus the study of Alavi et al. (2020) at z ∼ 1.3, and Bridge et al. (2010) at z ∼ 0.7. We have also added the data from the composite SEDs of the LAE sample of Fletcher et al. (2019) at z = 3.1, who report LyC-detected and nondetected subsamples.

Fig. 16. Ionizing-to-nonionizing observed flux ratio (f LyC/f1500) out in units of Fν as a function of $M_{1500}^{\mathrm{obs}}$ for high- and low-z galaxies. At high redshift, we only show measurements determined from stacks. LyC detections from Pahl et al. (2021) and Fletcher et al. (2019) at z ∼ 3 are plotted in black and blue symbols. Circles show detections and downward pointing triangles are 1σ upper limits that are color coded by sample (Pahl et al. 2021; Alavi et al. 2020; Bian & Fan 2020; Fletcher et al. 2019; Naidu et al. 2018; Grazian et al. 2017; Rutkowski et al. 2017; Japelj et al. 2017). The high-z data probe z ∼ 2.5 − 4, except for the study of Alavi et al. (2020) at z ∼ 1.3, and Bridge et al. (2010) at z ∼ 0.7. The low-z data from the present study (LzLCS and literature) are plotted in light gray and red, as in the previous figure.

Interestingly, the Lyman decrement observed both from the KLCS (Steidel et al. 2018; Pahl et al. 2021) and the LAE sample of Fletcher et al. (2019) at z ∼ 3 are very similar, and comparable to the observed values (f _LyC/f₁₅₀₀) _out ∼ 0.1 at z ∼ 0.3 and at same UV magnitudes, $M_{\text{1500}}^{\text{obs}}$ ∼ −19.5 to −21. The upper limits found for galaxies with $M_{\text{1500}}^{\text{int}}$ ≲ −21 indicate small Lyman decrements that are compatible with low-redshift observations. On the UV-faint end ($M_{\text{1500}}^{\text{obs}}$ ≳ −20), the high- and low-z observations also appear fairly consistent with each other¹⁰. Further studies are required to better understand what determines the differences found between the various high- and low-z samples and to determine representative population averages among low-redshift galaxies. In summary, the high-z observations show a consensus of decreasing Lyman decrement towards brighter UV magnitudes, although the overall correlation is weak. This trend is also tentatively shared by the LzLCS sample.

7.4. Main assumptions, limitations, and caveats

Among the possible sources of uncertainty in our analysis are (1) the choice of the dust-attenuation law, (2) the assumed gas and dust geometry, (3) the adopted picket-fence model, and (4) the definition and meaning of the LyC escape fraction.

(1) Assumptions on the dust-attenuation law. As previously described (Sect. 4), the choice of the attenuation curve primarily impacts the derived stellar population parameters and UV attenuation. In particular, X_i coefficients (see Eq. (1)) and the dust-attenuation parameters (E_B−V) can both change, while average stellar ages and metallicities remain consistent (same light-fractions, see Chisholm et al. 2019).

A steeper attenuation law (like the SMC curve) yields a lower UV attenuation and hence systematically higher escape fractions $f_{\mathrm{esc}}^{\mathrm{abs}}$ compared to the flatter R16 attenuation law adopted as default here. Typically, $f_{\mathrm{esc}}^{\mathrm{abs}}$ is found to be higher by a factor of between 1.05 and 1.50 from high to low LyC escape fractions (see Fig. 4). The intrinsic, that is, extinction-corrected UV magnitudes ($M_{1500}^{\mathrm{int}}$) are therefore also fainter for the SMC law. On the other hand, line equivalent widths and residual fluxes are essentially independent of the attenuation law, because the stellar continuum level is similarly modeled for the different laws.

From the present UV data alone, it is difficult to distinguish between the two attenuation laws used here. A joint analysis of the full SED from the UV to the optical, including an independent attenuation measurement such as the Balmer decrement, would be necessary to constrain the attenuation law (see e.g., Shivaei et al. 2020). The offsets found between the color excess E_B−V derived from the UV spectra and the Balmer decrement shown in Fig. 3 indicate that an attenuation law steeper than the R16 curve is probably appropriate for a fraction of the galaxies in our sample. For the compact, high-O₃₂ LyC emitters detected prior to the LzLCS, an SMC-like or even steeper curve in some cases had been found by Izotov et al. (2016b, 2018a,b).

Also, from the analysis of 218 z = 1.4 − 2.6 SF galaxies with both UV rest-frame and Balmer decrement measurements plus known gas-phase metallicities, Shivaei et al. (2020) suggest that “low-metallicity” galaxies have a steeper, SMC-like attenuation curve, whereas galaxies with 12 + log(O/H) ≳ 8.5 are better described by the R16 law. In any case, no consensus exists so far as to the most appropriate UV attenuation law, or even as to the number of dust components that should be considered (Charlot & Fall 2000), and significant variations may exist between individual galaxies, different galaxy types, and so on.

We primarily use the R16 as the default here because it has been observationally constrained down to the Lyman limit, and thus avoids the need for extrapolations of the SMC law to the spectral range covered by the HST/COS spectra used here. Furthermore the R16 law is often used for high-redshift studies (e.g., Reddy et al. 2016b; Steidel et al. 2018), and may thus simplify comparisons.

(2) Assumptions on the gas and dust geometry. The so-called picket-fence model (Sect. 3) proposes two simple ways to describe the relative distribution of dust with respect to the gas in the ISM. Either the dust is distributed in a homogeneous foreground layer (uniform dust screen, A), or it resides within the same clouds as the gas (clumpy geometry, B).

Both UV attenuation and the geometrical covering fraction depend on the adopted geometry, as amply discussed in Vasei et al. (2016), Reddy et al. (2016b), Gazagnes et al. (2018), Steidel et al. (2018). Overall, higher values of E_B−V and C_f are found in the clumpy geometry (geometry B, referred to as the holes scenario by Steidel et al. 2018), as seen for example in Gazagnes et al. (2018), Chisholm et al. (2018), Steidel et al. (2018). This is simply explained by the fact that, in this case, a fraction (1 − C_f) of UV continuum light is directly transmitted, that is, unaltered, to the observer. To fit the observed UV slope and the depth of the absorption lines, higher E_B−V and C_f are therefore needed.

Despite the differences for E_B−V and C_f, Chisholm et al. (2018) and Gazagnes et al. (2020) showed that, to first order, the determination of $f_{\mathrm{esc}}^{\mathrm{abs}}$ from the UV absorption lines is independent of the assumed dust distribution between these two picket-fence geometries. However, for reasons that are not completely clear, Steidel et al. (2018) obtain somewhat lower average absolute LyC escape fractions when assuming the “holes” geometry.

Our analysis of the H I and low-ionization UV absorption lines in Sect. 5.2 clearly shows that the observed behavior favors a clumpy dust geometry, i.e., geometry B, as also suggested by earlier studies (see e.g., Du et al. 2021). However, for our spectral fits, we assume the foreground dust-screen scenario (A), which is not the favored scenario. This assumption is primarily driven by practical reasons and simplicity, because the fits of the stellar continuum (where ISM absorption lines are masked) can then be done independently of C_f, and the absorption lines can subsequently be analyzed. Indeed, in this geometry, the covering fraction is simply related to the residual flux of the (saturated) lines via R = 1 − C_f (Gazagnes et al. 2018).

As long as the two assumptions yield almost identical LyC escape fractions, this apparent inconsistency is not relevant. However, it must be borne in mind that values such as the UV attenuation (quantified by E_B−V) and the covering fraction depend on geometrical assumptions, and they should be handled with care in comparisons between different studies.

(3) Validity of the picket-fence model. In this paper, we make use of the picket-fence model to derive simple analytical solutions of the radiative transfer equations, which can be used to fit the observed galaxy spectra using SEDs of integrated stellar populations and ISM absorption lines (cf., Sect. 3, Vasei et al. 2016; Reddy et al. 2016b). With the assumptions that the “channels” between the pickets have a negligible column density and the absorption lines are saturated in the ISM clouds (pickets), this model then allows one to determine the geometric covering fraction of the absorbing gas from the residual flux of these lines.

Although this simple model works well to describe the observations and predict the average LyC continuum escape fractions, as shown in this paper and earlier studies (Gazagnes et al. 2018; Steidel et al. 2018; Chisholm et al. 2018), we remind the reader that this model is a drastic over-simplification of the complexity observed in resolved nearby galaxies and predicted by state-of-the-art high-resolution galaxy simulations. In realistic situations, the geometry may differ from the simple assumptions of a single UV-emitting source behind a picket-fence ISM and a uniform dust screen made here. Different star-forming regions, whose emission contributes to the integrated spectra, may be “covered” by varying column densities of gas and dust, and each of those may also show substructure. Furthermore, the kinematics of the gas may be complex. How important these effects are for the diagnostics used in the present paper is nevertheless impossible to estimate. In addition, some galaxies, for example, the most compact ones, may also be better approximated by a picket-fence model than others.

The recent work of Mauerhofer et al. (2021), who simulate a z ∼ 3 low-mass galaxy at high resolution, and performed mock observations of UV absorption lines, may provide some guidance towards answering these questions. Indeed, these simulations illustrate the nontrivial definition of the covering fraction and the complex relations between this quantity, the residual line fluxes, and the escape fraction of LyC photons. Despite this, the simulations reveal correlations between residual fluxes and LyC escape, also with significant scatter, similar to those found in our analysis. In any case, going beyond the analytical picket-fence model for the analysis of spectroscopic observations will require the development of new methods and data of higher quality (in S/N and resolution).

(4) Definition and meaning of the LyC escape fraction. As described in Sect. 3, the LyC escape fraction $f_{\mathrm{esc}}^{\mathrm{abs}}$ derived in this paper relies on comparing the observed LyC flux in a limited window shortward of the Lyman edge to the expected intrinsic LyC flux. In our analysis, we measure the LyC typically from 850 to 900 Å in the COS spectra (Flury et al. 2022a), and the intrinsic LyC flux is predicted by the fits to the nonionizing UV spectrum. This method is the most basic and most commonly used method for “direct” determination of the LyC escape fraction, that is, when a direct measurement of the LyC flux is available. It has been used by almost all existing studies of this kind, allowing for meaningful comparisons of $f_{\mathrm{esc}}^{\mathrm{abs}}$ at all redshifts (Izotov et al. 2016a; Steidel et al. 2018).

The method can be applied to both spectroscopic observations and to imaging where two or more filters are used to probe the LyC flux and the nonionizing UV. When spectra are used, the LyC flux has so far always been measured over a small wavelength range just shortward of 912 Å. Imaging may probe a wider spectral range in the LyC, depending on the filter widths used; for HST/WFC3 imaging of the strong LyC leaker Ion2, e.g., the LyC is probed over ∼730–890 Å restframe (Vanzella et al. 2016).

One question that presents itself in regard to this method is to what extent a measure of such a monochromatic escape fraction reflects the total escape fraction of all H-ionizing photons, including those at shorter wavelengths. For example, in their numerical simulations of escaping LyC radiation, Gnedin et al. (2008) find that the total escape fraction of ionizing photons below 912 Å is approximately a factor of 1.25 higher than the fraction escaping close to the Lyman edge. Furthermore, McCandliss & O’Meara (2017) show the relation between these two escape fractions for slabs of variable H I column density and fixed ionization fraction, finding fairly large enhancements (a factor of 2 and higher for “Lyman edge escape fractions” < 0.32, comparable to our sample). However, such a simplified geometry and the transmission curves computed in their work cannot be used to predict the resulting spectrum of density-bounded H II regions (cf., below, and Inoue 2010).

Generally, if both (i) the escape fraction of LyC photons is wavelength-independent and (ii) the only emission process at the measured spectral range is stellar emission, the ratio of the observed to intrinsic stellar SED will provide a correct measure of the total escape fraction of LyC photons. Therefore, both of these conditions need to be examined. In fact, the escape fraction is to first order wavelength-independent if the ISM can be described by a picket-fence model, where the “leaking” regions have a negligible H I column density and the total escape fraction is governed by the geometrical covering factor (see also Inoue 2010). The detailed analyses in this paper and earlier studies of LyC-emitting galaxies (Gazagnes et al. 2018; Steidel et al. 2018) have shown clear evidence in favor of the picket-fence model, although a low residual column density in the “escape channels” cannot be excluded. The second condition (ii) might not always be validated, because free-bound emission of hydrogen in ionized nebulae emits shortward of the Lyman edge, and may therefore contribute a non-negligible additional emission source close to the Lyman edge, as discussed for example by Inoue (2010). Theoretically, the strength of this emission primarily depends on the mean energy of the ionizing photons (i.e., on the detailed SED) and on the electron temperature of the ionized gas, and might lead to a “Lyman bump” in extreme cases (Inoue 2010)¹¹ To the best of our knowledge, this emission has so far not been identified, although it is expected. If present, our monochromatic determinations of $f_{\mathrm{esc}}^{\mathrm{abs}}$ could overestimate the true total escape fraction.

A second, independent method to determine the LyC escape fraction has also been used for the LzLCS survey and the literature samples; see Flury et al. (2022a,b). This method uses the de-reddened Hβ flux, which is proportional to the total number of ionizing photons absorbed (i.e., nonescaping) in the galaxy, and compares it to the observed LyC flux, using theoretical SEDs to relate the total ionizing photon flux to the monochromatic flux shortward of the Lyman edge. As discussed by Izotov et al. (2018b), this method therefore measures the total ionizing photon escape fraction. However, if strong free-bound emission is present, this method might also overestimate $f_{\mathrm{esc}}^{\mathrm{abs}}$. Comparing this method to the “monochromatic” method, Izotov et al. (2018b) find good agreement overall and no systematic differences (cf., Izotov et al. 2018a,b, 2021). Similarly, for the larger sample analyzed here, we do not find systematic differences between the two methods, although the scatter is fairly large (Flury et al. 2022a), and possibly driven by uncertainties related to the reddening corrections (both in the optical and UV). Indeed, as shown above (see Fig. 4), the $f_{\mathrm{esc}}^{\mathrm{abs}}$ values derived here solely from the UV spectrum depend quite strongly on the assumed attenuation law, making this probably the largest source of uncertainty of our analysis.

Finally, as it is known that theoretical SEDs of star-forming galaxies have a significant lack of ionizing photons at energies above 54 eV, which are required to explain nebular lines of He II frequently observed in low-metallicity SF galaxies (see e.g., Stasińska et al. 2015; Schaerer et al. 2019), one might wonder whether or not these photons could alter our analysis and play a role in the ionizing photon budget of hydrogen. The answer appears to be no, primarily because the number of photons emitted above the He II ionization potential (λ < 228 Å) is small compared to those ionizing hydrogen. Indeed, their relative ratio, Q₂/Q_O, can be determined empirically from the observed He II and H recombination lines and is typically ∼(0.2 − 2)×10⁻² in low-z galaxies with nebular He IIλ4686 detections (see Schaerer et al. 2019). Their contribution to the ionization of hydrogen is therefore negligible. From model SEDs, it is well known that the bulk of the H-ionizing photons are emitted relatively close to the Lyman edge; for example ∼60%−70% of the photons in the LyC are emitted longward of the He I ionization potential (24.6 eV, or 504 Å) at young ages (see e.g., Nakajima et al. 2018; Xiao et al. 2018). Although all photons with energies > 13.6 eV can ionize hydrogen, the majority are emitted at relatively low energies. In any case, methods using Hβ (or other H recombination lines) to determine the LyC escape fraction inherently probe the total H-ionizing photon budget.

8. Summary and conclusions

The Low-Redshift Lyman Continuum Survey (LzLCS, Flury et al. 2022a) recently observed 66 low-z (z ∼ 0.3) star-forming galaxies with HST/COS to measure the LyC in a large sample of galaxies for the first time. We augmented the LzLCS data, adding 23 previous LyC observations taken with COS (Izotov et al. 2016a,b, 2018a,b, 2021, and Wang et al. 2019), and the available LyC, UV, and optical spectra analyzed in a homogenous manner. The sample provides unique insight into the LyC emission, the LyC escape fraction, and many related properties; it covers galaxies with a wide range of physical properties (UV magnitude, metallicity, and others) allowing us to empirically study how the escape of LyC emission varies from galaxy-to-galaxy and, in particular, how this phenomenon is related to the ISM properties of these sources.

To study the ISM properties of these low-redshift galaxies, we analyzed the COS UV spectra in the framework of a picket-fence model with a uniform dust screen. First, we fit the UV stellar continuum (see Sect. 3.1) to constrain the massive star content of these galaxies. From our nonparametric fits, we derived direct properties such us the stellar ages and metallicities, UV dust attenuations (E_B−V), intrinsic UV magnitudes ($M_{1500}^{\mathrm{int}}$), and absolute LyC escape fractions ($f_{\mathrm{esc}}^{\mathrm{abs}}$). The stellar population properties will be discussed in a separate paper (Chisholm et al., in prep.).

The derived LyC escape fractions range up to 60% (85%) for the LzLCS (published) sample. Approximately 45% of the sample is not detected in the LyC, with 1σ upper limits of $f_{\text{esc}}^{\text{abs}}$ < 0.01. The absolute escape fraction strongly depends on the dust-attenuation parameter (E_B−V), with the strongest LyC emitters only found with low UV attenuation. We also study the impact of the assumed attenuation law on $f_{\mathrm{esc}}^{\mathrm{abs}}$: adopting the SMC extinction law instead of the attenuation law from Reddy et al. (2016a) leads to an increase in $f_{\mathrm{esc}}^{\mathrm{abs}}$ by ∼10%−50%, and the largest differences appear for galaxies with the lowest LyC escape fractions (Sect. 4).

With the stellar continuum fits in hand, we then measured the equivalent widths (W_λ) and residual fluxes (R_λ) of different interstellar absorption lines of H I (Lyman series) and metallic low ionization state (LIS) lines, such as Si IIλλ1190,1193, Si IIλ1260, O Iλ1302, and C IIλ1334 (see Sect. 3.2). With the picket-fence model, the residual fluxes allow us to infer a geometrical covering fraction of the UV source by the gas. Our analysis of the absorption lines and other galaxy properties has led us to the following main conclusions:

• Lyman continuum emitters (leakers) have weak but likely saturated H I and LIS lines (Sect. 5.1), and changes in W_λ are mainly driven by the neutral gas covering fraction (C_f(H I)). Metallicity plays a secondary role (Sect. 5.2).

• We find that leakers have high H I and LIS residual fluxes (R_λ) and low UV attenuations (A_UV). This indicates that leakers tend to host a dust-poor ISM, with low H I column density channels intermixed between higher column density regions (Sect. 5.2). Our large sample confirms the finding that these channels enable the escape of ionizing photons, as shown in other recent studies (Reddy et al. 2016b; Gazagnes et al. 2018; Steidel et al. 2018).

• Leakers have large Lyα equivalent widths (W_Lyα), and this correlates strongly with the H I and LIS equivalent widths, showing that the absorption lines trace the same gas that regulates the radiative transfer of Lyα photons. W_Lyα is therefore partly driven by the neutral gas covering fraction C_f(H I) (Sect. 6.1), as also found in earlier studies (e.g., Gazagnes et al. 2020; Pahl et al. 2020.

• We find H I and LIS covering fractions varying from ∼0.5 to 1 and ∼0.2 to 1, respectively. As in Reddy et al. (2016b) and Gazagnes et al. (2018), the LIS covering fraction is on average smaller than that of H I, and the two quantities are strongly correlated. This probably indicates that H I and metals are not evenly distributed within the host galaxy, but that they trace each other (Sect. 6.2).

• The LyC photon escape fraction ($f_{\mathrm{esc}}^{\mathrm{abs}}$) strongly correlates with both H I and LIS equivalent widths and residual fluxes of the lines (Sect. 6.3), and also with the UV attenuation (Sect. 6.4). We derive empirical calibrations to estimate $f_{\mathrm{esc}}^{\mathrm{abs}}$ from UV absorption line measurements.

• Finally, we show that simultaneous UV absorption line and attenuation measurements can, in general, predict the escape fractions of galaxies (Sect. 7.1). We apply our methods to literature measurements of UV LIS lines in z ∼ 4 − 6 galaxies, estimating modestly low LyC escape fractions of $f_{\text{esc}}^{\text{abs}}$ ≲ 0.1 for relatively bright, $M_{\text{1500}}^{\text{obs}}$ ≲ −21, high-redshift galaxies (Sect. 7.3).

In short, our detailed quantitative analysis of the UV spectra of 89 low-z galaxies shows that LyC leakage is a complex phenomenon (see Sect. 7.2), which mainly depends on the gas distribution (covering fractions) and dust content (UV attenuation) of the host galaxy. We demonstrate that the H I Lyman series lines shortward of Lyα and metallic absorption lines can be used to statistically predict the escape fraction of ionizing photons, as proposed by Chisholm et al. (2018). Although several uncertainties and open questions remain, our study demonstrates how LIS UV absorption lines can be used to study and constrain the escape of ionizing photons. These diagnostics will be invaluable and key in upcoming studies of galaxies up to the Epoch of Reionization, with the advent of the James Webb Space Telescope (JWST) and other 30-meter class telescopes.

Acknowledgments

A.S.L. and D.S. acknowledge support from Swiss National Science Foundation. S.R.F and A.E.J. acknowledge support from NASA through grant HST-GO-15626 from the Space Telescope Science Institute. R.A. acknowledges support from ANID Fondecyt Regular 1202007. M.T. acknowledges support from the NWO grant 016.VIDI.189.162 (“ODIN”). This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are associated with program GO 15626 (P.I. Jaskot). Additional work was obtained from the data archive at the Space Telescope Science Institute from HST proposals 13744, 14635, 15341, and 15639.

References

Appendix A: Fits and absorption lines results

Summary tables with the main products from our UV–SED fits (+errors) and the absorption line measurements (+errors) are shown in this appendix for both the LzLCS and literature samples. Tables A.1 and A.3 include the E_B−V average stellar age and metallicities, inferred $f_{\mathrm{esc}}^{\mathrm{abs}}$, and intrinsic (dust-free) and observed UV absolute magnitudes (in AB system). Tables A.2 and A.4 include the average H I and LIS equivalent widths, residual fluxes, and derived $f_{\mathrm{esc}}^{\mathrm{abs},\mathrm{HI}}$ and $f_{\mathrm{esc}}^{\mathrm{abs},\mathrm{LIS}}$. Sources are identified and sorted by its coordinates, and the strong/weak/nonleaker categories are differentiated according to the LzLCS criteria.

Appendix B: Absorption lines – saturation

The ratio of two absorption lines of the same element and ionization state can indicate whether the species in question is saturated or not. Let f_i, λ_i be the oscillator strength and central rest-frame wavelength of the line in question, where the 2(1) subscript represents the stronger (weaker) transition. If both transitions are saturated, the equivalent width ratio W₂/W₁ can be approximated by Draine (2011)

$$\begin{array}{c}\hfill {\displaystyle \frac{W_2}{W_1}}\approx {[1+{\displaystyle \frac{ln(f_2{\lambda }_2/f_1{\lambda }_1)}{ln({\tau }_1/ln2)}}]}^{(1/2)},\end{array}$$(B.1)

where τ₁ is the optical depth of the stronger transition, being τ₁ ≫ 1 the optically thick limit. If, on the opposite, both lines behave under the optically thin regime (τ₁ ≪ 1), the equivalent width doublet-ratio is simply:

$$\begin{array}{c}\hfill {\displaystyle \frac{W_2}{W_1}}={\displaystyle \frac{f_2{\lambda }_2}{f_1{\lambda }_1}}.\end{array}$$(B.2)

In practice, we assume τ₁ > 1 for saturation, and a smooth transition between the two regimes. To test the saturation of H I in our galaxy spectra, we compare the equivalent width of H Iλ1026 to H Iλ950. The use of lines pairs which are close in wavelength ensures that both lines are in the same optical regime simultaneously.

The observations are shown in Figure B.1, together with the regions delimiting the optically thick (saturation) and optically thin (linear) regimes, respectively, following the equations above. As stressed in the text, most of the sources are compatible with saturation for H I within the 1σ errors, and only a low number of them fall in the optically thin limit. Several sources are found to show a higher equivalent width in H Iλ950 than H Iλ1026, which is not expected from simple curve-of-growth analysis. The exact reasons for this behavior are not known but could be due to multiple effects, such as systematic errors in the continuum placement in the blue part of the spectrum, blends with other lines such as O I, and other causes. In any case, in most of the cases, they are compatible with the saturated region at 1σ. Spectra with higher resolution and S/N, such as the subset of galaxies analyzed by Gazagnes et al. (2018), would be required to better examine the question of saturation.

Fig. B.1. Absorption line H I equivalent width comparisons for studying the condition of saturation. The optically thin (light-gray) and optically thick (dark-gray) regions are delimited using the approximations given in Draine (2011). As in the entire document, purple symbols represent the leakers while the black symbols show the nonleakers. Red squares correspond to the Izotov et al. (2016a,b, 2018a,b, 2021) and Wang et al. (2019) galaxies.

Appendix C: Absorption lines – systematic errors

Here, we carefully account for the impact of the resolution and S/N in the measured covering fraction. The spectrograph resolution tends to make the absorption lines wider but less deep. Together with the noise and sensitivity –and other possible instrumental effects– this can lead to a systematically overestimated measurement of the residual flux.

We use absorption line simulations to study how the resolution (R) and S/N both impact our measurements of the residual flux. We first simulate Lyδ, Lyβ, Si IIλ1260, and C IIλ1334 line intensities (I_λ) in a foreground dust-screen picket-fence model:

$$\begin{array}{c}\hfill I_{\lambda }=C_f(\lambda )\times exp(-{\sigma }_{\lambda }{N}_{\lambda })+1-C_f(\lambda ),\end{array}$$(C.1)

where C_f(λ) represents the covering fraction of the line in question, i.e., the fraction of sight lines for which the transition is optically thick around the host galaxy. C_f(λ) can be inferred from the line depth as C_f(λ) = 1 − R(λ) for the dust-screen picket-fence model. The τ_λ = σ_λN_λ product is usually known as the optical depth of the line, and the λ subscripts belong to the transition in question. The line cross-section, σ_λ, which shapes the absorption profile, is given by a usual Voigt function so that:

$$\begin{array}{c}\hfill {\sigma }_{\lambda }={\displaystyle \frac{\sqrt{\pi }{e}^2{f}_{\lambda }}{m_{e}c}}{\displaystyle \frac{\lambda }{b}}\mathrm{Voigt}(\lambda ,{\mathrm{A}}_{\lambda },\mathrm{b}).\end{array}$$(C.2)

For the Voigt(λ,A_λ,b) function, we use the numerical approximation described in Smith et al. (2015). The Doppler broadening parameter is input through b (km s⁻¹). Finally, oscillator strengths f_λ and Einstein coefficients A_λ are taken from the NIST database¹². We do not consider any velocity gradient effect (outflows) in the mock realizations.

For our purposes, we assume a Doppler parameter of b = 125 km/s, and typical column densities in SF galaxies of N_λ = 10¹⁹, 10¹⁹, 10¹⁷, 10¹⁷ cm⁻² for Lyδ, Lyβ, Si IIλ1260, and C IIλ1334 lines, respectively. We also explore a wide range of (N, b) parameters and conclude that results do not change drastically as soon as the lines are in the saturation region of the curve of growth. The covering fraction is fixed to C_f(λ) = 0.85 throughout the processes.

The mock lines are convolved with a Gaussian kernel to match the COS/G140L nominal resolution (R = 1000) and interpolated into the rest-frame COS/G140L wavelength spacing (0.5621/(1+z) Å. where z = 0.3). A set of Gaussian S/N are then randomly added to every pixel. Finally, the residual flux is measured in the same way as in the original real spectra, and related to the covering fraction by Eq. (4) (uniform dust-screen). This process is repeated 500 times for every S/N value, and the median covering fraction is compared to the input value.

This method is summarized in Fig. C.1, which shows the covering fraction relative difference (in percentage) as a function of the input S/N per pixel in the simulations. The four colored lines corresponds to the different lines considered. The upward arrows at the bottom correspond to the measured S/N for the LzLCS sample around the H Iλ950 line, for comparison (the S/N has been measured at the continuum level close to the line in question). In general, the observed continuum S/N increases at shorter wavelengths, i.e., it is larger for Lyδ and Lyβ, but systematically lower for Si IIλ1260 and C IIλ1334.

Fig. C.1. Estimation of the systematic error on the observed covering fraction as a function of the spectral S/N per pixel, for different simulated lines (see legend). A resolution of R = 1000 is assumed. The upward arrows at the bottom correspond to the measured S/N for the LzLCS sample around the H Iλ950 line. Horizontal dashed line limits the 10% representative value.

In all cases, and as indicated in the text, the systematic error on the covering fraction spans between 5% and 20% for these four representative lines (with a typical value of 10%). These values are within the original error bars we are reporting for single covering fraction measurements. We therefore leave our initial lines measurements and errors unaltered, and we prefer not to correct the C_f(λ) values by this systematic offset.

Appendix D: The LINMIX package

The LINMIX estimator (Kelly 2007) is a Python-based Bayesian fitting code that allows the user to model a two-dimensional data set (x, y) with a linear regression, accounting for errors on both variables and intrinsic random scatter, and with the capability of including censored (upper or lower limits) data.

LINMIX assumes that measured data points (x, y) can be sampled from a two-dimensional Gaussian distribution P₁ with mean (ξ, η), and covariance matrix determined from measurement uncertainties. It also assumes that the measurement errors in the x- and y-directions are normally distributed and uncorrelated.

The “true” value of the η variable is drawn from another Gaussian distribution P₂, with mean α + βη, and variance ${\sigma }_y^2$, where β describes the slope, α the intercept of the straight line, and ${\sigma }_y^2$ is the intrinsic y−dispersion (also normally distributed). Finally, the “true” value of the independent variable ξ is sampled from a weighted sum of K Gaussian distributions P₃ (we choose K = 2). The distributions P₁, P₂, P₃ are convolved hierarchically to compute the full likelihood of obtaining the observed data (x, y) given the parameters of the fit: α, β, ${\sigma }_y^2$.

Assuming uniform prior distributions for the three parameters: α, β ∈ ( − ∞, ∞), ${\sigma }_y^2$ ∈ [0, ∞); a Markov Chain Monte-Carlo (MCMC) algorithm is used to sample from the posterior distributions until the convergence of the chains (×4 chains, 5000 to 10 000 iterations each). For a practical example, see e.g., Kennicutt & De Los Reyes (2021). A PYTHON version of the LINMIX code can be found in https://linmix.readthedocs.io/en/latest/index.html.

Appendix E: Main results at a glance

This appendix is dedicated to illustrating the global behavior of the absolute LyC photon escape fraction ($f_{\mathrm{esc}}^{\mathrm{abs}}$) with the main UV spectral properties measured in this work, namely: R(H I), R(LIS), W_H I, W_LIS and E_B−V. These results are presented in Figs. 11–13, but here we focus on the median values of the working sample as a take- home message. To do this, we compute the median value of R(H I), R(LIS), $W_{\text{H}\text{I}}$, W_LIS and E_B−V in three regimes of $f_{\mathrm{esc}}^{\mathrm{abs}}$, specifically at $f_{\text{esc}}^{\text{abs}}$ ∈ [0%−5%), [5%−10%), [10%−100%). We additionally apply a bootstrap + Monte Carlo algorithm to the number of sources at each $f_{\mathrm{esc}}^{\mathrm{abs}}$ interval allowing us to compute the 0.16 and 0.84 quantiles at each bin. The results are summarized in Table E.1.

In Fig. E.1, the running medians are plotted as a function of $f_{\mathrm{esc}}^{\mathrm{abs}}$ (in black), together with the individual LzLCS+literature measurements (gray and red dots in the background). The reader may find a detailed discussion about these results in Sect. 6. A proper analysis of the median results of the sample will be investigated in future LzLCS publications, by the use of stacked spectra.