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

# Constraining Reionization with Ly$\alpha$ Damping-wing Absorption in Galaxy Spectra: A Machine Learning Model Based on Reionization Simulations

Hyunbae Park¹, Intae Jung², Hidenobu Yajima¹, Jenny G. Sorce³,⁴,⁵, Paul R. Shapiro⁶, Kyungjin Ahn⁷, Pierre Ocvirk⁸, Romain Teyssier⁹, Gustavo Yepes¹⁰,¹¹, Ilian T. Iliev¹², and Joseph S. W. Lewis¹³

¹Center for Computational Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan; hyunbae.park@gmail.com

²Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

³Université de Lille, CNRS, Centrale Lille, UMR 9189 CRISTAL, F-59000 Lille, France

⁴Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405, Orsay, France

⁵Leibniz-Institut für Astrophysik, An der Sternwarte 16, 14482 Potsdam, Germany

⁶Department of Astronomy, University of Texas, Austin, TX 78712-1083, USA

⁷Department of Earth Sciences, Chosun University, Gwangju 61452, Republic of Korea

⁸Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS UMR 7550, 11 rue de l'Université, 67000 Strasbourg, France

⁹Institute for Theoretical Physics, University of Zurich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland

¹⁰Departamento de Física Teórica M-8, Universidad Autónoma de Madrid, Cantoblanco, 28049, Madrid, Spain

¹¹Centro de Investigación Avanzada en Física Fundamental (CIAFF), Universidad Autónoma de Madrid, 28049 Madrid, Spain

¹²Astronomy Centre, Department of Physics & Astronomy, Pevensey III Building, University of Sussex, Falmer, Brighton, BN1 9QH, UK

¹³Institut d'Astrophysique de Paris, UMR 7095, CNRS, UPMC Univ. Paris VI, 98 bis boulevard Arago, 75014 Paris, France

Received 2024 October 9; revised 2025 February 10; accepted 2025 March 11; published 2025 April 10

## Abstract

Recently, NIRSpec PRISM/CLEAR observations by JWST have begun providing rest-frame UV continuum measurements of galaxies at $z \gtrsim 7$, revealing signatures of Ly$\alpha$ damping-wing (DW) absorption by the intergalactic medium (IGM). We develop a methodology to constrain the global ionization fraction of the IGM ($Q_{HII}$) using low-resolution spectra, employing the random forest classification (RFC) method. We construct mock spectra using the simulated galaxies and the IGM from the Cosmic Dawn II simulation and train RFC models to estimate $Q_{HII}$ at the redshift of the source and to detect the presence of a damped Ly$\alpha$ absorber. We find that individual galaxy spectra with spectral bins between 1220 and 1270 Å and with signal-to-noise ratios (SNR) >20 can place tight constraints on $Q_{HII}$, provided the UV continuum is accurately modeled. This method is particularly effective for the early phase of reionization ($Q_{HII} < 50\%$), when the IGM opacity is high in the DW. As a demonstration, we apply our model to existing NIRSpec PRISM/CLEAR spectra, placing upper bounds of $Q_{HII} = 59.6\%$, 5.6\%, and 18.5\% at $z = 7.7, 9.4$, and 10.6, respectively, with 68\% confidence, though several modeling uncertainties remain to be discussed. These constraints favor late-starting reionization models, where $\gtrsim 80\%$ of the IGM is ionized after $z = 8$. We conclude that high-SNR observations of carefully selected targets around $z \sim 7-9$ can effectively constrain reionization models.

Unified Astronomy Thesaurus concepts: Reionization (1383); Early universe (435); Origin of the universe (1186); Random Forests (1935); Radiative transfer simulations (1967)

## 1. Introduction

One of the pivotal epochs in the quest to understand cosmic evolution is the reionization era, during which early galaxies and black holes ionized the neutral hydrogen in intergalactic space. The ionizing efficiency of these sources is a key unknown in star formation physics, essential for understanding the subsequent evolution of the cosmos. In recent years, advances in both observational and computational capabilities have enabled us to probe this era in diverse ways.

Galactic radiation near Ly$\alpha$ is subject to attenuation by the neutral intergalactic medium (IGM) and is, therefore, considered an effective probe of the H II regions created by early galaxies (see reviews by, e.g., M. Ouchi et al. 2020). The equivalent width (EW) of the Ly$\alpha$ emission line has been widely used in recent studies to estimate the IGM opacity to Ly$\alpha$, and considerable progress has been made in understanding the ionization structure in the IGM up to $z \sim 8$ with Ly$\alpha$ EW statistics and Ly$\alpha$

luminosity function (E. Vanzella et al. 2011; Y. Ono et al. 2012; T. Shibuya et al. 2012; T. Treu et al. 2012, 2013; S. L. Finkelstein et al. 2013; M. A. Schenker et al. 2014; A. Zitrin et al. 2015; M. Song et al. 2016; I. Jung et al. 2018, 2020, 2022; C. A. Mason et al. 2018; L. Pentericci et al. 2018; A. Hoag et al. 2019; Y. Tilvi et al. 2020; L. R. Whitler et al. 2020; A. M. Morales et al. 2021; P. Bolan et al. 2022; R. Endsley & D. P. Stark 2022; R. L. Larson et al. 2022; I. G. B. Wold et al. 2022; S. Bruton et al. 2023; T. Morishita et al. 2023; M. Nakane et al. 2024).

A challenge in utilizing Ly$\alpha$ EW is that both the intrinsic emission strength and IGM opacity can vary significantly from galaxy to galaxy and sight line to sight line due to the complexity of the internal structure of galaxies and density/velocity fluctuations in the IGM (e.g., A. Hutter et al. 2014; H. Park et al. 2021; A. Smith et al. 2022a). The sight-line variation causes a large scatter in the relationship between Ly$\alpha$ EW and global ionization fraction, requiring a large number of galaxies to accurately constrain the global ionization fraction. Additionally, the emergent Ly$\alpha$ emission strength before the IGM attenuation is highly uncertain and is expected to depend on the viewing angle, further increasing the scatter (A. Smith et al. 2022b; J. Blaizot et al. 2023). For instance, the JADES and CEERS surveys by the James Webb

<figure><img src="image_2.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.
Space Telescope (JWST; J. P. Gardner et al. 2006) found significantly different fractions of strong Lyα emitters (LAEs; EW > 25 Å) from 25 galaxies each (20% versus 52%). This discrepancy can be explained by the statistical variation in Lyα opacity calculated from numerically simulated IGM (L. Napolitano et al. 2024). It has been suggested that the number of galaxies required to constrain IGM opacity to within 10% accuracy is larger than 100 (H. Park et al. 2021).

In this work, we focus on the Lyα damping-wing (DW) absorption in the continuum of Lyman break galaxies (LBGs) as an alternative probe that is free from some of the limitations of the Lyα EW statistics. When the Universe was highly neutral, the Lyα DW can produce a significant absorption feature up to ~100 Å redward of the resonance (J. Miralda-Escudé 1998). The DW feature from the high-z Universe was typically searched for in bright active galactic nuclei (AGNs) spectra, as their strong flux provides high signal-to-noise ratio (SNR) measurements (D. J. Mortlock et al. 2011; E. Bañados et al. 2018; F. B. Davies et al. 2018; F. Wang et al. 2020; J. Yang et al. 2020; J. F. Hennawi et al. 2024; T. Kist et al. 2024). Despite the lower flux, DWs from LBGs is potentially useful for constraining the reionization era as they are much more numerous than quasars, and their narrow emission lines, no broader than a few hundred kilometers per second, allow easier and more accurate fitting of the continuum component to extract the IGM transmissivity (K. E. Heintz et al. 2025).

A clean measurement of galactic continuum around Lyα is challenging for galaxies at $z > 6$ with ground-based telescopes, as their spectra mostly fall outside the visible window of Earth's atmosphere. It is, therefore, advantageous to utilize space instruments such as JWST to acquire damping-wing absorption data from the reionization era. Recent observations by NIRSpec/PRISM on JWST have yielded clean measurements of Lyα DW from high-z galaxies (P. Arrabal Haro et al. 2023; A. J. Bunker et al. 2023; K. Boyett et al. 2024; J. Witstok et al. 2025), and several studies have produced constraints on the reionization history using these data (E. Curtis-Lake et al. 2023; T. Y.-Y. Hsiao et al. 2023; K. E. Heintz et al. 2025; H. Umeda et al. 2024; C. A. Mason et al. 2025). A significant portion of these high-z observations exhibits a strong damped Lyα absorber (DLA) component added to the IGM component, presenting a challenge in separating the IGM opacity.

In this study, we develop a methodology to constrain the reionization history using LBG spectra, combining a reionization simulation with a machine learning (ML) algorithm. Cosmological simulations capture the complexity of density/velocity structures in the IGM, providing a realistic database of mock DW spectra. Using reionization simulations is advantageous over purely analytic models, as simulated ionization fields yield the possible range of H II bubble sizes around the source galaxy at a given global IGM ionization fraction (H. Chen 2024; L. C. Keating et al. 2024), providing a prior that can tighten the constraint.

The rest of the paper is as follows. In Section 2, we describe how we prepare the mock Lyα damping-wing spectra from a reionization simulation. In Section 3, we present our ML model for constraining the volume-averaged ionization fraction of the IGM, $Q_{HII}$, from input galactic spectra. In Section 4, we apply our ML model to existing JWST observation data to examine our methodology. In Section 5, we summarize and discuss our results.

## 2. Mock Lyα Damping-wing Spectra

### 2.1. Cosmic Dawn II Simulation

The Lyα opacity of the IGM is influenced by density and velocity fluctuations on subMpc scales, making numerical simulations essential for accurately modeling the complexity of individual sight lines. To this end, we utilize data from the Cosmic Dawn II (CoDa II) simulation to generate mock Lyα damping-wing (DW) spectra from galaxies during the reionization era. CoDa II is a cosmological simulation designed to replicate the formation of galaxies and the evolution of the IGM during cosmic reionization, in a cubic box of 94.4 Mpc or $64h^{-1}$ Mpc. The simulation computes the density, velocity, ionization fraction, and gas temperature on a 4096³ mesh, resolving structures with the cell size of 23.6 kpc. Dark matter is represented using 4096³ collisionless particles. Subgrid prescriptions are used to convert gas into stars, represented by collisionless star particles, and to estimate the ionizing radiation produced by these stars. These prescriptions are calibrated to reproduce the observed galaxy luminosity function at $z \sim 6$ (R. J. Bouwens et al. 2015, 2017; S. L. Finkelstein et al. 2015; H. Atek et al. 2018) and match existing constraints on reionization history, with reionization ending at $z = 6.2$. The simulation fully couples radiative transfer calculations with hydrodynamics, allowing for a realistic modeling of H II regions expanding from star-forming sites and their feedback on the IGM. The cosmology parameters used in the simulation are $\Omega_{m,0} = 0.307$, $\Omega_{b,0} = 0.048$, and $h = 0.678$, which are based on Planck Collaboration et al. (2014). Further details of the simulation can be found in P. Ocvirk et al. (2020).

CoDa II provides a robust data set for calculating Lyα opacity, capturing hundreds of H II bubbles from star-forming galaxies while resolving the IGM down to pressure-smoothing scales uniformly across the simulation volume. The Lyα opacity is sensitive to compact structures below the atomic-cooling mass ($\lesssim 10^8 M_{\odot}$; H. Park et al. 2021; A. Smith et al. 2022a; H. Park et al. 2024) and requires a spatial resolution element of 20 kpc or smaller for accurate convergence (Z. Lukić et al. 2015). The CoDa II simulation's cell size of 23.6 kpc is sufficiently small to resolve these small-scale structures. However, we note that the simulation box size of 94.4 Mpc is limited in sampling rare, bright galaxies with $M_{UV} \lesssim -22$ at $z \gtrsim 6$.

For this study, we use snapshots at $z = 6, 6.5, 7, 8$, and 10 to compute the mock spectra. The corresponding global ionization fractions in these snapshots are $Q_{HII} = 100\%, 80\%, 50\%, 13.2\%, \text{and } 1.6\%$, respectively. Although we use a limited number of snapshots for this proof-of-concept study, it is preferable to utilize more snapshots in practical applications to more accurately capture the redshift evolution of Lyα opacity. In future work, we plan to employ the CoDa III simulation data, which offers higher resolution at 8192³ and updated reionization prescriptions that better align with recent constraints,<fnref n="14" /> and to incorporate more than 20 snapshots to enhance the robustness of our results.

### 2.2. Lyα Opacity Calculation

To calculate the Lyα opacity by the IGM, we follow the general workflow outlined in H. Park et al. (2021) to conduct

<footnote n="14">¹⁴ The main paper for CoDa III is in preparation. However, some details of CoDa III can be found in J. S. W. Lewis et al. (2022) and P. Ocvirk et al. (2024).</footnote>
the calculation. The major difference in this work is that the calculation is done for a wider range of wavelengths (1216–1317 Å) to model the DW feature in low-resolution spectra, while H. Park et al. (2021) covered only from 1212–1220 Å to model Lyα emission lines.

The Lyα opacity is obtained by integrating the product of HI density and the Lyα cross section along the line of sight:

$$ \tau_{\alpha}(\nu_e) = \int_{r_{\min}}^{r_{\max}} n_{\text{HI}}(\mathbf{r}) \sigma_{\alpha}(T(\mathbf{r}), \nu(\mathbf{r})) a(z) dr, \quad (1) $$

where **r** is the proper position vector originating from the source galaxy, *a* is the cosmic scaling factor, $\nu_e$ is the frequency in the rest frame of the source, $n_{\text{HI}}$ is the neutral hydrogen density, and *T* is the gas temperature. Additionally,

$$ \nu(\mathbf{r}) = \nu_e \left[ 1 - \frac{\hat{\mathbf{r}} \cdot \mathbf{v}_{\text{pe}}(\mathbf{r}) + a(z)H(z)r}{c} \right] $$

accounts for the peculiar velocity of the IGM ($\mathbf{v}_{\text{pe}}$) and the cosmic redshift of photons ($aHr$) to calculate the Lyα cross section in the IGM frame. We integrate from the virial radius of the source galaxy (i.e., $r_{\min} = r_{\text{vir}}$) to the end of reionization ($r_{\max} = r(z=6)$) to calculate the IGM opacity.<fnref n="15" /> In the case of the $z=6$ snapshot, where reionization has ended in the simulation (i.e., $Q_{\text{HII}} = 100\%$), we simply set $\tau_{\alpha} = 0$ for all of the samples.

For other snapshots, we first integrate up to $r = 24 h^{-1} \text{Mpc}$ using the density, velocity, and temperature fields in the snapshot, assuming the fields remain static during the propagation of the radiation. This assumption is justified by the fact that the light travel time over that distance (e.g., 13.5 Myr at $z=7$) is much smaller than the Hubble time (e.g., ~800 Myr at $z=7$) or the duration of reionization (~500 Myr between $z=10$ and 6). We then evaluate the rest of the integral using the globally averaged H II fraction interpolated from averaged values across ~50 snapshots:

$$ \tau_{\alpha,\text{LS}} = \int_{24h^{-1}\text{Mpc}}^{r(z=6)} (1 - Q_{\text{HII}}(z)) \bar{n}_{\text{H}}(z) \sigma_{\alpha} a(z) dr. \quad (2) $$

This approach ignores fluctuations in $n_{\text{HI}}$ at distances beyond $24 h^{-1} \text{Mpc}$, which should not significantly affect the DW opacity. This large-scale opacity, $\tau_{\alpha,\text{LS}}$, contributes up to ~0.1 in optical depth (see H. Park et al. 2021 and A. Smith et al. 2022a for detailed analyses). For this term, we also account for the asymmetry in the DW cross section due to contributions from higher Lyman series damping wings (H.-W. Lee 2013):

$$ \sigma_{\alpha,\text{Lee}}(d\nu) = \sigma_{\alpha}(d\nu)(1 - 1.792 d\nu/\nu_{\alpha}), \quad (3) $$

where $\nu_{\alpha} = 2.47 \times 10^{15}$ Hz is the Lyα frequency.

We calculate the Lyα transmission, $T = \exp(-\tau_{\alpha})$, for 2000 sight lines from each of the 100 brightest galaxies of each snapshot, across 1000 equally spaced Lyα velocity offsets between $d\nu_{\alpha} = 1000$ and $25,000 \text{ km s}^{-1}$, corresponding to wavelengths between 1220 and 1317 Å. In total, we produce $2000 \times 100 \times 5 = 10^6$ transmission curves from the five snapshots.

A subset of the transmission curves is shown in Figure 1. The transmission shows a stronger dependence on $Q_{\text{HII}}$ below

<figure><img src="image_3.png" /><figcaption>**Figure 1.** A subset of Lyα transmission curves calculated from the CoDa II simulation data. Red, yellow, green, blue, and magenta curves are from the snapshots with the global ionization fraction of $Q_{\text{HII}} = 100\%$, 80\%, 50\%, 13.2\%, and 1.6\%, respectively.</figcaption></figure>

50%, with well-segregated transmission for $Q_{\text{HII}} = 50\%$, 13.2%, and 1.6% and no boundary cases—instances where the transmission values are near the overlap between different $Q_{\text{HII}}$ levels. Thus, these cases are easily classifiable, visually. However, the $Q_{\text{HII}}$ dependence becomes much weaker for higher global ionization cases (80% and 100%) as the transmission converges to 100%. Although these cases remain segregated from one another with certain boundaries, it is clear that noise from realistic observations (to be introduced below) would make the classification challenging. A similar trend can be seen from Figure 4 of H. Chen (2024), which shows results based on a different reionization simulation.

## 2.3. Mock Damping-wing Spectra

We process the Lyα transmission curves to generate realistic mock galaxy spectra near the Lyα resonance. LBGs typically exhibit near power-law rest-frame UV continua, where the normalization and slope can be accurately fitted given sufficient spectral coverage at wavelengths longer than the Lyα resonance (up to ~3000 Å). The transmission curve can then be recovered by normalizing the spectrum with the fitted continuum extrapolated to Lyα. Thus, the Lyα transmission curves calculated from the simulation effectively serve as normalized spectral energy distributions (SEDs) from observations.

Given that the continuum SEDs of high-redshift galaxies are generally very faint, low-resolution spectrographs, such as NIRSpec/PRISM, are typically used to measure the continuum with sufficient SNR for each spectral bin. To model such low-resolution spectra, we coarsen the transmission curves into eight bins, each with a width of 3000 km s⁻¹ in the rest frame, corresponding to the spectral bin size of such low-resolution instruments. We randomly select an example from the $Q_{\text{HII}} = 50\%$ cases to visually demonstrate how an SED is modified through the following processes. While this specific

<footnote n="15">¹⁵ We assume matter inside the virial radius belongs to the circumgalactic medium and interstellar medium.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Binned transmission curves (i.e., normalized SEDs) with added random Gaussian photon noise only (upper panels) and those with additional DLA opacity in addition to the photon noise (lower panels). The left, middle, and right panels show samples with photon noise of 1%, 5%, and 20% of the continuum flux, respectively. The black solid line is a randomly selected example from the $Q_{HII}$ = 50% samples without DLA opacity, and the three black dashed lines represent cases where a DLA opacity of $N_{HI}$ = $5 \times 10^{20}$, $1.5 \times 10^{21}$, and $1 \times 10^{22} \text{cm}^{-2}$ is added to the example case.</figcaption></figure>

case has a relatively higher flux than the average of the $Q_{HII}$ = 50% cases, it remains within the typical range, falling at the 24th percentile.

Next, we add random Gaussian noise to simulate different SNR levels. The standard deviation of the noise is set to 1%, 5%, and 20% of the continuum flux in each bin, corresponding to SNRs at continuum of 100, 20, and 5, respectively. A large fraction of JWST NIRSpec spectra are detector noise-limited.<fnref n="16" /> Thus, we assume SNR scales linearly with flux, with constant noise in the wavelength range where the flux is significantly reduced by the DW opacity. While a continuum SNR of 100 would be feasible only for rare bright objects, an SNR of 20 is achievable for moderately UV-bright sources with $M_{UV} \lesssim -22$ at $z \sim 7-8$ with a 1 hr integration of NIRSpec PRISM observations, and for exceptionally bright sources such as GN-z11, even at $z > 10$. The noise-added spectra are illustrated in the upper panels of Figure 2. With 1% noise (upper-left panel), it is relatively easy to visually determine that the example belongs to the $Q_{HII}$ = 50% samples. However, this classification becomes more challenging with 5% noise (upper-middle panel) and nearly impossible with 20% noise (upper-right panel).

In addition to photon noise, we simulate samples with additional DLA absorption. More than half of galactic

continuum spectra from $z \gtrsim 6$ appear to have extra absorption by a DLA with a column density of up to $\sim 10^{23} \text{cm}^{-2}$, in addition to the absorption by the IGM (T. Y.-Y. Hsiao et al. 2023; K. E. Heintz et al. 2025; H. Umeda et al. 2024). The DLA component is generally considered to be dense gas clumps in the ISM or CGM associated with the source galaxy. There have also been instances where a DLA appears to be associated with a remote absorber with a redshift offset of $\Delta z \sim 0.5$ from the source (C. Terp et al. 2024), though such cases appear to be rare. Since the spatial resolution of CoDa II is insufficient for resolving structures in the CGM and ISM, we opt to artificially add DLA opacity to the IGM transmission curves.

Based on recent measurements (e.g., K. E. Heintz et al. 2025), we simulate additional DLA components by duplicating the existing transmission curves and randomly assigning column densities between $N_{HI} = 2 \times 10^{20}$ and $1 \times 10^{23} \text{cm}^{-2}$, drawn uniformly. We assume the DLA is a pointlike absorber associated with the source galaxy. We apply the DLA opacity to the transmission before adding noise to avoid suppressing the noise, which should remain constant, by the DLA opacity.

The “DLA-added” samples are merged with the original “DLA-free” samples. A subset of the DLA-added mock spectra is shown in the lower panels of Figure 2. We also duplicate the example case into three additional cases, where weak, moderate

<footnote n="16">¹⁶ https://jwst-docs.stsci.edu/jwst-near-infrared-spectrograph/nirspec-observing-strategies/nirspec-background-recommended-strategies</footnote>
and strong DLA opacity with $N_{HI} = 5 \times 10^{20}$, $1.5 \times 10^{21}$, and $1 \times 10^{22} \text{ cm}^{-2}$ are added to sample spectra. The DLA-added cases present a much greater challenge in determining $Q_{HII}$ at all noise levels, as the attenuation by DLA opacity removes the boundary between different $Q_{HII}$ values on the transmission.

## 2.4. IGM versus DLA Opacity in the Data Space

Then, how do we deal with the contamination by DLA opacity? The difference in wavelength dependence can help us distinguish the DLA opacity from the IGM opacity. The $Ly\alpha$ cross section in DW scales as the inverse square of the wavelength offset. Since the DLA opacity comes from a pointlike absorber, it scales in the same way as the cross section itself does (i.e., $\tau_{DLA} \propto dv_\alpha^{-2}$) while the large-scale IGM opacity comes from the optical depth integrated from the edge of an H II region to a semi-infinite distance (several hundred Mpc to the end of EoR), thereby approximately scaling as the inverse of the offset (i.e., $\tau_{IGM} \propto dv_\alpha^{-1}$).<fnref n="17" />

Motivated by this fact, we make a scatter plot of transmissions at two different $Ly\alpha$ velocity offsets, $dv_\alpha = 3500$ and $5000 \text{ km s}^{-1}$, corresponding to 1226 and 1238 Å in wavelength in Figure 3 to distinguish DLA-added and DLA-free spectra based on the wavelength dependence of the opacity. Those two quantities plotted, $T(1226 \text{ Å})$ and $T(1238 \text{ Å})$, correspond to the two data points with the shortest wavelengths in Figure 2. The example cases from Figure 2 are shown as black symbols to show their location in the data space.

The 1% noise case shown in Figure 3(a) clearly reveals the distinct distribution of DLA-added versus DLA-free cases in the data space. As marked by two the dashed arrows, the DLA-added cases have a relatively higher flux at 1238 Å when the flux at 1226 Å is similar due to the steeper decline of DLA opacity with increasing wavelength. Also, the DLA-added cases have a much wider range of flux inherited from the large dynamic range of the DLA opacity.

The black square in Figure 3(a) is a $Q_{HII} = 50\%$ case with a DLA opacity of $N_{HI} = 1.5 \times 10^{21} \text{ cm}^{-2}$, and it exemplifies the difference between DLA-added versus DLA-free groups. It has $T(1226 \text{ Å}) = 0.5$, which falls within the range occupied by DLA-free samples with $Q_{HII} = 1.6\%$ (magenta dots), $0.45 \lesssim T(1226 \text{ Å}) \lesssim 0.65$. However, its flux in the second bin, $T(1238 \text{ Å}) = 0.86$, is outside the range of that group, $0.75 \lesssim T(1238 \text{ Å}) \lesssim 0.82$, thereby allowing us to determine that it does not belong to the DLA-free $Q_{HII} = 1.6\%$ cases.

However, the difference in the wavelength dependence between DLA and IGM opacity does not appear large enough to allow for an easy determination of $Q_{HII}$ for DLA-added samples, which have both opacity components. Unlike the DLA-free sample, which forms compact and distinct clusters by $Q_{HII}$ in the data space, DLA-added samples form extended trails with much shorter distances between different $Q_{HII}$ values. Moreover, the trails converge to a single line with increasing DLA opacity, as $T(1226 \text{ Å})$ is attenuated to near zero, exemplified by the black triangle symbol in Figure 3(a), a case with $Q_{HII} = 50\%$ and $N_{HI} \gtrsim 10^{22} \text{ cm}^{-2}$. Obviously, it would be nearly impossible to determine $Q_{HII}$ for these strong-DLA samples.

Increasing the noise level increases scatter on the data space, mixing samples with different labels. In the 5% noise case shown in Figure 3(b), DLA-added samples are mostly mixed with DLA-free cases at $T(1226 \text{ Å}) > 0.5$. The example DLA-added cases with $N_{HI} = 5 \times 10^{20}$ (square) and $1.5 \times 10^{21} \text{ cm}^{-2}$ (diamond) are now visually indistinguishable from DLA-free cases with $Q_{HII} = 13.2$ (blue dots) and $1.6\%$ (magenta dots), respectively. Also, DLA-free cases with $Q_{HII} = 50\%, 80\%,$ and $100\%$ are largely mixed with each other due to the small difference in the IGM opacity. When the noise level is increased to 20% (Figure 3(c)), the grouping of samples by $Q_{HII}$ is mostly disrupted due to the increased scatter, severely hindering visual classification.

We note that other combinations of spectral bins at longer wavelengths also exhibit similar trends but are more vulnerable to noise as the $Ly\alpha$ opacity is weaker for those bins. Thus, the first two spectral bins are most useful for determining $Q_{HII}$ and DLA presence.

# 3. ML Model for Constraining $Q_{HII}$

## 3.1. Random Forest Classifier

We use the Random Forest Classifier (hereafter RFC) to create an ML model that constrains $Q_{HII}$ from input galactic rest-frame UV spectra with DW feature. RFC is a supervised ML algorithm widely used for classification and regression tasks involving labeled data. It is an ensemble learning algorithm, which combines the predictions of multiple base estimators to improve predictive accuracy and reduce the risk of overfitting (L. Breiman 2001).

Each base estimator in RFC is a decision tree (see Figure 4 for an example). A decision tree is a flowchart-like structure where each internal node represents a test on a feature, each branch represents the outcome of the test, and each leaf node represents a class label (in classification) or a continuous value (in regression). The paths from the root to the leaf represent classification rules.

An individual decision tree is limited in modeling complex boundaries in data sets when used alone. By aggregating the results of many such decision trees, the random forest algorithm significantly enhances generalizability and robustness. This ensemble approach leverages the diversity of the individual trees to make more accurate and stable predictions.

## 3.2. Model Training

In our study, the input data for our ML model comprises $10^6$ mock $Ly\alpha$ DW spectra, each consisting of eight flux data bins, with labels corresponding to the $Q_{HII}$ values. We utilize the “RandomForestClassifier” class from the “sklearn.ensemble” module in the SciPy package to train an RFC model. Although the “RandomForestRegressor” class would typically be more appropriate for predicting continuous $Q_{HII}$ values, we opted for the classifier due to the discrete nature of our labels, limited to five distinct $Q_{HII}$ values in this work, and for easier visualization of the results. “RandomForestRegressor” essentially uses an ensemble of decision trees similar to RFCs for regression tasks. Thus, the underlying principles of our RFC model remain fundamentally aligned with those of the regressor.

As well as the RFC for determining $Q_{HII}$ (Q-RFC, hereafter), we also train another RFC to determine whether the spectrum has a DLA opacity component (DLA-RFC, hereafter). The DLA-RFC is trained on the same training data, but with another

<footnote n="17">¹⁷ In reality, the IGM opacity would scale differently from $dv_\alpha^{-1}$ because the mean H I density falls as reionization progresses and the Universe expands.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Scatter plot of $\mathcal{T}(1238\AA)$ vs. $\mathcal{T}(1226\AA)$ or transmission in the first vs. second shortest wavelength bins for the mock spectra shown in Figure 2. The circular dots are for the DLA-free cases and the minus symbols are for the DLA-added cases. Red, yellow, green, blue, and magenta colors correspond to $Q_{HII} = 100\%, 80\%, 50\%, 13.2\%,$ and $1.6\%$, respectively. The star symbol denotes the example DLA-free case from Figure 2, and the diamond, square, and triangle symbols denote the DLA-added versions of the example case with $N_{HII} = 5 \times 10^{20}$, $1.5 \times 10^{21}$, and $1 \times 10^{22} \text{ cm}^{-2}$, respectively. Panels (a), (b), and (c) are the results for the flux noise levels of 1%, 5%, and 20%, respectively.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Example decision tree for predicting $Q_{HII}$ based on Ly$\alpha$ transmission at $\lambda = 1226$ Å. This is a simplified case using data from only one wavelength, while a practical model would utilize all eight wavelengths. The tree starts at depth 1 with the first decision node, branching into subsequent nodes based on transmission thresholds. Each decision is made by evaluating whether $T(1226\AA)$ exceeds a specific threshold.</figcaption></figure>

set of binary labels containing whether the spectrum was added an additional DLA opacity. As demonstrated in Figure 3, DLA-added cases appear more difficult to determine their $Q_{HII}$. Thus, we aim for improving the classification accuracy by sorting out DLA-added cases using this DLA-RFC.

The Q-RFC can be trained both with and without the DLA-added mock spectra. We use the Q-RFC trained with the DLA-added spectra in Section 3.3 to estimate the Q-RFC accuracy for both the DLA-added and DLA-free sample. In Section 4, we use the Q-RFC trained with the DLA-free spectra only, as we will not apply our model to DLA-added samples based on the analysis from Section 3.3 except for one exceptionally high-SNR case analyzed in Section 4.4.

In addition to our main eight-parameter RFC models, we also generate two-parameter models using transmission values for the two shortest wavelengths, $T(1238\AA)$ and $T(1226\AA)$, to visualize the decision boundaries in the transmission-transmission plots of Figure 3. Although the two-parameter models yield less-accurate classification results and are not used in the main analysis, the decision boundaries visualized on a 2D plane for different noise levels provide useful insights into how the model is affected by noise. The decision boundary from this 2D data space offers a useful perspective on the full 8D data space because these two bins provide the most constraining information in practice due to the much weaker Ly$\alpha$ opacity in the other six bins with longer wavelengths. These six bins are more susceptible to noise and provide weaker constraints, thus only offering a “minority report” that slightly improves the model’s accuracy.

The number of data points per sample defines the dimensionality of the data space on which the RFC will train. For continuous quantities such as SED, it is often advisable to reduce data dimensionality using methods like principal component analysis (PCA) to avoid the complexities of training ML models in high-dimensional spaces. However, in our case, the mock spectra are already binned into just eight wavelength bins, making further compression unnecessary. We experimented with PCA to reduce the data dimensions even further, but this did not enhance the classification results.

Therefore, we use the Ly$\alpha$ DW spectra in their original binned form without any additional processing.

For the hyperparameters of RFC, we set the number of base estimators to 30 (i.e., “n_estimators = 30”) and the maximum depth of individual estimator to 20 (i.e., “max_depth = 20”). Setting these parameters too high can result in overfitted models, which create unnecessarily complex decision boundaries and increase computational time without meaningful gains in classification accuracy. This choice of hyperparameters leads to overfitted models in some cases, giving a significantly higher accuracy for training data set than for the test set, but the computation time for querying the RFCs stays below a millisecond and does not consume a significant computing time. Thus, we apply the above-mentioned parameter choice ($n_estimators = 30$ & $max_depth = 20$) across all of the cases for consistency.

The accuracy of the model is defined by the fraction of correctly classified cases from the test data set. Owing to the large number of our mock spectra ($10^6$), the accuracy of the models converges after training with only 5% of the entire mock samples ($5 \times 10^4$). Thus, we use the remaining 95% as the test set to calculate the accuracy of the RFCs.

As will be shown below, a mock spectrum with a weak DLA ($N_{HII} \lesssim 10^{21} \text{cm}^{-2}$) is much more prone to misidentification to a DLA-free case with a lower $Q_{HII}$ value than its correct $Q_{HII}$ label than stronger DLA cases. This implies that both the DLA-RFC and Q-RFC accuracies are subject to the assumed probability distribution of the column density of added DLA, and the accuracy calculated in this work should not be considered an absolute value.

### 3.3. Classification Results for the Test Samples

In Figure 5, we visualize the decision boundaries of our 2D RFC models for the three noise levels (1%, 5%, and 20%) as colored backgrounds overlaid on the transmission-transmission plots from Figure 3. The decision boundaries are generated by querying the trained 2D Q-RFC and DLA-RFC models on a $1000 \times 1000$ grid in the 2D data space.<fnref n="18" /> The hue of the background represents the decision boundary for $Q_{HII}$ based on the Q-RFC model, while the brightness indicates the DLA-RFC results, with brighter/darker regions corresponding to DLA-added/DLA-free classifications, respectively.

Each panel displays the classification accuracy values for DLA-RFC and Q-RFC, derived from the primary 8D models, in the upper-left corner. The Q-RFC accuracy is shown separately for DLA-free and DLA-added groups, where we use the DLA-RFC result—rather than the true DLA presence label—to separate these groups. This approach provides a conservative estimate of the Q-RFC accuracy for DLA-free samples, assuming that the true DLA presence is unknown and must be determined by the DLA-RFC. Consequently, any misclassified DLA-free samples would lower the classification accuracy compared to a scenario where the true DLA presence is known, as Q-RFC tends to be less accurate for DLA-added samples.

Figure 6 shows the Q-RFC results for $5 \times 5$ combinations of input and predicted $Q_{HII}$ values as colored histograms, separately for the DLA-free and DLA-added groups. The x-labels represent the predicted $Q_{HII}$ values, and the color of the bar indicates the true $Q_{HII}$ value of the input spectra. The bar

<footnote n="18">¹⁸ As noted above, these 2D models and decision boundaries are used for visualization purposes only.</footnote>
(a) Noise level = 1 %

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

(b) Noise level = 5 %

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

(c) Noise level = 20 %

<figure><img src="image_9.png" /><figcaption>**Figure 5.** The 2D RFC results for the first and second spectral bins overlaid to the Figure 3 as colored background. The hue of the background visualizes the Q-RFC results, where red, yellow, green, blue, and magenta corresponds to the same $Q_{HII}$ values as the data points do: 100%, 80%, 50%, 13.2%, and 1.6%, respectively. The brightness of the colors visualizes the DLA-RFC results, where the brighter/darker background marks the region where the samples are classified DLA-added/DLA-free.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 6. Histogram of Q-RFC prediction results by the $Q_{HII}$ values of the input sample, 100%, 80%, 50%, 13.2%, and 1.6% shown as red, yellow, green, blue, and magenta histograms, respectively. The predicted $Q_{HII}$ value on the x-axis is colored in the same way so that the histogram matching the x-label color is given the fraction of corrected classified for sample (i.e., Q-RFC accuracy) for the corresponding $Q_{HII}$ value. The DLA-free (upper panels) and DLA-added cases (lower panels) are shown for three different noise levels, 1%, 5%, and 20%, in the left, middle, and right panels, respectively. Here, the DLA-free and DLA-added samples are based on the DLA-RFC results, not on the true labels.</figcaption></figure>

whose color matches the x-label color shows the correctly classified fraction (i.e., classification accuracy) of the Q-RFC for the corresponding $Q_{HII}$ label.

### 3.3.1. One Percent Noise; SNR = 100

The RFC results for 1% noise or continuum SNR of 100 are shown in Figure 5(a). The Q-RFC correctly classifies 90% of the DLA-free samples, but is less accurate (64%) for the DLA-added samples. The DLA-RFC performs with high accuracy, achieving a 97% success rate in distinguishing DLA-added from DLA-free samples.

The Q-RFC decision boundary effectively separates different $Q_{HII}$ groups of the DLA-free samples, owing to the strong clustering of these groups in the data space. It is also effective for the DLA-added samples with $Q_{HII} = 1.6\%$, 13.2%, and 50% samples with $\mathcal{T}(1226 \AA) > 0.2$, as indicated by the bright magenta, blue, and green regions that distinguish these samples. The DLA-free example case (star symbol) and the cases with moderate ($1.5 \times 10^{21} \text{cm}^{-2}$; square symbol) and weak ($5 \times 10^{20} \text{cm}^{-2}$; diamond symbols) DLAs are all correctly classified as having $Q_{HII} = 50\%$.

Misclassification occurs mostly in regions with fragmented decision boundaries, where the RFC struggles to find clear separations between labels and overfits to the training data. The boundary between $Q_{HII} = 100\%$ (red) and 80% (yellow) is highly

fragmented for both DLA-free and DLA-added samples, as the 1% noise is comparable to the opacity difference between these $Q_{HII}$ values. The decision boundary is also poorly defined at $\mathcal{T}(1226 \AA) < 0.2$ and $\mathcal{T}(1238 \AA) < 0.6$ (lower-left corner in the figure), a region populated by strong DLA samples. As pointed out in Section 2.4, classification is highly ineffective for these samples, across all noise levels, as the IGM and DLA opacity become highly degenerate, causing the samples to converge to the same region in the data space. The example case with a strong DLA ($N_{HI} = 10^{22} \text{cm}^{-2}$; triangle symbol) is located in this region and is incorrectly classified as having $Q_{HII} = 100\%$.

The DLA-RFC decision boundary between the dark and bright background is well defined, reflected in the 97% classification accuracy. Misclassifications involve weak-DLA cases with $N_{HI} \sim 2 \times 10^{20} \text{cm}^{-2}$. The DLA-added examples with $N_{HI} = 5 \times 10^{20}$, $1.5 \times 10^{21}$, and $1 \times 10^{22} \text{cm}^{-2}$ are all correctly classified at this noise levels.

The Q-RFC prediction statistics in Figure 6 also show that the Q-RFC performs better for DLA-free samples (Figure 6(a)) compared to DLA-added sample (Figure 6(d)). Misclassification primarily occurs between adjacent $Q_{HII}$ values with unclear decision boundaries. In the DLA-free samples (Figure 6(a)), most misclassifications involve confusing around 20% of $Q_{HII} = 80\%$ and 100% samples with each other, with about 3% of the $Q_{HII} = 50\%$ samples misclassified as $Q_{HII} = 80\%$.
Conversely, the $Q_{HII}$ = 13.2% and 1.6% samples are nearly perfectly classified. While the accuracy decreases for DLA-added samples, it remains relatively high for $Q_{HII}$ = 13.2% and 1.6%, indicating that low-resolution DW spectra are more effective at constraining reionization history during the early stages ($Q_{HII}$ ≲ 50%) than during the later stages.

These results demonstrate the impressive performance of RFC models in constraining $Q_{HII}$ using LBG spectra observed with SNRs of 100 or higher. However, achieving such high SNRs from reionization-era galaxies is feasible only for rare bright objects, making it more practical to consider scenarios with higher noise levels.

### 3.3.2. Five Percent Noise; SNR = 20

An SNR of 20 is achievable for moderately bright galaxies with $M_{UV}$ ≲ −22 at $z$ ~ 7 with ~1 hr exposure by NIRSpec PRISM observation. When the noise level is increased to 5% (see Figure 5(b)), the increased mixing between different labels results in significantly more fragmented decision boundaries. The Q-RFC accuracy drops to 67% for DLA-free samples and 39% for DLA-added samples. Given that a complete random classification would yield a 20% accuracy, the Q-RFC is particularly ineffective for DLA-added samples at this noise level. Meanwhile, the DLA-RFC maintains relatively high accuracy at 90.3%.

Even with the increased noise, the Q-RFC still identifies well-defined boundaries between $Q_{HII}$ = 1.6% (dark magenta) and 13.2% (dark blue), and between $Q_{HII}$ = 13.2% and 50% (dark green) for DLA-free samples. However, for DLA-added samples across all $Q_{HII}$ values and DLA-free samples with $Q_{HII}$ ≳ 50%, the boundaries are poorly defined due to noise surpassing the opacity differences. The DLA-free example case (star symbol) with a true $Q_{HII}$ value of 50% is incorrectly classified as $Q_{HII}$ = 80%. Moreover, the DLA-added examples with $N_{HI}$ = 5 × 10²⁰ (diamond symbol), 1.5 × 1 × 10²¹ (square symbol), and 10²² cm⁻² (triangle symbol) are misclassified as $Q_{HII}$ = 13.2%, 1.6%, and 1.6%, respectively.

In case of the DLA-RFC, the samples with $T(1226\AA) < 0.5$ are classified as DLA-added with nearly 100% accuracy. The triangle symbol (a case with $N_{HI}$ = 10²² cm⁻²) falls into this category and is correctly identified as DLA-added with 100% confidence. However, DLA-added samples with $T(1226\AA) > 0.5$ are largely mixed with DLA-free samples, as the noise far exceeds the difference in wavelength dependence between DLA and IGM opacity. Many DLA-added samples with $T(1226\AA) > 0.5$ are misclassified as DLA-free cases with a lower $Q_{HII}$ than its true value. The $Q_{HII}$ = 50% examples with a weak ($N_{HI}$ = 5 × 10²⁰; diamond symbol) and a moderate DLA ($N_{HI}$ = 1.5 × 10²¹; square symbol), which are misclassified as DLA-free cases with $Q_{HII}$ = 13.2%, fall into this category. The fraction of such confused cases is only 10% in this case, but it could be higher if the fraction of weak DLAs is higher than assumed in this work.

For the DLA-RFC, samples with $T(1226\AA) < 0.5$ are classified as DLA-added with nearly 100% accuracy, including the triangle symbol (representing $N_{HI}$ = 10²² cm⁻²), which is correctly identified with 100% confidence. However, DLA-added samples with $T(1226\AA) > 0.5$ are frequently confused with DLA-free samples due to the noise overwhelming the differences in the wavelength dependence of DLA and IGM opacity, and are often misclassified as DLA-free cases with a lower $Q_{HII}$ than their true values. For example, the $Q_{HII}$ = 50% cases with a weak ($N_{HI}$ = 5 × 10²⁰; diamond symbol) and moderate DLA

($N_{HI}$ = 1.5 × 10²¹; square symbol) are misclassified as DLA-free cases with $Q_{HII}$ = 13.2%. Although only 10% of the DLA-added samples fall into this category, this fraction could be higher if weak DLAs are more prevalent than assumed.

The Q-RFC prediction statistics in Figure 6(e) highlight that DLA-added samples are much less effective in constraining $Q_{HII}$, except for $Q_{HII}$ = 1.6%, which still achieves 50% accuracy. Samples with other $Q_{HII}$ values slightly exceed 20% accuracy, with $Q_{HII}$ = 13.2% reaching 40%. This is because the highly fragmented decision boundary still roughly captures the statistical differences in the data distribution. However, given the low accuracy and the lack of clear visual confirmation, we do not consider the Q-RFC effective in this regime. Figure 6(b) indicates that the Q-RFC can reach meaningful accuracy if DLA-added cases are filtered using the DLA-RFC, with accuracy particularly high for $Q_{HII}$ = 13.2% and 1.6%, exceeding 85%. The accuracy is slightly lower for $Q_{HII}$ = 50%, 80%, and 100% cases, but still significantly higher than the random baseline of 20%.

### 3.3.3. Twenty Percent Noise; SNR = 5

In the 20% noise case, shown in Figure 5(c), the Q-RFC accuracy degrades to 34% for DLA-free samples and 30% for DLA-added samples (see Figure 5(c)). The highly fragmented decision boundaries between different colors illustrate the ineffectiveness of the model. The Q-RFC prediction statistics in Figures 6(c) and (f) also demonstrate that the Q-RFC is unable to correctly determine $Q_{HII}$ regardless of DLA presence.

The DLA-RFC establishes an approximate boundary between DLA-free and DLA-added samples at $T(1226\AA)$ ~ 0.5 (indicated by the bright versus dark background in Figure 5(c)) and correctly determine DLA presence for more than 80% of the samples. Figure 7 shows that most of the correctly classified DLA-added cases (i.e., true positives) are distributed in the low-transmission region in the data space. Given that DLA-RFC performs a binary classification between DLA-free versus DLA-added, its 80% accuracy in this noise level is modest at best. However, it can still serve as a screening tool for higher-SNR observations.

## 4. Application to Existing Observations

We train and apply our RFC models to existing high-redshift spectra to assess their effectiveness in constraining $Q_{HII}$ at the redshift of the source galaxy, thereby providing insights into reionization models. To ensure the RFC models are trained accurately for the specific conditions of each spectrum, we generate mock Ly$\alpha$ DW spectra tailored to each input spectrum. This involves matching the wavelength range of the spectral bins, the SNR of the flux data in each bin, and the spectral resolution of the instrument.

First, we smooth the Ly$\alpha$ transmission curves from the CoDa II simulation using a Gaussian line-spread function,

$$f_{\text{LSF}}(dv_{\alpha}) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{dv_{\alpha}^2}{2\sigma}\right), \quad (4)$$

where we set $\sigma = c/(2.355R)$ so that the FWHM is given by the speed of light, c, divided by the spectral resolution,¹⁹R.

<footnote n="19">¹⁹ The spectral resolution of PRISM can be obtained from https://jwst-docs.stsci.edu/jwst-near-infrared-spectrograph/nirspec-instrumentation/nirspec-dispersers-and-filters.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Histogram of total DLA samples (black) compared to the correctly classified DLA samples by our DLA-RFC model (blue) with 20% noise as a function of $T(1226\AA)$, the transmission in the shortest wavelength bin. The gray histogram represents the limit of a completely random 50–50 selection.</figcaption></figure>

Next, we rebin the smoothed transmission curves to match the wavelength range of the spectral bins and add Gaussian noise based on the SNR of each bin. We use spectral bins whose entire extent falls within the range of $dv_\alpha = 1000–25,000\,km\,s^{-1}$ to sample the wavelength range where DW opacity is significant, while avoiding modeling uncertainties related to the velocity offset of the Lyman break. We then train both the DLA-RFC and Q-RFC models using these tailored mock spectra to determine the presence of a DLA and the value of $Q_{HII}$ for each spectrum.

The continuum in the galaxy spectrum is fitted using a power-law between 1350 and 2600 Å, and the intrinsic UV continuum level is extrapolated to 1215.67 Å to estimate the DW opacity of the IGM. Additionally, we subtract any Ly$\alpha$ emission-line leakage into other spectral bins caused by the low-resolution spectroscopy, provided the exact location and strength of the line are known from higher-resolution spectroscopy. As with the smoothing of the transmission curves, we apply the line-spread function from Equation (4) to estimate the leakage across spectral bins.

We inspect publicly available spectra from the JWST Advanced Deep Extragalactic Survey<fnref n="20" /> (JADES, GTO 1180, GTO 1181, PI: D. Eisenstein, GTO 1210, PI: N. Lützgendorf, GO 3215, PI: D. Eisenstein & R. Maiolino; D. J. Eisenstein et al. 2023; A. J. Bunker et al. 2024; F. D'Eugenio et al. 2025) and the Cosmic Evolution Early Release Science program<fnref n="21" /> (ERS 1345, PI: S. Finkelstein; S. L. Finkelstein et al. 2024) to identify suitable spectra for applying our model. We select galaxies with continuum SNRs of 10 or above, along with measured Ly$\alpha$ line offsets and EW. As demonstrated below, accurate Ly$\alpha$ line measurement is crucial for our analysis, as leakage from the Ly$\alpha$ line can significantly contaminate flux in other spectral bins. We then use our DLA-RFC to select three galaxies that are likely DLA-free (GN-z11, CEERS-1026, and JADES-12637) for analysis. Additionally, we select JADES-58975 for analysis despite the presence of apparent DLA

opacity, due to its high continuum SNR of 40. The Jupyter notebook detailing these procedures is publicly available on GitHub.<fnref n="22" />

### 4.1. GN-z11

GN-z11 is an extremely bright galaxy with $M_{UV} = -21.5$, spectroscopically confirmed at $z > 10$ by JWST. It exhibits a clear DW feature in its NIRSpec PRISM spectrum (A. J. Bunker et al. 2023). This galaxy hosts an AGN with a broad Ly$\alpha$ line extended to $dv_\alpha \approx 1000\,km\,s^{-1}$ (R. Maiolino et al. 2024b, 2024a). The emission line is centered at $dv_\alpha = 555\,km\,s^{-1}$ with EW of 18 Å. Two publicly available spectra of GN-z11 from JADES were combined by averaging them, weighting each spectrum by the inverse of its SNR. The SNR of the combined spectrum was calculated using the Pythagorean summation of the SNRs from each observation. We use eight spectral bins at $dv_\alpha = 3400, 5900, 8500, 11,100, 13,700, 16,200, 18,800$, and $21,400\,km\,s^{-1}$ for the analysis.

Figure 8 illustrates the workflow used to estimate $Q_{HII}$ at the redshift of GN-z11 ($z = 10.6$). Panel (a) shows the raw SED near the rest-frame Ly$\alpha$ and the power-law fit to the continuum. Panel (b) displays the eight spectral bins near Ly$\alpha$ used to constrain $Q_{HII}$. We subtract the leakage from the Ly$\alpha$ emission line based on the spectral resolution of PRISM around the rest-frame Ly$\alpha$ ($R = 32$). The emission-line subtraction results in a significant change in flux, with the normalized flux being altered by approximately 0.35 in the shortest wavelength bin with $dv_\alpha \approx 4000\,km\,s^{-1}$. Panel (c) shows that the combined SNR of the two GN-z11 spectra is close to 20, comparable to the 5% noise cases in Section 3.

Panel (d) presents the mock spectra for different $Q_{HII}$ values alongside the GN-z11 spectrum. At the first spectral bin ($dv_\alpha = 4000\,km\,s^{-1}$ or $\lambda_{rest} = 1229\AA$), the flux of GN-z11 is slightly below $Q_{HII} = 1.6\%$ samples after the emission-line subtraction (0.44) while it is above $Q_{HII} = 13.2\%$ samples before subtraction (0.78). It is evident that this flux change from the emission-line subtraction would significantly impact the classification result for GN-z11.

Panel (e) shows the mock and GN-z11 spectra in the data space of the two shortest wavelength bins at 1229 and 1239 Å, with the results of the 2D Q-RFC and DLA-RFC models represented by the colored background, similar to in Figure 5. As described in Section 3.3, all of the data points are used for the final classification. GN-z11 is classified as having $Q_{HII} = 1.6\%$ after the emission-line subtraction (star symbol) and $Q_{HII} = 50\%$ before subtraction. The data point before emission-line subtraction is highly offset from the mock data points, indicating that it is invalid for classification without subtraction. Regarding DLA presence, GN-z11 is slightly outside the decision boundary for DLA-free samples due to its slightly lower flux compared to typical $Q_{HII} = 1.6\%$ cases. The DLA-RFC model determines GN-z11 has a 60% probability of containing a DLA. However, its location in the data space is still populated by other DLA-free $Q_{HII} = 1.6\%$ samples, suggesting that GN-z11 could still be DLA-free. Given that prediction probability depends on the DLA fraction and the assumed DLA column density distribution in the mock samples, it remains inconclusive whether GN-z11 indeed contains a DLA.

Panel (f) illustrates the “prediction probability” for each $Q_{HII}$ value, based on the number of votes from individual decision trees in the Q-RFC for the emission-line subtracted spectrum.

<footnote n="20">²⁰ https://jades-survey.github.io</footnote>

<footnote n="21">²¹ https://ceers.github.io</footnote>

<footnote n="22">²² https://github.com/hcosmosb/RFC4LBGDW</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 8.** This figure summarizes the workflow of constraining $Q_{HII}$ at $z = 10.6$ with the spectrum of GN-z11. Panel (a): the raw spectral energy distribution (SED) at the rest frame is shown as the black line, and the power-law curve, which best fits the data between 1350 and 2600 Å, is shown as the blue line. The fitting range is enclosed by two yellow circular dots. Panel (b): Lyα transmission as the function of $dv_α$ between 0 and 25,000 km s$^{-1}$, obtained by normalizing the SED with the fitting function. The dashed and solid lines represent the normalized flux before and after correcting for the leakage from the emission line, respectively. Panel (c): signal-to-noise ratio (SNR) of each data point in panel (b). Panel (d): mock Lyα spectra from CoDa II generated according to the flux range of the data bins and SNR of the observed spectrum. The spectrum of GN-z11 is also shown as a black line. Panel (e): transmission-transmission plot of the two shortest wavelength bins, with the results of the 2D Q-RFC and DLA-RFC shown as colored backgrounds as in Figure 5. The GN-z11 spectrum after subtracting the emission-line contribution is shown as the black star symbol, and that before the subtraction is shown as the black square symbol. Panel (f): prediction probability of $Q_{HII}$ values for $z = 10.6$ from our Q-RFC model.</figcaption></figure>

We can regard the prediction probability for each label as the likelihood for $Q_{HII}$. The Q-RFC strongly favors $Q_{HII} = 1.6\%$ and 13.2% over other $Q_{HII}$ values, with prediction probabilities of 54% and 30%, respectively. Thus, we conclude $Q_{HII} < 31.6\%$ at $z=10.6$ with 85% confidence, where $31.6\% = [(13.2 + 50)/2]\%$ is the midpoint between the $Q_{HII} = 13.2$ and 50% bins. This result is unlikely to constrain reionization models effectively, as most predict the Universe was less ionized than 31.6% at this redshift. Instead, this outcome serves as a test confirming the validity of our RFC method.

## 4.2. CEERS-1029

CEERS-1029 is one of the targets of the CEERS program, with $M_{UV} = -21.53$ at $z = 8.61$, which was first discovered by R. L. Larson et al. (2022). The EW and velocity offset of the Lyα emission are measured to be 3 Å and $dv_α = 2100$ km s$^{-1}$, respectively (M. Tang et al. 2024). The PRISM spectrum has an SNR of approximately 15 in the continuum. We use four spectral bins at $dv_α = 5400, 9800, 14,300$, and $18,800$ km s$^{-1}$ for the analysis.

Figure 9 shows the raw flux data with the power-law fit, the normalized spectrum alongside the mock spectra for the target, a scatter plot of the two shortest wavelength bins, and the prediction probability of Q-RFC for the emission-line subtracted spectrum. The normalized flux in panel (d) show that subtracting the emission leakage results in only a minor change (few percent) in the shortest wavelength bin, due to the small EW of the Lyα emission.

In the data space shown in panel (e), CEERS-1029 falls in the region for DLA-free $Q_{HII} = 1.6\%$ cases, regardless of the emission-line subtraction. The Q-RFC model predicts a 90% probability for $Q_{HII} = 1.6\%$ and much lower probability for other $Q_{HII}$ values (panel (f)). We interpret this result as an upper bound of $Q_{HII} < 7.4\%$ at $z = 8.6$. This upper bound is quite tight for $z = 8.6$ and would rule out some of the early reionization models.

However, we note that the spectrum of CEERS-1029 may not have been precisely fitted by a power law. The raw flux in panel (a) shows a long-range fluctuation around the fit, suggesting there may be another component in the continuum not captured by the power law. The normalized spectrum of CEERS-1029 near Lyα in panel (d) appears flatter than those of the mock $Q_{HII} = 1.6\%$ cases, with the flux at the two longer
<figure><img src="image_13.png" /><figcaption>Figure 9. Same as panels (a), (d), (e), and (f) of Figure 8 for CEERS-1029 at z = 8.6.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 10. Same as panels (a), (d), (e), and (f) of Figure 8 for JADES-12637 at z = 7.7.</figcaption></figure>

wavelength bins being lower than the typical range of the mock spectra for $Q_{HII}$ = 1.6%. If the actual continuum is lower by, for example, 10%, than what the power-law fitting gave, the Q-RFC would have given significantly higher probabilities for $Q_{HII}$ = 13.2% or 50%.

### 4.3. JADES-12637

JADES-12637 is one of the targets of the JADES program, with $M_{UV} = -20.7$ at $z=7.7$. The Ly$\alpha$ equivalent width and $dv_\alpha$ of the Ly$\alpha$ emission are measured to be 33 Å and 277 km s$^{-1}$, respectively (M. Tang et al. 2024). The PRISM spectrum has an SNR of approximately 10 in the continuum. We use eight spectral bins at $dv_\alpha$ = 3300, 6000, 8700, 11,500, 14,200, 17,100, 19,900, and 22,800 km s$^{-1}$ for the analysis.

Figure 10 presents the raw flux with the power-law fit, the normalized spectrum of JADES-12637, and its corresponding mock spectra, a scatter plot of these spectra in the 2D data space of the two shortest wavelength bins, and the prediction probability for the Q-RFC model. The power-law fit in panel (a) appears to fit the data well. The normalized flux in panel (d) shows that subtracting the emission-line leakage, with EW = 33 Å, results in a significant change in the spectrum, nearly halving the flux in the first bin.

In the data space (panel (e)), JADES-12637 is positioned near the center of the DLA-free $Q_{HII}$ = 1.6% group after the emission-line subtraction (indicated by the black star symbol). Without the subtraction (black square symbol), it falls entirely outside all of the mock data points. In panel (f), the Q-RFC weakly prefers $Q_{HII}$ = 13.2% with a prediction probability of 35%, with probabilities for other $Q_{HII}$ values higher than 14%. The apparent inconsistency with panel (e), appearing to prefer $Q_{HII}$ = 1.6%, is due to the five data points at longer wavelengths, not described in panel (e), preferring higher $Q_{HII}$ due to their relatively higher flux. Combining probability

for the three lower $Q_{HII}$ values gives an upper bound of $Q_{HII} < 65\%$ with 71.2% confidence, where 65% is obtained by taking the midpoint between $Q_{HII}$ = 80% and 50%. The weakness of the constraint from this galaxy is mainly due to the low SNR of the observation. While high-SNR observations would likely improve the result, the constraint can also be tightened by combining results from other galaxies at similar redshifts, treating the prediction probability as likelihood.

### 4.4. JADES-58975

JADES-58975 is an $M_{UV} = -20.3$ galaxy at $z=9$, with a publicly available spectrum showing an exceptionally high continuum SNR of ~40. This galaxy is considered to be an obscured AGN (J. Scholtz et al. 2023). The Ly$\alpha$ emission of this galaxy has not been measured. For the analysis, we use four spectral bins at $dv_\alpha$ = 6400, 10,700, 15,100, and 19,500 km s$^{-1}$. One additional bin, at $dv_\alpha$ = 2100 km s$^{-1}$, is excluded because its wavelength range, $dv_\alpha$ ≈ 0–4200 km s$^{-1}$, includes uncertain flux around the Lyman break. Figure 11 presents the raw spectrum, normalized mock spectra, decision boundary in the 2D data space, and the $Q_{HII}$ prediction probability.

The spectrum is well fitted by a power law with small deviations (panel (a)). In panel (d), the normalized spectrum falls below all of the mock DLA-free spectra, indicating that this galaxy likely has DLA absorption. The DLA presence is further confirmed by its location in the data space (panel (e)), which is within the region populated by DLA-added mocks. The DLA-RFC detects DLA absorption with 100% confidence. JADES-58975 is well positioned in the $Q_{HII}$ = 1.6% region, with a 90% prediction probability for $Q_{HII}$ = 1.6% from the Q-RFC (panel (f)).

In Section 3.3, we have shown that the Q-RFC accuracy for galaxies with a DLA is generally low (<40%) for spectra with
<figure><img src="image_15.png" /><figcaption>Figure 11. Same as panels (a), (d), (e), and (f) of Figure 8 for JADES-58975 at z = 9.4.</figcaption></figure>

SNR=20, but $Q_{HII}$ = 1.6% samples are classified with relatively higher accuracy (~55%) even when DLA-added mocks are used (Figure 5(b)). Additionally, the continuum SNR of JADES-58975 is higher than 20, suggesting that the Q-RFC results for this galaxy can be considered reliable.

It is worth noting that the presence of a DLA makes strong Lyα emission, which we could not subtract for this galaxy, less likely, as DLA opacity significantly suppresses the flux. However, there are cases where galaxies exhibit both strong Lyα emission and DLA absorption (e.g., K. E. Heintz et al. 2025; J. Witstok et al. 2025), possibly due to Lyα emission originating from structures outside the DLA obscuration in some galaxies. Thus, it is possible that the results here may have been influenced by unsubtracted Lyα emission leakage. If such leakage exists and is subtracted, the resulting spectrum would show even lower flux, changing the Q-RFC result.

## 5. Summary and Discussion

Using the state-of-the-art CoDa II simulation of the reionization era, we have developed a methodology to constrain the reionization history through the Lyα DW absorption feature in the spectra of LBGs. Leveraging the simulated galaxies and HII map from CoDa II, we created approximately $10^6$ mock galactic Lyα DW spectra across five values of $Q_{HII}$ and trained RFCs to estimate the most likely value of $Q_{HII}$ and detect the presence of DLAs in input spectra.

We found that individual spectra can produce tight constraints on $Q_{HII}$ when the IGM is highly neutral ($Q_{HII} < 50\%$), provided that the rest-frame flux data at wavelengths between 1220 and 1270 Å have a continuum SNR of $\gtrsim 20$, and no DLA is present in the sight line. We applied our model to four LBG spectra with high SNR ($\gtrsim 10$) to confirm the model's effectiveness. During this analysis, we also showed that accurately subtracting the leakage from Lyα emission is crucial for low-resolution spectra, such as those obtained using the NIRSpec PRISM.

### 5.1. Reionization History Constraint

Figure 12 summarizes the reionization history constraints from Section 4 along with reionization history in the CoDa II simulation and that from the LATE START/LATE END, EARLY START/LATE END, and EARLY START/EARLY END models of C. Cain et al. (2025). The upper bounds from JADES-12637, JADES-58975, and GN-z11 are calculated for 68% confidence limit (CL), assuming the five $Q_{HII}$ values (100%, 80%, 50%, 13.2%, and 1.6 %) represent intervals [100, 90]%, [90, 65]%, [65, 31.6]%, [31.6, 7.4]%, and [7.4, 0]%, respectively, where

<figure><img src="image_16.png" /><figcaption>Figure 12. Reionization history constraints derived from our RFC model applied to JWST NIRSpec/PRISM data on Lyman break galaxies (LBGs) are shown as red arrows, indicating the upper bounds with 68% CL from JADES-12637, JADES-58975, and GN-z11, respectively. Constraints from another study using LBGs by E. Curtis-Lake et al. (2023) and H. Umeda et al. (2024) are shown as magenta error bars. Constraints from the DW feature in other observables are shown as a green error bar for GRBs (T. Totani et al. 2006) and blue error bars for QSOs (J. Schroeder et al. 2013; F. B. Davies et al. 2018; B. Greig et al. 2019; F. Wang et al. 2020). These error bars are for 68% CL except for the leftmost blue bar, with a thicker line from J. Schroeder et al. (2013), calculated for 95% CL. The reionization history of CoDa II is shown as a black solid line, and that of the LATE START/LATE END, EARLY START/LATE END, and EARLY START/EARLY END models from C. Cain et al. (2025) are shown as blue, black, and red dotted lines, respectively.</figcaption></figure>

we took midpoints between snapshots to determine the intervals. We assume the likelihood within the interval, as predicted by the Q-RFC, is uniformly distributed. The resulting upper bounds are 59.6%, 5.6%, and 18.5% at $z = 7.7, 9.4$, and $10.6$, respectively. We do not include the constraint from CEERS-1029 in the figure due to its seemingly inaccurate continuum fitting. Note that 68% CL constraints from other studies using the DW feature of Gamma-ray bursts (GRBs), QSOs, and LBGs are also shown.

While the constraint from GN-z11 does not rule out any model show here, the upper bound from JADES-58975 disfavors the EARLY START/LATE END and EARLY START/EARLY END models, where $Q_{HII}$ rises above 20% before $z = 10$. The constraint from JADES-12637 is weak, but disfavors the EARLY START/EARLY END scenario. The CoDa II and LATE START/LATE END models reach ≈15% ionization by $z = 8$ and exhibit similar reionization histories down to $z = 7$, both of which are consistent with the constraints from
this work and with most of the other DW constraints shown in the figure. The tight constraint from JADES-58975 highlights the potential effectiveness of high-SNR observations of DWs, analyzed using reionization simulations, which provides a prior on the possible HII bubble size distribution for a given $Q_{HII}$. However, the presence of a DLA and the absence of an $Ly\alpha$ emission-line measurement for this galaxy suggest higher uncertainty in the result. Excluding the constraint from this galaxy would also allow for the EARLY START/EARLY END scenario, in which ~40% of the IGM is reionized by $z=8$.

We note that numerous constraints from other $Ly\alpha$ observables (e.g., C. A. Mason et al. 2018; A. Hoag et al. 2019; I. Jung et al. 2020; Planck Collaboration et al. 2020; L. R. Whitler et al. 2020; A. M. Morales et al. 2021; P. Bolan et al. 2022; S. Bruton et al. 2023; T. Morishita et al. 2023; M. Nakane et al. 2024; M. Tang et al. 2024) are not shown in the figure. While those constraints are often inconsistent with each other, most disfavor scenarios in which reionization ends early, such as the EARLY START/EARLY END model (see also discussions in C. Cain et al. 2025).

We also acknowledge that the reionization history constraints from this work still have remaining uncertainties, which could not be quantified in the scope of this work. Therefore, the constraints presented here are preliminary and subject to refinement in future studies. In the following sections, we discuss these uncertainties (Section 5.2) and outline future prospects (Section 5.3).

## 5.2. Challenges and Possible Solutions

DLAs: In theory, DLA opacity is distinguishable from IGM opacity due to the difference in wavelength dependence. Determining DLA presence in a target galaxy can be done even with low-SNR (~5) spectra from NIRSpec PRISM-like instruments. However, constraining $Q_{HII}$ from DLA-contaminated spectra requires a high SNR (~100), which is only attainable from rare, bright objects such as JADES-58975 (see also discussions in M. Huberty et al. 2025). A practical approach involves initially identifying galaxies without DLA signatures using low-SNR ($\lesssim 20$) observations. Once DLA-free candidates are identified, they can be followed up with higher-SNR ($\gtrsim 20$) observations for input into our model. Additionally, analyzing multiple spectra from galaxies at similar redshifts can help exclude DLA-contaminated results and improve overall reliability.

Reionization model: The mock $Ly\alpha$ spectra in this work are based on the reionization model of the CoDa II simulation, which may be inaccurate if the assumed ionizing emissivity of galaxies deviates from reality. The HII bubble size distribution, a key factor influencing DW strength, depends on whether ionizing photons predominantly originate from rare, bright galaxies or numerous, faint galaxies (see, e.g., T.-Y. Lu et al. 2024). Another highly uncertain aspect of reionization is the ionizing efficiency of minihalos below the atomic-cooling threshold ($\lesssim 10^8 M_\odot$), which are unresolved in most cosmological simulations. These sources can dominate reionization at $z \gtrsim 10$ before the buildup of the Lyman-Werner (LW) background and lead to a substantially different reionization geometry (K. Ahn et al. 2012; K. Ahn & P. R. Shapiro 2021). After star formation is suppressed in these minihalos by the LW background, they can act as sinks for ionizing photons, altering the HII bubble size distribution in the later stages of reionization (Y. Mao et al. 2020; M. Bianco et al. 2021).

Furthermore, the DW strength depends not only on $Q_{HII}$ but also on redshift, which, in principle, requires multiple simulations with different reionization histories to model the possible range of DW strength at a given redshift. This aspect was not accounted for in this work to perform the calculation with a single simulation. Additionally, the simulation volume may be insufficient to capture rare, bright galaxies such as GN-z11. While observational uncertainties currently dominate, these modeling uncertainties need to be addressed as data quality improves.

*Continuum shape.* This work fits the galactic continuum with a simple power law, not accounting for uncertainty in its shape. For instance, CEERS-1029 shows deviations from a power-law continuum, potentially biasing our $Q_{HII}$ constraints. Highly star-forming galaxies during reionization can exhibit significant contributions from nebular continuum (A. J. Cameron et al. 2023; G. Roberts-Borsani et al. 2024), which can masquerade a DW near $Ly\alpha$, adding degeneracy in constraining the reionization history (H. Katz et al. 2024). A potential solution is to identify galaxies with nebular continuum signatures through detailed SED modeling and focus deep observations on galaxies less likely to have this component. As demonstrated, tight constraints can be achieved from a single galaxy with high-SNR data if there is no uncertainty in the continuum shape. Prioritizing galaxies with minimal continuum shape uncertainty and obtaining high-SNR data for these samples is recommended.

*$Ly\alpha$ emission.* Spectral bins closer to the $Ly\alpha$ resonance can provide stronger constraints on $Q_{HII}$ due to increased DW absorption near the resonance. However, these bins are more susceptible to contamination from $Ly\alpha$ emission. As shown in the cases of GN-z11 and JADES-12637, strong emissions with EW greater than 10 Å can render the analysis ineffective (A. J. Bunker et al. 2023). Therefore, accurately measuring the $Ly\alpha$ EW or setting stringent upper limits using high-resolution spectroscopy is crucial for this method to be effective.

## 5.3. Future Prospects

JWST will continue to observe galactic $Ly\alpha$ damping wings (DW) throughout its operational lifetime. The CAPERS program (P.I.: Mark Dickinson) on JWST aims to observe approximately 2000 objects, dedicating 3 hr of exposure per target with NIRSpec/PRISM. This extensive data set will be well suited for the methodologies developed in this work. CAPERS also plans to conduct medium-resolution spectroscopy on a subset of these targets, which we demonstrated to be crucial for accurately subtracting $Ly\alpha$ emission leakage. Additionally, EUCLID will be able to observe hundreds of $Ly\alpha$ DWs (Euclid Collaboration et al. 2022). By analyzing high-SNR spectra from selected targets at $z \sim 7-9$, we will be able to distinguish between different reionization scenarios with high confidence, enhancing our understanding of the star formation rate during this epoch.

Future work will involve improving our model's performance by incorporating the complete data set from the new CoDa III simulation, which features higher spatial resolution and a reionization history that better aligns with recent observational constraints. This will help reduce modeling uncertainties. In the updated model, we plan to utilize a larger number of snapshots and adopt regression models to refine the prediction of $Q_{HII}$, minimizing discretization effects. This approach will also enable us to regress the DLA column
density directly from the spectra, allowing us to construct detailed DLA statistics and explore the conditions of the CGM and ISM in the high-redshift Universe.

Acknowledgments

H.P. thanks Junhwan Choi and Charlot Mason for insightful comments on this work and Christopher Cain for providing his reionization model data. H.P. thanks the organizers of the Cosmic Dawn at High Latitude conference held in Stockholm, Sweden between 2024 June 24–28, which allowed for helpful discussion with the participants. G.Y. acknowledges Ministerio de Ciencia e Innovación (Spain) for partial financial support under research grant PID2021-122603NB-C21. K.A. is supported by NRF-2021R1A2C1095136 and RS-2022-00197685. H.Y. is supported by MEXT/JSPS KAKENHI grant No. 21H04489 and JST FOREST Program, grant No. JP-MJFR202Z.

ORCID iDs

Hyunbae Park https://orcid.org/0000-0002-7464-7857

Intae Jung https://orcid.org/0000-0003-1187-4240

Hidenobu Yajima https://orcid.org/0000-0002-1319-3433

Paul R. Shapiro https://orcid.org/0000-0002-0410-3045

Kyungjin Ahn https://orcid.org/0000-0003-3974-1239

Romain Teyssier https://orcid.org/0000-0001-7689-0933

Gustavo Yepes https://orcid.org/0000-0001-5031-7936

Ilian T. Iliev https://orcid.org/0000-0002-5174-1365

References

Ahn, K., Iliev, I. T., Shapiro, P. R., et al. 2012, ApJL, 756, L16

Ahn, K., & Shapiro, P. R. 2021, ApJ, 914, 44

Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, ApJL, 951, L22

Atek, H., Richard, J., Kneib, J.-P., & Schäerer, D. 2018, MNRAS, 479, 5184

Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Natur, 553, 473

Bianco, M., Iliev, I. T., Ahn, K., et al. 2021, MNRAS, 504, 2443

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

Bolan, P., Lemaux, B. C., Mason, C., et al. 2022, MNRAS, 517, 3263

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 811, 140

Bouwens, R. J., Oesch, P. A., Illingworth, G. D., Ellis, R. S., & Stefanon, M. 2017, ApJ, 843, 129

Boyett, K., Trenti, M., Leethochawalit, N., et al. 2024, NatAs, 8, 657

Breiman, L. 2001, Mach. Learn., 45, 5

Bruton, S., Lin, Y.-H., Scarlata, C., & Hayes, M. J. 2023, ApJL, 949, L40

Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2024, A&A, 690, A288

Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023, A&A, 677, A88

Cain, C., Lopez, G., D'Aloisio, A., et al. 2025, ApJ, 980, 83

Cameron, A. J., Katz, H., Rey, M. P., & Saxena, A. 2023, MNRAS, 523, 3516

Chen, H. 2024, MNRAS, 528, L33

Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, NatAs, 7, 622

Davies, F. B., Hennawi, J. F., Bañados, E., et al. 2018, ApJ, 864, 142

D'Eugenio, F., Cameron, A. J., Scholtz, J., et al. 2025, ApJS, 277, 4

Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv:2306.02465

Endsley, R., & Stark, D. P. 2022, MNRAS, 511, 6042

Euclid Collaboration, Scaramella, R., Amiaux, J., et al. 2022, A&A, 662, A112

Finkelstein, S. L., Papovich, C., Dickinson, M., et al. 2013, Natur, 502, 524

Finkelstein, S. L., Leung, G. C. K., Bagley, M. B., et al. 2024, ApJL, 969, L2

Finkelstein, S. L., Ryan, R. E., Jr., Papovich, C., et al. 2015, ApJ, 810, 71

Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, SSRv, 123, 485

Greig, B., Mesinger, A., & Bañados, E. 2019, MNRAS, 484, 5094

Heintz, K. E., Brammer, G. B., Watson, D., et al. 2025, A&A, 693, A60

Hennawi, J. F., Kist, T., Davies, F. B., & Tamanas, J. 2024, arXiv:2406.12070

Hoag, A., Bradač, M., Huang, K., et al. 2019, ApJ, 878, 12

Hsiao, T. Y.-Y., Coe, D., Abdurro'uf, et al. 2023, ApJL, 949, L34

Huberty, M., Scarlata, C., Hayes, M. J., & Gazagnes, S. 2025, arXiv:2501.13899

Hutter, A., Dayal, P., Partl, A. M., & Müller, V. 2014, MNRAS, 441, 2861

Jung, I., Finkelstein, S. L., Dickinson, M., et al. 2020, ApJ, 904, 144

Jung, I., Finkelstein, S. L., Livermore, R. C., et al. 2018, ApJ, 864, 103

Jung, I., Papovich, C., Finkelstein, S. L., et al. 2022, ApJ, 933, 87

Katz, H., Cameron, A. J., Saxena, A., et al. 2024, arXiv:2408.03189

Keating, L. C., Puchwein, E., Bolton, J. S., Haehnelt, M. G., & Kulkarni, G. 2024, MNRAS, 531, L34

Kist, T., Hennawi, J. F., & Davies, F. B. 2024, arXiv:2406.12071

Larson, R. L., Finkelstein, S. L., Hutchison, T. A., et al. 2022, ApJ, 930, 104

Lee, H.-W. 2013, ApJ, 772, 123

Lewis, J. S. W., Ocvirk, P., Sorce, J. G., et al. 2022, MNRAS, 516, 3389

Lu, T.-Y., Mason, C. A., Hutter, A., et al. 2024, MNRAS, 528, 4872

Lukić, Z., Stark, C. W., Nugent, P., et al. 2015, MNRAS, 446, 3697

Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2024a, A&A, 691, A145

Maiolino, R., Scholtz, J., Witstok, J., et al. 2024b, Natur, 627, 59

Mao, Y., Koda, J., Shapiro, P. R., et al. 2020, MNRAS, 491, 1600

Mason, C. A., Chen, Z., Stark, D. P., et al. 2025, arXiv:2501.11702

Mason, C. A., Treu, T., Dijkstra, M., et al. 2018, ApJ, 856, 2

Miralda-Escudé, J. 1998, ApJ, 501, 15

Morales, A. M., Mason, C. A., Bruton, S., et al. 2021, ApJ, 919, 120

Morishita, T., Roberts-Borsani, G., Treu, T., et al. 2023, ApJL, 947, L24

Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Natur, 474, 616

Nakane, M., Ouchi, M., Nakajima, K., et al. 2024, ApJ, 967, 28

Napolitano, L., Pentericci, L., Santini, P., et al. 2024, A&A, 688, A106

Ocvirk, P., Aubert, D., Sorce, J. G., et al. 2020, MNRAS, 496, 4087

Ocvirk, P., Lewis, J. S. W., Conaboy, L., et al. 2024, arXiv:2409.05946

Ono, Y., Ouchi, M., Mobasher, B., et al. 2012, ApJ, 744, 83

Ouchi, M., Ono, Y., & Shibuya, T. 2020, ARA&A, 58, 617

Park, H., Jung, I., Song, H., et al. 2021, ApJ, 922, 263

Park, H., Lukić, Z., Sexton, J., Alvarez, M. A., & Shapiro, P. R. 2024, ApJ, 969, 46

Pentericci, L., Vanzella, E., Castellano, M., et al. 2018, A&A, 619, A147

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A16

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6

Roberts-Borsani, G., Treu, T., Shapley, A., et al. 2024, ApJ, 976, 193

Schenker, M. A., Ellis, R. S., Konidaris, N. P., & Stark, D. P. 2014, ApJ, 795, 20

Scholtz, J., Maiolino, R., D'Eugenio, F., et al. 2023, arXiv:2311.18731

Schroeder, J., Mesinger, A., & Haiman, Z. 2013, MNRAS, 428, 3058

Shibuya, T., Kashikawa, N., Ota, K., et al. 2012, ApJ, 752, 114

Smith, A., Kannan, R., Garaldi, E., et al. 2022a, MNRAS, 512, 3243

Smith, A., Kannan, R., Tacchella, S., et al. 2022b, MNRAS, 517, 1

Song, M., Finkelstein, S. L., Livermore, R. C., et al. 2016, ApJ, 826, 113

Tang, M., Stark, D. P., Topping, M. W., Mason, C., & Ellis, R. S. 2024, ApJ, 975, 208

Terp, C., Heintz, K. E., Watson, D., et al. 2024, A&A, 690, A70

Tilvi, V., Malhotra, S., Rhoads, J. E., et al. 2020, ApJL, 891, L10

Totani, T., Kawai, N., Kosugi, G., et al. 2006, PASJ, 58, 485

Treu, T., Schmidt, K. B., Trenti, M., Bradley, L. D., & Stiavelli, M. 2013, ApJL, 775, L29

Treu, T., Trenti, M., Stiavelli, M., Auger, M. W., & Bradley, L. D. 2012, ApJ, 747, 27

Umeda, H., Ouchi, M., Nakajima, K., et al. 2024, ApJ, 971, 124

Vanzella, E., Pentericci, L., Fontana, A., et al. 2011, ApJL, 730, L35

Wang, F., Davies, F. B., Yang, J., et al. 2020, ApJ, 896, 23

Whitler, L. R., Mason, C. A., Ren, K., et al. 2020, MNRAS, 495, 3602

Witstok, J., Jakobsen, P., Maiolino, R., et al. 2025, Natur, 639, 897

Wold, I. G. B., Malhotra, S., Rhoads, J., et al. 2022, ApJ, 927, 36

Yang, J., Wang, F., Fan, X., et al. 2020, ApJ, 904, 26

Zitrin, A., Labbé, I., Belli, S., et al. 2015, ApJL, 810, L12