The VANDELS survey: the ionizing properties of star-forming galaxies at $3 \leq z \leq 5$ using deep rest-frame ultraviolet spectroscopy

To better understand the ionizing properties of galaxies in the EoR, we investigate deep, rest-frame ultraviolet (UV) spectra of $\simeq 500$ star-forming galaxies at $3 \leq z \leq 5$ selected from the public ESO-VANDELS spectroscopic survey. The absolute ionizing photon escape fraction ($f_{\rm esc}^{\rm abs}$) is derived by combining absorption line measurements with estimates of the UV attenuation. The ionizing production efficiency ($\xi_{ion}$) is calculated by fitting the far-UV (FUV) stellar continuum of the VANDELS galaxies. We find that the $f_{\rm esc}^{\rm abs}$ and $\xi_{ion}$ parameters increase towards low-mass, blue UV-continuum slopes and strong Ly$\alpha$ emitting galaxies, and both are just slightly higher-than-average for the UV-faintest galaxies in the sample. Potential Lyman Continuum Emitters (LCEs) and selected Lyman Alpha Emitters (LAEs) show systematically higher $\xi_{ion}$ ($\log \xi_{ion}$ (Hz\erg) $\approx 25.38, 25.41$) than non-LCEs and non-LAEs ($\log \xi_{ion}$ (Hz\erg) $\approx 25.18, 25.14$) at similar UV magnitudes. This indicates very young underlying stellar populations ($\approx 10~{\rm Myr}$) at relatively low metallicities ($\approx 0.2~{\rm Z_{\odot}}$). The FUV non-ionizing spectra of potential LCEs is characterized by very blue UV slopes ($\leq -2$), enhanced Ly$\alpha$ emission ($\leq -25$A), strong UV nebular lines (e.g., high CIV1550/CIII]1908 $\geq 0.75$ ratios), and weak absorption lines ($\leq 1$A). The latter suggests the existence of low gas-column-density channels in the interstellar medium which enables the escape of ionizing photons. By comparing our VANDELS results against other surveys in the literature, our findings imply that the ionizing budget in the EoR was likely dominated by UV-faint, low-mass and dustless galaxies.


INTRODUCTION
Several data sets (see Goto et al. 2021, and references therein) measuring the redshift ( ) evolution of the volume-averaged neutral ★ E-mail: alberto.saldanalopez@unige.ch hydrogen fraction provide evidence for the last major phase change underwent by the Universe, the Cosmic Reionization. Between = 9 and = 6 (Planck Collaboration et al. 2016), the number of ionizing photons emitted per unit time ( ion ) overcame the recombination rate (Γ HI ) of hydrogen atoms, ion ≥ Γ HI (Madau et al. 1999), so that the bulk of H neutral gas within the Intergalactic Medium (IGM) progressively transitioned to an ionized state (see Dayal & Ferrara 2018, for a review).
Yet, the sources mainly responsible for expelling such a vast amount of ionizing (or Lyman Continuum, LyC) photons to the IGM remain elusive (see Robertson 2022, for a review). Overall consensus exists about the minor contribution of Active Galactic Nuclei (AGN) to the ionizing budget during the Epoch of Reionization (EoR, Hassan et al. 2018;Kulkarni et al. 2019;Dayal et al. 2020), mainly because of their lower number density at early epochs (e.g., Matsuoka et al. 2018, but see Fontanot et al. (2012) and Cristiani et al. (2016) for a different viewpoint). This said, AGNs have played a key role in keeping the IGM ionized after the EoR (see e.g., Becker & Bolton 2013), whereas stars seem to provide a small contribution to the ionizing radiation budget at < 5 (Tanvir et al. 2019). Clearly, assuming that star-forming (SF) galaxies drove Cosmic Reionization shifts the current debate on whether more massive, UV-bright galaxies (Madau & Haardt 2015) or in opposite, low-mass, UV-faint counterparts (Robertson et al. 2013(Robertson et al. , 2015 dominated the ionizing emissivity at the EoR. On the one hand, the remarkable modeling efforts by Sharma et al. (2016) and Naidu et al. (2020), based on the evolution of the star-formation rate surface-density of galaxies, and the works by Naidu et al. (2022) and Matthee et al. (2022), based on the fraction of Ly emitters (LAEs) over time, support a late-completed and rapid Reionization purely dominated by more massive and moderately luminous sub− ★ UV systems. These results are compatible with the rapid reionization modeled by Mason et al. (2019), in which the ionizing emissivity was constrained from CMB optical depth and Ly forest dark pixel fraction data. On the other hand, semi-empirical models like those in Finkelstein et al. (2019), based on observational constraints on the UV luminosity function during the EoR, and independent cosmological hydrodynamical simulation such as Rosdahl et al. (2022), suggest an early-completed Reionization conducted primarily by low-mass and fainter galaxies (see also Trebitsch et al. 2022). The most recent measurements of the mean free path of ionizing photons by Becker et al. (2021) have shown a much shorter value than previously thought at 5 ≤ ≤ 6, supporting the rapid reionization scenario that, according to recent models, could still be conducted by the faintest galaxies (Cain et al. 2021).
In observations, the problem reduces to solving the equation for the ionizing emissivity of the average galaxy population ( ion ), i.e., the number of ionizing photons emitted per unit time and comoving volume (Robertson 2022): where abs esc stands for the absolute escape fraction of ionizing photons (i.e., the ratio between the number of escaping versus produced LyC photons by massive stars), and ion is the so-called ionizing photon production efficiency (ionizing photons generated per non-ionizing intrinsic UV luminosity). UV accounts for the non-ionizing UV luminosity density at the EoR, resulting from the integral of the UV luminosity function (UVLF, the number of galaxies per UV luminosity and comoving volume).
The UVLF (or UV , equivalently) of galaxies is relatively wellconstrained up to the very high-redshift Universe (Bouwens et al. 2015;Davidzon et al. 2017;Bouwens et al. 2021;Donnan et al. 2023;Finkelstein et al. 2023). At the EoR, some works observe a decrease in the UV luminosity density from = 8 to 10 Harikane et al. 2022), compatible with a fast build-up of the dark matter halo mass function at those redshifts while, in contrast, some others do not find such a suppression at all (McLeod et al. 2016;Livermore et al. 2017). So far, UV-bright galaxies are found to be several orders of magnitude less numerous than the bulk of the UV-faint detected galaxies (Atek et al. 2015), and therefore thought to play a minor role during Reionization (but see Marques-Chaves et al. 2022a), although a possible excess in the number of sources at the bright-end of the UVLF (Rojas-Ruiz et al. 2020), might make flip the argument towards an EoR governed by the "oligarchs".
The ionizing production efficiency of galaxies ( ion ) is more uncertain. Surveys targeting H emitters (HAEs, see Bouwens et al. 2016;Matthee et al. 2017a;Shivaei et al. 2018;Atek et al. 2022) and LAEs (Harikane et al. 2018;Nakajima et al. 2018a) at intermediate redshifts show, in general, higher production efficiencies for the low-mass and UV-faint galaxies (Prieto-Lyon et al. 2022), although with a huge scatter within which ion spans a wide range of values depending on the galaxy type (log ion (Hz/erg) = 24 − 26). Reassuringly, the overall ion evolution with galaxy properties at these redshifts is in line with the usual formalism by which a log ion (Hz/erg) ≥ 25.2 SFG population with constant SFH (over 100Myr) is able to fully reionize the Universe, assuming a fixed escape fraction of abs esc ≥ 5% (Robertson et al. 2013). Even so, given the stochastic nature of the LyC emission, linked to the predominantly bursty star-formation histories (SFHs) in low-mass SFGs (Muratov et al. 2015), these assumptions might not realistically apply anymore (see discussion in Atek et al. 2022).
Since the flux at LyC wavelengths will not be accessible at the EoR due to the increase of the IGM opacity (Inoue et al. 2014), indirect tracers of LyC radiation are needed. Thanks to the advent of the LzLCS (Flury et al. 2022a), in which both far-UV (FUV) and optical spectra are available for each galaxy, several abs esc diagnostics have been statistically tested for the first time in Flury et al. (2022b). Among them, those which properly account for the neutral gas and dust column densities as well as for geometrical effects (see Seive et al. 2022), so that LyC radiation escapes only along favoured, cleared sight-lines in the interstellar medium (ISM), remain the most promising proxies. In particular, the peak separation of the Ly line (Verhamme et al. 2017;Gazagnes et al. 2020;Naidu et al. 2022), the depth of the low-ionization state (LIS) UV absorption lines (Reddy et al. 2016b;Chisholm et al. 2018;Gazagnes et al. 2018;Saldana-Lopez et al. 2022) and the Mg doublet ratio (Henry et al. 2018;Chisholm et al. 2020;Xu et al. 2022) seem to closely probe the measured escape fraction of LzLCS galaxies (but see cautionary theoretical work by Mauerhofer et al. 2021;Katz et al. 2022).
In this work, we aim to indirectly study the ionizing properties of high-SFGs and their evolution along with the different galaxy-properties. For that, we make use of a sample of 500 deep, rest-frame ultraviolet spectra at 3 ≤ ≤ 5 drawn from the VAN-DELS survey (McLure et al. 2018b;Pentericci et al. 2018;Garilli et al. 2021). In particular, LyC absolute escape fractions ( abs esc ) are derived by measuring the depth of the absorption lines in combination with the UV attenuation, whilst ionizing production efficiencies ( ion ) are computed based on the best Spectral Energy Distribution (SED) fit to the FUV stellar continuum (based on the work by Chisholm et al. 2019). Our study complements ongoing efforts to understand the properties of LyC emitting galaxies at low (LzLCS, Flury et al. 2022a) and high redshifts (KLCS, Steidel et al. 2018), thereby setting the pathway to interpret the ionizing signatures of EoR-galaxies, whose number of detections has been dramatically boosted thanks to the high-quality performance of the first JWST observations.
The layout of this article is as follows. The VANDELS survey and sample definition are described in Sect. 2. The code for fitting the stellar SED of VANDELS spectra, and the methods for predicting individual ionizing efficiencies and escape fractions are outlined in Sect 3. The main results of this research, looking for correlations between abs esc and ion with the different galaxy properties are summarized in Sect. 4. Our results on the ionizing efficiency of high-galaxies are compared with different estimates in the lit- erature in Sect. 5, and finally the possible redshift evolution of the abs esc × ion product is discussed in Sect. 6, by comparing our values against state-of-the-art low-and high-surveys. We summarize our findings in Sect. 7.
Throughout this paper, a standard flat ΛCDM cosmology is used, with a matter density parameter Ω M = 0.3, a vacuum energy density parameter Ω Λ = 0.7, and a Hubble constant of 0 = 70 km s −1 Mpc −1 . All magnitudes are in AB system (Oke & Gunn 1983), and we adopt a solar metallicity value of 12 + log(O/H) = 8.69. All the stellar metallicities are quoted relative to the solar abundance (Z ) from Asplund et al. (2009), which has a composition by mass of Z = 0.014. Emission and absorption line equivalent widths (EWs) are given in rest-frame (unless stated otherwise), with positive (negative) EWs meaning lines seen in absorption (emission).

SAMPLE
We rely on rest-frame FUV spectroscopy in order to estimate the ionizing properties of high-SFGs. At 3 ≤ ≤ 5, the rest-frame UV spectrum of galaxies (1200 − 2000Å) is accessible from the ground through optical spectroscopy. The criteria by which we build our sample of high-rest-UV spectra (2.1), and the estimation of the main SED parameters (2.2) and other spectroscopic features (2.3) are described in detail in the following sections. Additionally, the survey properties of the two main comparison samples of LCEs in the literature are summarized (2.4).

Deep rest-frame UV spectra: the VANDELS survey
The VANDELS survey (final Data Release 4 in Garilli et al. 2021) is an ESO Public Spectroscopic Survey composed of around 2100 optical, high signal-to-noise ratio (S/N) and medium-resolution (R) spectra of galaxies at redshifts = 1−6.5. The VANDELS footprints are centered on two of the HST Cosmic Assembly Near-Infrared Deep Extragalactic Legacy Survey (CANDELS, Grogin et al. 2011;Koekemoer et al. 2011) fields, in particular the CDFS (Chandra Deep Field South, see Guo et al. 2013) and the UDS (UKIDSS Ultra Deep Survey, see Galametz et al. 2013), but covering a wider area. The primary targets were selected from the parent CDFS and UDS photometric catalogs attending to the quality of their photometric redshifts (see McLure et al. 2018b;Pentericci et al. 2018, for details).
Ultra deep, multi-object optical spectroscopy (at 4800 − 10000Å observed-frame) was conducted with the VIMOS instrument (Le Fèvre et al. 2003), at the Very Large Telescope (VLT), resulting in 1019 and 1068 sources observed in the CDFS and UDS, respectively, over 2087 measured spectroscopic reshifts in total. The final VANDELS catalog reaches a target selection completeness of 40% at AB = 25. VANDELS spectra have an unprecedented average S/N ≥ 7 per resolution element over 80% of the spectra, thanks to exposure times spanning from 20 up to 80 hours on source. The VIMOS -Full-With at Half Maximum-spectral resolution is ≡ /Δ ≈ 600 at 5500Å, with a wavelength dispersion of 2.5Å/pixel. For additional information about the survey design and data reduction, we encourage the reader to check McLure et al. (2018b); Pentericci et al. (2018) and Garilli et al. (2021) papers.
Our final sample comprises VANDELS-DR4 galaxies under the SFGs and LBGs categories only, whose spectroscopic redshifts fall in the 3 ≤ spec ≤ 5 range and have the highest reliability in the redshift measurement (i.e., flag=3/4). Given the VIMOS instrumental sensitivity, the redshift constraint ensures complete wavelength coverage from 1200 to 1600Å for all galaxies in the sample (in practice, from Ly 1216 to C 1550). We then compute the mean S/N per pixel over two regions free of absorption features, specifically 1350 − 1375Å and 1450 − 1475Å in the restframe (25Å). Then, we select the galaxies with S/N ≥ 2 in the two spectral windows simultaneously ( 50% of the remaining sample). Worth noticing is that the uncertainties of the DR4 spectra were systematically underestimated (see Garilli et al. 2021). For this reason, we applied individual correction factors to the error spectra provided by the VANDELS collaboration (Talia et al., in preparation). The average correction factor is ×1.4.
In summary, our VANDELS working sample includes 534 galaxies at 3 ≤ spec ≤ 5 with a median S/N ≈ 5 in the 1400Å continuum range, where 297 sources were observed in the CDFS versus 237 in the UDS field. Fig. 1 shows the observed −band (either HST/F160W or VISTA/H) apparent magnitude distribution as a function of the spectroscopic redshift for the selected sources, together with the entire VANDELS-DR4 sample. The median (16 ℎ and 84 ℎ percentiles) F160W and spec of the working sample cor- Our working sample is indicated through red filled symbols, and the large blue circles display the SFR running medians at each inter-quartile range of stellar mass. The golden dashed line follows the MS relation by Speagle et al. (2014), for comparison. Our sample fall ∼ 0.1dex. systematically above the MS at = 3.5, but it still probes the upper bound of the SFG population at these redshifts. respond to AB = 24.8 +0.8 −0.9 AB and spec = 3.56 +0.41 −0.36 , respectively.

Photometric properties of the VANDELS galaxies: stellar masses
The wide imaging coverage of VANDELS allow a robust characterization of the SED of every galaxy (Garilli et al. 2021). The physical integrated properties of the VANDELS-DR4 galaxies were obtained by SED fitting using the Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (B ) code 1 . B (see Carnall et al. 2018) inputs the spectroscopic redshifts measured by the VANDELS team (Pentericci et al. 2018) and all the available CANDELS (Galametz et al. 2013;Guo et al. 2013) plus ground-based photometry (McLure et al. 2018b) to run a Bayesian algorithm which samples the posterior of the SED parameters. It makes use of the 2016 updated version of the Bruzual & Charlot (2003) models, including the MILES stellar spectral library (Falcón-Barroso et al. 2011) and the stellar evolutionary tracks of Bressan et al. (2012) and Marigo et al. (2013). An exponentially declining SFH is assumed, with a minimum timescale of 10 Myr and a minimum age of 10 Myr. The stellar metallicity was fixed to 0.2 , coinciding with the average stellar metallicity found by Cullen et al. (2019) and Calabrò et al. (2021), probing a similar sample of galaxies in VANDELS. The dust attenuation was modeled using the Salim et al. (2018) prescription, and nebular emission was considered by adopting a constant ionization parameter of log( ) = −3 in the fits.
For our VANDELS-DR4 working sample, the median stellar mass and SFR are log( ★ /M ) = 9.34 +0.31 −0.39 and 1 B (Carnall et al. 2018) is a state-of-the-art P code for modeling galaxy spectra and fitting spectroscopic and photometric data, for details go to https://bagpipes.readthedocs.io/en/latest/. SFR/(M yr −1 ) = 14.46 +25.30 −8.47 , respectively, with 0.1M and 2.4M yr −1 typical uncertainties. The log SFR − log M ★ distribution for the selected sample lays 0.1 to 0.2 dex above the SF mainsequence (MS) relation of Speagle et al. (2014) at similar redshifts, with a scatter of ±1 dex in SFR, approximately (see McLure et al. 2018b). In Fig. 2, we explicitly show how the SFR running medians at the 25 ℎ , 50 ℎ and 75 ℎ stellar mass inter-percentile ranges fall above the Speagle et al. (2014) MS relation at = 3.5. This means that, at comparable stellar masses (i.e., log( ★ /M ) = 8.5 − 10.5), our sample is probing a slightly higher SFR regime than the bulk of SFGs, although it is still representative of the whole population given the intrinsic scatter of the MS at high-(see the discussion in Cullen et al. 2021). A similar behavior was found in Calabrò et al. (2022). Finally, our sample also suffers from luminosity (or stellar mass) incompleteness i.e., the highest redshifts objects ( 4) are biased towards the brightest AB absolute magnitudes because of the current flux-limiting nature of the VANDELS survey.

Spectral properties of the VANDELS galaxies: Ly equivalent widths
In this work, we make use of the methods and measurements provided by the VANDELS collaboration (Talia et al., in preparation) to compute the Ly equivalent width ( Ly ) of galaxies, either in emission or absorption. The Ly (H 1216) flux and equivalent width were measured by fitting a single-Gaussian profile to the Ly line using the dedicated code 2 . If the S/N of the line reach an input threshold, the script allows a ±1000 kms −1 offset of the line center respect to rest-frame value given by the spectroscopic redshift. The continuum level was defined by the Bruzual & Charlot (2003) best-fit template to the entire spectrum (which gives a first order correction for underlying stellar absorption), and the flux and equivalent width was measured over ±8000 kms −1 each side of the line peak (see Sect. 3.2 for the actual definition of the equivalent width). By convention, the fits report negative equivalent widths when the Ly line appears in emission, whereas a positive value is reported when the line is absorption dominated. Given the limited variety of Ly profiles observed in the VANDELS sample due to a medium ≈ 600 spectral resolution (see Kornei et al. 2010;Cullen et al. 2020, and methods therein), this approach only constitutes a first order estimation of Ly . We also note that, according to Talia et al. (in preparation), for 198 object ( 37% of the sample) a single Gaussian fit did not provide a good fit to the Ly line. Rather, two Gaussian components were needed. However, a more detailed calculation is out of the scope of this paper, and we finally adopt single Gaussian fits as our fiducial estimates of the Ly flux. Doing so, the median Ly equivalent width is Ly (Å) = 4.58 +19.01 −9.43 for our selected sample. According to the definition by Pentericci et al. (2009), 104 out of 534 of galaxies ( 20% of the sample) can be classified as LAEs with Ly ≥ −20Å (although other authors may adopt different criteria, e.g., see Stark et al. 2011;Nakajima et al. 2018a;Kusakabe et al. 2020).
2 is a C++ simple software that can be used to derive spectroscopic redshifts from 1D spectra and measure line properties (fluxes, velocity widths, offsets, etc.). More in https://github.com/cschreib/ slinefit, by Corentin Schreiber.

Comparison samples of LCEs
Asides from our VANDELS sample, primarily composed by SFGs at 3 ≤ ≤ 5, we will consider as benchmarks other samples of confirmed LCEs from the literature (see Sect. 6). The LCEs comparison samples are the following: The Low-Lyman Continuum Survey, or LzLCS, is a large HST programme (Flury et al. 2022a, PI: Jaskot, HST Project ID: 15626) targeting 66 SFGs at 0.22 ≤ ≤ 0.43, selected from SDSS and GALEX observations. Each galaxy was observed with the low-resolution COS/G140L grating, covering the LyC and the non-ionizing rest-frame FUV continuum (600 − 1450Å), with a spectral resolution of ≈ 1000 at 1100Å. In combination with these 66 galaxies, the LzLCS also includes a compilation of 23 archival sources drawn from the literature (Izotov et al. 2016b(Izotov et al. ,a, 2018aWang et al. 2019;Izotov et al. 2021). Out of 89 galaxies, 50 were detected in the LyC, with a median absolute escape fraction of abs esc = 0.04 +0.16 −0.02 , while the remaining 39 galaxies show strong abs esc upper limits, typically below 1%. The Keck Lyman Continuum Survey, or KLCS, is an extensive Keck/LRIS spectroscopic campaign (Steidel et al. 2018) to obtain deep rest-UV spectra of selected LBGs at 2.8 ≤ ≤ 3.4. The LRIS spectra of KLCS galaxies cover the LyC region and a wide FUV window longwards of the Ly line (880 − 1650Å), at a spectral resolution of ≈ 800 for LRIS-B and ≈ 1400 for LRIS-R. In the most recent HST imaging follow-up by Pahl et al. (2021Pahl et al. ( , 2023, they re-analyze a sample of 124 KLCS galaxies and look for foreground and low-redshift contamination that could potentially lie in projection within the angular extent of each LyC detected galaxy. In total, 13 individually-detected and 107 LyC undetected sources resulted in a sample-averaged escape fraction of abs esc = 0.06 ± 0.01.

METHODS
In this section, we first describe the F CUS code (Sect. 3.1), a customized P script to fit the stellar continuum of extragalactic UV spectra. Then, we present the methodology to indirectly infer the ionizing absolute escape fraction of galaxies ( abs esc ), by using the information provided by the UV absorption lines and the global SED fit to the non-ionizing UV continuum (Sect. 3.2). Finally, we discuss about the procedure for stacking individual VANDELS spectra, and the systematic effects included in the measurements due to the instrumental resolution and stacking (Sect. 3.3).

The F CUS code: F tting the stellar Continuum of U Spectra
With the purpose of quantifying the stellar continuum properties underlying the VANDELS spectra, we have created F CUS 3 . F CUS is a customized P script that stands for FItting the stellar Continuum of Uv Spectra and, in short, it returns an estimation of the galaxy light-weighted stellar age, metallicity and dust extinction as well as other secondary SED parameters by using a combination of best-fit stellar population templates. This methodology was widely described and previously tested in Chisholm et al. (2019, hereafter C19), and has been used in other papers such as Gazagnes et al. (2018Gazagnes et al. ( , 2020 and Saldana-Lopez et al. (2022, hereafter SL22). We refer the reader to those papers for similar approaches but slightly different applications of the current method.

Stellar population models
As described in SL22, the stellar continuum modeling was achieved by fitting every observed spectrum with a linear combination of multiple bursts of single-age and single-metallicity stellar population models. By default, F CUS inputs the fully theoretical -99 single-star models without stellar rotation (S99, Leitherer et al. 2011Leitherer et al. , 2014 using the Geneva evolution models (Meynet et al. 1994), and computed with the WM-B method (Pauldrach et al. 2001;Leitherer et al. 2010). The S99 models assume a Kroupa (2001) initial mass function (IMF) with a high-(low-)mass exponent of 2.3 (1.3), and a high-mass cutoff at 100 M . The spectral resolution of the S99 models is ≡ /Δ 2500, and it remains approximately constant at FUV wavelengths.
Four different metallicities (0.05, 0.2, 0.4 and 1 ) and ten ages for each metallicity (1,2,3,4,5,8,10,15,20 and 40 Myr) were chosen as a representative set of ×40 models for our high-UV spectra. In detail, ages were chosen to densely sample the stellar ages where the stellar continuum features appreciable change (see C19). For example, at older ages than 40 Myr, the FUV stellar continuum does not appreciably change such that a 40 Myr model looks very similar to a 100 Myr population. We caution the reader that the S99 high-resolution models do not densely sample the HR diagram at effective temperatures below 15,000K (Leitherer et al. 2010). The cooler models have artificially bluer continuum slopes, which makes older stellar populations appear slightly bluer than expected for their temperatures (for details, see Chisholm et al. 2022).
Finally, a nebular continuum was generated by self-consistently processing the stellar population synthesis models through the 17.0 code 4 (Ferland et al. 2017), assuming similar gasphase and stellar metallicities, an ionization parameter of log( ) = −2.5, and a volume hydrogen density of n = 100 cm −3 . The output nebular continua for each stellar population were added to the stellar models. The inclusion of the nebular continuum is only appreciable at the youngest ages of the bursts (≤ 10Myr) at wavelengths ≥ 1200Å, and produces redder stellar models than before, which has a pronounced effect on the fitted B−V of the young stellar populations. C19 tested the effect of different ionization parameters on the fitted stellar ages and metallicities, and found that the goodness of the fits does not significantly change for different log( ) values.

Fitting method
The observed spectra were first manually placed into the rest-frame by multiplying by the corresponding 1/(1 + ) factor. Both the spectra and the models were then normalized by the median flux within a wavelength interval free of stellar and ISM features (1350− 1370Å), and all the fits were performed in the same rest-wavelength range (1200 − 1925Å in the case of VANDELS spectra). Finally, the models were convolved by a Gaussian kernel to the instrumental resolution ( (VIMOS) ≈ 600).
First, the non-stellar features and the spectral regions that are affected by host-galaxy ISM absorptions and sky subtraction residuals had to be masked out. Bad pixels with zero flux were neither considered in the fit. When ran with F CUS, the UV stellar continuum ( ★ ( )) is fitted with a linear combination of multiple S99 models plus the nebular continuum. 99 single stellar population is usually not able to fully reproduce both the stellar features and the shape of the continuum of extragalactic spectra, such that a combination of instantaneous burst models is actually needed (C19).
Adopting a simple geometry where the dust is located as a uniform foreground screen surrounding the galaxy 5 , this results in: 2,3,4,5,8,10,15,20,40 ≡ 0.05, 0.2, 0.4, 1 where 10 −0.4 is the UV attenuation term, = B−V , and is given by the adopted dust-attenuation law. The linear coefficients determine the weight of single-stellar population within the fit, and the best fit is chosen through a non-linear 2 minimization algorithm with respect to the observed data ( 6 package, Newville et al. 2016). Errors were derived in a Monte-Carlo way, varying observed pixel fluxes by a Gaussian distribution whose mean is zero and standard deviation is the 1 −error of the flux in the same pixel, and re-fitting the continuum over iterations per spectrum (we chose = 100, enough realizations to sample the posterior continuum so that it approaches "Gaussianity" on each pixel). Fig. 3 shows the F CUS output for one of the galaxies in VAN-DELS (objID:CDFS017345, at = 3.61), with the observed (in black) and fitted stellar continuum (in red). The bottom panel of this figure additionally shows the distribution of the light-fractions at a given burst-age and metallicity for the CDFS017345 best-fit. This galaxy is dominated by a young and low-metallicity population ( 60% of the total light at 3Myr) with 4.80 Myr average light-weighted age and 0.18 Z light-weighted metallicity (see the following section). CDFS017345 is moderately attenuated, with a UV dust-attenuation parameter of B−V 0.12 mag. In Fig. 3, and together with other nebular (gold) and ISM lines (black) usually present in VANDELS spectra, we also highlight (in blue) the main stellar features that helped our algorithm to estimate the age and the metallicity of the stellar population: the N 1240 and C 1550 P-Cygni stellar wind profiles (see also O 1371 and Si 1402). Note that the inclusion of C 1550 requires a careful masking of the strong ISM component of this line, where only the blue wing of the P-Cygni was used in the fit. While N 1240 is mostly sensitive to the age of the burst -with the P-Cygni profile being more prominent at younger ages-, the C 1550 feature, on the other hand, changes accordingly to the metallicity (C19) of the stellar population -whose asymmetric profile is stronger at higher stellar metallicities-. Therefore, these two distinct spectral features partially break the age-★ degeneracy in the FUV, although the method is still affected by age-attenuation degeneracy effects. 5 The effect of a clumpy gas-to-dust geometry has been discussed in Gazagnes et al. (2018) and SL22. Regrettably, the current resolution and S/N of the VANDELS spectra is not enough to disentangle between the two geometries, so we assume the most simple scenario of a uniform screen of dust as a our fiducial case. 6    summary chart with the best-fit SED fitting parameters, i.e., reduced-2 , B−V (in mag., R16), light-weighted stellar age (in Myr) and metallicity (in Z ).

Results
Given the best-fit SED for the stellar continuum and the best-fit coefficients in Eq. 2, F CUS also provides a handful of other secondary SED parameters, that we described here.

Stellar age and metallicity
The light-weighted stellar age and metallicity ( ★ ) of the best-fit stellar population can be easily obtained using the weights as: where age , are the age and metallicity of -th synthetic stellar model. For instance, in the example of Fig. 3 Fig. 4). Compare to these works, our average metallicity is slightly higher: log( ★ /Z ) −0.63, although still compatible within the spread of the distributions.
Our shift in log ★ with respect to the previous works must be driven by use of instantaneous burst of star-formation instead of a continuous SFH in the theoretical stellar models. Our election of single-burst models is motivated by the fact that mixed-age populations usually do a better job when reproducing age-sensitive tracers (such us N 1240) of individual FUV galaxy-spectra (see Sect. 5.3 in C19). Moreover, by accounting for the dominance of young stellar populations, bursty star formation histories, which are expected to hold for intermediate-to-low-mass SFGs at high-redshift (e.g., Trebitsch et al. 2017), produce significantly more ionizing photons, at higher energies, than continuous star formation histories.
The effect that the choice of either constant SFRs or instantaneous bursts of star-formation have in the stellar age, metallicity and ionizing properties of the mean SFG galaxy population at 3 will be discussed in a future publication. We refer the reader to Cullen et al. (2019) and Calabrò et al. (2021) for more details on the stellar metallicity of VANDELS galaxies.
Dust-attenuation parameter: B−V The amount of dustattenuation is given by the 10 −0.4 term in Eq. 2, where = B−V is defined as the "UV attenuation", and the specific functional form for is determined by the dust-attenuation law. Therefore, the resulting best-fit values for B−V will differ depending on the assumed dust-law, and this will have important consequences for other related quantities which explicitly depend on B−V , like the absolute escape fraction ( abs esc ) as we will discuss later on (see Sect. 6).
With the goal of testing the influence of the dust-attenuation law in our results, we use two extreme cases for in SF galaxies (Shivaei et al. 2020): the Reddy et al. (2016a, hereafter R16) and Small Magellanic Cloud (Prevot et al. 1984;Bouchet et al. 1985;Gordon et al. 2003, hereafter SMC) attenuation curves. As widely discussed in the dedicated sections in SL22, and references therein, the use of R16 is, on the one hand, motivated by the fact that this law is one of the only properly defined below 1000Å by a significant number of galaxies. On the other hand, the use of steeper SMC-like curves has appeared to be more suitable for high-SFGs (Salim et al. 2018), low-metallicity starburst (Shivaei et al. 2020), and LCEs (Izotov et al. 2016a).
After performing the SED fits with F CUS using both the R16 and SMC prescriptions for , we find an average attenuation parameter of B−V = 0.22 +0.09 −0.08 mag. for R16 and B−V = 0.10 +0.04 −0.03 mag. when using SMC (see Fig. 4). This corresponds to UV attenuations of 1.82 and 2.83 mag. at 1600Å and 912Å, using R16 (1.25, 2.63 mag. for SMC). Although the resulting light-weighted stellar ages and metallicities are similar in both sets of fits (similar light-fractions), the UV attenuation term 10 −0.4 is slightly higher for SMC at all wavelengths, meaning that the fitted coefficients are slightly lower for SMC than R16 law, so that both fits can match the observations similarly (reduced-2 distributions are similar). Lower coefficients for SMC means that all absolute quantities derived from the intrinsic SEDs will also change accordingly.

UV-continuum slope (at 1500Å): 1500
spec The UV spectroscopic continuum slope at 1500Å, also called −slope, was computed by fitting a power-law of the form ∝ to every individual spectrum (Meurer et al. 1999). To do so, we take the average flux density in seven 15Å-width spectral windows, between 1275 − 1825Å (similar to the ones in Calzetti et al. 1994), and fit a linear relation to the log − log values using C -F 7 . For consistency, and in order to avoid inhomogeneous wavelength sampling for every source due to the different redshifts, corresponds to the continuum flux obtained from the best-fit SED model of the galaxy ( ★ ( )), which was described in the previous Sect. 3.1. This assumption imprints an intrinsic dependency of the −slope on the chosen attenuation law, since it modifies the shape of the dust-attenuated spectrum through the term. Thus, median of the −slope distribution for the VANDELS sample is 1500 spec = −1.34 +0.45 −0.42 for R16 and 1500 spec = −1.06 +0.54 −0.51 for SMC. When compared to the slopes obtained in Calabrò et al. (2021) for the same objects ( 1500 −1.80), our spectroscopic −slopes are overall redder (lower) than the former by 0.4, approximately (see distribution in Fig. 4). In Calabrò et al. (2021), the UV slopes were derived by fitting a power-law to the available multi-band VANDELS photometry, taking the photometric bands whose bandwidths lie inside the 1230-2750Å rest-frame wavelength range. We investigated the possible reasons that may reach to this discrepancy and conclude that it is mainly caused by the use of a wider range in wavelength with respect to the spectroscopic slopes. Although more effects will be discussed later in Sect. 6, there is always the possibility that the SED of galaxies in the rest-UV is not a perfect power law (Bouwens et al. 2012), a behavior that was previously observed in the VANDELS galaxies by Calabrò et al. (2021). Differential shifts between UV slopes derived from different methods (either photometry or spectroscopy) have been already reported in the literature (e.g., see Hathi et al. 2016).
Ionizing photon production efficiency: ion The ionizing photon production efficiency is defined as the total number of ioniz-7 C -F is a -based optimization algorithm included in the S P package for scientific analysis: in https://docs.scipy.org/doc/ scipy/reference/generated/scipy.optimize.curve_fit.html.
ing photons per unit time produced by a radiation field ( ) normalized by its intrinsic UV luminosity ( UV ). In practice (see C19): where 1500 is the intrinsic (dust-free) SED luminosity at 1500Å (in units) and is the number of ionizing photons produced by the best-fit stellar population. is calculated as the integral of the intrinsic SED below the Lyman limit ( < (H 0 ) = 912 Å, or equivalently > (H 0 ) = 13.6 eV): When averaged over the whole galaxy population and integrated over the full range of UV luminosities, ion can be used to compute the total emissivity of ionizing photons ( ion ) at any redshift (see Eq. 1).
For typical SF galaxies that can be described through a singleburst of star-formation, S99 templates predict an exponentially declining ion with age, whose log ion values span over 26−24 Hz/erg for 1 − 10Myr stellar populations at any metallicity, dropping dramatically towards ages older than 10 Myr (C19). In contrast, mixedage models can provide systematically higher ion for a given mean age than instantaneous burst models at the same evolutionary stage. Accordingly, our mixed-age fits give a median ion of log ion = 25.30 +0.28 −0.42 Hz/erg using the R16 law. ion does not strongly depend on the assumed attenuation law because it has been calculated as the ratio between two absolute quantities ( and UV ). Since the fitted stellar ages and metallicities are similar for both R16 and SMC (as they are primarily fixed by single spectral features), the shape of the intrinsic SED is also similar for both laws and therefore ion stays unaltered.
Intrinsic 900-to-1500Å flux ratio: (F 900 /F 1500 ) int The ionizing-to-nonionizing flux ratio, namely at 900-to-1500Å, depends on the physical properties of galaxies like the mean stellar age, metallicity, IMF and star formation history. Following C19, S99 single-burst bases set a limit of (F 900 /F 1500 ) int ≤ 2, exponentially declining with age down to negligible values at ages older than 20Myr, where int 900 ≈ 0. Again, a mixed-age population could in principle increase the flux ratio at intermediate ages with respect to a single-burst of star-formation. The median of the 900-to-1500Å flux ratio for the VANDELS sample results in 900 / 1500 int = 0.66 +0.41 −0.32 . This quantity does not show any dependency on the attenuation law for the same reasons as ion .
In photometric campaigns, when searching for LCEs at high-(e.g., Grazian et al. 2016Grazian et al. , 2017, authors usually relate the observed to the intrinsic 900-to-1500Å flux ratio in order to infer the relative LyC escape fraction of galaxies. Usually, these works only have access to the observed but not the intrinsic ratio, and they end up by assuming a constant value typically set by stellar population models predictions (e.g, Steidel et al. 2001). In the same works, a more common way to look at this quantity is indeed the 1500-to-900Å flux ratio but in units. Typically assumed values in the literature are 1500 / 900 int = 3 − 5 (see e.g., Guaita et al. 2016). We obtain a median value of 1500 / 900 int ≈ 4.2 for the VANDELS sample at 3 ≤ ≤ 5 (Fig. 4), in agreement with the previous studies. Fig. 4 offers a summary of the above-mentioned distributions: Age, ★ , B−V , 1500 spec , ion , and (F 900 /F 1500 ) int ; where the values resulting from our SED fits using either R16 or SMC law are shown in red and dark-blue, respectively. Error bars cover the 16 ℎ to 84 ℎ percentiles of each distribution, whose values have been quoted with respect to the median in the current section. spec , log ion , (F 900 /F 1500 ) int , measured R(LIS) and predicted log abs esc . Results derived using R16 or SMC attenuation laws are plotted in red and dark-blue, respectively. The error bars on top of each histogram encompass the 16 ℎ , 50 ℎ and 84 ℎ percentiles. The distributions of the mean residual flux of the LIS lines are shown in different colors because, by definition, this quantity does not depends on the attenuation law. Dashed vertical lines indicate log ion (Hz/erg) = 25.2 and abs esc = 5%.

The picket-fence model
The observation that neutral lines of hydrogen and other strong lowionization state (LIS) lines do not become black at minimum depth (or maximum optical depth) suggests a partial covering of the stellar continuum sources by the same cold, neutral and low-ionized gas (Heckman et al. 2001(Heckman et al. , 2011. In this scenario, it is expected that the residual flux of the absorption lines correlates with the fraction of LyC photons that escape the galaxy via uncovered channels. The commonly used 'picket-fence' model (Reddy et al. 2016b;Vasei et al. 2016) connecting the UV absorption features of a galaxy with the escaping ionizing radiation assumes a galaxy described by a patchy, ionization-bounded ISM (Zackrisson et al. 2013), where both the neutral and low-ionized enriched gas are distributed in highcolumn-density regions surrounding the ionizing radiation field (clumps). The fraction of sight-lines which are optically thick to the transition along these dense clouds is usually parametrized by the so-called covering fraction, ( ). Optically thick gas absorbs all of the continuum light at a given velocity whereas optically thin gas transmits all of the continuum. If the dust is homogeneously distributed as a foreground screen on top of the stars, the residual flux of the lines can be simply related to the gas covering fraction as follows: This simple model also assumes that the lines are described by a single gas component or, in other words, that all velocity components of the gas have the same covering fraction. Additionally, if one accounts for the dust attenuation within the galaxy ( B−V ), the escape fraction of ionizing photons ( abs, LIS esc ) can be predicted from the depth of the neutral and other LIS absorption lines using the following formulae (Reddy et al. 2016b;Gazagnes et al. 2018;Steidel et al. 2018): where the B−V is the UV dust-attenuation parameter measured according to the methods described in the previous Sect. 3.1, using a uniform screen geometry (see Eq. 2), and is the measured residual flux of the transition. The [ , ] coefficients correspond to the LIS to H lines residual flux conversion (see Gazagnes et al. (2018), Reddy et al. (2022) and SL22).
This methodology has been successfully validated against direct measurements of the absolute escape fraction of low-LCEs in Chisholm et al. (2018) and Gazagnes et al. (2020) and recently thanks to the advent of the LzLCS by SL22, where we pointed to this method as a good predictor for the real abs esc of galaxies at any redshift. The depth of the LIS lines have also been applied to predict the escape fraction of high-galaxy composites in Steidel et al. (2018) and Pahl et al. (2021), and some individual highspectra in Chisholm et al. (2018), with reasonable agreement with the observed abs esc values. However, this approach is not without of caveats (see Sect. 7.4 in SL22): the choice of the dust-attenuation law and the assumptions on the gas and dust geometry are some of the limitations of the picket-fence model. Indeed, as suggested by Mauerhofer et al. (2021), a 'picket-fence' gas/dust distribution is shown to be a very simplistic approximation to the real ISM geometry in state-of-theart galaxy simulations. In principle, all these assumptions and model dependencies can contribute to explain the scatter seen in the observed R(LIS) − abs esc relation (see SL22), but this model could still be a good approximation for unresolved high-studies, where the LyC emission coming from a single sight-line usually dominates.

Absorption line measurements
Our goal is to apply Eq. 6 to VANDELS spectra in order to indirectly estimate their ionizing escape fractions. To do so, we first measured the equivalent widths ( ) and residual fluxes ( ) for a set of ×4 UV LIS lines, namely Si 1260, O +Si 1302, C 1334 and Si 1527, which had simultaneous wavelength coverage in all the spectra at 3 ≤ ≤ 5. The equivalent width ( ) was then computed individually for every absorption line following Trainor et al. (2019): where is the observed spectral flux density and ★ ( ) is the modeled stellar continuum. The integration window (Δ ) was defined as ±1000kms −1 from the wavelength of minimum depth for the line in question. Then, the residual flux was measured as the median over a narrrow velocity interval of ±150kms −1 around the minimum flux of the line, or equivalently: One of the conditions of applicability of the picket-fence model requires that the column density of gas is large enough that the absorption lines are saturated in the curve of growth (i.e., optically thick limit). In order to test the condition of saturation, we performed a similar analysis to the one in SL22 (see also Trainor et al. 2015;Calabrò et al. 2022), comparing equivalent-width ratios for transitions of the same ion at different wavelengths. Most of the galaxies have a Si 1260-to-Si 1527 equivalent-width ratio which is compatible within the optically thick limits ( 1527 / 1260 ≤ 2.55) given by the theoretical curve-of-growth for these transitions (Draine 2011). We conclude that the ISM conditions are such that Si is optically thick (saturated), and we assume that C is also saturated for typical Si/C abundances (see discussion in Mauerhofer et al. 2021). Calabrò et al. (2022) independently found optically thick Si (log( SiII / cm −2 ) > 12.7) along a sample of VANDELS C emitters at similar redshifts, whereas higher ionization lines like Si 1400 were closer to the optically thin limit as they probe more rarefied gas.
The effect of resolution on the absorption lines When measuring absorption line properties from observed spectra, the low spectral resolution ( ) tends to make the lines broader and less deep, and the actual residual flux ( ) can be overestimated. To account for these effects we performed mock absorption line simulations. Considering the picket-fence model with a uniform dustscreen, we simulated Si 1260 absorption Voigt profiles assuming Gaussian distributions for the column density of Si ( SiII , cm −2 ), the Doppler broadening parameter ( , kms −1 ) and the line velocity shift ( , kms −1 ). The simulated spectra were then degraded to (VIMOS) ≈ 600 spectral resolution. Finally, the impact of the S/N was folded onto the simulations by assuming a constant S/N = 5 that matches the median of the S/N distribution of the VANDELS spectra. A set of ×100 simulations were run for every of the input covering fractions in the ( ) = 0−1 range, and the median measured residual flux of each distribution was then compared to the theoretical value of 1 − . See App. A for a more comprehensive view of our simulations.
As a result, a theoretical residual flux of (1 − ) = 0.1, 0.3 and 0.5 would be observed as R(LIS) = 0.25, 0.4 and 0.55 due to the smearing effect of the VIMOS instrumental resolution (i.e., 60%, 25% and 10% relative error), while a negligible correction will be applied at (1 − ) ≥ 0.6 because this regime is dominated by the S/N of the spectra. The resulting calibration (Eq. A3) was applied as a correction factor to the individually measured residual fluxes.

Indirect abs esc estimations
The mean residual flux of the LIS lines, R(LIS), is calculated by taking the average of the individual Si 1260, O +Si 1302, C 1334 and Si 1527) line's depths. If these lines were all saturated, their residual fluxes should be very similar and the mean is representative, allowing to gain S/N over the measurement. E.g, the average LIS residual flux of CDFS01735 (Fig. 3) is R(LIS) ≈ 0.69 (0.75 previous to the correction by resolution). Fig. 4 shows the histogram of the average R(LIS) distribution, characterized by R(LIS) = 0.31 +0.22 −0.23 (0.35 +0.19 −0.20 without resolution correction). Then, using the dust-attenuation parameter ( B−V ) provided by the F CUS SED fits, we applied Eq. 6 to every single galaxy in our sample. For the [ , ] coefficients in this equation, we chose the calibration by Gazagnes et al. (2018) i.e., [ , ] = [(0.63 ± 0.12), (0.44 ± 0.07)]. We remark that this calibration was obtained over a sample of emission line galaxies but, as stated in SL22, we do not expect this relation to change with galaxy type. As an example, the LyC escape fraction for CDFS017345 results in abs esc ≈ 16%, adopting the R16 attenuation law.
The predicted abs esc distribution following this method is also presented in Fig. 4. As expected, the inferred abs esc values show a clear dependency on the dust-attenuation law, with the median of the R16 distribution ×1.5 lower than the SMC law abs esc distribution. When using the R16 dust-attenuation law, the median ionizing absolute escape fraction, abs esc , for our VANDELS sample of 534 SFGs at 3 ≤ ≤ 5 is: while, when using the SMC attenuation law, it results in: abs esc = 0.03 ± 0.02 (0.04 ± 0.03), where the numbers in brackets correspond to the median abs esc previous to any resolution correction on the residual flux of the lines (see Eq. 6). As a comparison, Begley et al. (2022) obtained abs esc = 0.07 ± 0.02 for a similar sample of 3.5 VANDELS galaxies with deep VIMOS/U−band observations (see Fig. 4), that is, a factor of ×2 higher than our median value using the SMC law, although compatible within 1 uncertainty. Moreover, our average abs esc is in agreement with the extrapolations of the recent Trebitsch et al. (2022) model at 4. In Sect. 6, we will compare our results with the escape fraction derived from other surveys targeting LCEs in the literature, at lower and higher redshifts.

Composite spectra
Stacked spectra were built with the goal of increasing the S/N with respect to individual galaxy data, allowing us to clear out some of the underlying physical correlations between the different parameters in this study.  2022), all the individual spectra in the sample were first shifted into the rest-frame using the VANDELS spectroscopic redshift and then resampled onto a common wavelength range, specifically 1200 − 1925Å. According to the median redshift of the sample i.e., spec = 3.56, the resulting spectral binning was chosen to be 2.5Å/(1 + spec ) = 0.55Å (2.5 Å equals the VIMOS wavelength dispersion per resolution element). Before co-adding the spectra, they were normalized to the mean flux in the 1350 − 1370Å restframe interval. The final (normalized) flux at each wavelength was taken as the unweighted median of all the individual flux values after a regular 3 clipping in order to reject outliers and bad pixels. The uncertainty on the stacked spectrum was calculated via bootstrap resampling of the spectra included in the composite.
Guided by this scheme, we performed stacks in bins of UV magnitude ( 1500 ), UV intrinsic luminosity ( int 1500 ), UVcontinuum slope ( 1500 spec ) and Ly equivalent width ( Ly ). Concretely, we divided the sample according to the 25 ℎ , 50 ℎ and 75 ℎ percentiles (quartiles) of every property, resulting in four subsamples for each quantity sorted as Q1, Q2, Q3 and Q4. Then, the resulting stacks were processed through the F CUS code and all the secondary SED parameters described in Sect. 3.1 were obtained. App. B, Table B1 contains the main properties and inferred SED parameters ( abs esc and ion ) of the different composites.
The effect of stacking on the absorption lines Even though the equivalent widths of the H and metal lines do not change with the gas column density in an optically thick medium, the lines can be slightly broader or narrower depending on the gas thermal (Doppler) and turbulence velocities. More importantly, different galaxies may have different gas flow velocities which translates into red-or blueshifted line centers relative to the systemic velocity. Moreover, the use of the spectroscopic instead of the systemic redshift introduces an additional source of uncertainty in the position of minimum depth of the lines (see Llerena et al. 2022). All these effects contribute so that the UV-lines residual flux of the resulting galaxy composite can potentially be overestimated.
To correct for this smearing effect and making use of the simulations described in App. A, we randomly generated a set of = 100 simulated Si 1260 line profiles with different intrinsic gas properties (column densities, Doppler broadening, inflow/outflow velocities, etc.), but fixing the covering fraction (i.e., equivalent to one minus the residual flux). We then stacked the line profiles following the same methods and assumptions described in this section, where the effects of instrumental resolution ( (VIMOS) ≈ 600) and a constant S/N = 5 were also incorporated in the mock realizations. After that, we measured the depth of the composite Si 1260 line profile and compared it to the input value given by the covering fraction. According to the our simulations (Eq. A4), an input (1 − ) = 0.1, 0.3, 0.5 and 0.8 would require corrections factors as large as 70%, 40%, 25% and 10% due to the effect of stacking. Calabrò et al. (2022) showed that, although the bulk ISM velocity is globally in outflow, the average shift of the LIS lines is very close to the systemic velocity (−60 km/s), i.e., within the actual spectral resolution. Therefore, our simulations -which assume a single VIMOS resolution element (±150 km/s) as the standard deviation of the distribution for the velocity shift of the linesmight actually over-predict this effect. For that reason, we prefer not to apply such stacking corrections to our line measurements, but we encourage the reader to check App. A for more details.
It is also worth mentioning that Calabrò et al. (2022) do not find any correlation between the stellar mass or SFRs and the velocity shift of the ISM lines. Thence, we will not expect differential corrections on the residual flux when combining galaxies with very different masses or SFRs, a conclusion that can be extrapolated to other galaxy properties (e.g., UV magnitudes).

RESULTS
The following paragraphs describe on the main results of this paper. In Sect. 4.1, the global relations between the ionizing escape fractions and production efficiencies with different galaxy properties are shown on an individual galaxy-basis. In Sect. 4.2, the LCEs and non-LCEs composites are presented, and the differences in their non-ionizing rest-UV spectra are placed in the context of the physical ISM conditions which enable the ionizing radiation to escape. In Sect. 4.3, the global ionizing properties of the LCEs versus non-LCEs samples are discussed. Lastly, Sect. 4.4 summarizes the properties of previously known LCEs reported in the literature which were included in our sample.

The ionizing escape fraction and production efficiency dependence with observed galaxy properties
Fig. 5 investigates how our predicted ionizing absolute escape fraction ( abs esc ) and the ionizing photon production efficiency ( ion ) depends on galaxy properties. The individual VANDELS measurements are shown in the background together with the results from our stacking analysis on top (large symbols). Systematic differences in abs esc and ion due to the use of a shallower (R16, filled diamonds) or steeper (SMC, open diamonds) attenuation law are also explored. In general, the escape fraction values derived using the SMC law shift to slightly systematically higher escapes (by a factor of ×1.5) compared to the R16 ones.
A Kendall− statistics was applied to each pair of variables in order to figure out the strength ( ) and significance ( −value) of the correlation. For a sample of 500 objects, we consider correlations to be significant if . 10 −3 (3 ) and strong if | | 0.1. Coincidentally, all significant correlations studied in this work showed up to be strong correlations, so that in Fig. 5 we only indicate the coefficients for significant correlations (thick-framed panels), otherwise we write val. > 10 −3 .

Observed absolute UV magnitude
The predicted abs esc and the observed UV absolute magnitude ( 1500 ) for the VANDELS sample at 3 ≤ ≤ 5 do not show any clear correlation when considering the whole range of UV magnitudes. Tanvir et al. (2019) found no significant correlation of HI with galaxy UV luminosity across a wide range of redshifts. Assuming that the escape fraction is primarily regulated by the same optically thick neutral gas along the line-of-sight, as suggested by Eq. 6, a lack of correlation between abs esc and 1500 would therefore be expected, supporting our global trend.
However, fainter (low-mass) galaxies usually host lower gas and dust fractions compared to brighter systems (Trebitsch et al. 2017). Together with their low-gravitational potential and often bursty SFHs (Muratov et al. 2015), the expelling of gas and dust in a turbulent ISM (Kakiichi & Gronke 2021) favours the creation of holes through which the LyC photons can freely escape. In this scenario, a higher abs esc is naturally expected for the faintest galaxies (see Sect. 7.2 in SL22).
Even though the correlation is tentative, 1500 −20 and fainter VANDELS galaxies tend to have higher escape fractions (Fig. 5). This mild tendency has been reported through the study of galaxy composites at high-by the KLCS (Steidel et al. 2018;Pahl et al. 2021Pahl et al. , 2023 as well as observed at low-by the LzLCS (Flury et al. 2022a;Chisholm et al. 2022, and SL22). Contrarily, models such as Sharma et al. (2016)  . Scatter plots searching for correlations between the predicted ionizing photon escape fraction ( abs esc , first column), the ionizing production efficiency (log ion , second column) and the product of the two (log abs esc × ion , third column, see Eq. 1), versus different galaxy properties: observed and intrinsic UV absolute magnitudes ( 1500 , int 1500 ), UV-continuum slope at 1500Å from the best-fit SED ( 1500 spec ), and Ly equivalent width (−1 × Ly ). The Kendall ( ) correlation coefficients for individual R16 measurements (coloured dots in the background) are shown at the top but only for significant correlations (thick-framed panels). Results from our stacking analysis are displayed with filled (R16) and open (SMC) thick diamonds. Typical error bars for individual sources are plotted at the bottom part of each panel, and the arrows along the int 1500 panels measure the shift in int 1500 due to the use of the SMC dust-attenuation law. Grey-shaded regions mark assumed canonical values of abs esc ≥ 5% and log ion = 25.2 − 25.3 in classical Reionization models (Robertson et al. 2013).
ing escape fraction towards bright UV systems, in disagreement with our picket-fence formulation 8 . The ionizing production efficiency shows a very smooth (although non-significant) dependence on 1500 , with the stacked measurement at the faintest UV magnitude bin ( 1500 −20) being higher than the typically assumed canonical value for cosmic reionization of log ion (Hz/erg) = 25.2 (Robertson et al. 2013). Our results are in agreement with other studies in the literature, for example with the results by Bouwens et al. (2016) at similar UV magnitudes but, contrary to the former, our 1500 range is actually not wide enough to show a clear ion − 1500 trend. Even though, it has been shown that, at similar redshifts, the ion evolution with UV luminosity is a very smooth correlation (see Lam et al. 2019;Emami et al. 2020;Prieto-Lyon et al. 2022). A more complete picture on the ion − 1500 relationship will be given in Sect. 5.
In general, both abs esc and ion distributions with 1500 are characterized by a large scatter ( 1dex, e.g., 0.1% to 10% in abs esc ) and a weak evolution with 1500 , where only the faintest galaxies in the sample have tentatively higher values in abs esc (≥ 5%) and log ion (≥ 25.2 Hz/erg). Finally, the log abs esc × ion product preserves the overall evolution of abs esc with the UV magnitude.

Intrinsic (dust-free) UV luminosity
The intrinsic UV absolute magnitude for each galaxy ( int 1500 ) was calculated from the best-fit SED by taking the flux at 1500Å previous to attenuation by dust (i.e., the dust-free SED), and computing the AB magnitude via the usual distance modulus formulae.
On the one hand, abs esc versus int 1500 shows one of the strongest correlations of the present study ( ≈ 0.4), where the less intrinsically luminous galaxies have the highest abs esc . From our abs esc prescription (Eq. 6), this correlation is expected since the dustattenuation is by definition directly related to the escape fraction of galaxies, where the most UV-bright galaxies are attenuated (Finkelstein et al. 2012;Bouwens et al. 2014) and host substantially higher, more extended gas and dust reservoirs at the same time (see e.g., Ma et al. 2020, from the perspective of the FIRE-2 simulation). This behavior was previously reported in Begley et al. (2022), where they demonstrate how intrinsically UV-faint galaxies at 3.5 would require statistically higher escape fractions in order to reproduce the observed distribution of the ionizing-to-nonionizing flux ratio.
On the other hand, the ion versus int 1500 distribution appears completely flat irrespective of the int 1500 bin and the dustattenuation law. According to Eq. 4, the dust correction factor applied to any ion measurement would be: A <912 /A 1500 . This ratio does not strongly depend on the dust-law, thus no ion dependence on the attenuation will be expected neither. Therefore, the abs esc × ion product inherits the same dependence on the UV luminosity as abs esc does. When comparing the R16 and SMC stacks predictions (see Fig. 5), their abs esc × ion show the largest differences at the bright int 1500 end i.e., the more deviation appears for the more attenuated, redder galaxies. 8 This said, the hypothesis proposed in Naidu et al. (2020) of a reionization driven by −20 ≤ 1500 ≤ −18 galaxies (what they called the 'oligarchs'), cannot necessarily be ruled out by our observations because our galaxies do not span beyond UV ≥ −19 AB.

UV-continuum slope at 1500Å
Standing out as the strongest correlation of this study ( ≈ −0.6, see Fig. 5), the abs esc inversely scales with the best-fit UV-continuum slope at 1500Å ( 1500 spec ) so that the bluest galaxies in the sample emit a significantly larger fraction of ionizing photons than their redder counterparts. The "Dust Only" case (dotted line) follows the resulting Eq. 6 assuming there is no gas along the line of sight, so the dust is the only source of attenuation for LyC photons. This scenario, although provides a physical upper limit to the previous equation, does not expect to hold observationally. Because the dust and cold gas within the ISM are spatially correlated, galaxies which are more dust-attenuated will also have lower residual fluxes of the LIS lines, and they will deviate from the "Dust only" case towards low abs esc values, as we see in Fig. 5. This trend has been directly observed by Flury et al. (2022b) in the LzLCS survey, also lately investigated by Chisholm et al. (2022). Additionally, a decrease in abs esc with increasing UV colors was shown in Pahl et al. (2021) using KLCS galaxy composites at high-. Also interestingly, Begley et al. (2022) found that UV blue galaxies at 3.5 would require statistically higher escape fractions in order to reproduce the observed ionizing-to-nonionizing flux ratio distribution.
Regarding our indirect approach to predict abs esc (Eq. 6), a negative correlation between the escape fraction and 1500 spec must be expected by definition, since the UV slope is inherently linked to the dust-attenuation ( B−V ), so that redder SEDs usually means more attenuated galaxies (Meurer et al. 1999). For example, we obtain abs esc ≥ 5% for galaxies whose 1500 spec ≤ −1.5, decreasing to abs esc ≈ 0.1% at 1500 spec ≈ −0.5. As stated in Chisholm et al. (2022), this steep relation between the escape fraction and the UV slope is particularly important at the highest redshifts because (1) it can be easily used as an observational proxy to indirectly infer abs esc at the EoR (see Trebitsch et al. 2022), and (2) primordial/fainter galaxies are thought to host less dust than actual SF galaxies (Finkelstein et al. 2012;Bouwens et al. 2014;Cullen et al. 2023), which naturally explains why higher redshift galaxies emit more ionizing photons than their lower redshift counterparts (Finkelstein et al. 2019). As we will discuss in Sect. 6, our abs esc − 1500 spec results agree with the relation = 0.3 published by Chisholm et al. (2022) but extrapolated to redder UV slopes and lower escape fractions. This similarity suggests that the abs esc − 1500 spec does not strongly change with redshift. The ion parameter is also strongly correlated with 1500 spec , so that ion rapidly increases as the UV slope decreases, only giving log ion (Hz/erg) ≥ 25.2 for galaxies whose slopes are 1500 spec ≤ −1.5. Our results are in agreement with other studies in the literature, for example with the results by Bouwens et al. (2016) and Matthee et al. (2017a), but shifted to redder 1500 values. For more details on the ion − 1500 relationship and comparison with the literature, we refer to Sect. 5.
Qualitatively, the largest differences between the R16 and SMC derived abs esc values are found among the redder objects i.e., the most attenuated galaxies in the sample. However, both R16 and SMC derived ion values are similar irrespective of the dust-attenuation law.

Ly equivalent width
The relation between the properties of the Ly 1216 line and the escape and production efficiency of ionizing photons is also striking. This can be seen in Fig. 5, where the predicted abs esc monotonically raises as the Ly equivalent width also increases ( Ly ). For typ-ically assumed values of Ly ≤ −20Å in LAEs (e.g., Pentericci et al. 2009, but see Stark et al. (2011); Kusakabe et al. (2020)), we obtain abs esc ≥ 5%. This result is in concordance with previous results in the literature. As demonstrated by SL22 (but see also Gazagnes et al. 2020;Izotov et al. 2021), the strongest LCEs are usually among the strongest LAEs i.e., having the highest Ly equivalent widths.
With the exception of photon scattering and the effect of gas kinematics, the same mechanisms that regulates the escape of ionizing photons also regulate the leakage of Ly photons and therefore the shape and strength of the Ly line (e.g. Henry et al. 2015;Verhamme et al. 2015). Considering the aforementioned picket-fence model with a very simplistic assumption on the geometry of galaxies, these main mechanisms are (1) the neutral gas column density -whose optically-thick spatial distribution is parametrized in Eq. 6 via the covering fraction -and (2) the dust-attenuation term ( B−V ). In this way the Ly emission ( Ly ) and the escape of LyC photons ( abs esc ) would remain physically related (see Verhamme et al. 2017;Izotov et al. 2020;Flury et al. 2022b;Maji et al. 2022). In fact, studying the afterglow gas of a significant sample of Gamma Ray Burst (GRB) host galaxies, Vielfaure et al. (2021) found that the gas columns density of such gas was indeed optically thick, so that "the bulk of Ly photons produced by massive stars in the star-forming region hosting the GRB will be surrounded by these opaque lines of sights".
In order to illustrate the role of dust attenuation and the gas covering fraction in the escape of Ly and LyC photons, in Fig.  6 we plot the average LIS absorption profile for different stacks in bins of Ly equivalent width (see also Trainor et al. 2019). This plot shows the combined profile of several ISM LIS lines (Si 1260, O +Si 1302, C 1334, Si 1527), color-coded by the inferred UV dust-attenuation ( B−V ) and the predicted escape fraction for each composite ( abs esc ) is indicated in the legend. The resulting Ly profiles can also be seen in the inset, and finally Ly (in Å) is plotted against the LIS equivalent width of the combined LIS absorption profiles ( LIS , in Å).
As a result, as long as the dust-attenuation increases and the gas covering decreases (increasing residual fluxes), the resulting escape fractions and the Ly equivalent widths also decrease, from Ly values typically found in LAEs towards absorption, even damped, Ly profiles. Relatedly, the relation between the Ly and the strength of the LIS lines has been widely studied in individual and stacked rest-UV spectra of LBGs at = 2 − 5 (Shapley et al. 2003;Jones et al. 2012;Henry et al. 2015;Trainor et al. 2015;Du et al. 2018;Pahl et al. 2020), and our results agree with the overall picture described in these papers. On the one hand, the observed connection between the LIS equivalent width and the reddening can only be explained if the dust resides within the same clouds as the neutral gas in a clumpy ISM geometry (SL22). On the other hand, the dwindling of the Ly emission with the dust-attenuation is a a natural consequence of the Mass Metallicity Relation (MZR, see Maiolino & Mannucci 2019, for a review), and the empirical relation between the Ly strength and the stellar mass (Cullen et al. 2020).
A large Ly can also be ascribed to a boost in the production of ionizing photons (Nakajima et al. 2018a). The VANDELS high-Ly − ion relation can be found in the bottom-mid panel of Fig. 5, where a significantly higher log ion (Hz/erg) ≥ 25.2 is found for the strongest LAEs only, with a Ly ≈ −30Å. Thanks to JWST data in combination with HST photometry, an increasing log ion with Ly strength has also been shown for individual galaxies at the late-edge of the EoR (Ning et al. 2022). However, in the works by Cullen et al. (2020) and Reddy et al. (2022), ion has been demonstrated to be insufficient to solely account for the whole variation in Ly , rather "the covering fraction of optically-thick H (gas) appears to be the principal factor modulating the escape of Ly , with most of the Ly photons in down-the-barrel observations of galaxies escaping through low-column-density or ionized channels in the ISM", as stated in Reddy et al. (2022).
Similarly, our abs esc × ion − Ly relationship results from the convolution of the abs esc and ion dependence on the Ly strength, where the abs esc − Ly behavior plays the major role. The higher ion values for the objects with the highest Ly equivalent width in the sample (LAEs) naturally increase the abs esc × ion product (Matthee et al. 2022).

Stellar mass
The predicted abs esc strongly scales with the stellar mass ( ★ ) of the VANDELS galaxies so that low-mass galaxies have statistically higher escape fractions than more massive systems (Fig. 7). For instance, only log ★ ( ) < 9 VANDELS galaxies show cosmologically relevant escape fractions of abs esc ≥ 5% (Robertson et al. 2013, see grey-shaded area in Fig. 5), and these are actually 1dex higher than the average escape fraction of the most massive galaxies in the sample (log ★ ( ) ≥ 10). According to our picket-fence model, a physical abs esc − ★ connection is expected (Eq. 6) because more massive galaxies are intrinsically more attenuated than low-mass systems (see e.g., Finkelstein et al. 2012;McLure et al. 2018a;Fudamoto et al. 2020). As log ★ increases, the dust-attenuation increases and the line-depth decreases accordingly, and both yield progressively decreasing escape fractions ( abs esc ). This is consistent with the results by Reddy et al. (2022) where, using a joint modeling of the composite FUV and optical spectra of high-SFGs, they demonstrate how the H

Reddy+ law SMC law
Galaxies dominate the ion. budget at EoR (canonical ξ ion ) Figure 7. Relation between our derived abs esc and the stellar mass ( ★ , in ). The abs esc − ★ correlation is strong and significant for both R16 and SMC dust-laws, although systematically higher escapes are reported when using steeper, SMC-like attenuation laws. The layout is the same as in Fig.5. and the gas-enriched covering fraction decreases with the galaxy stellar mass.
So far, there is no consensus in simulations on whether more massive galaxies should emit a higher (e.g., Naidu et al. 2020) or lower (Rosdahl et al. 2022) fraction of their produced ionizing photons to the IGM compared to their less massive counterparts, with some of them even suggesting a turnover at intermediate masses (see Ma et al. 2020). In observations, while at low-there does not seem to be any clear relation between these two quantities (Izotov et al. 2021;Flury et al. 2022b Lately, Pahl et al. (2023) have also shown a negative and significant abs esc − log ★ trend in the 9 ≤ log( ★ / ) ≤ 10 stellar mass range, using 3 KLCS stacks. However, at similar redshifts, Begley et al. (2022) did not find any statistical distinction when splitting their sample into lower and higher masses with respect to the median, suggesting that ★ is at best a secondary indicator of the average abs esc . In the 8 ≤ log( ★ / ) ≤ 10 mass interval, our stacking analysis with VANDELS support the former scenario by which low-mass galaxies have higher escape fractions (similar to Ma et al. 2020;Pahl et al. 2023). In this vein, our picket-fence formalism put the physical picture suggested by models like Naidu et al. (2020) up against the ropes, since the latter yields a monotonic increase of the escape fraction with stellar mass, behaviour also disfavoured by other works in the recent literature (see compilation by Pahl et al. 2023).
Finally, whilst AGNs may not contribute directly to the ionizing photon budget at the EoR (e.g., Hassan et al. 2018), semianalytical simulations by Seiler et al. (2018) show that they could indirectly influence reionization by clearing out holes in the ionized regions, so that abs esc can be boosted after quasar wind events. In such scenario, the mean escape fraction peaks for intermediatemass galaxies, around log( ★ / ) ≈ 8. Regrettably, our sample is not able to probe this hypothesis, since it does not extend to lower stellar masses.

Non-ionizing rest-UV properties of potential LCEs at
3 ≤ ≤ 5 Having predicted the ionizing escape fraction ( abs esc ) for every galaxy allows us to split the sample into potential LCEs, with abs esc ≥ 0.05 (64 out of 534 galaxies, i.e., ≈ 10% of LyC detection rate), and non-LCEs with abs esc < 0.05 (470 galaxies, assuming R16). We then build composite spectra for LCEs and non-LCEs candidates, following the same method as described in Sect. 3.3. The main properties, fitted SED parameters and UV lines measurements are summarized in App. B, Table B2. Both LCEs and non-LCEs stacks are presented Fig. 8, together with a handful of insets and labels which highlight the differences in their non-ionizing FUV spectra.

UV-continuum slope at 1500Å
The UV slope for LCEs and non-LCEs are remarkably different: 1500 spec = −2.17 ± 0.03 for LCEs versus −1.23 ± 0.01 for non-LCEs composites (see Fig. 8). As seen in the previous section, the UV slope is by definition related to abs esc in the sense that it constitutes a proxy for dust attenuation at UV wavelengths , therefore favouring a picture in which bluer galaxies with lower levels of dust attenuation ( B−V ) display higher values of escape fraction (by construction, see Eq. 6). Having different levels of dust attenuation in LCEs and non-LCEs will partially influence the strength of other nebular emission lines in their UV spectra.
The intrinsically dustier nature of non-LCEs in opposite to the LCEs population has been explicitly reported in SL22 and further investigated in Chisholm et al. (2022) using LzLCS data at 0.3. Also recently, Pahl et al. (2023) showed a monotonic increase of abs esc with B−V from FUV SED fitting in KLCS composite spectra. Likewise, Pahl et al. (2021) and Begley et al. (2022) observationally demonstrated a decrease in abs esc with increasing UV colors at 3, using HST (KLCS) and ground-based (VANDELS) LyC imaging, respectively. All these different LyC data sets at low-and highwill be put together in Sect. 6.

Ly emission
LCEs clearly show a stronger Ly emission than non-LCEs ( Fig.  8): Ly = −29.71 ± 2.46Å against −3.62 ± 0.41Å for LCEs and non-LCEs stacks, respectively. Worth mentioning is that the Ly for the LCEs composite is compatible with the typical definition of LAEs in the literature (i.e., Ly ≤ −20Å see Pentericci et al. 2009), so that the LCEs composite is mostly composed of the Ly emitting galaxies in the VANDELS sample. The non-LCEs composite, however, shows little Ly emission compared to LCEs, with a damped Ly red-wing extending up to 1240Å. This suggests a higher column density of H gas beneath the bulk of the non-LCEs compared to the LCEs population, an unequivocal necessary condition for preventing the leak of Ly and LyC photons (see Henry et al. 2015). This is also compatible with the strong decrease of the Ly equivalent width (by a factor 9) between LCEs and non-LCE, which is stronger than the decrease of ion (a factor 4), and which is naturally explained by the impact of a higher column density and higher dust content on the Ly escape fraction (see e.g. Atek et al. 2014).
Observational evidences for an increase in Ly with increasing abs esc has been reported by a few recent studies: Flury et al. (2022b) through FUV spectroscopy of LzLCS sources at 0.3; Fletcher et al. (2019) and Pahl et al. (2021Pahl et al. ( , 2023

Other nebular emission lines
In Fig. 8, the C 1550, He 1640 and C ]1908 nebular lineprofiles are displayed for the LCEs (in blue) and non-LCEs composites (in red). The C 1550 and He 1640 lines for LCEs show slightly stronger emission and narrower profiles, while the equivalent width of the C ]1908 line is remarkably higher in the LCEs stacked spectra compared to the non-LCEs one. We obtain HeII = −1.11±0.25Å (−0.81±0.06Å) and CIII] = −2.15±0.67Å (−0.73 ± 0.13Å) for LCEs (non-LCEs).
Similarly, Naidu et al. (2022) stacked the UV spectra of potential LyC emitting candidates according to their Ly line properties, showing a very different He 1640 profile respect to the non-emitting candidates (probably attributed to sample selection). Contrarily, in the recent work by Marques-Chaves et al. (2022b) based on direct LyC measurements of low-LzLCS galaxies, the authors do not find any significant correlation between the LyC escape fraction and the spectral hardness, where the He 1640 intensity is fundamentally driven by changes in the metallicity rather than by abs esc . Therefore, the slightly higher observed equivalent widths in He 1640 and O ]1666 for the LCEs composite (see Llerena et al. 2022) may indicate a much dissimilar gas-phase metallicity between the LCEs and non-LCEs stacks than the stellar metallicities that are actually derived from our stellar population modeling (see Sect. 4.3).
Regarding the UV carbon nebular lines and based on the work by Schaerer et al. (2022a); Saxena et al. (2022b) and Mascia et al. (2023b), we compute the C 1550 over C ]1908 flux ratio, and we obtain C /C ] = 0.79 ± 0.21 for LCEs while C /C ] = 0.42 ± 0.20 for non-LCEs (the C 1550 doublet is blended at the current resolution). Schaerer et al. (2022a) empirically demonstrated that both C 1550 and C ]1908 lines tend to be stronger in the LzLCS galaxies with the highest abs esc , and proposed a C /C ] ≥ 0.75 threshold for significant LyC leakage which, once again, is in agreement with our stacking analysis. As we will see in the next section, a combined effect of escaping ionizing radiation (favoured neutral gas and dust geometry) with a higher ion parameter is responsible for enhancing the emission of high-ionization-state lines such as C 1550 and C ]1908 in LCEs. Fe 1608, Al 1670 nd Al 1858 doublet) are detected in our VAN-DELS composite spectra (Fig. 8).

Absorption lines and ISM properties
The N 1240 stellar-wind line in LCEs shows a characteristic P-Cygni profile that can be fully reproduced by stellar templates (excluding potential contamination by AGNs), and it indicates the presence of a young stellar population (≤ 5 Myr). Conversely, N 1240 is suppressed by the absorbing part of the damped Ly red-wing in non-LCEs, indicating a high H column density. Going to redder wavelengths, the Si 1393,1402 high-ionization state absorption lines, which share both a stellar and ISM origin, have the same equivalent width in the LCEs and non-LCEs composites, suggesting a homogeneous distribution of the diffuse high-ionized species in the ISM of LBGs.
The remaining set of LIS absorption lines in LCEs (Si 1260, O +Si 1302, C 1334, Si 1527, Fe 1608, Al 1670 and Al 1858 doublet) clearly show lower equivalent widths than in the non-LCEs stack. For example, the S 1260 and C 1334 in Fig. 8, SiII = 0.99 ± 0.13 (1.82 ± 0.05) for LCEs (non-LCEs), while CII = 0.64 ± 0.13 (1.74 ± 0.04) in the LCEs (non-LCEs) composite. Following the picket-fence model, the fact that strong LCEs have significantly weaker absorption lines, i.e., lower equivalent widths with higher residual fluxes, can be attributed to a lower covering fraction of the spatially co-existing neutral and low-ionized gas, so that the photons from a given transition escape through low column-density channels in the ISM, as well as LyC photons do (Gazagnes et al. 2018. The observed connection between LIS absorption lines and the strength of the Ly line is also observed in our galaxy composites, and it can be explained by the picket-fence formulation, in which the gas covering fraction along the line-of-sight primarily governs the escape of the absorbed LIS, Ly and LyC photons. This behavior is well summarized in Fig. 9, where the Ly is plotted against the measured (averaged) residual flux of the LIS lines, R(LIS), and the points are color-coded by B−V . Blue and red symbols correspond to our LCEs and non-LCEs galaxy composites, respectively.
As long as the line-of-sight covering fraction of gas decreases and the dust-attenuation increases (traced by higher R(LIS) and lower B−V values), the Ly equivalent width dramatically increases towards values typically found in LAEs. For being the case, the stronger LAEs, which has been previously suggested in this work to be a proxy for the bulk of the LCEs population, fall in the top-right part of the plot.
This trend matches the results by Steidel et al. (2018) in the framework of the KLCS, and echoes the main conclusions by Gazagnes et al. (2020), using individual LyC detected galaxies at low-. In the same Fig. 9, the median results for LCEs and non-LCEs in the LzLCS (Flury et al. 2022a), are indicated by white open hexagons (see also SL22). As we can see, the LzLCS results clearly agree with our VANDELS stacks, which reassure the use of the picket-fence approach for describing the ISM of LCEs. Additionally, the results by Jones et al. (2013) for individual LAEs at = 2 − 4 are shown on top of the VANDELS results, following the general trend.

Ionizing properties of potential LCEs at 3 ≤ ≤ 5
The ionizing properties of LCEs and non-LCEs can be accessed by extrapolating the non-ionizing best-fit SED to ionizing wavelengths. In particular, it is of our most interest to check whether potential differences in the ionizing-to-nonionizing flux ratio, (F 900 /F 1500 ) int , and production efficiency, ion , exist.
First, we input the LCEs and non-LCEs stacked spectra into the F CUS code in order to unveil the average properties of the

LCEs non-LCEs
Intrinsic ionizing SED Figure 11. The ionizing continuum of potential LCEs at 3 ≤ ≤ 5. Left: a comparison between the dust-attenuated SED for LCEs (blue) versus non-LCEs (orange) down to ionizing wavelengths, normalized at 1360Å. The inset shows the intrinsic ionizing spectrum (i.e., dust-free) of LCEs and non-LCEs (in blue and red, respectively). As written in the legend, the intrinsic 900-to-1500 flux ratio (in F units) and the ionizing photon production efficiency ( ion ) are higher for LCEs. The effect of dust-attenuation and neutral hydrogen absorption in both ionizing spectra is sketched through the grey and black arrows, where the blue and orange points represent the corresponding observed LyC flux once both absorptions have been accounted (a constant H cross-section with wavelength has been considered). Right: the AB magnitude distribution of the predicted observed flux in the LyC for the LCEs (blue) and non-LCEs (orange) samples, once multiplied by the average IGM transmission at = 4 of IGM = 0.3.
underlying stellar population of each galaxy sample. Then, we use the best-fit stellar continua to define different galaxy observables (see Sect. 3.1). As we can see, the fit to LCEs provides a lower stellar age (13 ± 3 Myr) than for non-LCEs (19 ± 1 Myr), tentatively suggesting younger stellar populations for LCEs. On the other side, stellar metallicities are very similar for both (Z ★ ≈ 0.2 Z ), and compatible with the median of the total sample. These values can be seen in Fig. 10, where the best-fit F CUS SEDs are superimposed to the LCEs and non-LCEs spectra.
The inferred dust-attenuation parameter is remarkably dissimilar for LCEs and non-LCEs, being B−V ≈ 0.07 and 0.23 mag., respectively. This is by construction, since the parent sample has been splitted based on the abs esc of individual VANDELS galaxies (Eq. 6), which is explicitly related to B−V . It reflects the expected less dusty nature of LCEs compared to the bulk of the LBG population, as we explained in the previous section.
The intrinsic (i.e., dust-free) ionizing spectra of the resulting best-fit SED for LCEs and non-LCEs is presented in the inset of Fig. 11 (left) through blue and red solid lines, correspondingly. The slight decrease in age and metallicity makes the intrinsic ionizingto-nonionizing flux ratio (900-to-1500Å, see C19) to grow up to (F 900 /F 1500 ) int = 0.64 in LCEs from 0.47 for non-LCEs. A different (intrinsic) 900-to-1500Å ratio due to an decrease in age also fosters a rise in the ionizing efficiency to log ion (Hz/erg) = 25.38 for LCEs, while it is 25.18 for the non-LCEs composites. In combination with the covering fraction of neutral gas, the ion parameter is partially responsible for shaping the strength of highionization-state lines such as C 1550 and C ]1908, and other nebular resonant lines such Ly .
Ultimately, we infer what the observed AB magnitude at LyC wavelengths would look like for the LCEs and non-LCEs populations, as a guidance for upcoming surveys targeting LyC emitting galaxies at similar FUV magnitudes. According to the fundamental definition of the absolute ionizing escape fraction of galaxies (see Izotov et al. 2016a, SL22), the observed flux close to the Lyman edge ( 900 ) can be simply predicted by doing: 900 = abs esc × int 900 , where int 900 is the intrinsic LyC flux or, in other words, the ionizing flux previous to any gas or dust absorption (in units). Since we have access to the intrinsic ionizing SED through the F CUS fits, and the abs esc has been predicted following Eq. 6, we compute 900 for every galaxy in the VANDELS sample. We then derive the AB apparent magnitude from the 900 synthetic flux measurements ( LyC , AB). The resulting distribution can be seen in Fig. 11 (right), where the blue (orange) histogram represent the LyC AB distribution for LCEs (non-LCEs), once we manually apply a IGM = 0.3 mean IGM transmission factor, i.e., the mean IGM at = 4 according to Inoue et al. (2014).
The consecutive effects of neutral H gas absorption and dust attenuation in the LCEs and non-LCEs ionizing spectra are represented in Fig. 11 (left) through the black and grey arrows (assuming a constant H cross-section with wavelength), whilst the blue and orange open circles denote the LyC flux measured at 870 − 910Å. After applying a IGM = 0.3 factor and translating into AB magnitudes, these two points would equal the median of the individual LyC AB magnitude distribution shown in the right panel of the same figure.
Clearly, the combined effects of a lower ionizing-tononionizing flux ratio plus lower escape fractions for non-LCEs makes their median output LyC flux dwindle 2mag. down with respect to a hypothetical LCEs population drawn from the same LBG sample. Therefore, a dedicated survey which attempts to detect the bulk of LCEs at 3 ≤ ≤ 5 would required a sensitivity such that it reaches, at least, LyC 31 AB at LyC wavelengths.
Using very deep VIMOS/U-band imaging of the CDFS field (1 limiting magnitude of 30.4 AB), Begley et al. (2022) estimated 10 to 15 LyC detections at 2 for sources whose redshift is 3.35 ≤ ≤ 3.95, assuming an average abs esc ≈ 0.07. Restricting to the CDFS only, we found 7 sources with predicted LyC ≤ 30 AB in the same redshift range. The differences can be attributed to our lower average escape fraction compared to the former work. This said, only two robust LyC emitters have been detected so far in the CDFS field (see following section). The stochasticity of the IGM opacity might reduce our expectations.

Confirmed LCEs in the VANDELS sample
In Fig. 9, we highlighted two LyC emitting galaxies included in our VANDELS sample which have been already published in the literature, in particular within the CDFS field (see Begley et al. 2022): CDFS012448 (Saxena et al. 2022a) and Ion1 Ji et al. 2020).
While both galaxies share a similar SED dust-attenuation ( B−V 0.2 mag.), we measure a much lower residual flux in CDFS012448 (R(LIS) 0.4) than in Ion1 ( 0.55), yielding to a predicted escape fraction of abs esc 3% versus 6%. It is worth stressing that our inferred abs esc value for Ion1 equals the reported value by Ji et al. (2020). In the case of CDFS012448, our abs esc is remarkably different to the one quoted in Saxena et al. (2022a) of ≈ 20%, possibly because the different approaches between papers. Our method was able to identify both galaxies as potential LCEs.
Additionally, these galaxies exhibit very different Ly properties between each other: while CDSFS012448 shows a relatively strong Ly in emission ( Ly −20Å), the Ion1 Ly profile appears in absorption ( Ly 1Å). In our VANDELS sample, 10% of the LCEs ( abs esc ≥ 5%) appear with Ly in absorption. This manifests, once again, the variety of Ly profiles that can be found among the LCEs galaxy population Flury et al. 2022b;Naidu et al. 2022) where, given the fact that Ly can be resonantly scattered whereas LyC cannot be scattered, both Ly and the ionizing emission do not necessarily come from the same spatial location within the host galaxy (e.g., Vanzella et al. 2012). A low IGM transmission could also, in principle, have killed out most of the Ly photons of Ion1.

THE PRODUCTION EFFICIENCY OF IONIZING PHOTONS ( ion ) IN HIGH-REDSHIFT GALAXIES
Here, we investigate the evolution of ion as a function of the UV absolute magnitude and the UV slope (see Bouwens et al. 2016;Chisholm et al. 2022), quantities which will be easily measured for a significant fraction of galaxies at the EoR thanks to upcoming gound-based (ELT, GMT, TMT, ...) and ongoing space-based facilities (JWST). Fig. 12 shows the ion parameter in stacks of the absolute UV magnitude for the VANDELS sample. As it was introduced in Sect. 4 and regardless of the attenuation law (see inset), the ionizing production efficiency monotonically, smoothly increases towards fainter UV magnitudes, so that log ion ≥ 25.2 Hz/erg (Robertson et al. 2013) at 1500 −20, i.e., the faint-end of our sample. Stacked results from other samples in the literature are also included in Fig. 12. The main comparison samples for our case of study are the MOSDEF galaxies from Shivaei et al. (2018) Bouwens et al. (2016) and Matthee et al. (2017a), both tracing photometrically-selected HAEs at similar 1500 than ours. Nevertheless, VANDELS samples a narrower range in 1500 so that the overall trend does not appear as clear as in the former works. At the faintest UV magnitude bin, our ion − 1500 relationship disagrees with the results by Shivaei et al. (2018) who report no log ion evolution with 1500 , a result probably induced by the different dust-attenuation corrections between studies (Balmer Decrements for MOSDEF versus global SED for the rest, see the discussion in Matthee et al. 2017a).
Due to their lower metallicities and burstiness of their SFHs, fainter (low-mass) galaxies are expected to have high ionizing production efficiencies (Ma et al. 2020). If the LyC escape fraction and production efficiency were intertwined (see below), the fact that fainter galaxies produce more ionizing photons per UV luminosity which, at the same time, would escape more easily than in highermass counterparts, supports an early and slow reionization scenario in which the ionizing budget would be dominated by the same lowmass, UV-faint galaxies (see Finkelstein et al. 2019;Trebitsch et al. 2022), thanks to their higher number density at the EoR (Bouwens et al. 2015(Bouwens et al. , 2021. Recently, Prieto-Lyon et al. (2022) used deep HST and JWST multiband photometry to estimate the H flux and thus the ionizing production efficiency extending to very faint UV magnitudes (−23 ≤ 1500 ≤ −15.5) and spanning a wide range in redshift ( = 3−7). Their results demonstrate that, although the general trend of increasing ion with increasing UV magnitude is very smooth, it holds at any 1500 (see also Lam et al. 2019;Emami et al. 2020). Additionally, our thorough compilation of ion measurements (Fig.  12) illustrates the scatter in this relation ( 1dex) and its dependence on the galaxy type, physical properties of the underlying galaxypopulation and methodology itself (with the dust corrections among the different samples playing a critical role in the scatter). While normal high-SFGs (Nakajima et al. 2018b;Lam et al. 2019) have modest production efficiencies, high-LAEs (Matthee et al. 2017b;Harikane et al. 2018;Nakajima et al. 2018a) and C ] emitters (Nakajima et al. 2018b) show much higher log ion values at similar UV magnitudes, since overall these galaxies need stronger ionizing spectra to produce such high excitation lines (see e.g., Reddy et al. 2018).
Intending to qualitatively illustrate this behavior (Fig. 13), we search for meaningful differences in ion between the bulk of the VANDELS LBG population and other sources like LCEs or LAEs, which would potentially boost ion above the canonical value. This is done by splitting the sample into LAEs and non-LAEs, where the the sources whose Ly > −20Å were considered as LAEs, following the definition by Pentericci et al. (2009, but see Stark et al. (2011Kusakabe et al. (2020)). A similar analysis as we did for the LCEs/non-LCEs composites was then performed for the LAEs/non-LAEs composites, and the results are described in the following lines.
Of particular interest is to search for ion differences within different galaxy subsamples. In Fig. 13 (top), the VANDELS sample is color-coded by potential LCEs ( abs esc ≥ 0.05, in blue) and non-LCEs ( abs esc < 0.05, in red), where approximately 60% of the LCEs sample (39 out of 64 LCEs candidates) has log ion ≥ 25.2 Hz/erg, and they uniformly populate the whole range in UV magnitude. As said in Sect. 4.3, the resulting ion for LCEs and non-LCEs are log ion (Hz/erg) ≈ 25.38 and 25.18, respectively. Compared with other values in the literature, our estimated ionizing efficiency for LCEs is consistent with the average of the LCEs sample of the LzLCS . The fact that our results at 3 ≤ ≤ 5 agree with low-LCEs sample supports the argument by which the LzLCS galaxies might be analogs of higher-redshift reionizers (Flury et al. 2022b), This has been lately demonstrated by Mascia et al. (2023a) through JWST spectroscopy of 29 moderately faint galaxies at 5 ≤ ≤ 8, whose optical line ratios, UV slopes,
Similarly, in Fig. 13 (bottom) we explore the ion − 1500 relationship for selected LAEs ( Ly > −20Å, in green) and non-LAEs ( Ly ≤ −20Å, in light-green). Although the 1500 distribution for both galaxy types equally span along the −axis of this figure, around 70% of the LAEs sample (75 out of 104 identified LAEs) falls above the log ion canonical value in the −axis. Our ion estimate for LAEs is in prefect agreement with the results by Harikane et al. (2018) at = 5 − 6, and it matches surprisingly well the global trend found by several high-LAEs surveys at 3 ≤ ≤ 7 (Nakajima et al. 2018a;Harikane et al. 2018;Matthee et al. 2017b, see also Ning et al. (2022)). Similar to Ly , the ionizing efficiency also shapes the emission of other nebular lines ) so that, for instance, higher-than-average ion values have been found among extreme high-[O ] emitters ) and local, high H equivalent width galaxy samples .
The escape of Ly photons and the production of strong UV nebular emission usually requires special ISM and radiation field conditions such as low metalliticy, low dust and neutral gas content, and a vast production of (high-energy) ionizing photons. In Fig. 13 one can visualize a boost in the ionizing production efficiency ( ion ) for high-LAEs' samples. At the same time, a rise in the ionizing photon flux may reduce the covering fraction or column density of neutral hydrogen, which may ease the escape of LyC and Ly photons (see Erb et al. 2014). Recently, Flury et al. (2022b) and Schaerer et al. (2022a) statistically demonstrated that the detection rate of LCEs is enhanced among the LAEs and C emitter galaxy-population (see also Saxena et al. 2022b;Mascia et al. 2023b), reiterating our results.
Moving forward with the discussion, Fig. 14 shows ion as function of the UV-continuum slope for the VANDELS sample when using either the R16 (left) or the SMC (right) attenuation law. Irrespective of the dust law, the ionizing efficiency rapidly decreases with the slope of the continuum, and the trend is more steep for the SMC law. Focusing on the stacks measurements, the canonical limit in ion (log ion ≥ 25.2 Hz/erg) is only achieved for galaxies whose slopes are bluer than 1500 ≤ −1.5. The ion − 1500 relation of a dustless, single-burst stellar population at increasing age is not able to fully reproduce the whole range of 1500 values, so that it would require "some" dust reddening to account for redder slopes (see horizontal arrow in Fig. 14, indicating the change in UV slope after adding 1 mag. of reddening). The nebular continuum would also contribute to increase the spectral slope of the youngest bursts of star formation.
Compared to the already reported relations by Bouwens et al. (2016) at = 4 − 5 and by Matthee et al. (2017a) at = 2 − 3, our results provide a similar slope but systematically shifted to redder 1500 spec . Several effects can contribute to this offset, but we mostly attribute it to the use of broad-band photometry to compute the UV slopes instead of direct spectral measurements, and the different wavelength range probed by multi-band photometry, usually broader than with spectroscopy. In contrast, the flattening of the ion − 1500 relation proposed by Shivaei et al. (2018) (Flury et al. 2022a) and high-LAEs (Matthee et al. 2017b;Harikane et al. 2018;Nakajima et al. 2018a, see Fig. 12 for symbols). The grey-shaded area marks the canonical log ion = 25.2 − 25.3 value given by a simple stellar population at constant SFR over 100Myr (Robertson et al. 2013), and the number of LCEs and LAEs above and below this limit is indicated in the right side of the plot. In conclusion, LCEs and LAEs have systematically higher ion than the bulk of the LBG population (see also Reddy et al. 2022).
for local, compact SFGs, but shifted to lower ion than ours. The new results by Prieto-Lyon et al. (2022, teal-dashed line in the plot) using combined HST+JWST photometry are broadly consistent with our estimations.
In summary, our results provide observational evidence for moderately UV-faint, low-mass and less dusty (un-obscured, blue UV colors) high-redshift galaxies most likely being able to drive the Cosmic Reionization at 6 ≤ ≤ 9 (Finkelstein et al. 2019;Chisholm et al. 2022;Trebitsch et al. 2022;Lin et al. 2023). According to recent JWST observations (e.g., see Endsley et al. 2022;Cullen et al. 2023;Mascia et al. 2023a), these properties seem to be common among the galaxy population at the reionization epoch.
These UV-faint galaxies would produce higher amounts of ionizing photons (through young bursts of star formation, possibly at low metallicities) and allow these photons to escape more easily due to their favoured ISM physical conditions (dustless, gas-less holes in the ISM, see Gazagnes et al. 2020, SL22). Due to their high number density at earlier epochs ( ★ UV ≈ −21 at = 6, see Bouwens et al. 2021), they might become the major contributors to the ionizing budget at the EoR (see Eq. 1).

THE IONIZING PROPERTIES OF GALAXIES: A DETAILED COMPARISON BETWEEN LOW-AND HIGH-REDSHIFT STUDIES
We now compare the ionizing properties of SFGs at 3 ≤ ≤ 5 with the estimated at other other redshifts and examine if and how they could evolve into the EoR. To do so, we put our abs esc and ion trends along with observed properties of galaxies, together with other spectroscopic surveys of LCEs in the literature. On the one hand, we use the LzLCS (Flury et al. 2022a), probing emission line galaxies at 0.3, as our "EoR-analogue sample" of galaxies. The more extreme nature of the LzLCS galaxies, characterized by high ionization parameters, high SFR surface-densities and strong nebular emission lines, makes it the perfect sample to compare with, since EoR galaxies are expected (Schaerer et al. 2016;Boyett et al. 2022) and seem to hold these properties (e.g., Endsley et al. 2022;Schaerer et al. 2022b;Mascia et al. 2023a;Lin et al. 2023). At 3, the KLCS (Steidel et al. 2018;Pahl et al. 2021Pahl et al. , 2023 is used for comparison, because it targets LBGs at similar UV magnitudes (and other properties) to the VANDELS galaxies. For both the LzLCS 9 and the KLCS direct LyC measurements are available. Besides drawing conclusions on the possible redshift evolution of the ionizing properties of galaxies, we will discuss the main systematic uncertainties and caveats of our methodology.
In Fig. 15, the relations between abs esc , ion and the abs esc × ion product with the UV magnitude, UV-continuum slope and Ly strength are compared for the three samples: VANDELS (at = 4, SMC), KLCS (at = 3) and LzLCS (at = 0.3). Particularly important is how the abs esc × ion product compares between the different surveys (third column of the figure).
Interestingly, the properties of the LzLCS and the VANDELS samples show overall fairly consistent trends, despite the redshift range difference and different methods (LzLCS represents direct abs esc measurements, and VANDELS predictions are based on absorption lines). The samples are also quite complementary. Extending over a wide range of UV magnitudes (−22 ≤ 1500 ≤ −18) and stellar masses (10 8 M ≤ ★ ≤ 10 10 M ), the escape fraction of LzLCS galaxies tentatively decreases with the UV brightness and stellar mass. At variance with this, VANDELS galaxies span a much narrower range in terms of UV magnitudes, and it is unable to reveal a clear abs esc correlation with 1500 . Again, the ionizing efficiencies ( ion ) of both low-and high-samples are quite comparable, except for the lack of some low ion objects and the presence of stronger Ly emitters in LzLCS.
LzLCS galaxies show systematically bluer UV slopes (−2.5 ≤ 1500 ≤ −1) than VANDELS galaxies, which translates into much higher abs esc values for the bluest galaxies in the LzLCS (up to 9 We note that the abs esc and ion for LzLCS galaxies were obtained by running our F CUS code in Flury et al. (2022a) and SL22. The abs esc was calculated as the ratio between the intrinsic ionizing flux predicted by F CUS and the observed flux at LyC wavelengths.  (2022) using HST+JWST data at = 3 − 7. The grey-shaded area marks the canonical log ion = 25.2 − 25.3 value given by a simple stellar population at constant SFR over 100Myr (Robertson et al. 2013). log ion decreases with the UV slope almost unanimously for all samples. abs esc ≈ 80%) than any of the galaxies within the VANDELS survey. However, both samples show a significant decrease in their ionizing properties ( abs esc and ion ) with the UV-continuum slope. Strikingly, our VANDELS results follow the same slope provided by the fit to LzLCS data (grey lines in Fig. 15, see Chisholm et al. 2022), but extrapolated to (≥ 0.5) redder UV slopes and lower values of abs esc and ion . This is interesting as it reassures the applicability of our method for inferring the escape fraction of galaxies (see Chisholm et al. (2018), SL22), since the abs esc quoted for the LzLCS is actually the observed escape fraction from measurements of the LyC flux. However, it is worth mentioning the shift in the intercept between the three surveys, with VANDELS and KLCS giving effectively higher abs esc × ion values at a given 1500 than the extrapolation by Chisholm et al. (2022). This is interesting by itself since (1) it might hint a possible redshift evolution of the abs esc × ion − 1500 relation between 0.3 (LzLCS) and 3 (KLCS, VANDELS) or (2), if the high-surveys probe intrinsically lower UV slopes, they may get a better handle on the actual abs esc × ion − 1500 trend in redder galaxies. In any case, the VANDELS abs esc and abs esc × ion relations with 1500 are compatible with Chisholm et al. (2022) within 1 , so none of these hypothesis can be confirmed nor ruled out.
Finally, due to the more extreme nature of the LzLCS objects, their abs esc − Ly and ion − Ly relationships are shifted to higher values of the Ly equivalent width respect to the VANDELS spectra. The intrinsically higher photon production efficiencies ( ion ) in LzLCS galaxies boost the strength of the Ly emission. Regarding the abs esc × ion correlation with Ly for the LzLCS, only LCEs (yellow circles, while downward triangles represent non-LCEs) provide log abs esc × ion (Hz/erg) ≥ 24.2 (see Robertson et al. 2013). Even though our VANDELS survey populate a similar parameters space in terms of 1500 , 1500 , Ly than the KLCS, the VANDELS abs esc values fall systematically below the KLCS points, the latter acting like an upper envelope in all the correlations. In or-der to understand these differences, we now carefully list the main limitations and caveats of our "picket-fence" modeling. There are two main assumptions on the applicability of Eq. 6: (1) the dustattenuation law and (2) the conversion between the metals and H covering fractions (see the full dedicated section in SL22).
First, the preference for a shallower (R16) or a steeper (SMC) attenuation curve can be tested by comparing our SED-derived attenuation versus independent measurements of the UV attenuation of galaxies. For example, one can translate IRX-excess measurements of high-galaxies into UV attenuation at 1600Å ( 1600 , in mag.), by using simple energy balance arguments (see Schaerer et al. 2013). In Fig. 16  Even if the attenuation law was known, abs esc differences may reach depending on the way the UV dust-attenuation ( B−V ) is measured. The dust-attenuation parameter at UV wavelengths is mainly constraint by the slope of the UV continuum (e.g., 1500 , see Chisholm et al. 2022), by comparing the observed 1500 against predictions of the intrinsic int 1500 from SED fitting, either from the spectra or from multi-band photometry. In our case, F CUS  Figure 15. The absolute photon escape fraction ( abs esc ), the ionizing production efficiency (log ion ) and the product of the two (log abs esc × ion , see Eq. 1) for different low-and high-samples. Our predictions for the VANDELS galaxies at 3 ≤ ≤ 5 are shown through dark-blue points in the background (SMC law), and the running locus of the stacking results is indicated via the thick blue line. For comparison, the red open symbols represent the latest KLCS (Steidel et al. 2018) composite results at ∼ 3 by Pahl et al. (2021), and the yellow circles show individual measurements at ∼ 0.3 from LzLCS (Flury et al. 2022a, triangles are 1 upper limits). The dashed grey lines draw the fit relations (median ±1 error) to the LzLCS points . Grey-shaded regions mark the canonical values of abs esc ≥ 5% and log ion = 25.2 − 25.3 in classical Reionization models (Robertson et al. 2013).
is able to estimate B−V based on the spectral shape of the UV continuum. However, the slope of the current spectral observations might be subject to uncertainty due to different factors (see Garilli et al. 2021): instrumental calibrations, reduction pipeline issues, flux losses because of the atmospheric dispersion and sky subtraction residuals. An alternative solution would be to use an independent determination of the UV slope. With this purpose, we use the 1500 measurements from Calabrò et al. (2021), taking all the photometric bands whose bandwidths lie entirely inside the 1230-2750Å restframe wavelength range, and the relation given in Chisholm et al. (2022) (Eq. 8, assuming int 1500 = −2.5) to convert this 1500 to B−V . Plugging these new dust-attenuation values into Eq. 6 pro-duces higher estimations of abs esc . The systematic deviation between the spectroscopic (F CUS) and photometric (Calabrò et al. 2021) UV slopes by 0.5 on average, corresponds to 1 mag difference in UV attenuation at LyC wavelengths ( 912 ), which leads to a ×2.5 higher photon escape fraction. As sketched in Fig. 17, the new abs esc values (pink pentagons) are in closer agreement with the measurements by Pahl et al. (2021) and Begley et al. (2022), and fully consistent with the physical scenario in which the escape fraction of galaxies decreases towards brighter systems.
Second, in Fig. 18 we study the influence that the adopted conversion between the metals and H covering fraction has on the predicted escape fractions. For comparison, we compute the NO resolution correction Figure 17. A representation of the potential systematic bias in abs esc due to the different values in the UV dust-attenuation that can be derived from dissimilar measurements of the −slope. The blue diamonds show our default escape fractions from F CUS, i.e., B−V is derived from spectroscopic measurements of the UV continuum slope ( spec ). Conversely, the pink pentagons display the abs esc values when B−V is inferred from the −slope measurements by Calabrò et al. (2021), using photometry only ( phot ). The average abs esc = 0.07 ± 0.02 reported by Begley et al. (2022) at 3.5 is indicated through the grey shaded-area, and the results by Pahl et al. (2021) are plotted through black squares. The effect of spectral resolution due to changes in the depths of the absorption lines is represented with the dashed line. predicted abs esc values for the KLCS sample adopting the "Screen Model" (see Table 9 in Steidel et al. 2018), which is identical to our default "picket-fence" geometry assuming a uniform slab of dust. We are not able to reconcile our escape fractions with the KLCS points, with our predictions lying a factor of ×2 below them. However, if one assumes a 1:1 metals-to-H covering fraction correspondence so that (LIS) = (HI), instead of the usual linear regression by which (HI > (LIS)) (see Gazagnes et al. (2018) and SL22), the resulting abs esc would agree with the estimate of the average abs esc by Begley et al. (2022, grey shaded area), which is now much more compatible with KLCS at the same time.
Finally, the resolution correction function applied to the residual flux of the LIS lines (App. A) might also affect our escape fraction estimates. In Fig. 17 and Fig. 18, we explicitly show the effect of our custom correction in the inferred abs esc (dashed blue line). In the more conservative case in which the spectral resolution is totally dismissed and non-considered in Eq. 6, abs esc would increase by a factor of ×1 − 1.5 at most: this is by far the effect with the least influence on the average abs esc among the ones considered in this manuscript.
In conclusion, different effects can potentially augment our average abs esc . Among them, the determination of the UV dustattenuation has the strongest impact in the escape fraction estimates. In App. C (Fig. C1), we plot the abs esc values resulting from plugging into Eq. 6 the dust-attenuation values derived from the 1500 measurements from Calabrò et al. (2021). As we can see, all VAN-DELS, KLCS and LzLCS abs esc correlations with physical quantities now converge to the same global trends, and they share a common behaviour by which abs esc increases towards fainter UV magnitudes, bluer UV slopes and stronger Ly emission. In view of these results, the use of 1500 from photometry as a proxy of the UV dust attenuation in combination with the depth of the absorption lines as a tracer of H -empty channels in the ISM of galaxies (Eq. 6), seem to be a robust way of predicting the escape fraction of galaxies at any redshift. Independently, there is also the chance of (1) adopting a very steep dust-attenuation law, whose extrapolation towards LyC wavelengths would be even higher than when using an SMClike law and (2) using a favourable metals-to-H covering fraction's conversion, in which the metals would spatially trace the neutral gas more closely than with the current assumed relation.
The remaining question is whether a universal relation between the ionizing properties of galaxies ( abs esc and ion ) and their physical properties ( 1500 , 1500 , Ly etc.) exists, in which the different surveys would populate different regions of the hypothetical universal trend or, contrarily, the observed differences between surveys (VANDELS, KLCS, LzLCS) indicates an evolution of these properties with redshift. Still, as far as our sample is concerned, the abs esc × ion dependency along with galaxy properties is mainly due to escape fraction variations, so that the redshift evolution of this product will ultimately follow the changes in the ISM geometry, gas and dust composition across cosmic time.
In any case, and regardless the caveats that we have just described, the LzLCS, KLCS and VANDELS results suggest a scenario in which moderately UV-faint, low-mass and dustless galaxies most likely dominate the SFG ionizing photon budget at any epoch. The stellar populations of EoR-galaxies would be characterized by strong radiation fields, with high ionizing production efficiencies (boosting the nebular emission), and whose ISM conditions would be such that they facilitate the escape of a copious amount of LyC photons to the IGM. This can be speculatively attributed to the stellar mass build-up and changes in the gas and dust properties of galaxies between the Cosmic Noon ( 3) and the EoR (6 ≤ ≤ 9). NO resolution correction C f (HI) = C f (LIS) Figure 18. A representation of the potential systematic bias in abs esc (blue diamonds) due to the conversion between the metals and neutral gas covering fractions (blue dotted line). For comparison we also show the predicted abs esc values for the KLCS stacks (red symbols), by applying our methodology to the published covering fractions and UV attenuations from Steidel et al. (2018). The layout is the same as in Fig. 17.

SUMMARY AND CONCLUSIONS
In this work, we fully exploit and highlight the ability of absorption line and rest-UV spectroscopy alone to decipher the ionising properties of high-z SFGs. Our novel approach make use of deep, rest-frame UV spectra from the VANDELS survey at 3 ≤ ≤ 5, to compute the ionizing absolute photon escape fraction ( abs esc ) and ionizing photon production efficiency ( ion ) of galaxies. abs esc has been derived by combining absorption line measurements with estimates of the UV attenuation (see Saldana-Lopez et al. 2022), while the ion parameter was computed by fitting the FUV stellar continuum of the VANDELS galaxies (following Chisholm et al. 2019).
In particular, we have searched for correlations between abs esc and ion along with different galaxy properties (UV magnitude, UV slope, Ly strength, etc.), and we thoroughly compared with independent literature estimates. We found that: • The predicted abs esc monotonically decreases with the stellar mass, the UV-continuum slope and the Ly equivalent width of the VANDELS galaxies (Fig. 5, 6, 7). We find a non-significant correlation between abs esc and the UV magnitude, although the faintest galaxies tentatively have higher escape fractions.
• The estimated ion statistically increases towards blue UVcontinuum slopes and strong Ly emitting galaxies (Fig. 12, 14), and it smoothly raises beyond the canonical value towards the UVfaintest galaxies in the sample.
Additionally, we constructed the average composite FUV spec-trum of LCEs at 3 ≤ ≤ 5 (Fig. 8, 10, 11), by stacking potential, individual emitters in the VANDELS survey (selected based on our predicted abs esc ≥ 5%), and explained their non-ionizing spectral properties in the framework of the ISM conditions which enable ionizing photons to escape. Our results show that the FUV spectra of typical high-LCEs would be characterized by: • Blue UV slopes ( 1500 spec −2) compared to non-LCEs, which directly translates into low UV attenuation ( B−V 0.1 mag.) and therefore low column density of dust along the line-of-sight.
• Enhanced Ly emission ( Ly −25Å) and strong UV nebular lines in contrast to the non-LCEs population, particularly high C 1908/C 1550 ratios ( 0.75). Together with the intrinsically higher ion parameter for LCEs, this indicates very young underlying stellar populations (≈ 10 Myr) at relatively low metallicities (≈ 0.2 Z ).
• Weak ( 1Å) ISM LIS absorption lines (e.g., Si 1260, C 1334), while the HIS absorption lines (e.g., Si 1400) are of similar strength as the bulk of non-LCEs. This, together with the low UV attenuation for LCEs, suggest the presence of dustless, low gas column-density channels in the ISM which favour the escape of ionzing photons.
Finally, we have compared our findings with other LyC surveys in the literature (Fig. 15), concretely the Keck Lyman Continuum Survey (or KLCS, Steidel et al. 2018;Pahl et al. 2021Pahl et al. , 2023, targeting 3 LBGs at similar UV magnitudes to VANDELS, and the Low-Redshift Lyman Continuum Survey (or LzlCS, Flury et al. 2022a,b), which targetted emission line galaxies at 0.3. VAN-DELS and LzLCS show overall fairly consistent trends, with LzLCS shifted to fainter UV magnitudes, bluer UV slopes and stronger Ly emission, and therefore their abs esc and ion are enhanced with respect to VANDELS galaxies. The escape fraction of KLCS galaxies fall above our estimates at all UV magnitudes. We mainly ascribed this discrepancy to the way the amount of UV attenuation is measured, and we propose to use 1500 measurements from photometry as an independent proxy of the UV dust attenuation of galaxies (Fig.  17). The dust-attenuation law (Fig. 16) and the metals-to-neutral-gas covering fraction conversion (Fig. 18) constitutes additional sources of uncertainty to the escape fraction.
Our joint analysis of the VANDELS, LzLCS and KLCS results shed light onto the fact that UV-faint, low-mass and dustless galaxies likely dominated the ionizing budget during the EoR. Their stellar populations would most likely be characterized by strong radiation fields, with high ionizing production efficiencies, and whose ISM conditions are such that they favour the escape of LyC photons. Recent JWST observations (e.g., Endsley et al. 2022;Schaerer et al. 2022b;Cullen et al. 2023;Lin et al. 2023;Mascia et al. 2023a), have revealed that blue UV slopes and strong emission lines also characterized the less massive and moderately faint galaxies at the EoR, supporting our results.
An increasing number of high-quality, FUV spectra at > 6 will be available in the mid-future thanks to upcoming ground-based (ELT, GMT, TMT, ...) and currently ongoing space-based facilities (JWST). This will give us the first insights into the properties of galaxies during the EoR. Nevertheless, according to our work, additional efforts are still needed in order to correctly decipher the physics underlying such vast data-sets. The two main questions that we think are urgent and clearly needed for future studies are: (1) if unique, what is the dust-attenuation law governing low-mass, low-metallicity galaxies? And (2) do metals trace the same spatial location as the neutral gas within the ISM of these galaxies?   Figure C1. The ionizing photon escape fraction ( abs esc ) versus the UV magnitude ( 1500 ), the UV slope at 1500Å ( 1500 ) and the Ly equivalent width ( Ly ). Blue dots represent the predicted abs esc values for VANDELS galaxies (Eq. 6, using SMC) when the dust-attenuation ( B−V ) is inferred from the 1500 measurement by Calabrò et al. (2021). Layout is the same as Fig. 15.