# Linking UV spectral properties of MUSE Ly α emitters at z 3 to Lyman continuum escape

I Kramarenko, J Kerutt, A Verhamme, P Oesch, L Barrufet, J Matthee, H

Kusakabe, I Goovaerts, T Thai

► To cite this version:

I Kramarenko, J Kerutt, A Verhamme, P Oesch, L Barrufet, et al.. Linking UV spectral properties of MUSE Ly α emitters at z 3 to Lyman continuum escape. Monthly Notices of the Royal Astronomical Society, 2024, 527 (4), pp.9853-9871. 10.1093/mnras/stad3853. hal-04505602

HAL Id: hal-04505602

https://hal.science/hal-04505602v1

Submitted on 15 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Linking UV spectral properties of MUSE Ly α emitters at $z \gtrsim 3$ to Lyman continuum escape

I. G. Kramarenko<sup>D, 1,2*</sup> J. Kerutt,<sup>1,3</sup> A. Verhamme,<sup>1,4</sup> P. A. Oesch<sup>D, 1,5</sup> L. Barrufet<sup>D, 1</sup> J. Matthee<sup>D, 2,6</sup>
H. Kusakabe,<sup>1</sup> I. Goovaerts<sup>7</sup> and T. T. Thai<sup>7,8</sup>

Accepted 2023 December 12. Received 2023 November 16; in original form 2023 May 12

## ABSTRACT

The physical conditions giving rise to high escape fractions of ionizing radiation (LyC $f_{esc}$) in star-forming galaxies – most likely protagonists of cosmic reionization – are not yet fully understood. Using the VLT/MUSE observations of ~1400 Ly α emitters at $2.9 < z < 6.7$, we compare stacked rest-frame UV spectra of candidates for LyC leakers and non-leakers selected based on their Ly α profiles. We find that the stacks of potential LyC leakers, i.e. galaxies with narrow, symmetric Ly α profiles with small peak separation, generally show (i) strong nebular O III]λ1666, [Si III]λ1883, and [C III]λ1907 + C III]λ1909 emission, indicating a high-ionization state of the interstellar medium (ISM); (ii) high equivalent widths of He IIλ1640 (~1 – 3 Å), suggesting the presence of hard ionizing radiation fields; (iii) Si II*λ1533 emission, revealing substantial amounts of neutral hydrogen off the line of sight; (iv) high C IVλλ1548, 1550 to [C III]λ1907 + C III]λ1909 ratios (C IV/C III] $\gtrsim$0.75), signalling the presence of low column density channels in the ISM. In contrast, the stacks with broad, asymmetric Ly α profiles with large peak separation show weak nebular emission lines, low He IIλ1640 equivalent widths ($\leq$1 Å), and low C IV/C III] ($\leq$0.25), implying low-ionization states and high-neutral hydrogen column densities. Our results suggest that C IV/C III] might be sensitive to the physical conditions that govern LyC photon escape, providing a promising tool for identification of ionizing sources among star-forming galaxies in the epoch of reionization.

**Key words:** galaxies: high-redshift – ISM: lines and bands – dark ages, reionization, first stars – ultraviolet: galaxies – ultraviolet: ISM.

## 1 INTRODUCTION

Ultraviolet (UV) and X-ray radiation emitted from the first stars and galaxies likely ionized the neutral hydrogen in the intergalactic medium (IGM) between $z \approx 12$ (e.g. Hinshaw et al. 2013; Planck Collaboration 2016) and $z \approx 6$ (e.g. Fan et al. 2006; McGreer, Mesinger & D'Odorico 2015). This last phase transition of the Universe is known as the Epoch of Reionization (EoR). The hydrogen reionization could have been potentially powered by active galactic nuclei (AGNs) due to their brightness and high escape fractions of ionizing radiation ($f_{esc}$). However, the AGN number density is likely to be insufficient to maintain this process at $z \gtrsim 6$ (e.g Parsa, Dunlop & McLure 2018; Kulkarni, Worseck & Hennawi 2019; Shen et al. 2020). It is currently believed that the most promising candidates

for the sources of the ionizing photons are young massive stars in star-forming galaxies (SFGs), although their contribution is still unclear, primarily due to the uncertainties associated with $f_{esc}$.

Recent studies have shown that $f_{esc}$ in SFGs at $z \gtrsim 6$ has to be on average at least 10–20 percent to successfully reproduce the observed constraints on the reionization history (e.g. Ouchi et al. 2009; Bouwens et al. 2015; Robertson et al. 2015; Khaire et al. 2016; Naidu et al. 2020). Direct measurements of $f_{esc}$, however, are severely hampered at such redshifts because ionizing radiation, or Lyman continuum (LyC, $\lambda < 912$ Å), is absorbed by the neutral hydrogen in the IGM (e.g. Madau 1995; Inoue et al. 2014). The current solution to this problem lies in the studies of LyC-leaking SFGs in the local Universe (e.g. Izotov et al. 2016a,b; Schaerer et al. 2018; Izotov et al. 2018a,b; Wang et al. 2019; Flury et al. 2022; Xu et al. 2022) and at intermediate redshifts ($1 \lesssim z \lesssim 4$; e.g. Vanzella et al. 2015; Shapley et al. 2016; Rivera-Thorsen et al. 2017; Vanzella et al. 2018; Marques-Chaves et al. 2021; Saxena et al. 2022a), often

<footnote n="1"><sup>1</sup>Department of Astronomy, Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland</footnote>

<footnote n="2"><sup>2</sup>Institute of Science and Technology Austria (ISTA), Am Campus 1, A-3400 Klosterneuburg, Austria</footnote>

<footnote n="3"><sup>3</sup>Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, NL-9700 AV Groningen, the Netherlands</footnote>

<footnote n="4"><sup>4</sup>University of Lyon, University of Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="5"><sup>5</sup>Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen, Jagtvej 128, København N, DK-2200, Denmark</footnote>

<footnote n="6"><sup>6</sup>Department of Physics, ETH Zürich, Wolfgang-Pauli-Strasse 27, Zürich CH-8093, Switzerland</footnote>

<footnote n="7"><sup>7</sup>Aix Marseille Université, CNRS, CNES, LAM (Laboratoire d'Astrophysique de Marseille), UMR 7326, IPhU, F-13388 Marseille, France</footnote>

<footnote n="8"><sup>8</sup>Department of Astrophysics, Vietnam National Space Center, Vietnam Academy of Science and Technology, 18, Hoang Quoc Viet, Nghia Do, Cau Giay, Ha Noi, Vietnam</footnote>

<footnote>*E-mail: im.kramarenko@gmail.com</footnote>
accompanied by spectral stacking with the aim of obtaining higher
signal-to-noise (S/N) data and averaging out the spatial variations of
the IGM transmission (e.g. Marchi et al. 2017; Steidel et al. 2018;
Meštrić et al. 2021). The low-z observations help to develop indirect
tracers of $f_{esc}$ which can then be extrapolated to the galaxies in the
EoR.

Historically, one of the first proposed indicators of LyC leakage
was the [O III]λλ4959,5007 to [O II]λλ3726,3729 ratio (O³²; Jaskot &
Oey 2013; Nakajima & Ouchi 2014). Recent observations suggest,
however, that the correlation between O₃ and $f_{esc}$ is weak (e.g. Naidu
et al. 2018; Izotov et al. 2018b), possibly due to the complex geometry
and the kinematics of the interstellar medium (ISM) modulating LyC
escape (Bassett et al. 2019; McKinney et al. 2019; Nakajima et al.
2020). A better proxy for $f_{esc}$ might be an ensemble of UV low-
ionization state (LIS) absorption lines (e.g. Si IIλ1260, O Iλ1302,
and C IIλ1334). Observationally, the neutral gas covering fraction
($f_{cov}$) inferred from the LIS lines correlates with the fraction of
LyC photons escaping through low column density channels in the
ISM (e.g. Reddy et al. 2016; Chisholm et al. 2018; Gazagnes et al.
2018; Saldana-Lopez et al. 2022), although simulations demonstrate
a substantial scatter in this relation (Mauerhofer et al. 2021). Alterna-
tively, the ISM conditions driving LyC escape could be probed with
resonant nebular emission lines sensitive to the radiation transfer
effects. These include the Mg IIλλ2796,2803 doublet (Feltre et al.
2018; Henry et al. 2018; Chisholm et al. 2020), the C IVλλ1548,1550
doublet (Schaerer et al. 2022; Saxena et al. 2022b), and the Lyman-
alpha (Ly α) line.

The Ly α–LyC relationship has been extensively explored in the literature. Verhamme et al. (2015) studied Ly α and LyC escape in the two following configurations of the ISM: (i) density-bound nebulae, and (ii) ionization-bounded nebulae with holes, also referred to as the ‘picket-fence’ model. They found that in the first case the Ly α line has a narrow profile and a small velocity offset with respect to systemic redshift ($v_{peak}^{red} \lesssim 150$ km s⁻¹), whereas in the second case the Ly α line is at systemic redshift. They proposed to use small Ly α peak separation ($v_{sep} \lesssim 300$ km s⁻¹) for double-peaked Ly α profiles as an indicator of low neutral hydrogen column density and thus potentially high $f_{esc}$. Dijkstra, Gronke & Venkatesan (2016) carried out radiative transfer simulations of the clumpy ISM and found that the LyC-emitting galaxies have narrower, more symmetric Ly α line profiles. Kimm et al. (2019) and Kakiichi & Gronke (2021) obtained similar results from radiation hydrodynamic simulations of turbulent molecular clouds. Finally, recent observational studies have shown that local SFGs with high $f_{esc}$ typically have a narrow Ly α line with small peak separation, confirming the theoretical expectations (Izotov et al. 2016a, 2018b, 2021, 2022; Verhamme et al. 2017; Flury et al. 2022).

These results demonstrate that the Ly α line profile could provide
robust predictions for $f_{esc}$, potentially even at high redshift. Ly α
profiles have been already used to select candidates for LyC leakers
and non-leakers from a representative sample of Ly α emitters
(LAEs) at $z \approx 2$ (Naidu et al. 2022). In this paper, we apply the same
technique to examine the possible relationship between LyC escape
and rest-frame UV spectral properties of ~1400 LAEs observed with
the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010)
at 2.9 < z < 6.7. We combine the 3D spectroscopic data from the
MUSE-Wide (Urrutia et al. 2019) and the MUSE Hubble Ultra Deep
Field (MUSE HUDF; Bacon et al. 2017) surveys and select a sample
of LAEs with reliable detections of the Ly α line (Section 2). In
Section 3, we use the observed properties of the Ly α line profile (e.g.
peak separation) to select groups of LAEs with potentially different
$f_{esc}$. In Section 4, we stack the individual spectra of our galaxies to

obtain high-S/N detections of the rest-frame UV lines, namely, high-
ionization nebular emission lines (e.g. C IVλλ1548,1550, He IIλ1640,
and [C III]λ1907 + C III]λ1909), ISM absorption lines (e.g. O Iλ1302
+ Si IIλ1304 and C IIλ1334), and fine-structure emission lines (e.g.
Si II*λ1533). In Section 5, we compare the stacked spectra of LyC-
leaker and non-leaker candidates to examine the ionizing properties
of galaxies and the physical conditions in their ISM indicative of
high $f_{esc}$. In Section 6, we present the summary of our main findings
and discuss the implications of our work for reionization studies.

## 2 DATA AND SAMPLE

### 2.1 Spectroscopic data from MUSE

For the purpose of this study, we use the spectroscopic data from the MUSE-Wide (Urrutia et al. 2019) and MUSE HUDF (Bacon et al. 2017; Data Release I) surveys taken during the guaranteed time observations (GTOs) of the MUSE consortium. The MUSE integral-field spectrograph is a powerful tool for spectroscopic diagnostics of emission-line galaxies thanks to a combination of a large field of view (1 × 1 arcmin²), high resolving power (ranging from R ≈ 1800 in the blue to ≈4000 in the red) and a large simultaneous spectral range (4750–9350 Å; Bacon et al. 2015). The GTO programmes take full advantage of these unique spectroscopic capabilities of the MUSE instrument, enabling studies of large, untargeted samples of faint (M<sub>UV</sub> down to ≈−16) LAEs at z > 3.

Being different in depths and sizes, MUSE-Wide and MUSE HUDF effectively complement each other, probing the Ly α-bright (L<sub>Ly α</sub> ≳ L*) and faint (L<sub>Ly α</sub> ≲ 0.1 L*) populations of LAEs, respectively. MUSE-Wide, a relatively wide and shallow component of the GTO ‘wedding cake’, covers 100 1 × 1 arcmin² fields at 1 hr observation time, among which 60 are in the CANDELS (Grogin et al. 2011; Koekemoer et al. 2011)/GOODS-South (GOODS-S; Giavalisco et al. 2004) field, and 23 are in the CANDELS/COSMOS (Scoville et al. 2007) field. The other eight pointings of MUSE-Wide are in the HUDF parallel fields. MUSE HUDF, the pencil-beam component of the GTO ‘wedding cake’, targets a much smaller area but for longer exposure times. This survey consists of nine fields of 10 hr depth in HUDF, completing the 100 fields of MUSE-Wide, plus a single UDF-10 field with 31 hr exposure time.

In this paper, we use the 1D spectra extracted from the reduced MUSE-Wide and MUSE HUDF data cubes by Schmidt et al. (2021). The untargeted search for emission line sources in the MUSE data cubes was performed using the LSDCAT tool (Herenz & Wisotzki 2017). Each of the detected sources was then assigned a subjective confidence C between zero (lowest confidence) and three (highest confidence) depending on the uncertainty of the line classification. The latter was carried out using the custom graphical user interface QTCLASSIFY (Kerutt 2017). Finally, the optimally extracted spectra were obtained with TDOSE (Schmidt et al. 2019) using morphological Hubble Space Telescope (HST) models as templates. To preserve the self-consistency of this procedure, we decided against incorporating the data from the MUSE eXtremely Deep Field (MXDF) survey added to the second data release of MUSE HUDF (Bacon et al. 2023).

### 2.2 Sample selection

We select MUSE LAEs for our analysis as follows. First, we only
consider LAEs with a confidence $C \ge 2$, i.e. sources in which
Ly α is detected at S/N > 5 and the line profile is compatible
with typical Ly α shapes (see Bacon et al. 2023 for a detailed
description of MUSE confidence levels). Next, we note that some spectra might represent the same galaxy in the total LAE sample due to a partial overlap between the MUSE-Wide fields, the UDF mosaic and the UDF-10 field. Among these duplicates, we consider spectra with the highest observation time only. Finally, we exclude 19 superpositions (spatially overlapping sources at different redshifts) and three galaxies classified as AGN in the 7 Ms Chandra Deep Field-South Survey catalogues (Luo et al. 2017).

Our final sample amounts to 1422 LAEs in the redshift range $2.92 < z < 6.64$, with 697 objects from the MUSE-Wide survey and 725 objects from the MUSE HUDF survey. Unlike Feltre et al. (2020), who carried out a similar spectroscopic analysis of MUSE LAEs at $z > 3$, we include in our sample objects at redshifts where the MUSE spectral range does not cover several important rest-frame UV lines (e.g. the [C III]λ1907 + C III]λ1909 doublet at $z > 3.9$). We find that the evolution of redshift coverage with wavelength does not affect our conclusions (see Section 4).

## 2.3 Sample properties

The distributions of the observed properties of our LAEs are shown in Fig. 1, with data taken from Kerutt et al. (2022). The redshifts (Fig. 1a) are measured using the Ly α line and have a median value of 3.89. We caution that redshift estimates based on Ly α emission have systematic uncertainties due to resonant scattering of Ly α photons by neutral gas in the ISM. High-neutral hydrogen column density leads to a significant shift of the peak of the Ly α line ($v_{\text{peak}}^{\text{red}}$), typically of the order of a few hundred km s$^{-1}$ (e.g. Shapley et al. 2003; Steidel et al. 2010; Song et al. 2014; Hashimoto et al. 2015; Muzahid et al. 2020; Matthee et al. 2021). However, having accurate redshift measurements is crucial for our study to correctly perform spectral stacking. We describe the approach that we used to recover systemic redshifts in Section 4.

As shown in Fig. 1(b), the absolute UV magnitudes of our LAEs span the range of $-21.9 < M_{\text{UV}} < -15.6$, with a median value of $-18.1$. Thanks to the inclusion of the MUSE HUDF data, our sample probes $\sim 1 - 2$ dex fainter magnitudes compared to the samples of LAEs selected using the narrow-band technique (e.g. Ouchi et al. 2008; Kashikawa et al. 2011; Zheng et al. 2014; Nakajima et al. 2018a; Matthee et al. 2021), or samples of Lyman-break galaxies (e.g. Shapley et al. 2003; Stark et al. 2010). Magnitudes as faint as $M_{\text{UV}} \sim -16$ or even fainter, down to $M_{\text{UV}} \sim -14$, have been previously achieved mostly in the studies of gravitationally lensed LAEs (e.g. Stark et al. 2014; de La Vieuville et al. 2020; Bouwens et al. 2022).

The distribution of the UV-continuum slopes ($\beta$) is shown in Fig. 1(c). These estimates are based on the HST observations of LAEs with at least two detections in the HST filter bands (856 objects, or 60 per cent of the total sample). The fact that most of our LAEs have blue UV slopes with a median $\beta = -2.1$ suggests low dust content, as has been previously reported in other LAE studies (e.g. Matthee et al. 2021). Therefore, we can safely neglect the effect of dust attenuation when analysing the spectra of our LAEs (Section 5).

## 2.4 Ly α line profile statistics

In this study, we select candidates for LyC-leakers and non-leakers using the properties of the Ly α line profile such as the peak separation, full width at half-maximum (FWHM) and asymmetry parameter. We show the distributions of the Ly α profile measurements in Fig. 1(d)–(f). We find that our galaxies exhibit a large variety of Ly α profiles in terms of the number of peaks, width, and

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Observed properties of LAEs studied in this work (data from Kerutt et al. 2022). The distributions are shown for the total sample (1422 LAEs; dark blue) and the sample of galaxies with double-peaked Ly α profiles (369 LAEs; cyan). (a)–(f): redshift, absolute UV magnitude, UV continuum slope, Ly α peak separation, Ly α FWHM, and Ly α asymmetry parameter. Dashed vertical lines indicate the median values of the distributions.</figcaption></figure>

skewness (see also fig. 5 from Kerutt et al. 2022). This implies a wide range of the physical conditions in the ISM that shape radiative transfer processes.

The peak separations – distances between the red and the blue peaks of the double-peaked Ly α profiles – range from $\approx 100$ to $\approx 1200$ km s$^{-1}$, with a median value of $\approx 470$ km s$^{-1}$ (Fig. 1d). We note that only 369 LAEs, or 26 per cent of the total sample, have resolved double-peaked Ly α profiles. By comparing the UV properties of these galaxies with that of the total sample, we find that the double-peaked Ly α sample is representative of the parent population (dark blue and cyan histograms in Fig. 1).

There might be several explanations for the lack of the double-peaked structure in most of the observed Ly α profiles. In some cases, low-S/N data and/or the presence of the IGM absorption hinder the detection of the blue peak (hereafter the blue bump), which is typically several times weaker than its red counterpart. Alternatively, the blue bump might be blended with the red peak due to the limited spectral resolution of the MUSE instrument, effectively resulting in a single-peaked Ly α profile. We discuss the limitations of using
the double-peaked Ly α sample for the analysis of LyC escape in Section 3.

In addition to the Ly α peak separation, we use the FWHM and asymmetry parameter of the red peak of the Ly α line in our selection of LyC-leaker candidates. These properties are measured in Kerutt et al. (2022) by fitting the Ly α profile with an asymmetric Gaussian function:

$$f(\lambda) = A \exp\left(-\frac{(\lambda - \lambda_0)^2}{2\sigma_{\text{asym}}^2}\right) + f_0, \quad (1)$$

where A is the amplitude, $\lambda_0$ is the wavelength of the red peak, $\sigma_{\text{asym}}$ is the asymmetric dispersion, and $f_0$ is the continuum level. The asymmetric dispersion is described by $\sigma_{\text{asym}} = a_{\text{asym}}(\lambda - \lambda_0) + d$, where $a_{\text{asym}}$ is the asymmetry parameter and $d$ is the typical width of the line. The asymmetry parameter is positive for the vast majority of the sample (87 per cent) and has a median value of +0.16, indicating that the Ly α line typically has a red wing (Fig. 1f). Negative values of $a_{\text{asym}}$ suggest the presence of a blue wing which might be the case for LAEs with an unresolved blue bump blended with the red peak. The FWHMs range from ≈60 to ≈440 km s⁻¹, with a median value of ≈220 km s⁻¹ (Fig. 1e). The FWHMs are corrected for the wavelength-dependent line spread function (LSF) of MUSE (Bacon et al. 2017) under the assumption that the Ly α line profile can be approximated by a Gaussian (see Kerutt et al. 2022 for the details).

## 3 SELECTING LYC-LEAKER CANDIDATES

In this section, we classify our LAEs as potential LyC leakers and non-leakers using the properties of the Ly α line profile described in Section 2 (peak separation, FWHM and asymmetry parameter). We briefly discuss the physics behind the relationship between the Ly α line profile and LyC escape to motivate our selection criteria. In addition, we estimate the median $f_{\text{esc}}$ for each resulting group of LAEs wherever possible. At the end of this section, we examine the global UV properties ($M_{\text{UV}}$, $\beta$) of LyC-leaker and non-leaker candidates, and make a comparison with the literature to place the results of our classification in a broader context.

The accuracy of our classification procedure depends on uncertainties of the observed Ly α properties, i.e. an LAE might be attributed to a wrong group if its Ly α peak separation, FWHM and/or asymmetry parameter have high statistical uncertainties. Thus, on the one hand, we follow a strategy of excluding LAEs whose observed Ly α properties are highly uncertain. On the other hand, we keep the subsamples large enough to achieve high-S/N values of the stacked spectra to be able to detect weaker rest-frame UV lines compared to Ly α (Sections 4 and 5).

### 3.1 Ly α peak separation

The Ly α peak separation is considered as one of the most reliable tracers of $f_{\text{esc}}$ in low-redshift galaxies (e.g. Izotov et al. 2021). In the density-bounded scenario, small peak separation indicates that Ly α photons escape from an almost fully ionized ISM with little scattering thanks to low column densities of neutral hydrogen (e.g. Verhamme et al. 2015). The same physical conditions allow LyC photons escape from the ISM without being absorbed by neutral gas.

The limited number of LAEs with resolved blue bumps (see Section 2.4) prevents us from having more than two groups selected based on the peak separation. After excluding 40 objects with $\sigma_{v_{\text{sep}}} > 200$ km s⁻¹, we obtain a sample of 329 LAEs with resolved double-peaked Ly α line profiles (23 per cent of the total sample).

From this sample, we select two groups of roughly the same size with LAEs having $v_{\text{sep}} < 463$ km s⁻¹ (165 objects; potential leakers) and $v_{\text{sep}} > 463$ km s⁻¹ (164 objects; potential non-leakers), respectively.

Izotov et al. (2018b) collected $v_{\text{sep}}$ measurements for a sample of low-redshift LyC-leaking galaxies and found the following empirical relationship between $f_{\text{esc}}$ and $v_{\text{sep}}$:

$$f_{\text{esc}}(v_{\text{sep}}) = \frac{3.23 \cdot 10^4}{v_{\text{sep}}^2} - \frac{1.05 \cdot 10^2}{v_{\text{sep}}} + 0.095, \quad (2)$$

where $v_{\text{sep}}$ is in km s⁻¹. We use this equation to estimate the expected median LyC $f_{\text{esc}}$ and obtain values of 6.0 per cent and 1.1 per cent for the low- and high-$v_{\text{sep}}$ subsamples, respectively. Therefore, our classification of a galaxy as a LyC leaker based on the Ly α peak separation allows lower values of $f_{\text{esc}}$ compared to the thresholds usually adopted in other studies ($f_{\text{esc}} > 20$ per cent; Naidu et al. 2022; or $f_{\text{esc}} > 10$ per cent; Schaerer et al. 2022).

We note that the Ly α line profile can be affected by the IGM (and the CGM) absorption, especially at high redshift where the Universe becomes more neutral. Taking into account the stochastic nature of this effect, we caution that the resulting $v_{\text{sep}}$ subsamples might be ‘contaminated’ by LAEs observed along the sightlines with suppressed Ly α transmission. Importantly, both theoretical and observational works show that the neutral hydrogen in the IGM predominantly attenuates the blue part of the Ly α line (e.g. Laursen, Sommer-Larsen & Razoumov 2011; Hayes et al. 2021). This highlights the importance of other properties of the Ly α profile (FWHM, asymmetry parameter) which are measured from the red peak of the Ly α line and therefore less sensitive to the IGM absorption.

### 3.2 Ly α FWHM

Similarly to the Ly α peak separation, the Ly α FWHM traces the neutral hydrogen column density in the ISM. For the density-bounded geometry, a narrower Ly α line profile suggests less scattering of Ly α photons, higher fractions of ionized gas in the ISM and higher $f_{\text{esc}}$. We select 964 LAEs with $\sigma_{\text{FWHM}_{\text{Ly α}}} < 100$ km s⁻¹ and split them into four $\text{FWHM}_{\text{Ly α}}$ quartiles (Q1–Q4) of 241 LAEs each as follows: (Q1) $\text{FWHM}_{\text{Ly α}} < 161$, (Q2) $\text{FWHM}_{\text{Ly α}} \in [161, 208]$, (Q3) $\text{FWHM}_{\text{Ly α}} \in [208, 260]$, and (Q4) $\text{FWHM}_{\text{Ly α}} > 260$ km s⁻¹. We find that 58 objects (24 per cent) in the group of LyC-leaker candidates (Q1, the lowest FWHMs) have Ly α line profiles with two resolved peaks. Most of these galaxies (39 LAEs, or 67 per cent) fall into the low-$v_{\text{sep}}$ subsample (Fig. 2, top), reflecting the fact that $\text{FWHM}_{\text{Ly α}}$ correlates with $v_{\text{sep}}$ since both depend on the neutral hydrogen column density as discussed above.

Expanding shell models predict that $\text{FWHM}_{\text{Ly α}}$ could be approximated by half of the peak separation (Verhamme et al. 2018). This allows us to use equation (2) once again, this time to estimate $f_{\text{esc}}$ for the $\text{FWHM}_{\text{Ly α}}$ subsamples. We caution that the relation between $\text{FWHM}_{\text{Ly α}}$ and $v_{\text{sep}}$ holds as long as the medium enables at least some degree of Ly α photon scattering, otherwise the Ly α profile would have a single peak at systemic velocity. However, this should be the case for most SFGs since the minimal column densities required for a substantial number of scattering events are extremely low (for instance, the Ly α forest is observed at column densities as low as $n_{\text{HI}} \sim 10^{13}$ cm⁻²).

We find that the expected median $f_{\text{esc}}$ ranges from 1.1 per cent in the highest $\text{FWHM}_{\text{Ly α}}$ quartile (Q4) to 15 per cent in the lowest $\text{FWHM}_{\text{Ly α}}$ quartile (Q1). Therefore, the strongest leakers selected by the Ly α FWHM are likely to have higher escape fractions than
### (i) double-peaked Ly$\alpha$ sample

<figure><img src="image_3.png" /></figure>

### (ii) total sample

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Venn diagram showing the samples of LyC-leaker candidates selected as having low Ly$\alpha$ peak separation ($<463$ km s$^{-1}$; cyan), FWHM ($<161$ km s$^{-1}$; dark blue), and asymmetry parameter ($<0.19$; red). *Top*: LAEs with double-peaked Ly$\alpha$ profiles. *Bottom*: Total sample. The estimates of the median $f_{\text{esc}}$ are based on the fit from Izotov et al. (2018b). All three properties of the Ly$\alpha$ line act as tracers of $f_{\text{esc}}$ due to their high sensitivity to the neutral hydrogen column density – a key parameter of the ISM controlling escape of ionizing photons.</figcaption></figure>

the LAEs with the smallest peak separations ($f_{\text{esc}} \approx 6$ per cent). None the less, even the lowest FWHM$_{Ly\alpha}$ quartile (Q1) has a median $f_{\text{esc}}$ smaller than the ones observed in the individual LyC leakers at $2 \lesssim z \lesssim 4$ ($\gtrsim 20$ per cent; see Izotov et al. 2021 for a review). Most likely, we are probing lower escape fractions because our sample is free of a selection bias typical for the individual detections which favour galaxies with the highest LyC fluxes (see Steidel et al. 2018 for a discussion).

## 3.3 Ly$\alpha$ asymmetry parameter

The asymmetry of the Ly$\alpha$ line profile provides an alternative way of estimating the optical depth of neutral gas in the ISM (e.g. Erb et al. 2014). More symmetric profiles might suggest the presence of an unresolved blue bump blended with the red peak of Ly$\alpha$ (see Section 2.4), or Ly$\alpha$ emission at systemic velocity in cases where Ly$\alpha$ (and LyC) photons escape through a porous ISM with little scattering (e.g. Verhamme et al. 2015). However, quantifying the asymmetry of Ly$\alpha$ profiles is challenging at relatively low-spectral resolutions ($R \lesssim 3000$) and high skewness (Childs & Stanway 2018). After applying a moderate cut-off on the $a_{\text{asym}}$ uncertainty ($\sigma_{a_{\text{asym}}} < 0.07$), we are left with 554 LAEs, or 39 per cent of the total sample. From these LAEs, we select two groups with $a_{\text{asym}} < 0.19$ (potential leakers) and $a_{\text{asym}} > 0.19$ (potential non-leakers), both comprising of 277 objects.

The quantitative comparison between $a_{\text{asym}}$ and $f_{\text{esc}}$ is missing from the literature. If we use FWHM$_{Ly\alpha}$ to estimate the median $f_{\text{esc}}$, we obtain values of 2.0 per cent (low-$a_{\text{asym}}$ subsample) and 2.8 per cent (high-$a_{\text{asym}}$ subsample) that do not comply with our classification – LAEs with more symmetric Ly$\alpha$ profiles are expected to have higher $f_{\text{esc}}$. Furthermore, we find that only 81 objects (29 per cent) from the low-$a_{\text{asym}}$ subsample are part of the low-$v_{\text{sep}}$ and -FWHM$_{Ly\alpha}$ subsamples (Fig. 2, bottom). This could indicate that $a_{\text{asym}}$, on the one hand, and $v_{\text{sep}}$ and FWHM$_{Ly\alpha}$, on the other hand, are sensitive to different geometries of the ISM (namely, the ‘picket-fence’ and the density-bounded geometries) which could both drive LyC escape but not necessarily under the same physical conditions. We discuss different ISM geometries further in Section 5 when examining the stacked spectra of our LAEs.

## 3.4 Connection to the global UV properties

Finally, we compare the global UV properties ($M_{\text{UV}}$, $\beta$) of potential LyC leakers and non-leakers. We find that $M_{\text{UV}}$ becomes fainter with increasing $f_{\text{esc}}$ inferred from equation (2), although the correlation is generally weak. The FWHM$_{Ly\alpha}$ subsamples show the largest variations in $M_{\text{UV}}$, with magnitudes ranging from $M_{\text{UV}} \approx -17.7$ for leakers (FWHM$_{Ly\alpha}$ Q1) to $M_{\text{UV}} \approx -19.0$ for non-leakers (FWHM$_{Ly\alpha}$ Q4). This trend is consistent with the cosmological radiation-hydrodynamical simulations from Rosdahl et al. (2022), who demonstrated that $f_{\text{esc}}$ peaks at $M_{\text{UV}} \approx -17$ (see fig. 9 in their paper). Having deeper data e.g. from lensing field observations would be important to test their predictions for a declining $f_{\text{esc}}$ at magnitudes below $M_{\text{UV}} \approx -17$.

We also find that on average, the LyC-leaker candidates are bluer. For instance, the low- and high-$v_{\text{sep}}$ subsamples have median $\beta = -2.16$ and $\beta = -1.98$, respectively. This result is in agreement with Chisholm et al. (2022) who reported a $6\sigma$ inverse correlation between the $\beta$ slope and $f_{\text{esc}}$ based on the Low-redshift Lyman continuum survey (LzLCS) observations of 89 SFGs at $z \approx 0.3$. We caution, however, that the dynamic range of $\beta$ slopes across our samples of leakers and non-leakers is small ($\Delta\beta \lesssim 0.2$). The lack of significant variations in $\beta$ possibly indicates that the indirect tracers of LyC based on the Ly$\alpha$ line profile are not as sensitive to the dust content as the $\beta$ slope.

# 4 SPECTRAL STACKING

## 4.1 Method

Most of the rest-frame UV lines, both in emission and absorption, are too faint to be detected in the individual spectra of our LAEs. To address this problem, we perform spectral stacking to increase the S/N of the observed spectra.

First, we convert wavelengths from air into vacuum conditions<fnref n="1" /> and shift the spectra to the rest frame. We use systemic redshifts recovered from the Ly$\alpha$ line profiles following the method described in Verhamme et al. (2018, hereafter V18). They propose two diagnostics derived from the spectroscopic observations of LAEs with accurate systemic redshift measurements. In the case of the

<footnote n="1">¹http://www.astro.uu.se/valdwiki/Air-to-vacuum%20conversion</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: Median-stacked spectra of the total sample of MUSE LAEs described in Section 2. Grey shaded regions represent the 1σ noise level. The stacked spectra reveal various nebular emission lines (pink), ISM absorption lines (green), and fine-structure transitions (yellow-green) typically observed in high-redshift and local metal-poor SFGs. Dashed vertical lines mark the rest-frame wavelengths of the spectral lines. Right: Fraction of LAEs included in the stack (black solid line) and redshift coverage (light orange area) as a function of wavelength. The median redshift is shown by the dashed orange line.</figcaption></figure>

low- and high-$v_{\text{sep}}$ subsamples, we use equation (1) of V18 which relates $v_{\text{peak}}^{\text{red}}$ to the Ly$\alpha$ peak separation. For other subsamples, we use the empirical correlation between $v_{\text{peak}}^{\text{red}}$ and the Ly$\alpha$ FWHM (equation 2 of V18). We discuss possible caveats associated with this approach in Section 4.3.

The rest-frame spectra of galaxies can have a varying wavelength sampling depending on their redshift. Therefore, we resample the spectra using the flux-conserving SPECTRES tool (Carnall 2017). We adopt a target sampling of 0.25 Å corresponding to a rest-frame sampling of a spectrum observed by MUSE at z = 4. We then

create the median-stacked spectra and estimate uncertainties of the spectral flux densities using the standard deviation of 200 bootstrap replications. We apply the median statistics instead of the mean or the weighted average to ensure that the stacked spectra are not dominated by the few brightest sources in our sample.

## 4.2 Rest-frame UV lines revealed by the stacked spectra

The median-stacked spectra for the total LAE sample are shown in Fig. 3. We apply a median filter with a window size of 100 Å to
fit the continuum and obtain an average S/N of 3.1 per spectral bin (0.25 Å) in the wavelength range 1250 Å < λ < 1900 Å. This value significantly exceeds a typical S/N of an individual LAE spectrum (≤0.1), demonstrating the efficiency of spectral stacking.

In the same figure, we mark the rest-frame wavelengths of various nebular emission lines, ISM absorption lines, and fine-structure transitions which are commonly observed in the spectra of high-redshift and local metal-poor SFGs (see Feltre et al. 2020 and references therein). We measure the line EWs by fitting Gaussian functions to the continuum-subtracted spectra sliced to ±500 km s⁻¹ regions around the lines. For the fitting procedure, we use the AS-TROPY implementation of the Levenberg-Marquardt optimization algorithm, constraining the FWHM within the range of 50–300 km s⁻¹. If the fit is unsuccessful, we measure the EW non-parametrically by summing up the flux in the same spectral region where the fit is performed. In the case of the [C III]λ1907 + C III]λ1909 doublet (hereafter C III), we use a sum of two Gaussian profiles with a tied peak separation of 2.05 Å. We fit Gaussian functions to the rest of the line doublets separately because their components do not overlap with each other. In the case of the ISM absorption lines, we shift the spectral regions by −200 km s⁻¹ to take into account absorption line velocity offsets due to large-scale gas outflows (e.g. Shapley et al. 2003). Finally, we estimate the EW uncertainties using a Monte Carlo approach by repeating the measurements on the spectra perturbed with noise 1000 times.

Among the high-ionization emission lines, we report a S/N > 3 detection of the C IVλλ1548,1551 resonant doublet (hereafter C IV), collisionally excited O III]λλ1661,1666 and C III] doublets, and the He IIλ1640 line (hereafter He II; the EWs are reported in Table A1). We measure the EW of only the O III]λ1666 component of the doublet (hereafter O III), because the weaker bluer component (O III]λ1661) is not detected in our stacks with S/N > 3. Additionally, we find [Si III]λ1883 + Si III]λ1892 nebular emission in some of the stacks (Section 5). We limit ourselves to the analysis of the [Si III]λ1883 component (hereafter [Si III]) because it always dominates the total flux of the doublet.

While O III], C III], and [Si III] represent exclusively nebular emission, the origins of He II are likely more complex. In particular, strong and dense stellar winds of Wolf–Rayet (W–R) stars give rise to broad He II features observed along with nebular He II (e.g. Nanayakkara et al. 2019). The C IV line also includes other components in addition to nebular emission, i.e. emission from stellar winds in O and B stars characterized by P Cygni profiles, and ISM absorption (e.g. Berg et al. 2018). The S/N of our stacked spectra is insufficient to study the He II and C IV line profiles in detail in order to constrain contributions from different sources of emission (and absorption). Therefore, we report only the total flux of He II and C IV obtained from fitting Gaussian profiles.

Finally, we compare our stacked spectra with the full stack of 220 LAEs from the MUSE HUDF survey shown in fig. 3 of Feltre et al. (2020). Their spectra reveal a similar collection of nebular emission lines (e.g. O III], He II, and C III]) and ISM absorption features (e.g. C IIλ1335 and Si IVλ1403). However, they obtain a higher S/N per spectral pixel (S/N ~ 5 versus S/N ~ 3) despite using a similar wavelength sampling (0.3 Å versus 0.25 Å), likely because their sample only includes the deep MUSE HUDF data. On the other hand, by including the MUSE-Wide LAEs in our sample, we are probing both bright and faint LAEs, thus expanding the parameter space. This enables a more comprehensive analysis of the galaxy properties related to LyC escape.

## 4.3 The impact of redshift uncertainties on measuring EWs

Empirical relations used to recover systemic redshifts provide lower accuracy than the direct estimates based on the shifts of non-resonant lines. If systemic redshifts are measured with high statistical uncertainties, line fluxes are spread over a large wavelength range in the spectral stacks, resulting in a low-S/N and making the assessment of the line properties more difficult. If the spectra are stacked using the median statistics, the line EWs could in addition be biased towards lower values. For example Feltre et al. (2020) stacked copies of idealized spectra with a line centre shifted according to equation (2) of V18 and found that the EW is underestimated by 15–20 per cent. To investigate these effects, we select 24 LAEs with detected (S/N > 3) non-resonant emission lines, i.e. collisionally excited nebular emission lines and the He II line. We compare two spectral stacks obtained using (i) the ‘true’ systemic redshifts based on the non-resonant line velocity offsets with respect to the line rest wavelengths, and (ii) systemic redshift estimates given by equation (2) of V18. We find that if systemic redshifts are recovered using the Ly α line properties (equation 2 of V18), the EW measurements have similar absolute uncertainties but the line profiles are 1.2–1.6 times broader (Fig. 4, top panel) and the line S/Ns are typically lower by a few tens of per cent (Fig. 4, middle panel). The bottom panel in Fig. 4 shows that the EWs of C IV, He II, and C III] are underestimated by 5–20 per cent, in agreement with the simulation results from Feltre et al. (2020). However, we find that the [Si III] EW is underestimated by as much as $44_{-27}^{+39}$ per cent, whilst the O III] EW is higher than the ‘true’ EW by $15_{-42}^{+78}$ per cent. Taking into account high statistical uncertainties of these measurements, larger samples of LAEs with known systemic redshifts are required to provide a more accurate estimate of the EW bias. Nevertheless, we caution that the EWs reported in this paper might in some cases be underestimated by up to ≈20 per cent.

## 4.4 Smaller redshift coverage at longer wavelengths

The effective number of objects used to compute the median stack varies as a function of wavelength since our LAEs are observed at different redshifts (right panels in Fig. 3). Starting from 1422 sources (100 per cent) at the Ly α rest-frame wavelength (1215.67 Å), this number gradually decreases towards longer wavelengths, reaching 721 (51 per cent) at 1900 Å (at λ > 1900 Å the S/N quickly drops down). At the same time, the maximum possible redshift of an object included in the stack decreases from $z_{max}$ = 6.7 at 1215 Å to $z_{max}$ = 3.9 at 1900 Å. Despite having smaller redshift coverage at longer wavelengths, we can expect that e.g. the properties of the C III] line observed in LAEs at 2.9 < z < 3.9 would be representative for LAEs at 4 ≲ z ≲ 6 as well. This assumption is based on the observational evidence that LAEs at 2 ≲ z ≲ 6 share a similar distribution of several fundamental properties including sizes (e.g. Malhotra et al. 2012), UV slopes, Ly α EWs, and scale length parameters (e.g. Santos et al. 2020), and Ly α line profiles corrected for the IGM absorption (Hayes et al. 2021). In addition, the observed LAE luminosity function does not evolve significantly at such redshifts (e.g. Cassata et al. 2011; Herenz et al. 2019). We also compare the full stack and the stack of LAEs at 2.9 < z < 3.9, and find that the C IV/C III] ratio differs by only ~6 per cent. Finally, we note that the median redshift does not change significantly over the same wavelength range, decreasing from $z_{median}$ = 3.9 at 1215 Å to $z_{median}$ = 3.4 at 1900 Å. Throughout the rest of this paper, we assume that the evolution of redshift coverage with wavelength has minimal impact on the observed line properties.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison between the properties of the rest-frame UV emission lines detected in the median-stacked spectra of 24 LAEs with at least one non-resonant line detection at S/N > 3. Before stacking, the individual spectra are shifted to the rest frame using systemic redshifts estimated from (i) velocity offsets of non-resonant lines (white), and (ii) the empirical relation between the Ly α velocity offset and the Ly α FWHM (equation 2 of V18; dark blue). Top to Bottom: Line FWHMs, S/N ratios, and EWs. Vertical error bars mark the 16th and the 84th percentiles of the distributions of the measured quantities. Dashed horizontal lines in the top panel indicate the bounded constraints applied to the line FWHMs in the fitting procedure. The empirical method used to estimate systemic redshifts results, on average, in line broadening, a decrease in the S/N and EW underestimation.</figcaption></figure>

## 5 IONIZING PROPERTIES AND PHYSICAL CONDITIONS IN THE ISM OF LYC-LEAKER CANDIDATES

In this section, we investigate the rest-frame UV spectral properties of our LyC-leaker candidates selected as having narrow, symmetric Ly α profiles with small peak separation (Section 3). We apply the stacking technique described in Section 4 to each group of galaxies. In Fig. 5, we show the median-stacked spectra for both low- and high-v<sub>sep</sub> subsamples, with a focus on the spectral regions around the rest-frame UV lines. We also show the Gaussian profiles fitted to the spectral lines with a S/N > 3 detection. Similar figures but for the groups of different Ly α FWHMs and a<sub>asym</sub> are presented in Appendix B.

The EWs of nebular emission lines, LIS absorption lines and fine-structure emission lines detected with S/N > 3 in at least one of the stacks are listed in Table A1 and plotted in Fig. 6. Apart from Ly α, one of the strongest emission lines we observe is C III, with the total EW ranging from ~3 Å (FWHM<sub>Ly α</sub>, Q4) to ~8 Å (the low-v<sub>sep</sub> sub-

sample; see column 12 in Table A1). We resolve the double-peaked profile of C III] in all the stacks except for FWHM<sub>Ly α</sub>, Q1, for which we fit a single peak. Next, we report a S/N > 3 detection of the C IV line in the low-v<sub>sep</sub> and -a<sub>asym</sub> groups, and the two lowest FWHM<sub>Ly α</sub> quartiles (Q1 and Q2), with the EW as high as 14 ± 2 Å (FWHM<sub>Ly α</sub>, Q1; see column 8 in Table A1). We do not find significant (>2σ) C IV absorption in any of the stacks, suggesting that the contribution from the stellar winds of OB stars to C IV is negligible compared to the nebular emission from ionized gas. Finally, the EWs of He II, O III], and [Si III] take values between ~1 and ~4 Å, with an exception of [Si III] in the high-a<sub>asym</sub> subsample in which the EW is consistent with zero (see columns 9–11 in Table A1). [Si III] is detected with S/N > 3 only in the low-a<sub>asym</sub> stack, whereas He II – also in the two v<sub>sep</sub> stacks, and O III] – in both the low- and high-a<sub>asym</sub> stacks, the high-v<sub>sep</sub> stack and the two intermediate FWHM<sub>Ly α</sub> quartiles.

By measuring the line EWs and the line ratios, we qualitatively compare the properties of the ionizing sources (Sections 5.1–5.3) and the physical conditions in the ISM (Sections 5.4 and 5.5) of the LyC-leaker and non-leaker candidates. We discuss the relationships between the production of ionizing photons, nebular emission, ISM absorption, and Ly α and LyC escape.

### 5.1 Overture: source(s) of the ionizing radiation

Determining the source(s) of ionizing photons is one of the main goals in reionization studies. We argue that AGNs cannot dominate the ionizing radiation from our LAEs. First, we carry out the spectroscopic diagnostics based on the line ratios which trace the presence of an AGN. We find that log(C III]/He II) > 0 and log(O III]/He II) > −0.5 (Fig. 7), in contradiction with the pure AGN photoionization models (Feltre, Charlot & Gutkin 2016, Fig. A1). Second, we detect moderate strengths of C III] (EW ~ 3 – 8 Å; see Table A1, column 12) and C IV (EW ~ 0 – 14 Å; ibid., column 8), which are not compatible with the AGN scenario either (EW > 20 and EW > 12 Å, respectively; Nakajima et al. 2018b). Furthermore, we observe a lack of N v emission (second row in Figs 5, B1, and B2) – a sign of the AGN activity due to a very high photon energy required for ionization of this element (>77.4 eV; e.g. Sobral et al. 2018). Thus, we suggest that young, metal-poor stars have a dominant contribution to the ionizing photon budget of our LAEs, with no exception for the LyC-leaker candidates.

However, additional ionizing sources might be needed to explain the He II EWs (~1–3 Å; see Table A1, column 9). Although we cannot disentangle the nebular and stellar components of He II to examine its possible origins, we note that the current photoionization models are unable to fully reproduce the observed He II emission (e.g. Nanayakkara et al. 2019). Various mechanisms accounting for the missing ionizing photons are proposed to resolve this tension, including emission from Pop III stars, radiative shocks, increased stellar rotation, production of stripped stars, and emission from X-ray binaries and ultraluminous X-ray sources (e.g. Eldridge & Stanway 2012; Izotov, Thuan & Privon 2012; Smith, Götberg & de Mink 2018; Schaerer, Fragos & Izotov 2019; Simmonds, Schaerer & Verhamme 2021).

### 5.2 Nebular emission lines: the ionizing photon budget

In the framework of the density-bounded model of the ISM, high escape fractions of LyC are achieved in a nearly completely ionized medium (e.g. Nakajima & Ouchi 2014). Therefore, efficient LyC escape suggests high ionizing fluxes from the regions of intense star formation. To test this hypothesis, we inspect the nebular emission
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Median-stacked spectra of LAEs with low (<463 km s⁻¹; blue) and high (>463 km s⁻¹; green) Ly α peak separations. Each panel represents a spectral region around one of the rest-frame UV lines (or groups of lines) whose rest-frame wavelengths are indicated by vertical pink (nebular emission), green (ISM absorption), and yellow-green (fine-structure transition) dashed lines. Shaded regions represent the 1σ noise level computed via the bootstrap method. Horizontal dashed lines mark the continuum level. Gaussian profiles fitted to the spectral lines are shown for S/N > 3 detections (thick solid lines).</figcaption></figure>

lines, whose properties provide insight into the ionising photon budget of a galaxy. More specifically, strong nebular emission lines indicate low, subsolar gas-phase metallicities, and high-ionization parameters (e.g. Jaskot & Ravindranath 2016; Steidel et al. 2016; Senchyna et al. 2017) which imply the presence of young stellar populations producing copious amounts of ionizing photons. We find higher EWs of nebular O III], [Si III], and C III] predominantly in the low-$v_{\text{sep}}$ and -FWHM<sub>Ly α</sub> stacks (Fig. 6h1–h2, i1–i2, and j1–j2), supporting the idea that LyC leakers have an increased production rate of ionizing photons compared to non-leakers (see also Naidu et al. 2022). The EWs of Ly α, C IV, and He II show even more compelling trend (ibid., a1–a2, f1–f2, and g1–g2). However, we caution that (i) He II might have a significant stellar component in addition to the nebular (see Section 5.1), (ii) Ly α and C IV are resonant lines, a priori sensitive to the radiation transfer effects which could modulate their EWs alongside with LyC<sub>esc</sub> (see Section 5.5).

The LAE-selected sample gives us the opportunity to use the Ly α line itself to select galaxies with substantially different production rates of ionizing photons. We create four Ly α EW and luminosity quartiles using the data from Kerutt et al. (2022). We find that the EWs of all the nebular emission lines increase dramatically with the Ly α EW (Fig. C1, third column), in good agreement with previous studies (e.g. Shapley et al. 2003; Stark et al. 2014; Feltre et al. 2020). This correlation strongly suggests that high-redshift galaxies with higher Ly α EWs have higher ionizing photon production efficiency (ξ<sub>ion</sub>), implying stellar populations with lower metallicities and younger ages as demonstrated by Maseda et al. (2020) for a sample of continuum-faint LAEs at z ≈ 4–5. The stacks of LAEs with different Ly α luminosities, on the other hand, have nebular emission lines of similar strengths (Fig. C1, fourth column), which could be explained by a variety of metallicites and ages presented at fixed Ly α luminosity.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** EWs of rest-frame UV lines detected in the stacked spectra of candidates for LyC-leakers and non-leakers. Markers show the EWs measured for each group of LAEs. Histograms show the distributions of the Ly α line properties. Left column: Ly α (lavender; a1–a3), ISM absorption lines (sea foam green; b1–d3), and the SiII*λ1533 fine-structure transition (moss green; e1–e3). Right column: C IV (lavender; f1–f3), He II (yellow; g1–g3), and collisionally excited nebular emission lines (pink; h1–j3). The colour-coding of markers follows that of the stacked spectra of respective subsamples in Figs 5, B1, and B2. The distributions of the Ly α line properties are comprised of objects included in the stacks at respective wavelengths; the number of objects is shown in the top right corners of panels.</figcaption></figure>

## 5.3 He II: hardness of the ionizing spectrum

Whether or not a hard ionizing spectrum is a necessary condition for efficient LyC escape remains a subject of debate (Marques-Chaves et al. 2022; Naidu et al. 2022; Schaerer et al. 2022; Enders, Bomans & Wittje 2023). We search for evidence of a harder ionizing spectrum in the stacks of LyC-leaker candidates by inspecting He II emission. We find that the He II EW is typically a few times higher among potential LyC leakers (Fig. 6, g1–g3), suggesting an elevated production rate of He⁺-ionizing photons with energies >54.4 eV. In addition, we report lower C III]/He II and O III]/He II ratios in the low-υ$_{sep}$ and -a$_{asym}$ stacks (Fig. 7) which possibly indicates

a harder ionizing spectrum, given that the ionization energies for C III] and O III] are much lower compared to that of He II (24.4 and 35.1 eV, respectively). Similar results were obtained by Naidu et al. (2022) who studied composite spectra of 26 LAEs at z ~ 2 from the X-SHOOTER Ly α survey (XLS-z2). They detected prominent narrow He II (and C IV) emission in the ‘high escape’ stack (i.e. galaxies with low υ$_{sep}$ or high fraction of the Ly α flux at nearly systemic velocity), and only lower-ionization lines (e.g. C III] and O III]) in the ‘low escape’ stack (i.e. galaxies with high-υ$_{sep}$ and low fraction of the Ly α flux at nearly systemic velocity).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** [C III]λ1907 + C III]λ1909/He IIλ1640 and O III]λ1666/He IIλ1640 flux ratios for the LAE subsamples described in Section 3. The observed line ratios are consistent with the photoionization models of SFGs (e.g. Feltre, Charlot & Gutkin 2016), which rules out a significant contribution of AGNs to the nebular emission detected in the stacks. Background histograms are as in Fig. 6.</figcaption></figure>

At low-redshift ($z \sim 0.3 - 0.4$), Schaerer et al. (2022) studied rest-frame UV spectra of eight LyC emitters and found that the galaxies with $f_{esc} > 10$ per cent show strong He II emission with the rest-frame EWs ranging from 3 to 8 Å. They argued that such high values are primarily due to elevated $\xi_{ion}$, and that galaxies with high $f_{esc}$ do not exhibit exceptionally hard ionizing spectra compared to other galaxies at similar metallicities. Recently, Marques-Chaves et al. (2022) studied a large sample of SFGs from the HST Low-z Lyman continuum survey and found no correlation between He IIλ4686/Hβλ4861 and $f_{esc}$, arriving at the same conclusions. Our results suggest that LyC leakers indeed have higher $\xi_{ion}$ (Section 5.2), however we cannot rule out the possibility that higher $f_{esc}$ also accompanies harder radiation fields, especially in the light of our C III]/He II and O III]/He II measurements.

## 5.4 LIS absorption lines and Si II*: the ISM geometry

Interestingly, the EWs of O III], C III], and even He II (but not Ly α, C IV, and [Si III]) remain constant within $1\sigma$ uncertainties if we select LyC-leaker candidates using $a_{asym}$ (Fig. 6h3, j3, and g3). We suggest that the specifics of the ISM geometry could break the link between nebular emission and $f_{esc}$. Unlike $v_{sep}$ and FWHM$_{Ly\alpha}$ which likely probe the density-bounded geometry, the asymmetry of the Ly α line profile might be more sensitive to the ‘picket-fence’ geometry, i.e. a clumpy ISM with open lines of sight (LOS). In this scenario, high $f_{esc}$ does not require constant support by a strong ionizing radiation field because LyC photons escape freely through empty channels in the ISM.

For the ‘picket-fence’ geometry, a canonical parameter measuring the fraction of LOS covered by dense clumps of neutral hydrogen which block LyC radiation is $f_{cov}$, or the neutral gas covering fraction. Low $f_{cov}$ (or, equivalently, high $f_{esc}$) likely manifests itself in residual fluxes in the cores of LIS absorption lines (e.g. Heckman et al. 2011; Saldana-Lopez et al. 2022). In agreement with such interpretation of absorption line measurements, we find that the EWs of O Iλ1302 + Si IIλ1304 and Si IIλ1527 are lower for the stacks with more symmetric Ly α profiles (Fig. 6b3 and d3). The EW of C IIλ1335 is about -0.9 Å for both stacks (Fig. 6c3), but some additional flux

bluer to the fitted Gaussian is probably lost in the high-$a_{asym}$ stack (see B2). Conversely, the EWs of the LIS absorption lines detected in the $v_{sep}$ and FWHM$_{Ly\alpha}$ stacks do not show any clear trend (Fig. 6b1–b2, c1–c2, d1–d2), suggesting that $f_{cov}$ is loosely coupled to $v_{sep}$ and FWHM$_{Ly\alpha}$.

We caution that the LIS absorption line depths are also affected by the metallicity and the kinematics of the gas, and the infilling by resonant emission (e.g. Vasei et al. 2016). Moreover, the LIS diagnostics is hampered by the low-S/N data. Such effects could explain the lack of the correlation between the LIS absorption line EWs and the Ly α EW (Fig. C1b1–d1) which was reported by previous studies (e.g. Shapley et al. 2003; Jones, Stark & Ellis 2012). This result suggests that a higher sensitivity is required for more rigorous analyses of the LIS absorption lines.

Both the density-bounded and ‘picket-fence’ geometries represent idealized physical models of the ISM. To bring them into a unified picture, we examine the overall distribution of neutral gas in our galaxies. The neutral material off the LOS can be probed with fluorescent non-resonant emission lines such as Si II*λ1533 (hereafter Si II*; e.g. Jaskot & Oey 2014). We detect Si II* emission at a significance level of $>2\sigma$ in all of the $v_{sep}$ and FWHM$_{Ly\alpha}$ stacks, with the highest Si II* EW in the lowest FWHM$_{Ly\alpha}$ quartile (Fig. 6e1–e2). Strong Si II* in the stacks of LyC-leaker candidates suggests large amounts of H I off the LOS, contradicting the ‘classical’ density-bounded scenario which assumes a largely isotropic, almost completely ionized ISM. Ly α (and LyC) photon escape in an anisotropic medium could be attributed to the presence of highly ionized channels of low-optical depth (e.g. Zackrisson, Inoue & Jensen 2013; Rivera-Thorsen et al. 2017). In this case, the relative transparency of these channels determines $v_{sep}$ and FWHM$_{Ly\alpha}$ (Jaskot et al. 2019; McKinney et al. 2019; Kakiichi & Gronke 2021). The distribution of high column density gas, in turn, might be more relevant to $f_{cov}$ in the context of the ‘picket-fence’ geometry (Jaskot et al. 2019), which we associate with $a_{asym}$. In good agreement with this scenario, the Si II* EW in the low-$a_{asym}$ stack is consistent with zero, indicating that a relatively low number of dense molecular clouds imply $f_{cov} < 1$.

This unified picture of the ISM in the context of LyC escape is also supported by high-resolution spectroscopic observations of known LyC leakers at $z \sim 2 - 4$. Most of these galaxies including the Sunburst Arc ($z = 2.4$, Rivera-Thorsen et al. 2017), Ion2 ($z = 3.2$, Vanzella et al. 2015), and Ion3 ($z \approx 4.0$, Vanzella et al. 2018) show complex Ly α profiles, often with a central peak at systemic redshift in addition to conventional red and blue peaks (see also Naidu et al. 2022). The simultaneous observation of Ly α photons that travel through the ISM with and without resonant scattering is naturally explained by the model where the ‘picket-fence’ and the density bound geometries are mixed together as described above.

## 5.5 C IV: the ISM opacity

In addition to being a tracer of hard radiation fields similarly to the He II line, C IV undergoes resonant scattering in the ISM, probing the density of high-ionization gas. This makes nebular C IV emission a proxy for both the production and escape of LyC (e.g. Berg et al. 2019). Thanks to these unique properties, the C IV line has the potential to become one of the most reliable tools for identification of LyC leakers. We note that the same information about the LyC production and escape is also imprinted in Ly α emission, with the only difference being that the Ly α profile is sensitive to the column density of neutral gas, whereas the C IV profile is sensitive to the column density of high-ionization gas. However, the Ly α
<figure><img src="image_10.png" /><figcaption>**Figure 8.** C IV λλ1548,1550 to [C III] λ1907 + C III] λ1909 flux ratio for the LAE subsamples described in Section 3. The observed line ratios measured for the LyC-leaker candidates (purple and dark blue) are consistent with C IV/C III] ratios measured by Schaerer et al. (2022) for a sample of low-redshift (z ~ 0.3 – 0.4) galaxies with *f*<sub>esc</sub> > 10 per cent (C IV/C III]; shaded area). The rapid growth of C IV/C III] with increasing *f*<sub>esc</sub> suggests that this line ratio could serve as an alternative indirect tracer of *f*<sub>esc</sub> at redshifts where Ly α is significantly attenuated by the IGM (z ~ 6). Background histograms are as in Fig. 6.</figcaption></figure>

transmission declines rapidly at z > 6 due to the increasing opacity of the neutral IGM (e.g. Gronke et al. 2021), which emphasizes the importance of C IV as a stand alone probe for $f_{\text{esc}}$ in the EoR.

The LyC diagnostics based on C IV usually involves comparing it with a reference nebular emission line unaffected by resonant scattering. In particular, recent studies have explored the relationship between LyC leakage and the C IV/C III] ratio. Based on the spectroscopic observations of eight local (z ~ 0.3 – 0.4) SFGs with LyC measurements, Schaerer et al. (2022) proposed the following criterion for classifying a galaxy as a strong LyC leaker ($f_{\text{esc}}$ > 0.1): C IV/C III] > 0.75. At higher redshifts (z = 3.1 – 4.6), Saxena et al. (2022b) found that three out of five C IV emitters with $f_{\text{esc}}$ ≈ 0.05 – 0.3 inferred from the LIS absorption lines have similar values of C IV/C III] (≥0.75).

Fig. 8 shows C IV/C III] measured from our stacked spectra of LAEs. We find that C IV/C III] ≥ 0.75 in all the stacks of LyC-leaker candidates (low-$v_{\text{sep}}$, FWHM<sub>Ly α</sub> Q1 and Q2, and low-$a_{\text{asym}}$), in good agreement with the results from Schaerer et al. (2022). Conversely, C IV/C III] is lower than ~ 0.4 in the stacks of potential non-leakers (high-$v_{\text{sep}}$, FWHM<sub>Ly α</sub> Q3 and Q4, and high-$a_{\text{asym}}$). The most striking difference can be seen between the FWHM<sub>Ly α</sub> quartiles, where the C IV to C III] ratio drops down from C IV/C III] = 1.0<sup>+0.3</sup><sub>-0.2</sub> (FWHM<sub>Ly α</sub>, Q1) to zero (FWHM<sub>Ly α</sub>, Q4).

We argue that the observed C IV/C III] ratios in the stacks of potential leakers result not only from more intense ionizing radiation fields, but also from the increased transparency of the ISM. First, we compare the C IV and He II emission in the high-$v_{\text{sep}}$ and -FWHM<sub>Ly α</sub> stacks (Fig. 6f1–f2, g1–g2). We find that the C IV EW is consistent with zero but the He II EW is not, even though the ionization energy of C IV (47.9 eV) is lower than that of He II (54.4 eV). The absence of C IV emission from potential non-leakers thus indicates that C IV experiences a large optical depth from high-ionization gas. Second, we visually inspect the C IV velocity offsets. We find that C IV in the stacks with potential leakers is at the systemic velocity, while in the case of non-leakers, the C IV line is generally redshifted (Figs 5, B1, and B2). The non-zero velocity offsets are a clear signature of resonant scattering of C IV by high-ionization gas, which ultimately leads to a decrease in the observed C IV/C III] ratio. Therefore, we conclude that the ISM opacity effects are likely to play an important role in modulating the observed C IV emission.

## 6 SUMMARY

In this work, we select LyC-leaker candidates from a sample of 1422 MUSE LAEs using theoretically and empirically motivated criteria based on the Ly α peak separation ($v_{\text{sep}}$, FWHM<sub>Ly α</sub>, and asymmetry parameter ($a_{\text{asym}}$). We perform spectral stacking and obtain high-S/N detections of rest-frame UV emission and absorption lines containing valuable information on the ionizing properties and the physical conditions in the ISM of our galaxies. By comparing the stacked spectra of potential LyC leakers and non-leakers, we find the following:

(i) The stacks of LyC-leaker candidates generally show strong nebular emission, revealing an extreme ionization state of the ISM (Section 5.2). This highlights the importance of young, metal-poor stars in creating low column density channels in the ISM through which LyC photons escape.

(ii) Strong He II emission (EW ~ 1 – 3 Å) is typical among LyC-leaker candidates, implying a significant production rate of ionizing photons with energies >54.4 eV (Section 5.3). While elevated $\xi_{\text{ion}}$ is likely to play an important role in boosting the He II EW, lower C III]/He II and O III]/He II in the stacks of potential LyC leakers suggest that stronger He II might also be a consequence of a harder ionizing spectrum possibly associated with higher $f_{\text{esc}}$.

(iii) The LIS absorption lines are generally weaker for the low-$a_{\text{asym}}$ stacks, suggesting that the asymmetry of the Ly α line profile depends on the distribution of high column density gas and hence traces $f_{\text{cov}}$ in the framework of the ‘picket-fence’ model of the ISM. Conversely, the $v_{\text{sep}}$ and FWHM<sub>Ly α</sub> stacks demonstrate a scatter of the LIS absorption line depths, possibly indicating that these Ly α properties are more sensitive to the density-bounded geometry. Despite having a different impact on the Ly α line profile, together, the ‘picket-fence’ and the density-bounded ISM models provide a coherent physical picture of LyC escape (Section 5.4).

(iv) Si II* emission is detected in all the $v_{\text{sep}}$ and FWHM<sub>Ly α</sub> stacks, implying the presence of neutral gas off the LOS even in potential LyC leakers (Section 5.4). This indicates that LyC photons escape from a highly anisotropic ISM through ionized channels surrounded by a high column density medium.

(v) High C IV/C III] ratios (>0.75) are common among the LAEs with potentially high $f_{\text{esc}}$ (≥0.1; Section 5.5). The C IV line profile and the comparison between the EWs of C IV and He II indicate that
elevated C IV/C III] partly arises from the ISM opacity effects related to LyC escape.

Our work suggests that the synergy between the extreme ionizing fields produced by young stellar populations and the patchy ISM riddled with low column density channels creates an ideal physical environment for efficient escape of LyC photons. The C IV/C III] ratio provides the best illustration to our conclusions, being sensitive to both the ionizing photon production efficiency and the opacity of the ISM (Fig. 8). Importantly, unlike the Ly α line, C IV is not affected by changes in the H I content of the IGM which makes it a very promising tool for identification and analysis of LyC-leaking galaxies even in the EoR. Calibrating the relationship between C IV/C III] and $f_{esc}$ at high redshift ($z \gtrsim 3$) will likely become possible in the future thanks to blue-optimized integral-field spectrographs such as BlueMUSE (Richard et al. 2019), which will probe LyC emission from sources at $z \sim 3-4$. Subsequently, C IV/C III] and other spectral properties of SFGs indicative of high $f_{esc}$ could be measured at $z \gtrsim 6$ with JWST and the new generation spectrographs such as MOONS/VLT, PFS/Subaru or HARMONI/ELT, opening new paths for indirect identification of ionizing sources in the EoR.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for the constructive feedback that helped to improve the manuscript. We thank Michael Maseda, Daniel Schaerer, Charlotte Simmonds, and Rashmi Gottumukkala for useful comments and productive discussions. We also thank the organizers and participants of the 24th MUSE Science Busy Week in Leiden. IGK acknowledges an Excellence Master Fellowship granted by the Faculty of Science of the University of Geneva. This work has received funding from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00072, as well as from the Swiss National Science Foundation (SNSF) through project grant number 200020_207349 and SNSF Professorship grant number 190079. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant number 140. This paper is based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes A-0289(B), 095.A-0010(A), 096.A-0045(A), 096.A-0045(B), 094.A-0205, 095.A-0240, 096.A-0090, 097.A-0160, and 098.A-0017. We made extensive use of several open-source software packages and we are thankful to the respective authors for sharing their work: NUMPY (Harris et al. 2020), ASTROPY (Astropy Collaboration 2022), MATPLOTLIB (Hunter 2007), IPYTHON (Perez & Granger 2007), and TOPCAT (Taylor 2005).

## DATA AVAILABILITY

Raw MUSE data cubes are available from the ESO Science Archive Facility.<fnref n="2" /> The MUSE-Wide data from the first data release, which covers 44 fields in GOODS-S/CDFS, are available on the MUSE-Wide project page.<fnref n="3" /> The MUSE HUDF data from the second data release are available via the Advanced MUSE Data products (AMUSED) web interface.<fnref n="4" /> The stacked spectra in a machine-readable format and the spectral line measurements will be shared upon request to the corresponding author.

## REFERENCES

* Astropy Collaboration, 2022, ApJ, 935, 167

* Bacon R. et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, SPIE Conf. Ser. Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. p. 773508

* Bacon R. et al., 2015, A&A, 575, A75

* Bacon R. et al., 2017, A&A, 608, A1

* Bacon R. et al., 2023, A&A, 670, A4

* Bassett R. et al., 2019, MNRAS, 483, 5223

* Berg D. A., Erb D. K., Auger M. W., Pettini M., Brammer G. B., 2018, ApJ, 859, 164

* Berg D. A., Chisholm J., Erb D. K., Pogge R., Henry A., Olivier G. M., 2019, ApJ, 878, L3

* Bouwens R. J., Illingworth G. D., Oesch P. A., Caruana J., Holwerda B., Smit R., Wilkins S., 2015, ApJ, 811, 140

* Bouwens R. J., Illingworth G., Ellis R. S., Oesch P., Stefanon M., 2022, ApJ, 940, 55

* Carnall A. C., 2017, preprint (arXiv:1705.05165)

* Cassata P. et al., 2011, A&A, 525, A143

* Childs H. J. T., Stanway E. R., 2018, MNRAS, 480, 1938

* Chisholm J. et al., 2018, A&A, 616, A30

* Chisholm J., Prochaska J. X., Schaerer D., Gazagnes S., Henry A., 2020, MNRAS, 498, 2554

* Chisholm J. et al., 2022, MNRAS, 517, 5104

* de La Vieuville G. et al., 2020, A&A, 644, A39

* Dijkstra M., Gronke M., Venkatesan A., 2016, ApJ, 828, 71

* Eldridge J. J., Stanway E. R., 2012, MNRAS, 419, 479

* Enders A. U., Bomans D. J., Wittje A., 2023, A&A, 672, A11

* Erb D. K. et al., 2014, ApJ, 795, 33

* Fan X. et al., 2006, AJ, 132, 117

* Feltre A., Charlot S., Gutkin J., 2016, MNRAS, 456, 3354

* Feltre A. et al., 2018, A&A, 617, A62

* Feltre A. et al., 2020, A&A, 641, A118

* Flury S. R. et al., 2022, ApJ, 930, 126

* Gazagnes S., Chisholm J., Schaerer D., Verhamme A., Rigby J. R., Bayliss M., 2018, A&A, 616, A29

* Giavalisco M. et al., 2004, ApJ, 600, L93

* Grogin N. A. et al., 2011, ApJS, 197, 35

* Gronke M. et al., 2021, MNRAS, 508, 3697

* Harris C. R. et al., 2020, Nature, 585, 357

* Hashimoto T. et al., 2015, ApJ, 812, 157

* Hayes M. J., Runnholm A., Gronke M., Scarlata C., 2021, ApJ, 908, 36

* Heckman T. M. et al., 2011, ApJ, 730, 5

* Henry A., Berg D. A., Scarlata C., Verhamme A., Erb D., 2018, ApJ, 855, 96

* Herenz E. C., Wisotzki L., 2017, A&A, 602, A111

* Herenz E. C. et al., 2019, A&A, 621, A107

* Hinshaw G. et al., 2013, ApJS, 208, 19

* Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

* Inoue A. K., Shimizu I., Iwata I., Tanaka M., 2014, MNRAS, 442, 1805

* Izotov Y. I., Thuan T. X., Privon G., 2012, MNRAS, 427, 1229

* Izotov Y. I., Schaerer D., Thuan T. X., Worseck G., Guseva N. G., Orlitová I., Verhamme A., 2016a, MNRAS, 461, 3683

* Izotov Y. I., Orlitová I., Schaerer D., Thuan T. X., Verhamme A., Guseva N. G., Worseck G., 2016b, Nature, 529, 178

* Izotov Y. I., Schaerer D., Worseck G., Guseva N. G., Thuan T. X., Verhamme A., Orlitová I., Fricke K. J., 2018a, MNRAS, 474, 4514

* Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Thuan T. X., Fricke Verhamme A., Orlitová I., 2018b, MNRAS, 478, 4851

* Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Chisholm J., Thuan T. X., Fricke K. J., Verhamme A., 2021, MNRAS, 503, 1734

* Izotov Y. I., Chisholm J., Worseck G., Guseva N. G., Schaerer D., Prochaska J. X., 2022, MNRAS, 515, 2864

* Jaskot A. E., Oey M. S., 2013, ApJ, 766, 91

* Jaskot A. E., Oey M. S., 2014, ApJ, 791, L19

* Jaskot A. E., Ravindranath S., 2016, ApJ, 833, 136

* Jaskot A. E., Dowd T., Oey M. S., Scarlata C., McKinney J., 2019, ApJ, 885, 96

<footnote n="2">²http://archive.eso.org</footnote>

<footnote n="3">³https://musewide.aip.de</footnote>

<footnote n="4">⁴https://amused.univ-lyon1.fr</footnote>
Jones T., Stark D. P., Ellis R. S., 2012, ApJ, 751, 51

Kakiichi K., Gronke M., 2021, ApJ, 908, 30

Kashikawa N. et al., 2011, ApJ, 734, 119

Kerutt J., 2017, Astrophysics Source Code Library, record ascl:1703.0111,

Kerutt J. et al., 2022, A&A, 659, A183

Khaire V., Srianand R., Choudhury T. R., Gaikwad P., 2016, MNRAS, 457,
4051

Kimm T., Blaizot J., Garel T., Michel-Dansac L., Katz H., Rosdahl J.,
Verhamme A., Haehnelt M., 2019, MNRAS, 486, 2215

Koekemoer A. M. et al., 2011, ApJS, 197, 36

Kulkarni G., Worseck G., Hennawi J. F., 2019, MNRAS, 488, 1035

Laursen P., Sommer-Larsen J., Razoumov A. O., 2011, ApJ, 728, 52

Luo B. et al., 2017, ApJS, 228, 2

Madau P., 1995, ApJ, 441, 18

Malhotra S., Rhoads J. E., Finkelstein S. L., Hathi N., Nilsson K., McLinden
E., Pirzkal N., 2012, ApJ, 750, L36

Marchi F. et al., 2017, A&A, 601, A73

Marques-Chaves R., Schaerer D., Álvarez-Márquez J., Colina L., Dessauges-Zavadsky M., Pérez-Fournon I., Saldana-Lopez A., Verhamme A., 2021, MNRAS, 507, 524

Marques-Chaves R. et al., 2022, A&A, 663, L1

Maseda M. V. et al., 2020, MNRAS, 493, 5120

Matthee J. et al., 2021, MNRAS, 505, 1382

Mauerhofer V., Verhamme A., Blaizot J., Garel T., Kimm T., Michel-Dansac L., Rosdahl J., 2021, A&A, 646, A80

McGreer I. D., Mesinger A., D'Odorico V., 2015, MNRAS, 447, 499

McKinney J. H., Jaskot A. E., Oey M. S., Yun M. S., Dowd T., Lowenthal J.
D., 2019, ApJ, 874, 52

Meštrić U., Ryan-Weber E. V., Cooke J., Bassett R., Prichard L. J., Rafelski
M., 2021, MNRAS, 508, 4443

Muzahid S. et al., 2020, MNRAS, 496, 1013

Naidu R. P., Forrest B., Oesch P. A., Tran K.-V. H., Holden B. P., 2018,
MNRAS, 478, 791

Naidu R. P., Tacchella S., Mason C. A., Bose S., Oesch P. A., Conroy C.,
2020, ApJ, 892, 109

Naidu R. P. et al., 2022, MNRAS, 510, 4582

Nakajima K., Ouchi M., 2014, MNRAS, 442, 900

Nakajima K., Fletcher T., Ellis R. S., Robertson B. E., Iwata I., 2018a,
MNRAS, 477, 2098

Nakajima K. et al., 2018b, A&A, 612, A94

Nakajima K., Ellis R. S., Robertson B. E., Tang M., Stark D. P., 2020, ApJ,
889, 161

Nanayakkara T. et al., 2019, A&A, 624, A89

Ouchi M. et al., 2008, ApJS, 176, 301

Ouchi M. et al., 2009, ApJ, 706, 1136

Parsa S., Dunlop J. S., McLure R. J., 2018, MNRAS, 474, 2904

Perez F., Granger B. E., 2007, Comput. Sci. Eng., 9, 21

Planck Collaboration XLVII et al., 2016, A&A, 596, A108

Reddy N. A., Steidel C. C., Pettini M., Bogosavljević M., Shapley A. E.,
2016, ApJ, 828, 108

Richard J. et al., 2019, preprint (arXiv:1906.01657)

Rivera-Thorsen T. E. et al., 2017, A&A, 608, L4

Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJ, 802, L19

Rosdahl J. et al., 2022, MNRAS, 515, 2386

Saldana-Lopez A. et al., 2022, A&A, 663, A59

Santos S. et al., 2020, MNRAS, 493, 141

Saxena A. et al., 2022a, MNRAS, 511, 120

Saxena A. et al., 2022b, MNRAS, 517, 1098

Schaerer D., Izotov Y. I., Nakajima K., Worseck G., Chisholm J., Verhamme
A., Thuan T. X., de Barros S., 2018, A&A, 616, L14

Schaerer D., Fragos T., Izotov Y. I., 2019, A&A, 622, L10

Schaerer D. et al., 2022, A&A, 658, L11

Schmidt K. B. et al., 2019, A&A, 628, A91

Schmidt K. B. et al., 2021, A&A, 654, A80

Scoville N. et al., 2007, ApJS, 172, 1

Senchyna P. et al., 2017, MNRAS, 472, 2608

Shapley A. E., Steidel C. C., Pettini M., Adelberger K. L., 2003, ApJ, 588, 65

Shapley A. E., Steidel C. C., Strom A. L., Bogosavljević M., Reddy N. A.,
Siana B., Mostardi R. E., Rudie G. C., 2016, ApJ, 826, L24

Shen X., Hopkins P. F., Faucher-Giguère C.-A., Alexander D. M., Richards
G. T., Ross N. P., Hickox R. C., 2020, MNRAS, 495, 3252

Simmonds C., Schaerer D., Verhamme A., 2021, A&A, 656, A127

Smith N., Göteborg Y., de Mink S. E., 2018, MNRAS, 475, 772

Sobral D. et al., 2018, MNRAS, 477, 2817

Song M. et al., 2014, ApJ, 791, 3

Stark D. P., Ellis R. S., Chiu K., Ouchi M., Bunker A., 2010, MNRAS, 408, 1628

Stark D. P. et al., 2014, MNRAS, 445, 3200

Steidel C. C., Erb D. K., Shapley A. E., Pettini M., Reddy N., Bogosavljević
M., Rudie G. C., Rakic O., 2010, ApJ, 717, 289

Steidel C. C., Strom A. L., Pettini M., Rudie G. C., Reddy N. A., Trainor R.
F., 2016, ApJ, 826, 159

Steidel C. C., Bogosavljević M., Shapley A. E., Reddy N. A., Rudie G. C.,
Pettini M., Trainor R. F., Strom A. L., 2018, ApJ, 869, 123

Taylor M. B., 2005, in Shopbell P., Britton M., Ebert R. eds, ASP Conf. Ser. Vol. 347, Astronomical Data Analysis Software and Systems XIV, Astron. Soc. Pac., San Francisco, p. 29

Urrutia T. et al., 2019, A&A, 624, A141

Vanzella E. et al., 2015, A&A, 576, A116

Vanzella E. et al., 2018, MNRAS, 476, L15

Vasei K. et al., 2016, ApJ, 831, 38

Verhamme A., Orlitová I., Schaerer D., Hayes M., 2015, A&A, 578, A7

Verhamme A., Orlitová I., Schaerer D., Izotov Y., Worseck G., Thuan T. X.,
Guseva N., 2017, A&A, 597, A13

Verhamme A. et al., 2018, MNRAS, 478, L60 (V18)

Wang B., Heckman T. M., Leitherer C., Alexandroff R., Borthakur S.,
Overzier R. A., 2019, ApJ, 885, 57

Xu X. et al., 2022, ApJ, 933, 202

Zackrisson E., Inoue A. K., Jensen H., 2013, ApJ, 777, 39

Zheng Z.-Y., Wang J.-X., Malhotra S., Rhoads J. E., Finkelstein S. L.,
Finkelstein K., 2014, MNRAS, 439, 1101

APPENDIX A: EQUIVALENT WIDTHS OF REST-FRAME UV LINES DETECTED IN THE STACKED SPECTRA OF MUSE LAES
<table><caption>**Table A1.** EWs of rest-frame UV lines detected in the stacked spectra of MUSE LAEs.</caption><thead><tr><th>Sample<br>(1)</th><th>Σ<br>(2)</th><th>Lyman α<br>(3)</th><th>OI + Si IIλ1302,4<br>(4)</th><th>C IIλ1335<br>(5)</th><th>Si IIλ1527<br>(6)</th><th>Si II*λ1533<br>(7)</th><th>C IVλλ1548,51<br>(8)</th><th>He IIλ1640<br>(9)</th><th>O IIλ1666<br>(10)</th><th>[Si IIλ1883<br>(11)]</th><th>C IIIλλ1907,9<br>(12)</th></tr></thead><tbody><tr><td>All</td><td>1422</td><td>83.6<sup>+1.3</sup><sub>-1.5</sub></td><td>-0.4<sup>+0.1</sup><sub>-0.1</sub></td><td>-0.8<sup>+0.2</sup><sub>-0.2</sub></td><td>-0.6<sup>+0.2</sup><sub>-0.2</sub></td><td>1.2<sup>+0.2</sup><sub>-0.2</sub></td><td>2.2<sup>+0.5</sup><sub>-0.5</sub></td><td>1.1<sup>+0.3</sup><sub>-0.3</sub></td><td>1.9<sup>+0.3</sup><sub>-0.3</sub></td><td>1.0<sup>+0.5</sup><sub>-0.5</sub></td><td>5.7<sup>+0.8</sup><sub>-0.7</sub></td></tr><tr><td>υ<sub>sep</sub> &lt; 463 km s<sup>-1</sup></td><td>165</td><td>100.2<sup>+5.3</sup><sub>-4.7</sub></td><td>-0.7<sup>+0.7</sup><sub>-0.8</sub></td><td>-0.7<sup>+0.8</sup><sub>-0.8</sub></td><td>-1.9<sup>+0.6</sup><sub>-0.6</sub></td><td>1.0<sup>+0.5</sup><sub>-0.5</sub></td><td>4.0<sup>+1.2</sup><sub>-1.1</sub></td><td>2.9<sup>+1.0</sup><sub>-1.1</sub></td><td>2.3<sup>+0.9</sup><sub>-0.8</sub></td><td>0.9<sup>+0.8</sup><sub>-0.4</sub></td><td>8.1<sup>+2.3</sup><sub>-2.2</sub></td></tr><tr><td>υ<sub>sep</sub> &gt; 463 km s<sup>-1</sup></td><td>164</td><td>61.9<sup>+2.4</sup><sub>-2.4</sub></td><td>-0.8<sup>+0.3</sup><sub>-0.3</sub></td><td>-1.2<sup>+0.3</sup><sub>-0.2</sub></td><td>0.7<sup>+0.5</sup><sub>-0.5</sub></td><td>1.4<sup>+0.3</sup><sub>-0.3</sub></td><td>0.6<sup>+0.8</sup><sub>-0.7</sub></td><td>1.0<sup>+0.3</sup><sub>-0.3</sub></td><td>1.3<sup>+0.4</sup><sub>-0.4</sub></td><td>1.6<sup>+0.7</sup><sub>-0.6</sub></td><td>5.1<sup>+0.9</sup><sub>-1.0</sub></td></tr><tr><td>FWHM<sub>Lyα</sub> (Q1)</td><td>241</td><td>162.5<sup>+9.3</sup><sub>-10.1</sub></td><td>0.5<sup>+1.3</sup><sub>-1.4</sub></td><td>-1.8<sup>+0.7</sup><sub>-0.9</sub></td><td>-1.6<sup>+0.5</sup><sub>-0.7</sub></td><td>5.5<sup>+1.9</sup><sub>-1.8</sub></td><td>14.1<sup>+2.5</sup><sub>-2.3</sub></td><td>3.4<sup>+1.1</sup><sub>-1.2</sub></td><td>3.8<sup>+1.8</sup><sub>-1.0</sub></td><td>3.1<sup>+1.3</sup><sub>-1.1</sub></td><td>6.9<sup>+1.8</sup><sub>-1.4</sub></td></tr><tr><td>FWHM<sub>Lyα</sub> (Q2)</td><td>241</td><td>144.6<sup>+5.5</sup><sub>-5.5</sub></td><td>0.7<sup>+1.0</sup><sub>-0.9</sub></td><td>0.3<sup>+1.0</sup><sub>-1.0</sub></td><td>-1.1<sup>+0.4</sup><sub>-0.6</sub></td><td>1.9<sup>+0.7</sup><sub>-0.6</sub></td><td>6.4<sup>+1.4</sup><sub>-1.4</sub></td><td>1.3<sup>+0.6</sup><sub>-0.5</sub></td><td>5.1<sup>+1.0</sup><sub>-1.0</sub></td><td>1.1<sup>+0.6</sup><sub>-0.5</sub></td><td>7.1<sup>+1.5</sup><sub>-2.4</sub></td></tr><tr><td>FWHM<sub>Lyα</sub> (Q3)</td><td>241</td><td>105.4<sup>+5.2</sup><sub>-5.6</sub></td><td>-0.6<sup>+0.3</sup><sub>-0.2</sub></td><td>-1.2<sup>+0.5</sup><sub>-0.4</sub></td><td>-1.8<sup>+0.8</sup><sub>-0.8</sub></td><td>1.3<sup>+0.5</sup><sub>-0.5</sub></td><td>1.6<sup>+0.6</sup><sub>-0.6</sub></td><td>1.6<sup>+0.6</sup><sub>-0.5</sub></td><td>2.1<sup>+0.7</sup><sub>-0.6</sub></td><td>2.4<sup>+1.0</sup><sub>-0.9</sub></td><td>7.8<sup>+1.3</sup><sub>-1.4</sub></td></tr><tr><td>FWHM<sub>Lyα</sub> (Q4)</td><td>241</td><td>51.2<sup>+1.8</sup><sub>-1.9</sub></td><td>-0.8<sup>+0.2</sup><sub>-0.2</sub></td><td>-0.9<sup>+0.2</sup><sub>-0.2</sub></td><td>-1.0<sup>+0.3</sup><sub>-0.3</sub></td><td>0.6<sup>+0.4</sup><sub>-0.4</sub></td><td>-0.3<sup>+0.5</sup><sub>-0.5</sub></td><td>0.6<sup>+0.3</sup><sub>-0.3</sub></td><td>0.7<sup>+0.3</sup><sub>-0.2</sub></td><td>1.2<sup>+0.5</sup><sub>-0.4</sub></td><td>2.8<sup>+0.7</sup><sub>-0.9</sub></td></tr><tr><td>a<sub>asym</sub> &lt; 0.19</td><td>277</td><td>87.3<sup>+2.9</sup><sub>-2.8</sub></td><td>-0.2<sup>+0.4</sup><sub>-0.3</sub></td><td>-0.9<sup>+0.2</sup><sub>-0.9</sub></td><td>-0.7<sup>+0.4</sup><sub>-0.4</sub></td><td>0.3<sup>+0.4</sup><sub>-0.4</sub></td><td>2.6<sup>+0.6</sup><sub>-0.6</sub></td><td>1.4<sup>+0.3</sup><sub>-0.4</sub></td><td>2.4<sup>+0.5</sup><sub>-0.5</sub></td><td>1.5<sup>+0.4</sup><sub>-0.5</sub></td><td>5.4<sup>+0.7</sup><sub>-1.0</sub></td></tr><tr><td>a<sub>asym</sub> &gt; 0.19</td><td>277</td><td>68.8<sup>+2.1</sup><sub>-2.2</sub></td><td>-0.6<sup>+0.2</sup><sub>-0.2</sub></td><td>-0.9<sup>+0.3</sup><sub>-0.3</sub></td><td>-1.1<sup>+0.3</sup><sub>-1.0</sub></td><td>1.2<sup>+0.5</sup><sub>-0.4</sub></td><td>1.2<sup>+0.5</sup><sub>-0.6</sub></td><td>0.8<sup>+0.3</sup><sub>-0.3</sub></td><td>2.2<sup>+0.6</sup><sub>-0.8</sub></td><td>-0.8<sup>+1.0</sup><sub>-0.9</sub></td><td>5.7<sup>+1.5</sup><sub>-1.2</sub></td></tr><tr><td>EW<sub>Lyα</sub> (Q1)</td><td>262</td><td>34.3<sup>+1.5</sup><sub>-1.5</sub></td><td>-0.3<sup>+0.4</sup><sub>-0.4</sub></td><td>-1.0<sup>+0.2</sup><sub>-0.2</sub></td><td>-0.7<sup>+0.3</sup><sub>-0.3</sub></td><td>1.6<sup>+0.5</sup><sub>-0.5</sub></td><td>0.7<sup>+0.4</sup><sub>-0.5</sub></td><td>0.4<sup>+0.2</sup><sub>-0.2</sub></td><td>0.6<sup>+0.3</sup><sub>-0.3</sub></td><td>0.8<sup>+0.6</sup><sub>-0.4</sub></td><td>1.2<sup>+1.9</sup><sub>-0.8</sub></td></tr><tr><td>EW<sub>Lyα</sub> (Q2)</td><td>262</td><td>70.2<sup>+2.5</sup><sub>-2.8</sub></td><td>-0.7<sup>+0.3</sup><sub>-0.3</sub></td><td>-0.7<sup>+0.3</sup><sub>-0.3</sub></td><td>-1.2<sup>+0.4</sup><sub>-1.0</sub></td><td>0.7<sup>+0.3</sup><sub>-0.3</sub></td><td>1.9<sup>+0.4</sup><sub>-0.7</sub></td><td>0.2<sup>+0.7</sup><sub>-0.7</sub></td><td>2.3<sup>+0.6</sup><sub>-0.7</sub></td><td>1.3<sup>+0.4</sup><sub>-0.4</sub></td><td>3.4<sup>+0.7</sup><sub>-0.7</sub></td></tr><tr><td>EW<sub>Lyα</sub> (Q3)</td><td>262</td><td>101.3<sup>+3.8</sup><sub>-3.7</sub></td><td>-0.9<sup>+0.3</sup><sub>-0.4</sub></td><td>-1.3<sup>+0.5</sup><sub>-0.5</sub></td><td>-1.3<sup>+0.6</sup><sub>-0.6</sub></td><td>2.0<sup>+0.7</sup><sub>-0.7</sub></td><td>2.2<sup>+1.2</sup><sub>-0.7</sub></td><td>0.5<sup>+0.4</sup><sub>-0.3</sub></td><td>2.5<sup>+0.6</sup><sub>-0.6</sub></td><td>0.9<sup>+0.5</sup><sub>-0.4</sub></td><td>9.9<sup>+2.0</sup><sub>-2.2</sub></td></tr><tr><td>EW<sub>Lyα</sub> (Q4)</td><td>262</td><td>176.6<sup>+4.1</sup><sub>-4.7</sub></td><td>1.5<sup>+1.0</sup><sub>-1.0</sub></td><td>-1.5<sup>+0.5</sup><sub>-0.6</sub></td><td>-0.9<sup>+1.1</sup><sub>-1.2</sub></td><td>0.1<sup>+1.1</sup><sub>-1.0</sub></td><td>6.4<sup>+1.8</sup><sub>-1.8</sub></td><td>1.5<sup>+0.7</sup><sub>-0.6</sub></td><td>4.1<sup>+1.4</sup><sub>-1.2</sub></td><td>3.9<sup>+3.4</sup><sub>-1.6</sub></td><td>25.0<sup>+4.9</sup><sub>-5.9</sub></td></tr><tr><td>L<sub>Lyα</sub> (Q1)</td><td>356</td><td>91.1<sup>+2.2</sup><sub>-1.9</sub></td><td>-1.7<sup>+0.8</sup><sub>-0.9</sub></td><td>1.3<sup>+0.9</sup><sub>-0.9</sub></td><td>0.9<sup>+0.9</sup><sub>-1.0</sub></td><td>2.1<sup>+0.6</sup><sub>-0.6</sub></td><td>3.5<sup>+1.0</sup><sub>-1.1</sub></td><td>1.0<sup>+0.7</sup><sub>-0.7</sub></td><td>2.4<sup>+0.6</sup><sub>-0.6</sub></td><td>1.7<sup>+1.2</sup><sub>-0.9</sub></td><td>4.9<sup>+1.3</sup><sub>-1.3</sub></td></tr><tr><td>L<sub>Lyα</sub> (Q2)</td><td>355</td><td>86.5<sup>+2.4</sup><sub>-2.4</sub></td><td>-0.8<sup>+0.2</sup><sub>-0.2</sub></td><td>-1.0<sup>+0.3</sup><sub>-0.5</sub></td><td>-1.1<sup>+0.7</sup><sub>-0.7</sub></td><td>0.7<sup>+0.4</sup><sub>-0.4</sub></td><td>1.6<sup>+0.6</sup><sub>-0.6</sub></td><td>0.6<sup>+0.5</sup><sub>-0.2</sub></td><td>2.6<sup>+0.6</sup><sub>-0.6</sub></td><td>0.8<sup>+0.4</sup><sub>-0.4</sub></td><td>5.8<sup>+1.2</sup><sub>-1.2</sub></td></tr><tr><td>L<sub>Lyα</sub> (Q3)</td><td>355</td><td>86.2<sup>+2.3</sup><sub>-2.0</sub></td><td>-0.9<sup>+0.4</sup><sub>-0.3</sub></td><td>-1.0<sup>+0.3</sup><sub>-0.3</sub></td><td>-0.7<sup>+0.2</sup><sub>-0.4</sub></td><td>1.5<sup>+0.5</sup><sub>-0.5</sub></td><td>3.2<sup>+1.0</sup><sub>-0.8</sub></td><td>0.8<sup>+0.2</sup><sub>-0.8</sub></td><td>1.1<sup>+0.6</sup><sub>-0.5</sub></td><td>2.5<sup>+1.3</sup><sub>-0.9</sub></td><td>7.3<sup>+0.8</sup><sub>-1.1</sub></td></tr><tr><td>L<sub>Lyα</sub> (Q4)</td><td>355</td><td>91.6<sup>+2.1</sup><sub>-2.1</sub></td><td>0.2<sup>+0.4</sup><sub>-0.5</sub></td><td>-0.8<sup>+0.3</sup><sub>-0.5</sub></td><td>-1.2<sup>+0.3</sup><sub>-0.4</sub></td><td>1.0<sup>+0.4</sup><sub>-0.4</sub></td><td>2.3<sup>+0.6</sup><sub>-0.6</sub></td><td>1.2<sup>+0.5</sup><sub>-0.5</sub></td><td>2.4<sup>+0.4</sup><sub>-0.4</sub></td><td>1.1<sup>+0.6</sup><sub>-0.5</sub></td><td>4.4<sup>+1.2</sup><sub>-1.4</sub></td></tr></tbody></table>

Notes. Columns are as follows: (1) sample description, (2) sample size, (3) Ly α EW, (4–6) EWs of ISM absorption lines (O IIλ1302 + Si IIλ1304, C IIλ1335, and Si IIλ1527), (7) EW of the Si II*λ1533 fine-structure emission line, (8–12) EWs of high-ionization nebular emission lines (C IVλλ1548,1551, He IIλ1640, O IIIλ1666, [Si IIλ1883, and [C IIIλ1907 + C IIIλ1909]). All EWs are given for the rest-frame spectra and expressed in Å.
# APPENDIX B: STACKED SPECTRA OF POTENTIAL LYC LEAKERS AND NON-LEAKERS

## B1 Subsamples of different Ly α FWHMs

## B2 Subsamples of different Ly α asymmetry parameters

<figure><img src="image_11.png" /><figcaption>Figure B1. Similar to Fig. 5, but for the subsamples of LAEs with different Ly α FWHMs.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure B1. continued</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure B2.** Similar to Fig. 5, but for the subsamples of LAEs with low ($<$0.19; blue) and high ($>$0.19; green) Ly$\alpha$ asymmetry parameters.</figcaption></figure>
APPENDIX C: EQUIVALENT WIDTHS OF REST-FRAME UV LINES FOR DIFFERENT LY α EWS AND LUMINOSITIES

<figure><img src="image_14.png" /><figcaption>**Figure C1.** Similar to Fig. 6, but for the subsamples of LAEs with different Ly α EWS and luminosities.</figcaption></figure>

This paper has been typeset from a TeX/LA TeX file prepared by the author.