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

A. Saldana-Lopez, D. Schaerer, J. Chisholm, A. Calabrò, L. Pentericci, A.

Saxena, R. Amorín, A.C. Carnall, F. Fontanot, J.P.U. Fynbo, et al.

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

## ► To cite this version:

A. Saldana-Lopez, D. Schaerer, J. Chisholm, A. Calabrò, L. Pentericci, et al.. The VANDELS survey: the ionizing properties of star-forming galaxies at $3 \le z \le 5$ using deep rest-frame ultraviolet spectroscopy. Monthly Notices of the Royal Astronomical Society, 2023, 522 (4), pp.6295-6325. 10.1093/mnras/stad1283. hal-03862325

HAL Id: hal-03862325

https://hal.science/hal-03862325

Submitted on 16 Feb 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The VANDELS survey: the ionizing properties of star-forming galaxies at $3 \le z \le 5$ using deep rest-frame ultraviolet spectroscopy

A. Saldana-Lopez, D. Schaerer, J. Chisholm, A. Calabrò, L. Pentericci, A.

Saxena, R. Amorín, A.C. Carnall, F. Fontanot, J.P.U. Fynbo, et al.

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

## ► To cite this version:

A. Saldana-Lopez, D. Schaerer, J. Chisholm, A. Calabrò, L. Pentericci, et al.. The VANDELS survey: the ionizing properties of star-forming galaxies at $3 \le z \le 5$ using deep rest-frame ultraviolet spectroscopy. Monthly Notices of the Royal Astronomical Society, 2023, 522 (4), pp.6295-6325. 10.1093/mnras/stad1283. hal-03862325

HAL Id: hal-03862325

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

Submitted on 16 Feb 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The VANDELS survey: the ionizing properties of star-forming galaxies at $3 \le z \le 5$ using deep rest-frame ultraviolet spectroscopy

A. Saldana-Lopez<sup>1,★</sup> D. Schaerer,<sup>1,2</sup> J. Chisholm<sup>1,3</sup> A. Calabrò,<sup>4</sup> L. Pentericci,<sup>4</sup> F. Cullen<sup>1,5</sup>

A. Saxena,<sup>6,7</sup> R. Amorín<sup>1,8</sup> A. C. Carnall<sup>1,5</sup> F. Fontanot<sup>1,9</sup> J. P. U. Fynbo,<sup>10,11</sup> L. Guaita,<sup>12</sup>

N. P. Hathi<sup>1,13</sup> P. Hibon<sup>1,14</sup> Z. Ji,<sup>15</sup> D. J. McLeod,<sup>5</sup> E. Pompei<sup>14</sup> and G. Zamorani<sup>16</sup>

Accepted 2023 April 24. Received 2023 April 24; in original form 2022 November 21

## ABSTRACT

The physical properties of Epoch of Reionization (EoR) galaxies are still poorly constrained by observations. To better understand the ionizing properties of galaxies in the EoR, we investigate deep, rest-frame ultraviolet (UV) spectra of ≈500 star-forming galaxies at $3 \le z \le 5$ selected from the public ESO-VANDELS spectroscopic survey. The absolute ionizing photon escape fraction ($f_{esc}^{abs}$, i.e. the ratio of *leaking* against *produced* ionizing photons) is derived by combining absorption line measurements with estimates of the UV attenuation. The ionizing production efficiency ($\xi_{ion}$, i.e. the number of ionizing photons produced per non-ionizing UV luminosity) is calculated by fitting the far-UV (FUV) stellar continuum of the VANDELS galaxies. We find that the $f_{esc}^{abs}$ and $\xi_{ion}$ parameters increase towards low-mass, blue UV-continuum slopes and strong Ly$\alpha$ emitting galaxies, and both are slightly higher-than-average for the UV-faintest galaxies in the sample. Potential Lyman Continuum Emitters (LCEs, $f_{esc}^{abs} \ge 5$ per cent) and selected Lyman Alpha Emitters (LAEs, $W_{Ly\alpha} \le -20$ Å) show systematically higher $\xi_{ion}$ (log $\xi_{ion}$(Hz erg$^{-1}$) ≈ 25.38, 25.41) than non-LCEs and non-LAEs (log $\xi_{ion}$(Hz erg$^{-1}$) ≈ 25.18, 25.14) at similar UV magnitudes. This indicates very young underlying stellar populations (≈10 Myr) at relatively low metallicities (≈0.2 Z$_{⊙}$). The FUV non-ionizing spectra of potential LCEs is characterized by blue UV slopes (≤−2), enhanced Ly$\alpha$ emission (≤−25 Å), strong UV nebular lines (e.g. high C IV 1550/C III 1908 ≥0.75 ratios), and weak absorption lines (≤1 Å). 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.

**Key words:** dust, extinction – galaxies: high-redshift – galaxies: ISM – galaxies: stellar content – dark ages, reionization, first stars – ultraviolet: galaxies.

## 1 INTRODUCTION

Several data sets (see Goto et al. 2021, and references therein) measuring the redshift ($z$) evolution of the volume-averaged neutral hydrogen fraction provide evidence for the last major phase

change underwent by the Universe, the Cosmic Reionization. Between $z=9$ and $z=6$ (Planck Collaboration XIII 2016), the number of ionizing photons emitted per unit time ($\dot{N}_{ion}$) overcame the recombination rate ($\Gamma_{H1}$) of hydrogen atoms, $\dot{N}_{ion} \ge \Gamma_{H1}$ (Madau, Haardt & Rees 1999), so that the bulk of H1 neutral gas within the Intergalactic Medium (IGM) *progressively* transitioned to an ionized state (see Dayal & Ferrara 2018, for a review).

<footnote n="1"><sup>1</sup>Department of Astronomy, University of Geneva, 51 Chemin Pegasi, CH-1290 Versoix, Switzerland</footnote>

<footnote n="2"><sup>2</sup>CNRS, IRAP, 14 Avenue E. Belin, F-31400 Toulouse, France</footnote>

<footnote n="3"><sup>3</sup>Department of Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712-1205, USA</footnote>

<footnote n="4"><sup>4</sup>INAF – Osservatorio Astronomico di Roma, via Frascati 33, I-00078 Monteporzio Catone, Italy</footnote>

<footnote n="5"><sup>5</sup>Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="6"><sup>6</sup>Sub-department of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="7"><sup>7</sup>Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="8"><sup>8</sup>Instituto de Investigación Multidisciplinar en Ciencia y Tecnología, Universidad de La Serena, Raúl Bitrán, 1305 La Serena, Chile</footnote>

<footnote n="9"><sup>9</sup>INAF – Trieste Observatory, via GB Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="10"><sup>10</sup>Cosmic Dawn Center (DAWN), Jagtvej 128, DK-2200 Copenhagen N, Denmark</footnote>

<footnote n="11"><sup>11</sup>Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, DK-2100 Copenhagen Ø, Denmark</footnote>

<footnote n="12"><sup>12</sup>Departamento de Ciencias Físicas, Facultad de Ciencias Exactas, Universidad Andres Bello, Fernandez Concha 700, Las Condes, Santiago, Chile</footnote>

<footnote n="13"><sup>13</sup>Space Telescope Science Institute, Baltimore, MD 21218, USA</footnote>

<footnote n="14"><sup>14</sup>European Southern Observatory, Alonso de Córdova 3107, Vitacura, Santiago de Chile, Chile</footnote>

<footnote n="15"><sup>15</sup>Department of Astronomy, University of Massachusetts, Amherst, 710 North Pleasant Street, Amherst, MA 01003-9305, USA</footnote>

<footnote n="16"><sup>16</sup>INAF – Astrophysics and Space Science Observatory, Via Piero Gobetti 93/3, I-40129 Bologna, Italy</footnote>

<footnote>*E-mail: alberto.saldanalopez@unige.ch</footnote>
# The VANDELS survey: the ionizing properties of star-forming galaxies at $3 \le z \le 5$ using deep rest-frame ultraviolet spectroscopy

A. Saldana-Lopez<sup>1,★</sup> D. Schaerer,<sup>1,2</sup> J. Chisholm<sup>1,3</sup> A. Calabrò,<sup>4</sup> L. Pentericci,<sup>4</sup> F. Cullen<sup>1,5</sup>

A. Saxena,<sup>6,7</sup> R. Amorín<sup>1,8</sup> A. C. Carnall<sup>1,5</sup> F. Fontanot<sup>1,9</sup> J. P. U. Fynbo,<sup>10,11</sup> L. Guaita,<sup>12</sup>

N. P. Hathi<sup>1,13</sup> P. Hibon<sup>1,14</sup> Z. Ji,<sup>15</sup> D. J. McLeod,<sup>5</sup> E. Pompei<sup>14</sup> and G. Zamorani<sup>16</sup>

Accepted 2023 April 24. Received 2023 April 24; in original form 2022 November 21

## ABSTRACT

The physical properties of Epoch of Reionization (EoR) galaxies are still poorly constrained by observations. To better understand the ionizing properties of galaxies in the EoR, we investigate deep, rest-frame ultraviolet (UV) spectra of ≈500 star-forming galaxies at $3 \le z \le 5$ selected from the public ESO-VANDELS spectroscopic survey. The absolute ionizing photon escape fraction ($f_{esc}^{abs}$, i.e. the ratio of *leaking* against *produced* ionizing photons) is derived by combining absorption line measurements with estimates of the UV attenuation. The ionizing production efficiency ($\xi_{ion}$, i.e. the number of ionizing photons produced per non-ionizing UV luminosity) is calculated by fitting the far-UV (FUV) stellar continuum of the VANDELS galaxies. We find that the $f_{esc}^{abs}$ and $\xi_{ion}$ parameters increase towards low-mass, blue UV-continuum slopes and strong Ly$\alpha$ emitting galaxies, and both are slightly higher-than-average for the UV-faintest galaxies in the sample. Potential Lyman Continuum Emitters (LCEs, $f_{esc}^{abs} \ge 5$ per cent) and selected Lyman Alpha Emitters (LAEs, $W_{Ly\alpha} \le -20$ Å) show systematically higher $\xi_{ion}$ (log $\xi_{ion}$ (Hz erg$^{-1}$) ≈ 25.38, 25.41) than non-LCEs and non-LAEs (log $\xi_{ion}$ (Hz erg$^{-1}$) ≈ 25.18, 25.14) at similar UV magnitudes. This indicates very young underlying stellar populations (≈10 Myr) at relatively low metallicities (≈0.2 Z$_{⊙}$). The FUV non-ionizing spectra of potential LCEs is characterized by blue UV slopes (≤-2), enhanced Ly$\alpha$ emission (≤-25 Å), strong UV nebular lines (e.g. high C IV 1550/C III 1908 ≥0.75 ratios), and weak absorption lines (≤1 Å). 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.

**Key words:** dust, extinction – galaxies: high-redshift – galaxies: ISM – galaxies: stellar content – dark ages, reionization, first stars – ultraviolet: galaxies.

## 1 INTRODUCTION

Several data sets (see Goto et al. 2021, and references therein) measuring the redshift ($z$) evolution of the volume-averaged neutral hydrogen fraction provide evidence for the last major phase

change underwent by the Universe, the Cosmic Reionization. Between $z=9$ and $z=6$ (Planck Collaboration XIII 2016), the number of ionizing photons emitted per unit time ($\dot{N}_{ion}$) overcame the recombination rate ($\Gamma_{H1}$) of hydrogen atoms, $\dot{N}_{ion} \ge \Gamma_{H1}$ (Madau, Haardt & Rees 1999), so that the bulk of H1 neutral gas within the Intergalactic Medium (IGM) *progressively* transitioned to an ionized state (see Dayal & Ferrara 2018, for a review).

<footnote n="1"><sup>1</sup>Department of Astronomy, University of Geneva, 51 Chemin Pegasi, CH-1290 Versoix, Switzerland</footnote>

<footnote n="2"><sup>2</sup>CNRS, IRAP, 14 Avenue E. Belin, F-31400 Toulouse, France</footnote>

<footnote n="3"><sup>3</sup>Department of Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712-1205, USA</footnote>

<footnote n="4"><sup>4</sup>INAF – Osservatorio Astronomico di Roma, via Frascati 33, I-00078 Monteporzio Catone, Italy</footnote>

<footnote n="5"><sup>5</sup>Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="6"><sup>6</sup>Sub-department of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="7"><sup>7</sup>Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="8"><sup>8</sup>Instituto de Investigación Multidisciplinar en Ciencia y Tecnología, Universidad de La Serena, Raúl Bitrán, 1305 La Serena, Chile</footnote>

<footnote n="9"><sup>9</sup>INAF – Trieste Observatory, via GB Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="10"><sup>10</sup>Cosmic Dawn Center (DAWN), Jagtvej 128, DK-2200 Copenhagen N, Denmark</footnote>

<footnote n="11"><sup>11</sup>Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, DK-2100 Copenhagen Ø, Denmark</footnote>

<footnote n="12"><sup>12</sup>Departamento de Ciencias Físicas, Facultad de Ciencias Exactas, Universidad Andres Bello, Fernandez Concha 700, Las Condes, Santiago, Chile</footnote>

<footnote n="13"><sup>13</sup>Space Telescope Science Institute, Baltimore, MD 21218, USA</footnote>

<footnote n="14"><sup>14</sup>European Southern Observatory, Alonso de Córdova 3107, Vitacura, Santiago de Chile, Chile</footnote>

<footnote n="15"><sup>15</sup>Department of Astronomy, University of Massachusetts, Amherst, 710 North Pleasant Street, Amherst, MA 01003-9305, USA</footnote>

<footnote n="16"><sup>16</sup>INAF – Astrophysics and Space Science Observatory, Via Piero Gobetti 93/3, I-40129 Bologna, Italy</footnote>

<footnote>*E-mail: alberto.saldanalopez@unige.ch</footnote>
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 (AGNs) to the ionizing budget during the Epoch of Reionization (EoR; Hassan et al. 2018; Kulkarni, Worseck & Hennawi 2019; Dayal et al. 2020), mainly because of their lower number density at early epochs (e.g. Matsuoka et al. 2018, but see Fontanot, Cristiani & Vanzella (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 $z < 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, 2015) dominated the ionizing emissivity at the EoR.

On the one hand, the remarkable modelling 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$\alpha$ emitters (LAEs) over time, support a late-completed and rapid reionization purely dominated by more massive and moderately luminous sub-$M_{\text{UV}}^*$ systems. These results are compatible with the rapid reionization modelled by Mason et al. (2019), in which the ionizing emissivity was constrained from CMB optical depth and Ly$\alpha$ 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 \le z \le 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 ($\dot{n}_{\text{ion}}$), i.e. the number of ionizing photons emitted per unit time and comoving volume (Robertson 2022):

$$ \dot{n}_{\text{ion}} = f_{\text{esc}}^{\text{abs}} \xi_{\text{ion}} \rho_{\text{UV}}, \quad (1) $$

where $f_{\text{esc}}^{\text{abs}}$ 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 $\xi_{\text{ion}}$ is the so-called ionizing photon production efficiency (ionizing photons generated per non-ionizing intrinsic UV luminosity). $\rho_{\text{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 $\rho_{\text{UV}}$, equivalently) of galaxies is relatively well-constrained 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 $z = 8$ to 10 (Oesch et al. 2018; 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, McLure & Dunlop 2016; Livermore, Finkelstein & Lotz 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’.

Our paradigm of the early Universe is rapidly changing thanks to the JWST. The JWST early release science and observations have just made possible the photometric selection (e.g. Atek et al. 2023; Furtak et al. 2023; Labbé 2023) and spectroscopic characterization (Arellano-Córdova et al. 2022; Brinchmann 2022; Schaerer et al. 2022b; Tacchella et al. 2022; Trussler et al. 2022; Carnall et al. 2023; Curti et al. 2023; Trump et al. 2023) of galaxies at the EoR. Joint efforts combining both photometry and spectroscopy of some of the first JWST programmes (see e.g. Harikane et al. 2023; Isobe et al. 2023; Nakajima et al. 2023) have shown as the best approach so far to study the properties of EoR galaxies in the context of galaxy evolution, with surveys such as GLASS (Castellano et al. 2022a, b; Mascia et al. 2023a; Nanayakkara et al. 2023; Santini et al. 2023), UNCOVER (Bezanson et al. 2022; Weaver et al. 2023), CEERS (Endsley et al. 2022; Topping et al. 2022; Finkelstein et al. 2023; Fujimoto et al. 2023; Tang et al. 2023), or JADES (Cameron et al. 2023; Saxena et al. 2023; Curtis-Lake et al. 2023; Robertson et al. 2023).

The ionizing production efficiency of galaxies ($\xi_{\text{ion}}$) is more uncertain. Surveys targeting Ly$\alpha$ 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. 2023), although with a huge scatter within which $\xi_{\text{ion}}$ spans a wide range of values depending on the galaxy type (log $\xi_{\text{ion}}(\text{Hz erg}^{-1}) = 24-26$). Reassuringly, the overall $\xi_{\text{ion}}$ evolution with galaxy properties at these redshifts is in line with the usual formalism by which a log $\xi_{\text{ion}}(\text{Hz erg}^{-1}) \ge 25.2$ SFG population with constant SFH (over 100 Myr) is able to fully reionize the Universe, assuming a fixed escape fraction of $f_{\text{esc}}^{\text{abs}} \ge 5$ per cent (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).

Finally, the LyC absolute escape fraction ($f_{\text{esc}}^{\text{abs}}$) is by far the most unknown parameter in equation (1) (Faisst 2016). Starting a decade ago, the search for LyC emitters (LCEs), targeting Lyman Break Galaxies (LBGs) at $z = 0.7-4$ through expensive imaging (Vanzella et al. 2010, 2012, 2015; Mostardi et al. 2015; de Barros et al. 2016; Grazian et al. 2016, 2017; Japelj et al. 2017; Micheva et al. 2017; Rutkowski et al. 2017; Naidu et al. 2018; Alavi et al. 2020; Bian & Fan 2020; Begley et al. 2022) and spectroscopic campaigns (Steidel, Pettini & Adelberger 2001; Shapley et al. 2006; Bridge et al. 2010; Marchi et al. 2017; Meštrić et al. 2021; Prichard et al. 2022) remained unsuccessful, where most of the estimates for global escape fraction of the SF galaxy-population relied on stacking and $f_{\text{esc}}^{\text{abs}}$ upper limits. The systematic searches at $z \simeq 3$ by the Lyman Continuum Escape Survey (LACES; Fletcher et al. 2019; Nakajima et al. 2020) and other surveys by Saxena et al. (2022b) and Rivera-Thorsen, Hayes & Melinder (2022) have provided the first statistical correlations between individual escape fraction measurements and diverse galaxy properties. Other remarkable individual detections are Vanzella et al. (2012, 2015, 2018), de Barros et al. (2016), Shapley et al. (2016), Saha et al. (2020), Marques-Chaves et al. (2021, 2022a).

However, our knowledge about LCEs and their physical properties is mainly due to the *Keck Lyman Continuum Survey* (KLCS; Steidel
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 (AGNs) to the ionizing budget during the Epoch of Reionization (EoR; Hassan et al. 2018; Kulkarni, Worseck & Hennawi 2019; Dayal et al. 2020), mainly because of their lower number density at early epochs (e.g. Matsuoka et al. 2018, but see Fontanot, Cristiani & Vanzella (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 $z < 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, 2015) dominated the ionizing emissivity at the EoR.

On the one hand, the remarkable modelling 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$\alpha$ emitters (LAEs) over time, support a late-completed and rapid reionization purely dominated by more massive and moderately luminous sub-$M_{\text{UV}}^*$ systems. These results are compatible with the rapid reionization modelled by Mason et al. (2019), in which the ionizing emissivity was constrained from CMB optical depth and Ly$\alpha$ 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 \le z \le 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 ($\dot{n}_{\text{ion}}$), i.e. the number of ionizing photons emitted per unit time and comoving volume (Robertson 2022):

$$ \dot{n}_{\text{ion}} = f_{\text{esc}}^{\text{abs}} \xi_{\text{ion}} \rho_{\text{UV}}, \quad (1) $$

where $f_{\text{esc}}^{\text{abs}}$ 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 $\xi_{\text{ion}}$ is the so-called ionizing photon production efficiency (ionizing photons generated per non-ionizing intrinsic UV luminosity). $\rho_{\text{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 $\rho_{\text{UV}}$, equivalently) of galaxies is relatively well-constrained 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 $z = 8$ to 10 (Oesch et al. 2018; 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, McLure & Dunlop 2016; Livermore, Finkelstein & Lotz 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’.

Our paradigm of the early Universe is rapidly changing thanks to the JWST. The JWST early release science and observations have just made possible the photometric selection (e.g. Atek et al. 2023; Furtak et al. 2023; Labbé 2023) and spectroscopic characterization (Arellano-Córdova et al. 2022; Brinchmann 2022; Schaerer et al. 2022b; Tacchella et al. 2022; Trussler et al. 2022; Carnall et al. 2023; Curti et al. 2023; Trump et al. 2023) of galaxies at the EoR. Joint efforts combining both photometry and spectroscopy of some of the first JWST programmes (see e.g. Harikane et al. 2023; Isobe et al. 2023; Nakajima et al. 2023) have shown as the best approach so far to study the properties of EoR galaxies in the context of galaxy evolution, with surveys such as GLASS (Castellano et al. 2022a, b; Mascia et al. 2023a; Nanayakkara et al. 2023; Santini et al. 2023), UNCOVER (Bezanson et al. 2022; Weaver et al. 2023), CEERS (Endsley et al. 2022; Topping et al. 2022; Finkelstein et al. 2023; Fujimoto et al. 2023; Tang et al. 2023), or JADES (Cameron et al. 2023; Saxena et al. 2023; Curtis-Lake et al. 2023; Robertson et al. 2023).

The ionizing production efficiency of galaxies ($\xi_{\text{ion}}$) is more uncertain. Surveys targeting Ly$\alpha$ 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. 2023), although with a huge scatter within which $\xi_{\text{ion}}$ spans a wide range of values depending on the galaxy type (log $\xi_{\text{ion}}$(Hz erg$^{-1}$) = 24–26). Reassuringly, the overall $\xi_{\text{ion}}$ evolution with galaxy properties at these redshifts is in line with the usual formalism by which a log $\xi_{\text{ion}}$(Hz erg$^{-1}$) $\ge 25.2$ SFG population with constant SFH (over 100 Myr) is able to fully reionize the Universe, assuming a fixed escape fraction of $f_{\text{esc}}^{\text{abs}} \ge 5$ per cent (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).

Finally, the LyC absolute escape fraction ($f_{\text{esc}}^{\text{abs}}$) is by far the most unknown parameter in equation (1) (Faisst 2016). Starting a decade ago, the search for LyC emitters (LCEs), targeting Lyman Break Galaxies (LBGs) at $z = 0.7–4$ through expensive imaging (Vanzella et al. 2010, 2012, 2015; Mostardi et al. 2015; de Barros et al. 2016; Grazian et al. 2016, 2017; Japelj et al. 2017; Micheva et al. 2017; Rutkowski et al. 2017; Naidu et al. 2018; Alavi et al. 2020; Bian & Fan 2020; Begley et al. 2022) and spectroscopic campaigns (Steidel, Pettini & Adelberger 2001; Shapley et al. 2006; Bridge et al. 2010; Marchi et al. 2017; Meštrić et al. 2021; Prichard et al. 2022) remained unsuccessful, where most of the estimates for global escape fraction of the SF galaxy-population relied on stacking and $f_{\text{esc}}^{\text{abs}}$ upper limits. The systematic searches at $z \simeq 3$ by the Lyman Continuum Escape Survey (LACES; Fletcher et al. 2019; Nakajima et al. 2020) and other surveys by Saxena et al. (2022b) and Rivera-Thorsen, Hayes & Melinder (2022) have provided the first statistical correlations between individual escape fraction measurements and diverse galaxy properties. Other remarkable individual detections are Vanzella et al. (2012, 2015, 2018), de Barros et al. (2016), Shapley et al. (2016), Saha et al. (2020), Marques-Chaves et al. (2021, 2022a).

However, our knowledge about LCEs and their physical properties is mainly due to the *Keck Lyman Continuum Survey* (KLCS; Steidel
et al. 2018; Pahl et al. 2021, 2023) at $z \simeq 3$, and the pioneering work by Izotov et al. (2016b, a, 2018a, b, 2021, 2022) and the recent Low-Redshift Lyman Continuum Survey (LzLCS; Flury et al. 2022a) at $z \simeq 0.3$. Particularly, compactness, high ionization parameters (traced by optical line ratios), high star-formation rate (SFR) surface densities, strong Ly $\alpha$ emission, and low dust-attenuation (traced by the UV slope) seem to characterize the strongest LCEs at low-$z$ (see also Wang et al. 2019; Izotov et al. 2020), with ionizing photon production efficiencies analogued to $z \simeq 6$ galaxies (Schaerer et al. 2016; Chisholm et al. 2022). Interestingly, and according to the analysis of new JWST observations by Mascia et al. (2023a), these characteristics resemble the properties of galaxies at the EoR (see also Endsley et al. 2022; Schaerer et al. 2022b; Cullen et al. 2023; Lin et al. 2023).

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 $f_{\text{esc}}^{\text{abs}}$ 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 $\alpha$ 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 II 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-$z$ SFGs and their evolution along with the different galaxy properties. For that, we make use of a sample of $\simeq 500$ deep, rest-frame ultraviolet spectra at $3 \le z \le 5$ drawn from the VANDELS survey (McLure et al. 2018b; Pentericci et al. 2018; Garilli et al. 2021). In particular, LyC absolute escape fractions ($f_{\text{esc}}^{\text{abs}}$) are derived by measuring the depth of the absorption lines in combination with the UV attenuation, whilst ionizing production efficiencies ($\xi_{\text{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 Section 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 Section 3. The main results of this research, looking for correlations between $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$ with the different galaxy properties are summarized in Section 4. Our results on the ionizing efficiency of high-$z$ galaxies are compared with different estimates in the literature in Section 5, and finally the possible redshift evolution of the $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ product is discussed in Section 6, by comparing our values against state-of-the-art low- and high-$z$ surveys. We summarize our findings in Section 7.

Throughout this paper, a standard flat $\Lambda$CDM cosmology is used, with a matter density parameter $\Omega_M = 0.3$, a vacuum energy density

parameter $\Omega_\Lambda = 0.7$, and a Hubble constant of $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. All magnitudes are in AB system (Oke & Gunn 1983), and we adopt a solar metallicity value of $12 + \log (\text{O/H})_\odot = 8.69$. All the stellar metallicities are quoted relative to the solar abundance ($Z_\odot$) from Asplund et al. (2009), which has a composition by mass of $Z_\odot = 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).

## 2 SAMPLE

We rely on rest-frame FUV spectroscopy in order to estimate the ionizing properties of high-$z$ SFGs. At $3 \le z \le 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-$z$ 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).

### 2.1 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 $z = 1-6.5$. The VANDELS footprints are centred 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 catalogues 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–10 000 Å 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 redshifts in total. The final VANDELS catalogue reaches a target selection completeness of 40 per cent at $i_{\text{AB}} = 25$. VANDELS spectra have an unprecedented average S/N $\ge 7$ per resolution element over 80 per cent of the spectra, thanks to exposure times spanning from 20 up to 80 h on source. The VIMOS – Full-Width at Half Maximum – spectral resolution is $R \equiv \lambda/\Delta\lambda \approx 600$ at 5500 Å, with a wavelength dispersion of 2.5 Å per 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.

According to the selection criteria, the galaxies in VANDELS can be classified into passive/quiescent ($H_{\text{AB}} \le 22.5$, $1 \le z_{\text{phot}} < 2.5$), bright SFGs ($i_{\text{AB}} \le 25$, $2.4 \le z_{\text{phot}} < 5.5$), LBGs ($25 \le H_{\text{AB}} < 27$, $i_{\text{AB}} \le 27.5$, $3 \le z_{\text{phot}} < 6.5$), and a smaller subset of AGN candidates (Garilli et al. 2021). For non-AGN type galaxies, the spectroscopic redshifts are ‘flagged’ as a function of the reliability of the measurement as marginal ($\simeq 40$ per cent), fair ($\simeq 80$ per cent), or robust ($\ge 95$ per cent) reliability, corresponding to flag = 1, 2, 3/4, respectively (see Le Fèvre et al. 2015; Garilli et al. 2021, for more details). Based on a subset of VANDELS galaxies with C III]1908 detection, Llerena et al. (2022) found that the systemic redshifts are slightly larger than the spectroscopic redshifts, with a mean difference of 0.002.
et al. 2018; Pahl et al. 2021, 2023) at $z \simeq 3$, and the pioneering work by Izotov et al. (2016b, a, 2018a, b, 2021, 2022) and the recent Low-Redshift Lyman Continuum Survey (LzLCS; Flury et al. 2022a) at $z \simeq 0.3$. Particularly, compactness, high ionization parameters (traced by optical line ratios), high star-formation rate (SFR) surface densities, strong Ly $\alpha$ emission, and low dust-attenuation (traced by the UV slope) seem to characterize the strongest LCEs at low-$z$ (see also Wang et al. 2019; Izotov et al. 2020), with ionizing photon production efficiencies analogued to $z \simeq 6$ galaxies (Schaerer et al. 2016; Chisholm et al. 2022). Interestingly, and according to the analysis of new JWST observations by Mascia et al. (2023a), these characteristics resemble the properties of galaxies at the EoR (see also Endsley et al. 2022; Schaerer et al. 2022b; Cullen et al. 2023; Lin et al. 2023).

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 $f_{\text{esc}}^{\text{abs}}$ 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 $\alpha$ 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 II 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-$z$ SFGs and their evolution along with the different galaxy properties. For that, we make use of a sample of $\simeq 500$ deep, rest-frame ultraviolet spectra at $3 \le z \le 5$ drawn from the VANDELS survey (McLure et al. 2018b; Pentericci et al. 2018; Garilli et al. 2021). In particular, LyC absolute escape fractions ($f_{\text{esc}}^{\text{abs}}$) are derived by measuring the depth of the absorption lines in combination with the UV attenuation, whilst ionizing production efficiencies ($\xi_{\text{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 Section 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 Section 3. The main results of this research, looking for correlations between $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$ with the different galaxy properties are summarized in Section 4. Our results on the ionizing efficiency of high-$z$ galaxies are compared with different estimates in the literature in Section 5, and finally the possible redshift evolution of the $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ product is discussed in Section 6, by comparing our values against state-of-the-art low- and high-$z$ surveys. We summarize our findings in Section 7.

Throughout this paper, a standard flat $\Lambda$CDM cosmology is used, with a matter density parameter $\Omega_M = 0.3$, a vacuum energy density

parameter $\Omega_\Lambda = 0.7$, and a Hubble constant of $H_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)_\odot = 8.69$. All the stellar metallicities are quoted relative to the solar abundance ($Z_\odot$) from Asplund et al. (2009), which has a composition by mass of $Z_\odot = 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).

## 2 SAMPLE

We rely on rest-frame FUV spectroscopy in order to estimate the ionizing properties of high-$z$ SFGs. At $3 \le z \le 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-$z$ 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).

### 2.1 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 $z = 1-6.5$. The VANDELS footprints are centred 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 catalogues 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–10 000 Å 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 redshifts in total. The final VANDELS catalogue reaches a target selection completeness of 40 per cent at $i_{\text{AB}} = 25$. VANDELS spectra have an unprecedented average S/N $\ge 7$ per resolution element over 80 per cent of the spectra, thanks to exposure times spanning from 20 up to 80 h on source. The VIMOS–Full-Width at Half Maximum–spectral resolution is $R \equiv \lambda/\Delta\lambda \approx 600$ at 5500 Å, with a wavelength dispersion of 2.5 Å per 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.

According to the selection criteria, the galaxies in VANDELS can be classified into passive/quiescent ($H_{\text{AB}} \le 22.5$, $1 \le z_{\text{phot}} < 2.5$), bright SFGs ($i_{\text{AB}} \le 25$, $2.4 \le z_{\text{phot}} < 5.5$), LBGs ($25 \le H_{\text{AB}} < 27$, $i_{\text{AB}} \le 27.5$, $3 \le z_{\text{phot}} < 6.5$), and a smaller subset of AGN candidates (Garilli et al. 2021). For non-AGN type galaxies, the spectroscopic redshifts are ‘flagged’ as a function of the reliability of the measurement as marginal ($\simeq 40$ per cent), fair ($\simeq 80$ per cent), or robust ($\ge 95$ per cent) reliability, corresponding to flag = 1, 2, 3/4, respectively (see Le Fèvre et al. 2015; Garilli et al. 2021, for more details). Based on a subset of VANDELS galaxies with C III]1908 detection, Llerena et al. (2022) found that the systemic redshifts are slightly larger than the spectroscopic redshifts, with a mean difference of 0.002.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** H-band AB magnitude versus spectroscopic redshift ($z_{\text{spec}}$) for the VANDELS-DR4 (grey open circles, 2087 galaxies). Our working sample of 534 galaxies is highlighted with red filled symbols, and the vertical dashed line indicates the minimum redshift at which LyC imaging is available through the VIMOS/U-band. The top right and bottom histograms show the $H_{\text{AB}}$ and $z_{\text{spec}}$ distributions of both the original VANDELS-DR4 and working samples (shaded grey and red, respectively). Error bars comprise the 16th, 50th, and 84th percentiles of the distributions.</figcaption></figure>

Our final sample comprises VANDELS-DR4 galaxies under the SFGs and LBGs categories only, whose spectroscopic redshifts fall in the $3 \le z_{\text{spec}} \le 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$\alpha$1216 to C iv1550). We then compute the mean S/N per pixel over two regions free of absorption features, specifically 1350–1375 Å and 1450–1475 Å in the rest-frame (25 Å). Then, we select the galaxies with S/N $\ge$ 2 in the two spectral windows simultaneously ($\sim$50 per cent 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 $\times 1.4$.

In summary, our VANDELS working sample includes 534 galaxies at $3 \le z_{\text{spec}} \le 5$ with a median S/N $\approx 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 H-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 (16th and 84th percentiles) F160W and $z_{\text{spec}}$ of the working sample correspond to $\langle H_{\text{AB}} \rangle = 24.8_{-0.9}^{+0.8}$ AB and $\langle z_{\text{spec}} \rangle = 3.56_{-0.36}^{+0.41}$, respectively.

## 2.2 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

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Star-forming main sequence (MS) diagram i.e. log SFR – log M$_{\star}$, for the VANDELS-DR4 (grey open circles). Our working sample is indicated through red filled symbols, and the large blue circles display the SFR running medians at each interquartile 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 $z = 3.5$, but it still probes the upper bound of the SFG population at these redshifts.</figcaption></figure>

Inference and Parameter Estimation (BAGPIPES) code.¹BAGPIPES (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 time-scale of 10 Myr and a minimum age of 10 Myr. The stellar metallicity was fixed to 0.2Z$_{\odot}$, 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 modelled using the Salim, Boquien & Lee (2018) prescription, and nebular emission was considered by adopting a constant ionization parameter of log (U) = -3 in the fits.

For our VANDELS-DR4 working sample, the median stellar mass and SFR are $\langle \log(M_{\star}/M_{\odot}) \rangle = 9.34_{-0.39}^{+0.31}$ and $\langle \text{SFR}/(M_{\odot} \text{ yr}^{-1}) \rangle = 14.46_{-8.47}^{+25.30}$, respectively, with $0.1M_{\odot}$ and $2.4M_{\odot} \text{ yr}^{-1}$ typical uncertainties. The log SFR – log M$_{\star}$ distribution for the selected sample lays 0.1 to 0.2 dex above the SF main-sequence (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 25th, 50th, and 75th stellar mass interpercentile ranges fall above the Speagle et al. (2014) MS relation at $z = 3.5$. This means that, at comparable stellar masses (i.e. log (M$_{\star}$/M$_{\odot}$) = 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-z (see the discussion in Cullen et al. 2021). A similar

<footnote n="1">¹BAGPIPES (Carnall et al. 2018) is a state-of-the-art PYTHON code for modelling galaxy spectra and fitting spectroscopic and photometric data, for details go to https://bagpipes.readthedocs.io/en/latest/.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** H-band AB magnitude versus spectroscopic redshift ($z_{\text{spec}}$) for the VANDELS-DR4 (grey open circles, 2087 galaxies). Our working sample of 534 galaxies is highlighted with red filled symbols, and the vertical dashed line indicates the minimum redshift at which LyC imaging is available through the VIMOS/U-band. The top right and bottom histograms show the $H_{\text{AB}}$ and $z_{\text{spec}}$ distributions of both the original VANDELS-DR4 and working samples (shaded grey and red, respectively). Error bars comprise the 16th, 50th, and 84th percentiles of the distributions.</figcaption></figure>

Our final sample comprises VANDELS-DR4 galaxies under the SFGs and LBGs categories only, whose spectroscopic redshifts fall in the $3 \le z_{\text{spec}} \le 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$\alpha$1216 to C iv1550). We then compute the mean S/N per pixel over two regions free of absorption features, specifically 1350–1375 Å and 1450–1475 Å in the rest-frame (25 Å). Then, we select the galaxies with S/N $\ge$ 2 in the two spectral windows simultaneously ($\sim$50 per cent 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 $\times 1.4$.

In summary, our VANDELS working sample includes 534 galaxies at $3 \le z_{\text{spec}} \le 5$ with a median S/N $\approx 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 H-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 (16th and 84th percentiles) F160W and $z_{\text{spec}}$ of the working sample correspond to $\langle H_{\text{AB}} \rangle = 24.8_{-0.9}^{+0.8}$ AB and $\langle z_{\text{spec}} \rangle = 3.56_{-0.36}^{+0.41}$, respectively.

## 2.2 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

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Star-forming main sequence (MS) diagram i.e. log SFR – log M$_{\star}$, for the VANDELS-DR4 (grey open circles). Our working sample is indicated through red filled symbols, and the large blue circles display the SFR running medians at each interquartile 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 $z = 3.5$, but it still probes the upper bound of the SFG population at these redshifts.</figcaption></figure>

Inference and Parameter Estimation (BAGPIPES) code.¹BAGPIPES (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 time-scale of 10 Myr and a minimum age of 10 Myr. The stellar metallicity was fixed to 0.2Z$_{\odot}$, 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 modelled using the Salim, Boquien & Lee (2018) prescription, and nebular emission was considered by adopting a constant ionization parameter of log (U) = -3 in the fits.

For our VANDELS-DR4 working sample, the median stellar mass and SFR are $\langle \log(M_{\star}/M_{\odot}) \rangle = 9.34_{-0.39}^{+0.31}$ and $\langle \text{SFR}/(M_{\odot} \text{yr}^{-1}) \rangle = 14.46_{-8.47}^{+25.30}$, respectively, with $0.1M_{\odot}$ and $2.4M_{\odot} \text{yr}^{-1}$ typical uncertainties. The log SFR – log M$_{\star}$ distribution for the selected sample lays 0.1 to 0.2 dex above the SF main-sequence (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 25th, 50th, and 75th stellar mass interpercentile ranges fall above the Speagle et al. (2014) MS relation at $z = 3.5$. This means that, at comparable stellar masses (i.e. log (M$_{\star}/M_{\odot}$) = 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-z (see the discussion in Cullen et al. 2021). A similar

<footnote n="1">¹BAGPIPES (Carnall et al. 2018) is a state-of-the-art PYTHON code for modelling galaxy spectra and fitting spectroscopic and photometric data, for details go to https://bagpipes.readthedocs.io/en/latest/.</footnote>
behaviour was found in Calabrò et al. (2022). Finally, our sample also suffers from luminosity (or stellar mass) incompleteness, i.e. the highest redshifts objects ($z \gtrsim 4$) are biased towards the brightest AB absolute magnitudes because of the current flux-limiting nature of the VANDELS survey.

## 2.3 Spectral properties of the VANDELS galaxies: Ly$\alpha$ 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$\alpha$ equivalent width ($W_{Ly\alpha}$) of galaxies, either in emission or absorption.

The Ly$\alpha$ (H1216) flux and equivalent width were measured by fitting a single-Gaussian profile to the Ly$\alpha$ line using the dedicated SLINEFIT code.<fnref n="2" /> If the S/N of the line reach an input threshold, the script allows a $\pm 1000$ kms$^{-1}$ offset of the line centre respect to rest-frame value given by the spectroscopic redshift. The continuum level was defined by the Bruzual & Charlot (2003) best-fitting template to the entire spectrum (which gives a first order correction for underlying stellar absorption), and the flux and equivalent width was measured over $\pm 8000$ kms$^{-1}$ each side of the line peak (see Section 3.2 for the actual definition of the equivalent width). By convention, the fits report negative equivalent widths when the Ly$\alpha$ line appears in emission, whereas a positive value is reported when the line is absorption dominated. Given the limited variety of Ly$\alpha$ profiles observed in the VANDELS sample due to a medium $R \approx 600$ spectral resolution (see Kornei et al. 2010; Cullen et al. 2020, and methods therein), this approach only constitutes a first order estimation of $W_{Ly\alpha}$. We also note that, according to Talia et al. (in preparation), for 198 object ($\simeq 37$ per cent of the sample) a single Gaussian fit did not provide a good fit to the Ly$\alpha$ 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$\alpha$ flux. Doing so, the median Ly$\alpha$ equivalent width is $\langle W_{Ly\alpha} (\AA) \rangle = 4.58_{-9.43}^{+19.01}$ for our selected sample. According to the definition by Pentericci et al. (2009), 104 out of 534 of galaxies ($\simeq 20$ per cent of the sample) can be classified as LAEs with $W_{Ly\alpha} \ge -20$ Å (although other authors may adopt different criteria, e.g. see Stark, Ellis & Ouchi 2011; Nakajima et al. 2018a; Kusakabe et al. 2020).

## 2.4 Comparison samples of LCEs

Asides from our VANDELS sample, primarily composed by SFGs at $3 \le z \le 5$, we will consider as benchmarks other samples of confirmed LCEs from the literature (see Section 6). The LCEs comparison samples are the following:

**The Low-$z$ 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 \le z \le 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 $R \approx 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, a, 2018a,

Wang et al. 2019; Izotov et al. 2021). Out of 89 galaxies, 50 were detected in the LyC, with a median absolute escape fraction of $f_{esc}^{abs} = 0.04_{-0.02}^{+0.16}$, while the remaining 39 galaxies show strong $f_{esc}^{abs}$ upper limits, typically below 1 per cent.

**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 \le z \le 3.4$. The LRIS spectra of KLCS galaxies cover the LyC region and a wide FUV window longwards of the Ly$\alpha$ line (880–1650 Å), at a spectral resolution of $R \approx 800$ for LRIS-B and $R \approx 1400$ for LRIS-R. In the most recent HST imaging follow-up by Pahl et al. (2021, 2023), they re-analyse 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 $f_{esc}^{abs} = 0.06 \pm 0.01$.

# 3 METHODS

In this section, we first describe the FICUS code (Section 3.1), a customized PYTHON 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 ($f_{esc}^{abs}$), by using the information provided by the UV absorption lines and the global SED fit to the non-ionizing UV continuum (Section 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 (Section 3.3).

## 3.1 The FICUS code: Fitting the stellar continuum of UV spectra

With the purpose of quantifying the stellar continuum properties underlying the VANDELS spectra, we have created FICUS.<fnref n="3" /> FICUS is a customized PYTHON 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-fitting 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. (2018, 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.

### 3.1.1 Stellar population models

As described in SL22, the stellar continuum modelling 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, FICUS inputs the fully theoretical STARBURST99 single-star models without stellar rotation (S99; Leitherer et al. 2011, 2014) using the Geneva evolution models (Meynet et al. 1994), and computed with the WM-BASIC method (Pauldrach, Hoffmann & Lennon 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_{\odot}$. The spectral resolution of the S99 models is $R \equiv \lambda/\Delta\lambda \simeq 2500$, and it remains approximately constant at FUV wavelengths.

<footnote n="2">²SLINEFIT 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/cschreiber/slinefit, by Corentin Schreiber.</footnote>

<footnote n="3">³The FICUS code is publicly available and can be cloned from the author's GITHUB repository: https://github.com/asalda/FiCUS.git.</footnote>
behaviour was found in Calabrò et al. (2022). Finally, our sample also suffers from luminosity (or stellar mass) incompleteness, i.e. the highest redshifts objects ($z \gtrsim 4$) are biased towards the brightest AB absolute magnitudes because of the current flux-limiting nature of the VANDELS survey.

## 2.3 Spectral properties of the VANDELS galaxies: Ly$\alpha$ 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$\alpha$ equivalent width ($W_{Ly\alpha}$) of galaxies, either in emission or absorption.

The Ly$\alpha$ (H1216) flux and equivalent width were measured by fitting a single-Gaussian profile to the Ly$\alpha$ line using the dedicated SLINEFIT code.<fnref n="2" /> If the S/N of the line reach an input threshold, the script allows a $\pm 1000 \text{ kms}^{-1}$ offset of the line centre respect to rest-frame value given by the spectroscopic redshift. The continuum level was defined by the Bruzual & Charlot (2003) best-fitting template to the entire spectrum (which gives a first order correction for underlying stellar absorption), and the flux and equivalent width was measured over $\pm 8000 \text{ kms}^{-1}$ each side of the line peak (see Section 3.2 for the actual definition of the equivalent width). By convention, the fits report negative equivalent widths when the Ly$\alpha$ line appears in emission, whereas a positive value is reported when the line is absorption dominated. Given the limited variety of Ly$\alpha$ profiles observed in the VANDELS sample due to a medium $R \approx 600$ spectral resolution (see Kornei et al. 2010; Cullen et al. 2020, and methods therein), this approach only constitutes a first order estimation of $W_{Ly\alpha}$. We also note that, according to Talia et al. (in preparation), for 198 object ($\simeq 37$ per cent of the sample) a single Gaussian fit did not provide a good fit to the Ly$\alpha$ 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$\alpha$ flux. Doing so, the median Ly$\alpha$ equivalent width is $\langle W_{Ly\alpha} (\AA) \rangle = 4.58_{-9.43}^{+19.01}$ for our selected sample. According to the definition by Pentericci et al. (2009), 104 out of 534 of galaxies ($\simeq 20$ per cent of the sample) can be classified as LAEs with $W_{Ly\alpha} \ge -20 \AA$ (although other authors may adopt different criteria, e.g. see Stark, Ellis & Ouchi 2011; Nakajima et al. 2018a; Kusakabe et al. 2020).

## 2.4 Comparison samples of LCEs

Asides from our VANDELS sample, primarily composed by SFGs at $3 \le z \le 5$, we will consider as benchmarks other samples of confirmed LCEs from the literature (see Section 6). The LCEs comparison samples are the following:

**The Low-$z$ 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 \le z \le 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 $R \approx 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, a, 2018a,

Wang et al. 2019; Izotov et al. 2021). Out of 89 galaxies, 50 were detected in the LyC, with a median absolute escape fraction of $f_{esc}^{abs} = 0.04_{-0.02}^{+0.16}$, while the remaining 39 galaxies show strong $f_{esc}^{abs}$ upper limits, typically below 1 per cent.

**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 \le z \le 3.4$. The LRIS spectra of KLCS galaxies cover the LyC region and a wide FUV window longwards of the Ly$\alpha$ line (880–1650 Å), at a spectral resolution of $R \approx 800$ for LRIS-B and $R \approx 1400$ for LRIS-R. In the most recent HST imaging follow-up by Pahl et al. (2021, 2023), they re-analyse 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 $f_{esc}^{abs} = 0.06 \pm 0.01$.

# 3 METHODS

In this section, we first describe the FICUS code (Section 3.1), a customized PYTHON 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 ($f_{esc}^{abs}$), by using the information provided by the UV absorption lines and the global SED fit to the non-ionizing UV continuum (Section 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 (Section 3.3).

## 3.1 The FICUS code: Fitting the stellar continuum of UV spectra

With the purpose of quantifying the stellar continuum properties underlying the VANDELS spectra, we have created FICUS.<fnref n="3" /> FICUS is a customized PYTHON 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-fitting 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. (2018, 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.

### 3.1.1 Stellar population models

As described in SL22, the stellar continuum modelling 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, FICUS inputs the fully theoretical STARBURST99 single-star models without stellar rotation (S99; Leitherer et al. 2011, 2014) using the Geneva evolution models (Meynet et al. 1994), and computed with the WM-BASIC method (Pauldrach, Hoffmann & Lennon 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_{\odot}$. The spectral resolution of the S99 models is $R \equiv \lambda/\Delta\lambda \simeq 2500$, and it remains approximately constant at FUV wavelengths.

<footnote n="2">²SLINEFIT 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/cschreiber/slinefit, by Corentin Schreiber.</footnote>

<footnote n="3">³The FICUS code is publicly available and can be cloned from the author's GITHUB repository: https://github.com/asalda/FiCUS.git.</footnote>
Four different metallicities (0.05, 0.2, 0.4, and 1 Z⊙) 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-z 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 000 K (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 CLOUDY v17.0 code<fnref n="4" /> (Ferland et al. 2017), assuming similar gas-phase and stellar metallicities, an ionization parameter of log(U) = −2.5, and a volume hydrogen density of nH = 100 cm⁻³. 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 (≤10 Myr) at wavelengths ≥1200 Å, and produces redder stellar models than before, which has a pronounced effect on the fitted E<sub>B-V</sub> 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(U) values.

### 3.1.2 Fitting method

The observed spectra were first manually placed into the rest-frame by multiplying by the corresponding 1/(1+z) 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 (R(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 FICUS, the UV stellar continuum (F*(λ)) is fitted with a linear combination of multiple S99 models plus the nebular continuum. F<sub>ij</sub><sup>S99</sup> 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,<fnref n="5" /> this results in:

$$ 
\begin{aligned}
F^*(\lambda) &= 10^{-0.4k_\lambda E_{B-V}} \sum_{i,j} X_{ij} F_{ij}^{S99}, \\
& i \equiv 1, 2, 3, 4, 5, 8, 10, 15, 20, 40 \text{ Myr} \\
& j \equiv 0.05, 0.2, 0.4, 1 \text{ Z}_\odot
\end{aligned} 
\quad (2) 
$$

where 10⁻⁰.⁴Aλ is the UV attenuation term, Aλ = kλE<sub>B-V</sub>, and kλ is given by the adopted dust-attenuation law. The X<sub>ij</sub> linear coefficients determine the weight of single-stellar population within the fit, and the best fit is chosen through a non-linear χ² minimization algorithm with respect to the observed data (LMFIT<fnref n="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 N iterations per spectrum (we chose N = 100, enough realizations to sample the posterior continuum so that it approaches ‘Gaussianity’ on each pixel).

Fig. 3 shows the FICUS output for one of the galaxies in VANDELS (objID: CDFS017345, at z = 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 per cent of the total light at 3 Myr) 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 E<sub>B-V</sub> ≈ 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 vλ1240 and C ivλ1550 P-Cygni stellar wind profiles (see also O vλ1371 and Si ivλ1402). Note that the inclusion of C iv1550 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 vλ1240 is mostly sensitive to the age of the burst – with the P-Cygni profile being more prominent at younger ages – the C ivλ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-Z* degeneracy in the FUV, although the method is still affected by age-attenuation degeneracy effects.

### 3.1.3 Results

Given the best-fitting SED for the stellar continuum and the best-fitting X<sub>i</sub> coefficients in equation (2), FICUS also provides a handful of other secondary SED parameters that we described here.

**Stellar age and metallicity:** The light-weighted stellar age and metallicity (Z<sub>*</sub>) of the best-fitting stellar population can be easily obtained using the X<sub>k</sub> weights as:

$$ \text{Age} = \frac{\sum_k X_k \text{age}_k}{\sum_k X_k}, \quad \text{Z}^* = \frac{\sum_k X_k Z_k}{\sum_k X_k} \quad (3) $$

where age<sub>k</sub>, z<sub>k</sub> are the age and metallicity of k-th synthetic stellar model. For instance, in the example of Fig. 3, the light-weighted age is 4.80 Myr and the light-weighted metallicity is 0.18 Z⊙. The median stellar age and metallicity of the VANDELS sample are ⟨Age⟩ = 8.5<sup>+10.1</sup><sub>-4.4</sub> Myr and ⟨Z<sub>*</sub>/Z⊙⟩ = 0.23<sup>+0.18</sup><sub>-0.12</sub>, respectively. Cullen et al. (2019) and Calabrò et al. (2021) estimate ⟨log (Z<sub>*</sub>/Z⊙)⟩ ≈ −0.8 from independent methods (see white vertical line in Fig. 4). Compared to these works, our average metallicity is slightly higher: ⟨log (Z<sub>*</sub>/Z⊙)⟩ ≈ −0.63, although still compatible within the spread of the distributions.

<footnote n="4">⁴https://trac.nublado.org/</footnote>

<footnote n="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.</footnote>

<footnote n="6">⁶A PYTHON version of the LMFIT package can be found at https://lmfit.github.io/lmfit-py/.</footnote>
Four different metallicities (0.05, 0.2, 0.4, and 1 Z⊙) 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-z 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 000 K (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 CLOUDY v17.0 code<fnref n="4" /> (Ferland et al. 2017), assuming similar gas-phase and stellar metallicities, an ionization parameter of log(U) = −2.5, and a volume hydrogen density of nH = 100 cm⁻³. 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 (≤10 Myr) at wavelengths ≥1200 Å, and produces redder stellar models than before, which has a pronounced effect on the fitted E<sub>B-V</sub> 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(U) values.

### 3.1.2 Fitting method

The observed spectra were first manually placed into the rest-frame by multiplying by the corresponding 1/(1 + z) 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 (R(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 FICUS, the UV stellar continuum (F*(λ)) is fitted with a linear combination of multiple S99 models plus the nebular continuum. F<sub>ij</sub><sup>S99</sup> 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,<fnref n="5" /> this results in:

$$ 
\begin{aligned}
F^*(\lambda) &= 10^{-0.4k_{\lambda} E_{B-V}} \sum_{i,j} X_{ij} F_{ij}^{S99}, \\
&i \equiv 1, 2, 3, 4, 5, 8, 10, 15, 20, 40 \text{ Myr} \\
j \equiv 0.05, 0.2, 0.4, 1 \text{ Z}_{\odot}
\end{aligned} 
\quad (2) 
$$

where 10⁻⁰.⁴Aλ is the UV attenuation term, Aλ = kλE<sub>B-V</sub>, and kλ is given by the adopted dust-attenuation law. The X<sub>ij</sub> linear coefficients determine the weight of single-stellar population within the fit, and the best fit is chosen through a non-linear χ² minimization algorithm with respect to the observed data (LMFIT<fnref n="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 N iterations per spectrum (we chose N = 100, enough realizations to sample the posterior continuum so that it approaches ‘Gaussianity’ on each pixel).

Fig. 3 shows the FICUS output for one of the galaxies in VANDELS (objID: CDFS017345, at z = 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 per cent of the total light at 3 Myr) 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 E<sub>B-V</sub> ≈ 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 vλ1240 and C ivλ1550 P-Cygni stellar wind profiles (see also O vλ1371 and Si ivλ1402). Note that the inclusion of C iv1550 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 vλ1240 is mostly sensitive to the age of the burst – with the P-Cygni profile being more prominent at younger ages – the C ivλ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-Z* degeneracy in the FUV, although the method is still affected by age-attenuation degeneracy effects.

### 3.1.3 Results

Given the best-fitting SED for the stellar continuum and the best-fitting X<sub>i</sub> coefficients in equation (2), FICUS also provides a handful of other secondary SED parameters that we described here.

**Stellar age and metallicity:** The light-weighted stellar age and metallicity (Z<sub>*</sub>) of the best-fitting stellar population can be easily obtained using the X<sub>k</sub> weights as:

$$ \text{Age} = \frac{\sum_k X_k \text{age}_k}{\sum_k X_k}, \quad \text{Z}^* = \frac{\sum_k X_k Z_k}{\sum_k X_k} \quad (3) $$

where age<sub>k</sub>, z<sub>k</sub> are the age and metallicity of k-th synthetic stellar model. For instance, in the example of Fig. 3, the light-weighted age is 4.80 Myr and the light-weighted metallicity is 0.18 Z⊙. The median stellar age and metallicity of the VANDELS sample are ⟨Age⟩ = 8.5<sup>+10.1</sup><sub>-4.4</sub> Myr and ⟨Z<sub>*</sub>/Z⊙⟩ = 0.23<sup>+0.18</sup><sub>-0.12</sub>, respectively. Cullen et al. (2019) and Calabrò et al. (2021) estimate ⟨log (Z<sub>*</sub>/Z⊙)⟩ ≈ −0.8 from independent methods (see white vertical line in Fig. 4). Compared to these works, our average metallicity is slightly higher: ⟨log (Z<sub>*</sub>/Z⊙)⟩ ≈ −0.63, although still compatible within the spread of the distributions.

<footnote n="4">⁴https://trac.nublado.org/</footnote>

<footnote n="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.</footnote>

<footnote n="6">⁶A PYTHON version of the LMFIT package can be found at https://lmfit.github.io/lmfit-py/.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top: FICUS SED fitting results for the CDFS017345 VANDELS galaxy. The observed and error spectra are displayed in black and light green. The best-fitting stellar continuum is overplotted in red, while the spectral regions masked during the fit are shown in grey. The spectrum, shown in $F_\lambda$ units, has been normalized over 1350–1370 Å. The most prominent stellar features, nebular and ISM absorption lines are indicated with blue, gold, and black vertical lines at the top part of the figure. Bottom left-hand panel: light-fractions as a function of age and metallicity for the best-fitting composite stellar population model. Bottom right-hand panel: summary chart with the best-fitting SED fitting parameters, i.e. reduced-$\chi^2$, $E_{B-V}$ (in mag., R16), light-weighted stellar age (in Myr) and metallicity (in $Z_\odot$).</figcaption></figure>

Our shift in log $Z_*$ 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 Nv1240) of individual FUV galaxy-spectra (see Section 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 $z \simeq 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: $E_{B-V}$:** The amount of dust-attenuation is given by the $10^{-0.4A_\lambda}$ term in equation (2), where $A_\lambda = k_\lambda E_{B-V}$ is defined as the ‘UV attenuation’, and the specific functional form for $k_\lambda$ is determined by the dust-attenuation law. Therefore, the resulting best-fitting values for $E_{B-V}$ will differ depending on the assumed dust-law, and this will have important consequences for other related quantities which explicitly depend on $E_{B-V}$, like the absolute escape fraction ($f_{esc}^{abs}$) as we will discuss later on (see Section 6).

With the goal of testing the influence of the dust-attenuation law in our results, we use two extreme cases for $k_\lambda$ 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-z SFGs (Salim et al. 2018), low-metallicity starburst (Shivaei et al. 2020), and LCEs (Izotov et al. 2016a).

After performing the SED fits with FICUS using both the R16 and SMC prescriptions for $k_\lambda$, we find an average attenuation parameter of $\langle E_{B-V} \rangle = 0.22_{-0.08}^{+0.09}$ mag. for R16 and $\langle E_{B-V} \rangle = 0.10_{-0.03}^{+0.04}$ 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.4A_\lambda}$ is slightly higher for SMC at all wavelengths, meaning that the fitted $X_i$ coefficients are slightly lower for SMC than R16 law, so that both fits can match the observations similarly (reduced-$\chi^2$ distributions are similar). Lower $X_i$ coefficients for SMC means that all absolute quantities derived from the intrinsic SEDs will also change accordingly.

**UV-continuum slope (at 1500 Å): $\beta_{spec}^{1500}$:** The UV spectroscopic continuum slope at 1500 Å, also called $\beta$-slope, was computed by
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top: FICUS SED fitting results for the CDFS017345 VANDELS galaxy. The observed and error spectra are displayed in black and light green. The best-fitting stellar continuum is overplotted in red, while the spectral regions masked during the fit are shown in grey. The spectrum, shown in $F_\lambda$ units, has been normalized over 1350–1370 Å. The most prominent stellar features, nebular and ISM absorption lines are indicated with blue, gold, and black vertical lines at the top part of the figure. Bottom left-hand panel: light-fractions as a function of age and metallicity for the best-fitting composite stellar population model. Bottom right-hand panel: summary chart with the best-fitting SED fitting parameters, i.e. reduced-$\chi^2$, $E_{B-V}$ (in mag., R16), light-weighted stellar age (in Myr) and metallicity (in $Z_\odot$).</figcaption></figure>

Our shift in log $Z_*$ 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 Nv1240) of individual FUV galaxy-spectra (see Section 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 $z \simeq 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: $E_{B-V}$:** The amount of dust-attenuation is given by the $10^{-0.4A_\lambda}$ term in equation (2), where $A_\lambda = k_\lambda E_{B-V}$ is defined as the ‘UV attenuation’, and the specific functional form for $k_\lambda$ is determined by the dust-attenuation law. Therefore, the resulting best-fitting values for $E_{B-V}$ will differ depending on the assumed dust-law, and this will have important consequences for other related quantities which explicitly depend on $E_{B-V}$, like the absolute escape fraction ($f_{esc}^{abs}$) as we will discuss later on (see Section 6).

With the goal of testing the influence of the dust-attenuation law in our results, we use two extreme cases for $k_\lambda$ 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-$z$ SFGs (Salim et al. 2018), low-metallicity starburst (Shivaei et al. 2020), and LCEs (Izotov et al. 2016a).

After performing the SED fits with FICUS using both the R16 and SMC prescriptions for $k_\lambda$, we find an average attenuation parameter of $\langle E_{B-V} \rangle = 0.22_{-0.08}^{+0.09}$ mag. for R16 and $\langle E_{B-V} \rangle = 0.10_{-0.03}^{+0.04}$ 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.4A_\lambda}$ is slightly higher for SMC at all wavelengths, meaning that the fitted $X_i$ coefficients are slightly lower for SMC than R16 law, so that both fits can match the observations similarly (reduced-$\chi^2$ distributions are similar). Lower $X_i$ coefficients for SMC means that all absolute quantities derived from the intrinsic SEDs will also change accordingly.

**UV-continuum slope (at 1500 Å): $\beta_{spec}^{1500}$:** The UV spectroscopic continuum slope at 1500 Å, also called $\beta$-slope, was computed by
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Histograms showing the distribution of diverse secondary products resulting from our FICUS SED fitting to the VANDELS spectra. From top left-hand to bottom right-hand panel: light-weighted stellar Age, log Z, $E_{B-V}$, $\beta_{spec}^{1500}$, log $\xi_{ion}$, ($F_{\lambda 900}/F_{\lambda 1500}$)$_{int}$, measured R(LIS) and predicted log $f_{esc}^{abs}$. 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 16th, 50th, and 84th percentiles. The distributions of the mean residual flux of the LIS lines are shown in different colours because, by definition, this quantity does not depend on the attenuation law. Dashed vertical lines indicate log $\xi_{ion}$(Hz erg$^{-1}$) = 25.2 and $f_{esc}^{abs} = 5$ per cent.</figcaption></figure>

fitting a power law of the form $F_{\lambda} \propto \lambda^{\beta}$ to every individual spectrum (Meurer, Heckman & Calzetti 1999). To do so, we take the average flux density in seven 15 Å-width spectral windows, between 1275 and 1825 Å (similar to the ones in Calzetti, Kinney & Storchi-Bergmann 1994), and fit a linear relation to the log $F_{\lambda} - \log \lambda$

values using CURVE-FIT.<fnref n="7" /> For consistency, and in order to avoid inhomogeneous wavelength sampling for every source due to the different redshifts, $F_{\lambda}$ corresponds to the continuum flux obtained from the best-fitting SED model of the galaxy ($F^{*}(\lambda)$), which was described in the previous Section 3.1. This assumption imprints an intrinsic dependence of the $\beta$-slope on the chosen attenuation law, since it modifies the shape of the dust-attenuated spectrum through the $k_{\lambda}$ term. Thus, median of the $\beta$-slope distribution for the VANDELS sample is $\langle \beta_{spec}^{1500} \rangle = -1.34_{-0.42}^{+0.45}$ for R16 and $\langle \beta_{spec}^{1500} \rangle = -1.06_{-0.51}^{+0.54}$ for SMC. When compared to the slopes obtained in Calabrò et al. (2021) for the same objects ($\langle \beta_{1500} \rangle \simeq -1.80$), our spectroscopic $\beta$-slopes are overall redder (lower) than the former by $\simeq 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 multiband 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 Section 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 behaviour 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: $\xi_{ion}$:** The ionizing photon production efficiency is defined as the total number of ionizing photons per unit time produced by a radiation field ($Q$) normalized by its intrinsic UV luminosity ($L_{UV}$). In practice (see C19):

$$ \xi_{ion} [\text{Hz erg}^{-1}] = \frac{Q [\text{s}^{-1}]}{L_{\text{v1500}} [\text{erg s}^{-1} \text{ Hz}^{-1}]}, \quad (4) $$

where $L_{\text{v1500}}$ is the intrinsic (dust-free) SED luminosity at 1500 Å (in $F_{\lambda}$ units) and $Q$ is the number of ionizing photons produced by the best-fitting stellar population. $Q$ is calculated as the integral of the intrinsic SED below the Lyman limit ($\lambda < \lambda(H^0) = 912$ Å, or equivalently $E > E(H^0) = 13.6$ eV): $Q = \int_{<\lambda(H^0)} F^{*}(\lambda)d\lambda$. When averaged over the whole galaxy population and integrated over the full range of UV luminosities, $\xi_{ion}$ can be used to compute the total emissivity of ionizing photons ($\bar{n}_{ion}$) at any redshift (see equation 1).

For typical SF galaxies that can be described through a single-burst of star-formation, S99 templates predict an exponentially declining $\xi_{ion}$ with age, whose log $\xi_{ion}$ values span over 26–24 Hz erg$^{-1}$ for 1–10 Myr stellar populations at any metallicity, dropping dramatically towards ages older than 10 Myr (C19). In contrast, mixed-age models can provide systematically higher $\xi_{ion}$ for a given mean age than instantaneous burst models at the same evolutionary stage. Accordingly, our mixed-age fits give a median $\xi_{ion}$ of $\langle \log \xi_{ion} \rangle = 25.30_{-0.42}^{+0.28}$ Hz erg$^{-1}$ using the R16 law. $\xi_{ion}$ does not strongly depend on the assumed attenuation law because it has been calculated as the ratio between two absolute quantities ($Q$ and $L_{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 $\xi_{ion}$ stays unaltered.

<footnote n="7">⁷CURVE-FIT is a PYTHON-based optimization algorithm included in the SciPy package for scientific analysis: in https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html.⁷</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Histograms showing the distribution of diverse secondary products resulting from our FICUS SED fitting to the VANDELS spectra. From top left-hand to bottom right-hand panel: light-weighted stellar Age, log Z*, E<sub>B-V</sub>, β<sub>spec</sub><sup>1500</sup>, log ξ<sub>ion</sub>, (F<sub>λ900</sub>/F<sub>λ1500</sub>)<sub>int</sub>, measured R(LIS) and predicted log f<sub>esc</sub><sup>abs</sup>. 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 16th, 50th, and 84th percentiles. The distributions of the mean residual flux of the LIS lines are shown in different colours because, by definition, this quantity does not depend on the attenuation law. Dashed vertical lines indicate log ξ<sub>ion</sub>(Hz erg<sup>-1</sup>) = 25.2 and f<sub>esc</sub><sup>abs</sup> = 5 per cent.</figcaption></figure>

fitting a power law of the form $F_{\lambda} \propto \lambda^{\beta}$ to every individual spectrum (Meurer, Heckman & Calzetti 1999). To do so, we take the average flux density in seven 15 Å-width spectral windows, between 1275 and 1825 Å (similar to the ones in Calzetti, Kinney & Storchi-Bergmann 1994), and fit a linear relation to the log $F_{\lambda} - \log \lambda$

values using CURVE-FIT.<fnref n="7" /> For consistency, and in order to avoid inhomogeneous wavelength sampling for every source due to the different redshifts, $F_{\lambda}$ corresponds to the continuum flux obtained from the best-fitting SED model of the galaxy ($F^{*}(\lambda)$), which was described in the previous Section 3.1. This assumption imprints an intrinsic dependence of the β-slope on the chosen attenuation law, since it modifies the shape of the dust-attenuated spectrum through the $k_{\lambda}$ term. Thus, median of the β-slope distribution for the VANDELS sample is $\langle \beta_{spec}^{1500} \rangle = -1.34_{-0.42}^{+0.45}$ for R16 and $\langle \beta_{spec}^{1500} \rangle = -1.06_{-0.51}^{+0.54}$ for SMC. When compared to the slopes obtained in Calabrò et al. (2021) for the same objects ($\langle \beta_{1500} \rangle \simeq -1.80$), our spectroscopic β-slopes are overall redder (lower) than the former by $\simeq 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 multiband 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 Section 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 behaviour 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: ξ<sub>ion</sub>:** The ionizing photon production efficiency is defined as the total number of ionizing photons per unit time produced by a radiation field (*Q*) normalized by its intrinsic UV luminosity (*L*<sub>UV</sub>). In practice (see C19):

$$ \xi_{\text{ion}} [\text{Hz erg}^{-1}] = \frac{Q [\text{s}^{-1}]}{L_{\text{v1500}} [\text{erg s}^{-1} \text{ Hz}^{-1}]}, \quad (4) $$

where $L_{v1500}$ is the intrinsic (dust-free) SED luminosity at 1500 Å (in $F_v$ units) and $Q$ is the number of ionizing photons produced by the best-fitting stellar population. $Q$ is calculated as the integral of the intrinsic SED below the Lyman limit ($\lambda < \lambda(H^0) = 912$ Å, or equivalently $E > E(H^0) = 13.6$ eV): $Q = \int_{<\lambda(H^0)} F^*(\lambda)d\lambda$. When averaged over the whole galaxy population and integrated over the full range of UV luminosities, $\xi_{ion}$ can be used to compute the total emissivity of ionizing photons ($\bar{n}_{ion}$) at any redshift (see equation 1).

For typical SF galaxies that can be described through a single-burst of star-formation, S99 templates predict an exponentially declining $\xi_{ion}$ with age, whose log $\xi_{ion}$ values span over 26–24 Hz erg<sup>-1</sup> for 1–10 Myr stellar populations at any metallicity, dropping dramatically towards ages older than 10 Myr (C19). In contrast, mixed-age models can provide systematically higher $\xi_{ion}$ for a given mean age than instantaneous burst models at the same evolutionary stage. Accordingly, our mixed-age fits give a median $\xi_{ion}$ of $\langle \log \xi_{ion} \rangle = 25.30_{-0.42}^{+0.28}$ Hz erg<sup>-1</sup> using the R16 law. $\xi_{ion}$ does not strongly depend on the assumed attenuation law because it has been calculated as the ratio between two absolute quantities ($Q$ and $L_{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 $\xi_{ion}$ stays unaltered.

<footnote n="7">⁷CURVE-FIT is a PYTHON-based optimization algorithm included in the SCIPY package for scientific analysis: in https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html.⁷</footnote>
**Intrinsic 900-to-1500 Å flux ratio: ($F_{\lambda 900}/F_{\lambda 1500}$)int:** The ionizing-to-non-ionizing 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_{\lambda 900}/F_{\lambda 1500}$)int $\le 2$, exponentially declining with age down to negligible values at ages older than 20 Myr, where $F_{\lambda 900}^{\text{int}} \approx 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 $\langle F_{\lambda 900}/F_{\lambda 1500} \rangle_{\text{int}} = 0.66_{-0.32}^{+0.41}$. This quantity does not show any dependence on the attenuation law for the same reasons as $\xi_{\text{ion}}$.

In photometric campaigns, when searching for LCEs at high-z (e.g. Grazian et al. 2016, 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 $F_{\text{v}}$ units. Typically assumed values in the literature are $\langle F_{\text{v1500}}/F_{\text{v900}} \rangle_{\text{int}} = 3-5$ (see e.g. Guaita et al. 2016). We obtain a median value of $\langle F_{\text{v1500}}/F_{\text{v900}} \rangle_{\text{int}} \approx 4.2$ for the VANDELS sample at $3 \le z \le 5$ (Fig. 4), in agreement with the previous studies.

Fig. 4 offers a summary of the above-mentioned distributions: Age, $Z_{\text{*}}$, $E_{\text{B-V}}$, $\beta_{\text{spec}}^{\text{1500}}$, $\xi_{\text{ion}}$, and ($F_{\lambda 900}/F_{\lambda 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 16th to 84th percentiles of each distribution, whose values have been quoted with respect to the median in the current section.

## 3.2 Predicting ionizing escape fractions using UV absorption lines

### 3.2.1 The picket-fence model

The observation that neutral lines of hydrogen and other strong low-ionization 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, 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, Inoue & Jensen 2013), where both the neutral and low-ionized enriched gas are distributed in high-column-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, $C_f(\lambda_i)$. 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:

$$R_{\lambda_i} = 1 - C_f(\lambda_i). \qquad (5)$$

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 ($E_{\text{B-V}}$), the escape fraction of ionizing photons ($f_{\text{esc}}^{\text{abs, LIS}}$) 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):

$$f_{\text{esc}}^{\text{abs, LIS}} = 10^{-0.4k_{912}E_{\text{B-V}}} \times (a \times R_{\lambda_i} + b), \qquad (6)$$

where the $E_{\text{B-V}}$ is the UV dust-attenuation parameter measured according to the methods described in the previous Section 3.1, using a uniform screen geometry (see equation 2), and $R_{\lambda_i}$ is the measured residual flux of the $\lambda_i$ transition. The [a, b] coefficients correspond to the LIS to H I 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-z 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 $f_{\text{esc}}^{\text{abs}}$ of galaxies at any redshift. The depth of the LIS lines have also been applied to predict the escape fraction of high-z galaxy composites in Steidel et al. (2018) and Pahl et al. (2021), and some individual high-z spectra in Chisholm et al. (2018), with reasonable agreement with the observed $f_{\text{esc}}^{\text{abs}}$ values.

However, this approach is not without of caveats (see Section 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-the-art galaxy simulations. In principle, all these assumptions and model dependencies can contribute to explain the scatter seen in the observed R(LIS) – $f_{\text{esc}}^{\text{abs}}$ relation (see SL22), but this model could still be a good approximation for unresolved high-z studies, where the LyC emission coming from a single sight-line usually dominates.

### 3.2.2 Absorption line measurements

Our goal is to apply equation (6) to VANDELS spectra in order to indirectly estimate their ionizing escape fractions. To do so, we first measured the equivalent widths ($W_{\lambda}$) and residual fluxes ($R_{\lambda}$) for a set of $\times 4$ UV LIS lines, namely Si II1260, O I+Si II1302, C II1334 and Si II1527, which had simultaneous wavelength coverage in all the spectra at $3 \le z \le 5$. The equivalent width ($W_{\lambda_i}$) was then computed individually for every absorption line following Trainor et al. (2019):

$$W_{\lambda_i} = \int_{\Delta\lambda_i} (1 - f_{\lambda}/F^{*}(\lambda)) d\lambda, \qquad (7)$$

where $f_{\lambda}$ is the observed spectral flux density and $F^{*}(\lambda)$ is the modelled stellar continuum. The integration window ($\Delta\lambda_i$) was defined as $\pm 1000 \text{ kms}^{-1}$ from the wavelength of minimum depth for the line in question. Then, the residual flux was measured as the median over a narrow velocity interval of $\pm 150 \text{ kms}^{-1}$ around the minimum flux of the line, or equivalently:

$$R_{\lambda_i} = \langle f_{\lambda_i} / F^{*} \rangle. \qquad (8)$$

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
**Intrinsic 900-to-1500 Å flux ratio: ($F_{\lambda900}/F_{\lambda1500}$)int:** The ionizing-to-non-ionizing 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_{\lambda900}/F_{\lambda1500}$)int ≤ 2, exponentially declining with age down to negligible values at ages older than 20 Myr, where $F_{\lambda900}^{int} \approx 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 $\langle F_{\lambda900}/F_{\lambda1500} \rangle_{int} = 0.66_{-0.32}^{+0.41}$. This quantity does not show any dependence on the attenuation law for the same reasons as $\xi_{ion}$.

In photometric campaigns, when searching for LCEs at high-z (e.g. Grazian et al. 2016, 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 $F_{\nu}$ units. Typically assumed values in the literature are $\langle F_{\nu1500}/F_{\nu900} \rangle_{int} = 3-5$ (see e.g. Guaita et al. 2016). We obtain a median value of $\langle F_{\nu1500}/F_{\nu900} \rangle_{int} \approx 4.2$ for the VANDELS sample at $3 \le z \le 5$ (Fig. 4), in agreement with the previous studies.

Fig. 4 offers a summary of the above-mentioned distributions: Age, $Z_{\alpha}$, $E_{B-V}$, $\beta_{spec}^{1500}$, $\xi_{ion}$, and ($F_{\lambda900}/F_{\lambda1500}$)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 16th to 84th percentiles of each distribution, whose values have been quoted with respect to the median in the current section.

## 3.2 Predicting ionizing escape fractions using UV absorption lines

### 3.2.1 The picket-fence model

The observation that neutral lines of hydrogen and other strong low-ionization 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, 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, Inoue & Jensen 2013), where both the neutral and low-ionized enriched gas are distributed in high-column-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, $C_f(\lambda_i)$. 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:

$$R_{\lambda_i} = 1 - C_f(\lambda_i). \qquad (5)$$

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 ($E_{B-V}$), the escape fraction of ionizing photons ($f_{esc}^{abs, LIS}$) 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):

$$f_{esc}^{abs, LIS} = 10^{-0.4k_{912}E_{B-V}} \times (a \times R_{\lambda_i} + b), \qquad (6)$$

where the $E_{B-V}$ is the UV dust-attenuation parameter measured according to the methods described in the previous Section 3.1, using a uniform screen geometry (see equation 2), and $R_{\lambda_i}$ is the measured residual flux of the $\lambda_i$ transition. The [a, b] coefficients correspond to the LIS to H I 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-z 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 $f_{esc}^{abs}$ of galaxies at any redshift. The depth of the LIS lines have also been applied to predict the escape fraction of high-z galaxy composites in Steidel et al. (2018) and Pahl et al. (2021), and some individual high-z spectra in Chisholm et al. (2018), with reasonable agreement with the observed $f_{esc}^{abs}$ values.

However, this approach is not without of caveats (see Section 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-the-art galaxy simulations. In principle, all these assumptions and model dependencies can contribute to explain the scatter seen in the observed R(LIS) – $f_{esc}^{abs}$ relation (see SL22), but this model could still be a good approximation for unresolved high-z studies, where the LyC emission coming from a single sight-line usually dominates.

### 3.2.2 Absorption line measurements

Our goal is to apply equation (6) to VANDELS spectra in order to indirectly estimate their ionizing escape fractions. To do so, we first measured the equivalent widths ($W_{\lambda}$) and residual fluxes ($R_{\lambda}$) for a set of × 4 UV LIS lines, namely Si II1260, O I+Si II1302, C II1334 and Si II1527, which had simultaneous wavelength coverage in all the spectra at $3 \le z \le 5$. The equivalent width ($W_{\lambda_i}$) was then computed individually for every absorption line following Trainor et al. (2019):

$$W_{\lambda_i} = \int_{\Delta\lambda_i} (1 - f_{\lambda}/F^{*}(\lambda)) d\lambda, \qquad (7)$$

where $f_{\lambda}$ is the observed spectral flux density and $F^{*}(\lambda)$ is the modelled stellar continuum. The integration window ($\Delta\lambda_i$) was defined as $\pm 1000 \text{ kms}^{-1}$ from the wavelength of minimum depth for the line in question. Then, the residual flux was measured as the median over a narrow velocity interval of $\pm 150 \text{ kms}^{-1}$ around the minimum flux of the line, or equivalently:

$$R_{\lambda_i} = \langle f_{\lambda_i} / F^{*} \rangle. \qquad (8)$$

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 II 1260-to-Si II 1527 equivalent-width ratio which is compatible within the optically thick limits ($W_{1527}/W_{1260} \le 2.55$) given by the theoretical curve-of-growth for these transitions (Draine 2011). We conclude that the ISM conditions are such that Si II is optically thick (saturated), and we assume that C II is also saturated for typical Si/C abundances (see discussion in Mauerhofer et al. 2021). Calabrò et al. (2022) independently found optically thick Si II ($\log(N_{SiII}/\text{cm}^{-2}) > 12.7$) along a sample of VANDELS C III emitters at similar redshifts, whereas higher ionization lines like Si IV 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 ($R$) tends to make the lines broader and less deep, and the actual residual flux ($R_{\lambda i}$) can be overestimated. To account for these effects we performed mock absorption line simulations. Considering the picket-fence model with a uniform dust-screen, we simulated Si II 1260 absorption Voigt profiles assuming Gaussian distributions for the column density of Si II ($N_{SiII}$, cm$^{-2}$), the Doppler broadening parameter ($b$, kms$^{-1}$), and the line velocity shift ($v$, kms$^{-1}$). The simulated spectra were then degraded to $R$ (VIMOS) $\approx 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 $\times 100$ simulations were run for every of the input covering fractions in the $C_f(\lambda_i) = 0-1$ range, and the median measured residual flux of each distribution was then compared to the theoretical value of $1-C_f$. See Appendix A for a more comprehensive view of our simulations.

As a result, a theoretical residual flux of $(1-C_f) = 0.1, 0.3$, and 0.5 would be observed as $R(\text{LIS}) = 0.25, 0.4$, and 0.55 due to the smearing effect of the VIMOS instrumental resolution (i.e. 60 per cent, 25 per cent, and 10 per cent relative error), while a negligible correction will be applied at $(1-C_f) \ge 0.6$ because this $C_f$ regime is dominated by the S/N of the spectra. The resulting calibration (equation A3) was applied as a correction factor to the individually measured residual fluxes.

### 3.2.3 Indirect $f_{\text{esc}}^{\text{abs}}$ estimations

The mean residual flux of the LIS lines, $R(\text{LIS})$, is calculated by taking the average of the individual Si II 1260, O I + Si II 1302, C II 1334, and Si II 1527 line's depths. If these lines were all saturated, their residual fluxes should be very similar and the mean is representative, allowing us to gain S/N over the measurement. E.g. the average LIS residual flux of CDFS01735 (Fig. 3) is $R(\text{LIS}) \approx 0.69$ (0.75 previous to the correction by resolution).

Fig. 4 shows the histogram of the average $R(\text{LIS})$ distribution, characterized by $\langle R(\text{LIS}) \rangle = 0.31_{-0.23}^{+0.22}$ ($0.35_{-0.20}^{+0.19}$ without resolution correction). Then, using the dust-attenuation parameter ($E_{B-V}$) provided by the FICUS SED fits, we applied equation (6) to every single galaxy in our sample. For the [a, b] coefficients in this equation, we chose the calibration by Gazagnes et al. (2018) i.e. $[a, b] = [(0.63 \pm 0.12), (0.44 \pm 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 $f_{\text{esc}}^{\text{abs}} \approx 16$ per cent, adopting the R16 attenuation law.

The predicted $f_{\text{esc}}^{\text{abs}}$ distribution following this method is also presented in Fig. 4. As expected, the inferred $f_{\text{esc}}^{\text{abs}}$ values show a clear dependence on the dust-attenuation law, with the median of the R16

distribution $\times 1.5$ lower than the SMC law $f_{\text{esc}}^{\text{abs}}$ distribution. When using the R16 dust-attenuation law, the median ionizing absolute escape fraction, $\langle f_{\text{esc}}^{\text{abs}} \rangle$, for our VANDELS sample of 534 SFGs at $3 \le z \le 5$ is:

$$ \langle f_{\text{esc}}^{\text{abs}} \rangle = 0.02 \pm 0.01 \ (0.03 \pm 0.02), $$

while, when using the SMC attenuation law, it results in:

$$ \langle f_{\text{esc}}^{\text{abs}} \rangle = 0.03 \pm 0.02(0.04 \pm 0.03), $$

where the numbers in brackets correspond to the median $f_{\text{esc}}^{\text{abs}}$ previous to any resolution correction on the residual flux of the lines (see equation 6). As a comparison, Begley et al. (2022) obtained $\langle f_{\text{esc}}^{\text{abs}} \rangle = 0.07 \pm 0.02$ for a similar sample of $z \approx 3.5$ VANDELS galaxies with deep VIMOS/U-band observations (see Fig. 4), that is, a factor of $\times 2$ higher than our median value using the SMC law, although compatible within $1\sigma$ uncertainty. Moreover, our average $f_{\text{esc}}^{\text{abs}}$ is in agreement with the extrapolations of the recent Trebitsch et al. (2022) model at $z \approx 4$. In Section 6, we will compare our results with the escape fraction derived from other surveys targeting LCEs in the literature, at lower and higher redshifts.

## 3.3 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.

Following Cullen et al. (2019), Calabrò et al. (2021), Llerena et al. (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. $\langle z_{\text{spec}} \rangle = 3.56$, the resulting spectral binning was chosen to be $2.5\AA/(1 + (z_{\text{spec}})) = 0.55\AA$ ($2.5\AA$ equals the VIMOS wavelength dispersion per resolution element). Before co-adding the spectra, they were normalized to the mean flux in the 1350–1370 Å rest-frame interval. The final (normalized) flux at each wavelength was taken as the unweighted median of all the individual flux values after a regular $3\sigma$ 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 ($M_{1500}$), UV intrinsic luminosity ($M_{1500}^{\text{int}}$), UV-continuum slope ($\beta_{\text{spec}}^{1500}$), and Ly$\alpha$ equivalent width ($W_{Ly\alpha}$). Concretely, we divided the sample according to the 25th, 50th, and 75th percentiles (quartiles) of every property, resulting in four sub-samples for each quantity sorted as Q1, Q2, Q3, and Q4. Then, the resulting stacks were processed through the FICUS code and all the secondary SED parameters described in Section 3.1 were obtained. Appendix B, Table B1 contains the main properties and inferred SED parameters ($f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$) of the different composites.

**The effect of stacking on the absorption lines:** Even though the equivalent widths of the H I 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 blue-shifted line centres 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
galaxies have a Si II 1260-to-Si II 1527 equivalent-width ratio which is compatible within the optically thick limits ($W_{1527}/W_{1260} \le 2.55$) given by the theoretical curve-of-growth for these transitions (Draine 2011). We conclude that the ISM conditions are such that Si II is optically thick (saturated), and we assume that C II is also saturated for typical Si/C abundances (see discussion in Mauerhofer et al. 2021). Calabrò et al. (2022) independently found optically thick Si II ($\log(N_{SiII}/\text{cm}^{-2}) > 12.7$) along a sample of VANDELS C III emitters at similar redshifts, whereas higher ionization lines like Si IV 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 ($R$) tends to make the lines broader and less deep, and the actual residual flux ($R_{\lambda i}$) can be overestimated. To account for these effects we performed mock absorption line simulations. Considering the picket-fence model with a uniform dust-screen, we simulated Si II 1260 absorption Voigt profiles assuming Gaussian distributions for the column density of Si II ($N_{SiII}$, cm$^{-2}$), the Doppler broadening parameter ($b$, kms$^{-1}$), and the line velocity shift ($v$, kms$^{-1}$). The simulated spectra were then degraded to $R$ (VIMOS) $\approx 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 $\times 100$ simulations were run for every of the input covering fractions in the $C_f(\lambda_i) = 0-1$ range, and the median measured residual flux of each distribution was then compared to the theoretical value of $1-C_f$. See Appendix A for a more comprehensive view of our simulations.

As a result, a theoretical residual flux of $(1-C_f) = 0.1, 0.3$, and 0.5 would be observed as $R(\text{LIS}) = 0.25, 0.4$, and 0.55 due to the smearing effect of the VIMOS instrumental resolution (i.e. 60 per cent, 25 per cent, and 10 per cent relative error), while a negligible correction will be applied at $(1-C_f) \ge 0.6$ because this $C_f$ regime is dominated by the S/N of the spectra. The resulting calibration (equation A3) was applied as a correction factor to the individually measured residual fluxes.

### 3.2.3 Indirect $f_{\text{esc}}^{\text{abs}}$ estimations

The mean residual flux of the LIS lines, $R(\text{LIS})$, is calculated by taking the average of the individual Si II 1260, O I + Si II 1302, C II 1334, and Si II 1527 line's depths. If these lines were all saturated, their residual fluxes should be very similar and the mean is representative, allowing us to gain S/N over the measurement. E.g. the average LIS residual flux of CDFS01735 (Fig. 3) is $R(\text{LIS}) \approx 0.69$ (0.75 previous to the correction by resolution).

Fig. 4 shows the histogram of the average $R(\text{LIS})$ distribution, characterized by $\langle R(\text{LIS}) \rangle = 0.31_{-0.23}^{+0.22}$ ($0.35_{-0.20}^{+0.19}$ without resolution correction). Then, using the dust-attenuation parameter ($E_{B-V}$) provided by the FICUS SED fits, we applied equation (6) to every single galaxy in our sample. For the [a, b] coefficients in this equation, we chose the calibration by Gazagnes et al. (2018) i.e. $[a, b] = [(0.63 \pm 0.12), (0.44 \pm 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 $f_{\text{esc}}^{\text{abs}} \approx 16$ per cent, adopting the R16 attenuation law.

The predicted $f_{\text{esc}}^{\text{abs}}$ distribution following this method is also presented in Fig. 4. As expected, the inferred $f_{\text{esc}}^{\text{abs}}$ values show a clear dependence on the dust-attenuation law, with the median of the R16

distribution $\times 1.5$ lower than the SMC law $f_{\text{esc}}^{\text{abs}}$ distribution. When using the R16 dust-attenuation law, the median ionizing absolute escape fraction, $\langle f_{\text{esc}}^{\text{abs}} \rangle$, for our VANDELS sample of 534 SFGs at $3 \le z \le 5$ is:

$$ \langle f_{\text{esc}}^{\text{abs}} \rangle = 0.02 \pm 0.01 (0.03 \pm 0.02), $$

while, when using the SMC attenuation law, it results in:

$$ \langle f_{\text{esc}}^{\text{abs}} \rangle = 0.03 \pm 0.02(0.04 \pm 0.03), $$

where the numbers in brackets correspond to the median $f_{\text{esc}}^{\text{abs}}$ previous to any resolution correction on the residual flux of the lines (see equation 6). As a comparison, Begley et al. (2022) obtained $\langle f_{\text{esc}}^{\text{abs}} \rangle = 0.07 \pm 0.02$ for a similar sample of $z \approx 3.5$ VANDELS galaxies with deep VIMOS/U-band observations (see Fig. 4), that is, a factor of $\times 2$ higher than our median value using the SMC law, although compatible within $1\sigma$ uncertainty. Moreover, our average $f_{\text{esc}}^{\text{abs}}$ is in agreement with the extrapolations of the recent Trebitsch et al. (2022) model at $z \approx 4$. In Section 6, we will compare our results with the escape fraction derived from other surveys targeting LCEs in the literature, at lower and higher redshifts.

## 3.3 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.

Following Cullen et al. (2019), Calabrò et al. (2021), Llerena et al. (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. $\langle z_{\text{spec}} \rangle = 3.56$, the resulting spectral binning was chosen to be $2.5\AA/(1 + (z_{\text{spec}})) = 0.55\AA$ ($2.5\AA$ equals the VIMOS wavelength dispersion per resolution element). Before co-adding the spectra, they were normalized to the mean flux in the 1350–1370 Å rest-frame interval. The final (normalized) flux at each wavelength was taken as the unweighted median of all the individual flux values after a regular $3\sigma$ 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 ($M_{1500}$), UV intrinsic luminosity ($M_{1500}^{\text{int}}$), UV-continuum slope ($\beta_{\text{spec}}^{1500}$), and Ly$\alpha$ equivalent width ($W_{Ly\alpha}$). Concretely, we divided the sample according to the 25th, 50th, and 75th percentiles (quartiles) of every property, resulting in four sub-samples for each quantity sorted as Q1, Q2, Q3, and Q4. Then, the resulting stacks were processed through the FICUS code and all the secondary SED parameters described in Section 3.1 were obtained. Appendix B, Table B1 contains the main properties and inferred SED parameters ($f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$) of the different composites.

**The effect of stacking on the absorption lines:** Even though the equivalent widths of the H I 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 blue-shifted line centres 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 Appendix A, we randomly generated a set of $N = 100$ simulated $\mathrm{Si\text{i}I}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 ($R$) (VIMOS) ≈ 600) and a constant $S/N = 5$ were also incorporated in the mock realizations. After that, we measured the depth of the composite $\mathrm{Si\text{i}I}1260$ line profile and compared it to the input value given by the covering fraction. According to the our simulations (equation A4), an input ($1 - C_f$) = 0.1, 0.3, 0.5, and 0.8 would require corrections factors as large as 70 per cent, 40 per cent, 25 per cent, and 10 per cent 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 lines – might actually overpredict this effect. For that reason, we prefer not to apply such stacking corrections to our line measurements, but we encourage the reader to check Appendix 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).

# 4 RESULTS

The following paragraphs describe on the main results of this paper. In Section 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 Section 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 Section 4.3, the global ionizing properties of the LCEs versus non-LCEs samples are discussed. Lastly, Section 4.4 summarizes the properties of previously known LCEs reported in the literature which were included in our sample.

## 4.1 The ionizing escape fraction and production efficiency dependence with observed galaxy properties

Fig. 5 investigates how our predicted ionizing absolute escape fraction ($f_{\text{esc}}^{\text{abs}}$) and the ionizing photon production efficiency ($\xi_{\text{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 $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{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-$\tau$ statistics was applied to each pair of variables in order to figure out the strength ($\tau$) and significance ($p$-value)

of the correlation. For a sample of ≈500 objects, we consider correlations to be significant if $p_{\text{val.}} \le 10^{-3}$ (3$\sigma$) and strong if $|\tau| \gtrsim 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 $\tau$ coefficients for significant correlations (thick-framed panels), otherwise we write $p_{\text{val.}} > 10^{-3}$.

### 4.1.1 Observed absolute UV magnitude

The predicted $f_{\text{esc}}^{\text{abs}}$ and the observed UV absolute magnitude ($M_{1500}$) for the VANDELS sample at $3 \le z \le 5$ do not show any clear correlation when considering the whole range of UV magnitudes. Tanvir et al. (2019) found no significant correlation of $N_{\text{H}\text{i}}$ 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 equation (6), a lack of correlation between $f_{\text{esc}}^{\text{abs}}$ and $M_{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 $f_{\text{esc}}^{\text{abs}}$ is naturally expected for the faintest galaxies (see Section 7.2 in SL22).

Even though the correlation is tentative, $M_{1500} \gtrsim -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-$z$ by the KLCS (Steidel et al. 2018; Pahl et al. 2021, 2023) as well as observed at low-$z$ by the LzLCS (Chisholm et al. 2022; Flury et al. 2022a, and SL22). Contrarily, models such as Sharma et al. (2016) and Naidu et al. (2020) give an increasing escape fraction towards bright $M_{\text{UV}}$ systems, in disagreement with our picket-fence formulation.<fnref n="8" />

The ionizing production efficiency shows a very smooth (although non-significant) dependence on $M_{1500}$, with the stacked measurement at the faintest UV magnitude bin ($M_{1500} \gtrsim -20$) being higher than the typically assumed canonical value for cosmic reionization of $\log \xi_{\text{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 $M_{1500}$ range is actually not wide enough to show a clear $\xi_{\text{ion}} - M_{1500}$ trend. Even though, it has been shown that, at similar redshifts, the $\xi_{\text{ion}}$ evolution with UV luminosity is a very smooth correlation (see Lam et al. 2019; Emami et al. 2020; Prieto-Lyon et al. 2023). A more complete picture on the $\xi_{\text{ion}} - M_{1500}$ relationship will be given in Section 5.

In general, both $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$ distributions with $M_{1500}$ are characterized by a large scatter (≈1 dex, e.g. 0.1 per cent to 10 per cent in $f_{\text{esc}}^{\text{abs}}$) and a weak evolution with $M_{1500}$, where only the faintest galaxies in the sample have tentatively higher values in $f_{\text{esc}}^{\text{abs}}$ (≥ 5 per cent) and $\log \xi_{\text{ion}}$ (≥ 25.2 Hz erg⁻¹). Finally, the $\log f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ product preserves the overall evolution of $f_{\text{esc}}^{\text{abs}}$ with the UV magnitude.

<footnote n="8">⁸This said, the hypothesis proposed in Naidu et al. (2020) of a reionization driven by $-20 \le M_{1500} \le -18$ galaxies (what they called the ‘oligarchs’), cannot necessarily be ruled out by our observations because our galaxies do not span beyond $M_{\text{UV}} \ge -19$ AB.</footnote>
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 Appendix A, we randomly generated a set of $N = 100$ simulated $\mathrm{Si\texttt{II}}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 ($R(\mathrm{VIMOS}) \approx 600$) and a constant $S/N = 5$ were also incorporated in the mock realizations. After that, we measured the depth of the composite $\mathrm{Si\texttt{II}}1260$ line profile and compared it to the input value given by the covering fraction. According to the our simulations (equation A4), an input ($1 - C_f = 0.1, 0.3, 0.5$, and $0.8$ would require corrections factors as large as 70 per cent, 40 per cent, 25 per cent, and 10 per cent 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 \text{ km s}^{-1}$), i.e. within the actual spectral resolution. Therefore, our simulations – which assume a single VIMOS resolution element ($\pm 150 \text{ km s}^{-1}$) as the standard deviation of the distribution for the velocity shift of the lines – might actually overpredict this effect. For that reason, we prefer not to apply such stacking corrections to our line measurements, but we encourage the reader to check Appendix 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).

# 4 RESULTS

The following paragraphs describe on the main results of this paper. In Section 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 Section 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 Section 4.3, the global ionizing properties of the LCEs versus non-LCEs samples are discussed. Lastly, Section 4.4 summarizes the properties of previously known LCEs reported in the literature which were included in our sample.

## 4.1 The ionizing escape fraction and production efficiency dependence with observed galaxy properties

Fig. 5 investigates how our predicted ionizing absolute escape fraction ($f_{\text{esc}}^{\text{abs}}$) and the ionizing photon production efficiency ($\xi_{\text{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 $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{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 $\times 1.5$) compared to the R16 ones.

A Kendall-$\tau$ statistics was applied to each pair of variables in order to figure out the strength ($\tau$) and significance ($p$-value)

of the correlation. For a sample of $\sim 500$ objects, we consider correlations to be significant if $p_{\text{val.}} \le 10^{-3}$ ($3\sigma$) and strong if $|\tau| \gtrsim 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 $\tau$ coefficients for significant correlations (thick-framed panels), otherwise we write $p_{\text{val.}} > 10^{-3}$.

### 4.1.1 Observed absolute UV magnitude

The predicted $f_{\text{esc}}^{\text{abs}}$ and the observed UV absolute magnitude ($M_{1500}$) for the VANDELS sample at $3 \le z \le 5$ do not show any clear correlation when considering the whole range of UV magnitudes. Tanvir et al. (2019) found no significant correlation of $N_{\text{H}\text{I}}$ 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 equation (6), a lack of correlation between $f_{\text{esc}}^{\text{abs}}$ and $M_{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 $f_{\text{esc}}^{\text{abs}}$ is naturally expected for the faintest galaxies (see Section 7.2 in SL22).

Even though the correlation is tentative, $M_{1500} \gtrsim -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-$z$ by the KLCS (Steidel et al. 2018; Pahl et al. 2021, 2023) as well as observed at low-$z$ by the LzLCS (Chisholm et al. 2022; Flury et al. 2022a, and SL22). Contrarily, models such as Sharma et al. (2016) and Naidu et al. (2020) give an increasing escape fraction towards bright $M_{\text{UV}}$ systems, in disagreement with our picket-fence formulation.<fnref n="8" />

The ionizing production efficiency shows a very smooth (although non-significant) dependence on $M_{1500}$, with the stacked measurement at the faintest UV magnitude bin ($M_{1500} \gtrsim -20$) being higher than the typically assumed canonical value for cosmic reionization of $\log \xi_{\text{ion}}(\text{Hz erg}^{-1}) = 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 $M_{1500}$ range is actually not wide enough to show a clear $\xi_{\text{ion}} - M_{1500}$ trend. Even though, it has been shown that, at similar redshifts, the $\xi_{\text{ion}}$ evolution with UV luminosity is a very smooth correlation (see Lam et al. 2019; Emami et al. 2020; Prieto-Lyon et al. 2023). A more complete picture on the $\xi_{\text{ion}} - M_{1500}$ relationship will be given in Section 5.

In general, both $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$ distributions with $M_{1500}$ are characterized by a large scatter ($\simeq 1$ dex, e.g. 0.1 per cent to 10 per cent in $f_{\text{esc}}^{\text{abs}}$) and a weak evolution with $M_{1500}$, where only the faintest galaxies in the sample have tentatively higher values in $f_{\text{esc}}^{\text{abs}}$ ($\ge 5$ per cent) and $\log \xi_{\text{ion}}$ ($\ge 25.2 \text{ Hz erg}^{-1}$). Finally, the $\log f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ product preserves the overall evolution of $f_{\text{esc}}^{\text{abs}}$ with the UV magnitude.

<footnote n="8">⁸This said, the hypothesis proposed in Naidu et al. (2020) of a reionization driven by $-20 \le M_{1500} \le -18$ galaxies (what they called the ‘oligarchs’), cannot necessarily be ruled out by our observations because our galaxies do not span beyond $M_{\text{UV}} \ge -19$ AB.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Scatter plots searching for correlations between the predicted ionizing photon escape fraction ($f_{\text{esc}}^{\text{abs}}$, first column), the ionizing production efficiency (log $\xi_{\text{ion}}$, second column), and the product of the two (log $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$, third column, see equation 1), versus different galaxy properties: observed and intrinsic UV absolute magnitudes ($M_{1500}$, $M_{1500}^{\text{int}}$), UV-continuum slope at 1500 Å from the best-fitting SED ($\beta_{\text{spec}}^{1500}$), and Ly$\alpha$ equivalent width ($-1 \times W_{\text{Ly}\alpha}$). The Kendall ($\tau$) 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 $M_{1500}^{\text{int}}$ panels measure the shift in $M_{1500}^{\text{int}}$ due to the use of the SMC dust-attenuation law. Grey-shaded regions mark assumed canonical values of $f_{\text{esc}}^{\text{abs}} \ge 5$ per cent and log $\xi_{\text{ion}} = 25.2 - 25.3$ in classical Reionization models (Robertson et al. 2013).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Scatter plots searching for correlations between the predicted ionizing photon escape fraction ($f_{\text{esc}}^{\text{abs}}$, first column), the ionizing production efficiency (log $\xi_{\text{ion}}$, second column), and the product of the two (log $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$, third column, see equation 1), versus different galaxy properties: observed and intrinsic UV absolute magnitudes ($M_{1500}$, $M_{1500}^{\text{int}}$), UV-continuum slope at 1500 Å from the best-fitting SED ($\beta_{\text{spec}}^{1500}$), and Ly$\alpha$ equivalent width ($-1 \times W_{\text{Ly}\alpha}$). The Kendall ($\tau$) 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 $M_{1500}^{\text{int}}$ panels measure the shift in $M_{1500}^{\text{int}}$ due to the use of the SMC dust-attenuation law. Grey-shaded regions mark assumed canonical values of $f_{\text{esc}}^{\text{abs}} \ge 5$ per cent and log $\xi_{\text{ion}} = 25.2 - 25.3$ in classical Reionization models (Robertson et al. 2013).</figcaption></figure>
## 4.1.2 Intrinsic (dust-free) UV luminosity

The intrinsic UV absolute magnitude for each galaxy ($M_{1500}^{int}$) was calculated from the best-fitting 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, $f_{esc}^{abs}$ versus $M_{1500}^{int}$ shows one of the strongest correlations of this study ($\tau \approx 0.4$), where the less intrinsically luminous galaxies have the highest $f_{esc}^{abs}$. From our $f_{esc}^{abs}$ prescription (equation 6), this correlation is expected since the dust-attenuation 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 behaviour was previously reported in Begley et al. (2022), where they demonstrate how intrinsically UV-faint galaxies at $z \simeq 3.5$ would require statistically higher escape fractions in order to reproduce the observed distribution of the ionizing-to-non-ionizing flux ratio.

On the other hand, the $\xi_{ion}$ versus $M_{1500}^{int}$ distribution appears completely flat irrespective of the $M_{1500}^{int}$ bin and the dust-attenuation law. According to equation (4), the dust correction factor applied to any $\xi_{ion}$ measurement would be: $A_{<912}/A_{1500}$. This ratio does not strongly depend on the dust law, thus no $\xi_{ion}$ dependence on the attenuation will be expected neither. Therefore, the $f_{esc}^{abs} \times \xi_{ion}$ product inherits the same dependence on the UV luminosity as $f_{esc}^{abs}$ does. When comparing the R16 and SMC stacks predictions (see Fig. 5), their $f_{esc}^{abs} \times \xi_{ion}$ show the largest differences at the bright $M_{1500}^{int}$ end i.e. the more deviation appears for the more attenuated, redder galaxies.

## 4.1.3 UV-continuum slope at 1500 Å

Standing out as the strongest correlation of this study ($\tau \approx -0.6$, see Fig. 5), the $f_{esc}^{abs}$ inversely scales with the best-fitting UV-continuum slope at 1500 Å ($\beta_{spec}^{1500}$) 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 equation (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 $f_{esc}^{abs}$ 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 $f_{esc}^{abs}$ with increasing UV colours was shown in Pahl et al. (2021) using KLCS galaxy composites at high-$z$. Also interestingly, Begley et al. (2022) found that UV blue galaxies at $z \simeq 3.5$ would require statistically higher escape fractions in order to reproduce the observed ionizing-to-non-ionizing flux ratio distribution.

Regarding our indirect approach to predict $f_{esc}^{abs}$ (equation 6), a negative correlation between the escape fraction and $\beta_{spec}^{1500}$ must be expected by definition, since the UV slope is inherently linked to the dust-attenuation ($E_{B-V}$), so that redder SEDs usually means more attenuated galaxies (Meurer et al. 1999). For example, we obtain $f_{esc}^{abs} \ge 5$ per cent for galaxies whose $\beta_{spec}^{1500} \le -1.5$, decreasing to $f_{esc}^{abs} \approx 0.1$ per cent at $\beta_{spec}^{1500} \approx -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 $f_{esc}^{abs}$ at the EoR (see Trebitsch et al. 2022), and (2) primordial/f 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 Section 6, our $f_{esc}^{abs} - \beta_{spec}^{1500}$ results agree with the relation $z = 0.3$ published by Chisholm et al. (2022) but extrapolated to redder UV slopes and lower escape fractions. This similarity suggests that the $f_{esc}^{abs} - \beta_{spec}^{1500}$ does not strongly change with redshift.

The $\xi_{ion}$ parameter is also strongly correlated with $\beta_{spec}^{1500}$, so that $\xi_{ion}$ rapidly increases as the UV slope decreases, only giving $\log \xi_{ion}(\text{Hz erg}^{-1}) \ge 25.2$ for galaxies whose slopes are $\beta_{spec}^{1500} \le -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 $\beta_{1500}$ values. For more details on the $\xi_{ion} - \beta_{1500}$ relationship and comparison with the literature, we refer to Section 5.

Qualitatively, the largest differences between the R16 and SMC derived $f_{esc}^{abs}$ values are found among the redder objects, i.e. the most attenuated galaxies in the sample. However, both R16 and SMC derived $\xi_{ion}$ values are similar irrespective of the dust-attenuation law.

## 4.1.4 Ly$\alpha$ equivalent width

The relation between the properties of the Ly$\alpha$ 1216 line and the escape and production efficiency of ionizing photons is also striking. This can be seen in Fig. 5, where the predicted $f_{esc}^{abs}$ monotonically raises as the Ly$\alpha$ equivalent width also increases ($W_{Ly\alpha}$). For typically assumed values of $W_{Ly\alpha} \le -20$ Å in LAEs (e.g. Pentericci et al. 2009, but see Stark et al. 2011; Kusakabe et al. 2020), we obtain $f_{esc}^{abs} \ge 5$ per cent. 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$\alpha$ 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$\alpha$ photons and therefore the shape and strength of the Ly$\alpha$ 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 equation (6) via the covering fraction – and (2) the dust-attenuation term ($E_{B-V}$). In this way the Ly$\alpha$ emission ($W_{Ly\alpha}$) and the escape of LyC photons ($f_{esc}^{abs}$) 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$\alpha$ 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$\alpha$ and LyC photons, in Fig. 6 we plot the average LIS absorption profile for different stacks in bins of Ly$\alpha$ equivalent width (see also Trainor et al. 2019). This plot shows the combined profile of several ISM LIS lines (Si II 1260,
## 4.1.2 Intrinsic (dust-free) UV luminosity

The intrinsic UV absolute magnitude for each galaxy ($M_{1500}^{int}$) was calculated from the best-fitting 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, $f_{esc}^{abs}$ versus $M_{1500}^{int}$ shows one of the strongest correlations of this study ($\tau \approx 0.4$), where the less intrinsically luminous galaxies have the highest $f_{esc}^{abs}$. From our $f_{esc}^{abs}$ prescription (equation 6), this correlation is expected since the dust-attenuation 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 behaviour was previously reported in Begley et al. (2022), where they demonstrate how intrinsically UV-faint galaxies at $z \simeq 3.5$ would require statistically higher escape fractions in order to reproduce the observed distribution of the ionizing-to-non-ionizing flux ratio.

On the other hand, the $\xi_{ion}$ versus $M_{1500}^{int}$ distribution appears completely flat irrespective of the $M_{1500}^{int}$ bin and the dust-attenuation law. According to equation (4), the dust correction factor applied to any $\xi_{ion}$ measurement would be: $A_{<912}/A_{1500}$. This ratio does not strongly depend on the dust law, thus no $\xi_{ion}$ dependence on the attenuation will be expected neither. Therefore, the $f_{esc}^{abs} \times \xi_{ion}$ product inherits the same dependence on the UV luminosity as $f_{esc}^{abs}$ does. When comparing the R16 and SMC stacks predictions (see Fig. 5), their $f_{esc}^{abs} \times \xi_{ion}$ show the largest differences at the bright $M_{1500}^{int}$ end i.e. the more deviation appears for the more attenuated, redder galaxies.

## 4.1.3 UV-continuum slope at 1500 Å

Standing out as the strongest correlation of this study ($\tau \approx -0.6$, see Fig. 5), the $f_{esc}^{abs}$ inversely scales with the best-fitting UV-continuum slope at 1500 Å ($\beta_{spec}^{1500}$) 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 equation (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 $f_{esc}^{abs}$ 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 $f_{esc}^{abs}$ with increasing UV colours was shown in Pahl et al. (2021) using KLCS galaxy composites at high-$z$. Also interestingly, Begley et al. (2022) found that UV blue galaxies at $z \simeq 3.5$ would require statistically higher escape fractions in order to reproduce the observed ionizing-to-non-ionizing flux ratio distribution.

Regarding our indirect approach to predict $f_{esc}^{abs}$ (equation 6), a negative correlation between the escape fraction and $\beta_{spec}^{1500}$ must be expected by definition, since the UV slope is inherently linked to the dust-attenuation ($E_{B-V}$), so that redder SEDs usually means more attenuated galaxies (Meurer et al. 1999). For example, we obtain $f_{esc}^{abs} \ge 5$ per cent for galaxies whose $\beta_{spec}^{1500} \le -1.5$, decreasing to $f_{esc}^{abs} \approx 0.1$ per cent at $\beta_{spec}^{1500} \approx -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 $f_{esc}^{abs}$ at the EoR (see Trebitsch et al. 2022), and (2) primordial/f 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 Section 6, our $f_{esc}^{abs} - \beta_{spec}^{1500}$ results agree with the relation $z = 0.3$ published by Chisholm et al. (2022) but extrapolated to redder UV slopes and lower escape fractions. This similarity suggests that the $f_{esc}^{abs} - \beta_{spec}^{1500}$ does not strongly change with redshift.

The $\xi_{ion}$ parameter is also strongly correlated with $\beta_{spec}^{1500}$, so that $\xi_{ion}$ rapidly increases as the UV slope decreases, only giving $\log \xi_{ion}(\text{Hz erg}^{-1}) \ge 25.2$ for galaxies whose slopes are $\beta_{spec}^{1500} \le -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 $\beta_{1500}$ values. For more details on the $\xi_{ion} - \beta_{1500}$ relationship and comparison with the literature, we refer to Section 5.

Qualitatively, the largest differences between the R16 and SMC derived $f_{esc}^{abs}$ values are found among the redder objects, i.e. the most attenuated galaxies in the sample. However, both R16 and SMC derived $\xi_{ion}$ values are similar irrespective of the dust-attenuation law.

## 4.1.4 Ly$\alpha$ equivalent width

The relation between the properties of the Ly$\alpha$ 1216 line and the escape and production efficiency of ionizing photons is also striking. This can be seen in Fig. 5, where the predicted $f_{esc}^{abs}$ monotonically raises as the Ly$\alpha$ equivalent width also increases ($W_{Ly\alpha}$). For typically assumed values of $W_{Ly\alpha} \le -20$ Å in LAEs (e.g. Pentericci et al. 2009, but see Stark et al. 2011; Kusakabe et al. 2020), we obtain $f_{esc}^{abs} \ge 5$ per cent. 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$\alpha$ 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$\alpha$ photons and therefore the shape and strength of the Ly$\alpha$ 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 equation (6) via the covering fraction – and (2) the dust-attenuation term ($E_{B-V}$). In this way the Ly$\alpha$ emission ($W_{Ly\alpha}$) and the escape of LyC photons ($f_{esc}^{abs}$) 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$\alpha$ 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$\alpha$ and LyC photons, in Fig. 6 we plot the average LIS absorption profile for different stacks in bins of Ly$\alpha$ equivalent width (see also Trainor et al. 2019). This plot shows the combined profile of several ISM LIS lines (Si II 1260,
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Combined low-ionization-state (LIS) absorption profile for the stacks in bins of Ly α equivalent width ($W_{Lyα}$, see values in Appendix B), normalized by the local continuum and constructed by averaging the profiles of individual LIS absorption lines (Si II1260, O I+Si II1302, C II1334, Si II1527). The composites are colour coded by the dust-attenuation parameter ($E_{B-V}$) derived from the FICUS SED fits, and the predicted escape fraction for each composite ($f_{esc}^{abs}$) is indicated in the legend. The insets show the resulting continuum-normalized Ly α profile for the stacked spectra (left-hand panel), and the evolution of the Ly α equivalent width as function of the LIS equivalent width in absorption (right-hand panel).</figcaption></figure>

O I+Si II1302, C II1334, Si II1527), colour-coded by the inferred UV dust-attenuation ($E_{B-V}$) and the predicted escape fraction for each composite ($f_{esc}^{abs}$) is indicated in the legend. The resulting Ly α profiles can also be seen in the inset, and finally $W_{Lyα}$ (in Å) is plotted against the LIS equivalent width of the combined LIS absorption profiles ($W_{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 $W_{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 $z=2-5$ (Shapley et al. 2003; Jones, Stark & Ellis 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 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 $W_{Lyα}$ can also be ascribed to a boost in the production of ionizing photons (Nakajima et al. 2018a). The VANDELS high-z $W_{Lyα} - \xi_{ion}$ relation can be found in the bottom-mid panel of Fig. 5, where a significantly higher log $\xi_{ion}$ (Hz erg⁻¹) ≥ 25.2 is found for the strongest LAEs only, with a $W_{Lyα} \approx -30$ Å. Thanks to JWST data in combination with HST photometry, an increasing log $\xi_{ion}$ with Ly α strength has also been shown for individual galaxies at the late-edge of the EoR (Ning et al. 2023). However, in the works by Cullen et al. (2020) and Reddy et al. (2022), $\xi_{ion}$ has been demonstrated to be insufficient to solely account for the whole variation in $W_{Lyα}$,

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Relation between our derived $f_{esc}^{abs}$ and the stellar mass ($M_*$, in $M_⊙$). The $f_{esc}^{abs} - M_*$ 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.</figcaption></figure>

rather ‘the covering fraction of optically thick H I (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 $f_{esc}^{abs} \times \xi_{ion} - W_{Lyα}$ relationship results from the convolution of the $f_{esc}^{abs}$ and $\xi_{ion}$ dependence on the Ly α strength, where the $f_{esc}^{abs} - W_{Lyα}$ behaviour plays the major role. The higher $\xi_{ion}$ values for the objects with the highest Ly α equivalent width in the sample (LAEs) naturally increase the $f_{esc}^{abs} \times \xi_{ion}$ product (Matthee et al. 2022).

### 4.1.5 Stellar mass

The predicted $f_{esc}^{abs}$ strongly scales with the stellar mass ($M_*$) of the VANDELS galaxies so that low-mass galaxies have statistically higher escape fractions than more massive systems (Fig. 7). For instance, only $log M_*(M_⊙) < 9$ VANDELS galaxies show cosmologically relevant escape fractions of $f_{esc}^{abs} \ge 5$ per cent (Robertson et al. 2013, see grey-shaded area in Fig. 5), and these are actually ≃1 dex higher than the average escape fraction of the most massive galaxies in the sample ($log M_*(M_⊙) \ge 10$).

According to our picket-fence model, a physical $f_{esc}^{abs} - M_*$ connection is expected (equation 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 M_*$ increases, the dust-attenuation increases and the line-depth decreases accordingly, and both yield progressively decreasing escape fractions ($f_{esc}^{abs}$). This is consistent with the results by Reddy et al. (2022) where, using a joint modelling of the composite FUV and optical spectra of high-z SFGs, they demonstrate how the H I 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-z there does not seem to be any clear relation between these two quantities (Izotov et al. 2021; Flury
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Combined low-ionization-state (LIS) absorption profile for the stacks in bins of Ly α equivalent width ($W_{Lyα}$, see values in Appendix B), normalized by the local continuum and constructed by averaging the profiles of individual LIS absorption lines (Si II1260, O I+Si II1302, C II1334, Si II1527). The composites are colour coded by the dust-attenuation parameter ($E_{B-V}$) derived from the FICUS SED fits, and the predicted escape fraction for each composite ($f_{esc}^{abs}$) is indicated in the legend. The insets show the resulting continuum-normalized Ly α profile for the stacked spectra (left-hand panel), and the evolution of the Ly α equivalent width as function of the LIS equivalent width in absorption (right-hand panel).</figcaption></figure>

O I+Si II1302, C II1334, Si II1527), colour-coded by the inferred UV dust-attenuation ($E_{B-V}$) and the predicted escape fraction for each composite ($f_{esc}^{abs}$) is indicated in the legend. The resulting Ly α profiles can also be seen in the inset, and finally $W_{Lyα}$ (in Å) is plotted against the LIS equivalent width of the combined LIS absorption profiles ($W_{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 $W_{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 $z=2-5$ (Shapley et al. 2003; Jones, Stark & Ellis 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 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 $W_{Lyα}$ can also be ascribed to a boost in the production of ionizing photons (Nakajima et al. 2018a). The VANDELS high-z $W_{Lyα} - \xi_{ion}$ relation can be found in the bottom-mid panel of Fig. 5, where a significantly higher log $\xi_{ion}$ (Hz erg⁻¹) ≥ 25.2 is found for the strongest LAEs only, with a $W_{Lyα} \approx -30$ Å. Thanks to JWST data in combination with HST photometry, an increasing log $\xi_{ion}$ with Ly α strength has also been shown for individual galaxies at the late-edge of the EoR (Ning et al. 2023). However, in the works by Cullen et al. (2020) and Reddy et al. (2022), $\xi_{ion}$ has been demonstrated to be insufficient to solely account for the whole variation in $W_{Lyα}$,

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Relation between our derived $f_{esc}^{abs}$ and the stellar mass ($M_*$, in $M_⊙$). The $f_{esc}^{abs} - M_*$ 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.</figcaption></figure>

rather ‘the covering fraction of optically thick H I (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 $f_{esc}^{abs} \times \xi_{ion} - W_{Lyα}$ relationship results from the convolution of the $f_{esc}^{abs}$ and $\xi_{ion}$ dependence on the Ly α strength, where the $f_{esc}^{abs} - W_{Lyα}$ behaviour plays the major role. The higher $\xi_{ion}$ values for the objects with the highest Ly α equivalent width in the sample (LAEs) naturally increase the $f_{esc}^{abs} \times \xi_{ion}$ product (Matthee et al. 2022).

### 4.1.5 Stellar mass

The predicted $f_{esc}^{abs}$ strongly scales with the stellar mass ($M_*$) of the VANDELS galaxies so that low-mass galaxies have statistically higher escape fractions than more massive systems (Fig. 7). For instance, only $log M_*(M_⊙) < 9$ VANDELS galaxies show cosmologically relevant escape fractions of $f_{esc}^{abs} \ge 5$ per cent (Robertson et al. 2013, see grey-shaded area in Fig. 5), and these are actually ≃1 dex higher than the average escape fraction of the most massive galaxies in the sample ($log M_*(M_⊙) \ge 10$).

According to our picket-fence model, a physical $f_{esc}^{abs} - M_*$ connection is expected (equation 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 M_*$ increases, the dust-attenuation increases and the line-depth decreases accordingly, and both yield progressively decreasing escape fractions ($f_{esc}^{abs}$). This is consistent with the results by Reddy et al. (2022) where, using a joint modelling of the composite FUV and optical spectra of high-z SFGs, they demonstrate how the H I 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-z there does not seem to be any clear relation between these two quantities (Izotov et al. 2021; Flury
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Composite spectra (main plot) for LCEs (blue) and non-LCEs samples (orange) normalized at 1360 Å (grey band). The position of the main stellar-wind lines (Nv1240, Ov1371, Siiv1400 doublet, Civ1550), nebular lines (Lyα1216, He II1640, O iii]1666, C iii]1908) and ISM lines (Si ii1260, O i+Si ii1302, C iii1334, Si ii1527, Fe ii1608, Al ii1670, and Al iii1858 doublet) are marked in dark-blue, gold, and black labels, respectively. The $W_{Lyα}$ (in Å) and $β_{spec}^{1500}$ values are also indicated in the inset. The top panel shows the number of objects included in each composite as a function of wavelength, and the bottom panel zooms in some of the lines, highlighting fundamental differences in their non-ionizing spectra (equivalent width for each line are listed, see the text).</figcaption></figure>

et al. 2022b), even when spanning a wide range of stellar masses, a gradual increase in $f_{esc}^{abs}$ towards lower masses have been recently reported at high-z in Fletcher et al. (2019), Saxena et al. (2022b), for individual detections.

Lately, Pahl et al. (2023) have also shown a negative and significant $f_{esc}^{abs} - \log M_{\star}$ trend in the $9 \le \log(M_{\star}/M_{\odot}) \le 10$ stellar mass range, using $z \simeq 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 $M_{\star}$ is at best a secondary indicator of the average $f_{esc}^{abs}$. In the $8 \le \log(M_{\star}/M_{\odot}) \le 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), semi-analytical simulations by Seiler et al. (2018) show that they could indirectly influence reionization by clearing out holes in the ionized regions, so that $f_{esc}^{abs}$ can be boosted after quasar wind events. In such scenario, the mean escape fraction peaks for intermediate-mass galaxies, around

$\log(M_{\star}/M_{\odot}) \approx 8$. Regrettably, our sample is not able to probe this hypothesis, since it does not extend to lower stellar masses.

## 4.2 Non-ionizing rest-UV properties of potential LCEs at $3 \le z \le 5$

Having predicted the ionizing escape fraction ($f_{esc}^{abs}$) for every galaxy allows us to split the sample into potential LCEs, with $f_{esc}^{abs} \ge 0.05$ (64 out of 534 galaxies, i.e. $\approx 10$ per cent of LyC detection rate), and non-LCEs with $f_{esc}^{abs} < 0.05$ (470 galaxies, assuming R16). We then build composite spectra for LCEs and non-LCEs candidates, following the same method as described in Section 3.3. The main properties, fitted SED parameters and UV lines measurements are summarized in Appendix 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.

### 4.2.1 UV-continuum slope at 1500 Å

The UV slope for LCEs and non-LCEs are remarkably different: $\beta_{spec}^{1500} = -2.17 \pm 0.03$ for LCEs versus $-1.23 \pm 0.01$ for non-LCEs composites (see Fig. 8). As seen in the previous section, the UV slope is by definition related to $f_{esc}^{abs}$ in the sense that it constitutes
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Composite spectra (main plot) for LCEs (blue) and non-LCEs samples (orange) normalized at 1360 Å (grey band). The position of the main stellar-wind lines (Nv1240, Ov1371, Siiv1400 doublet, Civ1550), nebular lines (Lyα1216, He II1640, Oiii1666, Ciii1908) and ISM lines (Si ii1260, OI+Si ii1302, C iii1334, Si ii1527, Fe ii1608, Al ii1670, and Al iii1858 doublet) are marked in dark-blue, gold, and black labels, respectively. The $W_{Lyα}$ (in Å) and $β_{spec}^{1500}$ values are also indicated in the inset. The top panel shows the number of objects included in each composite as a function of wavelength, and the bottom panel zooms in some of the lines, highlighting fundamental differences in their non-ionizing spectra (equivalent width for each line are listed, see the text).</figcaption></figure>

et al. 2022b), even when spanning a wide range of stellar masses, a gradual increase in $f_{esc}^{abs}$ towards lower masses have been recently reported at high-z in Fletcher et al. (2019), Saxena et al. (2022b), for individual detections.

Lately, Pahl et al. (2023) have also shown a negative and significant $f_{esc}^{abs} - \log M_{\star}$ trend in the $9 \le \log(M_{\star}/M_{\odot}) \le 10$ stellar mass range, using $z \simeq 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 $M_{\star}$ is at best a secondary indicator of the average $f_{esc}^{abs}$. In the $8 \le \log(M_{\star}/M_{\odot}) \le 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), semi-analytical simulations by Seiler et al. (2018) show that they could indirectly influence reionization by clearing out holes in the ionized regions, so that $f_{esc}^{abs}$ can be boosted after quasar wind events. In such scenario, the mean escape fraction peaks for intermediate-mass galaxies, around

$\log(M_{\star}/M_{\odot}) \approx 8$. Regrettably, our sample is not able to probe this hypothesis, since it does not extend to lower stellar masses.

## 4.2 Non-ionizing rest-UV properties of potential LCEs at $3 \le z \le 5$

Having predicted the ionizing escape fraction ($f_{esc}^{abs}$) for every galaxy allows us to split the sample into potential LCEs, with $f_{esc}^{abs} \ge 0.05$ (64 out of 534 galaxies, i.e. $\approx 10$ per cent of LyC detection rate), and non-LCEs with $f_{esc}^{abs} < 0.05$ (470 galaxies, assuming R16). We then build composite spectra for LCEs and non-LCEs candidates, following the same method as described in Section 3.3. The main properties, fitted SED parameters and UV lines measurements are summarized in Appendix 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.

### 4.2.1 UV-continuum slope at 1500 Å

The UV slope for LCEs and non-LCEs are remarkably different: $\beta_{spec}^{1500} = -2.17 \pm 0.03$ for LCEs versus $-1.23 \pm 0.01$ for non-LCEs composites (see Fig. 8). As seen in the previous section, the UV slope is by definition related to $f_{esc}^{abs}$ in the sense that it constitutes