Comparing the VANDELS Sample to a Zoom-in Radiative Hydrodynamical Simulation: Using the Si ii and C ii Line Spectra as Tracers of Galaxy Evolution and Lyman Continuum Leakage

We compare mock ultraviolet C ii and Si ii absorption and emission line features generated using a ∼109 M ⊙ virtual galaxy with observations of 131 z ∼ 3 galaxies from the vandels survey. We find that the mock spectra reproduce reasonably well a large majority (83%) of the vandels spectra (χ 2 < 2), but do not resemble the most massive objects (⪆1010 M ⊙), which exhibit broad absorption features. Interestingly, the best-matching mock spectra originate from periods of intense star formation in the virtual galaxy, where its luminosity is 4 times higher than in periods of relative quiescence. Furthermore, for each galaxy, we predict the Lyman continuum (LyC) escape fractions ( fesc(pred)LyC ) using the environment of the virtual galaxy. We derive an average fesc(pred)LyC of 0.01 ± 0.02, consistent with other estimates from the literature. The fesc(pred)LyC are tightly correlated with the Lyα escape fractions and highly consistent with observed empirical trends. Additionally, galaxies with larger fesc(pred)LyC exhibit bluer β slopes, more Lyα flux, and weaker low-ionization absorption lines. Building upon the good agreement between fesc(pred)LyC and observationally established LyC diagnostics, we examine the LyC leakage mechanisms in the simulation. We find that LyC photon leakage is enhanced in directions where the observed flux dominantly emerges from compact regions depleted of neutral gas and dust, mirroring the scenario inferred from observational data. In general, this study further highlights the potential of high-resolution radiation hydrodynamics simulations in analyzing UV absorption and emission line features and providing valuable insights into the LyC leakage of star-forming galaxies.


INTRODUCTION
Ultraviolet (UV) absorption and emission line features are formidable tracers of the Interstellar and Circumgalactic Media (ISM, CGM) properties within starforming galaxies.These UV features provide direct insights into key gas properties, including its metallicity, geometry, and kinematics.Such parameters are fundamental in tracking the evolving environment of galaxies over cosmic epochs, forming a cornerstone for the establishment of precise galaxy evolution models.In particular, the analysis of UV features enables us to explore properties intrinsically linked to the baryon cycle of galaxies, quantify the metal content of neutral gas, and infer the metal enrichment of the ISM (e.g., Heckman et al. 2000;James et al. 2014;Heckman et al. 2015;Shapley et al. 2003;Erb et al. 2012;Rubin et al. 2011;Wofford et al. 2013;Chisholm et al. 2015aChisholm et al. , 2017;;Finley et al. 2017;Steidel et al. 2018;Wang et al. 2020;Xu et al. 2022).
UV absorption lines originating from low-ionization states (LIS) of metals, such as C II or Si II, hold particular relevance in the context of reionization investigations.These absorption features, characterized by substantial residual flux (i.e., the flux remaining at the line's maximum depth) or relatively low equivalent widths (EWs), have been identified in various galaxy samples exhibiting ionizing photon leakage (referred to as Lyman Continuum Emitters or LCEs, Reddy et al. 2016;Gazagnes et al. 2018;Saldana-Lopez et al. 2022).Indeed, these features directly probe the absence of dense neutral clouds or the presence of low-column density channels through which ionizing photons can escape the galaxy's ISM (e.g., Heckman et al. 2001;Alexandroff et al. 2015;Reddy et al. 2016;Gazagnes et al. 2018;Chisholm et al. 2018;Steidel et al. 2018;Saldana-Lopez et al. 2022).This property underscores the potential of LIS absorption lines as promising diagnostic tools for predicting the escape fraction of ionizing photons (f LyC esc , where LyC represents the Lyman Continuum photons with λ < 912 Å) in the high-redshift Universe, where direct measurements are hindered by the opacity of the intergalactic medium (IGM).
While UV LIS absorption and emission lines offer substantial potential for exploring the characteristics of a galaxy's ISM and constraining the LyC photons escaping star-forming galaxies, conventional analyses and diagnostics of these features often rely on spectrally averaged measurements, such as equivalent width or residual flux, which merely capture a fraction of the rich informational potential.In practice, the emergent profiles of UV absorption and emission lines can exhibit a high degree of complexity resulting from the interplay of multiple gas clouds.The ISM geometry and complexity are critical keys to understanding the physical processes that have transformed the galaxy environment and might have enabled the leakage of LyC photons (e.g., Jaskot et al. 2019;Rivera-Thorsen et al. 2017).To capture the vital features that are directly related to certain properties such as the LyC escape, one must imperatively advance our analytical methodologies to fully exploit the amount of information available within these features.
Over the past decades, semi-analytical techniques building upon outflow and inflow geometrical models have shown promising results for interpreting simultaneously and consistently the resonant and fluorescent features of low and high ionization states of metals (e.g., Scarlata & Panagia 2015;Carr et al. 2018;Carr & Scarlata 2022;Carr et al. 2023;Li et al. 2024).Simultaneously, the development and progress of hydrodynamical simulation techniques have made a significant step in the development of new frameworks for the interpretation of spectroscopic observations.Building upon Mauerhofer et al. (2021), whose authors used a zoom-in high-resolution radiation hydrodynamics (RHD) simulation post-processed with radiativetransfer techniques, Blaizot et al. (2023) highlighted that the mock Lyman−α (Lyα) profiles produced in such a simulation closely resemble some of the notoriously complex Lyα profiles observed in the low−z universe.Gazagnes et al. (2023) showed that the same simulation produced mock C II and Si II spectra (including the fluorescent C II* and Si II* features) which could well reproduce 38 out of 45 observations of low−z galaxies from the COS Legacy Archive Survey (classy, Berg et al. 2022).This strong correspondence between UV absorption and emission lines from star-forming galaxies and simulated spectra generated from a single virtual galaxy has opened a promising avenue for interpreting galaxy observations.It suggests that we can now construct physically meaningful models that faithfully replicate realistic environments and help us connect observations to the actual physical galaxy properties responsible for the emergent spectral properties (see also Choustikov et al. 2023;Katz et al. 2023).
In this work, we extend our exploration of how mock observations from high-resolution simulations can offer a novel and insightful approach to connect spectroscopic observations to the underlying physical processes governing galaxy evolution.Our investigation dives into two distinct aspects.First, while prior work by Blaizot et al. (2023) and Gazagnes et al. (2023) focused on observations at redshifts around z < 0.5, we assess the simulation's reliability in replicating the C II and Si II spectra of a more extensive galaxy sample at 3.35 < z < 3.95, drawn from the ultra-deep VANDELS spectroscopic survey (McLure et al. 2018;Pentericci et al. 2018;Garilli et al. 2021).Secondly, we investigate the origin of the similarities between observed and simulated spectra and employ the simulation as a predictive tool for Lyman Continuum escape (f LyC esc ), comparing these predictions with both direct and indirect f LyC esc constraints from Begley et al. (2022) and Saldana-Lopez et al. (2023), as well as with the Lyα escape fraction measurements from Begley et al. (2024).One objective of this paper is to further illustrate the potential of high-resolution RHD simulations in the interpretation of UV spectroscopic observations and in the extraction of meaningful LyC diagnostics, an aspect that holds particular relevance in the current landscape of high-redshift observations.This paper is organized as follows: Section 2 describes the simulation and observation data sets used in this work.Section 3 compares the C II λ1334+C II* λ1335 and Si II λ1260+Si II* λ1265 spectra in the simulation and in the observations.Section 4 discusses the origin of the resemblance between the mock and observed spectra in the context of the virtual galaxy properties.Section 5 presents f LyC esc predictions made using the simulation and explores their consistency with previous LyC constraints and the shape of UV spectral features.Finally, Section 6 summarizes our findings.

DATA
In this work, we compare the C II λ1334+C II* λ1335 and Si II λ1260+Si II* λ1265 observations from vandels to mock spectra generated from the simulated galaxy of Mauerhofer et al. (2021).Section 2.1 describes the simulation and Section 2.2 details the VANDELS observations.

Zoom-in simulation
The RHD simulation used in the work is identical to the one used in Gazagnes et al. (2023) where we analyzed the C II λ1334+C II* λ1335 and Si II λ1260+Si II* λ1265 observations from the classy sample (Berg et al. 2022).Therefore, we only present here the most important characteristics of the simulation and refer the readers to Section 2 of Gazagnes et al. (2023) for a complete and thorough description.
The zoom-in simulation run was originally presented in Mauerhofer et al. (2021), and was obtained with the adaptive mesh refinement code Ramses-RT (Teyssier 2002;Rosdahl et al. 2013;Rosdahl & Teyssier 2015) which builds upon the physics of cooling, supernova feedback, and star-formation from sphinx (Rosdahl et al. 2018).The initial conditions were generated using music (Hahn & Abel 2013) and chosen such that the resulting galaxy has a stellar mass of M ⋆ ∼ 10 9 M ⊙ at z = 3, with a maximum cell resolution of 14 pc around this redshift.We selected 75 outputs with 10 Myrs timesteps between z = 4.19 to z = 3.00.During this period, the halo mass increases from 10 10.6 to 10 10.75 M ⊙ , the stellar mass (M ⋆ ) increases from 10 8.8 to 10 9.4 M ⊙ , both the averaged stellar and gas metallicity increase from 0.20Z ⊙ to 0.42Z ⊙ , and the star-formation rate (SFR) oscillates between ∼ 1.5 and 4.0 M ⊙ yr −1 , with two main SFR peaks around 1.64 and 2.04 Gyrs (the evolution of the SFR is shown and discussed in Section 4.2).
To produce the LIS mock spectra, we first compute the number densities of C II and Si II using the following procedure: we derive the densities of carbon and silicon in each simulation cell assuming solar abundance ratios, and then compute their ionization fractions with the customizable chemical network code krome (Grassi et al. 2014), using the recombination rates from (Badnell 2006), the collisional ionization rates from (Voronov 1997), the photoionization cross-sections from (Verner et al. 1996), and assuming chemical equilibrium for both ions.We also include the depletion fraction of metals on dust as a function of gas-phase metallicity from De Cia et al. ( 2016) and Konstantopoulou et al. (2022) and apply the correction formula cell by cell in the simulation.
The ground-state energy level of the C II and Si II ions is split into two spin states having an energy differ-ence of 0.0079 eV for C II and 0.0356 eV for Si II.These fine-structure levels produce two additional channels (referred to as fluorescent channels, noted as C II* and Si II*) of photon absorption and emission at 1335.66 Å and 1335.71Å, respectively, for C II, and at 1264.73 Å and 1265.02Å, respectively, for Si II.At typical ISM temperatures, nearly all of the C and Si will reside in the ground state.Under these conditions, large radiative excitation through the 1260 and 1334 Å channels leads to radiative de-excitation through either the resonant or the fluorescent channel.Since the densities of photons and particles are low, the fine structure state will radiatively decay back down to the ground state.In other words, this means that while the ground state might be optically thick to radiation, the fluorescent line can remain optically thin to radiative excitation such that fluorescent photons are not radiatively trapped as easily as in the resonant channel.
Similarly to Mauerhofer et al. (2021) and Gazagnes et al. (2023), we consider both the resonant and fluorescent channels of excitation and de-excitation when generating the LIS mock spectra.The final simulated line spectra are created using the radiative transfer postprocessing code rascas (Michel-Dansac et al. 2020) using the following procedure: a total of one million photon packets sample the spectra of stellar populations contained in the star particles, within a wavelength range of ± 10 Å around the C II λ1334 and Si II λ1260 lines.Photon packets are emitted from star particles and travel through the simulation grid.In each cell the optical depth is the sum of all channels of interaction: The optical depth of each line is the product of the ion column density in each cell and the ion cross section which is defined by Eq (2) in Mauerhofer et al. (2021).The line parameters are taken from the NIST atomic database (Kramida et al. 2018) 1 .We compute the turbulent velocity of the gas (v turb ) accounting for the gas density and velocity in the neighboring cells, as in the star-formation recipe of e.g., Trebitsch et al. (2021).τ dust is set following the implementation of Katz et al. (2022) where the dust-to-metal ratio is a double power-law function of the metallicity and is set to 0 in cells with temperatures above 10 5 K (inspired from the results of Rémy-Ruyer et al. 2014).
We implement the radiative transfer post-processing method within a spherical volume with a radius equivalent to three times the galaxy's virial radius (R vir ).Within this volume, we generate a set of 22,500 mock spectra for C II and Si II lines.These 22,500 mock spectra correspond to 300 distinct lines of sight of observations applied to each of the 75 simulation outputs.The 300 directions, obtained with HealPix, are uniformly sampled across the sphere and are the same for each simulation output.We employ an aperture size of 1" diameter (∼7 kpc at z = 4.2 and 7.9 kpc at z = 3).

Observations
In this work, we compare the mock spectra described in Section 2.1 to observations from the vandels survey (McLure et al. 2018;Pentericci et al. 2018;Garilli et al. 2021).The vandels survey was conducted using the Visible Multi-Object Spectrograph (VIMOS) on the ESO Very Large Telescope (VLT), featuring a median spectral resolution of R = 580.All vandels observations were taken with the MR grism+GG475 order sorting filter with a 1" slit width and a minimum slit length of 7".This extensive survey comprises deep spectra of 2087 galaxies, covering observed wavelengths ranging from 4,800 Å to 10,000 Å.The survey encompasses regions in the Chandra Deep Field South (CDFS) and UKIDSS Ultra Deep Survey (UDS) fields.Notably, the majority of objects included in the vandels dataset consist of star-forming galaxies within the redshift range of 2.4 ≤ z ≤ 7.0 such that VIMOS provides coverage of their rest-frame far-ultraviolet (FUV) spectra.
For this work, we focus on the analysis of the C II λ1334 and Si II λ1260 absorption lines (including the nearby C II* and Si II* fluorescent features).We pre-select a first sample of 965 galaxies at z > 2.8 to ensure that both the Si II and C II features are covered.We note that all the selected galaxies have high-quality spectroscopic redshifts (vandels quality flags of 3, 4, or 9 which means the redshifts are reliable at > 95%, Garilli et al. 2021).Nevertheless, in some cases, the redshift estimations are based on either Lyman−α or ISM lines so the exact systemic redshift might differ by a few hundred km s −1 (less or equal to the spectroscopic resolution element).We further detailed how we accounted for possible line shifts in Section 3.2 where we directly compare the vandels and mock spectra line profiles.
From this first sample of 965 galaxies with both Si II and C II observed, we select spectra that have a signalto-noise ratio (S/N) of at least 5 per resolution element.We chose a high S/N cut because vandels spectra have relatively low resolution, so the number of spectral bins that can be directly compared between the simulated  and observed spectra is small.Selecting high S/N spectra ensures that observations have reduced noise features, reducing the chance of biasing our comparisons.We visually checked each vandels spectrum to ensure there was no interloper contamination around the C II and Si II lines.We ended up with a final set of 131 galaxies.Figure 1 illustrates the redshift, M ⋆ , SFR, and absolute AB magnitude at 1250 Å (M 1250 ) ranges of our final vandels sub-sample alongside those of the virtual galaxy, considering its evolution from z = 4.19 to z = 3.00.In Gazagnes et al. (2023), we compared the current simulation to a sample of z < 0.18 galaxies with a broad range of stellar mass (∼ 10 6 − 10 10 M ⊙ ).
In contrast, the vandels sample better matches the redshift range of the virtual object, and its range of M ⋆ is slightly larger than that of our virtual galaxy (∼ 10 9 − 10 10.5 M ⊙ ).
The SFR and M 1250 of the vandels galaxies are significantly larger than our virtual galaxy.Yet, in Gazagnes et al. (2023), we showed that the simulated spectra, derived from the same virtual galaxy, successfully replicated observations from galaxies within the mass range of ∼ 10 6 − 10 9.5 M ⊙ and SFRs between 10 −2 to 10 1.58 M ⊙ yr −1 .Hence, the lack of an exact match between the properties of the virtual galaxy and observed counterparts does not present a significant obstacle to comparing the LIS line spectra.We discuss further this aspect in Section 4. In general, the selected vandels galaxies, which cover substantially larger z, M ⋆ , SFR, and magnitudes as compared to classy, is an excellent sample to test the conclusions from our previous study.

COMPARING THE C II AND Si II SPECTRA IN VANDELS AND IN THE SIMULATION
Here, we compare the simulated C II λ1334+C II* λ1335 and Si II λ1260+Si II* λ1265 spectra and the vandels observations.We first explore the overlap in the lines' equivalent widths in Section 3.1, and then identify the mock spectra that best resemble the observed line profiles of each vandels galaxy in Section 3.2.

Overview of C II and Si II measurements
In this section, we measure the EWs of the C II λ1334 and Si II λ1260 spectra to compare the overall consistency of the lines' properties in vandels and in the virtual galaxy.The spectral resolution of the vandels spectra is relatively low, causing the C II λ1334 and C II* λ1335 features (separated by ∼1.2 Å, equivalent to half a resolution element) to blend.For comparison, the separation of the Si II λ1260 and Si II* λ1265 features is approximately 4.5 Å apart, equivalent to roughly two resolution elements.To simplify the comparison of the C II and Si II line spectra, we choose to measure the combined EW(C II λ1334 + C II* λ1335) and EW(Si II λ1260 + Si II* λ1265).
The measurement method is as follows: we normalize the Si II and C II simulated spectra using the median of the flux using a feature-free continuum interval around 1255 Å and 1331 Å, respectively and convolve them to the resolution of the vandels spectra (R∼580).We derive the final EW measurements by integrating the continuum-normalized flux in an interval of ± 2000 km s −1 around λ = 1334.53Å and λ = 1260.42Å, respectively.Spectra dominated by absorption (emission) have positive (negative) EW values.We note that there exists a S II absorption line at 1259.52 Å, blended with Si II λ1260, which is not included in our radiative transfer model.In Gazagnes et al. (2023), we showed that its contribution, resolved in higher-resolution spectra, is typically much weaker than Si II λ1260 and has a negligible impact on the line measurements.
Figure 2 compares the EW(C II λ1334 + C II* λ1335) and EW(Si II λ1260 + Si II* λ1265) values in both the simulation and the vandels observations.The C II EW values scale strongly with the Si II EWs, highlighting that the properties of both ions are closely related.Further, the correlation and scatter in the EW values of the simulation align well with observational data.We observe a slight shift towards higher EWs (above 2.5 Å) in the observations, with 33% of the vandels spectra measurements falling outside the range of simulated values.This difference is likely attributable to the generally larger stellar masses of the vandels galaxies compared to our simulated galaxy, and these objects typically have broad absorption profiles (> 1000 km s −1 in width) unseen in the virtual galaxy (Gazagnes et al. 2023).
The scatter in the simulation's EW values arises because we look at several independent lines of sight, each revealing diverse ISM conditions.This results in a fluctuation of EW values around typical peaks for both ions, as seen in the 1D distribution.Similarly, this line-ofsight variability might also contribute to the scatter observed in real-world data, as each observation captures a unique perspective of the targeted galaxy (see also the discussion in the context of comparing the simulation to the classy observations, Gazagnes et al. 2023).
In the following section, we dive deeper into the comparison of the simulated and observed spectra and analyze the similitude of the line spectral shapes in both datasets.

Finding the best-matching mock spectra
Comparing the EW(C II λ1334 + C II* λ1335) and EW(Si II λ1260 + Si II* λ1265) only offers limited insights into the resemblance of the line profiles.This is because we combine the properties of both resonant and fluorescent features, overlooking the specific shape of the lines.A direct comparison of the shape of UV absorption line profiles is crucial for gaining insights into the geometry of the absorbing gas (see e.g., Scarlata & Panagia 2015;Carr et al. 2018;Carr & Scarlata 2022), in particular for comparing how distinct models reproduce the data (e.g clumpy CGN versus homogeneous wind model, Li et al. 2024).In this section, we directly assess the spectral similitude between the simulated and vandels observations.To achieve this, we look for the best-matching mock spectra from the simulation for each vandels spectrum.We employ the same minimization procedure as outlined in Gazagnes et al. (2023).The mock and observed profiles are normalized using the median flux within featurefree intervals located around 1256 Å and 1331 Å for the Si II and C II spectra, respectively.Then, the 22,500 mock spectra are convolved to the resolution of the vandels observations (R∼580).For each mock spectrum, we compute the reduced χ2 value: , where (2) F obs ion and F sim ion are the observed and simulated flux in the wavelength regions of each ion, respectively.σ obs ion represents the error on the observed flux, and n obs ion is the number of wavelength bins used to compare both spectra.
In our analysis, we address the redshift uncertainty (vandels redshift flag of 3, 4, or, 9, meaning their redshift is reliable at > 95%, see Section 2.2) in our fits by initially conducting a cross-correlation between the observed spectra and each mock spectrum, allowing for the observations to be shifted by at maximum one spectral resolution element.This cross-correlation process determines the optimal shift required to align the spectral features.However, we acknowledge that this approach tends to cancel out any actual differences in the line velocity offsets between the mock and observed spectra.To assess the potential impact of this pre-alignment step, we conducted an independent experiment wherein we identified the best-matching mock spectrum both with and without including the cross-correlation step.While we observed that the typical χ2 of the best match was larger and some best-matching spectra did not align as well when excluding the cross-correlation step, we find that this does not have a significant impact on the overall results and subsequent discussions presented in the following sections.We elaborate on the outcomes of these tests in Appendix A, highlighting that employing an alternative strategy yields similar results.
In Gazagnes et al. (2023), we selected a unique, best-matching mock spectrum for each observed galaxy.Here, for each vandels object, we also look at the distribution of the χ2 for all 22,500 mock spectra to explore the diversity of spectral shape from the simulation that closely mimics a given observation.This approach enables us to further assess the simulation's accuracy in replicating a specific observation by evaluating the number of mock spectra that closely align with it.It also facilitates the examination of common features or characteristics among all mock spectra that well match with the observation, an aspect that we discuss in detail in Section 5.1.
In Figure 3, we present the matching outcomes for three galaxies drawn from our selected vandels sample.We chose representative galaxies covering/spanning a range of stellar masses, specifically 10 9.2 , 10 9.6 , and 10 10.6 M ⊙ .For each galaxy, we show the best-matching mock spectrum (i.e., with the lowest χ2 value).The best-matching spectrum replicates reasonably well the line profiles in the two lower-mass galaxies ( χ2 < 1.5).However, they fall short of capturing the complete width of the broad absorption lines observed in the 10 10.6 M ⊙ galaxy ( χ2 > 3).This result reaffirms the findings from Gazagnes et al. (2023); Blaizot et al. (2023), emphasizing that mock line spectra generated from a 10 9.2 M ⊙ virtual galaxy do not provide an accurate match for the most massive galaxies, where the ISM gas environment may be significantly different (e.g., more turbulent or larger outflows).Figure 4 shows the overall distribution of χ2 for the whole sample.The top panel presents the distribution of the χ2 value of the best-matching spectrum for each object.37% of the spectra are well reproduced ( χ2 < 1), and 83% reasonably matched ( χ2 < 2).These numbers are slightly lower than for the classy sample of z ∼ 0 galaxies (60% with χ2 < 1, 89% with χ2 < 2).These results may be explained by the presence of more massive galaxies in the selected vandels sample whose spectra are more complex to match using the current simulation.We discuss further the influence of galaxy properties such as M ⋆ and SFR on the agreement between the best matches and observations in the next section.
4. CONNECTING THE BEST-MATCHED SPECTRA WITH GALAXY PROPERTIES Section 3.2 showed that mock spectra, derived from a single virtual galaxy, reasonably well replicate ( χ2 < 2) 83% of the selected vandels galaxies exhibiting a diverse range of line profiles and characteristics.These results echo those previously discussed in the context of the comparison with the low-z classy sample in Gazagnes et al. (2023).In that study, we emphasized that monitoring the virtual galaxy's evolution over 690 million years and exploring it from various viewing angles enabled us to encounter a wide range of ISM configurations characterized by evolving properties.In this section, we further explore the relationship between the best matches and galaxy properties, examining which factors determine the agreement between simulated and observed line spectra.Section 4.1 explores the influence of M ⋆ and SFR for finding the best match to a given observation, and Section 4.2 delves deeper into the role of star formation bursts into producing LIS spectral features that resemble the observations.4.1.The influence of M ⋆ and SFR when matching simulations and observations In this section, we discuss the influence of M ⋆ and SFR on the matching quality.Due to the lack of precise overlap between M ⋆ and SFR between the simulation and vandels (see Fig 1), we investigate whether a greater number of mock spectra closely match a given observation when the corresponding vandels galaxy has a stellar mass and SFR more similar to those of the virtual galaxy.To do so, we analyze the distribution of the average number of mock spectra with χ2 values less than 1, 2, and 3, respectively, per observed galaxy, across bins of increasing mass and SFR.
The top panel of Figure 5 shows the distribution of the average number of mock spectra with χ2 values less than 1, 2, and 3 per galaxy across increasing mass bins.The first bin (M ⋆ < 10 9.3 M ⊙ ) directly aligns with the stellar mass range of the virtual galaxy (∼ 10 8.9 to 10 9.3 M ⊙ ).Interestingly, this bin exhibits the highest average number of mock spectra with χ2 values less than 1, 2, and ) > 1 0 1 .9 10 0 10 1 10 2

Counts
Figure 5. Top: histogram presenting the average number of mock spectra with χ2 < 1 (black), χ2 < 2 (gray), and χ2 < 3 (white) per observed galaxy, in bins of stellar mass.The most massive vandels galaxies (∼> 10 10 ) are not well replicated by the mock spectra for our virtual galaxy.Bottom: Same but for the SFR.
3, indicating a greater resemblance between these mock spectra from the simulation and the observations within the same mass range.Specifically, for galaxies with M ⋆ < 10 9.8 M ⊙ , there are, on average, 5 to 8 mock spectra with χ2 < 1, while none are found for galaxies in the two largest stellar mass bins (∼> 10 10 M ⊙ ).These results echo Gazagnes et al. (2023); Blaizot et al. (2023), indicating that the current simulation may not accurately represent galaxies with M ⋆ > 10 10 M ⊙ .Furthermore, we observe a declining trend in the average number of mock spectra with χ2 < 2 and χ2 < 3 as a function of stellar mass.This trend suggests that as galaxies increase in mass, it becomes increasingly challenging to find multiple mock spectra closely resembling the observed data.
The bottom panel of Figure 5 focuses on the SFR.We note that in vandels, no galaxy exhibits an SFR consistent with the SFR of our virtual galaxy (∼ 10 0 to 10 0.6 M ⊙ yr −1 ).This panel shows that, similarly to the analysis for stellar mass, the bins representing the most extreme star-forming galaxies yield no matches with χ2 < 1.Given that these galaxies are also the most massive, this outcome aligns with the observations from the top panel.We do not observe a declining trend in the average number of mock spectra with χ2 < 1 from the lowest SFR bin to the largest SFR bin.However, we do discern a slight decrease in the average number of χ2 < 2 and χ2 < 3, suggesting that, similarly to M ⋆ , albeit to a lesser extent, it becomes increasingly challenging to find a spectrum that adequately reproduces an observation as the galaxy exhibits a higher SFR.
In general, our results reveal the existence of certain dependencies on both M ⋆ and SFR.Specifically, the spectra of the most massive or star-forming galaxies in vandels are not well matched in the simulation.Conversely, a greater number of mock spectra resemble observations with properties similar to those of the virtual galaxies.However, we also find that an exact match between observation and simulation properties is not strictly necessary to find mock spectra resembling real galaxies.This outcome likely stems from the extensive sampling of 300 viewing angles over a 690 Myr period, which unveils rarer ISM and CGM configurations that may be more prevalent in galaxies with distinct properties.
Previous studies revealed noticeable yet scattered dependencies of LIS absorption line properties on M⋆ and SFR (e.g., the EW and velocity offsets; Rupke et al. 2005;Chisholm et al. 2015bChisholm et al. , 2016)).Our findings may provide important clues to understanding the scatter observed in typical line property to stellar mass trends in those studies.Given that galaxies are not idealized symmetric objects, there exists a distribution of ISM and CGM configurations observable depending on the viewing angle.The peak and deviation of this distribution depend on galaxy properties, but our measurements only capture a unique point in this distribution contingent upon our line of sight which may deviate from this peak.
Overall, Fig 5 suggests that the fits we obtain are genuine rather than coincidental, i.e. that the scatter in the mock spectra reflects physical variations in time and line-of-sight that are bounded in a range loosely set by the mass of the system.Achieving good agreement between observed and simulated spectra doesn't necessarily imply an exact match between the virtual galaxy environment and the object considered, but it suggests that the observed galaxy falls within a certain dynamical range that can be replicated with a given simulation.In turn, analyzing the most massive objects requires the inclusion of virtual galaxies with properties closer to the observed range (Rosdahl et al. 2018, e.g. using those from the sphinx simulations).While this is an important effort for delineating the range of spectral properties that single-galaxy simulations can accurately reproduce as a function of their properties, this effort extends beyond the scope of the present paper. .The evolution of the (black) and AB magnitude at 1250 Å (red) of the virtual galaxy over the 75 time steps used to produce the mock C II and Si II line profiles (690 Myrs period).The magnitude curve is derived using the median and standard deviation from the 300 line-of-sight measurements at each time step.At the top, we present a histogram illustrating the distribution of time steps where the 2,790 best-matching mock spectra are found (taking 30 best mocks for each of the 93 vandels galaxies selected in Section 3.2).Overall, the best-matching spectra dominantly originate from luminosity peaks closely preceding the two most intense star-formation bursts.
An important and novel aspect of our analysis involves examining the average distribution of χ2 values across multiple mock spectra for each galaxy.This strategy allows us to establish a more robust understanding of the range of mock spectra that closely match the observed spectrum, as opposed to relying solely on a single best-matching spectrum.In the subsequent sections, we build upon this distribution to interpret our results, focusing on a smaller subset of vandels galaxies that have been reasonably well replicated by the simulation.Specifically, we consider the 93 vandels galaxies with at least 30 mock spectra with χ2 < 2 as suitable candidates for further investigating the origin of the good agreement between simulation and observation.This selection criterion is derived from the insights presented in the lower panel of Figure 4, which demonstrates that galaxies within the two highest stellar mass bins fail to meet this threshold.
In this section, we explore further the connection between the virtual galaxy star-formation history and the origin of the good agreement between observations and simulations.To do so, we look at the distribution of time steps at which each of the best mock spectra was generated and investigate whether they originate from specific stages in the galaxy's evolutionary timeline.To do so, we extract the corresponding 2,790 best-matching mock spectra (30 mock spectra × 93 well-matched vandels galaxies).
In Figure 6, we plot the evolution of the SFR of the virtual galaxy over the period considered in this work.We measure the AB magnitude of the virtual galaxy taking the flux at 1250 Å (M 1250 Å) for each output and direction and plot the evolution of the median and standard deviation values from 300 directions.The distribution of the best-matching spectra is very tightly correlated with the evolution of M 1250 Å, and most of the best mocks originate from bright luminosity peaks.Interestingly, these peaks typically precede strong starformation bursts where the galaxy is at its peak of star formation.However, we note that the SFR of the virtual galaxy at time t is calculated as an average over the preceding 100 Myrs.As a result, the luminosity peaks may coincide with SFR peaks, albeit with a delay in its influence on the 100-Myrs averaged curve.
We find that 85 out of 93 selected vandels galaxies exhibit at least one best-matching mock spectrum originating from either the 1.53 Gyrs or the 2.02 Gyrs time steps, corresponding to the two highest peaks in the histogram distribution near the peaks of SFR.This finding has two implications.Firstly, it indicates that a comparable success rate in matching spectra to observations could have been achieved by solely considering these two time steps.Secondly, given the diversity of line properties exhibited by the 93 vandels galaxies (shown in Fig 2 ), it implies that the LIS profiles generated at these times encompass a wide range of variability such that two distinct viewing angles likely have significantly different line profiles.We note again that, while a strong agreement between observed and simulated spectra during these two bursts is notable, it does not systematically imply an exact match between the environments of the observed and simulated objects.
The fact that the best-matching mock spectra aren't randomly distributed in time but rather linked to specific star-formation phases is a fundamental result as it suggests that these phases have a direct impact on the LIS metal line spectra.We highlighted in Gazagnes et al. (2023) that one of the main differences in the line properties was that the absorption lines are significantly more redshifted during periods of relative quiescence as compared to star-formation events (see Figure 10 in particular).This suggests that the gas kinematics may be one powerful tracer of the current star-formation activity (see also Trebitsch et al. 2017;Rosdahl et al. 2018).
As detailed in Section 4.1, although the virtual galaxy properties do not precisely match those of the selected vandels sample, we still identify mock spectra closely resembling the observations.This means that, in the context of understanding how LIS metal line spectra connect to the dominant physical processes, the relative evolution of the SFR of a galaxy is likely more relevant than its absolute fluctuations at a certain time.One can also note the four-fold variation in brightness, closely tied to the SFR, of the virtual galaxy.Hence, at high z, where galaxies may exhibit bursty behavior due to ample neutral hydrogen reservoirs, the observation of bright luminous objects may be systematically biased toward objects experiencing intense star-formation bursts.However, these objects may not be as massive as their luminosity would imply since they are caught in a very specific phase of their evolution.
To conclude on a rather general note, Figure 6 has critical implications for understanding existing observational biases inherent in current high-z surveys.vandels is a UV continuum selected sample, meaning that it is implicitly biased towards UV-bright galaxies (Garilli et al. 2021) likely in a phase of star-formation (a bias that we further accentuated by selecting the highest SNR objects).Interestingly, when matching mock and observed spectra, we recover the selection function of the observed survey because the best mock spectra correspond to times when the virtual galaxy is also UV-bright.This means that the spectral properties of such galaxies are significantly influenced by their current stage in the evolutionary cycle, differing quite notably from those in less luminous or less active star-forming phases.

PREDICTING THE LEAKAGE OF IONIZING PHOTONS IN vandels
In the high-z universe, the opacity of the IGM prevents any robust measurements of the LyC escape fraction from direct observations of the flux below 912 Å.Nevertheless, the study of low−z LCEs has enabled us to determine indirect diagnostics to f LyC esc (e.g., O 32 , Lyα Mg II, β slopes, or LIS covering fractions, see Jaskot & Oey 2013;Verhamme et al. 2017;Izotov et al. 2023;Chisholm et al. 2020Chisholm et al. , 2022;;Flury et al. 2022a;Gazagnes et al. 2018;Reddy et al. 2016;Saldana-Lopez et al. 2023).Although empirical trends between f LyC esc and these diagnostics show a certain scatter, they are already used to put the first constraints on the escape fractions of reionization-era star-forming galaxies (e.g Choustikov et al. 2023;Witten et al. 2023).Here, for each vandels galaxy in our sample, we predict f LyC esc using the simulation and compare these predictions with previous constraints from Begley et al. (2022), Begley et al. (2024), andSaldana-Lopez et al. (2023).Section 5.1 details our simulation-based approach to predict f LyC esc for each vandels galaxy, Section 5.2 explores the consistency of these predictions with other known LyC diagnostics, and Section 5.3 investigates what makes LIS metal lines good tracers of the LyC leakage.The RHD simulation used in this work includes a selfconsistent propagation of the ionizing photons through radiative transfer, meaning that we can follow for each line of sight the amount of ionizing photons escaping the galaxy.To determine f LyC esc in the simulation, we take the approach of (Mauerhofer et al. 2021): for each of the 300 directions and 75 outputs, we generate mock spectra of the virtual galaxy between 10 Å and 912 Å, such that each ionizing photon produced encounter an optical depth that is a sum of contributions from H, He, He+, and dust.We then compute the total number of ionizing photons going out of the virial radius and di-vide it by the intrinsic number produced to get f LyC esc .In this work, we only focus on the escape fraction at 900 Å, determined by averaging the spectra between 890 Å and 910 Å, to maintain consistency with LyC constraints in the existing literature.We note that for the current virtual galaxy, in the absence of hard ionizing sources, the escape fraction at 900 Å closely relates to the total escape fraction of ionizing photons (see Mauerhofer et al. 2021, for more details on this aspect).
Using the detailed strategy, we determine the LyC escape fraction corresponding to each of the 22,500 mock spectra generated.Then, to extract a f LyC esc prediction for each of the 93 vandels best-matched objects, we take the median and standard deviation of the 30 f LyC esc values corresponding to the 30 best-matching mock spectra.To differentiate our predictions from direct LyC measurements, we refer to these predicted f LyC esc as f LyC esc (pred) .Figure 7 presents the distribution of the 93 f LyC esc (pred) measurements.The vast majority of values fall below 1% escape, and only two objects have f LyC esc (pred) > 10%.The average f LyC esc (pred) for the 93 galaxies is 0.01 ± 0.02.While lower, this estimate falls within 2σ of the 0.07±0.02escape fraction reported by Begley et al. (2022) who directly constrained f LyC esc using ultra-deep, ground-based, U-band imaging.Interestingly, our average predicted escape fraction is highly consistent with Saldana-Lopez et al. ( 2023) who determined an average f LyC esc of 0.02±0.01using an indirect approach combining the residual flux of LIS absorption lines (including Si II λ1260, 1302, 1527, O I λ1302, and C II λ1334) and the expected dust attenuation at 912 Å (Eq.( 5) from Chisholm et al. 2018).The consistency of both approaches is interesting, particularly because the predictions derived in this work do not require an explicit forecast of the dust attenuation at 912 Å.We further discuss this point in Section 5.3.
Using RHD simulations to predict the escape fraction from the simulation is a promising avenue for determining f LyC esc of actual observations (see also Katz et al. 2020Katz et al. , 2022Katz et al. , 2023;;Choustikov et al. 2023).In this study, the primary limitation stems from relying on a single zoomin simulation to forecast f LyC esc .This approach inherently restricts the range of potential f LyC esc values to those observed within the simulation itself, with the maximum f LyC esc observed in our virtual galaxy being around 42% (Mauerhofer et al. 2021).Further work is necessary to corroborate the simulation's predictions by comparing them with direct f LyC esc measurements from lower-redshift samples.However, we demonstrate in the following subsection that the simulation-based predictions are in alignment with the anticipated intensity of LyC leakage for each object, as inferred from their Lyα and UV spectral properties.

Testing the consistency of the predicted f LyC esc
Here we assess if the f LyC esc (pred) derived in the previous section correlates with other established LyC diagnostics.

Comparing to the Lyα escape fractions
The Lyα line is a powerful tracer for the LyC leakage (e.g., Verhamme et al. 2017;Steidel et al. 2018;Izotov et al. 2023), as Lyα photons are affected by the same neutral hydrogen gas and dust extinction regulating the escape of LyC photons.Begley et al. (2024) investigated the connection between the escape fraction of Lyα (f Lyα esc ) and f LyC esc within a sample of 152 vandels objects, 11 of which are also part of the 93 well-matched vandels galaxies.These 11 objects span M ⋆ values from 10 8.8 to 10 9.8 M ⊙ and SFR values between 10 0.8 to 10 1.8 M ⊙ yr −1 .Further, the median χ 2 of their 30 best mock spectra spans values from 0.7 to 1.9.This indicates that these 11 objects represent the broader wellmatched vandels sample examined in this study and do not exhibit bias toward specific groups of objects that are "exceptionally" well-matched or possess a particular set of galaxy properties.
Figure 8 compares the predicted LyC escape fractions and the f Lyα esc values measured in Begley et al. ( 2024) for these objects.We observe a significant correlation between f LyC esc (pred) and the measured f Lyα esc .Notably, this trend is highly consistent with two empirical f LyC escf Lyα esc relationships, one derived in Begley et al. (2024) using vandels galaxies, and the other established in Izotov et al. (2023) using a sample of z < 0.5 LyC emitters.This consistency is remarkable as it shows that our LyC escape fraction predictions, derived from comparing LIS metal line spectra in observations and RHD simulations, adhere to established Lyα-LyC observational trends.This result underscores the consistency of both Lyα and LIS metal lines as indicators of f LyC esc .Further, it brings confidence in the f LyC esc (pred) , suggesting these predictions may be close to the actual LyC escape fractions of these objects.

5.2.2.
Comparing UV spectral differences Begley et al. (2024) investigated the spectral difference of composite spectra created using bins of f Lyα esc in vandels.Here, we explore a similar exercise but producing composite stacks in three bins of f LyC esc (pred) : below 1%, between 1% and 10%, and larger or equal to 10% respectively.We use these three f LyC esc (pred) bins to represent categories of non, weak, or relatively strong leakers, respectively.As discussed in Section 5.1, only two objects have f LyC esc ≥ 10%.To increase the number of objects in those bins, we also include four galaxies whose predicted escape fractions are consistent with a strong leakage (i.e.f LyC esc (pred) +1σ > 10%).Figure 9 present these three stacks, with four panels zooming in on the Lyα profiles, LIS absorption lines (from Si II, O I, and C II), and the high-energy UV features C IV and He II, tracing hard ionizing radiation field above 47.9 and 54.4 eV respectively.For each of the stacks, we measure the β slope between 1250 and 1700 Å (selecting feature-free intervals in that range) by fitting a line such that F obs ∝ λ β .β is a robust indicator of the LyC leakage in low-z LyC leakers (Chisholm et al. 2022) because the UV slope strongly depends on the dust extinction which also regulates the escape of ionizing photons.β is especially powerful because its measurement can be relatively simply extracted from spectra of high-z objects with JWST (e.g., Cullen et al. 2023a,b).
Examining the distinctions among the three composite spectra, we observe several trends.Stacks with higher predicted f LyC esc exhibit bluer β values, an increase in Lyα flux, and weaker LIS absorption lines, respectively.These observations align with the physical picture that lower levels of dust and absorbing gas facilitate and enhance the escape of both LyC and Lyα photons (Gazagnes et al. 2020;Begley et al. 2024).We note that the measured β slopes in the weak and relatively strong leakers are "redder" than what is typically expected from low-z LyC emitters.Using the β to f LyC esc relation from Chisholm et al. (2022) derived from the LzLCS sample (Flury et al. 2022b), a β slope of −1.35, and −1.84, would correspond to f LyC esc of ∼ 1%, and ∼ 3%, respectively, hence lower than what is expected based on the simulation-based predictions.
Regarding the high-ionization UV features, we observe no difference in the C IV absorption or He II emission.However, we do observe a slightly larger C IV P-Cygni emission flux in the f LyC esc > 10% stack, which could be either due to a younger stellar population (Chisholm et al. 2020) or enhanced collisionally-excited C IV emission in a high-ionization galaxy environment (Berg et al. 2019).Prominent C IV and He II emission are sometimes observed in extreme LyC leakers (f LyC esc > 20% or f Lyα esc > 50%, see Schaerer et al. 2022;Naidu et al. 2022) with large intrinsic production of ionizing photons.None of the f LyC esc (pred) derived in this work are consistent with the presence of an extreme leakage2 .We note that Marques-Chaves et al. (2022) explored the correlation between strong LyC leakage and spectral hardness in a sample of 89 galaxies and found that extreme leakers do not have harder ionizing spectra than non-leakers at similar metallicity.Hence, this aspect may also explain the absence of differences seen for higher-energy lines in the three stacks.
Overall, this section highlights that the predicted LyC escape fractions derived from this simulation align closely with empirical LyC trends observed in observational studies, reinforcing our confidence in the physical accuracy of the current RHD simulation.This means we can use the virtual galaxy environment to elucidate the physical processes responsible for the LyC escape (see also Katz et al. 2020Katz et al. , 2022Katz et al. , 2023;;Choustikov et al. 2023), and we explore that aspect in the next section.

What makes LIS lines good tracers of f LyC esc
In the previous Section, we highlighted that the f LyC esc (pred) values, extracted from the simulation from the simple comparison of the LIS line spectra in the observations and simulation, consistently scale with independent LyC diagnostics such as Lyα and the β slopes, and align well with direct and indirect LyC constraints established by previous studies (Begley et al. 2022;Saldana-Lopez et al. 2023).This result is significant as it reinforces both theoretical and empirical evidence suggesting that LIS metal line spectra are reliable indicators of LyC leakage (e.g.Heckman et al. 2001;Alexandroff et al. 2015;Gazagnes et al. 2018;Saldana-Lopez et al. 2023).We use three bins of non, weak, and (relatively) strong leakers using cuts at 1 and 10% leakage.We measure for each stack the β slope using the observed continuum flux between 1250 and 1700 Å, and the four panels on the bottom zoom on the Lyα profiles, LIS absorption lines, and high-energy C IV and He II features.
Yet, this outcome brings us to the following question: what makes the LIS metal lines good tracers of the LyC leakage?The tight f LyC esc empirical trends found in low-z galaxies over the past suggested that LyC photons escape through relatively dust-free, low neutral gas column-density channels (Gazagnes et al. 2018;Chisholm et al. 2018;Saldana-Lopez et al. 2023;Reddy et al. 2016).Further, Flury et al. (2022a) found a significant correlation between larger LyC leakage and the objects' compactness, which may support that stellar feedback more efficiently clear escape channels in such environment (Cen 2020).Hence, LIS metal lines may trace the LyC leakage because they indirectly probe the presence of compact, dust-free, and neutral-gas-free regions through which LyC photons can escape.
We explore here the complex relation between f LyC esc , LIS metal line EWs, dust attenuation, and compactness.To do so, we generate Integral Field Unit cubes for each of the 300 lines of sight for a single simulation output of the galaxy at t = 2 Gyrs (the time-step where most of the best mock spectra are found).Each data cube comprises 64 × 64 spatial pixels, with each pixel corresponding to a physical size of 120 × 120 pc, covering a total area of approximately 7.7 by 7.7 kpc on the entire cube face.
For each of these 300 viewing directions, we extract three essential properties: (i) the radius at which the integrated flux accounts for more than 75% of the total flux (referred to as r 75 ), (ii) the luminosity-weighted dust attenuation, calculated by comparing the observed Figure 10.f LyC esc as a function of the average EW of the Si II and C II absorption lines (noted EW(LIS)) along 300 lines of sight for a single simulation output of the virtual galaxy at t = 2 Gyrs.Squares and circles differentiate sightlines where less or more than 75% of the intrinsic flux escapes the galaxy, respectively (using the luminosity weighted F obs /Fint ratio).Each point is color-coded by r75, which is the radius at which the integrated flux accounts for more than 75% of the total escaping flux in that direction.
and intrinsic flux and taking the luminosity-weighted average across all pixels, and (iii) the EW of the absorption lines for Si II and C II (denoted as EW(LIS)) along that specific line of sight.We opted for the use of r 75 as a measure of compactness instead of the halflight radius r 50 primarily because in numerous viewing directions, a single pixel accounts for 50% of the total flux.This makes it challenging to discern meaningful variations when relying only on r 50 .
Figure 10 shows the connection between f LyC esc in the simulation, EW(LIS), r75, and the dust attenuation.Directions with high f LyC esc exhibit low EW(LIS), minimal dust attenuation, and a smaller r75.This aligns with the expected LyC-leakage picture: in our virtual galaxy, LyC photons likely escape through areas that are spatially confined (indicated by a low r 75 ), have a low neutral gas amount (low EWs of ions tracing the low-ionized gas) and are dust depleted (higher Fobs/F int ratio).
Interestingly, these three characteristics (lack of neutral gas and dust and spatial compactness) tend to occur together.Directions with low EW(LIS) usually have lower dust attenuation and a smaller r 75 .Therefore, this explains why LIS line profiles are good tracers of the f LyC esc values: they indirectly reflect all the factors influencing LyC leakage.We note that this is not always true, a few directions have only one or two of these criteria met (e.g., low r 75 and dust attenuation but high EW(LIS), which results in lower f LyC esc .These cases might explain the observed scatter in LyC diagnostic trends since LyC diagnostics usually trace only one of these characteristics.
In general, Figure 10 supports a physical picture where LyC photons escape through compact regions, depleted of neutral gas and dust.However, further work is needed to refine this model and to more thoroughly evaluate the reliability and the physics behind existing LyC diagnostics.This effort will be the topic of a separate paper.

SUMMARY
In this study, we compared mock UV C II and Si II absorption and emission line spectra simulated from a single ∼ 10 9 M ⊙ virtual galaxy to the observed spectra of 131 z ∼ 3 galaxies of stellar mass between ∼ 10 9 and 10 10.5 from the vandels survey (McLure et al. 2018;Pentericci et al. 2018;Garilli et al. 2021).We assessed the degree of resemblance between the simulated and observed spectra and investigated the potential of the LIS metal line spectra as tracers of galaxy properties and Lyman Continuum leakage within star-forming galaxies.In particular, we established a simulation-based f LyC esc estimation method which enabled us to predict the strength of the LyC leakage in the selected vandels sample.We summarize our findings as follows: • The correlation and dispersion in the distribution of equivalent widths for the C II and Si II line spectra closely match the vandels observations.We note a slight shift towards higher EWs in the observations, with 33% of the vandels measurements falling outside the range of simulated values.This difference is attributable to the generally larger stellar masses of the selected vandels galaxies compared to the virtual galaxy (Section 3.1).
• For each vandels galaxy, we determine the similarity ( χ2 ) between the mock spectra and the observed absorption and emission line profiles.We find for 37% of the vandels spectra a most spectrum with excellent agreement ( χ2 < 1), while 83% are reasonably well matched (best mock spectrum with χ2 < 2).Notably, we observed a dependence on stellar mass, with the simulated spectra failing to reproduce the most massive galaxies (M⋆ > 10 10.1 M⊙, Section 4.1) due to the presence of broader absorption profiles.This finding further supports the conclusions from Gazagnes et al. (2023); Blaizot et al. (2023), highlighting the limitations of the current virtual galaxy in representing the environments of the most massive galaxies.
• The best-matching mock spectra originate from time steps when the virtual galaxy's luminosity and SFR are at or close to a peak (Section 4).These luminosity peaks are four times larger than during weaker star formation phases.The influence of intense star formation in explaining the good agreement between simulated and observed spectra is consistent with the nature of the selected vandels galaxies as UV-bright, relatively star-forming objects.This aspect emphasizes that star-formation phases have a direct impact on the LIS metal line spectra.In general, this outcome reminds us of the strong observational biases of current surveys that disproportionately favor galaxies in a star-bursting phase at any given stellar mass (Figure 6).
• We predicted the f LyC esc for each galaxy in the selected vandels sample using the virtual galaxy environment (f LyC esc (pred) , Section 5.1).We derived an average predicted f LyC esc (pred) of 0.01 ± 0.02, within 2σ of the value 0.07±0.02reported by Begley et al. (2022), and directly consistent with the indirect constraint from Saldana-Lopez et al. ( 2023) who used the residual flux of LIS absorption lines and the expected dust attenuation at 912 Å.
• The simulation-based f LyC esc (pred) tightly correlate with the Lyα escape fractions from Begley et al. (2024) and are in remarkable agreement with the empirical f LyC esc -f Lyα esc relations from Izotov et al. (2023) and Begley et al. (2024) (Section 5.2.1, Figure 8).We produced composite spectra in bins of f LyC esc (pred) (Section 5.2.2) and show that the stack comprising the higher f LyC esc (pred) galaxies exhibit bluer β values, an In Section 3.2, we implemented a matching strategy to find, for each vandels spectrum analyzed in this work, the best-matching mock spectra in the simulation.We discussed that, although the selected vandels objects have a redshift reliable at > 95%, a few of these redshift estimates rely on ISM features that can be kinematicsdependent and therefore result in slight offsets (typically within a few hundred km s −1 at most) compared to the systemic redshift typically estimated using optical nebular emission lines.In our main analysis, we addressed this redshift uncertainty by initially conducting a cross-correlation between the observed spectra and each mock spectrum to determine the optimal spectral shift required to align the spectral features.Yet we noted that this approach nullifies actual differences in the line velocity offsets between the mock and observed spectra.
To assess the potential impact of this pre-alignment step, we conducted an independent experiment, detailed here, where we identified the best-matching mock spectra both with and without including this "redshift uncertainty".In Figure 11, we present the main differences of using either approach.The top panel shows the distribution of χ2 values and highlights that a few more galaxies have χ2 < 1.5 when we first pre-align the mock and observed spectra.This result is not unexpected as by performing a pre-alignment, we already find the spectral shift that provides the minimal difference between both spectra.The middle panel reproduces Figure 1 from Section 4.2 and compares whether we observe any differences in the time-steps origin of the 30 best-matching spectra for each vandels object when we do or do not account for potential redshift uncertainties.The two distributions are highly similar, suggesting that the conclusions drawn in Section 4.2 are unaffected by the choice of either method.Finally, the bottom panel compares the differences in the predicted f LyC esc , using the methodology detailed in Section 5, when one or the other approach.The predicted values are relatively consistent with each other and choosing one or the other approach does not impact the overall predicted average LyC escape fraction (0.01±0.02 in both cases).
Overall, Figure 11 emphasizes that choosing a bestmatching approach that does not include a prealignment step would have had no impact on the main conclusions drawn in this work.The comparison of the best-match χ2 for the two different cases.As expected, allowing for a spectral shift when finding the best matching spectra leads to typically better χ2 .Middle: The differences in the distribution of the best-matching spectra with and without z uncertainty.We do not find a significant difference in both distributions.Right: Comparing the predicted f LyC esc when using the two different approaches.We do observe small variations of the f LyC esc (pred) , however, these variations are within the uncertainties and do not change the interpretation of galaxies which could be categorized as non, weak, or strong leakers.Overall, these there panels support that choosing one or the other approach has no impact on the main conclusions drawn in this work.
Figure1.Histograms comparing the range of redshift, stellar mass, SFR, and absolute AB magnitude at 1250 Å for the 131 vandels galaxies selected in this work (z > 2.8 and S/N > 5, in blue) and for the virtual galaxy (yellow).For the latter, the ranges shown are determined by the minimum and maximum values for each property over the 690 Myrs period used to generate the mock spectra.

Figure 2 .
Figure 2. The comparison of the EW(C II λ1334 + C II* λ1335) and EW(Si II λ1260 + Si II* λ1265) in the simulation (black contour) and in the vandels observations (circles).The vandels points are color-coded by their stellar mass (the virtual galaxy has M⋆ within ∼ 10 9 − 10 10.5 M⊙), and we show in the bottom right the typical uncertainty on the measurements.The black contour includes the entire range of the mock spectra measurements.The vandels galaxy sample exhibits a marginally wider range of EW values, featuring notably larger values in comparison to the simulation.

Figure 3 .
Figure3.The best-matching spectra of the Si II λ1260+Si II* λ1265 (top panels) and C II λ1334+C II* λ1335 (bottom panels) line profiles for three vandels galaxies of different M⋆ (shown the top left corner of the top panels) using the procedure detailed in Section 3.2.The black and orange spectra are the observed flux and best-matching spectrum (lowest χ2 values, indicated in the top left corner of the bottom panels), respectively.The vertical dotted lines indicate the expected locations of each transition.The higher mass galaxy, which presents broad absorption profiles, is not well matched by the lowest χ2 spectrum.

4. 2 .
Figure6.The evolution of the (black) and AB magnitude at 1250 Å (red) of the virtual galaxy over the 75 time steps used to produce the mock C II and Si II line profiles (690 Myrs period).The magnitude curve is derived using the median and standard deviation from the 300 line-of-sight measurements at each time step.At the top, we present a histogram illustrating the distribution of time steps where the 2,790 best-matching mock spectra are found (taking 30 best mocks for each of the 93 vandels galaxies selected in Section 3.2).Overall, the best-matching spectra dominantly originate from luminosity peaks closely preceding the two most intense star-formation bursts.

5. 1 .Figure 7 .
Figure 7.The distribution of f LyCesc (pred)  for the vandels sample analyzed in this work.The f LyC esc (pred) is determined by taking the average f LyC esc value across the 30 best mock spectra for each of the 93 selected vandels objects.The average f LyC esc (pred) for this sample is 0.01 ± 0.02.

Figure 8 .
Figure 8.The comparison of the f LyC esc (pred) derived using the simulation and the Lyα escape fractions measured in Begley et al. (2024) for the 11 vandels objects that are present in both samples.The red and blue dotted lines show the empirical relations found in Izotov et al. (2023) (using a z < 0.5 LyC emitter sample) and in Begley et al. (2024) (using the vandels sample), respectively.

Figure 9 .
Figure9.Composite stacks of vandels spectra based on the f LyC esc (pred) from this work.We use three bins of non, weak, and (relatively) strong leakers using cuts at 1 and 10% leakage.We measure for each stack the β slope using the observed continuum flux between 1250 and 1700 Å, and the four panels on the bottom zoom on the Lyα profiles, LIS absorption lines, and high-energy C IV and He II features.

Figure 11 .
Figure11.Exploring the impact of accounting for potential redshift uncertainties when finding the best-matching mock spectra.Top: The comparison of the best-match χ2 for the two different cases.As expected, allowing for a spectral shift when finding the best matching spectra leads to typically better χ2 .Middle: The differences in the distribution of the best-matching spectra with and without z uncertainty.We do not find a significant difference in both distributions.Right: Comparing the predicted f LyC