# 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

Simon Gazagnes, Fergus Cullen, Valentin Mauerhofer, Ryan Begley, Danielle Berg, Jeremy Blaizot, John Chisholm, Thibault Garel, Floriane Leclercq, Ross J. Mclure, et al.

► To cite this version:

Simon Gazagnes, Fergus Cullen, Valentin Mauerhofer, Ryan Begley, Danielle Berg, et al. 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. The Astrophysical Journal, 2024, 969, 10.3847/1538-4357/ad47a4. insu-04838862

HAL Id: insu-04838862

https://insu.hal.science/insu-04838862v1

Submitted on 15 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

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

CrossMark

# 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

Simon Gazagnes¹, Fergus Cullen², Valentin Mauerhofer³, Ryan Begley², Danielle Berg¹, Jeremy Blaizot⁴, John Chisholm¹, Thibault Garel⁵, Floriane Leclercq¹, Ross J. McLure², and Anne Verhamme⁵

¹Department of Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712-1205, USA; gazagnes@utexas.edu

²Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK

³Kapteyn Astronomical Institute, University of Groningen, P.O Box 800, 9700 AV Groningen, The Netherlands

⁴Centre de Recherche Astrophysique de Lyon, UMR5574, F-69230, Saint-Genis-Laval, France

⁵Department of Astronomy, University of Geneva, 51 Chemin Pegasi, 1290 Versoix, Switzerland

Received 2024 January 30; revised 2024 April 28; accepted 2024 May 2; published 2024 July 1

## Abstract

We compare mock ultraviolet C II and Si II absorption and emission line features generated using a ~$10^9 M_\odot$ virtual galaxy with observations of 131 $z \sim 3$ galaxies from the VANDELS survey. We find that the mock spectra reproduce reasonably well a large majority (83%) of the VANDELS spectra ($\chi^2 < 2$), but do not resemble the most massive objects ($\gtrsim 10^{10} M_\odot$), 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 ($f_{\text{esc (pred)}}^{\text{LyC}}$) using the environment of the virtual galaxy. We derive an average $f_{\text{esc (pred)}}^{\text{LyC}}$ of 0.01 ± 0.02, consistent with other estimates from the literature. The $f_{\text{esc (pred)}}^{\text{LyC}}$ are tightly correlated with the Ly$\alpha$ escape fractions and highly consistent with observed empirical trends. Additionally, galaxies with larger $f_{\text{esc (pred)}}^{\text{LyC}}$ exhibit bluer $\beta$ slopes, more Ly$\alpha$ flux, and weaker low-ionization absorption lines. Building upon the good agreement between $f_{\text{esc (pred)}}^{\text{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.

*Unified Astronomy Thesaurus concepts: Ultraviolet astronomy (1736); Interstellar medium (847); Starburst galaxies (1570); Galaxy evolution (594)*

## 1. Introduction

Ultraviolet (UV) absorption and emission line features are formidable tracers of the interstellar and circumgalactic media (ISM and CGM) properties within star-forming 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; Shapley et al. 2003; Rubin et al. 2011; Erb et al. 2012; Wofford et al. 2013; James et al. 2014; Chisholm et al. 2015; Heckman et al. 2015; Chisholm 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 (LyC) 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; Chisholm et al. 2018; Gazagnes 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_{\text{esc}}^{\text{LyC}}$, where LyC represents the Lyman continuum photons with $\lambda < 912 \AA$) 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,

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
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., Rivera-Thorsen et al. 2017;
Jaskot et al. 2019). To capture the vital features that are directly
related to certain properties such as the LyC escape, it is
imperative that we 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 consis-
tently 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
the work of Mauerhofer et al. (2021), who used a zoom-in
high-resolution radiation hydrodynamics (RHD) simulation
postprocessed with radiative-transfer techniques, Blaizot et al.
(2023) highlighted that the mock 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 environ-
ments and help us connect observations to the actual physical
galaxy properties responsible for the emergent spectral proper-
ties (see also Katz et al. 2023a; Choustikov et al. 2024).

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 evol-
ution. 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).
Second, we investigate the origin of the similarities between
observed and simulated spectra and employ the simulation as a
predictive tool for LyC escape ($f_{\text{esc}}^{\text{LyC}}$), comparing these
predictions with both direct and indirect $f_{\text{esc}}^{\text{LyC}}$ 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 interpreta-
tion 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_{\text{esc}}^{\text{LyC}}$ 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.

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

### 2.1. 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 observa-
tions 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_* \sim 10^9 M_\odot$ at
$z=3$, with a maximum cell resolution of 14 pc around this
redshift. We selected 75 outputs with 10 Myr time steps
between $z=4.19$ and $z=3.00$. During this period, the halo
mass increases from $10^{10.6}$ to $10^{10.75} M_\odot$, the stellar mass ($M_*$)
increases from $10^{8.8}$ to $10^{9.4} M_\odot$, both the averaged stellar and
gas metallicity increase from 0.20 $Z_\odot$ to 0.42 $Z_\odot$, and the star
formation rate (SFR) oscillates between $\sim 1.5$ and $4.0 M_\odot \text{ yr}^{-1}$,
with two main SFR peaks around 1.64 and 2.04 Gyr (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), and 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 difference 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:

$$ \tau_{\text{cell}} = \tau_{\text{CII}}^{1334.53} + \tau_{\text{CII}*}^{1335.66} + \tau_{\text{CII}*}^{1335.71} + \tau_{\text{SiII}}^{1260.42} + \tau_{\text{SiII}*}^{1264.73} + \tau_{\text{SiII}*}^{1265.02} + \tau_{\text{dust}}. \quad (1) $$

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 Equation (2) in Mauerhofer et al. (2021). The line parameters are taken from the NIST atomic database (Kramida et al. 2023).<fnref n="6" /> We compute the turbulent velocity of the gas ($v_{\text{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). $\tau_{\text{dust}}$ is set following the implementation of Katz et al. (2023b), 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 postprocessing method within a spherical volume with a radius equivalent to 3 times the galaxy's virial radius ($R_{\text{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$).

## 2.2. 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 4800 Å to 10000 Å. 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 data set consist of star-forming galaxies within the redshift range of $2.4 \le z \le 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 preselect 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 $Ly\alpha$ or ISM lines, so the exact systemic redshift might differ by a few hundred km s$^{-1}$ (less than or equal to the spectroscopic resolution element). We further detail how we account 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 signal-to-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\AA}$) ranges of our final VANDELS subsample 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_\odot$). 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_\odot$).

The SFR and $M_{1250\AA}$ 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_\odot$ and SFRs between $10^{-2}$ to $10^{1.58} M_\odot$ 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, comprise an excellent sample to test the conclusions from our previous study.

<footnote n="6">⁶ https://www.nist.gov/pml/atomic-spectra-database</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** 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 Myr period used to generate the mock spectra.</figcaption></figure>

### 3. Comparing the C II and Si II Spectra in VANDELS and in the Simulation

Here, we compare the simulated C II $\lambda$1334 + C II* $\lambda$1335 and Si II $\lambda$1260 + Si II* $\lambda$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.

#### 3.1. Overview of C II and Si II Measurements

In this section, we measure the EWs of the C II $\lambda$1334 and Si II $\lambda$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 $\lambda$1334 and C II* $\lambda$1335 features (separated by ~1.2 Å, equivalent to half a resolution element) to blend. For comparison, the separation of the Si II $\lambda$1260 and Si II* $\lambda$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 $\lambda$1334 + C II* $\lambda$1335) and EW(Si II $\lambda$1260 + Si II* $\lambda$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⁻¹ around λ = 1334.53 Å and λ = 1260.42 Å, respectively. Spectra dominated by absorption (emission) have positive (negative) EW values. We note that there exists an S II absorption line at 1259.52 Å, blended with Si II $\lambda$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

<figure><img src="image_5.png" /><figcaption>**Figure 2.** The comparison of the EW(C II $\lambda$1334 + C II* $\lambda$1335) and EW(Si II $\lambda$1260 + Si II* $\lambda$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_\odot$), 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.</figcaption></figure>

much weaker than Si II $\lambda$1260 and has a negligible impact on the line measurements.

Figure 2 compares the EW(C II $\lambda$1334 + C II* $\lambda$1335) and EW(Si II $\lambda$1260 + Si II* $\lambda$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 toward 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⁻¹ 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-of-sight 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, presented in 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 data sets.

#### 3.2. Finding the Best-matching Mock Spectra

Comparing the EW(C II $\lambda$1334 + C II* $\lambda$1335) and EW(Si II $\lambda$1260 + Si II* $\lambda$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
<figure><img src="image_6.png" /><figcaption>Figure 3. 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 in 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 χ² 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 χ² spectrum.</figcaption></figure>

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 feature-free 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 χ² value:

$$ 
\begin{align*}
\chi^2 &= \frac{\chi_{C \text{ II}}^2 + \chi_{\text{Si II}}^2}{n_{C \text{ II}}^{\text{obs}} + n_{\text{Si II}}^{\text{obs}}}, \quad \text{where} \\
\chi_{C \text{ II}}^2 &= \sum \left( \frac{F_{C \text{ II}}^{\text{obs}} - F_{C \text{ II}}^{\text{sim}}}{\sigma_{C \text{ II}}^{\text{obs}}} \right)^2 \quad \text{and} \\
\chi_{\text{Si II}}^2 &= \sum \left( \frac{F_{\text{Si II}}^{\text{obs}} - F_{\text{Si II}}^{\text{sim}}}{\sigma_{\text{Si II}}^{\text{obs}}} \right)^2. \tag{2}
\end{align*}
$$

$F_{\text{ion}}^{\text{obs}}$ and $F_{\text{ion}}^{\text{sim}}$ are the observed and simulated flux in the wavelength regions of each ion, respectively. $\sigma_{\text{ion}}^{\text{obs}}$ represents the error on the observed flux, and $n_{\text{ion}}^{\text{obs}}$ 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 prealignment 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 χ² 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, 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 χ² for all 22,500 mock spectra, to explore the diversity of spectral shapes 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⁹.², 10⁹.⁶, and 10¹⁰.⁶ M☉. For each galaxy, we show the best-matching mock spectrum (i.e., with the lowest χ² value). The best-matching spectrum replicates reasonably well the line profiles in the two lower-mass galaxies (χ² < 1.5). However, they fall short of capturing the complete width of the broad absorption lines observed in the 10¹⁰.⁶ M☉ galaxy (χ² > 3). This result reaffirms the findings from Blaizot et al. (2023) and Gazagnes et al. (2023), emphasizing that mock line spectra generated from a 10⁹.² M☉ virtual galaxy do not
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Histogram of the best-match $\tilde{\chi}^2$ for the 131 VANDELS spectra analyzed in this work. Of these spectra, 37% are well reproduced ($\tilde{\chi}^2 < 1$) and 83% are reasonably matched ($\tilde{\chi}^2 < 2$).</figcaption></figure>

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 $\tilde{\chi}^2$ for the whole sample. The top panel presents the distribution of the $\tilde{\chi}^2$ value of the best-matching spectrum for each object. Of these spectra, 37% are well reproduced ($\tilde{\chi}^2 < 1$) and 83% are reasonably matched ($\tilde{\chi}^2 < 2$). These numbers are slightly lower than for the CLASSY sample of $z \sim 0$ galaxies (60% with $\tilde{\chi}^2 < 1$, 89% with $\tilde{\chi}^2 < 2$). These results may be explained by the presence of more massive galaxies in the selected VANDELS sample whose complex spectra are more difficult 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 ($\tilde{\chi}^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 in 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 Figure 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 an SFR more similar to those of the virtual galaxy. To do so, we analyze the distribution of the average number of mock spectra with $\tilde{\chi}^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 $\tilde{\chi}^2$ values less than 1, 2, and 3 per galaxy across increasing mass bins. The first bin ($M_* < 10^{9.3} M_\odot$) directly aligns with the stellar mass range of the virtual galaxy ($\sim 10^{8.9}-10^{9.3} M_\odot$). Interestingly, this bin exhibits the highest average number of mock spectra with $\tilde{\chi}^2$ values less than 1, 2, and 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_\odot$, there are, on average, five to eight mock spectra with $\tilde{\chi}^2 < 1$, while none are found for galaxies in the two largest stellar mass bins ($\sim > 10^{10} M_\odot$). These results echo those of Blaizot et al. (2023) and Gazagnes et al. (2023), indicating that the current simulation may not accurately represent galaxies with $M_* > 10^{10} M_\odot$. Furthermore, we observe a declining trend in the average number of mock spectra with $\tilde{\chi}^2 < 2$ and $\tilde{\chi}^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 ($\sim 10^0-10^{0.6} M_\odot \text{ 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 $\tilde{\chi}^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 $\tilde{\chi}^2 < 1$ from the lowest SFR bin to the largest SFR bin. However, we do discern a slight decrease in the average number of $\tilde{\chi}^2 < 2$ and $\tilde{\chi}^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. 2015, 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
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Top: histogram presenting the average number of mock spectra with $\bar{\chi}^2 < 1$ (black), $\bar{\chi}^2 < 2$ (gray), and $\bar{\chi}^2 < 3$ (white) per observed galaxy, in bins of stellar mass. The most massive VANDELS galaxies ($\sim > 10^{10}$) are not well replicated by the mock spectra for our virtual galaxy. Bottom: same, but for the SFR.</figcaption></figure>

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, Figure 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 does not 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

<figure><img src="image_9.png" /><figcaption>**Figure 6.** The evolution of the SFR (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 Myr 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 2790 best-matching mock spectra are found (taking the 30 best mocks for each of the 93 VANDELS galaxies selected in Section 3.2). Overall, the best-matching spectra predominantly originate from luminosity peaks closely preceding the two most intense star formation bursts.</figcaption></figure>

simulations can accurately reproduce as a function of their properties, this effort extends beyond the scope of the present paper.

## 4.2. Best-matching Spectra Originate from Luminosity Peaks Close to Star Formation Bursts

An important and novel aspect of our analysis involves examining the average distribution of $\bar{\chi}^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 $\bar{\chi}^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 2790 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\AA}$) 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\AA}$, and most
of the best mocks originate from bright luminosity peaks. Interestingly, these peaks typically precede strong star formation 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 Myr. As a result, the luminosity peaks may coincide with SFR peaks, albeit with a delay in its influence on the 100 Myr 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 Gyr or the 2.02 Gyr time step, corresponding to the two highest peaks in the histogram distribution near the peaks of SFR. This finding has two implications. First, it indicates that a comparable success rate in matching spectra to observations could have been achieved by solely considering these two time steps. Second, given the diversity of line properties exhibited by the 93 VANDELS galaxies (shown in Figure 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 are not 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 fourfold 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 toward UV-bright galaxies (Garilli et al. 2021) likely in a phase of star formation (a bias that we further accentuated by selecting the highest-S/N 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.

## 5. 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_{\text{esc}}^{\text{LyC}}$ (e.g., $O_{32}$, $Ly\alpha$ Mg II, $\beta$ slopes, or LIS covering fractions; see Jaskot & Oey 2013; Reddy et al. 2016; Verhamme et al. 2017; Gazagnes et al. 2018; Chisholm et al. 2020; Flury et al. 2022a; Chisholm et al. 2022; Saldana-Lopez et al. 2023; Izotov et al. 2024). Although empirical trends between $f_{\text{esc}}^{\text{LyC}}$ 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_{\text{esc}}^{\text{LyC}}$ using the simulation and compare these predictions with previous constraints from Begley et al. (2022, 2024) and Saldana-Lopez et al. (2023). Section 5.1 details our simulation-based approach to predict $f_{\text{esc}}^{\text{LyC}}$ 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.

### 5.1. Extracting Simulation-based $f_{\text{esc}}^{\text{LyC}}$ Predictions

The RHD simulation used in this work includes a self-consistent 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_{\text{esc}}^{\text{LyC}}$ 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 encounters 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 divide it by the intrinsic number produced to get $f_{\text{esc}}^{\text{LyC}}$. 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 an $f_{\text{esc}}^{\text{LyC}}$ prediction for each of the 93 VANDELS best-matched objects, we take the median and standard deviation of the 30 $f_{\text{esc}}^{\text{LyC}}$ values corresponding to the 30 best-matching mock spectra. To differentiate our predictions from direct LyC measurements, we refer to these predicted $f_{\text{esc}}^{\text{LyC}}$ as $f_{\text{esc}}^{\text{pred}}$.

Figure 7 presents the distribution of the 93 $f_{\text{esc}}^{\text{LyC}}$ measurements. The vast majority of values fall below 1% escape, and only two objects have $f_{\text{esc}}^{\text{LyC}}$ > 10%. The
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The distribution of $f_{\text{esc (pred)}}^{\text{LyC}}$ for the VANDELS sample analyzed in this work. The $f_{\text{esc (pred)}}^{\text{LyC}}$ is determined by taking the average $f_{\text{esc}}^{\text{LyC}}$ value across the 30 best mock spectra for each of the 93 selected VANDELS objects. The average $f_{\text{esc (pred)}}^{\text{LyC}}$ for this sample is $0.01 \pm 0.02$.</figcaption></figure>

average $f_{\text{esc (pred)}}^{\text{LyC}}$ for the 93 galaxies is $0.01 \pm 0.02$. While lower, this estimate falls within $2\sigma$ of the $0.07 \pm 0.02$ escape fraction reported by Begley et al. (2022), who directly constrained $f_{\text{esc}}^{\text{LyC}}$ using ultradeep, ground-based, U-band imaging. Interestingly, our average predicted escape fraction is highly consistent with that of Saldana-Lopez et al. (2023), who determined an average $f_{\text{esc}}^{\text{LyC}}$ of $0.02 \pm 0.01$ using 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 Å (Equation (5) from Chisholm et al. 2018). The consistency of both approaches is interesting, in particular 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 the $f_{\text{esc}}^{\text{LyC}}$ of actual observations (see also Katz et al. 2020, 2023a, 2023b; Choustikov et al. 2023). In this study, the primary limitation stems from relying on a single zoom-in simulation to forecast $f_{\text{esc}}^{\text{LyC}}$. This approach inherently restricts the range of potential $f_{\text{esc}}^{\text{LyC}}$ values to those observed within the simulation itself, with the maximum $f_{\text{esc}}^{\text{LyC}}$ 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_{\text{esc}}^{\text{LyC}}$ 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$\alpha$ and UV spectral properties.

## 5.2. Testing the Consistency of the Predicted $f_{\text{esc}}^{\text{LyC}}$

Here, we assess if the $f_{\text{esc (pred)}}^{\text{LyC}}$ derived in the previous section correlates with other established LyC diagnostics.

### 5.2.1. Comparing to the Ly$\alpha$ Escape Fractions

The Ly$\alpha$ line is a powerful tracer for the LyC leakage (e.g., Verhamme et al. 2017; Steidel et al. 2018; Izotov et al. 2024), as Ly$\alpha$ photons are affected by the same neutral hydrogen gas and dust extinction regulating the escape of LyC photons.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** The comparison of the $f_{\text{esc (pred)}}^{\text{LyC}}$ derived using the simulation and the Ly$\alpha$ 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. (2024; using a $z < 0.5$ LyC emitter sample) and in Begley et al. (2024; using the VANDELS sample), respectively.</figcaption></figure>

Begley et al. (2024) investigated the connection between the escape fraction of Ly$\alpha$ ($f_{\text{esc}}^{\text{Ly}\alpha}$) and $f_{\text{esc}}^{\text{LyC}}$ 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_{\odot}$ and SFR values between $10^{0.8}$ and $10^{1.8} M_{\odot}$ yr$^{-1}$. Further, the median $\chi^2$ of their 30 best mock spectra spans values from 0.7 to 1.9. This indicates that these 11 objects represent the broader well-matched 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_{\text{esc}}^{\text{Ly}\alpha}$ values measured in Begley et al. (2024) for these objects. We observe a significant correlation between $f_{\text{esc (pred)}}^{\text{LyC}}$ and the measured $f_{\text{esc}}^{\text{Ly}\alpha}$. Notably, this trend is highly consistent with two empirical $f_{\text{esc}}^{\text{LyC}}-f_{\text{esc}}^{\text{Ly}\alpha}$ relationships, one derived in Begley et al. (2024) using VANDELS galaxies, and the other established in Izotov et al. (2024) 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$\alpha$ LyC observational trends. This result underscores the consistency of both Ly$\alpha$ and LIS metal lines as indicators of $f_{\text{esc}}^{\text{LyC}}$. Further, it gives us confidence in the $f_{\text{esc (pred)}}^{\text{LyC}}$, 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_{\text{esc}}^{\text{Ly}\alpha}$ in VANDELS. Here, we explore a similar exercise, but producing composite stacks in three bins of $f_{\text{esc (pred)}}^{\text{LyC}}$: below 1%, between 1% and 10%, and larger than or equal to 10%, respectively. We use these three $f_{\text{esc (pred)}}^{\text{LyC}}$ bins to represent categories of non, weak, or relatively strong leakers, respectively. As discussed in Section 5.1, only two objects have $f_{\text{esc}}^{\text{LyC}} \ge 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_{\text{esc (pred)}}^{\text{LyC}}+1\sigma > 10\%$).

Figure 9 presents these three stacks, with four panels zooming in on the Ly$\alpha$ profiles, the LIS absorption lines (from Si II, O I, and C II), and the high-energy UV features C IV and
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Composite stacks of VANDELS spectra based on the $f_{\text{esc}}^{\text{LyC}}$ (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 $\beta$ slope using the observed continuum flux between 1250 and 1700 Å, and the four panels on the bottom zoom in on the Ly$\alpha$ profiles, LIS absorption lines, and high-energy C IV and He II features.</figcaption></figure>

He II, tracing the hard ionizing radiation field above 47.9 and 54.4 eV respectively. For each of the stacks, we measure the $\beta$ slope between 1250 and 1700 Å (selecting feature-free intervals in that range) by fitting a line such that $F_{\text{obs}} \propto \lambda^{\beta}$. $\beta$ 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. $\beta$ is especially powerful because its measurement can be relatively simply extracted from spectra of high-$z$ objects with JWST (e.g., Cullen et al. 2023, 2024).

Examining the distinctions among the three composite spectra, we observe several trends. Stacks with higher predicted $f_{\text{esc}}^{\text{LyC}}$ exhibit bluer $\beta$ values, an increase in Ly$\alpha$ 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$\alpha$ photons (Gazagnes et al. 2020; Begley et al. 2024). We note that the measured $\beta$ slopes in the weak and relatively strong leakers are “redder” than what is typically expected from low-$z$ LyC emitters. Using the $\beta$ to $f_{\text{esc}}^{\text{LyC}}$ relation from Chisholm et al. (2022) derived from the LzLCS sample (Flury et al. 2022b), $\beta$ slopes of -1.35 and -1.84 would correspond to $f_{\text{esc}}^{\text{LyC}}$ of ~1% and ~3%, respectively, and hence would be 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_{\text{esc}}^{\text{LyC}}$ > 10% stack, which could be due to either 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_{\text{esc}}^{\text{LyC}}$ > 20% or $f_{\text{esc}}^{\text{Ly}\alpha}$ > 50%; see Naidu et al. 2022; Schaeerer et al. 2022) with large intrinsic production of ionizing photons. None of the $f_{\text{esc}}^{\text{LyC}}$ derived in this work are consistent with the

presence of an extreme leakage.<fnref n="7" /> 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 nonleakers 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. 2020, 2023a, 2023b; Choustikov et al. 2023), and we explore that aspect in the next section.

### 5.3. What Makes LIS Lines Good Tracers of $f_{\text{esc}}^{\text{LyC}}$

In the previous section, we highlighted that the $f_{\text{esc}}^{\text{LyC}}$ 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$\alpha$ and the $\beta$ slopes, and they 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).

Yet, this outcome brings us to the following question: what makes the LIS metal lines good tracers of the LyC leakage? The tight $f_{\text{esc}}^{\text{LyC}}$ empirical trends found in low-$z$ galaxies over

<footnote n="7">⁷ We note that Saxena et al. (2020) studied the properties of 51 He II emitters in VANDELS but none of these objects are present in our current sample.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** $f_{\text{esc}}^{\text{LyC}}$ 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$ Gyr. Squares and circles differentiate sightlines where less or more than 75% of the intrinsic flux escapes the galaxy, respectively (using the luminosity-weighted $F_{\text{obs}}/F_{\text{int}}$ ratio). Each point is color-coded by $r_{75}$, which is the radius at which the integrated flux accounts for more than 75% of the total escaping flux in that direction.</figcaption></figure>

the past suggested that LyC photons escape through relatively dust-free, low neutral gas column density channels (Reddy et al. 2016; Chisholm et al. 2018; Gazagnes et al. 2018; Saldana-Lopez et al. 2023). Further, Flury et al. (2022a) found a significant correlation between larger LyC leakage and the objects’ compactness, which may support the notion that stellar feedback more efficiently clears escape channels in such environments (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_{\text{esc}}^{\text{LyC}}$, 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$ Gyr (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 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 half-light 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_{\text{esc}}^{\text{LyC}}$ in the simulation, EW(LIS), $r_{75}$, and the dust attenuation. Directions with high $f_{\text{esc}}^{\text{LyC}}$ exhibit low EW(LIS), minimal dust attenuation, and a smaller $r_{75}$. 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 $F_{\text{obs}}/F_{\text{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_{\text{esc}}^{\text{LyC}}$ 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_{\text{esc}}^{\text{LyC}}$. 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 its reliability and the physics behind existing LyC diagnostics. This effort will be the topic of a separate paper.

## 6. 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_{\odot}$ virtual galaxy to the observed spectra of 131 $z \sim 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 LyC leakage within star-forming galaxies. In particular, we established a simulation-based $f_{\text{esc}}^{\text{LyC}}$ estimation method that enabled us to predict the strength of the LyC leakage in the selected VANDELS sample. We summarize our findings as follows:

1. 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 toward 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).

2. For each VANDELS galaxy, we determine the similarity ($\tilde{\chi}^2$) between the mock spectra and the observed absorption and emission line profiles. We find for 37% of the VANDELS spectra a mock spectrum with excellent agreement ($\tilde{\chi}^2 < 1$), while 83% are reasonably well matched (best mock spectrum with $\tilde{\chi}^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_{\odot}$; Section 4.1) due to the presence of broader absorption profiles. This finding further supports the conclusions from Blaizot et al. (2023) and Gazagnes et al. (2023), highlighting the limitations of the current virtual galaxy in representing the environments of the most massive galaxies.

3. 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 4 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).

4. We predicted the $f_{\text{esc}}^{\text{LyC}}$ for each galaxy in the selected VANDELS sample using the virtual galaxy environment ($f_{\text{esc}}^{\text{LyC}}$; see Section 5.1). We derived an average predicted $f_{\text{esc}}^{\text{LyC}}$ of $0.01 \pm 0.02$, within $2\sigma$ of the value $0.07 \pm 0.02$ reported 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 Å.

5. The simulation-based $f_{\text{esc}}^{\text{LyC}}$ tightly correlate with the Ly$\alpha$ escape fractions from Begley et al. (2024) and are in remarkable agreement with the empirical $f_{\text{esc}}^{\text{LyC}} - f_{\text{esc}}^{\text{Ly}\alpha}$ relations from Izotov et al. (2024) and Begley et al. (2024) (Section 5.2.1, Figure 8). We produced composite spectra in bins of $f_{\text{esc}}^{\text{LyC}}$ (Section 5.2.2) and show that the stack comprising the higher $f_{\text{esc}}^{\text{LyC}}$ galaxies exhibit bluer $\beta$ values, an increase in Ly$\alpha$ flux, and weaker LIS absorption lines. We do not find significant differences in high-energy lines such as C IV and He II, which may suggest that the current sample does not hold extreme LyC leakage candidates or that the relation between $f_{\text{esc}}^{\text{LyC}}$ and spectral hardness is more complex than expected.

6. Building upon the good agreement between the simulation-based $f_{\text{esc}}^{\text{LyC}}$ and well-established LyC empirical trends, we explored what makes LIS lines good tracers of $f_{\text{esc}}^{\text{LyC}}$ in the simulation (Section 5.3). We find that directions with larger $f_{\text{esc}}^{\text{LyC}}$ exhibit lower LIS absorption line EW, lower dust attenuation, and compact spatial profiles (smaller $r_{75}$). These results align with the expected LyC-leakage picture derived from observations where LyC photons likely escape through spatially confined areas depleted of neutral gas and dust. The consistent co-occurrence of these three characteristics (low neutral gas, low dust, and spatial compactness) may explain why LIS-based LyC diagnostics provide relatively robust indirect estimates into the $f_{\text{esc}}^{\text{LyC}}$ of galaxies.

This study extends the success of using a zoom-in simulation approach to analyze UV absorption and emission lines, as seen in recent studies (Mauerhofer et al. 2021; Blaizot et al. 2023; Gazagnes et al. 2023). It reaffirms the potential of these methods in understanding the factors contributing to various and intricate line profiles in galaxies. Additionally, it underscores their relevance in the context of high-redshift observations (see also Choustikov et al. 2023; Katz et al. 2023a) and their potential to provide valuable insights, including estimates of the LyC escape fraction, in the current epoch of reionization-era observations with JWST.

## Acknowledgments

The authors would like to thank the referee for comments that helped improve the quality of this manuscript. S.G. is grateful for the support enabled by the Harlan J. Smith McDonald fellowship, which made this project possible. F.C. acknowledges the support from a UKRI Frontier Research Guarantee Grant (PI: Cullen; grant reference EP/X021025/1). T.G. and A.V. are supported by the professorship grant SNFProf_PP00P2_211023 from the Swiss National Foundation for Scientific Research. V.M. acknowledges support from the NWO grant 016.VIDI.189.162 ("ODIN").

*Software: Astropy* (Astropy Collaboration et al. 2013, 2018, 2022), NumPy (Harris et al. 2020), pandas (McKinney 2010), SciPy (Virtanen et al. 2020), RASCAS (Michel-Dansac et al. 2020).

## Data Availability

The VANDELS survey is a European Southern Observatory Public Spectroscopic Survey. The full spectroscopic data set, together with photometric catalogs and derived quantities, is available from http://vandels.inaf.it, as well as from the ESO archive https://www.eso.org/qi. Any secondary product and/or data underlying this article will be shared on reasonable request to the corresponding author.

## Appendix

### Accounting for a Redshift Uncertainty When Finding the Best-matching Spectra

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 kinematics-dependent and therefore result in slight offsets (typically within a few hundred km s⁻¹ 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 prealignment 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 $\tilde{\chi}^2$ values and highlights that a few more galaxies have $\tilde{\chi}^2 < 1.5$ when we first prealign the mock and observed spectra. This result is not unexpected, as by performing a prealignment, 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-step 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
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Exploring the impact of accounting for potential redshift uncertainties when finding the best-matching mock spectra. Top: The comparison of the best-match $\tilde{\chi}^2$ for the two different cases. As expected, allowing for a spectral shift when finding the best-matching spectra typically leads to better $\tilde{\chi}^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 either distribution. Right: Comparing the predicted $f_{\text{esc}}^{\text{LyC}}$ when using the two different approaches. We do observe small variations of the $f_{\text{esc}}^{\text{LyC}}$, however, these variations are within the uncertainties and do not change the interpretation of galaxies that could be categorized as non, weak, or strong leakers. Overall, these panels support that choosing one or the other approach has no impact on the main conclusions drawn in this work.</figcaption></figure>

panel compares the differences in the predicted $f_{\text{esc}}^{\text{LyC}}$, using the methodology detailed in Section 5, when following 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 \pm 0.02$, in both cases).

Overall, Figure 11 emphasizes that choosing a best-matching approach that does not include a prealignment step would have had no impact on the main conclusions drawn in this work.

## ORCID iDs

Simon Gazagnes https://orcid.org/0000-0002-5659-4974

Fergus Cullen https://orcid.org/0000-0002-3736-476X

Valentin Mauerhofer https://orcid.org/0000-0003-0595-9483

Ryan Begley https://orcid.org/0000-0003-0629-8074

Danielle Berg https://orcid.org/0000-0002-4153-053X

Jeremy Blaizot https://orcid.org/0000-0003-1609-7911

John Chisholm https://orcid.org/0000-0002-0302-2577

Thibault Garel https://orcid.org/0000-0002-9613-9044

Floriane Leclercq https://orcid.org/0000-0002-6085-5073

Anne Verhamme https://orcid.org/0000-0002-2201-1865

## References

Alexandroff, R. M., Heckman, T. M., Borthakur, S., Overzier, R., & Leitherer, C. 2015, ApJ, 810, 104

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Badnell, N. R. 2006, ApJS, 167, 334

Begley, R., Cullen, F., McLure, R. J., et al. 2022, MNRAS, 513, 3510

Begley, R., Cullen, F., McLure, R. J., et al. 2024, MNRAS, 527, 4040

Berg, D. A., Chisholm, J., Erb, D. K., et al. 2019, ApJL, 878, L3

Berg, D. A., James, B. L., King, T., et al. 2022, ApJS, 261, 41

Blaizot, J., Garel, T., Verhamme, A., et al. 2023, MNRAS, 523, 3749

Carr, C., Michel-Dansac, L., Blaizot, J., et al. 2023, ApJ, 952, 88

Carr, C., & Scarlata, C. 2022, ApJ, 939, 11

Carr, C., Scarlata, C., Panagia, N., & Henry, A. 2018, ApJ, 860, 143

Cen, R. 2020, ApJL, 889, L22

Chisholm, J., Gazagnes, S., Schaerer, D., et al. 2018, A&A, 616, A30

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

Chisholm, J., Saldana-Lopez, A., Flury, S., et al. 2022, MNRAS, 517, 5104

Chisholm, J., Tremonti, C. A., Leitherer, C., et al. 2015, ApJ, 811, 149

Chisholm, J., Tremonti, C. A., Leitherer, C., & Chen, Y. 2017, MNRAS, 469, 4831

Chisholm, J., Tremonti, C. A., Leitherer, C., Chen, Y., & Wofford, A. 2016, MNRAS, 457, 3133

Choustikov, N., Katz, H., Saxena, A., et al. 2024, MNRAS, 529, 3751

Cullen, F., McLeod, D. J., McLure, R. J., et al. 2024, MNRAS, 531, 997

Cullen, F., McLure, R. J., McLeod, D. J., et al. 2023, MNRAS, 520, 14

De Cia, A., Ledoux, C., Mattsson, L., et al. 2016, A&A, 596, A97

Erb, D. K., Quider, A. M., Henry, A. L., & Martin, C. L. 2012, ApJ, 759, 26

Finley, H., Bouché, N., Contini, T., et al. 2017, A&A, 605, A118

Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022a, ApJ, 930, 126

Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022b, ApJS, 260, 1

Garilli, B., McLure, R., Pentericci, L., et al. 2021, A&A, 647, A150

Gazagnes, S., Chisholm, J., Schaerer, D., et al. 2018, A&A, 616, A29

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

Gazagnes, S., Mauerhofer, V., Berg, D. A., et al. 2023, ApJ, 952, 164

Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386

Hahn, O., & Abel, T., 2013 MUSIC: MUIti-Scale Initial Conditions, Astrophysics Source Code Library, ascl:1311.011

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Heckman, T. M., Alexandroff, R. M., Borthakur, S., Overzier, R., & Leitherer, C. 2015, ApJ, 809, 147

Heckman, T. M., Lehnert, M. D., Strickland, D. K., & Armus, L. 2000, ApJS, 129, 493

Heckman, T. M., Sembach, K. R., Meurer, G. R., et al. 2001, ApJ, 558, 56

Izotov, Y. I., Thuan, T. X., Guseva, N. G., et al. 2024, MNRAS, 527, 281

James, B. L., Aloisi, A., Heckman, T., Sohn, S. T., & Wolfe, M. A. 2014, ApJ, 795, 109

Jaskot, A. E., Dowd, T., Oey, M. S., Scarlata, C., & McKinney, J. 2019, ApJ, 885, 27

Jaskot, A. E., & Oey, M. S. 2013, ApJ, 766, 91

Katz, H., Durovčíková, D., Kimm, T., et al. 2020, MNRAS, 498, 164

Katz, H., Rosdahl, J., Kimm, T., et al. 2023a, OJAp, 6, 44

Katz, H., Saxena, A., Rosdahl, Y., et al. 2023b, MNRAS, 518, 270
Konstantopoulou, C., De Cia, A., Krogager, J.-K., et al. 2022, A&A, 666, A12

Kramida, A., Ralchenko, Y., Reader, J. & NIST ASD Team, 2023 NIST Atomic Spectra Database, version 5.11, doi:10.18434/T4W30F

Li, Z., Gronke, M., & Steidel, C. C. 2024, MNRAS, 529, 444

Marques-Chaves, R., Schaerer, D., Amorín, R. O., et al. 2022, A&A, 663, L1

Mauerhofer, V., Verhamme, A., Blaizot, J., et al. 2021, A&A, 646, A80

McKinney, W. 2010, in SciPy 2010: Proc. 9th Python in Science Conference, 56

McLure, R. J., Pentericci, L., Cimatti, A., et al. 2018, MNRAS, 479, 25

Michel-Dansac, L., Blaizot, J., Garel, T., et al. 2020, A&A, 635, A154

Naidu, R. P., Matthee, J., Oesch, P. A., et al. 2022, MNRAS, 510, 4582

Pentericci, L., McLure, R. J., Garilli, B., et al. 2018, A&A, 616, A174

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

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31

Rivera-Thorsen, T. E., Dahle, H., Gronke, M., et al. 2017, A&A, 608, L4

Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188

Rosdahl, J., Katz, H., Blaizot, J., et al. 2018, MNRAS, 479, 994

Rosdahl, J., & Teyssier, R. 2015, MNRAS, 449, 4380

Rubin, K. H. R., Prochaska, J. X., Ménard, B., et al. 2011, ApJ, 728, 55

Rupke, D. S., Veilleux, S., & Sanders, D. B. 2005, ApJS, 160, 115

Saldana-Lopez, A., Schaerer, D., Chisholm, J., et al. 2022, A&A, 663, A59

Saldana-Lopez, A., Schaerer, D., Chisholm, J., et al. 2023, MNRAS, 522, 6295

Saxena, A., Pentericci, L., Mirabelli, M., et al. 2020, A&A, 636, A47

Scarlata, C., & Panagia, N. 2015, ApJ, 801, 43

Schaerer, D., Izotov, Y. I., Worseck, G., et al. 2022, A&A, 658, L11

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

Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123

Teyssier, R. 2002, A&A, 385, 337

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

Trebitsch, M., Dubois, Y., Volonteri, M., et al. 2021, A&A, 653, A154

Verhamme, A., Orlitová, I., Schaerer, D., et al. 2017, A&A, 597, A13

Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Voronov, G. S. 1997, ADNDT, 65, 1

Wang, B., Heckman, T. M., Zhu, G., & Norman, C. A. 2020, ApJ, 894, 149

Witten, C. E. C., Laporte, N., & Katz, H. 2023, ApJ, 944, 61

Wofford, A., Leitherer, C., & Salzer, J. 2013, ApJ, 765, 118

Xu, X., Heckman, T., Henry, A., et al. 2022, ApJ, 933, 222