a proxy for dust attenuation at UV wavelengths (Chisholm et al. 2022), therefore favouring a picture in which bluer galaxies with lower levels of dust attenuation ($E_{B-V}$) display higher values of escape fraction (by construction, see equation 6). Having different levels of dust attenuation in LCEs and non-LCEs will partially influence the strength of other nebular emission lines in their UV spectra.

The intrinsically dustier nature of non-LCEs in opposite to the LCEs population has been explicitly reported in SL22 and further investigated in Chisholm et al. (2022) using LzLCS data at $z \simeq 0.3$. Also recently, Pahl et al. (2023) showed a monotonic increase of $f_{esc}^{abs}$ with $E_{B-V}$ from FUV SED fitting in KLCS composite spectra. Likewise, Pahl et al. (2021) and Begley et al. (2022) observationally demonstrated a decrease in $f_{esc}^{abs}$ with increasing UV colours at $z \simeq 3$, using HST (KLCS) and ground-based (VANDELS) LyC imaging, respectively. All these different LyC data sets at low- and high-$z$ will be put together in Section 6.

### 4.2.2 Ly$\alpha$ emission

LCEs clearly show a stronger Ly$\alpha$ emission than non-LCEs (Fig. 8): $W_{Ly\alpha} = -29.71 \pm 2.46$ Å against $-3.62 \pm 0.41$ Å for LCEs and non-LCEs stacks, respectively. Worth mentioning is that the $W_{Ly\alpha}$ for the LCEs composite is compatible with the typical definition of LAEs in the literature (i.e. $W_{Ly\alpha} \le -20$ Å see Pentericci et al. 2009), so that the LCEs composite is mostly composed of the Ly$\alpha$ emitting galaxies in the VANDELS sample. The non-LCEs composite, however, shows little Ly$\alpha$ emission compared to LCEs, with a damped Ly$\alpha$ red-wing extending up to 1240 Å. This suggests a higher column density of H I gas beneath the bulk of the non-LCEs compared to the LCEs population, an unequivocal necessary condition for preventing the leak of Ly$\alpha$ and LyC photons (see Henry et al. 2015). This is also compatible with the strong decrease of the Ly$\alpha$ equivalent width (by a factor $\simeq 9$) between LCEs and non-LCE, which is stronger than the decrease of $\xi_{ion}$ (a factor $\simeq 4$), and which is naturally explained by the impact of a higher column density and higher dust content on the Ly$\alpha$ escape fraction (see e.g. Atek et al. 2014).

Observational evidences for an increase in $W_{Ly\alpha}$ with increasing $f_{esc}^{abs}$ has been reported by a few recent studies: Flury et al. (2022b) through FUV spectroscopy of LzLCS sources at $z \simeq 0.3$; Fletcher et al. (2019) and Pahl et al. (2021, 2023) through LyC imaging of LCEs at $z \simeq 3$ (the LACES and KLCS surveys, respectively); and by Begley et al. (2022), through the statistical measurement of the average $f_{esc}^{abs}$ at $z \simeq 3.5$ (VANDELS).

### 4.2.3 Other nebular emission lines

In Fig. 8, the C IV 1550, He II 1640, and C III 1908 nebular line-profiles are displayed for the LCEs (in blue) and non-LCEs composites (in red). The C IV 1550 and He II 1640 lines for LCEs show slightly stronger emission and narrower profiles, while the equivalent width of the C III 1908 line is remarkably higher in the LCEs stacked spectra compared to the non-LCEs one. We obtain $W_{HeII} = -1.11 \pm 0.25$ Å ($-0.81 \pm 0.06$ Å) and $W_{CIII} = -2.15 \pm 0.67$ Å ($-0.73 \pm 0.13$ Å) for LCEs (non-LCEs).

Similarly, Naidu et al. (2022) stacked the UV spectra of potential LyC emitting candidates according to their Ly$\alpha$ line properties, showing a very different He II 1640 profile respect to the non-emitting candidates (probably attributed to sample selection). Contrarily, in the recent work by Marques-Chaves et al. (2022b) based on direct LyC measurements of low-$z$ LzLCS galaxies, the authors do not find any significant correlation between the LyC escape

fraction and the spectral hardness, where the He II 1640 intensity is fundamentally driven by changes in the metallicity rather than by $f_{esc}^{abs}$. Therefore, the slightly higher observed equivalent widths in He II 1640 and O III 1666 for the LCEs composite (see Llerena et al. 2022) may indicate a much dissimilar gas-phase metallicity between the LCEs and non-LCEs stacks than the stellar metallicities that are actually derived from our stellar population modelling (see Section 4.3).

Regarding the UV carbon nebular lines and based on the work by Schaerer et al. (2022a), Saxena et al. (2022a), and Mascia et al. (2023b), we compute the C IV 1550 over C III 1908 flux ratio, and we obtain C IV/C III] = 0.79 ± 0.21 for LCEs while C IV/C III] = 0.42 ± 0.20 for non-LCEs (the C IV 1550 doublet is blended at the current resolution). Schaerer et al. (2022a) empirically demonstrated that both C IV 1550 and C III 1908 lines tend to be stronger in the LzLCS galaxies with the highest $f_{esc}^{abs}$, and proposed a C IV/C III] ≥ 0.75 threshold for significant LyC leakage which, once again, is in agreement with our stacking analysis. As we will see in the next section, a combined effect of escaping ionizing radiation (favoured neutral gas and dust geometry) with a higher $\xi_{ion}$ parameter is responsible for enhancing the emission of high-ionization-state lines such as C IV 1550 and C III 1908 in LCEs.

### 4.2.4 Absorption lines and ISM properties

Several stellar-wind (N v 1240, Si IV 1400, C IV 1550) and other ISM UV absorption lines (Si II 1260, O I + Si II 1302, C II 1334, Si II 1527, Fe II 1608, Al II 1670, and Al III 1858 doublet) are detected in our VANDELS composite spectra (Fig. 8).

The N v 1240 stellar-wind line in LCEs shows a characteristic P-Cygni profile that can be fully reproduced by stellar templates (excluding potential contamination by AGNs), and it indicates the presence of a young stellar population ($\le 5$ Myr). Conversely, N v 1240 is suppressed by the absorbing part of the damped Ly$\alpha$ red-wing in non-LCEs, indicating a high H I column density. Going to redder wavelengths, the Si IV 1393, 1402 high-ionization state absorption lines, which share both a stellar and ISM origin, have the same equivalent width in the LCEs and non-LCEs composites, suggesting a homogeneous distribution of the diffuse high-ionized species in the ISM of LBGs.

The remaining set of LIS absorption lines in LCEs (Si II 1260, O I + Si II 1302, C II 1334, Si II 1527, Fe II 1608, Al II 1670, and Al III 1858 doublet) clearly show lower equivalent widths than in the non-LCEs stack. For example, the S II 1260 and C II 1334 in Fig. 8, $W_{SiII} = 0.99 \pm 0.13$ (1.82 ± 0.05) for LCEs (non-LCEs), while $W_{CII} = 0.64 \pm 0.13$ (1.74 ± 0.04) in the LCEs (non-LCEs) composite. Following the picket-fence model, the fact that strong LCEs have significantly weaker absorption lines, i.e. lower equivalent widths with higher residual fluxes, can be attributed to a lower covering fraction of the spatially co-existing neutral and low-ionized gas, so that the photons from a given transition escape through low column-density channels in the ISM, as well as LyC photons do (Gazagnes et al. 2018, 2020, SL22).

The observed connection between LIS absorption lines and the strength of the Ly$\alpha$ line is also observed in our galaxy composites, and it can be explained by the picket-fence formulation, in which the gas covering fraction along the line-of-sight primarily governs the escape of the absorbed LIS, Ly$\alpha$, and LyC photons. This behaviour is well summarized in Fig. 9, where the $W_{Ly\alpha}$ is plotted against the measured (averaged) residual flux of the LIS lines, R(LIS), and the points are colour-coded by $E_{B-V}$. Blue and red symbols correspond to our LCEs and non-LCEs galaxy composites, respectively.
a proxy for dust attenuation at UV wavelengths (Chisholm et al. 2022), therefore favouring a picture in which bluer galaxies with lower levels of dust attenuation ($E_{B-V}$) display higher values of escape fraction (by construction, see equation 6). Having different levels of dust attenuation in LCEs and non-LCEs will partially influence the strength of other nebular emission lines in their UV spectra.

The intrinsically dustier nature of non-LCEs in opposite to the LCEs population has been explicitly reported in SL22 and further investigated in Chisholm et al. (2022) using LzLCS data at $z \simeq 0.3$. Also recently, Pahl et al. (2023) showed a monotonic increase of $f_{esc}^{abs}$ with $E_{B-V}$ from FUV SED fitting in KLCS composite spectra. Likewise, Pahl et al. (2021) and Begley et al. (2022) observationally demonstrated a decrease in $f_{esc}^{abs}$ with increasing UV colours at $z \simeq 3$, using HST (KLCS) and ground-based (VANDELS) LyC imaging, respectively. All these different LyC data sets at low- and high-$z$ will be put together in Section 6.

## 4.2.2 Ly$\alpha$ emission

LCEs clearly show a stronger Ly$\alpha$ emission than non-LCEs (Fig. 8): $W_{Ly\alpha} = -29.71 \pm 2.46$ Å against $-3.62 \pm 0.41$ Å for LCEs and non-LCEs stacks, respectively. Worth mentioning is that the $W_{Ly\alpha}$ for the LCEs composite is compatible with the typical definition of LAEs in the literature (i.e. $W_{Ly\alpha} \le -20$ Å see Pentericci et al. 2009), so that the LCEs composite is mostly composed of the Ly$\alpha$ emitting galaxies in the VANDELS sample. The non-LCEs composite, however, shows little Ly$\alpha$ emission compared to LCEs, with a damped Ly$\alpha$ red-wing extending up to 1240 Å. This suggests a higher column density of H I gas beneath the bulk of the non-LCEs compared to the LCEs population, an unequivocal necessary condition for preventing the leak of Ly$\alpha$ and LyC photons (see Henry et al. 2015). This is also compatible with the strong decrease of the Ly$\alpha$ equivalent width (by a factor $\simeq 9$) between LCEs and non-LCE, which is stronger than the decrease of $\xi_{ion}$ (a factor $\simeq 4$), and which is naturally explained by the impact of a higher column density and higher dust content on the Ly$\alpha$ escape fraction (see e.g. Atek et al. 2014).

Observational evidences for an increase in $W_{Ly\alpha}$ with increasing $f_{esc}^{abs}$ has been reported by a few recent studies: Flury et al. (2022b) through FUV spectroscopy of LzLCS sources at $z \simeq 0.3$; Fletcher et al. (2019) and Pahl et al. (2021, 2023) through LyC imaging of LCEs at $z \simeq 3$ (the LACES and KLCS surveys, respectively); and by Begley et al. (2022), through the statistical measurement of the average $f_{esc}^{abs}$ at $z \simeq 3.5$ (VANDELS).

## 4.2.3 Other nebular emission lines

In Fig. 8, the C IV 1550, He II 1640, and C III 1908 nebular line-profiles are displayed for the LCEs (in blue) and non-LCEs composites (in red). The C IV 1550 and He II 1640 lines for LCEs show slightly stronger emission and narrower profiles, while the equivalent width of the C III 1908 line is remarkably higher in the LCEs stacked spectra compared to the non-LCEs one. We obtain $W_{HeII} = -1.11 \pm 0.25$ Å ($-0.81 \pm 0.06$ Å) and $W_{CIII} = -2.15 \pm 0.67$ Å ($-0.73 \pm 0.13$ Å) for LCEs (non-LCEs).

Similarly, Naidu et al. (2022) stacked the UV spectra of potential LyC emitting candidates according to their Ly$\alpha$ line properties, showing a very different He II 1640 profile respect to the non-emitting candidates (probably attributed to sample selection). Contrarily, in the recent work by Marques-Chaves et al. (2022b) based on direct LyC measurements of low-$z$ LzLCS galaxies, the authors do not find any significant correlation between the LyC escape

fraction and the spectral hardness, where the He II 1640 intensity is fundamentally driven by changes in the metallicity rather than by $f_{esc}^{abs}$. Therefore, the slightly higher observed equivalent widths in He II 1640 and O III 1666 for the LCEs composite (see Llerena et al. 2022) may indicate a much dissimilar gas-phase metallicity between the LCEs and non-LCEs stacks than the stellar metallicities that are actually derived from our stellar population modelling (see Section 4.3).

Regarding the UV carbon nebular lines and based on the work by Schaerer et al. (2022a), Saxena et al. (2022a), and Mascia et al. (2023b), we compute the C IV 1550 over C III 1908 flux ratio, and we obtain C IV/C III] = 0.79 ± 0.21 for LCEs while C IV/C III] = 0.42 ± 0.20 for non-LCEs (the C IV 1550 doublet is blended at the current resolution). Schaerer et al. (2022a) empirically demonstrated that both C IV 1550 and C III 1908 lines tend to be stronger in the LzLCS galaxies with the highest $f_{esc}^{abs}$, and proposed a C IV/C III] ≥ 0.75 threshold for significant LyC leakage which, once again, is in agreement with our stacking analysis. As we will see in the next section, a combined effect of escaping ionizing radiation (favoured neutral gas and dust geometry) with a higher $\xi_{ion}$ parameter is responsible for enhancing the emission of high-ionization-state lines such as C IV 1550 and C III 1908 in LCEs.

## 4.2.4 Absorption lines and ISM properties

Several stellar-wind (N v 1240, Si IV 1400, C IV 1550) and other ISM UV absorption lines (Si II 1260, O I+Si II 1302, C II 1334, Si II 1527, Fe II 1608, Al II 1670, and Al III 1858 doublet) are detected in our VANDELS composite spectra (Fig. 8).

The N v 1240 stellar-wind line in LCEs shows a characteristic P-Cygni profile that can be fully reproduced by stellar templates (excluding potential contamination by AGNs), and it indicates the presence of a young stellar population ($\le 5$ Myr). Conversely, N v 1240 is suppressed by the absorbing part of the damped Ly$\alpha$ red-wing in non-LCEs, indicating a high H I column density. Going to redder wavelengths, the Si IV 1393, 1402 high-ionization state absorption lines, which share both a stellar and ISM origin, have the same equivalent width in the LCEs and non-LCEs composites, suggesting a homogeneous distribution of the diffuse high-ionized species in the ISM of LBGs.

The remaining set of LIS absorption lines in LCEs (Si II 1260, O I+Si II 1302, C II 1334, Si II 1527, Fe II 1608, Al II 1670, and Al III 1858 doublet) clearly show lower equivalent widths than in the non-LCEs stack. For example, the S II 1260 and C II 1334 in Fig. 8, $W_{SiII} = 0.99 \pm 0.13$ (1.82 ± 0.05) for LCEs (non-LCEs), while $W_{CII} = 0.64 \pm 0.13$ (1.74 ± 0.04) in the LCEs (non-LCEs) composite. Following the picket-fence model, the fact that strong LCEs have significantly weaker absorption lines, i.e. lower equivalent widths with higher residual fluxes, can be attributed to a lower covering fraction of the spatially co-existing neutral and low-ionized gas, so that the photons from a given transition escape through low column-density channels in the ISM, as well as LyC photons do (Gazagnes et al. 2018, 2020, SL22).

The observed connection between LIS absorption lines and the strength of the Ly$\alpha$ line is also observed in our galaxy composites, and it can be explained by the picket-fence formulation, in which the gas covering fraction along the line-of-sight primarily governs the escape of the absorbed LIS, Ly$\alpha$, and LyC photons. This behaviour is well summarized in Fig. 9, where the $W_{Ly\alpha}$ is plotted against the measured (averaged) residual flux of the LIS lines, R(LIS), and the points are colour-coded by $E_{B-V}$. Blue and red symbols correspond to our LCEs and non-LCEs galaxy composites, respectively.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The relation between the Ly α equivalent width (−1 × $W_{Lyα}$, in Å) and the measured mean residual flux of the LIS lines (R(LIS)). The points are colour-coded by the dust-attenuation parameter ($E_{B-V}$, in mag.). The blue and red circles indicate the position of the LCEs and non-LCEs composite spectra, once corrected by the effect of resolution on the depth of the lines (see Appendix A). For comparison, the empty hexagons represent the measurements of LCEs and non-LCEs stacks from the LzLCS (Flury et al. 2022a, $z = 0.2 - 0.4$), and the black crosses show the individual measurements by Jones et al. (2013, $z = 2 - 4$). Arrows point to already confirmed LCEs in our sample: Ion1 by Vanzella et al. (2012) and CDFS012448 by Saxena et al. (2022b).</figcaption></figure>

As long as the line-of-sight covering fraction of gas decreases and the dust-attenuation increases (traced by higher R(LIS) and lower $E_{B-V}$ values), the Ly α equivalent width dramatically increases towards values typically found in LAEs. For being the case, the stronger LAEs, which has been previously suggested in this work to be a proxy for the bulk of the LCEs population, fall in the top right part of the plot.

This trend matches the results by Steidel et al. (2018) in the framework of the KLCS, and echoes the main conclusions by Gazagnes et al. (2020), using individual LyC detected galaxies at low-$z$. In the same Fig. 9, the median results for LCEs and non-LCEs in the LzLCS (Flury et al. 2022a), are indicated by white open hexagons (see also SL22). As we can see, the LzLCS results clearly agree with our VANDELS stacks, which reassure the use of the picket-fence approach for describing the ISM of LCEs. Additionally, the results by Jones et al. (2013) for individual LAEs at $z = 2 - 4$ are shown on top of the VANDELS results, following the general trend.

## 4.3 Ionizing properties of potential LCEs at $3 \le z \le 5$

The ionizing properties of LCEs and non-LCEs can be accessed by extrapolating the non-ionizing best-fitting SED to ionizing wavelengths. In particular, it is of our most interest to check whether potential differences in the ionizing-to-non-ionizing flux ratio, ($F_{\lambda,900}/F_{\lambda,1500}$)int, and production efficiency, $\xi_{ion}$, exist.

First, we input the LCEs and non-LCEs stacked spectra into the FICUS code in order to unveil the average properties of the underlying stellar population of each galaxy sample. Then, we use the best-fitting stellar continua to define different galaxy observables (see Section 3.1). As we can see, the fit to LCEs provides a lower stellar age ($13 \pm 3$ Myr) than for non-LCEs ($19 \pm 1$ Myr), tentatively suggesting younger stellar populations for LCEs. On the other side, stellar metallicities are very similar

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Main FICUS SED fitting results for the LCEs (top) and non-LCEs (bottom) composite spectra: $E_{B-V}$, stellar age, and metallicity (see the legend). The best-fitting SED models for LCEs (blue) and non-LCEs (orange) are plotted on top of each composite (in grey), which is normalized at 1360 Å. The main stellar, nebular, and ISM lines are labelled through dark-blue, gold, and black labels, as in Fig. 8.</figcaption></figure>

for both ($Z_{*} \approx 0.2 Z_{\odot}$), and compatible with the median of the total sample. These values can be seen in Fig. 10, where the best-fitting FICUS SEDs are superimposed to the LCEs and non-LCEs spectra.

The inferred dust-attenuation parameter is remarkably dissimilar for LCEs and non-LCEs, being $E_{B-V} \approx 0.07$ and $0.23$ mag., respectively. This is by construction, since the parent sample has been split based on the $f_{esc}^{abs}$ of individual VANDELS galaxies (equation 6), which is explicitly related to $E_{B-V}$. It reflects the expected less dusty nature of LCEs compared to the bulk of the LBG population, as we explained in the previous section.

The intrinsic (i.e. dust-free) ionizing spectra of the resulting best-fitting SED for LCEs and non-LCEs is presented in the inset of Fig. 11 (left-hand panel) through blue and red solid lines, correspondingly. The slight decrease in age and metallicity makes the intrinsic ionizing-to-non-ionizing flux ratio (900-to-1500 Å, see C19) to grow up to ($F_{\lambda,900}/F_{\lambda,1500}$)int = 0.64 in LCEs from 0.47 for non-LCEs. A different (intrinsic) 900-to-1500 Å ratio due to an decrease in age also fosters a rise in the ionizing efficiency to $log \xi_{ion}$ (Hz erg$^{-1}$) = 25.38 for LCEs, while it is 25.18 for the non-LCEs composites. In combination with the covering fraction of neutral gas, the $\xi_{ion}$ parameter is partially responsible for shaping the strength of high-ionization-state lines such as C IV 1550 and C III]1908, and other nebular resonant lines such Ly α.

Ultimately, we infer what the observed AB magnitude at LyC wavelengths would look like for the LCEs and non-LCEs populations, as a guidance for upcoming surveys targeting LyC emitting galaxies at similar FUV magnitudes. According to the fundamental definition of the absolute ionizing escape fraction of galaxies (see Izotov et al. 2016a, SL22), the observed flux close to the Lyman edge ($F_{\lambda,900}$) can be simply predicted by doing: $F_{\lambda,900} = f_{esc}^{abs} \times F_{\lambda,900}^{int}$,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The relation between the Ly α equivalent width (−1 × $W_{Lyα}$, in Å) and the measured mean residual flux of the LIS lines (R(LIS)). The points are colour-coded by the dust-attenuation parameter ($E_{B-V}$, in mag.). The blue and red circles indicate the position of the LCEs and non-LCEs composite spectra, once corrected by the effect of resolution on the depth of the lines (see Appendix A). For comparison, the empty hexagons represent the measurements of LCEs and non-LCEs stacks from the LzLCS (Flury et al. 2022a, $z = 0.2 - 0.4$), and the black crosses show the individual measurements by Jones et al. (2013, $z = 2 - 4$). Arrows point to already confirmed LCEs in our sample: Ion1 by Vanzella et al. (2012) and CDFS012448 by Saxena et al. (2022b).</figcaption></figure>

As long as the line-of-sight covering fraction of gas decreases and the dust-attenuation increases (traced by higher R(LIS) and lower $E_{B-V}$ values), the Ly α equivalent width dramatically increases towards values typically found in LAEs. For being the case, the stronger LAEs, which has been previously suggested in this work to be a proxy for the bulk of the LCEs population, fall in the top right part of the plot.

This trend matches the results by Steidel et al. (2018) in the framework of the KLCS, and echoes the main conclusions by Gazagnes et al. (2020), using individual LyC detected galaxies at low-$z$. In the same Fig. 9, the median results for LCEs and non-LCEs in the LzLCS (Flury et al. 2022a), are indicated by white open hexagons (see also SL22). As we can see, the LzLCS results clearly agree with our VANDELS stacks, which reassure the use of the picket-fence approach for describing the ISM of LCEs. Additionally, the results by Jones et al. (2013) for individual LAEs at $z = 2 - 4$ are shown on top of the VANDELS results, following the general trend.

## 4.3 Ionizing properties of potential LCEs at $3 \le z \le 5$

The ionizing properties of LCEs and non-LCEs can be accessed by extrapolating the non-ionizing best-fitting SED to ionizing wavelengths. In particular, it is of our most interest to check whether potential differences in the ionizing-to-non-ionizing flux ratio, ($F_{\lambda,900}/F_{\lambda,1500}$)int, and production efficiency, $\xi_{ion}$, exist.

First, we input the LCEs and non-LCEs stacked spectra into the FICUS code in order to unveil the average properties of the underlying stellar population of each galaxy sample. Then, we use the best-fitting stellar continua to define different galaxy observables (see Section 3.1). As we can see, the fit to LCEs provides a lower stellar age (13 ± 3 Myr) than for non-LCEs (19 ± 1 Myr), tentatively suggesting younger stellar populations for LCEs. On the other side, stellar metallicities are very similar

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Main FICUS SED fitting results for the LCEs (top) and non-LCEs (bottom) composite spectra: $E_{B-V}$, stellar age, and metallicity (see the legend). The best-fitting SED models for LCEs (blue) and non-LCEs (orange) are plotted on top of each composite (in grey), which is normalized at 1360 Å. The main stellar, nebular, and ISM lines are labelled through dark-blue, gold, and black labels, as in Fig. 8.</figcaption></figure>

for both ($Z_{*} \approx 0.2 Z_{\odot}$), and compatible with the median of the total sample. These values can be seen in Fig. 10, where the best-fitting FICUS SEDs are superimposed to the LCEs and non-LCEs spectra.

The inferred dust-attenuation parameter is remarkably dissimilar for LCEs and non-LCEs, being $E_{B-V} \approx 0.07$ and 0.23 mag., respectively. This is by construction, since the parent sample has been split based on the $f_{esc}^{abs}$ of individual VANDELS galaxies (equation 6), which is explicitly related to $E_{B-V}$. It reflects the expected less dusty nature of LCEs compared to the bulk of the LBG population, as we explained in the previous section.

The intrinsic (i.e. dust-free) ionizing spectra of the resulting best-fitting SED for LCEs and non-LCEs is presented in the inset of Fig. 11 (left-hand panel) through blue and red solid lines, correspondingly. The slight decrease in age and metallicity makes the intrinsic ionizing-to-non-ionizing flux ratio (900-to-1500 Å, see C19) to grow up to ($F_{\lambda,900}/F_{\lambda,1500}$)int = 0.64 in LCEs from 0.47 for non-LCEs. A different (intrinsic) 900-to-1500 Å ratio due to an decrease in age also fosters a rise in the ionizing efficiency to $log \xi_{ion}$ (Hz erg⁻¹) = 25.38 for LCEs, while it is 25.18 for the non-LCEs composites. In combination with the covering fraction of neutral gas, the $\xi_{ion}$ parameter is partially responsible for shaping the strength of high-ionization-state lines such as C IV 1550 and C III]1908, and other nebular resonant lines such Ly α.

Ultimately, we infer what the observed AB magnitude at LyC wavelengths would look like for the LCEs and non-LCEs populations, as a guidance for upcoming surveys targeting LyC emitting galaxies at similar FUV magnitudes. According to the fundamental definition of the absolute ionizing escape fraction of galaxies (see Izotov et al. 2016a, SL22), the observed flux close to the Lyman edge ($F_{\lambda,900}$) can be simply predicted by doing: $F_{\lambda,900} = f_{esc}^{abs} \times F_{\lambda,900}^{int}$,
<figure><img src="image_13.png" /><figcaption>**Figure 11.** The ionizing continuum of potential LCEs at $3 \le z \le 5$. Left-hand panel: a comparison between the dust-attenuated SED for LCEs (blue) versus non-LCEs (orange) down to ionizing wavelengths, normalized at 1360 Å. The inset shows the intrinsic ionizing spectrum (i.e. dust-free) of LCEs and non-LCEs (in blue and red, respectively). As written in the legend, the intrinsic 900-to-1500 flux ratio (in $F_\lambda$ units) and the ionizing photon production efficiency ($\xi_{ion}$) are higher for LCEs. The effect of dust-attenuation and neutral hydrogen absorption in both ionizing spectra is sketched through the grey and black arrows, where the blue and orange points represent the corresponding observed LyC flux once both absorptions have been accounted (a constant H I cross-section with wavelength has been considered). Right-hand panel: the AB magnitude distribution of the predicted observed flux in the LyC for the LCEs (blue) and non-LCEs (orange) samples, once multiplied by the average IGM transmission at $\langle z \rangle = 4$ of $\langle T_{IGM} \rangle = 0.3$.</figcaption></figure>

where $F_{\lambda_{900}}^{\text{int}}$ is the intrinsic LyC flux or, in other words, the ionizing flux previous to any gas or dust absorption (in $F_\lambda$ units).

Since we have access to the intrinsic ionizing SED through the FICUS fits, and the $f_{esc}^{\text{abs}}$ has been predicted following equation (6), we compute $F_{2,900}$ for every galaxy in the VANDELS sample. We then derive the AB apparent magnitude from the $F_{2,900}$ synthetic flux measurements ($m_{\text{LYC, AB}}$). The resulting distribution can be seen in Fig. 11 (right-hand panel), where the blue (orange) histogram represent the LyC AB distribution for LCEs (non-LCEs), once we manually apply a $\langle T_{\text{IGM}} \rangle = 0.3$ mean IGM transmission factor, i.e. the mean $T_{\text{IGM}}$ at $z = 4$ according to Inoue et al. (2014).

The consecutive effects of neutral H I gas absorption and dust attenuation in the LCEs and non-LCEs ionizing spectra are represented in Fig. 11 (left-hand panel) through the black and grey arrows (assuming a constant H I cross-section with wavelength), whilst the blue and orange open circles denote the LyC flux measured at 870–910 Å. After applying a $\langle T_{\text{IGM}} \rangle = 0.3$ factor and translating into AB magnitudes, these two points would equal the median of the individual LyC AB magnitude distribution shown in the right-hand panel of the same figure.

Clearly, the combined effects of a lower ionizing-to-non-ionizing flux ratio plus lower escape fractions for non-LCEs makes their median output LyC flux dwindle $\simeq 2$ mag down with respect to a hypothetical LCEs population drawn from the same LBG sample. Therefore, a dedicated survey which attempts to detect the bulk of LCEs at $3 \le z \le 5$ would require a sensitivity such that it reaches, at least, $m_{\text{LYC}} \simeq 31$ AB at LyC wavelengths.

Using very deep VIMOS/U-band imaging of the CDFS field (1σ limiting magnitude of 30.4 AB), Begley et al. (2022) estimated 10 to 15 LyC detections at $2\sigma$ for sources whose redshift is $3.35 \le z \le 3.95$, assuming an average $f_{esc}^{\text{abs}} \approx 0.07$. Restricting to the CDFS only, we found seven sources with predicted $m_{\text{LYC}} \le 30$ AB in the same redshift range. The differences can be attributed to our lower average escape fraction compared to the former work. This said, only

two robust LyC emitters have been detected so far in the CDFS field (see following section). The stochasticity of the IGM opacity might reduce our expectations.

## 4.4 Confirmed LCEs in the VANDELS sample

In Fig. 9, we highlighted two LyC emitting galaxies included in our VANDELS sample which have been already published in the literature, in particular within the CDFS field (see Begley et al. 2022): CDFS012448 (Saxena et al. 2022b) and Ion1 (Vanzella et al. 2012; Ji et al. 2020).

While both galaxies share a similar SED dust-attenuation ($E_{B-V} \simeq 0.2$ mag), we measure a much lower residual flux in CDFS012448 (R(LIS) $\simeq 0.4$) than in Ion1 ($\simeq 0.55$), yielding to a predicted escape fraction of $f_{esc}^{\text{abs}} \simeq 3$ per cent versus $\simeq 6$ per cent. It is worth stressing that our inferred $f_{esc}^{\text{abs}}$ value for Ion1 equals the reported value by Ji et al. (2020). In the case of CDFS012448, our $f_{esc}^{\text{abs}}$ is remarkably different to the one quoted in Saxena et al. (2022b) of $\approx 20$ per cent, possibly because the different approaches between papers. Our method was able to identify both galaxies as potential LCEs.

Additionally, these galaxies exhibit very different Ly$\alpha$ properties between each other: while CDSFS012448 shows a relatively strong Ly$\alpha$ in emission ($W_{Ly\alpha} \simeq -20$ Å), the Ion1 Ly$\alpha$ profile appears in absorption ($W_{Ly\alpha} \simeq 1$ Å). In our VANDELS sample, $\simeq 10$ per cent of the LCEs ($f_{esc}^{\text{abs}} \ge 5$ per cent) appear with Ly$\alpha$ in absorption. This manifests, once again, the variety of Ly$\alpha$ profiles that can be found among the LCEs galaxy population (Izotov et al. 2020; Flury et al. 2022b; Naidu et al. 2022) where, given the fact that Ly$\alpha$ can be resonantly scattered whereas LyC cannot be scattered, both Ly$\alpha$ and the ionizing emission do not necessarily come from the same spatial location within the host galaxy (e.g. Vanzella et al. 2012). A low IGM transmission could also, in principle, have killed out most of the Ly$\alpha$ photons of Ion1.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** The ionizing continuum of potential LCEs at $3 \le z \le 5$. Left-hand panel: a comparison between the dust-attenuated SED for LCEs (blue) versus non-LCEs (orange) down to ionizing wavelengths, normalized at 1360 Å. The inset shows the intrinsic ionizing spectrum (i.e. dust-free) of LCEs and non-LCEs (in blue and red, respectively). As written in the legend, the intrinsic 900-to-1500 flux ratio (in $F_\lambda$ units) and the ionizing photon production efficiency ($\xi_{ion}$) are higher for LCEs. The effect of dust-attenuation and neutral hydrogen absorption in both ionizing spectra is sketched through the grey and black arrows, where the blue and orange points represent the corresponding observed LyC flux once both absorptions have been accounted (a constant H I cross-section with wavelength has been considered). Right-hand panel: the AB magnitude distribution of the predicted observed flux in the LyC for the LCEs (blue) and non-LCEs (orange) samples, once multiplied by the average IGM transmission at $\langle z \rangle = 4$ of $\langle T_{IGM} \rangle = 0.3$.</figcaption></figure>

where $F_{\lambda_{900}}^{\text{int}}$ is the intrinsic LyC flux or, in other words, the ionizing flux previous to any gas or dust absorption (in $F_\lambda$ units).

Since we have access to the intrinsic ionizing SED through the FICUS fits, and the $f_{esc}^{\text{abs}}$ has been predicted following equation (6), we compute $F_{2,900}$ for every galaxy in the VANDELS sample. We then derive the AB apparent magnitude from the $F_{2,900}$ synthetic flux measurements ($m_{LyC}$, AB). The resulting distribution can be seen in Fig. 11 (right-hand panel), where the blue (orange) histogram represent the LyC AB distribution for LCEs (non-LCEs), once we manually apply a $\langle T_{IGM} \rangle = 0.3$ mean IGM transmission factor, i.e. the mean $T_{IGM}$ at $z = 4$ according to Inoue et al. (2014).

The consecutive effects of neutral H I gas absorption and dust attenuation in the LCEs and non-LCEs ionizing spectra are represented in Fig. 11 (left-hand panel) through the black and grey arrows (assuming a constant H I cross-section with wavelength), whilst the blue and orange open circles denote the LyC flux measured at 870–910 Å. After applying a $\langle T_{IGM} \rangle = 0.3$ factor and translating into AB magnitudes, these two points would equal the median of the individual LyC AB magnitude distribution shown in the right-hand panel of the same figure.

Clearly, the combined effects of a lower ionizing-to-non-ionizing flux ratio plus lower escape fractions for non-LCEs makes their median output LyC flux dwindle $\simeq 2$ mag down with respect to a hypothetical LCEs population drawn from the same LBG sample. Therefore, a dedicated survey which attempts to detect the bulk of LCEs at $3 \le z \le 5$ would require a sensitivity such that it reaches, at least, $m_{LyC} \simeq 31$ AB at LyC wavelengths.

Using very deep VIMOS/U-band imaging of the CDFS field (1σ limiting magnitude of 30.4 AB), Begley et al. (2022) estimated 10 to 15 LyC detections at $2\sigma$ for sources whose redshift is $3.35 \le z \le 3.95$, assuming an average $f_{esc}^{\text{abs}} \approx 0.07$. Restricting to the CDFS only, we found seven sources with predicted $m_{LyC} \le 30$ AB in the same redshift range. The differences can be attributed to our lower average escape fraction compared to the former work. This said, only

two robust LyC emitters have been detected so far in the CDFS field (see following section). The stochasticity of the IGM opacity might reduce our expectations.

## 4.4 Confirmed LCEs in the VANDELS sample

In Fig. 9, we highlighted two LyC emitting galaxies included in our VANDELS sample which have been already published in the literature, in particular within the CDFS field (see Begley et al. 2022): CDFS012448 (Saxena et al. 2022b) and Ion1 (Vanzella et al. 2012; Ji et al. 2020).

While both galaxies share a similar SED dust-attenuation ($E_{B-V} \simeq 0.2$ mag), we measure a much lower residual flux in CDFS012448 (R(LIS) $\simeq 0.4$) than in Ion1 ($\simeq 0.55$), yielding to a predicted escape fraction of $f_{esc}^{\text{abs}} \simeq 3$ per cent versus $\simeq 6$ per cent. It is worth stressing that our inferred $f_{esc}^{\text{abs}}$ value for Ion1 equals the reported value by Ji et al. (2020). In the case of CDFS012448, our $f_{esc}^{\text{abs}}$ is remarkably different to the one quoted in Saxena et al. (2022b) of $\approx 20$ per cent, possibly because the different approaches between papers. Our method was able to identify both galaxies as potential LCEs.

Additionally, these galaxies exhibit very different Ly$\alpha$ properties between each other: while CDSFS012448 shows a relatively strong Ly$\alpha$ in emission ($W_{Ly\alpha} \simeq -20$ Å), the Ion1 Ly$\alpha$ profile appears in absorption ($W_{Ly\alpha} \simeq 1$ Å). In our VANDELS sample, $\simeq 10$ per cent of the LCEs ($f_{esc}^{\text{abs}} \ge 5$ per cent) appear with Ly$\alpha$ in absorption. This manifests, once again, the variety of Ly$\alpha$ profiles that can be found among the LCEs galaxy population (Izotov et al. 2020; Flury et al. 2022b; Naidu et al. 2022) where, given the fact that Ly$\alpha$ can be resonantly scattered whereas LyC cannot be scattered, both Ly$\alpha$ and the ionizing emission do not necessarily come from the same spatial location within the host galaxy (e.g. Vanzella et al. 2012). A low IGM transmission could also, in principle, have killed out most of the Ly$\alpha$ photons of Ion1.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** log $\xi_{\text{ion}}$ as a function of the observed absolute UV magnitude. The red thick symbols show the results of our $3 \le z \le 5$ VANDELS composite spectra in bins of UV magnitude when using the R16 attenuation law (SMC dust-law results are shown in the inset). Stacked results from other literature samples at various redshifts are represented via open white symbols (see the legend), in particular for low-$z$ LCEs (Flury et al. 2022a), high-$z$ LAEs (Matthee et al. 2017b; Harikane et al. 2018; Nakajima et al. 2018a; Maseda et al. 2020), and normal high-$z$ SFGs (Nakajima et al. 2018b; Lam et al. 2019). Our main comparison sources (studies targeting LBGs at UV magnitudes similar to the VANDELS galaxies) are the MOSDEF galaxies at $z = 2-3$ (Shivaei et al. 2018, orange squares), the extensive HST+IRAC campaigns at $z = 4-5$ by Bouwens et al. (2016, blue diamonds), and the HAEs sample at $z = 2-3$ by Matthee et al. (2017a, grey crosses). The individual LAEs and C IV emitters measurements from Nakajima et al. (2016) at $z = 3-4$ and Stark et al. (2015) at $z = 7$ are plotted in coloured circles. The teal dashed line follows the extrapolation to brighter magnitudes of the recent HST+JWST results at $z = 3-7$ (Prieto-Lyon et al. 2023). The grey-shaded area (Maseda et al. 2020) marks the canonical log $\xi_{\text{ion}} = 25.2-25.3$ value given by a simple stellar population at constant SFR over 100 Myr (Robertson et al. 2013). For reference, $M_{\text{UV}}^{*} \approx -21$ at $z = 6$ (Bouwens et al. 2021). Overall, log $\xi_{\text{ion}}$ monotonically increases with $M_{1500}$.</figcaption></figure>

# 5 THE PRODUCTION EFFICIENCY OF IONIZING PHOTONS ($\xi_{\text{ION}}$) IN HIGH-REDSHIFT GALAXIES

Here, we investigate the evolution of $\xi_{\text{ion}}$ as a function of the UV absolute magnitude and the UV slope (see Bouwens et al. 2016; Chisholm et al. 2022), quantities which will be easily measured for a significant fraction of galaxies at the EoR thanks to upcoming ground-based (ELT, GMT, TMT,...) and ongoing space-based facilities (JWST).

Fig. 12 shows the $\xi_{\text{ion}}$ parameter in stacks of the absolute UV magnitude for the VANDELS sample. As it was introduced in Section 4 and regardless of the attenuation law (see the inset), the ionizing production efficiency monotonically, smoothly increases towards fainter UV magnitudes, so that log $\xi_{\text{ion}} \ge 25.2 \text{ Hz erg}^{-1}$ (Robertson et al. 2013) at $M_{1500} \gtrsim -20$, i.e. the faint-end of our sample. Stacked results from other samples in the literature are also included in Fig. 12. The main comparison samples for our case of study are the MOSDEF galaxies from Shivaei et al. (2018) at $z = 2-3$, the extensive HST+IRAC searches by Bouwens et al. (2016) at $z = 4-5$, and the HAEs sample of Matthee et al. (2017a) at $z = 2-3$. Our results are qualitatively in agreement with Bouwens et al. (2016) and Matthee et al. (2017a), both tracing photometrically selected HAEs at similar $M_{1500}$ than ours. Nevertheless, VANDELS samples a narrower range in $M_{1500}$ so that the overall trend does not appear as clear as in the former works. At the faintest UV magnitude bin, our $\xi_{\text{ion}} - M_{1500}$ relationship disagrees with the results by Shivaei et al. (2018) who report no log $\xi_{\text{ion}}$ evolution with $M_{1500}$, a result probably induced by the different dust-attenuation corrections between studies

(Balmer Decrement for MOSDEF versus global SED for the rest, see the discussion in Matthee et al. 2017a).

Due to their lower metallicities and burstiness of their SFHs, fainter (low-mass) galaxies are expected to have high ionizing production efficiencies (Ma et al. 2020). If the LyC escape fraction and production efficiency were intertwined (see below), the fact that fainter galaxies produce more ionizing photons per UV luminosity which, at the same time, would escape more easily than in higher mass counterparts, supports an early and slow reionization scenario in which the ionizing budget would be dominated by the same low-mass, UV-faint galaxies (see Finkelstein et al. 2019; Trebitsch et al. 2022), thanks to their higher number density at the EoR (Bouwens et al. 2015, 2021).

Recently, Prieto-Lyon et al. (2023) used deep HST and JWST multiband photometry to estimate the H$\alpha$ flux and thus the ionizing production efficiency extending to very faint UV magnitudes ($-23 \le M_{1500} \le -15.5$) and spanning a wide range in redshift ($z = 3-7$). Their results demonstrate that, although the general trend of increasing $\xi_{\text{ion}}$ with increasing UV magnitude is very smooth, it holds at any $M_{1500}$ (see also Lam et al. 2019; Emami et al. 2020). Additionally, our thorough compilation of $\xi_{\text{ion}}$ measurements (Fig. 12) illustrates the scatter in this relation ($\sim 1$ dex) and its dependence on the galaxy type, physical properties of the underlying galaxy-population and methodology itself (with the dust corrections among the different samples playing a critical role in the scatter). While normal high-$z$ SFGs (Nakajima et al. 2018b; Lam et al. 2019) have modest production efficiencies, high-$z$ LAEs (Matthee et al. 2017b; Harikane et al. 2018; Nakajima et al. 2018a) and C III emitters
<figure><img src="image_14.png" /><figcaption>**Figure 12.** log $\xi_{\text{ion}}$ as a function of the observed absolute UV magnitude. The red thick symbols show the results of our $3 \le z \le 5$ VANDELS composite spectra in bins of UV magnitude when using the R16 attenuation law (SMC dust-law results are shown in the inset). Stacked results from other literature samples at various redshifts are represented via open white symbols (see the legend), in particular for low-$z$ LCEs (Flury et al. 2022a), high-$z$ LAEs (Matthee et al. 2017b; Harikane et al. 2018; Nakajima et al. 2018a; Maseda et al. 2020), and normal high-$z$ SFGs (Nakajima et al. 2018b; Lam et al. 2019). Our main comparison sources (studies targeting LBGs at UV magnitudes similar to the VANDELS galaxies) are the MOSDEF galaxies at $z = 2-3$ (Shivaei et al. 2018, orange squares), the extensive HST+IRAC campaigns at $z = 4-5$ by Bouwens et al. (2016, blue diamonds), and the HAEs sample at $z = 2-3$ by Matthee et al. (2017a, grey crosses). The individual LAEs and C IV emitters measurements from Nakajima et al. (2016) at $z = 3-4$ and Stark et al. (2015) at $z = 7$ are plotted in coloured circles. The teal dashed line follows the extrapolation to brighter magnitudes of the recent HST+JWST results at $z = 3-7$ (Prieto-Lyon et al. 2023). The grey-shaded area (Maseda et al. 2020) marks the canonical log $\xi_{\text{ion}} = 25.2-25.3$ value given by a simple stellar population at constant SFR over 100 Myr (Robertson et al. 2013). For reference, $M_{\text{UV}}^{*} \approx -21$ at $z = 6$ (Bouwens et al. 2021). Overall, log $\xi_{\text{ion}}$ monotonically increases with $M_{1500}$.</figcaption></figure>

## 5 THE PRODUCTION EFFICIENCY OF IONIZING PHOTONS ($\xi_{\text{ION}}$) IN HIGH-REDSHIFT GALAXIES

Here, we investigate the evolution of $\xi_{\text{ion}}$ as a function of the UV absolute magnitude and the UV slope (see Bouwens et al. 2016; Chisholm et al. 2022), quantities which will be easily measured for a significant fraction of galaxies at the EoR thanks to upcoming ground-based (ELT, GMT, TMT,...) and ongoing space-based facilities (JWST).

Fig. 12 shows the $\xi_{\text{ion}}$ parameter in stacks of the absolute UV magnitude for the VANDELS sample. As it was introduced in Section 4 and regardless of the attenuation law (see the inset), the ionizing production efficiency monotonically, smoothly increases towards fainter UV magnitudes, so that log $\xi_{\text{ion}} \ge 25.2 \text{ Hz erg}^{-1}$ (Robertson et al. 2013) at $M_{1500} \gtrsim -20$, i.e. the faint-end of our sample. Stacked results from other samples in the literature are also included in Fig. 12. The main comparison samples for our case of study are the MOSDEF galaxies from Shivaei et al. (2018) at $z = 2-3$, the extensive HST+IRAC searches by Bouwens et al. (2016) at $z = 4-5$, and the HAEs sample of Matthee et al. (2017a) at $z = 2-3$. Our results are qualitatively in agreement with Bouwens et al. (2016) and Matthee et al. (2017a), both tracing photometrically selected HAEs at similar $M_{1500}$ than ours. Nevertheless, VANDELS samples a narrower range in $M_{1500}$ so that the overall trend does not appear as clear as in the former works. At the faintest UV magnitude bin, our $\xi_{\text{ion}} - M_{1500}$ relationship disagrees with the results by Shivaei et al. (2018) who report no log $\xi_{\text{ion}}$ evolution with $M_{1500}$, a result probably induced by the different dust-attenuation corrections between studies

(Balmer Decrement for MOSDEF versus global SED for the rest, see the discussion in Matthee et al. 2017a).

Due to their lower metallicities and burstiness of their SFHs, fainter (low-mass) galaxies are expected to have high ionizing production efficiencies (Ma et al. 2020). If the LyC escape fraction and production efficiency were intertwined (see below), the fact that fainter galaxies produce more ionizing photons per UV luminosity which, at the same time, would escape more easily than in higher mass counterparts, supports an early and slow reionization scenario in which the ionizing budget would be dominated by the same low-mass, UV-faint galaxies (see Finkelstein et al. 2019; Trebitsch et al. 2022), thanks to their higher number density at the EoR (Bouwens et al. 2015, 2021).

Recently, Prieto-Lyon et al. (2023) used deep HST and JWST multiband photometry to estimate the H$\alpha$ flux and thus the ionizing production efficiency extending to very faint UV magnitudes ($-23 \le M_{1500} \le -15.5$) and spanning a wide range in redshift ($z = 3-7$). Their results demonstrate that, although the general trend of increasing $\xi_{\text{ion}}$ with increasing UV magnitude is very smooth, it holds at any $M_{1500}$ (see also Lam et al. 2019; Emami et al. 2020). Additionally, our thorough compilation of $\xi_{\text{ion}}$ measurements (Fig. 12) illustrates the scatter in this relation ($\sim 1$ dex) and its dependence on the galaxy type, physical properties of the underlying galaxy-population and methodology itself (with the dust corrections among the different samples playing a critical role in the scatter). While normal high-$z$ SFGs (Nakajima et al. 2018b; Lam et al. 2019) have modest production efficiencies, high-$z$ LAEs (Matthee et al. 2017b; Harikane et al. 2018; Nakajima et al. 2018a) and C III emitters
<figure><img src="image_15.png" /><figcaption>**Figure 13.** log $\xi_{\text{ion}}$ versus $M_{1500}$ colour-coded for different galaxy samples. Top: potential LCEs versus non-LCEs (in blue and red, respectively). Bottom: selected LAEs and non-LAEs (in green and light-green). The corresponding measurements for LCEs/non-LCEs and LAEs/non-LAEs composite spectra are represented through thick filled symbols, and the thick black line represents our stacked measurements in bins of UV magnitude. Stacked results from other literature samples are represented via open white symbols, in particular for low-z LCEs (Flury et al. 2022a) and high-z LAEs (Matthee et al. 2017b; Harikane et al. 2018; Nakajima et al. 2018a, see Fig. 12 for symbols). The grey-shaded area marks the canonical log $\xi_{\text{ion}} = 25.2 - 25.3$ value given by a simple stellar population at constant SFR over 100 Myr (Robertson et al. 2013), and the number of LCEs and LAEs above and below this limit is indicated in the right-hand side of the plot. In conclusion, LCEs and LAEs have systematically higher $\xi_{\text{ion}}$ than the bulk of the LBG population (see also Reddy et al. 2022).</figcaption></figure>

(Nakajima et al. 2018b) show much higher log $\xi_{\text{ion}}$ values at similar UV magnitudes, since overall these galaxies need stronger ionizing spectra to produce such high excitation lines (see e.g. Reddy et al. 2018).

Intending to qualitatively illustrate this behaviour (Fig. 13), we search for meaningful differences in $\xi_{\text{ion}}$ between the bulk of the VANDELS LBG population and other sources like LCEs or LAEs, which would potentially boost $\xi_{\text{ion}}$ above the canonical value. This is done by splitting the sample into LAEs and non-LAEs, where the sources whose $W_{\text{Ly}\alpha} > -20$ Å were considered as LAEs, following the definition by Pentericci et al. (2009), but see Stark et al. 2011, Kusakabe et al. 2020). A similar analysis as we did for the LCEs/non-

LCEs composites was then performed for the LAEs/non-LAEs composites, and the results are described in the following lines.

Of particular interest is to search for $\xi_{\text{ion}}$ differences within different galaxy subsamples. In Fig. 13 (top), the VANDELS sample is colour-coded by potential LCEs ($f_{\text{esc}}^{\text{abs}} \ge 0.05$, in blue) and non-LCEs ($f_{\text{esc}}^{\text{abs}} < 0.05$, in red), where approximately 60 per cent of the LCEs sample (39 out of 64 LCEs candidates) has log $\xi_{\text{ion}} \ge 25.2$ Hz erg$^{-1}$, and they uniformly populate the whole range in UV magnitude. As said in Section 4.3, the resulting $\xi_{\text{ion}}$ for LCEs and non-LCEs are log $\xi_{\text{ion}}$ (Hz erg$^{-1}) \approx 25.38$ and 25.18, respectively. Compared with other values in the literature, our estimated ionizing efficiency for LCEs is consistent with the average of the LCEs sample of the LzLCS (Chisholm et al. 2022). The fact that our results at $3 \le z \le 5$ agree with low-z LCEs sample supports the argument by which the LzLCS galaxies might be analogues of higher redshift reionizers (Flury et al. 2022b). This has been lately demonstrated by Mascia et al. (2023a) through JWST spectroscopy of 29 moderately faint galaxies at $5 \le z \le 8$, whose optical line ratios, UV slopes, compactness, masses and predicted $f_{\text{esc}}^{\text{abs}}$ were compatible with the strongest LCEs in LzLCS sample (see also Lin et al. 2023).

Similarly, in Fig. 13 (bottom) we explore the $\xi_{\text{ion}} - M_{1500}$ relationship for selected LAEs ($W_{\text{Ly}\alpha} > -20$ Å, in green) and non-LAEs ($W_{\text{Ly}\alpha} \le -20$ Å, in light-green). Although the $M_{1500}$ distribution for both galaxy types equally span along the x-axis of this figure, around 70 per cent of the LAEs sample (75 out of 104 identified LAEs) falls above the log $\xi_{\text{ion}}$ canonical value in the y-axis. Our $\xi_{\text{ion}}$ estimate for LAEs is in perfect agreement with the results by Harikane et al. (2018) at $z = 5 - 6$, and it matches surprisingly well the global trend found by several high-z LAEs surveys at $3 \le z \le 7$ (Matthee et al. 2017b; Harikane et al. 2018; Nakajima et al. 2018a, see also Ning et al. 2023). Similar to Ly $\alpha$, the ionizing efficiency also shapes the emission of other nebular lines (Reddy et al. 2018) so that, for instance, higher-than-average $\xi_{\text{ion}}$ values have been found among extreme high-z [O III] emitters (Tang et al. 2019) and local, high H $\beta$ equivalent width galaxy samples (Izotov et al. 2017).

The escape of Ly $\alpha$ photons and the production of strong UV nebular emission usually requires special ISM and radiation field conditions such as low metallicity, low dust, and neutral gas content, and a vast production of (high-energy) ionizing photons. In Fig. 13 one can visualize a boost in the ionizing production efficiency ($\xi_{\text{ion}}$) for high-z LAEs' samples. At the same time, a rise in the ionizing photon flux may reduce the covering fraction or column density of neutral hydrogen, which may ease the escape of LyC and Ly $\alpha$ photons (see Erb et al. 2014). Recently, Flury et al. (2022b) and Schaerer et al. (2022a) statistically demonstrated that the detection rate of LCEs is enhanced among the LAEs and C IV emitter galaxy-population (see also Saxena et al. 2022a; Mascia et al. 2023b), reiterating our results.

Moving forward with the discussion, Fig. 14 shows $\xi_{\text{ion}}$ as function of the UV-continuum slope for the VANDELS sample when using either the R16 (left-hand panel) or the SMC (right-hand panel) attenuation law. Irrespective of the dust law, the ionizing efficiency rapidly decreases with the slope of the continuum, and the trend is more steep for the SMC law. Focusing on the stacks measurements, the canonical limit in $\xi_{\text{ion}}$ (log $\xi_{\text{ion}} \ge 25.2$ Hz erg$^{-1}$) is only achieved for galaxies whose slopes are bluer than $\beta_{1500} \le -1.5$. The $\xi_{\text{ion}} - \beta_{1500}$ relation of a dustless, single-burst stellar population at increasing age is not able to fully reproduce the whole range of $\beta_{1500}$ values, so that it would require 'some' dust reddening to account for redder slopes (see horizontal arrow in Fig. 14, indicating the change in UV slope
<figure><img src="image_15.png" /><figcaption>**Figure 13.** log $\xi_{\text{ion}}$ versus $M_{1500}$ colour-coded for different galaxy samples. Top: potential LCEs versus non-LCEs (in blue and red, respectively). Bottom: selected LAEs and non-LAEs (in green and light-green). The corresponding measurements for LCEs/non-LCEs and LAEs/non-LAEs composite spectra are represented through thick filled symbols, and the thick black line represents our stacked measurements in bins of UV magnitude. Stacked results from other literature samples are represented via open white symbols, in particular for low-z LCEs (Flury et al. 2022a) and high-z LAEs (Matthee et al. 2017b; Harikane et al. 2018; Nakajima et al. 2018a, see Fig. 12 for symbols). The grey-shaded area marks the canonical log $\xi_{\text{ion}} = 25.2 - 25.3$ value given by a simple stellar population at constant SFR over 100 Myr (Robertson et al. 2013), and the number of LCEs and LAEs above and below this limit is indicated in the right-hand side of the plot. In conclusion, LCEs and LAEs have systematically higher $\xi_{\text{ion}}$ than the bulk of the LBG population (see also Reddy et al. 2022).</figcaption></figure>

(Nakajima et al. 2018b) show much higher log $\xi_{\text{ion}}$ values at similar UV magnitudes, since overall these galaxies need stronger ionizing spectra to produce such high excitation lines (see e.g. Reddy et al. 2018).

Intending to qualitatively illustrate this behaviour (Fig. 13), we search for meaningful differences in $\xi_{\text{ion}}$ between the bulk of the VANDELS LBG population and other sources like LCEs or LAEs, which would potentially boost $\xi_{\text{ion}}$ above the canonical value. This is done by splitting the sample into LAEs and non-LAEs, where the sources whose $W_{\text{Ly}\alpha} > -20$ Å were considered as LAEs, following the definition by Pentericci et al. (2009), but see Stark et al. 2011, Kusakabe et al. 2020). A similar analysis as we did for the LCEs/non-

LCEs composites was then performed for the LAEs/non-LAEs composites, and the results are described in the following lines.

Of particular interest is to search for $\xi_{\text{ion}}$ differences within different galaxy subsamples. In Fig. 13 (top), the VANDELS sample is colour-coded by potential LCEs ($f_{\text{esc}}^{\text{abs}} \ge 0.05$, in blue) and non-LCEs ($f_{\text{esc}}^{\text{abs}} < 0.05$, in red), where approximately 60 per cent of the LCEs sample (39 out of 64 LCEs candidates) has log $\xi_{\text{ion}} \ge 25.2$ Hz erg$^{-1}$, and they uniformly populate the whole range in UV magnitude. As said in Section 4.3, the resulting $\xi_{\text{ion}}$ for LCEs and non-LCEs are log $\xi_{\text{ion}}$ (Hz erg$^{-1}) \approx 25.38$ and 25.18, respectively. Compared with other values in the literature, our estimated ionizing efficiency for LCEs is consistent with the average of the LCEs sample of the LzLCS (Chisholm et al. 2022). The fact that our results at $3 \le z \le 5$ agree with low-z LCEs sample supports the argument by which the LzLCS galaxies might be analogues of higher redshift reionizers (Flury et al. 2022b). This has been lately demonstrated by Mascia et al. (2023a) through JWST spectroscopy of 29 moderately faint galaxies at $5 \le z \le 8$, whose optical line ratios, UV slopes, compactness, masses and predicted $f_{\text{esc}}^{\text{abs}}$ were compatible with the strongest LCEs in LzLCS sample (see also Lin et al. 2023).

Similarly, in Fig. 13 (bottom) we explore the $\xi_{\text{ion}} - M_{1500}$ relationship for selected LAEs ($W_{\text{Ly}\alpha} > -20$ Å, in green) and non-LAEs ($W_{\text{Ly}\alpha} \le -20$ Å, in light-green). Although the $M_{1500}$ distribution for both galaxy types equally span along the x-axis of this figure, around 70 per cent of the LAEs sample (75 out of 104 identified LAEs) falls above the log $\xi_{\text{ion}}$ canonical value in the y-axis. Our $\xi_{\text{ion}}$ estimate for LAEs is in perfect agreement with the results by Harikane et al. (2018) at $z = 5 - 6$, and it matches surprisingly well the global trend found by several high-z LAEs surveys at $3 \le z \le 7$ (Matthee et al. 2017b; Harikane et al. 2018; Nakajima et al. 2018a, see also Ning et al. 2023). Similar to Ly $\alpha$, the ionizing efficiency also shapes the emission of other nebular lines (Reddy et al. 2018) so that, for instance, higher-than-average $\xi_{\text{ion}}$ values have been found among extreme high-z [O III] emitters (Tang et al. 2019) and local, high H $\beta$ equivalent width galaxy samples (Izotov et al. 2017).

The escape of Ly $\alpha$ photons and the production of strong UV nebular emission usually requires special ISM and radiation field conditions such as low metallicity, low dust, and neutral gas content, and a vast production of (high-energy) ionizing photons. In Fig. 13 one can visualize a boost in the ionizing production efficiency ($\xi_{\text{ion}}$) for high-z LAEs' samples. At the same time, a rise in the ionizing photon flux may reduce the covering fraction or column density of neutral hydrogen, which may ease the escape of LyC and Ly $\alpha$ photons (see Erb et al. 2014). Recently, Flury et al. (2022b) and Schaerer et al. (2022a) statistically demonstrated that the detection rate of LCEs is enhanced among the LAEs and C IV emitter galaxy-population (see also Saxena et al. 2022a; Mascia et al. 2023b), reiterating our results.

Moving forward with the discussion, Fig. 14 shows $\xi_{\text{ion}}$ as function of the UV-continuum slope for the VANDELS sample when using either the R16 (left-hand panel) or the SMC (right-hand panel) attenuation law. Irrespective of the dust law, the ionizing efficiency rapidly decreases with the slope of the continuum, and the trend is more steep for the SMC law. Focusing on the stacks measurements, the canonical limit in $\xi_{\text{ion}}$ (log $\xi_{\text{ion}} \ge 25.2$ Hz erg$^{-1}$) is only achieved for galaxies whose slopes are bluer than $\beta_{1500} \le -1.5$. The $\xi_{\text{ion}} - \beta_{1500}$ relation of a dustless, single-burst stellar population at increasing age is not able to fully reproduce the whole range of $\beta_{1500}$ values, so that it would require 'some' dust reddening to account for redder slopes (see horizontal arrow in Fig. 14, indicating the change in UV slope
<figure><img src="image_16.png" /><figcaption>**Figure 14.** log $\xi_{ion}$ as function of the UV-continuum slope at 1500 Å ($\beta_{1500}$). The red (blue) dots show the results of our $3 \le z \le 5$ VANDELS individual spectra assuming R16 (SMC) dust-attenuation law. Thick symbols represent the results from our composites in bins of $\beta_{1500}$. Orange squares show the MOSDEF results at $z = 2-3$ by Shivaei et al. (2018), while blue diamonds and grey crosses correspond to the HAEs measurements in Bouwens et al. (2016) at $z = 4-5$ Matthee et al. (2017a) at $z = 2-3$, respectively. Open symbols were taken from Lam et al. (2019), and the teal dashed line corresponds to the results by Prieto-Lyon et al. (2023) using HST+JWST data at $z = 3-7$. The grey-shaded area marks the canonical log $\xi_{ion} = 25.2-25.3$ value given by a simple stellar population at constant SFR over 100 Myr (Robertson et al. 2013). log $\xi_{ion}$ decreases with the UV slope almost unanimously for all samples.</figcaption></figure>

after adding 1 mag of reddening). The nebular continuum would also contribute to increase the spectral slope of the youngest bursts of star formation.

Compared to the already reported relations by Bouwens et al. (2016) at $z = 4-5$ and by Matthee et al. (2017a) at $z = 2-3$, our results provide a similar slope but systematically shifted to redder $\beta_{spec}^{1500}$. Several effects can contribute to this offset, but we mostly attribute it to the use of broad-band photometry to compute the UV slopes instead of direct spectral measurements, and the different wavelength range probed by multiband photometry, usually broader than with spectroscopy. In contrast, the flattening of the $\xi_{ion} - \beta_{1500}$ relation proposed by Shivaei et al. (2018) at $\beta_{1500} > -1.6$ is not consistent with our findings. Izotov et al. (2017) reported a similar general trend of decreasing $\xi_{ion}$ with increasing UV slope for local, compact SFGs, but shifted to lower $\xi_{ion}$ than ours. The new results by Prieto-Lyon et al. (2023, teal-dashed line in the plot) using combined HST+JWST photometry are broadly consistent with our estimations.

In summary, our results provide observational evidence for moderately UV-faint, low-mass, and less dusty (un-obscured, blue UV colours) high-redshift galaxies most likely being able to drive the Cosmic Reionization at $6 \le z \le 9$ (Finkelstein et al. 2019; Chisholm et al. 2022; Trebitsch et al. 2022; Lin et al. 2023). According to recent JWST observations (e.g. see Endsley et al. 2022; Cullen et al. 2023; Mascia et al. 2023a), these properties seem to be common among the galaxy population at the reionization epoch. These UV-faint galaxies would produce higher amounts of ionizing photons (through young bursts of star formation, possibly at low metallicities) and allow these photons to escape more easily due to their favoured ISM physical conditions (dustless, gas-less holes in the ISM, see Gazagnes et al. 2020, SL22). Due to their high number density at earlier epochs ($M_{UV}^{*} \approx -21$ at $z = 6$, see Bouwens et al. 2021), they might become the major contributors to the ionizing budget at the EoR (see equation 1).

## 6 THE IONIZING PROPERTIES OF GALAXIES: A DETAILED COMPARISON BETWEEN LOW- AND HIGH-REDSHIFT STUDIES

We now compare the ionizing properties of SFGs at $3 \le z \le 5$ with the estimated at other other redshifts and examine if and how they could evolve into the EoR. To do so, we put our $f_{esc}^{abs}$ and $\xi_{ion}$ trends along with observed properties of galaxies, together with other spectroscopic surveys of LCEs in the literature. On the one hand, we use the LzLCS (Flury et al. 2022a), probing emission line galaxies at $z \simeq 0.3$, as our ‘EoR-analogue sample’ of galaxies. The more extreme nature of the LzLCS galaxies, characterized by high ionization parameters, high SFR surface-densities, and strong nebular emission lines, makes it the perfect sample to compare with, since EoR galaxies are expected (Schaerer et al. 2016; Boyett et al. 2022) and seem to hold these properties (e.g. Endsley et al. 2022; Schaerer et al. 2022b; Lin et al. 2023; Mascia et al. 2023a). At $z \simeq 3$, the KLCS (Steidel et al. 2018; Pahl et al. 2021, 2023) is used for comparison, because it targets LBGs at similar UV magnitudes (and other properties) to the VANDELS galaxies. For both the LzLCS<fnref n="9" /> and the KLCS direct LyC measurements are available. Besides drawing conclusions on the possible redshift evolution of the ionizing properties of galaxies, we will discuss the main systematic uncertainties and caveats of our methodology.

In Fig. 15, the relations between $f_{esc}^{abs}$, $\xi_{ion}$, and the $f_{esc}^{abs} \times \xi_{ion}$ product with the UV magnitude, UV-continuum slope, and Ly$\alpha$ strength are compared for the three samples: VANDELS (at $z = 4$, SMC), KLCS (at $z = 3$), and LzLCS (at $z = 0.3$). Particularly

<footnote n="9">⁹We note that the $f_{esc}^{abs}$ and $\xi_{ion}$ for LzLCS galaxies were obtained by running our FICUS code in Flury et al. (2022a) and SL22. The $f_{esc}^{abs}$ was calculated as the ratio between the intrinsic ionizing flux predicted by FICUS and the observed flux at LyC wavelengths.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** log $\xi_{ion}$ as function of the UV-continuum slope at 1500 Å ($\beta_{1500}$). The red (blue) dots show the results of our $3 \le z \le 5$ VANDELS individual spectra assuming R16 (SMC) dust-attenuation law. Thick symbols represent the results from our composites in bins of $\beta_{1500}$. Orange squares show the MOSDEF results at $z = 2-3$ by Shivaei et al. (2018), while blue diamonds and grey crosses correspond to the HAEs measurements in Bouwens et al. (2016) at $z = 4-5$ Matthee et al. (2017a) at $z = 2-3$, respectively. Open symbols were taken from Lam et al. (2019), and the teal dashed line corresponds to the results by Prieto-Lyon et al. (2023) using HST+JWST data at $z = 3-7$. The grey-shaded area marks the canonical log $\xi_{ion} = 25.2-25.3$ value given by a simple stellar population at constant SFR over 100 Myr (Robertson et al. 2013). log $\xi_{ion}$ decreases with the UV slope almost unanimously for all samples.</figcaption></figure>

after adding 1 mag of reddening). The nebular continuum would also contribute to increase the spectral slope of the youngest bursts of star formation.

Compared to the already reported relations by Bouwens et al. (2016) at $z = 4-5$ and by Matthee et al. (2017a) at $z = 2-3$, our results provide a similar slope but systematically shifted to redder $\beta_{spec}^{1500}$. Several effects can contribute to this offset, but we mostly attribute it to the use of broad-band photometry to compute the UV slopes instead of direct spectral measurements, and the different wavelength range probed by multiband photometry, usually broader than with spectroscopy. In contrast, the flattening of the $\xi_{ion} - \beta_{1500}$ relation proposed by Shivaei et al. (2018) at $\beta_{1500} > -1.6$ is not consistent with our findings. Izotov et al. (2017) reported a similar general trend of decreasing $\xi_{ion}$ with increasing UV slope for local, compact SFGs, but shifted to lower $\xi_{ion}$ than ours. The new results by Prieto-Lyon et al. (2023, teal-dashed line in the plot) using combined HST+JWST photometry are broadly consistent with our estimations.

In summary, our results provide observational evidence for moderately UV-faint, low-mass, and less dusty (un-obscured, blue UV colours) high-redshift galaxies most likely being able to drive the Cosmic Reionization at $6 \le z \le 9$ (Finkelstein et al. 2019; Chisholm et al. 2022; Trebitsch et al. 2022; Lin et al. 2023). According to recent JWST observations (e.g. see Endsley et al. 2022; Cullen et al. 2023; Mascia et al. 2023a), these properties seem to be common among the galaxy population at the reionization epoch. These UV-faint galaxies would produce higher amounts of ionizing photons (through young bursts of star formation, possibly at low metallicities) and allow these photons to escape more easily due to their favoured ISM physical conditions (dustless, gas-less holes in the ISM, see Gazagnes et al. 2020, SL22). Due to their high number density at earlier epochs ($M_{UV}^{*} \approx -21$ at $z = 6$, see Bouwens et al. 2021), they might become the major contributors to the ionizing budget at the EoR (see equation 1).

## 6 THE IONIZING PROPERTIES OF GALAXIES: A DETAILED COMPARISON BETWEEN LOW- AND HIGH-REDSHIFT STUDIES

We now compare the ionizing properties of SFGs at $3 \le z \le 5$ with the estimated at other other redshifts and examine if and how they could evolve into the EoR. To do so, we put our $f_{esc}^{abs}$ and $\xi_{ion}$ trends along with observed properties of galaxies, together with other spectroscopic surveys of LCEs in the literature. On the one hand, we use the LzLCS (Flury et al. 2022a), probing emission line galaxies at $z \simeq 0.3$, as our ‘EoR-analogue sample’ of galaxies. The more extreme nature of the LzLCS galaxies, characterized by high ionization parameters, high SFR surface-densities, and strong nebular emission lines, makes it the perfect sample to compare with, since EoR galaxies are expected (Schaerer et al. 2016; Boyett et al. 2022) and seem to hold these properties (e.g. Endsley et al. 2022; Schaerer et al. 2022b; Lin et al. 2023; Mascia et al. 2023a). At $z \simeq 3$, the KLCS (Steidel et al. 2018; Pahl et al. 2021, 2023) is used for comparison, because it targets LBGs at similar UV magnitudes (and other properties) to the VANDELS galaxies. For both the LzLCS<fnref n="9" /> and the KLCS direct LyC measurements are available. Besides drawing conclusions on the possible redshift evolution of the ionizing properties of galaxies, we will discuss the main systematic uncertainties and caveats of our methodology.

In Fig. 15, the relations between $f_{esc}^{abs}$, $\xi_{ion}$, and the $f_{esc}^{abs} \times \xi_{ion}$ product with the UV magnitude, UV-continuum slope, and Ly$\alpha$ strength are compared for the three samples: VANDELS (at $z = 4$, SMC), KLCS (at $z = 3$), and LzLCS (at $z = 0.3$). Particularly

<footnote n="9">⁹We note that the $f_{esc}^{abs}$ and $\xi_{ion}$ for LzLCS galaxies were obtained by running our FICUS code in Flury et al. (2022a) and SL22. The $f_{esc}^{abs}$ was calculated as the ratio between the intrinsic ionizing flux predicted by FICUS and the observed flux at LyC wavelengths.</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 15.** The absolute photon escape fraction ($f_{\text{esc}}^{\text{abs}}$), the ionizing production efficiency (log $\xi_{\text{ion}}$) and the product of the two (log $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$, see equation 1) for different low-$z$ and high-$z$ samples. Our predictions for the VANDELS galaxies at $3 \le z \le 5$ are shown through dark-blue points in the background (SMC law), and the running locus of the stacking results is indicated via the thick blue line. For comparison, the red open symbols represent the latest KLCS (Steidel et al. 2018) composite results at $z \sim 3$ by Pahl et al. (2021), and the yellow circles show individual measurements at $z \sim 0.3$ from LzLCS (Flury et al. 2022a, triangles are $1\sigma$ upper limits). The dashed grey lines draw the fit relations (median $\pm 1\sigma$ error) to the LzLCS points (Chisholm et al. 2022). Grey-shaded regions mark the 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>

important is how the $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ product compares between the different surveys (third column of the figure).

Interestingly, the properties of the LzLCS and the VANDELS samples show overall fairly consistent trends, despite the redshift range difference and different methods (LzLCS represents direct $f_{\text{esc}}^{\text{abs}}$ measurements, and VANDELS predictions are based on absorption lines). The samples are also quite complementary. Extending over a wide range of UV magnitudes ($-22 \le M_{\odot} \le -18$) and stellar masses ($10^8 M_{\odot} \le M_{\star} \le 10^{10} M_{\odot}$), the escape fraction of LzLCS galaxies tentatively decreases with the UV brightness and stellar mass. At variance with this, VANDELS galaxies span a much narrower range in terms of UV magnitudes, and it is unable to reveal

a clear $f_{\text{esc}}^{\text{abs}}$ correlation with $M_{1500}$. Again, the ionizing efficiencies ($\xi_{\text{ion}}$) of both low- and high-$z$ samples are quite comparable, except for the lack of some low $\xi_{\text{ion}}$ objects and the presence of stronger Ly $\alpha$ emitters in LzLCS.

LzLCS galaxies show systematically bluer UV slopes ($-2.5 \le \beta_{1500} \le -1$) than VANDELS galaxies, which translates into much higher $f_{\text{esc}}^{\text{abs}}$ values for the bluest galaxies in the LzLCS (up to $f_{\text{esc}}^{\text{abs}} \approx 80$ per cent) than any of the galaxies within the VANDELS survey. However, both samples show a significant decrease in their ionizing properties ($f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$) with the UV-continuum slope. Strikingly, our VANDELS results follow the same slope provided by the fit to LzLCS data (grey lines in Fig. 15, see Chisholm et al. 2022), but
<figure><img src="image_17.png" /><figcaption>**Figure 15.** The absolute photon escape fraction ($f_{\text{esc}}^{\text{abs}}$), the ionizing production efficiency (log $\xi_{\text{ion}}$) and the product of the two (log $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$, see equation 1) for different low-$z$ and high-$z$ samples. Our predictions for the VANDELS galaxies at $3 \le z \le 5$ are shown through dark-blue points in the background (SMC law), and the running locus of the stacking results is indicated via the thick blue line. For comparison, the red open symbols represent the latest KLCS (Steidel et al. 2018) composite results at $z \sim 3$ by Pahl et al. (2021), and the yellow circles show individual measurements at $z \sim 0.3$ from LzLCS (Flury et al. 2022a, triangles are $1\sigma$ upper limits). The dashed grey lines draw the fit relations (median $\pm 1\sigma$ error) to the LzLCS points (Chisholm et al. 2022). Grey-shaded regions mark the 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>

important is how the $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ product compares between the different surveys (third column of the figure).

Interestingly, the properties of the LzLCS and the VANDELS samples show overall fairly consistent trends, despite the redshift range difference and different methods (LzLCS represents direct $f_{\text{esc}}^{\text{abs}}$ measurements, and VANDELS predictions are based on absorption lines). The samples are also quite complementary. Extending over a wide range of UV magnitudes ($-22 \le M_{\odot} \le -18$) and stellar masses ($10^8 M_{\odot} \le M_{\star} \le 10^{10} M_{\odot}$), the escape fraction of LzLCS galaxies tentatively decreases with the UV brightness and stellar mass. At variance with this, VANDELS galaxies span a much narrower range in terms of UV magnitudes, and it is unable to reveal

a clear $f_{\text{esc}}^{\text{abs}}$ correlation with $M_{1500}$. Again, the ionizing efficiencies ($\xi_{\text{ion}}$) of both low- and high-$z$ samples are quite comparable, except for the lack of some low $\xi_{\text{ion}}$ objects and the presence of stronger Ly $\alpha$ emitters in LzLCS.

LzLCS galaxies show systematically bluer UV slopes ($-2.5 \le \beta_{1500} \le -1$) than VANDELS galaxies, which translates into much higher $f_{\text{esc}}^{\text{abs}}$ values for the bluest galaxies in the LzLCS (up to $f_{\text{esc}}^{\text{abs}} \approx 80$ per cent) than any of the galaxies within the VANDELS survey. However, both samples show a significant decrease in their ionizing properties ($f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$) with the UV-continuum slope. Strikingly, our VANDELS results follow the same slope provided by the fit to LzLCS data (grey lines in Fig. 15, see Chisholm et al. 2022), but
extrapolated to (≥0.5) redder UV slopes and lower values of $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$. This is interesting as it reassures the applicability of our method for inferring the escape fraction of galaxies (see Chisholm et al. 2018, SL22), since the $f_{\text{esc}}^{\text{abs}}$ quoted for the LzLCS is actually the observed escape fraction from measurements of the LyC flux. However, it is worth mentioning the shift in the intercept between the three surveys, with VANDELS and KLCS giving effectively higher $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ values at a given $\beta_{1500}$ than the extrapolation by Chisholm et al. (2022). This is interesting by itself since (1) it might hint a possible redshift evolution of the $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}} - \beta_{1500}$ relation between $z \simeq 0.3$ (LzLCS) and $z \simeq 3$ (KLCS, VANDELS) or (2), if the high-$z$ surveys probe intrinsically lower UV slopes, they may get a better handle on the actual $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}} - \beta_{1500}$ trend in redder galaxies. In any case, the VANDELS $f_{\text{esc}}^{\text{abs}}$ and $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ relations with $\beta_{1500}$ are compatible with Chisholm et al. (2022) within $1\sigma$, so none of these hypothesis can be confirmed nor ruled out.

Finally, due to the more extreme nature of the LzLCS objects, their $f_{\text{esc}}^{\text{abs}} - W_{\text{Ly}\alpha}$ and $\xi_{\text{ion}} - W_{\text{Ly}\alpha}$ relationships are shifted to higher values of the Ly $\alpha$ equivalent width respect to the VANDELS spectra. The intrinsically higher photon production efficiencies ($\xi_{\text{ion}}$) in LzLCS galaxies boost the strength of the Ly $\alpha$ emission. Regarding the $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ correlation with $W_{\text{Ly}\alpha}$ for the LzLCS, only LCEs (yellow circles, while downward triangles represent non-LCEs) provide $\log f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}(\text{Hz erg}^{-1}) \ge 24.2$ (see Robertson et al. 2013).

Even though our VANDELS survey populate a similar parameters space in terms of $M_{1500}$, $\beta_{1500}$, $W_{\text{Ly}\alpha}$ than the KLCS, the VANDELS $f_{\text{esc}}^{\text{abs}}$ values fall systematically below the KLCS points, the latter acting like an upper envelope in all the correlations. In order to understand these differences, we now carefully list the main limitations and caveats of our ‘picket-fence’ modelling. There are two main assumptions on the applicability of equation (6): (1) the dust-attenuation law and (2) the conversion between the metals and H I covering fractions (see the full dedicated section in SL22).

First, the preference for a shallower (R16) or a steeper (SMC) attenuation curve can be tested by comparing our SED-derived attenuation versus independent measurements of the UV attenuation of galaxies. For example, one can translate IRX-excess measurements of high-$z$ galaxies into UV attenuation at 1600 Å ($A_{1600}$, in mag), by using simple energy balance arguments (see Schaerer, de Barros & Sklias 2013). In Fig. 16 (see Cullen et al. 2018, for details), our VANDELS $A_{1600}$ versus stellar mass results – using both R16 (in red) and SMC (in blue) laws – are compared against the estimations by McLure et al. (2018a) at $z = 2 - 3$ and Fudamoto et al. (2020) at $z = 3 - 4$ from IRX-excess measurements (solid and dashed green lines, respectively). The VANDELS stacks measurements are represented through thick red and thick blue circles. If the more recent analysis of Fudamoto et al. (2020) is adopted, the SMC law, and hence higher $f_{\text{esc}}^{\text{abs}}$, are preferred. However, the McLure et al. (2018a) results at similar redshift would indicate an intermediate attenuation between R16 and SMC. The origin of these differences is difficult to understand and we currently have no way to distinguish which attenuation law is more appropriate. As discussed in Section 4, assuming a most favourable attenuation law like SMC, the $f_{\text{esc}}^{\text{abs}}$ would be $\times 1.5$ higher compared to the ones derived using R16.

Even if the attenuation law was known, $f_{\text{esc}}^{\text{abs}}$ differences may reach depending on the way the UV dust-attenuation ($E_{B\cdot V}$) is measured. The dust-attenuation parameter at UV wavelengths is mainly constraint by the slope of the UV continuum (e.g. $\beta_{1500}$, see Chisholm et al. 2022), by comparing the observed $\beta_{1500}$ against predictions of the intrinsic $\beta_{1500}^{\text{int}}$ from SED fitting, either from the spectra or from multiband photometry. In our case, FICUS is able

<figure><img src="image_18.png" /><figcaption>**Figure 16.** The UV attenuation at 1600 Å ($A_{1600}$, in mag) as function of the stellar mass ($M_*$, in $M_{\odot}$) for individual VANDELS galaxies (dots in the background) and our stacks measurements (big circles). The solid and dashed green thick-lines represent independent estimates by McLure et al. (2018a) and Fudamoto et al. (2020) from IRX-excess measurements at $z = 2 - 3$ and $z = 3 - 4$, respectively. Red symbols assume a R16 dust-attenuation law while the blue symbols assume an SMC attenuation curve.</figcaption></figure>

to estimate $E_{B\cdot V}$ based on the spectral shape of the UV continuum. However, the slope of the current spectral observations might be subject to uncertainty due to different factors (see Garilli et al. 2021): instrumental calibrations, reduction pipeline issues, flux losses because of the atmospheric dispersion and sky subtraction residuals. An alternative solution would be to use an independent determination of the UV slope. With this purpose, we use the $\beta_{1500}$ measurements from Calabrò et al. (2021), taking all the photometric bands whose bandwidths lie entirely inside the 1230–2750 Å rest-frame wavelength range, and the relation given in Chisholm et al. (2022) (equation 8, assuming $\beta_{1500}^{\text{int}} = -2.5$) to convert this $\beta_{1500}$ to $E_{B\cdot V}$. Plugging these new dust-attenuation values into equation (6) produces higher estimations of $f_{\text{esc}}^{\text{abs}}$. The systematic deviation between the spectroscopic (FICUS) and photometric (Calabrò et al. 2021) UV slopes by $\simeq 0.5$ on average, corresponds to $\simeq 1$ mag difference in UV attenuation at LyC wavelengths ($A_{912}$), which leads to a $\times 2.5$ higher photon escape fraction. As sketched in Fig. 17, the new $f_{\text{esc}}^{\text{abs}}$ values (pink pentagons) are in closer agreement with the measurements by Pahl et al. (2021) and Begley et al. (2022), and fully consistent with the physical scenario in which the escape fraction of galaxies decreases towards brighter systems.

Secondly, in Fig. 18 we study the influence that the adopted conversion between the metals and H I covering fraction has on the predicted escape fractions. For comparison, we compute the predicted $f_{\text{esc}}^{\text{abs}}$ values for the KLCS sample adopting the ‘Screen Model’ (see table 9 in Steidel et al. 2018), which is identical to our default ‘picket-fence’ geometry assuming a uniform slab of dust. We are not able to reconcile our escape fractions with the KLCS points, with our predictions lying a factor of $\times 2$ below them. However, if one assumes a 1:1 metals-to-H I covering fraction correspondence so that $C_f(\text{LIS}) = C_f(\text{H I})$, instead of the usual linear regression by which $C_f(\text{H I}) > C_f(\text{LIS})$) (see Gazagnes et al. 2018 and SL22), the resulting $f_{\text{esc}}^{\text{abs}}$ would agree with the estimate of the average $\langle f_{\text{esc}}^{\text{abs}} \rangle$
extrapolated to (≥0.5) redder UV slopes and lower values of $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$. This is interesting as it reassures the applicability of our method for inferring the escape fraction of galaxies (see Chisholm et al. 2018, SL22), since the $f_{\text{esc}}^{\text{abs}}$ quoted for the LzLCS is actually the observed escape fraction from measurements of the LyC flux. However, it is worth mentioning the shift in the intercept between the three surveys, with VANDELS and KLCS giving effectively higher $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ values at a given $\beta_{1500}$ than the extrapolation by Chisholm et al. (2022). This is interesting by itself since (1) it might hint a possible redshift evolution of the $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}} - \beta_{1500}$ relation between $z \simeq 0.3$ (LzLCS) and $z \simeq 3$ (KLCS, VANDELS) or (2), if the high-$z$ surveys probe intrinsically lower UV slopes, they may get a better handle on the actual $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}} - \beta_{1500}$ trend in redder galaxies. In any case, the VANDELS $f_{\text{esc}}^{\text{abs}}$ and $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ relations with $\beta_{1500}$ are compatible with Chisholm et al. (2022) within $1\sigma$, so none of these hypothesis can be confirmed nor ruled out.

Finally, due to the more extreme nature of the LzLCS objects, their $f_{\text{esc}}^{\text{abs}} - W_{\text{Ly}\alpha}$ and $\xi_{\text{ion}} - W_{\text{Ly}\alpha}$ relationships are shifted to higher values of the Ly $\alpha$ equivalent width respect to the VANDELS spectra. The intrinsically higher photon production efficiencies ($\xi_{\text{ion}}$) in LzLCS galaxies boost the strength of the Ly $\alpha$ emission. Regarding the $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ correlation with $W_{\text{Ly}\alpha}$ for the LzLCS, only LCEs (yellow circles, while downward triangles represent non-LCEs) provide $\log f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}(\text{Hz erg}^{-1}) \ge 24.2$ (see Robertson et al. 2013).

Even though our VANDELS survey populate a similar parameters space in terms of $M_{1500}$, $\beta_{1500}$, $W_{\text{Ly}\alpha}$ than the KLCS, the VANDELS $f_{\text{esc}}^{\text{abs}}$ values fall systematically below the KLCS points, the latter acting like an upper envelope in all the correlations. In order to understand these differences, we now carefully list the main limitations and caveats of our ‘picket-fence’ modelling. There are two main assumptions on the applicability of equation (6): (1) the dust-attenuation law and (2) the conversion between the metals and H I covering fractions (see the full dedicated section in SL22).

First, the preference for a shallower (R16) or a steeper (SMC) attenuation curve can be tested by comparing our SED-derived attenuation versus independent measurements of the UV attenuation of galaxies. For example, one can translate IRX-excess measurements of high-$z$ galaxies into UV attenuation at 1600 Å ($A_{1600}$, in mag), by using simple energy balance arguments (see Schaerer, de Barros & Sklias 2013). In Fig. 16 (see Cullen et al. 2018, for details), our VANDELS $A_{1600}$ versus stellar mass results – using both R16 (in red) and SMC (in blue) laws – are compared against the estimations by McLure et al. (2018a) at $z = 2 - 3$ and Fudamoto et al. (2020) at $z = 3 - 4$ from IRX-excess measurements (solid and dashed green lines, respectively). The VANDELS stacks measurements are represented through thick red and thick blue circles. If the more recent analysis of Fudamoto et al. (2020) is adopted, the SMC law, and hence higher $f_{\text{esc}}^{\text{abs}}$, are preferred. However, the McLure et al. (2018a) results at similar redshift would indicate an intermediate attenuation between R16 and SMC. The origin of these differences is difficult to understand and we currently have no way to distinguish which attenuation law is more appropriate. As discussed in Section 4, assuming a most favourable attenuation law like SMC, the $f_{\text{esc}}^{\text{abs}}$ would be $\times 1.5$ higher compared to the ones derived using R16.

Even if the attenuation law was known, $f_{\text{esc}}^{\text{abs}}$ differences may reach depending on the way the UV dust-attenuation ($E_{B\cdot V}$) is measured. The dust-attenuation parameter at UV wavelengths is mainly constraint by the slope of the UV continuum (e.g. $\beta_{1500}$, see Chisholm et al. 2022), by comparing the observed $\beta_{1500}$ against predictions of the intrinsic $\beta_{1500}^{\text{int}}$ from SED fitting, either from the spectra or from multiband photometry. In our case, FICUS is able

<figure><img src="image_18.png" /><figcaption>**Figure 16.** The UV attenuation at 1600 Å ($A_{1600}$, in mag) as function of the stellar mass ($M_*$, in $M_{\odot}$) for individual VANDELS galaxies (dots in the background) and our stacks measurements (big circles). The solid and dashed green thick-lines represent independent estimates by McLure et al. (2018a) and Fudamoto et al. (2020) from IRX-excess measurements at $z = 2 - 3$ and $z = 3 - 4$, respectively. Red symbols assume a R16 dust-attenuation law while the blue symbols assume an SMC attenuation curve.</figcaption></figure>

to estimate $E_{B\cdot V}$ based on the spectral shape of the UV continuum. However, the slope of the current spectral observations might be subject to uncertainty due to different factors (see Garilli et al. 2021): instrumental calibrations, reduction pipeline issues, flux losses because of the atmospheric dispersion and sky subtraction residuals. An alternative solution would be to use an independent determination of the UV slope. With this purpose, we use the $\beta_{1500}$ measurements from Calabrò et al. (2021), taking all the photometric bands whose bandwidths lie entirely inside the 1230–2750 Å rest-frame wavelength range, and the relation given in Chisholm et al. (2022) (equation 8, assuming $\beta_{1500}^{\text{int}} = -2.5$) to convert this $\beta_{1500}$ to $E_{B\cdot V}$. Plugging these new dust-attenuation values into equation (6) produces higher estimations of $f_{\text{esc}}^{\text{abs}}$. The systematic deviation between the spectroscopic (FICUS) and photometric (Calabrò et al. 2021) UV slopes by $\simeq 0.5$ on average, corresponds to $\simeq 1$ mag difference in UV attenuation at LyC wavelengths ($A_{912}$), which leads to a $\times 2.5$ higher photon escape fraction. As sketched in Fig. 17, the new $f_{\text{esc}}^{\text{abs}}$ values (pink pentagons) are in closer agreement with the measurements by Pahl et al. (2021) and Begley et al. (2022), and fully consistent with the physical scenario in which the escape fraction of galaxies decreases towards brighter systems.

Secondly, in Fig. 18 we study the influence that the adopted conversion between the metals and H I covering fraction has on the predicted escape fractions. For comparison, we compute the predicted $f_{\text{esc}}^{\text{abs}}$ values for the KLCS sample adopting the ‘Screen Model’ (see table 9 in Steidel et al. 2018), which is identical to our default ‘picket-fence’ geometry assuming a uniform slab of dust. We are not able to reconcile our escape fractions with the KLCS points, with our predictions lying a factor of $\times 2$ below them. However, if one assumes a 1:1 metals-to-H I covering fraction correspondence so that $C_f(\text{LIS}) = C_f(\text{H I})$, instead of the usual linear regression by which $C_f(\text{H I}) > C_f(\text{LIS})$) (see Gazagnes et al. 2018 and SL22), the resulting $f_{\text{esc}}^{\text{abs}}$ would agree with the estimate of the average $\langle f_{\text{esc}}^{\text{abs}} \rangle$
<figure><img src="image_19.png" /><figcaption>**Figure 17.** A representation of the potential systematic bias in $f_{\text{esc}}^{\text{abs}}$ due to the different values in the UV dust-attenuation that can be derived from dissimilar measurements of the β-slope. The blue diamonds show our default escape fractions from FICUS, i.e. $E_{\text{B-V}}$ is derived from spectroscopic measurements of the UV continuum slope ($\beta_{\text{spec}}$). Conversely, the pink pentagons display the $f_{\text{esc}}^{\text{abs}}$ values when $E_{\text{B-V}}$ is inferred from the β-slope measurements by Calabrò et al. (2021), using photometry only ($\beta_{\text{phot}}$). The average $\langle f_{\text{esc}}^{\text{abs}} \rangle = 0.07 \pm 0.02$ reported by Begley et al. (2022) at $z \simeq 3.5$ is indicated through the grey shaded-area, and the results by Pahl et al. (2021) are plotted through black squares. The effect of spectral resolution due to changes in the depths of the absorption lines is represented with the dashed line.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** A representation of the potential systematic bias in $f_{\text{esc}}^{\text{abs}}$ (blue diamonds) due to the conversion between the metals and neutral gas covering fractions (blue dotted line). For comparison we also show the predicted $f_{\text{esc}}^{\text{abs}}$ values for the KLCS stacks (red symbols), by applying our methodology to the published covering fractions and UV attenuations from Steidel et al. (2018). The layout is the same as in Fig. 17.</figcaption></figure>

by Begley et al. (2022, grey shaded area), which is now much more compatible with KLCS at the same time.

Finally, the resolution correction function applied to the residual flux of the LIS lines (Appendix A) might also affect our escape fraction estimates. In Figs 17 and 18, we explicitly show the effect of our custom correction in the inferred $f_{\text{esc}}^{\text{abs}}$ (dashed blue line). In the more conservative case in which the spectral resolution is totally dismissed and non-considered in equation (6), $f_{\text{esc}}^{\text{abs}}$ would increase by a factor of $\times 1 - 1.5$ at most: this is by far the effect with the least influence on the average $f_{\text{esc}}^{\text{abs}}$ among the ones considered in this paper.

In conclusion, different effects can potentially augment our average $f_{\text{esc}}^{\text{abs}}$. Among them, the determination of the UV dust-attenuation has the strongest impact in the escape fraction estimates. In Appendix C (Fig. C1), we plot the $f_{\text{esc}}^{\text{abs}}$ values resulting from plugging into equation (6) the dust-attenuation values derived from the $\beta_{1500}$ measurements from Calabrò et al. (2021). As we can see, all VANDELS, KLCS, and LzLCS $f_{\text{esc}}^{\text{abs}}$ correlations with physical quantities now converge to the same global trends, and they share a common behaviour by which $f_{\text{esc}}^{\text{abs}}$ increases towards fainter UV magnitudes, bluer UV slopes, and stronger Ly$\alpha$ emission. In view of these results, the use of $\beta_{1500}$ from photometry as a proxy of the UV dust attenuation in combination with the depth of the absorption lines as a tracer of H I-empty channels in the ISM of galaxies (equation 6), seem to be a robust way of predicting the escape fraction of galaxies at any redshift. Independently, there is also the chance of (1) adopting a very steep dust-attenuation law, whose $k_\lambda$ extrapolation towards LyC wavelengths would be even higher than when using an SMC-like law and (2) using a favourable metals-to-H I covering fraction's conversion, in which the metals would spatially trace the neutral gas more closely than with the current assumed relation.

The remaining question is whether a universal relation between the ionizing properties of galaxies ($f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$) and their physical properties ($M_{1500}$, $\beta_{1500}$, $W_{\text{Ly}\alpha}$ etc.) exists, in which the different surveys would populate different regions of the hypothetical universal trend or, contrarily, the observed differences between surveys (VANDELS, KLCS, LzLCS) indicates an evolution of these properties with redshift. Still, as far as our sample is concerned, the $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ dependence along with galaxy properties is mainly due to escape fraction variations, so that the redshift evolution of this product will ultimately follow the changes in the ISM geometry, gas, and dust composition across cosmic time.

In any case, and regardless the caveats that we have just described, the LzLCS, KLCS, and VANDELS results suggest a scenario in which moderately UV-faint, low-mass, and dustless galaxies most likely dominate the SFG ionizing photon budget at any epoch. The stellar populations of EoR-galaxies would be characterized by strong radiation fields, with high ionizing production efficiencies (boosting the nebular emission), and whose ISM conditions would be such that they facilitate the escape of a copious amount of LyC photons to the IGM. This can be speculatively attributed to the stellar mass build-up and changes in the gas and dust properties of galaxies between the Cosmic Noon ($z \simeq 3$) and the EoR ($6 \le z \le 9$).

# 7 SUMMARY AND CONCLUSIONS

In this work, we fully exploit and highlight the ability of absorption line and rest-UV spectroscopy alone to decipher the ionizing properties of high-z SFGs. Our novel approach make use of deep, rest-frame UV spectra from the VANDELS survey at $3 \le z \le 5$, to compute the ionizing absolute photon escape fraction ($f_{\text{esc}}^{\text{abs}}$) and ionizing photon production efficiency ($\xi_{\text{ion}}$) of galaxies. $f_{\text{esc}}^{\text{abs}}$ has been derived by
<figure><img src="image_19.png" /><figcaption>**Figure 17.** A representation of the potential systematic bias in $f_{\text{esc}}^{\text{abs}}$ due to the different values in the UV dust-attenuation that can be derived from dissimilar measurements of the β-slope. The blue diamonds show our default escape fractions from FICUS, i.e. $E_{\text{B-V}}$ is derived from spectroscopic measurements of the UV continuum slope ($\beta_{\text{spec}}$). Conversely, the pink pentagons display the $f_{\text{esc}}^{\text{abs}}$ values when $E_{\text{B-V}}$ is inferred from the β-slope measurements by Calabrò et al. (2021), using photometry only ($\beta_{\text{phot}}$). The average $\langle f_{\text{esc}}^{\text{abs}} \rangle = 0.07 \pm 0.02$ reported by Begley et al. (2022) at $z \simeq 3.5$ is indicated through the grey shaded-area, and the results by Pahl et al. (2021) are plotted through black squares. The effect of spectral resolution due to changes in the depths of the absorption lines is represented with the dashed line.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** A representation of the potential systematic bias in $f_{\text{esc}}^{\text{abs}}$ (blue diamonds) due to the conversion between the metals and neutral gas covering fractions (blue dotted line). For comparison we also show the predicted $f_{\text{esc}}^{\text{abs}}$ values for the KLCS stacks (red symbols), by applying our methodology to the published covering fractions and UV attenuations from Steidel et al. (2018). The layout is the same as in Fig. 17.</figcaption></figure>

by Begley et al. (2022, grey shaded area), which is now much more compatible with KLCS at the same time.

Finally, the resolution correction function applied to the residual flux of the LIS lines (Appendix A) might also affect our escape fraction estimates. In Figs 17 and 18, we explicitly show the effect of our custom correction in the inferred $f_{\text{esc}}^{\text{abs}}$ (dashed blue line). In the more conservative case in which the spectral resolution is totally dismissed and non-considered in equation (6), $f_{\text{esc}}^{\text{abs}}$ would increase by a factor of $\times 1 - 1.5$ at most: this is by far the effect with the least influence on the average $f_{\text{esc}}^{\text{abs}}$ among the ones considered in this paper.

In conclusion, different effects can potentially augment our average $f_{\text{esc}}^{\text{abs}}$. Among them, the determination of the UV dust-attenuation has the strongest impact in the escape fraction estimates. In Appendix C (Fig. C1), we plot the $f_{\text{esc}}^{\text{abs}}$ values resulting from plugging into equation (6) the dust-attenuation values derived from the $\beta_{1500}$ measurements from Calabrò et al. (2021). As we can see, all VANDELS, KLCS, and LzLCS $f_{\text{esc}}^{\text{abs}}$ correlations with physical quantities now converge to the same global trends, and they share a common behaviour by which $f_{\text{esc}}^{\text{abs}}$ increases towards fainter UV magnitudes, bluer UV slopes, and stronger Ly$\alpha$ emission. In view of these results, the use of $\beta_{1500}$ from photometry as a proxy of the UV dust attenuation in combination with the depth of the absorption lines as a tracer of H I-empty channels in the ISM of galaxies (equation 6), seem to be a robust way of predicting the escape fraction of galaxies at any redshift. Independently, there is also the chance of (1) adopting a very steep dust-attenuation law, whose $k_\lambda$ extrapolation towards LyC wavelengths would be even higher than when using an SMC-like law and (2) using a favourable metals-to-H I covering fraction's conversion, in which the metals would spatially trace the neutral gas more closely than with the current assumed relation.

The remaining question is whether a universal relation between the ionizing properties of galaxies ($f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$) and their physical properties ($M_{1500}$, $\beta_{1500}$, $W_{\text{Ly}\alpha}$ etc.) exists, in which the different surveys would populate different regions of the hypothetical universal trend or, contrarily, the observed differences between surveys (VANDELS, KLCS, LzLCS) indicates an evolution of these properties with redshift. Still, as far as our sample is concerned, the $f_{\text{esc}}^{\text{abs}} \times \xi_{\text{ion}}$ dependence along with galaxy properties is mainly due to escape fraction variations, so that the redshift evolution of this product will ultimately follow the changes in the ISM geometry, gas, and dust composition across cosmic time.

In any case, and regardless the caveats that we have just described, the LzLCS, KLCS, and VANDELS results suggest a scenario in which moderately UV-faint, low-mass, and dustless galaxies most likely dominate the SFG ionizing photon budget at any epoch. The stellar populations of EoR-galaxies would be characterized by strong radiation fields, with high ionizing production efficiencies (boosting the nebular emission), and whose ISM conditions would be such that they facilitate the escape of a copious amount of LyC photons to the IGM. This can be speculatively attributed to the stellar mass build-up and changes in the gas and dust properties of galaxies between the Cosmic Noon ($z \simeq 3$) and the EoR ($6 \le z \le 9$).

# 7 SUMMARY AND CONCLUSIONS

In this work, we fully exploit and highlight the ability of absorption line and rest-UV spectroscopy alone to decipher the ionizing properties of high-z SFGs. Our novel approach make use of deep, rest-frame UV spectra from the VANDELS survey at $3 \le z \le 5$, to compute the ionizing absolute photon escape fraction ($f_{\text{esc}}^{\text{abs}}$) and ionizing photon production efficiency ($\xi_{\text{ion}}$) of galaxies. $f_{\text{esc}}^{\text{abs}}$ has been derived by
combining absorption line measurements with estimates of the UV attenuation (see Saldana-Lopez et al. 2022), while the $\xi_{\text{ion}}$ parameter was computed by fitting the FUV stellar continuum of the VANDELS galaxies (following Chisholm et al. 2019).

In particular, we have searched for correlations between $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$ along with different galaxy properties (UV magnitude, UV slope, Ly $\alpha$ strength, etc.), and we thoroughly compared with independent literature estimates. We found that:

(i) The predicted $f_{\text{esc}}^{\text{abs}}$ monotonically decreases with the stellar mass, the UV-continuum slope, and the Ly $\alpha$ equivalent width of the VANDELS galaxies (Figs 5, 6, 7). We find a non-significant correlation between $f_{\text{esc}}^{\text{abs}}$ and the UV magnitude, although the faintest galaxies tentatively have higher escape fractions.

(ii) The estimated $\xi_{\text{ion}}$ statistically increases towards blue UV-continuum slopes and strong Ly $\alpha$ emitting galaxies (Figs 12 and 14), and it smoothly raises beyond the canonical value towards the UV-faintest galaxies in the sample.

(iii) Potential Lyman Continuum Emitters (LCEs) and selected Lyman Alpha Emitters (LAEs) show systematically higher $\xi_{\text{ion}}$ (log $\xi_{\text{ion}}$(Hz erg$^{-1}$) $\approx$ 25.38, 25.41, respectively) than non-LCEs and non-LAEs at similar UV magnitudes (log $\xi_{\text{ion}}$(Hz erg$^{-1}$) $\approx$ 25.18, 25.14), and our $\xi_{\text{ion}}$ values are in agreement with other LAEs surveys in the literature (Fig. 13).

Additionally, we constructed the average composite FUV spectrum of LCEs at $3 \le z \le 5$ (Figs 8, 10, and 11), by stacking potential, individual emitters in the VANDELS survey (selected based on our predicted $f_{\text{esc}}^{\text{abs}} \ge 5$ per cent), and explained their non-ionizing spectral properties in the framework of the ISM conditions which enable ionizing photons to escape. Our results show that the FUV spectra of typical high-$z$ LCEs would be characterized by:

(i) Blue UV slopes ($\beta_{\text{spec}}^{1500} \lesssim -2$) compared to non-LCEs, which directly translates into low UV attenuation ($E_{B-V} \lesssim 0.1$ mag) and therefore low column density of dust along the line-of-sight.

(ii) Enhanced Ly $\alpha$ emission ($W_{Ly\alpha} \lesssim -25$ Å) and strong UV nebular lines in contrast to the non-LCEs population, particularly high C III I 1908/C IV 1550 ratios ($\gtrsim 0.75$). Together with the intrinsically higher $\xi_{\text{ion}}$ parameter for LCEs, this indicates very young underlying stellar populations ($\approx 10$ Myr) at relatively low metallicities ($\approx 0.2$ Z$_{\odot}$).

(iii) Weak ($\lesssim 1$ Å) ISM LIS absorption lines (e.g. Si II 1260, C II 1334), while the HIS absorption lines (e.g. Si IV 1400) are of similar strength as the bulk of non-LCEs. This, together with the low UV attenuation for LCEs, suggest the presence of dustless, low gas column-density channels in the ISM which favour the escape of ionizing photons.

Finally, we have compared our findings with other LyC surveys in the literature (Fig. 15), concretely the *Keck Lyman Continuum Survey* (or KLCS; Steidel et al. 2018; Pahl et al. 2021, 2023), targeting $z \simeq 3$ LBGs at similar UV magnitudes to VANDELS, and the *Low-Redshift Lyman Continuum Survey* (or LzLCS; Flury et al. 2022a, b), which targetted emission line galaxies at $z \simeq 0.3$. VANDELS and LzLCS show overall fairly consistent trends, with LzLCS shifted to fainter UV magnitudes, bluer UV slopes and stronger Ly $\alpha$ emission, and therefore their $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$ are enhanced with respect to VANDELS galaxies. The escape fraction of KLCS galaxies fall above our estimates at all UV magnitudes. We mainly ascribed this discrepancy to the way the amount of UV attenuation is measured, and we propose to use $\beta_{1500}$ measurements from photometry as an independent proxy of the UV dust attenuation of galaxies (Fig. 17). The dust-attenuation law (Fig. 16) and the metals-to-neutral-gas

covering fraction conversion (Fig. 18) constitutes additional sources of uncertainty to the escape fraction.

Our joint analysis of the VANDELS, LzLCS, and KLCS results shed light onto the fact that UV-faint, low-mass, and dustless galaxies likely dominated the ionizing budget during the EoR. Their stellar populations would most likely be characterized by strong radiation fields, with high ionizing production efficiencies, and whose ISM conditions are such that they favour the escape of LyC photons. Recent JWST observations (e.g. Endsley et al. 2022; Schaerer et al. 2022b; Cullen et al. 2023; Lin et al. 2023; Mascia et al. 2023a), have revealed that blue UV slopes and strong emission lines also characterized the less massive and moderately faint galaxies at the EoR, supporting our results.

An increasing number of high-quality, FUV spectra at $z > 6$ will be available in the mid-future thanks to upcoming ground-based (ELT, GMT, TMT, ...) and currently ongoing space-based facilities (JWST). This will give us the first insights into the properties of galaxies during the EoR. Nevertheless, according to our work, additional efforts are still needed in order to correctly decipher the physics underlying such vast data sets. The two main questions that we think are urgent and clearly needed for future studies are: (1) if unique, what is the dust-attenuation law governing low-mass, low-metallicity galaxies? And (2) do metals trace the same spatial location as the neutral gas within the ISM of these galaxies?

## ACKNOWLEDGEMENTS

The authors thank the anonymous referee for providing useful comments, which have certainly improved the quality of this paper. The authors also thank Hakim Atek, Jorytt Matthee, and Irene Shivaei for kindly providing tabulated data from Atek et al. (2022), Matthee et al. (2017a), and Shivaei et al. (2018). ASL and DS acknowledge support from Swiss National Science Foundation. The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant No. 140. RA acknowledges support from ANID FONDECYT Regular 1202007. JPUF is supported by the DFF grant 1026-00066.

## DATA AVAILABILITY

The VANDELS Data Release 4 (DR4) is now publicly available and can be accessed using the VANDELS data base at http://vandels.in af.it/dr4.html, or through the ESO archives. Any secondary product and/or data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Alavi A. et al., 2020, ApJ, 904, 59

* Arellano-Córdova K. Z. et al., 2022, ApJ, 940, L23

* Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

* Atek H., Kunth D., Schaerer D., Mas-Hesse J. M., Hayes M., Östlin G., Kneib J.-P., 2014, A&A, 561, A89

* Atek H. et al., 2015, ApJ, 800, 18

* Atek H., Furtak L. J., Oesch P., van Dokkum P., Reddy N., Contini T., Illingworth G., Wilkins S., 2022, MNRAS, 511, 4464

* Atek H. et al., 2023, MNRAS, 519, 1201

* Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023

* Becker G. D., D'Aloisio A., Christenson H. M., Zhu Y., Worseck G., Bolton J. S., 2021, MNRAS, 508, 1853

* Begley R. et al., 2022, MNRAS, 513, 3510

* Bezanson R. et al., 2022, preprint (arXiv:2212.04026)

* Bian F., Fan X., 2020, MNRAS, 493, L65
combining absorption line measurements with estimates of the UV attenuation (see Saldana-Lopez et al. 2022), while the $\xi_{\text{ion}}$ parameter was computed by fitting the FUV stellar continuum of the VANDELS galaxies (following Chisholm et al. 2019).

In particular, we have searched for correlations between $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$ along with different galaxy properties (UV magnitude, UV slope, Ly $\alpha$ strength, etc.), and we thoroughly compared with independent literature estimates. We found that:

(i) The predicted $f_{\text{esc}}^{\text{abs}}$ monotonically decreases with the stellar mass, the UV-continuum slope, and the Ly $\alpha$ equivalent width of the VANDELS galaxies (Figs 5, 6, 7). We find a non-significant correlation between $f_{\text{esc}}^{\text{abs}}$ and the UV magnitude, although the faintest galaxies tentatively have higher escape fractions.

(ii) The estimated $\xi_{\text{ion}}$ statistically increases towards blue UV-continuum slopes and strong Ly $\alpha$ emitting galaxies (Figs 12 and 14), and it smoothly raises beyond the canonical value towards the UV-faintest galaxies in the sample.

(iii) Potential Lyman Continuum Emitters (LCEs) and selected Lyman Alpha Emitters (LAEs) show systematically higher $\xi_{\text{ion}}$ (log $\xi_{\text{ion}}$(Hz erg$^{-1}$) $\approx$ 25.38, 25.41, respectively) than non-LCEs and non-LAEs at similar UV magnitudes (log $\xi_{\text{ion}}$(Hz erg$^{-1}$) $\approx$ 25.18, 25.14), and our $\xi_{\text{ion}}$ values are in agreement with other LAEs surveys in the literature (Fig. 13).

Additionally, we constructed the average composite FUV spectrum of LCEs at $3 \le z \le 5$ (Figs 8, 10, and 11), by stacking potential, individual emitters in the VANDELS survey (selected based on our predicted $f_{\text{esc}}^{\text{abs}} \ge 5$ per cent), and explained their non-ionizing spectral properties in the framework of the ISM conditions which enable ionizing photons to escape. Our results show that the FUV spectra of typical high-$z$ LCEs would be characterized by:

(i) Blue UV slopes ($\beta_{\text{spec}}^{1500} \lesssim -2$) compared to non-LCEs, which directly translates into low UV attenuation ($E_{B-V} \lesssim 0.1$ mag) and therefore low column density of dust along the line-of-sight.

(ii) Enhanced Ly $\alpha$ emission ($W_{Ly\alpha} \lesssim -25$ Å) and strong UV nebular lines in contrast to the non-LCEs population, particularly high C III I 1908/C IV 1550 ratios ($\gtrsim 0.75$). Together with the intrinsically higher $\xi_{\text{ion}}$ parameter for LCEs, this indicates very young underlying stellar populations ($\approx 10$ Myr) at relatively low metallicities ($\approx 0.2$ Z$_{\odot}$).

(iii) Weak ($\lesssim 1$ Å) ISM LIS absorption lines (e.g. Si II 1260, C II 1334), while the HIS absorption lines (e.g. Si IV 1400) are of similar strength as the bulk of non-LCEs. This, together with the low UV attenuation for LCEs, suggest the presence of dustless, low gas column-density channels in the ISM which favour the escape of ionizing photons.

Finally, we have compared our findings with other LyC surveys in the literature (Fig. 15), concretely the *Keck Lyman Continuum Survey* (or KLCS; Steidel et al. 2018; Pahl et al. 2021, 2023), targeting $z \simeq 3$ LBGs at similar UV magnitudes to VANDELS, and the *Low-Redshift Lyman Continuum Survey* (or LzLCS; Flury et al. 2022a, b), which targetted emission line galaxies at $z \simeq 0.3$. VANDELS and LzLCS show overall fairly consistent trends, with LzLCS shifted to fainter UV magnitudes, bluer UV slopes and stronger Ly $\alpha$ emission, and therefore their $f_{\text{esc}}^{\text{abs}}$ and $\xi_{\text{ion}}$ are enhanced with respect to VANDELS galaxies. The escape fraction of KLCS galaxies fall above our estimates at all UV magnitudes. We mainly ascribed this discrepancy to the way the amount of UV attenuation is measured, and we propose to use $\beta_{1500}$ measurements from photometry as an independent proxy of the UV dust attenuation of galaxies (Fig. 17). The dust-attenuation law (Fig. 16) and the metals-to-neutral-gas

covering fraction conversion (Fig. 18) constitutes additional sources of uncertainty to the escape fraction.

Our joint analysis of the VANDELS, LzLCS, and KLCS results shed light onto the fact that UV-faint, low-mass, and dustless galaxies likely dominated the ionizing budget during the EoR. Their stellar populations would most likely be characterized by strong radiation fields, with high ionizing production efficiencies, and whose ISM conditions are such that they favour the escape of LyC photons. Recent JWST observations (e.g. Endsley et al. 2022; Schaerer et al. 2022b; Cullen et al. 2023; Lin et al. 2023; Mascia et al. 2023a), have revealed that blue UV slopes and strong emission lines also characterized the less massive and moderately faint galaxies at the EoR, supporting our results.

An increasing number of high-quality, FUV spectra at $z > 6$ will be available in the mid-future thanks to upcoming ground-based (ELT, GMT, TMT, ...) and currently ongoing space-based facilities (JWST). This will give us the first insights into the properties of galaxies during the EoR. Nevertheless, according to our work, additional efforts are still needed in order to correctly decipher the physics underlying such vast data sets. The two main questions that we think are urgent and clearly needed for future studies are: (1) if unique, what is the dust-attenuation law governing low-mass, low-metallicity galaxies? And (2) do metals trace the same spatial location as the neutral gas within the ISM of these galaxies?

## ACKNOWLEDGEMENTS

The authors thank the anonymous referee for providing useful comments, which have certainly improved the quality of this paper. The authors also thank Hakim Atek, Jorytt Matthee, and Irene Shivaei for kindly providing tabulated data from Atek et al. (2022), Matthee et al. (2017a), and Shivaei et al. (2018). ASL and DS acknowledge support from Swiss National Science Foundation. The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant No.140. RA acknowledges support from ANID FONDECYT Regular 1202007. JPUF is supported by the DFF grant 1026-00066.

## DATA AVAILABILITY

The VANDELS Data Release 4 (DR4) is now publicly available and can be accessed using the VANDELS data base at http://vandels.in af.it/dr4.html, or through the ESO archives. Any secondary product and/or data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Alavi A. et al., 2020, ApJ, 904, 59

* Arellano-Córdova K. Z. et al., 2022, ApJ, 940, L23

* Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

* Atek H., Kunth D., Schaerer D., Mas-Hesse J. M., Hayes M., Östlin G., Kneib J.-P., 2014, A&A, 561, A89

* Atek H. et al., 2015, ApJ, 800, 18

* Atek H., Furtak L. J., Oesch P., van Dokkum P., Reddy N., Contini T., Illingworth G., Wilkins S., 2022, MNRAS, 511, 4464

* Atek H. et al., 2023, MNRAS, 519, 1201

* Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023

* Becker G. D., D'Aloisio A., Christenson H. M., Zhu Y., Worseck G., Bolton J. S., 2021, MNRAS, 508, 1853

* Begley R. et al., 2022, MNRAS, 513, 3510

* Bezanson R. et al., 2022, preprint (arXiv:2212.04026)

* Bian F., Fan X., 2020, MNRAS, 493, L65
Bouchet P., Lequeux J., Maurice E., Prevot L., Prevot-Burnichon M. L., 1985, A&A, 149, 330

Bouwens R. J. et al., 2012, ApJ, 754, 83

Bouwens R. J. et al., 2014, ApJ, 793, 115

Bouwens R. J. et al., 2015, ApJ, 803, 34

Bouwens R. J., Smit R., Labbé I., Franx M., Caruana J., Oesch P., Stefanon M., Rasappu N., 2016, ApJ, 831, 176

Bouwens R. J. et al., 2021, AJ, 162, 47

Boyett K. N. K., Stark D. P., Bunker A. J., Tang M., Maseda M. V., 2022, MNRAS, 513, 4451

Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127

Bridge C. R. et al., 2010, ApJ, 720, 465

Brinchmann J., 2022, preprint (arXiv:2208.07467)

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Cain C., D'Aloisio A., Gangolli N., Becker G. D., 2021, ApJ, 917, L37

Calabrò A. et al., 2021, A&A, 646, A39

Calabrò A. et al., 2022, A&A, 667, A117

Calzetti D., Kinney A. L., Storchi-Bergmann T., 1994, ApJ, 429, 582

Cameron A. J. et al., 2023, preprint (arXiv:2302.04298)

Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379

Carnall A. C. et al., 2023, MNRAS, 518, L45

Castellano M. et al., 2022a, preprint (arXiv:2212.06666)

Castellano M. et al., 2022b, ApJ, 938, L15

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

Chisholm J., Rigby J. R., Bayliss M., Berg D. A., Dahle H., Gladders M., Sharon K., 2019, ApJ, 882, 182 (C19)

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

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

Cristiani S., Serrano L. M., Fontanot F., Vanzella E., Monaco P., 2016, MNRAS, 462, 2478

Cullen F. et al., 2018, MNRAS, 476, 3218

Cullen F. et al., 2019, MNRAS, 487, 2038

Cullen F. et al., 2020, MNRAS, 495, 1501

Cullen F. et al., 2021, MNRAS, 505, 903

Cullen F. et al., 2023, MNRAS, 520, 14

Curti M. et al., 2023, MNRAS, 518, 425

Curtis-Lake E. et al., 2023, Nature Astronomy

Davidzon I. et al., 2017, A&A, 605, A70

Dayal P., Ferrara A., 2018, Phys. Rep., 780, 1

Dayal P. et al., 2020, MNRAS, 495, 3065

de Barros S. et al., 2016, A&A, 585, A51

Donnan C. T. et al., 2023, MNRAS, 518, 6011

Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium, Princeton Series in Astrophysics. Princeton Univ. Press, Princeton

Du X. et al., 2018, ApJ, 860, 75

Emami N., Siana B., Alavi A., Gburek T., Freeman W. R., Richard J., Weisz D. R., Stark D. P., 2020, ApJ, 895, 116

Endsley R., Stark D. P., Whitler L., Topping M. W., Chen Z., Plat A., Chisholm J., Charlot S., 2022, preprint (arXiv:2208.14999)

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

Faisst A. L., 2016, ApJ, 829, 99

Falcón-Barroso J., Sánchez-Blázquez P., Vazdekis A., Ricciardelli E., Cardiel N., Cenarro A. J., Gorgas J., Peletier R. F., 2011, A&A, 532, A95

Ferland G. J. et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385

Finkelstein S. L. et al., 2012, ApJ, 756, 164

Finkelstein S. L. et al., 2019, ApJ, 879, 36

Finkelstein S. L. et al., 2023, ApJ, 946, L13

Fletcher T. J., Tang M., Robertson B. E., Nakajima K., Ellis R. S., Stark D. P., Inoue A., 2019, ApJ, 878, 87

Flury S. R. et al., 2022a, ApJS, 260, 1

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

Fontanot F., Cristiani S., Vanzella E., 2012, MNRAS, 425, 1413

Fudamoto Y. et al., 2020, MNRAS, 491, 4724

Fujimoto S. et al., 2023, preprint (arXiv:2301.09482)

Furtak L. J., Shuntov M., Atek H., Zitrin A., Richard J., Lehnert M. D., Chevallard J., 2023, MNRAS, 519, 3064

Galametz A. et al., 2013, ApJS, 206, 10

Garilli B. et al., 2021, A&A, 647, A150

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

Gazagnes S., Chisholm J., Schaerer D., Verhamme A., Izotov Y., 2020, A&A, 639, A85

Gordon K. D., Clayton G. C., Misselt K. A., Landolt A. U., Wolff M. J., 2003, ApJ, 594, 279 (SMC)

Goto H. et al., 2021, ApJ, 923, 229

Grazian A. et al., 2016, A&A, 585, A48

Grazian A. et al., 2017, A&A, 602, A18

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

Guaita L. et al., 2016, A&A, 587, A133

Guo Y. et al., 2013, ApJS, 207, 24

Harikane Y. et al., 2018, ApJ, 859, 84

Harikane Y. et al., 2022, ApJ, 929, 1

Harikane Y. et al., 2023, ApJS, 265, 5

Hassan S., Davé R., Mitra S., Finlator K., Ciardi B., Santos M. G., 2018, MNRAS, 473, 227

Hathi N. P. et al., 2016, A&A, 588, A26

Heckman T. M., Sembach K. R., Meurer G. R., Leitherer C., Calzetti D., Martin C. L., 2001, ApJ, 558, 56

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

Henry A., Scarlata C., Martin C. L., Erb D., 2015, ApJ, 809, 19

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

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

Isobe Y., Ouchi M., Nakajima K., Harikane Y., Ono Y., Xu Y., Zhang Y., Umeda H., 2023, preprint (arXiv:2301.06811)

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

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

Izotov Y. I., Guseva N. G., Fricke K. J., Henkel C., Schaerer D., 2017, MNRAS, 467, 4118

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

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

Izotov Y. I., Schaerer D., Worseck G., Verhamme A., Guseva N. G., Thuan T. X., Orlitová I., Fricke K. J., 2020, MNRAS, 491, 468

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

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

Japelj J. et al., 2017, MNRAS, 468, 389

Ji Z. et al., 2020, ApJ, 888, 109

Jones T., Stark D. P., Ellis R. S., 2012, ApJ, 751, 51

Jones T. A., Ellis R. S., Schenker M. A., Stark D. P., 2013, ApJ, 779, 52

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

Katz H. et al., 2022, MNRAS, 515, 4265

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

Kornei K. A., Shapley A. E., Erb D. K., Steidel C. C., Reddy N. A., Pettini M., Bogosavljević M., 2010, ApJ, 711, 693

Kroupa P., 2001, MNRAS, 322, 231

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

Kusakabe H. et al., 2020, A&A, 638, A12

Labbe I. et al., 2023, Nature, 616, 266

Lam D. et al., 2019, A&A, 627, A164

Le Fèvre O. et al., 2003, in Iye M., Moorwood A. F. M., eds, Proc. SPIE Conf. Ser. Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes. SPIE, Bellingham, p. 1670

Le Fèvre O. et al., 2015, A&A, 576, A79

Leitherer C., Ortiz Otálvaro P. A., Bresolin F., Kudritzki R.-P., Lo Faro B., Pauldrach A. W. A., Pettini M., Rix S. A., 2010, ApJS, 189, 309

Leitherer C. et al., 2011, Starburst99: Synthesis Models for Galaxies with Active Star Formation, record ascl:1104.003
Bouchet P., Lequeux J., Maurice E., Prevot L., Prevot-Burnichon M. L., 1985, A&A, 149, 330

Bouwens R. J. et al., 2012, ApJ, 754, 83

Bouwens R. J. et al., 2014, ApJ, 793, 115

Bouwens R. J. et al., 2015, ApJ, 803, 34

Bouwens R. J., Smit R., Labbé I., Franx M., Caruana J., Oesch P., Stefanon M., Rasappu N., 2016, ApJ, 831, 176

Bouwens R. J. et al., 2021, AJ, 162, 47

Boyett K. N. K., Stark D. P., Bunker A. J., Tang M., Maseda M. V., 2022, MNRAS, 513, 4451

Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127

Bridge C. R. et al., 2010, ApJ, 720, 465

Brinchmann J., 2022, preprint (arXiv:2208.07467)

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Cain C., D'Aloisio A., Gangolli N., Becker G. D., 2021, ApJ, 917, L37

Calabrò A. et al., 2021, A&A, 646, A39

Calabrò A. et al., 2022, A&A, 667, A117

Calzetti D., Kinney A. L., Storchi-Bergmann T., 1994, ApJ, 429, 582

Cameron A. J. et al., 2023, preprint (arXiv:2302.04298)

Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379

Carnall A. C. et al., 2023, MNRAS, 518, L45

Castellano M. et al., 2022a, preprint (arXiv:2212.06666)

Castellano M. et al., 2022b, ApJ, 938, L15

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

Chisholm J., Rigby J. R., Bayliss M., Berg D. A., Dahle H., Gladders M., Sharon K., 2019, ApJ, 882, 182 (C19)

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

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

Cristiani S., Serrano L. M., Fontanot F., Vanzella E., Monaco P., 2016, MNRAS, 462, 2478

Cullen F. et al., 2018, MNRAS, 476, 3218

Cullen F. et al., 2019, MNRAS, 487, 2038

Cullen F. et al., 2020, MNRAS, 495, 1501

Cullen F. et al., 2021, MNRAS, 505, 903

Cullen F. et al., 2023, MNRAS, 520, 14

Curti M. et al., 2023, MNRAS, 518, 425

Curtis-Lake E. et al., 2023, Nature Astronomy

Davidzon I. et al., 2017, A&A, 605, A70

Dayal P., Ferrara A., 2018, Phys. Rep., 780, 1

Dayal P. et al., 2020, MNRAS, 495, 3065

de Barros S. et al., 2016, A&A, 585, A51

Donnan C. T. et al., 2023, MNRAS, 518, 6011

Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium, Princeton Series in Astrophysics. Princeton Univ. Press, Princeton

Du X. et al., 2018, ApJ, 860, 75

Emami N., Siana B., Alavi A., Gburek T., Freeman W. R., Richard J., Weisz D. R., Stark D. P., 2020, ApJ, 895, 116

Endsley R., Stark D. P., Whitler L., Topping M. W., Chen Z., Plat A., Chisholm J., Charlot S., 2022, preprint (arXiv:2208.14999)

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

Faisst A. L., 2016, ApJ, 829, 99

Falcón-Barroso J., Sánchez-Blázquez P., Vazdekis A., Ricciardelli E., Cardiel N., Cenarro A. J., Gorgas J., Peletier R. F., 2011, A&A, 532, A95

Ferland G. J. et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385

Finkelstein S. L. et al., 2012, ApJ, 756, 164

Finkelstein S. L. et al., 2019, ApJ, 879, 36

Finkelstein S. L. et al., 2023, ApJ, 946, L13

Fletcher T. J., Tang M., Robertson B. E., Nakajima K., Ellis R. S., Stark D. P., Inoue A., 2019, ApJ, 878, 87

Flury S. R. et al., 2022a, ApJS, 260, 1

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

Fontanot F., Cristiani S., Vanzella E., 2012, MNRAS, 425, 1413

Fudamoto Y. et al., 2020, MNRAS, 491, 4724

Fujimoto S. et al., 2023, preprint (arXiv:2301.09482)

Furtak L. J., Shuntov M., Atek H., Zitrin A., Richard J., Lehnert M. D., Chevallard J., 2023, MNRAS, 519, 3064

Galametz A. et al., 2013, ApJS, 206, 10

Garilli B. et al., 2021, A&A, 647, A150

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

Gazagnes S., Chisholm J., Schaerer D., Verhamme A., Izotov Y., 2020, A&A, 639, A85

Gordon K. D., Clayton G. C., Misselt K. A., Landolt A. U., Wolff M. J., 2003, ApJ, 594, 279 (SMC)

Goto H. et al., 2021, ApJ, 923, 229

Grazian A. et al., 2016, A&A, 585, A48

Grazian A. et al., 2017, A&A, 602, A18

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

Guaita L. et al., 2016, A&A, 587, A133

Guo Y. et al., 2013, ApJS, 207, 24

Harikane Y. et al., 2018, ApJ, 859, 84

Harikane Y. et al., 2022, ApJ, 929, 1

Harikane Y. et al., 2023, ApJS, 265, 5

Hassan S., Davé R., Mitra S., Finlator K., Ciardi B., Santos M. G., 2018, MNRAS, 473, 227

Hathi N. P. et al., 2016, A&A, 588, A26

Heckman T. M., Sembach K. R., Meurer G. R., Leitherer C., Calzetti D., Martin C. L., 2001, ApJ, 558, 56

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

Henry A., Scarlata C., Martin C. L., Erb D., 2015, ApJ, 809, 19

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

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

Isobe Y., Ouchi M., Nakajima K., Harikane Y., Ono Y., Xu Y., Zhang Y., Umeda H., 2023, preprint (arXiv:2301.06811)

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

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

Izotov Y. I., Guseva N. G., Fricke K. J., Henkel C., Schaerer D., 2017, MNRAS, 467, 4118

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

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

Izotov Y. I., Schaerer D., Worseck G., Verhamme A., Guseva N. G., Thuan T. X., Orlitová I., Fricke K. J., 2020, MNRAS, 491, 468

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

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

Japelj J. et al., 2017, MNRAS, 468, 389

Ji Z. et al., 2020, ApJ, 888, 109

Jones T., Stark D. P., Ellis R. S., 2012, ApJ, 751, 51

Jones T. A., Ellis R. S., Schenker M. A., Stark D. P., 2013, ApJ, 779, 52

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

Katz H. et al., 2022, MNRAS, 515, 4265

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

Kornei K. A., Shapley A. E., Erb D. K., Steidel C. C., Reddy N. A., Pettini M., Bogosavljević M., 2010, ApJ, 711, 693

Kroupa P., 2001, MNRAS, 322, 231

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

Kusakabe H. et al., 2020, A&A, 638, A12

Labbe I. et al., 2023, Nature, 616, 266

Lam D. et al., 2019, A&A, 627, A164

Le Fèvre O. et al., 2003, in Iye M., Moorwood A. F. M., eds, Proc. SPIE Conf. Ser. Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes. SPIE, Bellingham, p. 1670

Le Fèvre O. et al., 2015, A&A, 576, A79

Leitherer C., Ortiz Otálvaro P. A., Bresolin F., Kudritzki R.-P., Lo Faro B., Pauldrach A. W. A., Pettini M., Rix S. A., 2010, ApJS, 189, 309

Leitherer C. et al., 2011, Starburst99: Synthesis Models for Galaxies with Active Star Formation, record ascl:1104.003
Leitherer C., Ekström S., Meynet G., Schaerer D., Agienko K. B., Levesque E. M., 2014, ApJS, 212, 14

Lin Y.-H. et al., 2023, preprint (arXiv:2303.04572)

Livermore R. C., Finkelstein S. L., Lotz J. M., 2017, ApJ, 835, 113

Llerena M. et al., 2022, A&A, 659, A16

Ma X., Quataert E., Wetzel A., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2020, MNRAS, 498, 2001

Madau P., Haardt F., 2015, ApJ, 813, L8

Madau P., Haardt F., Rees M. J., 1999, ApJ, 514, 648

Maiolino R., Mannucci F., 2019, A&AR, 27, 3

Maji M. et al., 2022, A&A, 663, A66

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

Marigo P., Bressan A., Nanni A., Girardi L., Pumo M. L., 2013, MNRAS, 434, 488

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

Marques-Chaves R. et al., 2022a, MNRAS, 517, 2972

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

Mascia S. et al., 2023a, A&A, 672, A155

Mascia S. et al., 2023b, preprint (arXiv:2301.09328)

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

Mason C. A., Naidu R. P., Tacchella S., Leja J., 2019, MNRAS, 489, 2669

Matsuoka Y. et al., 2018, ApJ, 869, 150

Matthee J., Sobral D., Best P., Khostovan A. A., Oteo I., Bouwens R., Röttgering H., 2017a, MNRAS, 465, 3637

Matthee J., Sobral D., Darvish B., Santos S., Mobasher B., Paulino-Afonso A., Röttgering H., Alegre L., 2017b, MNRAS, 472, 772

Matthee J. et al., 2022, MNRAS, 512, 5960

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

McLeod D. J., McLure R. J., Dunlop J. S., 2016, MNRAS, 459, 3812

McLure R. J. et al., 2018a, MNRAS, 476, 3991

McLure R. J. et al., 2018b, MNRAS, 479, 25

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

Meurer G. R., Heckman T. M., Calzetti D., 1999, ApJ, 521, 64

Meynet G., Maeder A., Schaller G., Schaerer D., Charbonnel C., 1994, A&AS, 103, 97

Micheva G., Iwata I., Inoue A. K., Matsuda Y., Yamada T., Hayashino T., 2017, MNRAS, 465, 316

Mostardi R. E., Shapley A. E., Steidel C. C., Trainor R. F., Reddy N. A., Siana B., 2015, ApJ, 810, 107

Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, MNRAS, 454, 2691

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

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

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

Nakajima K., Ellis R. S., Iwata I., Inoue A. K., Kusakabe H., Ouchi M., Robertson B. E., 2016, ApJ, 831, L9

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

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

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

Nakajima K., Ouchi M., Isobe Y., Harikane Y., Zhang Y., Ono Y., Umeda H., Oguri M., 2023, preprint (arXiv:2301.12825)

Nanayakkara T. et al., 2023, ApJ, 944, L1

Newville M., Stensitzki T., Allen D. B., Rawlik M., Ingargiola A., Nelson A., 2016, Lmfit: Non-Linear Least-Square Minimization and Curve-Fitting for Python, record ascl:1606.014

Ning Y., Cai Z., Jiang L., Lin X., Fu S., Spinoso D., 2023, ApJ, 944, L1

Oesch P. A., Bouwens R. J., Illingworth G. D., Labbé I., Stefanon M., 2018, ApJ, 855, 105

Oke J. B., Gunn J. E., 1983, ApJ, 266, 713

Pahl A. J., Shapley A., Faisst A. L., Capak P. L., Du X., Reddy N. A., Laursen P., Topping M. W., 2020, MNRAS, 493, 3194

Pahl A. J., Shapley A., Steidel C. C., Chen Y., Reddy N. A., 2021, MNRAS, 505, 2447

Pahl A. J., Shapley A., Steidel C. C., Reddy N. A., Chen Y., Rudie G. C., Strom A. L., 2023, MNRAS, 521, 3247

Pauldrach A. W. A., Hoffmann T. L., Lennon M., 2001, A&A, 375, 161

Pentericci L., Grazian A., Fontana A., Castellano M., Giallongo E., Salimbeni S., Santini P., 2009, A&A, 494, 553

Pentericci L. et al., 2018, A&A, 616, A174

Planck Collaboration XIII, 2016, A&A, 596, A108

Prevot M. L., Lequeux J., Maurice E., Prevot L., Rocca-Volmerange B., 1984, A&A, 132, 389

Prichard L. J. et al., 2022, ApJ, 924, 14

Prieto-Lyon G. et al., 2023, A&A, 672, A186

Reddy N. A., Steidel C. C., Pettini M., Bogosavljević M., 2016a, ApJ, 828, 107, (R16)

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

Reddy N. A. et al., 2018, ApJ, 869, 92

Reddy N. A. et al., 2022, ApJ, 926, 31

Rivera-Thorsen T. E., Hayes M., Melinder J., 2022, A&A, 666, A145

Robertson B. E., 2022, ARA&A, 60, 121

Robertson B. E. et al., 2013, ApJ, 768, 71

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

Robertson B. E. et al., 2023, Nature Astronomy

Rojas-Ruiz S., Finkelstein S. L., Bagley M. B., Stevans M., Finkelstein K. D., Larson R., Mechtley M., Diekmann J., 2020, ApJ, 891, 146

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

Rutkowski M. J. et al., 2017, ApJ, 841, L27

Saha K. et al., 2020, Nat. Astron., 4, 1185

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

Salim S., Boquien M., Lee J. C., 2018, ApJ, 859, 11

Santini P. et al., 2023, ApJ, 942, L27

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

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

Saxena A. et al., 2023, preprint (arXiv:2302.12805)

Schaerer D., de Barros S., Sklias P., 2013, A&A, 549, A4

Schaerer D., Izotov Y. I., Verhamme A., Orlitová I., Thuan T. X., Worseck G., Guseva N. G., 2016, A&A, 591, L8

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

Schaerer D., Marques-Chaves R., Barrufet L., Oesch P., Izotov Y. I., Naidu R., Guseva N. G., Brammer G., 2022b, A&A, 665, L4

Seiler J., Hutter A., Sinha M., Croton D., 2018, MNRAS, 480, L33

Seive T., Chisholm J., Leclercq F., Zeimann G., 2022, MNRAS, 515, 5556

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

Shapley A. E., Steidel C. C., Pettini M., Adelberger K. L., Erb D. K., 2006, ApJ, 651, 688

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

Sharma M., Theuns T., Frenk C., Bower R., Crain R., Schaller M., Schaye J., 2016, MNRAS, 458, L94

Shivaei I. et al., 2018, ApJ, 855, 42

Shivaei I. et al., 2020, ApJ, 899, 117

Smith A., Safranek-Shrader C., Bromm V., Milosavljević M., 2015, MNRAS, 449, 4336

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15

Stark D. P., Ellis R. S., Ouchi M., 2011, ApJ, 728, L2

Stark D. P. et al., 2015, MNRAS, 454, 1393

Steidel C. C., Pettini M., Adelberger K. L., 2001, ApJ, 546, 665

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

Tacchella S. et al., 2022, preprint (arXiv:2208.03281)

Tang M., Stark D. P., Chevallard J., Charlot S., 2019, MNRAS, 489, 2572

Tang M. et al., 2023, preprint (arXiv:2301.07072)
Leitherer C., Ekström S., Meynet G., Schaerer D., Agienko K. B., Levesque E. M., 2014, ApJS, 212, 14

Lin Y.-H. et al., 2023, preprint (arXiv:2303.04572)

Livermore R. C., Finkelstein S. L., Lotz J. M., 2017, ApJ, 835, 113

Llerena M. et al., 2022, A&A, 659, A16

Ma X., Quataert E., Wetzel A., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2020, MNRAS, 498, 2001

Madau P., Haardt F., 2015, ApJ, 813, L8

Madau P., Haardt F., Rees M. J., 1999, ApJ, 514, 648

Maiolino R., Mannucci F., 2019, A&AR, 27, 3

Maji M. et al., 2022, A&A, 663, A66

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

Marigo P., Bressan A., Nanni A., Girardi L., Pumo M. L., 2013, MNRAS, 434, 488

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

Marques-Chaves R. et al., 2022a, MNRAS, 517, 2972

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

Mascia S. et al., 2023a, A&A, 672, A155

Mascia S. et al., 2023b, preprint (arXiv:2301.09328)

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

Mason C. A., Naidu R. P., Tacchella S., Leja J., 2019, MNRAS, 489, 2669

Matsuoka Y. et al., 2018, ApJ, 869, 150

Matthee J., Sobral D., Best P., Khostovan A. A., Oteo I., Bouwens R., Röttgering H., 2017a, MNRAS, 465, 3637

Matthee J., Sobral D., Darvish B., Santos S., Mobasher B., Paulino-Afonso A., Röttgering H., Alegre L., 2017b, MNRAS, 472, 772

Matthee J. et al., 2022, MNRAS, 512, 5960

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

McLeod D. J., McLure R. J., Dunlop J. S., 2016, MNRAS, 459, 3812

McLure R. J. et al., 2018a, MNRAS, 476, 3991

McLure R. J. et al., 2018b, MNRAS, 479, 25

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

Meurer G. R., Heckman T. M., Calzetti D., 1999, ApJ, 521, 64

Meynet G., Maeder A., Schaller G., Schaerer D., Charbonnel C., 1994, A&AS, 103, 97

Micheva G., Iwata I., Inoue A. K., Matsuda Y., Yamada T., Hayashino T., 2017, MNRAS, 465, 316

Mostardi R. E., Shapley A. E., Steidel C. C., Trainor R. F., Reddy N. A., Siana B., 2015, ApJ, 810, 107

Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, MNRAS, 454, 2691

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

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

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

Nakajima K., Ellis R. S., Iwata I., Inoue A. K., Kusakabe H., Ouchi M., Robertson B. E., 2016, ApJ, 831, L9

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

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

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

Nakajima K., Ouchi M., Isobe Y., Harikane Y., Zhang Y., Ono Y., Umeda H., Oguri M., 2023, preprint (arXiv:2301.12825)

Nanayakkara T. et al., 2023, ApJ, 944, L1

Newville M., Stensitzki T., Allen D. B., Rawlik M., Ingargiola A., Nelson A., 2016, Lmfit: Non-Linear Least-Square Minimization and Curve-Fitting for Python, record ascl:1606.014

Ning Y., Cai Z., Jiang L., Lin X., Fu S., Spinoso D., 2023, ApJ, 944, L1

Oesch P. A., Bouwens R. J., Illingworth G. D., Labbé I., Stefanon M., 2018, ApJ, 855, 105

Oke J. B., Gunn J. E., 1983, ApJ, 266, 713

Pahl A. J., Shapley A., Faisst A. L., Capak P. L., Du X., Reddy N. A., Laursen P., Topping M. W., 2020, MNRAS, 493, 3194

Pahl A. J., Shapley A., Steidel C. C., Chen Y., Reddy N. A., 2021, MNRAS, 505, 2447

Pahl A. J., Shapley A., Steidel C. C., Reddy N. A., Chen Y., Rudie G. C., Strom A. L., 2023, MNRAS, 521, 3247

Pauldrach A. W. A., Hoffmann T. L., Lennon M., 2001, A&A, 375, 161

Pentericci L., Grazian A., Fontana A., Castellano M., Giallongo E., Salimbeni S., Santini P., 2009, A&A, 494, 553

Pentericci L. et al., 2018, A&A, 616, A174

Planck Collaboration XIII, 2016, A&A, 596, A108

Prevot M. L., Lequeux J., Maurice E., Prevot L., Rocca-Volmerange B., 1984, A&A, 132, 389

Prichard L. J. et al., 2022, ApJ, 924, 14

Prieto-Lyon G. et al., 2023, A&A, 672, A186

Reddy N. A., Steidel C. C., Pettini M., Bogosavljević M., 2016a, ApJ, 828, 107, (R16)

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

Reddy N. A. et al., 2018, ApJ, 869, 92

Reddy N. A. et al., 2022, ApJ, 926, 31

Rivera-Thorsen T. E., Hayes M., Melinder J., 2022, A&A, 666, A145

Robertson B. E., 2022, ARA&A, 60, 121

Robertson B. E. et al., 2013, ApJ, 768, 71

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

Robertson B. E. et al., 2023, Nature Astronomy

Rojas-Ruiz S., Finkelstein S. L., Bagley M. B., Stevans M., Finkelstein K. D., Larson R., Mechtley M., Diekmann J., 2020, ApJ, 891, 146

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

Rutkowski M. J. et al., 2017, ApJ, 841, L27

Saha K. et al., 2020, Nat. Astron., 4, 1185

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

Salim S., Boquien M., Lee J. C., 2018, ApJ, 859, 11

Santini P. et al., 2023, ApJ, 942, L27

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

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

Saxena A. et al., 2023, preprint (arXiv:2302.12805)

Schaerer D., de Barros S., Sklias P., 2013, A&A, 549, A4

Schaerer D., Izotov Y. I., Verhamme A., Orlitová I., Thuan T. X., Worseck G., Guseva N. G., 2016, A&A, 591, L8

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

Schaerer D., Marques-Chaves R., Barrufet L., Oesch P., Izotov Y. I., Naidu R., Guseva N. G., Brammer G., 2022b, A&A, 665, L4

Seiler J., Hutter A., Sinha M., Croton D., 2018, MNRAS, 480, L33

Seive T., Chisholm J., Leclercq F., Zeimann G., 2022, MNRAS, 515, 5556

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

Shapley A. E., Steidel C. C., Pettini M., Adelberger K. L., Erb D. K., 2006, ApJ, 651, 688

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

Sharma M., Theuns T., Frenk C., Bower R., Crain R., Schaller M., Schaye J., 2016, MNRAS, 458, L94

Shivaei I. et al., 2018, ApJ, 855, 42

Shivaei I. et al., 2020, ApJ, 899, 117

Smith A., Safranek-Shrader C., Bromm V., Milosavljević M., 2015, MNRAS, 449, 4336

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15

Stark D. P., Ellis R. S., Ouchi M., 2011, ApJ, 728, L2

Stark D. P. et al., 2015, MNRAS, 454, 1393

Steidel C. C., Pettini M., Adelberger K. L., 2001, ApJ, 546, 665

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

Tacchella S. et al., 2022, preprint (arXiv:2208.03281)

Tang M., Stark D. P., Chevallard J., Charlot S., 2019, MNRAS, 489, 2572

Tang M. et al., 2023, preprint (arXiv:2301.07072)
Tanvir N. R. et al., 2019, MNRAS, 483, 5380

Topping M. W., Stark D. P., Endsley R., Plat A., Whitler L., Chen Z., Charlot S., 2022, ApJ, 941, 153

Trainor R. F., Steidel C. C., Strom A. L., Rudie G. C., 2015, ApJ, 809, 89

Trainor R. F., Strom A. L., Steidel C. C., Rudie G. C., Chen Y., Theios R. L., 2019, ApJ, 887, 85

Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 470, 224

Trebitsch M. et al., 2022, preprint (arXiv:2212.06177)

Trump J. R. et al., 2023, ApJ, 945, 35

Trussler J. A. A. et al., 2022, preprint (arXiv:2207.14265)

Vanzella E. et al., 2010, ApJ, 725, 1011

Vanzella E. et al., 2012, ApJ, 751, 70

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

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

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

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

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

Vielfaure J. B. et al., 2021, A&A, 653, A83

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

Weaver J. R. et al., 2023, preprint (arXiv:2301.02671)

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

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

## APPENDIX A: THE IMPACT OF RESOLUTION AND STACKING ON THE DEPTH OF THE ABSORPTION LINES

Here we provide two calibrations to correct the observed line-depths measurements when (1) inherently degraded by the instrumental resolution and (2) artificially smoothed during usual stacking procedures. As we discussed in the main text, the spectrograph resolution makes the absorption lines wider but less deep, conserving the flux. Together with the noise, sensitivity, and other possible instrumental aberrations, this can lead to a systematically overestimated measurement of the residual flux.

To account for these effects, we simulate SiII 1260 line intensities ($I_\lambda$) assuming a foreground dust-screen geometry of the picket-fence model, describing the line through a single gas component as (Draine 2011):

$$I_\lambda = C_f(\lambda) \times \exp(-\sigma_\lambda N_\lambda) + 1 - C_f(\lambda), \quad (A1)$$

where $C_f(\lambda)$ represents the covering fraction of the line in question and the $\tau_\lambda = \sigma_\lambda N_\lambda$ product is usually known as the optical depth of the line. The line cross-section, $\sigma_\lambda$, which shapes the absorption profile, is given by the Voigt function so that:

$$\sigma_\lambda = \frac{\sqrt{\pi} e^2 f_\lambda}{m_e c} \frac{\lambda}{b} \text{Voigt}(\lambda, A_\lambda, b), \quad (A2)$$

where oscillator strengths $f_\lambda$ and Einstein coefficients $A_\lambda$ are taken from the NIST data base.<fnref n="10" /> For the Voigt function, Voigt($\lambda, A_\lambda, b$), we

use the numerical approximation described in Smith et al. (2015), and assume the following Gaussian distributions ($N(\mu, \sigma^2)$) of the input parameters:

(i) $N \sim N(16, 1) \text{ cm}^{-2}$

(ii) $b \sim N(125, 25) \text{ km s}^{-1}$

(iii) $v \sim N(0, 150) \text{ km s}^{-1}$, $\lambda = \lambda_0 + v/c$

where $\mu$ represents the median and $\sigma^2$ the standard deviation of the normal distribution. The gas column density distribution ($N, \text{cm}^{-2}$) is chosen so that the SiII 1260 equivalent width of the line falls within the optically thick limits of the curve-of-growth or, in other words, the line is always saturated. The line velocity shift ($v, \text{km s}^{-1}$) randomly changes to account for inflows or outflowing gas so that the line centre is blue- or red-shifted accordingly, i.e. $\lambda = \lambda_0 + v/c$, being $\lambda_0 = 1260.42 \AA$ for the SiII 1260 line. Finally, the gas thermal (or Doppler) broadening ($b, \text{km s}^{-1}$) varies around typical values for normal star-forming galaxies (Steidel et al. 2018). The input covering fraction is fixed for every set of simulations to $C_f(\lambda) = 0.3, 0.5, 0.65, 0.75, 0.85$, and 1, respectively.

After performing a set of $\times 100$ ideal Voigt profiles for each covering fraction, we introduce the effect of resolution by convolving the simulated spectra with a Gaussian kernel whose FWHM matches the instrumental resolution (a.k.a., $R$ (VIMOS) $\approx 600$). A constant S/N value ($S/N = 5$, the median of the VANDELS sample) is then added to every single spectra. At fixed $C_f(\lambda)$, the median measured (observed) residual flux of each mock distribution is compared to the theoretical line residual flux of $R_\lambda = 1 - C_f(\lambda)$ according to the picket-fence configuration, producing the desired calibration. Fig. A1 presents the resulting line profiles for a set of $\times 100$ simulations of the SiII 1260 line, all of them with an input covering fraction of $C_f(\lambda) = 0.85$ ($R_\lambda = 0.15$). The median of the distribution of individually measured $R_\lambda$ values is 0.3 (red line), showing that the observed line depth can statistically differ as much as 40 per cent with respect to the real value.

In the last step, we proceed to stack all the spectra contained on each simulation pack at fixed $C_f(\lambda)$, following the same methods described in the text. The co-added spectra for the $C_f(\lambda) = 0.85$ simulation is also shown in Fig. A1 (blue line), with a 65 per cent relative difference in residual flux with respect to the real input value.

To sum up, the resulting calibrations – fitting a linear polynomial to the measured versus theoretical $R_\lambda$ values – are:

$$R_{\text{SIII}} [\text{R} + \text{S/N}] = 0.68 \times (1 - C_f(\text{SiII})) + 0.20, \quad (A3)$$

and

$$R_{\text{SIII}} [\text{stacking}] = 0.73 \times (1 - C_f(\text{SiII})) + 0.29, \quad (A4)$$

when accounting for instrumental resolution plus S/N (equation A3), and spectral stacking (equation A4), respectively. Both relations are plotted in Fig. A2 through the dashed red and dashed blue lines. These calibrations were applied to individual residual flux measurements as well as when measuring the residual flux over composite spectra.

<footnote n="10">¹⁰NIST stands for National Institute of Standards and Technology: https://physics.nist.gov/PhysRefData/ASD/lines_form.html.</footnote>
Tanvir N. R. et al., 2019, MNRAS, 483, 5380

Topping M. W., Stark D. P., Endsley R., Plat A., Whitler L., Chen Z., Charlot S., 2022, ApJ, 941, 153

Trainor R. F., Steidel C. C., Strom A. L., Rudie G. C., 2015, ApJ, 809, 89

Trainor R. F., Strom A. L., Steidel C. C., Rudie G. C., Chen Y., Theios R. L., 2019, ApJ, 887, 85

Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 470, 224

Trebitsch M. et al., 2022, preprint (arXiv:2212.06177)

Trump J. R. et al., 2023, ApJ, 945, 35

Trussler J. A. A. et al., 2022, preprint (arXiv:2207.14265)

Vanzella E. et al., 2010, ApJ, 725, 1011

Vanzella E. et al., 2012, ApJ, 751, 70

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

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

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

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

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

Vielfaure J. B. et al., 2021, A&A, 653, A83

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

Weaver J. R. et al., 2023, preprint (arXiv:2301.02671)

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

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

## APPENDIX A: THE IMPACT OF RESOLUTION AND STACKING ON THE DEPTH OF THE ABSORPTION LINES

Here we provide two calibrations to correct the observed line-depths measurements when (1) inherently degraded by the instrumental resolution and (2) artificially smoothed during usual stacking procedures. As we discussed in the main text, the spectrograph resolution makes the absorption lines wider but less deep, conserving the flux. Together with the noise, sensitivity, and other possible instrumental aberrations, this can lead to a systematically overestimated measurement of the residual flux.

To account for these effects, we simulate SiII 1260 line intensities ($I_\lambda$) assuming a foreground dust-screen geometry of the picket-fence model, describing the line through a single gas component as (Draine 2011):

$$I_\lambda = C_f(\lambda) \times \exp(-\sigma_\lambda N_\lambda) + 1 - C_f(\lambda), \quad (A1)$$

where $C_f(\lambda)$ represents the covering fraction of the line in question and the $\tau_\lambda = \sigma_\lambda N_\lambda$ product is usually known as the optical depth of the line. The line cross-section, $\sigma_\lambda$, which shapes the absorption profile, is given by the Voigt function so that:

$$\sigma_\lambda = \frac{\sqrt{\pi} e^2 f_\lambda}{m_e c} \frac{\lambda}{b} \text{Voigt}(\lambda, A_\lambda, b), \quad (A2)$$

where oscillator strengths $f_\lambda$ and Einstein coefficients $A_\lambda$ are taken from the NIST data base.<fnref n="10" /> For the Voigt function, $\text{Voigt}(\lambda, A_\lambda, b)$, we

use the numerical approximation described in Smith et al. (2015), and assume the following Gaussian distributions ($N(\mu, \sigma^2)$) of the input parameters:

(i) $N \sim N(16, 1) \text{ cm}^{-2}$

(ii) $b \sim N(125, 25) \text{ km s}^{-1}$

(iii) $v \sim N(0, 150) \text{ km s}^{-1}$, $\lambda = \lambda_0 + v/c$

where $\mu$ represents the median and $\sigma^2$ the standard deviation of the normal distribution. The gas column density distribution ($N, \text{cm}^{-2}$) is chosen so that the SiII 1260 equivalent width of the line falls within the optically thick limits of the curve-of-growth or, in other words, the line is always saturated. The line velocity shift ($v, \text{km s}^{-1}$) randomly changes to account for inflows or outflowing gas so that the line centre is blue- or red-shifted accordingly, i.e. $\lambda = \lambda_0 + v/c$, being $\lambda_0 = 1260.42 \AA$ for the SiII 1260 line. Finally, the gas thermal (or Doppler) broadening ($b, \text{km s}^{-1}$) varies around typical values for normal star-forming galaxies (Steidel et al. 2018). The input covering fraction is fixed for every set of simulations to $C_f(\lambda) = 0.3, 0.5, 0.65, 0.75, 0.85$, and 1, respectively.

After performing a set of $\times 100$ ideal Voigt profiles for each covering fraction, we introduce the effect of resolution by convolving the simulated spectra with a Gaussian kernel whose FWHM matches the instrumental resolution (a.k.a., $R$ (VIMOS) $\approx 600$). A constant S/N value ($S/N = 5$, the median of the VANDELS sample) is then added to every single spectra. At fixed $C_f(\lambda)$, the median measured (observed) residual flux of each mock distribution is compared to the theoretical line residual flux of $R_\lambda = 1 - C_f(\lambda)$ according to the picket-fence configuration, producing the desired calibration. Fig. A1 presents the resulting line profiles for a set of $\times 100$ simulations of the SiII 1260 line, all of them with an input covering fraction of $C_f(\lambda) = 0.85$ ($R_\lambda = 0.15$). The median of the distribution of individually measured $R_\lambda$ values is 0.3 (red line), showing that the observed line depth can statistically differ as much as 40 per cent with respect to the real value.

In the last step, we proceed to stack all the spectra contained on each simulation pack at fixed $C_f(\lambda)$, following the same methods described in the text. The co-added spectra for the $C_f(\lambda) = 0.85$ simulation is also shown in Fig. A1 (blue line), with a 65 per cent relative difference in residual flux with respect to the real input value.

To sum up, the resulting calibrations – fitting a linear polynomial to the measured versus theoretical $R_\lambda$ values – are:

$$R_{\text{SIII}} [\text{R} + \text{S/N}] = 0.68 \times (1 - C_f(\text{SiII})) + 0.20, \quad (A3)$$

and

$$R_{\text{SIII}} [\text{stacking}] = 0.73 \times (1 - C_f(\text{SiII})) + 0.29, \quad (A4)$$

when accounting for instrumental resolution plus S/N (equation A3), and spectral stacking (equation A4), respectively. Both relations are plotted in Fig. A2 through the dashed red and dashed blue lines. These calibrations were applied to individual residual flux measurements as well as when measuring the residual flux over composite spectra.

<footnote n="10">¹⁰NIST stands for National Institute of Standards and Technology: https://physics.nist.gov/PhysRefData/ASD/lines_form.html.</footnote>
<figure><img src="image_21.png" /></figure>

**Figure A1.** The impact of resolution and stacking on the depth of the absorption lines. Left-hand panel: the coloured lines represent a set of ×100 mock profiles of the Si II 1260 line with an input covering fraction of $C_f = 0.85$ (equivalent to $R_{SiII} = 1 - C_f = 0.15$, delimited by the horizontal black line), and where the Si II gas column density, $N_{SiII}(\text{cm}^{-2})$, the Doppler broadening parameter, $b(\text{kms}^{-1})$, and the gas velocity, $v(\text{kms}^{-1})$, were varied following Gaussian distributions (see the text). Right-hand panel: Distribution of individual $R_{SiII}$ measurements for all the simulations once the lines have been degraded to a spectral resolution of $R(\text{VIMOS}) = 600$ (red line), including a constant Gaussian noise of median S/N = 5. The blue line indicates the resulting residual flux after co-adding all the spectra (stacking).

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

**Figure A2.** Calibration between the theoretical residual flux of the Si II 1260 line and the actual measured flux when accounting for (1, in red) the effect of a finite instrumental resolution and S/N ($R(\text{VIMOS}) \approx 600$, S/N = 5), and (2, in blue) the joint effect of resolution, S/N plus stacking a set of N = 100 simulated lines (see Fig. A1). The dashed-black line shows the 1:1 relation.

## APPENDIX B: STACK MEASUREMENTS AND FICUS SED RESULTS

As stated in Section 3.3, the VANDELS sample was divided according to the 25th, 50th, and 75th distribution percentiles (quartiles) of the different targeted quantities, resulting in four subsamples for each quantity which were quoted as Q1, Q2, Q3, and Q4. Then, stacked spectra in bins of UV magnitude ($M_{1500}$), UV intrinsic luminosity ($M_{\text{int}}^{\text{int}}$), UV continuum slope ($\beta_{\text{spec}}^{\text{1500}}$), and Ly $\alpha$ equivalent width ($W_{\text{Ly}\alpha}$) were built following the methodology described in the same section.

Table B1 includes the ionizing escape fractions ($f_{\text{esc}}^{\text{abs}}$) and production efficiencies ($\xi_{\text{ion}}$) derived for every composite spectra. It also contains the same results when using either the R16 or SMC dust-attenuation laws, for comparison. Table B2 summarizes the outputs of our FICUS SED fits and the main UV spectral measurements for the LCEs, non-LCEs, LAEs, and non-LAEs stacks.
<figure><img src="image_21.png" /></figure>

**Figure A1.** The impact of resolution and stacking on the depth of the absorption lines. Left-hand panel: the coloured lines represent a set of ×100 mock profiles of the Si II 1260 line with an input covering fraction of $C_f = 0.85$ (equivalent to $R_{SiII} = 1 - C_f = 0.15$, delimited by the horizontal black line), and where the Si II gas column density, $N_{SiII}(\text{cm}^{-2})$, the Doppler broadening parameter, $b(\text{kms}^{-1})$, and the gas velocity, $v(\text{kms}^{-1})$, were varied following Gaussian distributions (see the text). Right-hand panel: Distribution of individual $R_{SiII}$ measurements for all the simulations once the lines have been degraded to a spectral resolution of $R(\text{VIMOS}) = 600$ (red line), including a constant Gaussian noise of median S/N = 5. The blue line indicates the resulting residual flux after co-adding all the spectra (stacking).

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

**Figure A2.** Calibration between the theoretical residual flux of the Si II 1260 line and the actual measured flux when accounting for (1, in red) the effect of a finite instrumental resolution and S/N ($R(\text{VIMOS}) \approx 600$, S/N = 5), and (2, in blue) the joint effect of resolution, S/N plus stacking a set of N = 100 simulated lines (see Fig. A1). The dashed-black line shows the 1:1 relation.

## APPENDIX B: STACK MEASUREMENTS AND FICUS SED RESULTS

As stated in Section 3.3, the VANDELS sample was divided according to the 25th, 50th, and 75th distribution percentiles (quartiles) of the different targeted quantities, resulting in four subsamples for each quantity which were quoted as Q1, Q2, Q3, and Q4. Then, stacked spectra in bins of UV magnitude ($M_{1500}$), UV intrinsic luminosity ($M_{\text{int}}^{\text{int}}$), UV continuum slope ($\beta_{\text{spec}}^{\text{1500}}$), and Ly $\alpha$ equivalent width ($W_{\text{Ly}\alpha}$) were built following the methodology described in the same section.

Table B1 includes the ionizing escape fractions ($f_{\text{esc}}^{\text{abs}}$) and production efficiencies ($\xi_{\text{ion}}$) derived for every composite spectra. It also contains the same results when using either the R16 or SMC dust-attenuation laws, for comparison. Table B2 summarizes the outputs of our FICUS SED fits and the main UV spectral measurements for the LCEs, non-LCEs, LAEs, and non-LAEs stacks.
<table><caption>**Table B1.** Inferred ionizing absolute escape fractions ($f_{\text{esc}}^{\text{abs}}$) and production efficiencies ($\xi_{\text{ion}}$) for the VANDELS composites in this work.</caption><thead><tr><th>Composite ID</th><th>Bin range</th><th>Median</th><th>N<sub>gal</sub></th><th>f<sub>esc</sub><sup>abs</sup> (R16, per cent)</th><th>ξ<sub>ion</sub> (R16, Hz erg<sup>-1</sup>)</th><th>f<sub>esc</sub><sup>abs</sup> (SMC, per cent)</th><th>ξ<sub>ion</sub> (SMC, Hz erg<sup>-1</sup>)</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">Stacks in bins of observed absolute UV magnitude (AB)</td></tr><tr><td>M<sub>1500</sub> (Q1)</td><td>≤-20.94</td><td>-21.22</td><td>63</td><td>1.70 ± 0.30</td><td>25.21 ± 0.04</td><td>3.35 ± 0.47</td><td>25.18 ± 0.04</td></tr><tr><td>M<sub>1500</sub> (Q2)</td><td>-20.94, -20.53</td><td>-20.70</td><td>71</td><td>1.01 ± 0.23</td><td>25.26 ± 0.05</td><td>2.19 ± 0.40</td><td>25.20 ± 0.06</td></tr><tr><td>M<sub>1500</sub> (Q3)</td><td>-20.53, -20.16</td><td>-20.34</td><td>73</td><td>1.55 ± 0.34</td><td>25.30 ± 0.09</td><td>3.11 ± 0.62</td><td>25.31 ± 0.08</td></tr><tr><td>M<sub>1500</sub> (Q4)</td><td>>-20.16</td><td>-19.95</td><td>90</td><td>2.38 ± 0.33</td><td>25.42 ± 0.09</td><td>4.09 ± 0.56</td><td>25.43 ± 0.08</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">Stacks in bins of intrinsic UV luminosity (AB)</td></tr><tr><td>M<sub>1500</sub><sup>int</sup> (Q1)</td><td>≤-23.49</td><td>-23.91</td><td>74</td><td>0.43 ± 0.10</td><td>25.20 ± 0.03</td><td>1.39 ± 0.28</td><td>25.09 ± 0.04</td></tr><tr><td>M<sub>1500</sub><sup>int</sup> (Q2)</td><td>-23.49, -22.86</td><td>-23.16</td><td>76</td><td>1.02 ± 0.22</td><td>25.26 ± 0.08</td><td>2.45 ± 0.52</td><td>25.16 ± 0.08</td></tr><tr><td>M<sub>1500</sub><sup>int</sup> (Q3)</td><td>-22.86, -22.29</td><td>-22.61</td><td>71</td><td>2.20 ± 0.57</td><td>25.29 ± 0.09</td><td>3.78 ± 0.67</td><td>25.40 ± 0.08</td></tr><tr><td>M<sub>1500</sub><sup>int</sup> (Q4)</td><td>>-22.29</td><td>-21.85</td><td>71</td><td>5.64 ± 1.19</td><td>25.28 ± 0.05</td><td>7.21 ± 1.28</td><td>25.28 ± 0.05</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">Stacks in bins of UV continuum slope</td></tr><tr><td>β<sub>spec</sub><sup>1500</sup> (Q1)</td><td>≤-1.65</td><td>-1.82</td><td>69</td><td>6.37 ± 1.21</td><td>25.48 ± 0.05</td><td>9.31 ± 1.70</td><td>25.44 ± 0.06</td></tr><tr><td>β<sub>spec</sub><sup>1500</sup> (Q2)</td><td>-1.65, -1.34</td><td>-1.49</td><td>78</td><td>1.94 ± 0.46</td><td>25.27 ± 0.06</td><td>3.45 ± 0.56</td><td>25.28 ± 0.07</td></tr><tr><td>β<sub>spec</sub><sup>1500</sup> (Q3)</td><td>-1.34, -1.03</td><td>-1.20</td><td>74</td><td>0.97 ± 0.24</td><td>25.18 ± 0.04</td><td>2.16 ± 0.40</td><td>25.14 ± 0.07</td></tr><tr><td>β<sub>spec</sub><sup>1500</sup> (Q4)</td><td>>-1.03</td><td>-0.81</td><td>76</td><td>0.27 ± 0.09</td><td>25.12 ± 0.05</td><td>0.87 ± 0.20</td><td>24.99 ± 0.05</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">Stacks in bins of Ly α equivalent width (Å)</td></tr><tr><td>-1 × W<sub>Lyα</sub> (Q1)</td><td>≤-3.61</td><td>-5.86</td><td>71</td><td>0.60 ± 0.17</td><td>25.18 ± 0.02</td><td>1.44 ± 0.31</td><td>25.13 ± 0.04</td></tr><tr><td>-1 × W<sub>Lyα</sub> (Q2)</td><td>-3.61, 4.58</td><td>1.13</td><td>67</td><td>1.14 ± 0.23</td><td>25.14 ± 0.06</td><td>2.38 ± 0.40</td><td>25.14 ± 0.08</td></tr><tr><td>-1 × W<sub>Lyα</sub> (Q3)</td><td>4.58, 15.60</td><td>7.75</td><td>79</td><td>1.44 ± 0.37</td><td>25.13 ± 0.07</td><td>2.95 ± 0.55</td><td>25.11 ± 0.07</td></tr><tr><td>-1 × W<sub>Lyα</sub> (Q4)</td><td>>15.60</td><td>30.42</td><td>80</td><td>4.24 ± 0.92</td><td>25.41 ± 0.07</td><td>7.09 ± 1.27</td><td>25.44 ± 0.07</td></tr></tbody></table>

*Note.* Column 1: composite identifier. Column 2: interquartile range for each magnitude bin (units are in the spanning titles). Column 3: median of the interquartile range. Column 4: number of objects included in the stacks. Columns 5 and 6: estimated escape fractions (in per cent) and production efficiencies (in Hz erg⁻¹) when using R16 dust-attenuation law. Columns 7 and 8: estimated escape fractions (in per cent) and production efficiencies (in Hz erg⁻¹) when using SMC dust-attenuation law.

<table><caption>**Table B2.** Summary chart containing the main spectral measurements and FICUS SED fitting results for LCEs, non-LCEs, LAEs, and non-LAEs composites.</caption><thead><tr><th>Parameter <sup>a</sup></th><th>Units</th><th>LCEs<br>(f<sub>esc</sub><sup>abs</sup> ≥ 5per cent)</th><th>non-LCEs<br>(f<sub>esc</sub><sup>abs</sup> < 5per cent)</th><th>LAEs (W<sub>Lyα</sub> ≤ -20 Å)</th><th>non-LAEs (W<sub>Lyα</sub> > -20 Å)</th></tr></thead><tbody><tr><td>W<sub>Lyα</sub></td><td>Å</td><td>-29.71 ± 2.46</td><td>-3.62 ± 0.41</td><td>-33.45 ± 1.10</td><td>-1.63 ± 0.22</td></tr><tr><td>β<sub>spec</sub><sup>1500</sup></td><td>-</td><td>-2.17 ± 0.03</td><td>-1.23 ± 0.01</td><td>-1.68 ± 0.03</td><td>-1.22 ± 0.01</td></tr><tr><td>W<sub>LIS</sub></td><td>Å</td><td>0.86 ± 0.07</td><td>1.90 ± 0.02</td><td>1.04 ± 0.05</td><td>2.01 ± 0.03</td></tr><tr><td>R<sub>LIS</sub><sup>b</sup></td><td>-</td><td>0.73 ± 0.06</td><td>0.47 ± 0.02</td><td>0.64 ± 0.05</td><td>0.43 ± 0.03</td></tr><tr><td>W<sub>C IV 1550</sub></td><td>Å</td><td>-0.72 ± 0.18</td><td>-0.26 ± 0.05</td><td>-0.55 ± 0.10</td><td>-0.25 ± 0.05</td></tr><tr><td>W<sub>He II 1640</sub></td><td>Å</td><td>-1.11 ± 0.25</td><td>-0.81 ± 0.06</td><td>-1.16 ± 0.15</td><td>-0.76 ± 0.06</td></tr><tr><td>W<sub>C III 1908</sub></td><td>Å</td><td>-2.15 ± 0.67</td><td>-0.73 ± 0.13</td><td>-2.28 ± 0.31</td><td>-0.46 ± 0.14</td></tr><tr><td>C IV/C III</td><td>-</td><td>0.79 ± 0.21</td><td>0.42 ± 0.20</td><td>-</td><td>-</td></tr><tr><td>Age</td><td>Myr</td><td>13.60 ± 2.73</td><td>19.14 ± 0.92</td><td>13.63 ± 1.97</td><td>20.26 ± 0.90</td></tr><tr><td>Z<sub>*</sub></td><td>Z<sub>☉</sub></td><td>0.18 ± 0.02</td><td>0.21 ± 0.01</td><td>0.19 ± 0.02</td><td>0.21 ± 0.01</td></tr><tr><td>E<sub>B-V</sub></td><td>mag.</td><td>0.07 ± 0.01</td><td>0.23 ± 0.00</td><td>0.16 ± 0.00</td><td>0.24 ± 0.00</td></tr><tr><td>(F<sub>λ900</sub>/F<sub>λ1500</sub>)<sub>int</sub></td><td>-</td><td>0.64 ± 0.09</td><td>0.47 ± 0.02</td><td>0.69 ± 0.07</td><td>0.43 ± 0.02</td></tr><tr><td>ξ<sub>ion</sub></td><td>Hz erg<sup>-1</sup></td><td>25.38 ± 0.08</td><td>25.18 ± 0.03</td><td>25.41 ± 0.07</td><td>25.14 ± 0.03</td></tr></tbody></table>

*Notes.* Column 1: galaxy properties. Column 2: physical units.

Columns 3 and 4: list of measurements for the LCEs and non-LCEs stacks.

<sup>a</sup>All tabulated results use the R16 dust-attenuation law.

<sup>b</sup>Measured residual fluxes not corrected by the spectral resolution (see Appendix A).

## APPENDIX C: ALTERNATIVE ESCAPE FRACTION ESTIMATES

Here we offer complementary predictions for the ionizing photon escape fraction. $f_{\text{esc}}^{\text{abs}}$ is derived by combining absorption line measurements with estimates of the UV attenuation, in a similar manner as in the main text (equation 6). However, the dust-attenuation parameter ($E_{\text{B-V}}$) has now been estimated from the UV slope measurements of Calabrò et al. (2021), by fit-

ting a power law to the VANDELS photometry, requiring the whole bandwidths to reside between 1230 and 2750 Å rest-frame. equation (8) in Chisholm et al. (2022) was then used to convert $\beta_{1500}$ to $E_{\text{B-V}}$ assuming an SMC attenuation law. VANDELS, KLCS, and LzLCS $f_{\text{esc}}^{\text{abs}}$ correlations with physical quantities converge to the same global trends (see Fig. C1), with $f_{\text{esc}}^{\text{abs}}$ increasing towards fainter UV magnitudes, bluer UV slopes, and stronger Ly$\alpha$ emission. We refer to Section 6 for more details.
<table><caption>**Table B1.** Inferred ionizing absolute escape fractions ($f_{esc}^{abs}$) and production efficiencies ($\xi_{ion}$) for the VANDELS composites in this work.</caption><thead><tr><th>Composite ID</th><th>Bin range</th><th>Median</th><th>N<sub>gal</sub></th><th>f<sub>esc</sub><sup>abs</sup> (R16, per cent)</th><th>ξ<sub>ion</sub> (R16, Hz erg<sup>-1</sup>)</th><th>f<sub>esc</sub><sup>abs</sup> (SMC, per cent)</th><th>ξ<sub>ion</sub> (SMC, Hz erg<sup>-1</sup>)</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">Stacks in bins of observed absolute UV magnitude (AB)</td></tr><tr><td>M<sub>1500</sub> (Q1)</td><td>≤-20.94</td><td>-21.22</td><td>63</td><td>1.70 ± 0.30</td><td>25.21 ± 0.04</td><td>3.35 ± 0.47</td><td>25.18 ± 0.04</td></tr><tr><td>M<sub>1500</sub> (Q2)</td><td>-20.94, -20.53</td><td>-20.70</td><td>71</td><td>1.01 ± 0.23</td><td>25.26 ± 0.05</td><td>2.19 ± 0.40</td><td>25.20 ± 0.06</td></tr><tr><td>M<sub>1500</sub> (Q3)</td><td>-20.53, -20.16</td><td>-20.34</td><td>73</td><td>1.55 ± 0.34</td><td>25.30 ± 0.09</td><td>3.11 ± 0.62</td><td>25.31 ± 0.08</td></tr><tr><td>M<sub>1500</sub> (Q4)</td><td>>-20.16</td><td>-19.95</td><td>90</td><td>2.38 ± 0.33</td><td>25.42 ± 0.09</td><td>4.09 ± 0.56</td><td>25.43 ± 0.08</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">Stacks in bins of intrinsic UV luminosity (AB)</td></tr><tr><td>M<sub>1500</sub><sup>int</sup> (Q1)</td><td>≤-23.49</td><td>-23.91</td><td>74</td><td>0.43 ± 0.10</td><td>25.20 ± 0.03</td><td>1.39 ± 0.28</td><td>25.09 ± 0.04</td></tr><tr><td>M<sub>1500</sub><sup>int</sup> (Q2)</td><td>-23.49, -22.86</td><td>-23.16</td><td>76</td><td>1.02 ± 0.22</td><td>25.26 ± 0.08</td><td>2.45 ± 0.52</td><td>25.16 ± 0.08</td></tr><tr><td>M<sub>1500</sub><sup>int</sup> (Q3)</td><td>-22.86, -22.29</td><td>-22.61</td><td>71</td><td>2.20 ± 0.57</td><td>25.29 ± 0.09</td><td>3.78 ± 0.67</td><td>25.40 ± 0.08</td></tr><tr><td>M<sub>1500</sub><sup>int</sup> (Q4)</td><td>>-22.29</td><td>-21.85</td><td>71</td><td>5.64 ± 1.19</td><td>25.28 ± 0.05</td><td>7.21 ± 1.28</td><td>25.28 ± 0.05</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">Stacks in bins of UV continuum slope</td></tr><tr><td>β<sub>spec</sub><sup>1500</sup> (Q1)</td><td>≤-1.65</td><td>-1.82</td><td>69</td><td>6.37 ± 1.21</td><td>25.48 ± 0.05</td><td>9.31 ± 1.70</td><td>25.44 ± 0.06</td></tr><tr><td>β<sub>spec</sub><sup>1500</sup> (Q2)</td><td>-1.65, -1.34</td><td>-1.49</td><td>78</td><td>1.94 ± 0.46</td><td>25.27 ± 0.06</td><td>3.45 ± 0.56</td><td>25.28 ± 0.07</td></tr><tr><td>β<sub>spec</sub><sup>1500</sup> (Q3)</td><td>-1.34, -1.03</td><td>-1.20</td><td>74</td><td>0.97 ± 0.24</td><td>25.18 ± 0.04</td><td>2.16 ± 0.40</td><td>25.14 ± 0.07</td></tr><tr><td>β<sub>spec</sub><sup>1500</sup> (Q4)</td><td>>-1.03</td><td>-0.81</td><td>76</td><td>0.27 ± 0.09</td><td>25.12 ± 0.05</td><td>0.87 ± 0.20</td><td>24.99 ± 0.05</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">Stacks in bins of Ly α equivalent width (Å)</td></tr><tr><td>-1 × W<sub>Lyα</sub> (Q1)</td><td>≤-3.61</td><td>-5.86</td><td>71</td><td>0.60 ± 0.17</td><td>25.18 ± 0.02</td><td>1.44 ± 0.31</td><td>25.13 ± 0.04</td></tr><tr><td>-1 × W<sub>Lyα</sub> (Q2)</td><td>-3.61, 4.58</td><td>1.13</td><td>67</td><td>1.14 ± 0.23</td><td>25.14 ± 0.06</td><td>2.38 ± 0.40</td><td>25.14 ± 0.08</td></tr><tr><td>-1 × W<sub>Lyα</sub> (Q3)</td><td>4.58, 15.60</td><td>7.75</td><td>79</td><td>1.44 ± 0.37</td><td>25.13 ± 0.07</td><td>2.95 ± 0.55</td><td>25.11 ± 0.07</td></tr><tr><td>-1 × W<sub>Lyα</sub> (Q4)</td><td>>15.60</td><td>30.42</td><td>80</td><td>4.24 ± 0.92</td><td>25.41 ± 0.07</td><td>7.09 ± 1.27</td><td>25.44 ± 0.07</td></tr></tbody></table>

*Note.* Column 1: composite identifier. Column 2: interquartile range for each magnitude bin (units are in the spanning titles). Column 3: median of the interquartile range. Column 4: number of objects included in the stacks. Columns 5 and 6: estimated escape fractions (in per cent) and production efficiencies (in Hz erg⁻¹) when using R16 dust-attenuation law. Columns 7 and 8: estimated escape fractions (in per cent) and production efficiencies (in Hz erg⁻¹) when using SMC dust-attenuation law.

<table><caption>**Table B2.** Summary chart containing the main spectral measurements and FICUS SED fitting results for LCEs, non-LCEs, LAEs, and non-LAEs composites.</caption><thead><tr><th>Parameter <sup>a</sup></th><th>Units</th><th>LCEs<br>(f<sub>esc</sub><sup>abs</sup> ≥ 5per cent)</th><th>non-LCEs<br>(f<sub>esc</sub><sup>abs</sup> < 5per cent)</th><th>LAEs (W<sub>Lyα</sub> ≤ -20 Å)</th><th>non-LAEs (W<sub>Lyα</sub> > -20 Å)</th></tr></thead><tbody><tr><td>W<sub>Lyα</sub></td><td>Å</td><td>-29.71 ± 2.46</td><td>-3.62 ± 0.41</td><td>-33.45 ± 1.10</td><td>-1.63 ± 0.22</td></tr><tr><td>β<sub>spec</sub><sup>1500</sup></td><td>-</td><td>-2.17 ± 0.03</td><td>-1.23 ± 0.01</td><td>-1.68 ± 0.03</td><td>-1.22 ± 0.01</td></tr><tr><td>W<sub>LIS</sub></td><td>Å</td><td>0.86 ± 0.07</td><td>1.90 ± 0.02</td><td>1.04 ± 0.05</td><td>2.01 ± 0.03</td></tr><tr><td>R<sub>LIS</sub><sup>b</sup></td><td>-</td><td>0.73 ± 0.06</td><td>0.47 ± 0.02</td><td>0.64 ± 0.05</td><td>0.43 ± 0.03</td></tr><tr><td>W<sub>C IV 1550</sub></td><td>Å</td><td>-0.72 ± 0.18</td><td>-0.26 ± 0.05</td><td>-0.55 ± 0.10</td><td>-0.25 ± 0.05</td></tr><tr><td>W<sub>He II 1640</sub></td><td>Å</td><td>-1.11 ± 0.25</td><td>-0.81 ± 0.06</td><td>-1.16 ± 0.15</td><td>-0.76 ± 0.06</td></tr><tr><td>W<sub>C III 1908</sub></td><td>Å</td><td>-2.15 ± 0.67</td><td>-0.73 ± 0.13</td><td>-2.28 ± 0.31</td><td>-0.46 ± 0.14</td></tr><tr><td>C IV/C III</td><td>-</td><td>0.79 ± 0.21</td><td>0.42 ± 0.20</td><td>-</td><td>-</td></tr><tr><td>Age</td><td>Myr</td><td>13.60 ± 2.73</td><td>19.14 ± 0.92</td><td>13.63 ± 1.97</td><td>20.26 ± 0.90</td></tr><tr><td>Z<sub>*</sub></td><td>Z<sub>☉</sub></td><td>0.18 ± 0.02</td><td>0.21 ± 0.01</td><td>0.19 ± 0.02</td><td>0.21 ± 0.01</td></tr><tr><td>E<sub>B-V</sub></td><td>mag.</td><td>0.07 ± 0.01</td><td>0.23 ± 0.00</td><td>0.16 ± 0.00</td><td>0.24 ± 0.00</td></tr><tr><td>(F<sub>λ900</sub>/F<sub>λ1500</sub>)<sub>int</sub></td><td>-</td><td>0.64 ± 0.09</td><td>0.47 ± 0.02</td><td>0.69 ± 0.07</td><td>0.43 ± 0.02</td></tr><tr><td>ξ<sub>ion</sub></td><td>Hz erg<sup>-1</sup></td><td>25.38 ± 0.08</td><td>25.18 ± 0.03</td><td>25.41 ± 0.07</td><td>25.14 ± 0.03</td></tr></tbody></table>

*Notes.* Column 1: galaxy properties. Column 2: physical units.

Columns 3 and 4: list of measurements for the LCEs and non-LCEs stacks.

<sup>a</sup>All tabulated results use the R16 dust-attenuation law.

<sup>b</sup>Measured residual fluxes not corrected by the spectral resolution (see Appendix A).

## APPENDIX C: ALTERNATIVE ESCAPE FRACTION ESTIMATES

Here we offer complementary predictions for the ionizing photon escape fraction. $f_{esc}^{abs}$ is derived by combining absorption line measurements with estimates of the UV attenuation, in a similar manner as in the main text (equation 6). However, the dust-attenuation parameter ($E_{B-V}$) has now been estimated from the UV slope measurements of Calabrò et al. (2021), by fit-

ting a power law to the VANDELS photometry, requiring the whole bandwidths to reside between 1230 and 2750 Å rest-frame. equation (8) in Chisholm et al. (2022) was then used to convert $\beta_{1500}$ to $E_{B-V}$ assuming an SMC attenuation law. VANDELS, KLCS, and LzLCS $f_{esc}^{abs}$ correlations with physical quantities converge to the same global trends (see Fig. C1), with $f_{esc}^{abs}$ increasing towards fainter UV magnitudes, bluer UV slopes, and stronger Ly$\alpha$ emission. We refer to Section 6 for more details.
<figure><img src="image_23.png" /></figure>

**Figure C1.** The ionizing photon escape fraction ($f_{esc}^{abs}$) versus the UV magnitude ($M_{1500}$), the UV slope at 1500 Å ($\beta_{1500}$), and the Ly $\alpha$ equivalent width ($W_{Ly\alpha}$). Blue dots represent the predicted $f_{esc}^{abs}$ values for VANDELS galaxies (equation 6, using SMC) when the dust-attenuation ($E_{B-V}$) is inferred from the $\beta_{1500}$ measurement by Calabrò et al. (2021). Layout is the same as Fig. 15.

This paper has been typeset from a TEX/LATEX file prepared by the author.
<figure><img src="image_23.png" /></figure>

**Figure C1.** The ionizing photon escape fraction ($f_{esc}^{abs}$) versus the UV magnitude ($M_{1500}$), the UV slope at 1500 Å ($\beta_{1500}$), and the Ly $\alpha$ equivalent width ($W_{Ly\alpha}$). Blue dots represent the predicted $f_{esc}^{abs}$ values for VANDELS galaxies (equation 6, using SMC) when the dust-attenuation ($E_{B-V}$) is inferred from the $\beta_{1500}$ measurement by Calabrò et al. (2021). Layout is the same as Fig. 15.

This paper has been typeset from a TEX/LATEX file prepared by the author.