# The JWST/PASSAGE Survey: Testing Reionization Histories with JWST's First Unbiased Survey for Ly$\alpha$ Emitters at Redshifts 7.5–9.5

Axel Runnholm<fnref n="1" />, Matthew J. Hayes<fnref n="1" />, Vihang Mehta<fnref n="2" />, Matthew A. Malkan<fnref n="3" />, Claudia Scarlata<fnref n="4" />, Kalina V. Nedkova<fnref n="5" />,<fnref n="6" />, Marc Rafelski<fnref n="5" />,<fnref n="6" />, Benedetta Vulcani<fnref n="7" />, Mason Huberty<fnref n="4" />, E. Christian Herenz<fnref n="8" />, Anne Hutter<fnref n="9" />,<fnref n="10" />, Sean Bruton<fnref n="11" />, Ayan Acharyya<fnref n="7" />, Hakim Atek<fnref n="12" />, Ivano Baronchelli<fnref n="13" />, Andrew J. Battisti<fnref n="14" />,<fnref n="15" />,<fnref n="16" />, Maruša Bradač<fnref n="17" />,<fnref n="18" />, Andrew J. Bunker<fnref n="19" />, Y. Sophia Dai<fnref n="20" />, Clea Hannahs<fnref n="3" />, Farhanul Hasan<fnref n="6" />, Keunho J. Kim<fnref n="2" />, Nicha Leethochawalit²¹, Yu-Heng Lin<fnref n="2" />, Michael J. Rutkowski<fnref n="22" />, Alberto Saldana-Lopez<fnref n="1" />, Zahra Sattari<fnref n="2" />, and Xin Wang<fnref n="23" />,²⁴,<fnref n="25" />

## Abstract

Ly$\alpha$ emission is one of a few observable features of galaxies that can trace the neutral hydrogen content in the Universe during the Epoch of Reionization (EoR). To accomplish this, we need an efficient way to survey for Ly$\alpha$ emitters (LAEs) at redshifts beyond 7, requiring unbiased emission-line observations that are both sufficiently deep and wide to cover enough volume to detect them. Here we present results from PASSAGE—a pure-parallel JWST/NIRISS slitless spectroscopic survey to detect LAEs deep into the EoR, without the bias of photometric preselection. We identify four LAEs at 7.5 ≤ z ≤ 9.5 in four surveyed pointings and estimate the luminosity function (LF). We find that the LF does show a marked decrease compared to post-reionization measurements, but the change is a factor of ≲10, which is less than expected from theoretical calculations and simulations, as well as observational expectations from the pre-JWST literature. Modeling of the intergalactic medium and expected Ly$\alpha$ profiles implies that these galaxies reside in ionized bubbles of ≳2 physical Mpc. We also report that in the four fields we detect {3, 1, 0, 0} LAEs, which could indicate strong field-to-field variation in the LAE distribution, consistent with a patchy H I distribution at z ~ 8. We compare the recovered LAE number counts with expectations from simulations and discuss the potential implications for reionization and its morphology.

*Unified Astronomy Thesaurus concepts: Reionization (1383); Lyman-alpha galaxies (978); High-redshift galaxies (734)*

## 1. Introduction

The Epoch of Reionization (EoR) is a time of great astrophysical interest since it represents both a major phase transition of the Universe and the earliest stages of galaxy

evolution (see, e.g., B. E. Robertson 2022, for a review). Despite this, relatively little is certain about how the transition proceeded, especially at redshifts (z) beyond 7. The Planck experiment (Planck Collaboration et al. 2020) used electron scattering measurements to constrain the midpoint of the EoR to be z ~ 7.5. Observations of the Ly$\alpha$ forest and damping wing observed in quasar spectra have placed the end of reionization between redshift 5 and 6 (G. Kulkarni et al. 2019; L. C. Keating et al. 2020; S. E. I. Bosman et al. 2022; B. Greig et al. 2024; G. D. Becker et al. 2024; Y. Qin et al. 2024; Y. Zhu

<figure><img src="image_1.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.

<footnote n="1">¹Stockholm University, Department of Astronomy and Oskar Klein Centre for Cosmoparticle Physics, AlbaNova University Centre, SE-10691, Stockholm, Sweden; matthew.hayes@astro.su.se</footnote>

<footnote n="2">²IPAC, Mail Code 314-6, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA</footnote>

<footnote n="3">³University of California, Los Angeles, Department of Physics and Astronomy, 430 Portola Plaza, Los Angeles, CA 90095, USA</footnote>

<footnote n="4">⁴Minnesota Institute for Astrophysics, University of Minnesota, 116 Church Street SE, Minneapolis, MN 55455, USA</footnote>

<footnote n="5">⁵Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA</footnote>

<footnote n="6">⁶Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="7">⁷INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy</footnote>

<footnote n="8">⁸Inter-University Centre for Astronomy and Astrophysics (IUCAA), Pune University Campus, Pune 411 007, India</footnote>

<footnote n="9">⁹Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen N, Denmark</footnote>

<footnote n="10">¹⁰Cosmic Dawn Center (DAWN), Copenhagen, Denmark</footnote>

<footnote n="11">¹¹California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA</footnote>

<footnote n="12">¹²Institut d'Astrophysique de Paris, CNRS, Sorbonne Université, 98bis Boulevard Arago, 75014, Paris, France</footnote>

<footnote n="13">¹³INAF—Istituto di Radioastronomia, Via Gobetti 101 40129 Bologna, Italy</footnote>

<footnote n="14">¹⁴International Centre for Radio Astronomy Research (ICRAR), University of Western Australia, M468, 35 Stirling Highway, Crawley, WA 6009, Australia</footnote>

<footnote n="15">¹⁵Australian National University, Research School of Astronomy and Astrophysics, Canberra, ACT 2611, Australia</footnote>

<footnote n="16">¹⁶ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia</footnote>

<footnote n="17">¹⁷University of Ljubljana, Department of Mathematics and Physics, Jadranska ulica 19, SI-1000 Ljubljana, Slovenia</footnote>

<footnote n="18">¹⁸Department of Physics and Astronomy, University of California, Davis, 1 Shields Avenue, Davis, CA 95616, USA</footnote>

<footnote n="19">¹⁹Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="20">²⁰The National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100101, People's Republic of China</footnote>

<footnote n="22">²²Department of Physics and Astronomy, Minnesota State University, Mankato, Trafton North 158, Mankato, MN 56001, USA</footnote>

<footnote n="23">²³School of Astronomy and Space Science, University of Chinese Academy of Sciences (UCAS), Beijing 100049, People's Republic of China
²⁴National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, People's Republic of China</footnote>

<footnote n="25">²⁵Institute for Frontiers in Astronomy and Astrophysics, Beijing Normal University, Beijing 102206, People's Republic of China</footnote>

<footnote>Received 2024 December 19; revised 2025 February 26; accepted 2025 March 11; published 2025 April 30</footnote>
et al. 2024). However, the start of the EoR, as well as whether it proceeded rapidly or more gradually, is still debated (S. L. Finkelstein et al. 2019; R. P. Naidu et al. 2020; P. Bolan et al. 2022; G. D. Becker et al. 2024; L. Napolitano et al. 2024). In fact, the shape of the transition itself may have strong implications for the sources that powered it, in terms of power, number density, and spatial distribution, and even for their nature (Y. Qin & J. S. B. Wyithe 2025).

The neutral intergalactic medium (IGM) consists mostly of atomic hydrogen, which makes studying it directly extremely challenging. The Square Kilometer Array (SKA; P. E. Dewdney et al. 2009; L. Koopmans et al. 2015) will attempt to target 21 cm emission from this epoch; however, even with this new facility, the problem will remain a difficult one. Therefore, we must turn to additional, more indirect tracers to better constrain the reionization history. One such tracer is Lyα emission, which is an intrinsically strong emission line (~68% of ionizing photons result in Lyα emission under normal gas conditions) from recombining hydrogen gas (see, e.g., M. Dijkstra 2014; M. Ouchi et al. 2020). As a result of its resonant nature, Lyα radiation interacts strongly with neutral atomic hydrogen gas. Therefore, in principle, it can be used to trace the column of atomic gas in the IGM and hence the neutral fraction of the Universe.

This has been done using various metrics, including the evolution of the Lyα luminosity function (LF; N. Kashikawa et al. 2011; I. G. B. Wold et al. 2022), the volume-averaged escape fraction (M. Hayes et al. 2011), the Lyα emitter (LAE) fraction (D. P. Stark et al. 2010, 2011; L. Pentericci et al. 2011; Y. Ono et al. 2012; H. Kusakabe et al. 2022), and the equivalent width (EW) distribution (C. A. Mason et al. 2018, 2019; P. Bolan et al. 2022). Previous works, such as M. Ouchi et al. (2008) and E. C. Herenz et al. (2019), have shown that the Lyα LF remains relatively constant between redshifts 3 and 6, although some changes at the extreme ends of the LF have been noted (e.g., T. T. Thai et al. 2023). However, observing LAEs and, consequently, LFs at higher redshifts has proven challenging, but not impossible, using ground-based narrowband searches. Early attempts indicated a dearth of Lyα detections (J. G. Cuby et al. 2007; J. P. Willis et al. 2008), which on its own suggested a significant drop in the LAE density at redshifts above 6. More recently, R. Itoh et al. (2018) and I. G. B. Wold et al. (2022) detected small numbers of LAEs and confirmed the reduction of the Lyα LF at $z \sim 7$. A significant reduction in the Lyα LF at the same $z$ was also found by the SILVERRUSH collaboration (M. Ouchi et al. 2018; S. Kikuta et al. 2023; H. Umeda et al. 2025).

Nevertheless, a reduction of the Lyα LF is insufficient to conclude that the Lyα optical depth is increasing purely because of a neutral IGM. This could also be due to the galaxy population being different at higher redshifts. Therefore, IGM inferences have relied on various comparisons to the UV properties of the galaxy population at the same redshift. For example, M. Hayes et al. (2011), R. Itoh et al. (2018), and I. G. B. Wold et al. (2022) compare the implied Lyα luminosity density (LD) evolution with redshift to the UV LD evolution. R. Itoh et al. (2018) found that Lyα evolves faster—an indication of neutral IGM damping wing absorption—whereas I. G. B. Wold et al. (2022) conclude that there is no indication of this at $z \sim 7$, implying a fully ionized IGM. C. A. Mason et al. (2018, 2019) and P. Bolan et al. (2022), on the other hand, rely on modeling the Lyα EW distribution relation to the UV

LF and conclude that the neutral fraction, $x_{HI}$, is about 0.6 at redshift 7, whereas the results of C. A. Mason et al. (2025) indicate that the Universe reaches those $x_{HI}$ closer to $z \sim 9$. Clearly, there are large uncertainties and discrepancies in EoR progression inferences, and more observations are needed, especially at higher redshift, to put the conclusions on solid footing.

The high-sensitivity infrared capabilities of JWST have opened a new window into this and are enabling the detection of LAEs at unprecedented redshifts. For instance, M. Nakane et al. (2024) collected 53 galaxies with confirmed Lyα emission from GLASS (T. Treu et al. 2022), CEERS (S. L. Finkelstein et al. 2023), and JADES (D. J. Eisenstein et al. 2023; A. J. Bunker et al. 2024) above $z \ge 7$, and M. Tang et al. (2024) present a compilation of more than 210 LAEs at $z \ge 6.5$. Such data sets have already been used to estimate the sizes of ionized regions around the observed galaxies, finding necessary radii of 0.1–1 pMpc (J. Witstok et al. 2024, 2025b), as well as tentative growth of ionized regions with decreasing redshift (M. J. Hayes & C. Scarlata 2023). J. Witstok et al. (2025a) recently reported the detection of a strong Lyα line at an impressively high redshift of 13, which was wholly unexpected owing to the high IGM neutral fractions expected to efficiently scatter Lyα at these redshifts.

Using such observations, we are now able to start pushing EoR constraints from the end stages all the way to when the process began and determine the reionization history. Different sources of reionization have different implications for how the process proceeds; whether it is powered by a few bright galaxies (oligarchic reionization; M. Sharma et al. 2016; R. P. Naidu et al. 2020) or by emission from many faint galaxies (also referred to as democratic reionization; S. L. Finkelstein et al. 2019; H. Atek et al. 2024) will affect how gradual it is (Y. Qin & J. S. B. Wyithe 2025), as well as its topology. Measuring this well will require large unbiased samples of LAEs in the EoR—in this work we demonstrate the feasibility of obtaining such samples with JWST/NIRISS.

One limitation of current high-redshift Lyα samples is that they are derived from a small number of relatively small area surveys. This means that any inferences from them are potentially susceptible to cosmic variance (CV) and may be part of the explanation of why Lyα-based $x_{HI}$ determinations at the same redshifts may disagree with each other (see, e.g., L. Napolitano et al. 2024). The most effective way to overcome CV is to observe many uncorrelated fields rather than expanding a single contiguous field (S. Bruton et al. 2023b). This is the primary driver of the Cycle 1 JWST program Parallel Application of Slitless Spectroscopy for the Analysis of Galaxy Evolution (PASSAGE; Malkan et al., submitted), which uses pure parallel observations with NIRISS wide-field slitless spectroscopy (WFSS) to conduct an unbiased search for emission-line galaxies. In this work we develop a technique for pure emission-line selection in grism spectroscopy and use the deepest PASSAGE data to search for faint Lyα emission to build up a Lyα LF at $7.5 \le z \le 9.5$.

In Section 2 we motivate why WFSS is the most efficient way to perform an unbiased Lyα survey at $z > 7$. In Section 3 we describe the observations and data reduction, as well as the bespoke process of source detection, calibration, and validation. We describe the selected LAE candidates in Section 4 and discuss the implications of our observations for reionization in Section 5. Sections 6 and 7 give a broader discussion and
conclusions, respectively. We use a concordance cosmology ($\Omega_M = 0.3, \Omega_\Lambda = 0.7, H_0 = 70$) throughout.

## 2. Grism-based Emission-line Selection

Although LAEs are now being identified at $z > 7$, characterizing the population as a whole is not straightforward. The reason for this is that all nebular $Ly\alpha$ from these redshifts is identified by follow-up spectroscopy of continuum-selected galaxies. This introduces a bias, as the resulting LF is sampled as the UV LF convolved with the EW distribution (see, e.g., M. Dijkstra & J. S. B. Wyithe 2012; M. Gronke et al. 2015, for discussion) rather than being derived from $Ly\alpha$ emission directly. This makes it difficult, if not impossible, to quantify the sample completeness in $Ly\alpha$ luminosities.

Ideally, we should select LAEs directly by their $Ly\alpha$ emission. Many $Ly\alpha$ LF studies have done this either using IFU spectroscopy below redshift ~6.6 (e.g., E. C. Herenz et al. 2019) or at higher redshifts, using narrowband selection from the ground (see, e.g., J. G. Cuby et al. 2007; J. P. Willis et al. 2008). At $z \gtrsim 6.5$ these studies become very challenging and the reachable depths are heavily limited by the sky background —skylines, absorption, and IR background impact sensitivity and available $z$-range, which limits the observable volume. Though some campaigns, for instance, UltraVISTA (H. J. McCracken et al. 2012; P. Laursen et al. 2019), have targeted large sky areas to mitigate this, line sensitivity remains low. Furthermore, low-$z$ interlopers, such as [O II], [O III], and $H\alpha$ emitters, limit the purity of the samples. At Earth–Sun L2, JWST near-infrared (NIR) observations do not suffer from airglow background, which allows much greater depths to be reached in comparable short integration times. Additionally, the use of NIRISS’s WFSS capabilities lets us cover a much larger range in redshift in a single observation than a narrowband filter and provides better interloper rejection, due to the longer-wavelength coverage in which other emission lines can be identified. Together they enable us to conduct efficient, unbiased NIR emission-line surveys for the first time.

In this work we set out to derive the $7.5 \le z \le 9.5$ $Ly\alpha$ LF from NIRISS WFSS data. However, standard WFSS processing uses a direct detection image to do preselection before spectral extraction. This means that in many cases “emission-line galaxies” are actually continuum selected rather than line selected, and the samples suffer the same bias as discussed above. This becomes an even more pronounced issue when using parallel observations, where the direct image covering the emission line may be shallow or even missing. Therefore, we develop methods in this paper for selecting emission-line galaxies directly in the dispersed images, adopting the philosophy of H. Atek et al. (2010; see also M. B. Bagley et al. 2017; A. J. Battisti et al. 2024). We describe the procedure in detail in Section 3.4. This allows us to do a pure emission-line-selected $Ly\alpha$ survey at high $z$.

## 3. Observations and Data Processing

Here we describe the technical details of the observations, discuss how we extract spectra, and verify our LAE candidates. Very briefly, we do this by searching for emission lines in the GR150C grism data using the F115W filter and confirming that the same line exists in the GR150R grism direction. We then verify that there is no detection in deep ancillary imaging data blueward of the emission line and that no other emission lines

are present in the rest of the spectral region covered by the grism data, including observations in the F150W and F200W filters. We specifically search for single, symmetric emission lines since the characteristic $Ly\alpha$ spectral profile is unresolved by NIRISS and the [O III]$_{4959,5007}$ line commonly shows an asymmetric profile. Readers interested primarily in the results may skip to Section 4.

### 3.1. Observations

The main data set for this study comes from the PASSAGE survey (Malkan et al. 2025, submitted), which consists in total of 63 parallel pointings observed with NIRISS with various combinations of the F115W, F150W, and F200W filters and a wide range of depths. Here we search for high-redshift $Ly\alpha$-emitting sources, which are naturally faint in both continuum and the $Ly\alpha$ line. This places a limit on which fields can be used, since many observations are quite shallow, with less than 1.5 hr of exposure time per grism and filter combination. For this survey we require a minimum of 2 hr of exposure time in each dispersion direction in the F115W filter. Since we are working at close to the detection limit of the data, we require coverage in both NIRISS grism orients—GR150C and GR150R—to confirm the line detections and establish a wavelength calibration. Together, these requirements leave four initial parallel fields that are of interest for this work. In addition to the PASSAGE data, we also use archival NIRISS data in the same WFSS configurations from the NGDEEP program (M. B. Bagley et al. 2024), which covers the Hubble Ultra Deep Field.

In order to further validate single-line emitters as LAEs, we also need deep ancillary data that lie blueward of our F115W data, in order to verify that there is no emission blueward of the emission line—as expected for LAEs (see Section 3.5 for details). We search for such imaging in both Hubble Space Telescope (HST) archives and the Suprime-Cam Legacy Archive. For one field only very shallow Sloan Digital Sky Survey imaging exists, which means that we cannot exclude low-$z$ interlopers, and therefore we do not consider this field further. This leaves us with a total of three PASSAGE fields—Par28, Par50, and Par52 in the nomenclature of the APT file—and the NGDEEP field covering a total area of $\approx 19$ arcmin$^2$. A summary of the exposure times of the four fields is given in Table 1. The last two columns of Table 1 detail the ancillary data that we have for each of our four fields together with the observational depth in terms of both the $1\sigma$ sky noise and AB magnitude. We characterize the sky noise by placing $\sim 500$ apertures, each with a radius of 5 pixels, in sky regions of the image and calculating the rms of the measured fluxes.

### 3.2. Data Reduction

The data were obtained from the MAST and reduced by the standard pipeline to obtain rate files. The rate files were further processed using Grizli (G. Brammer 2019, 2023), to produce combined frames for direct imaging, and GR150R + GR150C dispersed images for each visit, using the Astro-Drizzle routine of DrizzlePac (S. Gonzaga et al. 2012; S. L. Hoffmann et al. 2021).

### 3.3. Emission-line Detection

We develop a source detection algorithm to identify galaxies by their line emission in the dispersed spectral images using the following procedure:
<table><caption>**Table 1**
Observational Log</caption><thead><tr><th rowspan="2">Parallel</th><th colspan="3">F115W</th><th colspan="3">F150W</th><th colspan="3">F200W</th><th colspan="3">Ancillary Data</th></tr><tr><th>Direct</th><th>GR150C</th><th>GR150R</th><th>Direct</th><th>GR150C</th><th>GR150R</th><th>Direct</th><th>GR150C</th><th>GR150R</th><th>Instrument</th><th>Filter</th><th>1σ Depth</th></tr></thead><tbody><tr><td>Par28</td><td>6442</td><td>34014</td><td>42517</td><td>5153</td><td>8503</td><td>17007</td><td>5153</td><td>6442</td><td>12884</td><td>HST/ACS</td><td>F814W</td><td>29.5</td></tr><tr><td>Par50</td><td>...</td><td>21645</td><td>21645</td><td>7344</td><td>...</td><td>14173</td><td>7344</td><td>...</td><td>14172</td><td>Suprime-Cam<sup>a</sup></td><td>R-band</td><td>28.5</td></tr><tr><td>Par52</td><td>...</td><td>21645</td><td>21645</td><td>7344</td><td>...</td><td>14173</td><td>7344</td><td>...</td><td>14172</td><td>Suprime-Cam<sup>a</sup></td><td>R-band</td><td>30.3</td></tr><tr><td>NGDEEP</td><td>5411</td><td>46962</td><td>46962</td><td>1803</td><td>20872</td><td>20872</td><td>1803</td><td>15654</td><td>15654</td><td>HST/ACS</td><td>F814W</td><td>32.1</td></tr></tbody></table>

Notes. Integration times are given in seconds, and limiting magnitudes in the AB system (J. B. Oke & J. E. Gunn 1983).

<sup>a</sup> Suprime-Cam data were retrieved from the Suprime-Cam Legacy Archive (S. Gwyn 2020) (https://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/scla/).

1. We create a source detection image by (1D) unsharp masking. We create a smoothed image by running a rolling median filter (width of 30 pixels) along each pixel row in the dispersion direction and subtracting it from the spectral image. The spectral traces are slightly curved in practice, but this curvature is significantly less than 1 pixel along the first order, so we do not account for this curvature. The median filtering efficiently removes continuum light but leaves emission lines as point sources in the subtracted image.

2. We run *Source Extractor* (E. Bertin & S. Arnouts 1996) as implemented in the Python library *sep* (K. Barbary 2016) on the median-filtered image and create a catalog of potential emission lines. We use a *thresh* of 3 and *minarea* of 5, which yields a very permissive candidate list.

3. Next, we create a catalog of bright sources in the direct image, again with *Source Extractor*. From this catalog we calculate the positions of expected zero orders, which we then cross-correlate with our emission-line catalog (from step 2). Sources that fall within a 5-pixel radius of an expected zero order are excluded from the emission-line catalog.

4. We then visually inspect the direct images, 2D spectra, and extracted 1D spectra of the emission-line catalog sources. We classify them as artifacts, low-z emission-line sources (usually with multiple emission lines such as [O II], [O III], or Hα), or single-line emitter candidates. With the very relaxed source extractor settings, we have approximately 800–1000 candidate sources per NIRISS field. Approximately 80%–90% of these are artifacts, commonly from poorly masked hot pixels, or artifacts produced by the median filtering close to the edges of bright spectra.

This procedure can in principle use any of the three NIRISS bands as the detection band, but we only use F115W since it covers Lyα emission between z ~ 7.5 and z ~ 9.5. Future work will extend this search to the longer wavelength bands.

## 3.4. Calibration and Extraction

### 3.4.1. Wavelength Calibration

Traditional grism extraction starts from detections in direct imaging and then uses this source position to model the source spectrum trace and the positions of surrounding objects to model the potential source contamination. However, we need to search for very faint line-only sources, with varying relative depths in the direct and dispersed images, and in some cases even with missing F115W direct images. Thus, not all sources detected in emission lines will be detected in the direct image.

We therefore use spectral images dispersed in both the R and C directions to validate sources and determine spatial location of sources for wavelength calibration. Our procedure is as follows:

1. From the emission-line detection location in the initial dispersed image we calculate all possible origin positions in the direct image. This equates to one position per assumed observed wavelength of the line. We refer to this region as the possible origin line (POL).

2. Given the POL in the direct image, we can calculate the possible locations at which the emission line must fall in the second grism orient.

3. We then perform source detection along this line in the second grism, allowing offsets of a maximum of 2.5 pixels on either side of the line to account for possible calibration offsets and shifts because of irregular morphologies.

4. Finally, we trace this second line detection back into the direct image and find the point at which the implied POLs intersect. This gives us the position that matches the line wavelengths and also the coordinate we use for the final trace calculations and extractions.

In order to model the trace, we use the *grismconf* formalism as implemented in the *GRIZLI* code (G. Brammer 2019) together with the wavelength calibration files derived by N. Pirzkal et al. (2024), who report a wavelength calibration uncertainty of less than half a pixel.

### 3.4.2. Spectral Extraction and Background Subtraction

We perform our spectral extractions in combined drizzled images produced using the *GRIZLI* preprocessing steps. The spectral extraction is done using an optimal extraction formalism, meaning that the spectrum is weighted in the cross-dispersion direction by a combination of the point-spread function (PSF) size and the size of the source. The source size is measured using a Gaussian fit in the cross-dispersion direction at the position of the detected emission line. However, we find that all LAE candidates are consistent with being unresolved or marginally resolved, so the extraction weighting is simply the PSF profile, which we approximate as a Gaussian with standard deviation 0.6 pixels.

Due to both contamination by nearby sources and an uneven background in the NIRISS data, local backgrounds can be significant. We therefore do a relatively simple local background subtraction. For each pixel along the predicted trace position we select two 3-pixel regions outside $2\sigma$ of the source width in the cross-dispersion direction. We then fit a straight line to these points to model the background at the trace pixel. This background is then subtracted from the data before the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Data for P50CLAE1. The top left panels show the G150C 2D and 1D spectra. The blue line indicates the central position of the trace, and the red ellipse highlights the approximate position of the emission line. The bottom row shows the same for the G150R grism. For this grism all three filters are shown in the columns. Extracted 1D spectra of the source are shown in $F_\lambda \times 10^{20}$. The top right panels show the direct image cutouts of the source from the F150W and F200W filters together with a red circle indicating the source coordinate used for extraction.</figcaption></figure>

weighted extraction. We find that this simplistic approach in general performs well, but if the contaminating source spectrum contains strong emission lines, it can lead to oversubtraction of the spectrum. We flag/remove such cases during visual inspection.

### 3.5. LAE Candidate Validation

In order to securely identify single-line emitters as real sources, we require that the emission line is detected in both grisms at $>2\sigma$ confidence, with line strengths and widths in both grisms consistent within approximately 50%. For parallels with deep F115W direct images (Par28 and NGDEEP), we also consider a direct source detection to lend further credence to the source, but for Par50 and Par52 we do not have direct F115W imaging.

Single-line emitters in our data most likely have two primary source populations: Hα emission from $z \sim 0.7$ and Lyα from $z \ge 7.5$. Other strong lines like [O II]$_{3727,29}$ and [O III]$_{5007}$ would also show other strong emission lines in the redder bands. Bright Hα emitters may also be distinguished by other lines, such as [S II]$_{6717,31}$, [S III]$_{9069,9531}$, and He I$_{1.087}$, but detection of these much fainter lines cannot be guaranteed. Differentiating these populations is crucial. The primary differentiation between the populations is the expected lack of any continuum emission blueward of a Lyα line. The faintness of the sources means that continuum is often undetected in the spectra. In order to attempt to lend confidence

to our candidates, we search for archival data from both HST and ground-based telescopes, primarily Subaru, in filters blueward of F115W. For all the parallels considered here we have collected deep Suprime-Cam imaging from Subaru from the Suprime-Cam Legacy Archive (S. Gwyn 2020), and for Par28 we also have deep HST imaging from the Cosmic Evolution Survey (COSMOS; A. M. Koekemoer et al. 2007; N. Scoville et al. 2007).

To quantify the depth of the ancillary imaging and confirm detections and nondetections of our LAE candidates, we need accurate noise level estimates. We obtain this for each ancillary image by measuring the fluxes in 3-pixel apertures placed in empty regions of the image and then calculating the rms of these measurements (right columns of Table 1). We then measure the flux at the expected source location in a matched aperture and require that flux to be below the $1\sigma$ sky level.

## 4. LAE Candidate Sample

After candidate selection and validation, we have reduced a sample of ~4000 initial selections to a final sample of four candidate line emitters with redshifts between 7.6 and 9.5. In this section we present measurements of these emitters, as well as the implied populations statistics.

The direct images, 2D spectra, and extracted 1D spectra are shown in Figures 1–4. First, we present how we characterize the sources together with some discussion of each individual
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Same as Figure 1, but for P50CLAE2.</figcaption></figure>

candidate, and then we discuss sample properties, such as completeness and implied LFs.

## 4.1. Source Characterization

To measure the line flux, we fit the extracted spectra with single Gaussian line profiles for each dispersion direction separately. While a simple background subtraction is done during the extraction, some residual background can remain, which we account for by simultaneously fitting a linear component. Line measurements are then averaged over the two grism orients to produce final fluxes, central wavelengths, and line widths.

Comparing the uncertainty frames produced by the pipeline with the observed pixel-to-pixel fluctuations, we find that the pipeline uncertainty in general appears to significantly underestimate the true uncertainty. For this reason we instead estimate the noise level of a given spectrum by measuring the standard deviation of all pixels, excluding a wavelength window around the emission line set to be $2\sigma_{\text{line}}$ wide, after subtracting a linear background to make sure no trends will skew the measurements. We then use a Monte Carlo (MC) procedure to estimate the final flux errors: we randomize the spectra by drawing flux perturbations from a Gaussian with mean zero and standard deviation equal to the fluctuations measured in the original spectrum, and we refit the emission line. This is repeated 1000 times, and the final errors are calculated as the standard deviation of the resulting distribution of flux measurements.

We do not detect continuum emission in the spectra of any of our LAEs. Despite this, we can still place limits on the EWs. We do this simply by considering a maximum continuum level equal to the local noise measured from the spectrum. The resulting lower limits on the EWs are shown in Figure 5. We find rest-frame EWs above ~30 Å for all sources. This is expected from the resolution of NIRISS and is comparable to the commonplace definitions of LAEs used in the literature (e.g., L. Pentericci et al. 2011; D. P. Stark et al. 2011; L. Napolitano et al. 2024).

### 4.1.1. P50CLAE1

P50CLAE1 is shown in Figure 1. This is the clearest candidate of the sample and has an unambiguous line detection at 11710 Å, corresponding to a redshift of 8.63. The line is narrow, with a line width of only 30 Å (corresponding to ≈ 770 km s, noting that 1 NIRISS pixel corresponds to ≈1200 km s at this redshift), and is detected at more than $12\sigma$ significance despite the noticeable contamination from a nearby source in the GR150R grism image. The source is also visible, although faint, in both the F150W and F200W images.

### 4.1.2. P50CLAE2

This source, shown in Figure 2, is the faintest of our candidates, has a redshift of 9.01, and has a rest-frame EW of 36.9 Å. The line falls close to strong contaminants in both dispersion directions but fortunately remains in a region where it can be detected. The source is undetected in the imaging in F150W and F200W (see top right panels of Figure 2). Such nondetections are to be expected, since the exposure times of
<figure><img src="image_4.png" /><figcaption>Figure 3. Same as Figure 1, but for P50CLAE3.</figcaption></figure>

the direct images are only one-third of the dispersed exposure time, and no strong emission lines are expected in the wavelength regions covered by F150W and F200W when the Lyα line falls in the F115W filter. This highlights the importance of searching for line emission in the dispersed images, as well as direct imaging.

#### 4.1.3. P50CLAE3

In Figure 3 we show the highest-redshift candidate in the sample, which, at $z = 9.48$, lies at the very edge of the F115W bandpass. The line is very weak but nevertheless corresponds to a rest-frame EW of ≈45 Å. We also do not detect the source in the F150W or F200W direct images, for the same reasons as described for P50CLAE2 above.

#### 4.1.4. P52CLAE1

This is the lowest-redshift source, shown in Figure 4, and the most ambiguous of our detections. The reason is that the source lies close to the edge of the field, which means that the F150W and F200W spectra fall outside the edge of the dispersed image. Additionally, the line detected is significantly broader than the other Lyα detections, and the source has a cross-dispersion size that is larger than the PSF by an additional $\sigma = 1.4$ pixels. The source remains undetected in F150W and F200W direct imaging.

### 4.2. Completeness Estimation

In order to determine accurate number densities and LFs (see Section 4.3), we need a precise understanding of our sample completeness. We estimate the completeness using a source

injection and recovery simulation following E. C. Herenz et al. (2019). The simulation mimics the source selection procedure as closely as possible, but, requiring the simulations to be automated, we therefore cannot include the visual inspection steps of the procedure.

Since all our candidates all are either unresolved or marginally resolved, we simplify our simulations by limiting them to point sources. We model the sources as pure emission-line sources i.e., a given line flux, an observed wavelength of the line, and zero continuum. The emission line is modeled as a Gaussian line profile of PSF width (FWHM = 1.42 pixels) in both the dispersion and cross-dispersion directions (corresponding to 0.094 in the spatial direction and 1700 km s⁻¹ in the dispersion direction, assuming the average redshift of our sample). From the generated synthetic 1D spectrum we compute the inverse of the optimal extraction weighting along a trace calculated for a randomly chosen point in the direct frame. This produces a 2D spectrum, which we then add to the science frame. When new source locations are chosen, we make sure that no sources are placed within 15 pixels of each other in order to avoid artificial crowding.

The new seeded frame is then passed to source detection, where we do median filtering and source extraction using the same settings as for the actual source detection. For each source that is recovered in the GR150C dispersion direction we then extract the GR150R spectrum and do line fitting on the known location of the simulated line, as we would on the real detections. The noise in the spectrum is estimated from the spectral regions outside a $2\sigma$ wide window around the emission line and used as input to a 500-iteration MC simulation of the line fitting. This allows us to estimate the error on the line
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Same as Figure 1, but for P52CLAE1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** EW limits for the four LAE candidates as a function of redshift. All points are lower $1\sigma$ limits since the continuum is undetected.</figcaption></figure>

detection. We then consider a source rediscovered in both dispersion directions if the fitted line has a higher signal-to-noise ratio than $2\sigma$.

We run these simulations for luminosities between $10^{41.8}$ and $10^{43.5}$ erg s$^{-1}$ in steps of $\log(L) = 0.15$ and $7.4 \le z \le 9.5$ in steps of 0.16, and for each combination we generate three images of 150 sources each. This sampling allows us to construct a grid in L and z that we can sample to get the completeness for a given source. This simulation procedure is done for each field independently.

We present the results of the completeness simulations in Figure 6. We note a few interesting features. The first is that the completeness for all fields reaches a maximum around 80%-90%. This is because we have to mask a 200-pixel-wide region along the border of the dispersed images, in which the zeroth-order images of sources that lie outside the direct image are located. These mimic emission lines, but we cannot model and exclude them using the direct image (see Section 3.3). This is not quite the same as standard completeness; however, the reduction in effective volume is the same when used in conjunction with the $C^{-}$ and $1/V_{max}$ methods (see Section 4.3). The second noteworthy feature are the relative depths of the fields. We note that NGDEEP is the deepest field, as expected, but that the depths of Par50, Par52, and Par28 are comparable, and Par28 may indeed be the least complete of the four fields, despite having the longest exposure time of the PASSAGE parallels. This is due to unfortunate bright star contamination of that field.

### 4.3. Luminosity Functions

One of the most fundamental properties of a galaxy population is its LF. The LAE LF has been found to be approximately constant up to $z \sim 6$ (M. Ouchi et al. 2008; E. C. Herenz et al. 2019), above which a precipitous drop occurs that is often interpreted as the impact of an increasingly neutral IGM (N. Kashikawa et al. 2011). To date, all LAEs above $z = 7.3$ have been selected by their continuum emission (e.g., as Lyman break galaxies (LBGs)) and similar. Here we instead have a sample, albeit small, of line-selected Ly$\alpha$-emitting galaxies. It is therefore interesting to examine the
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Completeness levels for each field as a function of redshift and luminosity ($\log_{10}(\text{erg s}^{-1})$). Trends like the decrease of completeness close to the edges of the filter (corresponding to low and high redshifts) can be clearly seen.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Cumulative $\text{Ly}\alpha$ LF derived with the $C^-$ method (blue points) and $1/V_{\text{max}}$ method (orange points). The black lines (dotted and dotted-dashed) are the LFs determined by E. C. Herenz et al. (2019; H+19) between redshifts ~3 and ~6, with the dotted-dashed line using a completeness function that considers LAEs as point sources and the dotted line considering them as spatially extended.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Differential LF. Blue points indicate data from this work using the binned $1/V_{\text{max}}$ estimator. The dotted line indicates the global LAE LF from E. C. Herenz et al. (2019) derived for redshifts 3–6.7 (see Figure 7), and the solid line is the E. C. Herenz et al. (2019) LF rescaled to fit our measurements. The dotted-dashed line indicates the LF from I. G. B. Wold et al. (2022). Green diamonds and red squares indicate LF measurements from A. B. Drake et al. (2017b) and A. B. Drake et al. (2017a), respectively, at $2.91 < z < 6.64$.</figcaption></figure>

implied LF of these galaxies and how it differs from other selection techniques.

We determine the cumulative LF $\Phi$ of our sample using two different methodologies: the $1/V_{\text{max}}$ method and the $C^-$ method. The $1/V_{\text{max}}$ method is a refinement of a simple number count per total observed volume and calculates the number density by considering the maximum volume within which a specific source could be detected, given the completeness simulation presented in the previous section. $C^-$, on the other hand, calculates the cumulative LF using a product of Dirac $\delta$ functions. We refer the reader to R. Johnston (2011) and E. C. Herenz et al. (2019) for more detailed discussions of these methods. We present the two resulting LFs in Figure 7, together with an LF of $3 \lesssim z \lesssim 6.6$ LAEs from E. C. Herenz et al. (2019), noting that our observed densities are significantly lower than those at lower $z$.

The small number of sources prevents us from deriving a full differential LF, $\phi$, usually parameterized by the Schechter function (P. Schechter 1976). However, we can examine the relative normalization by assuming the same Schechter shape parameters as redshift 6 and scaling the function to fit our points. We derive $\phi$ in two luminosity bins, with the sample divided below and above $10^{42.5} \text{ erg s}^{-1}$, resulting in two galaxies in each bin. We then rescale the Schechter fit from E. C. Herenz et al. (2019) to optimally fit our points by recalculating the normalization $\phi_*$.

The results are presented in Figure 8, and the rescaled $\log_{10}(\phi_*)$ value is -3.64, as compared to -2.71 from E. C. Herenz et al. (2019). We can see that the rescaled $\phi_*$ value lies significantly below the E. C. Herenz et al. (2019) result by a factor of ~10, indicating a considerable reduction in the volume density of LAEs. This downward evolution in the
LF is expected for a Universe that is significantly more neutral at $z = 8$ than it is at $z = 6$, and only Ly$\alpha$ emission from galaxies residing in large ionized regions will be observable.

## 5. Implications for Reionization

Here we examine the implications that our observations and LFs have for the neutral content of the Universe, paying specific attention to Ly$\alpha$ EWs, transmission through the IGM, and the comparison of empirical and simulated LFs.

### 5.1. Intrinsic Equivalent Widths

The EWs of LAEs have been used in many studies to place constraints on the reionization process (see, e.g., M. Nakane et al. 2024). While we only have lower limits on the EW of our emitters, we can still use these to estimate the limits of the rest-frame EW of the galaxy before IGM absorption. As in M. J. Hayes & C. Scarlata (2023), we refer to this as the “emergent” EW. In order to calculate the emergent EW of a source, we have to make a few key assumptions. First, we assume that the Universe is largely neutral at these redshifts but that the galaxies lie in ionized regions. Second, we assume an intrinsic Ly$\alpha$ spectral profile, for which we use the median stack of COS observations of low-$z$ Ly$\alpha$-emitting galaxies from M. J. Hayes et al. (2021) since the constituent galaxies in the stack are at sufficiently low redshift that IGM absorption is negligible. This low-$z$ profile is observationally indistinguishable from spectral profiles observed out to at least $z = 5.5$. Third, we require a radius of the ionized bubble that the galaxy resides in, which we vary between 0.5 and 15 pMpc.

In order to calculate the IGM attenuation at the redshifts of our Galaxy for a given ionized bubble size, we follow the methodology of J. Witstok et al. (2025a) and make use of the code lymana_absorption.<fnref n="26" /> We calculate the IGM transmission curve as a function of wavelength, which we use to calculate the fraction of transmitted light from our assumed Ly$\alpha$ profile. We then multiply our rest-frame EWs by this scaling factor to get the emergent EW. The results are shown in Figure 9.

We note that only moderately large ionized bubbles are required for the emergent EWs to be in line with theoretical bounds and observed LAE EWs. Requiring the maximum emergent EW to be less than 250 Å would only require bubbles ≳2 pMpc. This is broadly consistent with the radii reported by M. J. Hayes & C. Scarlata (2023) and H. Umeda et al. (2024). However, 250 Å is not a hard limit but simply a typical value for a zero-age stellar population with a Salpeter initial mass function (IMF), and galaxies at these redshifts likely do not conform to these assumptions. For instance, there have been several suggestions that as we approach very low metallicities the IMF may become top-heavy (see, e.g., A. J. Cameron et al. 2024; D. Nandal et al. 2024; M. J. Hayes et al. 2025; A. Hutter et al. 2025, for some examples). One should also note that this calculation assumes an escape fraction of Ly$\alpha$ of 1. A more realistic $f_{esc}$ assumption would require larger bubble sizes to remain below a given maximum theoretical EW. However, the observed emergent rest-frame EWs themselves imply a high $f_{esc}$ in these galaxies (e.g., D. Sobral & J. Matthee 2019).

<figure><img src="image_10.png" /><figcaption>Figure 9. Implied emergent EWs for our sources as a function of the size of the ionized bubble the galaxy resides in.</figcaption></figure>

### 5.2. Ly$\alpha$ Transmission Factor

Following I. G. B. Wold et al. (2022), we now try to quantify the evolution of the Ly$\alpha$ IGM transmission using the Ly$\alpha$ transmission factor, which is defined as

$$T_{Ly\alpha} = \frac{\rho_{Ly\alpha}^{z2}/\rho_{Ly\alpha}^{z1}}{\rho_{UV}^{z2}/\rho_{UV}^{z1}},$$

where $\rho_{Ly\alpha}^{z2}$ is the Ly$\alpha$ LD at a given redshift $z$ and $\rho_{UV}^{z2}$ is the UV LD. The luminosity densities are calculated as $\int L\phi(L)dL$ using the Schechter function fit of $\phi$, and we consider integration limits between $M_{UV} = -22$ and $-18$ mag and between $\log(L \text{ erg s}^{-1}) = 42.4$ and 45 for the UV LF and the Ly$\alpha$ LFs, respectively.

We compute the Ly$\alpha$ density at redshifts 7.5–9.5 from the rescaled Schechter function presented in Section 4.3, and we compute the $z \sim 6$ Ly$\alpha$ density from the original E. C. Herenz et al. (2019) fit. For the UV LF we use the redshift 5.9 function from R. J. Bouwens et al. (2021) as the low-redshift anchor point and R. Bouwens et al. (2023) for the high-$z$ measurement. We find that the implied $T_{Ly\alpha}$ is 0.27, i.e., the transmission at redshifts 7.5–9.5, dominated by the $z = 8.5–9.5$ end, is a quarter of that at redshift 6. This implies a significant evolution in the IGM neutral fraction between these redshifts, but not as strong as would be expected if reionization were late and rapid (e.g., R. P. Naidu et al. 2020). We note that a reduction in Ly$\alpha$ transmission of a factor of four is the same as the factor-of-four drop in Ly$\alpha$ fraction that M. Tang et al. (2024) estimate between redshifts 5.5 and 8 using Ly$\alpha$ follow-up of continuum-selected galaxies.

One should note two things about this method, however: The first is that a transmission factor does not uniquely correspond to a neutral fraction, since it would depend on the topology of the ionized regions in the IGM. Second, this method of calculating the transmission may be oversimplified by assuming a nonevolving EW distribution (although the elapsed

<footnote n="26">²⁶ https://github.com/joriswitstok/lymana_absorption</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Comparison of the differential LFs derived in this work (green circles) with those from various simulations. Left panel: simulations are taken from the ASTRAEUS simulation of A. Hutter et al. (2023; H+23), computed for a Salpeter IMF at various redshifts. Right panel: simulations are taken from A. M. Morales et al. (2021; M+21) using dotted lines and the analytical fit to the LF from S. Bruton et al. (2023a; B+23) using dotted-dashed lines. We also include versions of our LFs where we only consider P50CLAE1, the highest-significance source, which is shown as the black circle and solid line.</figcaption></figure>

cosmic time over this interval is only 350 Myr, and perhaps large changes in the emergent EW distribution should not be expected).

### 5.3. Comparisons to Expected Ly$\alpha$ LFs

Another way to use our findings to make inferences on the reionization history is to directly compare them to theoretical predictions. We compare our differential LFs with theoretical predictions in Figure 10.

In the left panel, we contrast the observed LFs with those of A. Hutter et al. (2023, 2025), computed using the ASTRAEUS simulation framework (A. Hutter et al. 2021), which includes both a semianalytic galaxy formation model and a seminumeric reionization simulation, assuming a Salpeter IMF. From the figure it is clear that our observed LF would best correspond to the simulated LF at redshift 7.3 and is quite discrepant from the expected level at our mean redshift of 8.6. One should note that there are some potential caveats in the simulations, such as the assumption that the emergent Ly$\alpha$ is a Gaussian profile with no velocity offset, which lower-redshift observations show to essentially never be the case (M. J. Hayes et al. 2023). A. Hutter et al. (2023) do show that assuming double-peaked Ly$\alpha$ makes only marginal difference to the LF. However, there is a possibility that more red-peak-dominated emergent line profiles, similar to those observed at lower z, would reduce the IGM sensitivity of the Ly$\alpha$ and boost the observed LF. The observed LF would then match simulated LFs at higher redshift, although we unfortunately cannot precisely estimate the redshift of best agreement.

A. Hutter et al. (2025) also present ASTRAEUS simulations using an evolving IMF, where dense gas produces a top-heavy IMF and enhances the number of massive-ionizing stars. Indeed, the IMF is not expected to be constant as we reach early cosmic times, when the ISM of star-forming becomes denser, although the exact nature of the evolution is hard to predict (e.g., A. J. Cameron et al. 2024). A consequence of such an IMF would be that lower-mass galaxies produce more ionizing and Ly$\alpha$ photons, resulting in different reionization

topology and Ly$\alpha$ LF (A. Hutter et al. 2025). We also compared our data with the evolving IMF simulations and found that our observations are well matched with the simulated LF at $z \simeq 8$. However, one should note here that the impact of the evolving IMF is relatively small at redshifts below 9 and becomes much more pronounced thereafter, which may also help in explaining the detection of a $z \sim 9.5$ source in our sample.

We now turn to the right panel of Figure 10, where we plot the LF predictions of S. Bruton et al. (2023a) at $z=7$ for three different values of $x_{H1}$ and LFs for $7 \le z \le 10$ for the simulated average $x_{H1}$ at those redshifts from A. M. Morales et al. (2021). We note that our data points lie in the upper span of the models presented, matching closely the line of S. Bruton et al. (2023a) for a neutral fraction of 0.5. Considering the A. M. Morales et al. (2021) LFs, our data lie between the $z=7$ and $z=8$ data points. Again, our points fall above the simulated LFs but are marginally consistent with the $z=8$ expectations at $1\sigma$.

### 5.4. Uncertainties on the Luminosity Function

There are three main contributors to the uncertainty on our LF determinations: Poisson/counting statistics in the low-number regime, the unknown interloper/contamination fraction, and CV. In this subsection we address each of these in turn.

We estimate the impact of counting statistics by making some simplifying assumptions: first, that the luminosity distribution within a LF bin is uniform; second, that objects do not scatter between bins; and third, that the bin widths are static. In order to estimate the counting errors, shown as error bars in Figures 8 and 10, we resample the number of sources in the bin from a Poisson distribution of mean 2. For each source we then assign a luminosity from a uniform distribution in log (L) with edges set by the bin edges. These new sources are then used to calculate $\phi$ using the $1/V_{max}$ methodology as before. We repeat this procedure 1000 times and show the 16th and 84th percentiles as the asymmetric error bars in the figures.
<table><caption>**Table 2**
Observed Lyα Characteristics</caption><thead><tr><th>Source<br>(1)</th><th>R.A.<br>(2)</th><th>Decl.<br>(3)</th><th>Redshift<br>(4)</th><th>Luminosity<sup>a</sup><br>(5)</th><th>Flux<sup>b</sup><br>(6)</th><th>Width<sup>c,d</sup><br>(7)</th><th>EW<sup>d,e</sup><br>(8)</th></tr></thead><tbody><tr><td>P50CLAE1</td><td>189.18718</td><td>62.05699</td><td>8.63 ± 0.004</td><td>3.4 ± 0.4</td><td>36.9 ± 4.4</td><td>30 ± 6</td><td>70.1</td></tr><tr><td>P50CLAE2</td><td>189.15394</td><td>62.06494</td><td>9.01 ± 0.01</td><td>1.5 ± 0.3</td><td>15.3 ± 3</td><td>52 ± 8</td><td>36.9</td></tr><tr><td>P50CLAE3</td><td>189.16159</td><td>62.07098</td><td>9.48 ± 0.02</td><td>2.8 ± 0.45</td><td>24.7 ± 4</td><td>47 ± 11</td><td>45.1</td></tr><tr><td>P52CLAE1</td><td>150.17060</td><td>2.03315</td><td>7.60 ± 0.01</td><td>8.4 ± 0.8</td><td>122.4 ± 12</td><td>88.5 ± 3</td><td>98.6</td></tr></tbody></table>

**Notes.**

a In units of $10^{42}$ erg s.

b In units of $10^{-19}$ erg s cm$^2$.

c $1\sigma$ line width.

d In units of Å.

e Rest-frame $1\sigma$ lower limit.

Estimating the impact of contamination is nontrivial, but to give a conservative measure, we omit all sources except for P50CLAE1, which is the highest-confidence detection, and recompute the LFs. This assumption is very conservative but useful for establishing limits of the implied LFs. The result is shown as the black circle and solid line in Figure 10. To compute the differential LF, we have to make a choice of the bin width, which we set to 0.5. The resulting point is the outcome of one $z = 8.6$ galaxy and now matches with the $z = 8$ curve from A. M. Morales et al. (2021), significantly reducing the tension with simulations. When comparing to the ASTRAEUS simulations of A. Hutter et al. (2023), this LF is consistent with expectations between $z \sim 8$ and 9.

CV describes the scatter in the measurements that is intrinsic to the Universe, encapsulating intrinsic dark matter halo variation in the volume, the bias, and the hydrogen-reionization process itself that determines the Lyα visibility. Pure CV uncertainties can be calculated from the dark matter density field and are largest in small contiguous fields, but both the large $\Delta z$ of our Lyα survey and the four independent pointings do work to reduce CV. Using the CV calculator of M. Trenti & M. Stiavelli (2008), we calculate an additional CV error budget of only 14% assuming dark matter halo masses above $10^{10} M_{\odot}$. In order to include the effects of reionization/IGM topology on the LAE counts, we need to turn to a simulation. We construct coarse light cones from the simulation of S. Bruton et al. (2023a), using the average reionization history updated for the JWST Cycle 1 LF in C. A. Mason et al. (2025, their Figure 9). We then select random locations in the cube, extract regions corresponding to the NIRISS field of view and redshift range $z = 7.5-9.5$, and count the number of LBGs with $M_{UV} < -19$ and LAEs with $L_{Ly\alpha} > 1.5 \times 10^{42}$ erg s$^{-1}$ (the least luminous object in our sample and an exact $5\sigma$ detection; Table 2). We repeat this for four pointings to mimic our survey, and we execute a 1000-realization MC simulation to predict the number of LBGs and LAEs. At the 16th and 84th percentiles, we expect to observe 4–10 LBGs and 0–3 LAEs. With four LAEs in our survey, we slightly exceed the expected values, as would also be predicted from our LFs. We develop these arguments and compare our result with the ASTRAEUS simulations in Sections 6.1 and 6.3.

# 6. Discussion

## 6.1. Implications of Recent Lyα Observations

Pre-JWST results suggested that finding significant numbers of Lyα-emitting galaxies at redshifts beyond 7 was going to be

very challenging with only a few detections reported (M. A. Schenker et al. 2012, 2014; A. Hoag et al. 2019; V. Tilvi et al. 2020; I. Jung et al. 2022). Large ground-based narrowband surveys, such as UltraVISTA, targeted $z \sim 8$ LAEs but found none, due to the difficulty of reaching the required flux limits against the bright atmospheric background in the NIR: the limiting flux of UltraVISTA was about $1.5 \times 10^{-17}$ erg s$^{-1}$ cm$^{-2}$, while our $z > 8$ detections have an average flux that is $\sim 6$ times fainter ($2.4 \times 10^{-18}$ erg s$^{-1}$ cm$^{-2}$). However, the situation has drastically changed after the launch of JWST, and detections now reach deep into the EoR. Detecting sources like GNz11 (A. J. Bunker et al. 2023), and now an LAE at $z = 13$ (J. Witstok et al. 2025a), was completely unexpected. M. Tang et al. (2024) used observations of the EGS, GOODS-N, GOODS-S, and A2744, covering 41 NIRSpec pointings (corresponding to total area coverage $\gtrsim 440$ arcmin$^2$), to detect 33 galaxies with Lyα emission above redshift 6.5, as well as confirm a $7 \le z \le 8$ galaxy overdensity in the EGS field. In this work we detect three LAEs with a similar $\Delta z$ to galaxies in the EGS overdensity, within a single NIRISS pointing of only $\sim 4.8$ arcmin$^2$.

The question that arises is whether the unexpected prevalence of Lyα at high redshift in fact has implications for models of the EoR. In a reionization scenario that is mainly powered by a few bright galaxies (as suggested by, e.g., R. P. Naidu et al. 2022, sometimes referred to as oligarchic reionization) we would expect to detect LAEs in a relatively clustered manner, as indeed we may be doing now, but we would also expect such sight lines to be relatively rare and that the LAEs detected should be dominated by bright sources. Given the total area currently surveyed, it appears that LAEs are more common than expected, and the LAEs we detect here are also unusually bright. On the other hand, Y. Qin & J. S. B. Wyithe (2025) suggest that a more distributed reionization process, where fainter galaxies play a more important role, will boost the likelihood of detecting fainter LAEs. This scenario is more likely to be in line with our findings.

A. Hutter et al. (2023) show that the LAEs in the ASTRAEUS simulation are dominated by more massive galaxies residing in larger ionized regions. We therefore conduct a test to see what the likelihood is of observing a conjunction of LAEs like the one we report here in a volume corresponding to a single NIRISS pointing. To do this, we analyze the population of LAEs in a $z = 9$ snapshot of the simulation: we select LAEs from the Salpeter IMF simulation with luminosities above $1.5 \times 10^{42}$ erg s—approximately the
lowest luminosity in our $8.5 \le z \le 9.5$ sample. Then, we bin the simulation volume in $120^3$ cells to obtain LAE number counts in each cell. For each cell we then count the numbers of LAEs within a volume corresponding to the $V_{max}$ volume estimate for Par50 centered on that cell.

We find that only ≈0.21% of such volumes have ≥3 sources in them. This suggests that the reionization history and topology produced in the Salpeter IMF version of ASTRAEUS are an unlikely representation of the $z=8.6$ Universe, when confronted with our NIRISS observations (under the assumption that all our sources are real). To investigate this further, we conduct the same experiment for the simulations with an evolving IMF. The resulting likelihood of ≈0.74% is three times higher than for the Salpeter IMF but nevertheless remains unlikely.

To draw fully quantitative conclusions regarding the implications of our observations for EoR morphology would require detailed modeling that is beyond the scope of this work, as well as follow-up of the sources. Nevertheless, these results, especially when taken in conjunction with other LAE observations, suggest that the picture of LAEs being very rare deep into the EoR may be incorrect. This may imply distributed reionization topologies and a boost in the observability of more common, fainter LAEs.

## 6.2. Field-to-field Variance

Early simulations by H. Jensen et al. (2013) assessed the possibility of using the field-to-field variation in Ly$\alpha$-emitting galaxies as a discriminator between reionization. Indeed, it has been found that at least one of the JWST deep fields, EGS, does cover a high-redshift overdensity (V. Tilvi et al. 2020; I. Jung et al. 2022; M. Tang et al. 2024). Another indication of significant field-to-field variance comes from L. Napolitano et al. (2024), who measured the Ly$\alpha$ fraction in the CEERS and JADES surveys using MSA spectroscopic data. Their results indicate large differences (up to a factor of 3) between the fields, depending on whether or not known overdensities are included. The results also appear to cast some doubt over the previously established increase of $X_{Ly\alpha}$ up to redshift 6 and subsequent drop, instead showing a large scatter of values across all redshifts. One should note, however, that the small slitlets of the NIRSpec MSA may lead to large slit losses of Ly$\alpha$, which is typically more extended than the UV emission (M. Hayes et al. 2013; L. Wisotzki et al. 2016; F. Leclercq et al. 2020; A. Runnholm et al. 2023).

The principal way to overcome and characterize CV is to observe many spatially uncorrelated fields (see, e.g., S. Bruton et al. 2023a, for Ly$\alpha$ studies), which indeed was one of the primary motivations for the PASSAGE project. The number counts of LAEs in our fields are quite uneven, with zero LAEs recovered in the deepest NGDEEP data but three found in one of the shallower parallel fields. The sample dispersion (≡$\frac{\text{Sample Mean}}{\text{Sample Variance}}$) is relatively large (1.5) compared to the expectation (1) for a Poisson process. We therefore test whether this is consistent with being simply a random draw from a uniform Poisson distribution or whether additional CV is required to explain the data.

To do this, we draw samples of four random variates from a Poisson distribution with the same mean as our observed sample (=1) and calculate the dispersion. We repeat this 10,000 times and then compare our observed dispersion to the distribution. We find that the probability of finding a dispersion

≥1.5 is 0.13, i.e., unlikely but not impossible. From this we conclude that while some impact of CV is probable, we do not require it to explain our observations. Characterizing CV and the deviations from simple Poissonian statistics in more detail requires more independent field observations. The Cycle 3 program #3383 (PI: Glazebrook) obtained pure parallel NIRISS data similar to PASSAGE; however, only two fields in this program have the required depth and F115W grism coverage to be usable for this purpose.

## 6.3. Improving $x_{HI}$ Precision with JWST/NIRISS

Obviously more data are required to make more precise estimates of the neutral fraction from Ly$\alpha$ LFs, and here we address the question of how many pointings would be needed. We noted in Section 5.3 that there are three sources of variation on the number counts: CV (including the overlaid state of the IGM/H II bubble size), interlopers/contamination, and Poisson uncertainty. Here we assess the Poisson errors directly in our data (e.g., assuming that the LF values are accurate but their error bars are larger than desired) and the impact of CV +reionization on Ly$\alpha$ visibility using updates to the simulations of S. Bruton et al. (2023a).

The right panel of Figure 10 shows that the high-luminosity green data point ($\log(L_{Ly\alpha}/\text{erg s}^{-1}) = 42.72$) provides rather strong constraints on the $x_{HI}$, because the simulated LFs of S. Bruton et al. (2023a) diverge at high L (see also W. Hu et al. 2019 and I. G. B. Wold et al. 2022). Taking Bruton et al.'s LFs as truth and assuming purely Poissonian errors, we estimate that reducing the error on $x_{HI}$ to $\Delta x_{HI} = 0.1$ would require only eight fields, or approximately doubling the survey compared to that presented here. However, the same simulations can also be used to address the CV error budget as in Section 5.4. With only eight fields, the inference would remain dominated by CV, with an estimated uncertainty of ≈70% on the LF data point, due mainly to the reionization topology.

However, the high-$L$ data point is strongly dependent on the luminosity of a single galaxy, P52CLAE1. This galaxy is ≈3 times more luminous than the other three sources and has by far the highest EW (see Table 2), making it the least representative object in our sample. The recovery of that galaxy may have been an improbable occurrence (as suggested by the LF alone and field-to-field variations at high luminosity; I. G. B. Wold et al. 2022), so we repeat the same calculation for the more common low-$L_{Ly\alpha}$ galaxies in the bin near $\log(L_{Ly\alpha}/\text{erg s}^{-1}) = 42.3$ (lower green data point). In this case the uncertainty of ±0.1 on $x_{HI}$ can be reached with ≈40 fields. With 40 random fields, the CV+reionization budget would add just 0.11 dex to the $\Delta\phi$ uncertainty in the LF and would become negligible (again following S. Bruton et al. 2023a). Such a survey would be a 10-fold increase in size compared to this study.

Finally, we address the issue of possible contamination in the same way, using the LF data point derived under the assumption that only one source is real (black point in the right panel of Figure 10). As an upper limit, only the higher error bound is finite, so we simply estimate the number of sources (of constant luminosity) required to bring this upper bound down to the level of $\Delta x_{HI} = 0.1$. To reach this level of precision, a similar NIRISS campaign would require ≈70 random, uncorrelated pointings. The CV contribution to the uncertainty on $\phi$ would again be negligible (0.08 dex).
## 7. Summary and Conclusions

We have conducted what is (to our knowledge) the first blind, unbiased survey for Lyα-emitting galaxies at $z > 7.5$ with JWST. We have conducted a WFSS survey using the NIRISS instrument, obtained under the PASSAGE program. We focus on three fields for which spectroscopic data have been obtained for at least 2 hr in both the row and column dispersion directions using the F115W filter, for which spectra in the F150W and F200W filters have also been obtained, and for which deep imaging exists at wavelengths blueward of the WFSS data. We further included data from the NGDEEP program in the Hubble Ultra Deep Field, which has obtained similar NIRISS data. Here we summarize the results presented in this work, as well as the main conclusions we draw.

1. We detect Lyα line emission from four LAEs in the redshift range $7.5 \le z \le 9.5$ across four surveyed NIRISS pointings. These are the highest-redshift Lyα-selected galaxies known to date. Their luminosities fall in the range $(1.5-8.4) \times 10^{42} \text{ erg s}^{-1}$.

2. We determine the lower limit of the rest-frame Lyα EWs to be 36–99 Å and investigate what these EWs imply for the emergent EW given a range of ionized bubble sizes. We conclude that relatively small bubbles of radius $\approx 2$ pMpc would yield physically probable emergent EWs.

3. We construct the cumulative and differential Lyα LFs and contrast them with the low-redshift results of E. C. Herenz et al. (2019). We find that the LF undergoes significant evolution between redshift 6 and our mean redshift of 8.6. However, this amounts to a factor of $\lesssim 10$, and contrasting this observed LF with theoretical expectations from A. Hutter et al. (2023), A. M. Morales et al. (2021), and S. Bruton et al. (2023a) shows the evolution to be significantly less than expected.

4. Contrasting the evolution of the Lyα LF to the UV LF from R. Bouwens et al. (2023) indicates that the Lyα LF undergoes significantly more change than the UV LF—demonstrating that it is unlikely that galaxy evolution is the main cause of Lyα LF change and that it is rather due to an increasingly neutral IGM.

5. One potential risk in making inferences on the average evolution of a galaxy population based on a small number of fields is the impact of CV on the results. We find large variations between our observed fields that may indicate that CV plays a role; however, having only four fields is not sufficient for us to conclude that the distribution is non-Poissonian.

6. We also estimate the increase in the number of observed JWST/NIRISS fields required to reach a precision of $\Delta x_{H\&I} = 0.1$. We find that this number ranges from $\sim 8$ to $\sim 70$ fields depending on whether we can reach the more constraining higher-luminosity sources, as well as the fraction of contamination in our current sample. This would be possible to achieve with future PASSAGE-like pure parallel JWST/NIRISS surveys.

## Acknowledgments

This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the

Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with programs #1571 and #2079. The specific observations analyzed can be accessed via doi: 10.17909/k1a0-qg30.

M.J.H. is supported by the Swedish Research Council (Vetenskapsrådet) and is Fellow of the Knut & Alice Wallenberg Foundation. M.B. acknowledges support from the ERC grant FIRSTLIGHT and Slovenian national research agency ARIS through grants N1-0238 and P1-0188. A.A. is supported by European Union–NextGenerationEU RFF M4C2 1.1 PRIN 2022 project 2022ZSL4BL INSIGHT. A.S.-L. acknowledges support from Knut and Alice Wallenberg Foundation. H.A. acknowledges support from CNES, focused on the JWST mission, and the Programme National Cosmology and Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, cofunded by CEA and CNES. A.H. acknowledges support by the VILLUM FONDEN under grant 37459. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. We thank the anonymous referee, whose thoughtful comments have enhanced the quality of the manuscript.

## ORCID iDs

Axel Runnholm https://orcid.org/0000-0002-1025-7569
Matthew J. Hayes https://orcid.org/0000-0001-8587-218X
Vihang Mehta https://orcid.org/0000-0001-7166-6035
Matthew A. Malkan https://orcid.org/0000-0001-6919-1237
Claudia Scarlata https://orcid.org/0000-0002-9136-8876
Kalina V. Nedkova https://orcid.org/0000-0001-5294-8002
Marc Rafelski https://orcid.org/0000-0002-9946-4731
Benedetta Vulcani https://orcid.org/0000-0003-0980-1499
Mason Huberty https://orcid.org/0009-0002-9932-4461
E. Christian Herenz https://orcid.org/0000-0002-8505-4678
Sean Bruton https://orcid.org/0000-0002-6503-5218
Ayan Acharyya https://orcid.org/0000-0003-4804-7142
Hakim Atek https://orcid.org/0000-0002-7570-0824
Ivano Baronchelli https://orcid.org/0000-0003-0556-2929
Andrew J. Battisti https://orcid.org/0000-0003-4569-2285
Maruša Bradač https://orcid.org/0000-0001-5984-0395
Andrew J. Bunker https://orcid.org/0000-0002-8651-9879
Y. Sophia Dai https://orcid.org/0000-0002-7928-416X
Clea Hannahs https://orcid.org/0009-0009-0868-8165
Farhanul Hasan https://orcid.org/0000-0002-0072-0281
Keunho J. Kim https://orcid.org/0000-0001-6505-0293
Nichla Leethochawalit https://orcid.org/0000-0003-4570-3159
Yu-Heng Lin https://orcid.org/0000-0001-8792-3091
Michael J. Rutkowski https://orcid.org/0000-0001-7016-5220
Alberto Saldana-Lopez https://orcid.org/0000-0001-8419-3062
Zahra Sattari https://orcid.org/0000-0002-0364-1159
Xin Wang https://orcid.org/0000-0002-9373-3865

## References

Atek, H., Labbé, I., Furtak, L. J., et al. 2024, Natur, 626, 975
Atek, H., Malkan, M., McCarthy, P., et al. 2010, ApJ, 723, 104
Bagley, M. B., Pirzkal, N., Finkelstein, S. L., et al. 2024, ApJL, 965, L6
Bagley, M. B., Scarlata, C., Henry, A., et al. 2017, ApJ, 837, 11
Barbary, K. 2016, JOSS, 1, 58
Battisti, A. J., Bagley, M. B., Rafelski, M., et al. 2024, MNRAS, 530, 894
Becker, G. D., Bolton, J. S., Zhu, Y., & Hashemi, S. 2024, MNRAS, 533, 1525
Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
Bolan, P., Lemaux, B. C., Mason, C., et al. 2022, MNRAS, 517, 3263
Bosman, S. E. I., Davies, F. B., Becker, G. D., et al. 2022, MNRAS, 514, 55
Bouwens, R., Illingworth, G., Oesch, P., et al. 2023, MNRAS, 523, 1009

Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47

Brammer, G., 2019 Grizli: Grism Redshift and Line Analysis Software, Astrophysics Source Code Library, ascl:1905.001

Brammer, G. 2023, Grizli, v1.9.11, Zenodo, doi:10.5281/zenodo.8370018

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

Bruton, S., Scarlata, C., Haardt, F., et al. 2023b, ApJ, 953, 29

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

Cameron, A. J., Katz, H., Witten, C., et al. 2024, MNRAS, 534, 523

Cuby, J. G., Hibon, P., Lidman, C., et al. 2007, A&A, 461, 911

Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W. 2009, IEEEP, 97, 1482

Dijkstra, M. 2014, PASA, 31, e040

Dijkstra, M., & Wyithe, J. S. B. 2012, MNRAS, 419, 3181

Drake, A. B., Garel, T., Wisotzki, L., et al. 2017a, A&A, 608, A6

Drake, A. B., Guiderdoni, B., Blaizot, J., et al. 2017b, MNRAS, 471, 267

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

Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13

Finkelstein, S. L., D'Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36

Gonzaga, S., Hack, W., Fruchter, A., & Mack, J. 2012, The DrizzlePac Handbook (Baltimore, MD: STScI)

Greig, B., Mesinger, A., Bañados, E., et al. 2024, MNRAS, 530, 3208

Gronke, M., Dijkstra, M., Trenti, M., & Wyithe, S. 2015, MNRAS, 449, 1284

Gwyn, S. 2020, in ASP Conf. Ser. 527, Astronomical Data Analysis Software and Systems XXIX, ed. R. Pizzo et al. (San Francisco, CA: ASP) 575

Hayes, M., Östlin, G., Schäerer, D., et al. 2013, ApJ, 765, L27

Hayes, M., Schäerer, D., Östlin, G., et al. 2011, ApJ, 730, 8

Hayes, M. J., Runnholm, A., Gronke, M., & Scarlata, C. 2021, ApJ, 908, 36

Hayes, M. J., Runnholm, A., Scarlata, C., Gronke, M., & Rivera-Thorsen, T. E. 2023, MNRAS, 520, 5903

Hayes, M. J., Saldana-Lopez, A., Citro, A., et al. 2025, ApJ, 982, 14

Hayes, M. J., & Scarlata, C. 2023, ApJL, 954, L14

Herenz, E. C., Wisotzki, L., Saust, R., et al. 2019, A&A, 621, A107

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

Hoffmann, S. L., Mack, J., Avila, R., et al. 2021, AAS Meeting, 238, 216.02

Hu, W., Wang, J., Zheng, Z.-Y., et al. 2019, ApJ, 886, 90

Hutter, A., Cueto, E. R., Dayal, P., et al. 2025, A&A, 694, A254

Hutter, A., Dayal, P., Yepes, G., et al. 2021, MNRAS, 503, 3698

Hutter, A., Trebitsch, M., Dayal, P., et al. 2023, MNRAS, 524, 6124

Itoh, R., Ouchi, M., Zhang, H., et al. 2018, ApJ, 867, 46

Jensen, H., Laursen, P., Mellema, G., et al. 2013, MNRAS, 428, 1366

Johnston, R. 2011, A&ARv, 19, 41

Jung, I., Finkelstein, S. L., Larson, R. L., et al. 2022, arXiv:2212.09850

Kashikawa, N., Shimasaku, K., Matsuda, Y., et al. 2011, ApJ, 734, 119

Keating, L. C., Weinberger, L. H., Kulkarni, G., et al. 2020, MNRAS, 491, 1736

Kikuta, S., Ouchi, M., Shibuya, T., et al. 2023, ApJS, 268, 24

Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196

Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14) (Trieste: SISSA), 1

Kulkarni, G., Keating, L. C., Haehnelt, M. G., et al. 2019, MNRAS, 485, L24

Kusakabe, H., Verhamme, A., Blaizot, J., et al. 2022, A&A, 660, A44

Laursen, P., Sommer-Larsen, J., Milvang-Jensen, B., Fynbo, J. P. U., & Razoumov, A. O. 2019, A&A, 627, A84

Leclercq, F., Bacon, R., Verhamme, A., et al. 2020, A&A, 635, A82

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

Mason, C. A., Fontana, A., Treu, T., et al. 2019, MNRAS, 485, 3947

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

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

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

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

Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109

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

Nandal, D., Sibony, Y., & Tsiatsiou, S. 2024, A&A, 688, A142

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

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

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

Ouchi, M., Harikane, Y., Shibuya, T., et al. 2018, PASJ, 70, S13

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

Ouchi, M., Shimasaku, K., Akiyama, M., et al. 2008, ApJS, 176, 301

Pentericci, L., Fontana, A., Vanzella, E., et al. 2011, ApJ, 743, 132

Pirzkal, N., Rothberg, B., Papovich, C., et al. 2024, ApJ, 969, 90

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

Qin, Y., Mesinger, A., Prelogović, D., et al. 2024, arXiv: 2412.00799

Qin, Y., & Wyithe, J. S. B. 2025, MNRAS Lett., 538, L16

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

Runnholm, A., Hayes, M. J., Lin, Y.-H., et al. 2023, MNRAS, 522, 4275

Schechter, P. 1976, ApJ, 203, 297

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

Schenker, M. A., Stark, D. P., Ellis, R. S., et al. 2012, ApJ, 744, 179

Scoville, N., Abraham, R. G., Aussel, H., et al. 2007, ApJS, 172, 38

Sharma, M., Theuns, T., Frenk, C., et al. 2016, MNRAS, 458, L94

Sobral, D., & Matthee, J. 2019, A&A, 623, A157

Stark, D. P., Ellis, R. S., Chiu, K., Ouchi, M., & Bunker, A. 2010, MNRAS, 408, 1628

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

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

Thai, T. T., Tuan-Anh, P., Pello, R., et al. 2023, A&A, 678, A139

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

Trenti, M., & Stiavelli, M. 2008, ApJ, 676, 767

Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110

Umeda, H., Ouchi, M., Kikuta, S., et al. 2025, ApJS, 277, 37

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

Willis, J. P., Courbin, F., Kneib, J. P., & Minniti, D. 2008, MNRAS, 384, 1039

Wisotzki, L., Bacon, R., Blaizot, J., et al. 2016, A&A, 587, A98

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

Witstok, J., Maiolino, R., Smit, R., et al. 2025b, MNRAS, 536, 27

Witstok, J., Smit, R., Saxena, A., et al. 2024, A&A, 682, A40

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

Zhu, Y., Becker, G. D., Bosman, S. E. I., et al. 2024, MNRAS, 533, L49