Linking UV spectral properties of MUSE Ly-alpha emitters at z>3 to Lyman continuum escape

The physical conditions giving rise to high escape fractions of ionizing radiation (LyC $f_{\rm{esc}}$) in star-forming galaxies - most likely protagonists of cosmic reionization - are not yet fully understood. Using the VLT/MUSE observations of ~1400 Ly$\alpha$ 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$\alpha$ profiles. We find that the stacks of potential LyC leakers, i.e. galaxies with narrow, symmetric Ly$\alpha$ profiles with small peak separation, generally show (i) strong nebular OIII]1666, [SiIII]1883, and [CIII]1907+CIII]1909 emission, indicating a high-ionization state of the interstellar medium (ISM); (ii) high equivalent widths of HeII1640 (~1-3 A), suggesting the presence of hard ionizing radiation fields; (iii) SiII*1533 emission, revealing substantial amounts of neutral hydrogen off the line of sight; (iv) high CIV1548,1550 to [CIII]1907+CIII]1909 ratios (CIV/CIII]>0.75), signalling the presence of low column density channels in the ISM. In contrast, the stacks with broad, asymmetric Ly$\alpha$ profiles with large peak separation show weak nebular emission lines, low HeII1640 equivalent widths (<1 A), and low CIV/CIII] (<0.25), implying low-ionization states and high-neutral hydrogen column densities. Our results suggest that CIV/CIII] 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.


INTRODUCTION
Ultraviolet and X-ray radiation emitted from the first stars and galaxies likely ionised the neutral hydrogen in the intergalactic medium (IGM) between  ≈ 12 (e.g.Hinshaw et al. 2013;Planck Collaboration et al. 2016) and  ≈ 6 (e.g.Fan et al. 2006;McGreer et al. 2015).This last phase transition of the Universe is known as the Epoch of Reionisation (EoR).The hydrogen reionisation could have been potentially powered by active galactic nuclei (AGNs) due to their brightness and high escape fractions of ionising radiation (  esc ).However, the AGN number density is likely to be insufficient to maintain this process at  ≳ 6 (e.g Parsa et al. 2018;Kulkarni et al. 2019;Shen et al. 2020).It is currently believed that the most promising candidates for the sources of the ionising photons are young massive ★ E-mail: im.kramarenko@gmail.comstars in star-forming galaxies (SFGs), although their contribution is still unclear, primarily due to the uncertainties associated with  esc .
Recent studies have shown that  esc in SFGs at  ≳ 6 has to be on average at least 10-20 per cent to successfully reproduce the observed constraints on the reionisation history (e.g.Ouchi et al. 2009;Robertson et al. 2015;Bouwens et al. 2015;Khaire et al. 2016;Naidu et al. 2020).Direct measurements of  esc , however, are severely hampered at such redshifts because ionising radiation, or Lyman continuum (LyC,  < 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. 2016aIzotov et al. ,b, 2018a,b;,b;Schaerer et al. 2018;Wang et al. 2019;Flury et al. 2022;Xu et al. 2022) and at intermediate redshifts (1 ≲  ≲ 4; e.g.Vanzella et al. 2015;Shapley et al. 2016;Vanzella et al. 2018;Rivera-Thorsen et al. 2017;Marques-Chaves et al. 2021;Saxena et al. 2022a), often 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- observations help to develop indirect tracers of  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 32 ;Jaskot & Oey 2013;Nakajima & Ouchi 2014).Recent observations suggest, however, that the correlation between O 32 and  esc is weak (e.g.Izotov et al. 2018b;Naidu et al. 2018), 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  esc might be an ensemble of UV low-ionisation state (LIS) absorption lines (e.g.Si ii1260, O i1302 and C ii1334).Observationally, the neutral gas covering fraction (  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;Gazagnes et al. 2018;Chisholm et al. 2018;Saldana-Lopez et al. 2022), although simulations demonstrate a substantial scatter in this relation (Mauerhofer et al. 2021).Alternatively, the ISM conditions driving LyC escape could be probed with resonant nebular emission lines sensitive to the radiation transfer effects.These include the Mg ii2796,2803 doublet (Feltre et al. 2018;Henry et al. 2018;Chisholm et al. 2020), the C iv1548,1550 doublet (Schaerer et al. 2022;Saxena et al. 2022b) and the Lymanalpha (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) ionisation-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 ( red peak ≲ 150 km s −1 ), whereas in the second case the Ly line is at systemic redshift.They proposed to use small Ly peak separation ( sep ≲ 300 km s −1 ) for doublepeaked Ly profiles as an indicator of low neutral hydrogen column density and thus potentially high  esc .Dĳkstra et al. (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  esc typically have a narrow Ly line with small peak separation, confirming the theoretical expectations (Izotov et al. 2016a;Verhamme et al. 2017;Izotov et al. 2018bIzotov et al. , 2021Izotov et al. , 2022;;Flury et al. 2022).
These results demonstrate that the Ly line profile could provide robust predictions for  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 Lyman-alpha emitters (LAEs) at  ≈ 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 <  < 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  esc .In section 4, we stack the individual spectra of our galaxies to obtain high S/N detections of the restframe UV lines, namely, high-ionisation nebular emission lines (e.g.C iv1548,1550, He ii1640 and [C iii]1907+C iii]1909), ISM absorption lines (e.g.O i1302+Si ii1304 and C ii1334) and finestructure 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 ionising properties of galaxies and the physical conditions in their ISM indicative of high  esc .In section 6, we present the summary of our main findings and discuss the implications of our work for reionisation studies.

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 integralfield 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 2 ), high resolving power (ranging from  ≈ 1800 in the blue to ≈ 4000 in the red) and a large simultaneous spectral range (4750-9350 Å; Bacon et al. 2015).The GTO programs take full advantage of these unique spectroscopic capabilities of the MUSE instrument, enabling studies of large, un-targeted samples of faint (M UV down to ≈ −16) LAEs at  > 3.
Being different in depths and sizes, MUSE-Wide and MUSE HUDF effectively complement each other, probing the Ly-bright (L Ly ≳ L * ) and faint (L Ly ≲ 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 2 fields at one hour observation time, among which 60 are in the CANDELS (Koekemoer et al. 2011;Grogin 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 ten hours depth in HUDF, completing the 100 fields of MUSE-Wide, plus a single UDF-10 field with 31 hours 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).

Sample selection
We select MUSE LAEs for our analysis as follows.First, we only consider LAEs with a confidence C ≥ 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 catalogs (Luo et al. 2017).
Our final sample amounts to 1422 LAEs in the redshift range 2.92 <  < 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  > 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  > 3.9).We find that the evolution of redshift coverage with wavelength does not affect our conclusions (see section 4).

Sample properties
The distributions of the observed properties of our LAEs are shown in Figure 1, with data taken from Kerutt et al. (2022).The redshifts (Figure 1, panel a) 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 ( red peak ), 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 Figure 1, panel b, the absolute UV magnitudes of our LAEs span the range of −21.9 < M UV < −15.6, with a median value of −18.1.Thanks to the inclusion of the MUSE HUDF data, our sample probes ∼ 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 UV ∼ −16 or even fainter, down to M UV ∼ −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 () is shown in Figure 1, panel 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  = −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).

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 Figure 1, panels d-f.We find that our galaxies exhibit a large variety of Ly profiles in terms of the number of peaks, width, and skewness (see  also Figure 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 ≈ 100 km s −1 to ≈ 1200 km s −1 , with a median value of ≈ 470 km s −1 (Figure 1, panel d).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 Figure 1).
There might be several explanations for the lack of the doublepeaked 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: where  is the amplitude,  0 is the wavelength of the red peak,  asym is the asymmetric dispersion and  0 is the continuum level.The asymmetric dispersion is described by  asym =  asym (− 0 )+, where  asym is the asymmetry parameter and  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 (Figure 1, panel f).Negative values of  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 km s −1 to ≈ 440 km s −1 , with a median value of ≈ 220 km s −1 (Figure 1, panel e).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).

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  esc for each resulting group of LAEs wherever possible.At the end of this section, we examine the global UV properties (M UV , ) 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. a 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 (section 4 and section 5).

Ly𝛼 peak separation
The Ly peak separation is considered as one of the most reliable tracers of  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 subsection 2.4) prevents us from having more than two groups selected based on the peak separation.After excluding 40 objects with From this sample, we select two groups of roughly the same size with LAEs having  sep < 463 km s −1 (165 objects; potential leakers) and  sep > 463 km s −1 (164 objects; potential non-leakers), respectively.Izotov et al. (2018b) collected  sep measurements for a sample of low-redshift LyC-leaking galaxies and found the following empirical relationship between  esc and  sep : where  sep is in km s −1 .We use this equation to estimate the expected median LyC  esc and obtain values of 6.0 per cent and 1.1 per cent for the low-and high- sep subsamples, respectively.Therefore, our classification of a galaxy as a LyC leaker based on the Ly peak separation allows lower values of  esc compared to the thresholds usually adopted in other studies (  esc > 20 per cent, Naidu et al. 2022; or  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  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 et al. 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.

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  esc .We select 964 LAEs with  FWHM Ly < 100 km s −1 and split them into four FWHM Ly quartiles (Q1-Q4) of 241 LAEs each as follows: and (Q4) FWHM Ly > 260 km s −1 .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- sep subsample (Figure 2, top), reflecting the fact that FWHM Ly correlates with  sep since both depend on the neutral hydrogen column density as discussed above.
Expanding shell models predict that FWHM 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  esc for the FWHM Ly subsamples.We caution that the relation between FWHM Ly and  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 star-forming galaxies 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 HI ∼ 10 13 cm −2 ).
We find that the expected median  esc ranges from 1.1 per cent in the highest FWHM Ly quartile (Q4) to 15 per cent in the lowest FWHM Ly quartile (Q1).Therefore, the strongest leakers selected by the Ly FWHM are likely to have higher escape fractions than the LAEs with the smallest peak separations (  esc ≈ 6 per cent).Nonetheless, even the lowest FWHM Ly quartile (Q1) has a median  esc smaller than the ones observed in the individual LyC leakers at 2 ≲  ≲ 4 (≳ 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).

Ly𝛼 asymmetry parameter
The asymmetry of the Ly line profile provides an alternative way of estimating the optical depth of neutral gas in the ISM (e.g. 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.

Connection to the global UV properties
Finally, we compare the global UV properties (M UV , ) of potential LyC leakers and non-leakers.We find that M UV becomes fainter with increasing  esc inferred from Equation 2, although the correlation is generally weak.The FWHM Ly subsamples show the largest variations in M UV , with magnitudes ranging from M UV ≈ −17.7 for leakers (FWHM Ly Q1) to M UV ≈ −19.0 for non-leakers (FWHM Ly Q4).This trend is consistent with the cosmological radiation-hydrodynamical simulations from Rosdahl et al. (2022), who demonstrated that  esc peaks at M UV ≈ −17 (see Figure 9 in their paper).Having deeper data e.g. from lensing field observations would be important to test their predictions for a declining  esc at magnitudes below M UV ≈ −17.
We also find that on average, the LyC-leaker candidates are bluer.For instance, the low-and high- sep subsamples have median  = −2.16 and  = −1.98,respectively.This result is in agreement with Chisholm et al. (2022) who reported a 6 inverse correlation between the  slope and  esc based on the Low-redshift Lyman Continuum Survey (LzLCS) observations of 89 SFGs at  ≈ 0.3.We caution, however, that the dynamic range of  slopes across our samples of leakers and non-leakers is small (Δ ≲ 0.2).The lack of significant variations in  possibly indicates that the indirect tracers of LyC based on the Ly line profile are not as sensitive to the dust content as the  slope.

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 1 and shift the spectra to the rest-frame.We use systemic redshifts recovered from the Ly 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 low-and high- sep subsamples, we use Eq. ( 1) of V18 which relates  red peak to the Ly peak separation.For other subsamples, we use the empirical correlation between  red peak and the Ly FWHM (Eq.( 2) of V18).We discuss possible caveats associated with this approach in subsection 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  = 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.

Rest-frame UV lines revealed by the stacked spectra
The median-stacked spectra for the total LAE sample are shown in Figure 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 highredshift and local metal-poor star-forming galaxies (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 −1 regions around the lines.For the fitting procedure, we use the astropy implementation of the the Levenberg-Marquardt optimisation algorithm, constraining the FWHM within the range of 50 − 300 km s −1 .If the fit is unsuccessful, we measure the EW nonparametrically 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 −1 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-ionisation emission lines, we report a S/N > 3 detection of the C iv1548,1551 resonant doublet (hereafter C iv), collisionally excited O iii]1661,1666 and C iii] doublets, and the He ii1640 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 1 http://www.astro.uu.se/valdwiki/Air-to-vacuum% 20conversion [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 characterised 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 Figure 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 ii1335 and Si iv1403).However, they obtain a higher S/N per spectral pixel (S/N ∼ 5 vs. S/N ∼ 3) despite using a similar wavelength sampling (0.3 Å vs. 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.

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 idealised spectra with a line center shifted according to Eq. ( 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 Eq. (2) of V18.We find that if systemic redshifts are recovered using the Ly line properties (Eq.(2) of V18), the EW measurements have similar absolute uncertainties but the line profiles are 1.2-1.6 times broader (Figure 4, top panel) and the line S/N-s are typically lower by a few tens of per cent (Figure 4, middle panel).The bottom panel in Figure 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 +39 −27 per cent, whilst the O iii] EW is higher than the "true" EW by 15 +78 −42 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.

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 Figure 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  max = 6.7 at 1215 Å to  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 <  < 3.9 would be representative for LAEs at 4 ≲  ≲ 6 as well.This assumption is based on the observational evidence that LAEs at 2 ≲  ≲ 6 share a similar distribution of several fundamental properties including sizes (e.g.Comparison between the properties of the rest-frame UV emission lines detected in the median-stacked spectra of 24 LAEs with at least one nonresonant 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 (Eq.(2) of V18; dark blue).Top to Bottom: Line FWHMs, S/N ratios and EWs.Vertical error bars mark the 16-th and the 84-th 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.
eters (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 <  < 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  median = 3.9 at 1215 Å to  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.

IONISING 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 Figure 5, we show the median-stacked spectra for both low-and high- sep 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  asym are presented in the appendix (section 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 Figure 6.Apart from Ly, one of the strongest emission lines we observe is C iii], with the total EW ranging from ∼ 3 Å (FWHM Ly , Q4) to ∼ 8 Å (the low- sep subsample; see column 12 in Table A1).We resolve the double-peaked profile of C iii] in all the stacks except for FWHM Ly , Q1, for which we fit a single peak.Next, we report a S/N > 3 detection of the C iv line in the low- sep and - asym groups, and the two lowest FWHM Ly quartiles (Q1 and Q2), with the EW as high as 14 ± 2 Å (FWHM Ly , 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 ionised 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- asym 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- asym stack, whereas He ii -also in the two  sep stacks, and O iii]in both the low-and high- asym stacks, the high- sep stack and the two intermediate FWHM Ly quartiles.
By measuring the line EWs and the line ratios, we qualitatively compare the properties of the ionising sources (subsection 5.1subsection 5.3) and the physical conditions in the ISM (subsection 5.4 and subsection 5.5) of the LyC-leaker and non-leaker candidates.We discuss the relationships between the production of ionising photons, nebular emission, ISM absorption, and Ly and LyC escape.

Overture: source(s) of the ionising radiation
Determining the source(s) of ionising photons is one of the main goals in reionisation studies.We argue that AGNs cannot dominate the ionising 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 (Figure 7), in contradiction with the pure AGN photoionisation models (Feltre et al. 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 Figure 5, Figure B1 and Figure B2) -a sign of the AGN activity due to a very high photon energy required for ionisation 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 ionising photon budget of our LAEs, with no exception for the LyC-leaker candidates.
However, additional ionising 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 photoionisation models are unable to fully reproduce the observed He ii emission (e.g.Nanayakkara et al. 2019).Various mechanisms accounting for the missing ionising 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 ultra-luminous X-ray sources (e.g.Izotov et al. 2012;Eldridge & Stanway 2012;Smith et al. 2018;Schaerer et al. 2019;Simmonds et al. 2021).

Nebular emission lines: the ionising photon budget
In the framework of the density-bounded model of the ISM, high escape fractions of LyC are achieved in a nearly completely ionised medium (e.g.Nakajima & Ouchi 2014).Therefore, efficient LyC escape suggests high ionising fluxes from the regions of intense star formation.To test this hypothesis, we inspect the nebular emission lines, whose properties provide insight into the ionising photon budget of a galaxy.More specifically, strong nebular emission lines indicate low, sub-solar gas-phase metallicities and high ionisation parameters (e.g.Steidel et al. 2016;Jaskot & Ravindranath 2016;Senchyna et al. 2017) which imply the presence of young stellar populations producing copious amounts of ionising photons.We find higher EWs of nebular O iii], [Si iii] and C iii] predominantly in the low- sep and -FWHM Ly stacks (Figure 6, panels h1-h2, i1-i2 and j1-j2), supporting the idea that LyC leakers have an increased production rate of ionising 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., panels 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 subsection 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  esc (see subsection 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 ionising 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 (Figure 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 ionising photon production efficiency ( ion ), implying stellar populations with lower metallicities and younger ages as demonstrated by Maseda et al. (2020) for a sample of continuum-faint LAEs at  ≈ 4 − 5.The stacks of LAEs with different Ly luminosities, on the other hand, have nebular emission lines of similar strengths (Figure C1, fourth column), which could be explained by a variety of metallicites and ages presented at fixed Ly luminosity.

He ii: hardness of the ionising spectrum
Whether or not a hard ionising spectrum is a necessary condition for efficient LyC escape remains a subject of debate (Naidu et   flux ratios for the LAE subsamples described in section 3. The observed line ratios are consistent with the photoionisation models of star-forming galaxies (e.g.Feltre et al. 2016), which rules out a significant contribution of AGNs to the nebular emission detected in the stacks.Background histograms are as in Figure 6.
We search for evidence of a harder ionising 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 (Figure 6, panels g1-g3), suggesting an elevated production rate of He + -ionising 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 - asym stacks (Figure 7) which possibly indicates a harder ionising spectrum, given that the ionisation energies for C iii] and O iii] are much lower compared to that of He ii (24.4 eV and 35.1 eV, respectively).Similar results were obtained by Naidu et al. (2022) who studied composite spectra of 26 LAEs at  ∼ 2 from the X-SHOOTER Ly survey (XLS-2).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-ionisation 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).At low redshift ( ∼ 0.3 − 0.4), Schaerer et al. (2022) studied restframe UV spectra of eight LyC emitters and found that the galaxies with  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  ion , and that galaxies with high  esc do not exhibit exceptionally hard ionising spectra compared to other galaxies at similar metallicities.Recently, Marques-Chaves et al. ( 2022) studied a large sample of SFGs from the HST Low- Lyman Continuum Survey and found no correlation between He ii4686/H4861 and  esc , arriving at the same conclusions.Our results suggest that LyC leakers indeed have higher  ion (subsection 5.2), however we cannot rule out the possibility that higher  esc also accompanies harder radiation fields, especially in the light of our C iii]/He ii and O iii]/He ii measurements.

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  uncertainties if we select LyC-leaker candidates using  asym (Figure 6, panels h3, j3 and g3).
We suggest that the specifics of the ISM geometry could break the link between nebular emission and  esc .Unlike  sep and FWHM Ly 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  esc does not require constant support by a strong ionising 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  cov , or the neutral gas covering fraction.Low  cov (or, equivalently, high  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 i1302+Si ii1304 and Si ii1527 are lower for the stacks with more symmetric Ly profiles (Figure 6, panels b3 and d3).The EW of C ii1335 is about −0.9 Å for both stacks (Figure 6, panel c3), but some additional flux bluer to the fitted Gaussian is probably lost in the high- asym stack (see B2). Conversely, the EWs of the LIS absorption lines detected in the  sep and FWHM Ly stacks do not show any clear trend (Figure 6, panels b1-b2, c1-c2, d1-d2), suggesting that  cov is loosely coupled to  sep and FWHM Ly .
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 (Figure C1, panels b1-d1) which was reported by previous studies (e.g.Shapley et al. 2003;Jones et al. 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 idealised 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 in all of the  sep and FWHM Ly stacks, with the highest Si ii * EW in the lowest FWHM Ly quartile (Figure 6, panels e1-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 ionised ISM.Ly (and LyC) photon escape in an anisotropic medium could be attributed to the presence of highly ionised channels of low optical depth (e.g.Zackrisson et al. 2013;Rivera-Thorsen et al. 2017).In this case, the relative transparency of these channels determines  sep and FWHM Ly (McKinney et al. 2019;Jaskot et al. 2019;Kakiichi & Gronke 2021).The distribution of high-column-density gas, in turn, might be more relevant to  cov in the context of the "picket-fence" geometry (Jaskot et al. 2019), which we associate with  asym .In good agreement with this scenario, the Si ii * EW in the low- asym stack is consistent with zero, indicating that a relatively low number of dense molecular clouds imply  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  ∼ 2−4.Most of these galaxies including the Sunburst Arc ( = 2.4, Rivera-Thorsen et al. 2017), Ion2 (𝑧 = 3.2, Vanzella et al. 2015) and Ion3 ( ≈ 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  2022) for a sample of low-redshift ( ∼ 0.3 − 0.4) galaxies with  esc > 10% (C iv/C iii] ≳ 0.75; shaded area).The rapid growth of C iv/C iii] with increasing  esc suggests that this line ratio could serve as an alternative indirect tracer of  esc at redshifts where Ly is significantly attenuated by the IGM ( ≳ 6).Background histograms are as in Figure 6.
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.

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-ionisation 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-ionisation gas.However, the Ly transmission declines rapidly at  > 6 due to the increasing opacity of the neutral IGM (e.g.Gronke et al. 2021), which emphasises the importance of C iv as a standalone probe for  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 ( ∼ 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 (  esc > 0.1): C iv/C iii] > 0.75.At higher redshifts ( = 3.1 − 4.6), Saxena et al. (2022b) found that three out of five C iv emitters with  esc ≈ 0.05 − 0.3 inferred from the LIS absorption lines have similar values of C iv/C iii] (≳ 0.75).
Figure 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- sep , FWHM Ly Q1 and Q2, and low- 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- sep , FWHM Ly Q3 and Q4, and high- asym ).The most striking difference can be seen between the FWHM Ly quartiles, where the C iv to C iii] ratio drops down from C iv/C iii] = 1.0 +0.3 −0.2 (FWHM Ly , Q1) to zero (FWHM Ly , Q4).We argue that the observed C iv/C iii] ratios in the stacks of potential leakers result not only from more intense ionising radiation fields, but also from the increased transparency of the ISM.First, we compare the C iv and He ii emission in the high- sep and -FWHM Ly stacks (Figure 6, panels f1-f2, g1-g2).We find that the C iv EW is consistent with zero but the He ii EW is not, even though the ionisation 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-ionisation 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 (Figure 5, Figure B1 and Figure B2).The non-zero velocity offsets are a clear signature of resonant scattering of C iv by high-ionisation 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.

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 ( sep ), full width at half maximum (FWHM Ly ) and asymmetry parameter ( asym ).We perform spectral stacking and obtain high-S/N detections of rest-frame UV emission and absorption lines containing valuable information on the ionising 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: • The stacks of LyC-leaker candidates generally show strong nebular emission, revealing an extreme ionisation state of the ISM (subsection 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.
• Strong He ii emission (EW ∼ 1 − 3 Å) is typical among LyCleaker candidates, implying a significant production rate of ionising photons with energies > 54.4 eV (subsection 5.3).While elevated  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 ionising spectrum possibly associated with higher  esc .
• The LIS absorption lines are generally weaker for the low- asym stacks, suggesting that the asymmetry of the Ly line profile depends on the distribution of high-column-density gas and hence traces  cov in the framework of the "picket-fence" model of the ISM.Conversely, the  sep and FWHM Ly 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 (subsection 5.4).
• Si ii * emission is detected in all the  sep and FWHM Ly stacks, implying the presence of neutral gas off the LOS even in potential LyC leakers (subsection 5.4).This indicates that LyC photons escape from a highly anisotropic ISM through ionised channels surrounded by a high-column-density medium.
• High C iv/C iii] ratios (> 0.75) are common among the LAEs with potentially high  esc (≳ 0.1; subsection 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 ionising 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 ionising photon production efficiency and the opacity of the ISM (Figure 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  esc at high redshift ( ≳ 3) will likely become possible in the future thanks to blue-optimised integral-field spectrographs such as BlueMUSE (Richard et al. 2019), which will probe LyC emission from sources at  ∼ 3−4.Subsequently, C iv/C iii] and other spectral properties of SFGs indicative of high  esc could be measured at  ≳ 6 with JWST and the new generation spectrographs such as MOONS/VLT, PFS/Subaru or HARMONI/ELT, opening new paths for indirect identification of ionising sources in the EoR.

Figure 2 .
Figure 2. Venn diagram showing the samples of LyC-leaker candidates selected as having low Ly 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 profiles.Bottom: Total sample.The estimates of the median  esc are based on the fit fromIzotov et al. (2018b).All three properties of the Ly line act as tracers of  esc due to their high sensitivity to the neutral hydrogen column density -a key parameter of the ISM controlling escape of ionising photons.

Figure 3 .
Figure3.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 star-forming galaxies.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.
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 nonresonant 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 (Eq.(2) of V18; dark blue).Top to Bottom: Line FWHMs, S/N ratios and EWs.Vertical error bars mark the 16-th and the 84-th 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.

Figure 5 .
Figure5.Median-stacked spectra of LAEs with low (< 463 km s −1 ; blue) and high (> 463 km s −1 ; 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).

Figure 6 .
Figure6.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; panels a1-a3), ISM absorption lines (sea foam green; panels b1-d3) and the Si ii * 1533 fine-structure transition (moss green; panels e1-e3).Right column: C iv (lavender; panels f1-f3), He ii (yellow; panels g1-g3) and collisionally excited nebular emission lines (pink; panels h1-j3).The color-coding of markers follows that of the stacked spectra of respective subsamples in Figure5, FigureB1and FigureB2.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.

Figure 7 .
Figure 7. [C iii]1907+Ciii]1909/He ii1640 and O iii]1666/He ii1640 flux ratios for the LAE subsamples described in section 3. The observed line ratios are consistent with the photoionisation models of star-forming galaxies (e.g.Feltre et al. 2016), which rules out a significant contribution of AGNs to the nebular emission detected in the stacks.Background histograms are as in Figure6.

Figure 8 .
Figure8.C iv1548,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 bySchaerer et al. (2022) for a sample of low-redshift ( ∼ 0.3 − 0.4) galaxies with  esc > 10% (C iv/C iii] ≳ 0.75; shaded area).The rapid growth of C iv/C iii] with increasing  esc suggests that this line ratio could serve as an alternative indirect tracer of  esc at redshifts where Ly is significantly attenuated by the IGM ( ≳ 6).Background histograms are as in Figure6.

Figure B1 .
Figure B1.Similar to Figure 5, but for the subsamples of LAEs with different Ly FWHMs.

From panel a-f: redshift
(Childs & Stanway 2018)cape through a porous ISM with little scattering (e.g.Verhamme et al. 2015).However, quantifying the asymmetry of Ly profiles is challenging at relatively low spectral resolutions ( ≲ 3000) and high skewness(Childs & Stanway 2018).After applying a moderate cutoff on the  asym uncertainty (  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  asym < 0.19 (potential leakers) and  asym > 0.19 (potential non-leakers), both comprising of 277 objects.The quantitative comparison between  asym and  esc is missing from the literature.If we use FWHM Ly to estimate the median  esc , we obtain values of 2.0 per cent (low- asym subsample) and 2.8 per cent (high- asym subsample) that do not comply with our classification -LAEs with more symmetric Ly profiles are expected to have higher  esc .Furthermore, we find that only 81 objects (29 per cent) from the low- asym subsample are part of the low- sep and -FWHM Ly subsamples (Figure2, bottom).This could indicate that  asym , on the one hand, and  sep and FWHM Ly , on the other hand, are sensitive to different geometries of the ISM (namely, Erb et al. 2014).More symmetric profiles might suggest the presence of an unresolved blue bump blended with the red peak of Ly (see subsection 2.4), or Ly emission at systemic velocity in cases where