A GLIMPSE into the very faint end of the Hβ+[O III]λλ4960,5008 luminosity function at z ∼ 7 – 9 behind Abell S1063

Damien Korber; Iryna Chemerynska; Lukas J. Furtak; Hakim Atek; Ryan Endsley; Daniel Schaerer; John Chisholm; Vasily Kokorev; Alberto Saldana-Lopez; Angela Adamo; Julian B. Muñoz; Pascal A. Oesch; Romain Meyer; Rui Marques-Chaves; Seiji Fujimoto

doi:10.1051/0004-6361/202556868

Home

All issues

Volume 708 (April 2026)

A&A, 708 (2026) A43

Full HTML

Open Access

Issue		A&A Volume 708, April 2026


Article Number		A43
Number of page(s)		24
Section		Extragalactic astronomy
DOI		https://doi.org/10.1051/0004-6361/202556868
Published online		31 March 2026

A&A, 708, A43 (2026)

A GLIMPSE into the very faint end of the Hβ+[O III]λλ4960,5008 luminosity function at z ∼ 7 – 9 behind Abell S1063

Damien Korber¹^★, Iryna Chemerynska³, Lukas J. Furtak⁴, Hakim Atek³, Ryan Endsley⁵, Daniel Schaerer¹^,2, John Chisholm⁵, Vasily Kokorev⁵, Alberto Saldana-Lopez⁶, Angela Adamo⁶, Julian B. Muñoz⁵, Pascal A. Oesch¹^,9, Romain Meyer¹, Rui Marques-Chaves¹ and Seiji Fujimoto⁷^,8

¹ Observatoire de Genève, Université de Genève, Chemin Pegasi 51, 1290, Versoix, Switzerland
² CNRS, IRAP, 14 Avenue E. Belin, 31400, Toulouse, France
³ Institut d’Astrophysique de Paris, CNRS, Sorbonne Université, 98bis Boulevard Arago, 75014, Paris, France
⁴ Department of Physics, Ben-Gurion University of the Negev, P.O. Box 653, Be’er-Sheva, 84105, Israel
⁵ Department of Astronomy, The University of Texas at Austin, Austin, TX, 78712, USA
⁶ Department of Astronomy, Oskar Klein Centre, Stockholm University, AlbaNova University Center, SE-106 91, Stockholm, Sweden
⁷ David A. Dunlap Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada
⁸ Dunlap Institute for Astronomy and Astrophysics, 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada
⁹ Cosmic Dawn Center (DAWN), Denmark, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, København N, DK-2200, Denmark

^★ Corresponding author: This email address is being protected from spambots. You need JavaScript enabled to view it.

Received: 14 August 2025
Accepted: 10 February 2026

Abstract

We used the ultra-deep GLIMPSE JWST/NIRCam survey to constrain the faint end of the [O III]+Hβ luminosity function (LF) down to 10³⁹erg s⁻¹ at z ∼ 7 − 9 behind the lensed Hubble Frontier Field galaxy cluster Abell S1063. We applied a spectral energy distribution fitting on a Lyman-break galaxy selected sample of 164 lensed galaxies and measured their combined Hβ+[O III]λλ4960, 5008 flux to build the emission line LF. We found a [O III]+Hβ LF with a faint-end slope (α = −1.78^+0.06_−0.06 for z = 7 − 8 and α = −1.55^+0.11_−0.11 for z = 8 − 9), which is flatter than the UV LF at similar redshifts (α ≤ −2) and suggests a lower number density of weak [O III]+Hβ emitting galaxies at fixed M_UV. We analysed several possible explanations: (i) a decrease in the [O III]+Hβ-to-UV ratio due to bursty star formation histories (SFHs), (ii) the effect of metallicity on the [O III]-to-Hβ ratio, or (iii) signs of a faint-end turnover in the UV LF. Under the assumption of an evolving [O III]-to-Hβ ratio, we separated the contribution of [O III]λ5008 and Hβ and obtained a flatter [O III]λ5008 LF (α = −1.66^+0.05_−0.05 for z = 7 − 8 and α = −1.45^+0.09_−0.10 for z = 8 − 9) but steeper Hβ LF ( $α = - 1 . 95_{- 0.08}^{+ 0.08}$ $Mathematical equation: $ \alpha=-1.95_{-0.08}^{+0.08} $$ for z = 7 − 8 and α = −1.68^+0.13_−0.14 for z = 8 − 9). The combination of a decreasing metallicity and bursty SFH can reconcile the observed differences between the UV and [O III]+Hβ LF. By converting this LF into the ionising photon-production rate Ṅ_ion, we show that galaxies with L_Hα ≥ 10³⁹ erg s⁻¹, that is, with a star formation rate (SFR) (Hα) ≥ 5 × 10⁻³ M_⊙ yr⁻¹) cause 31%−90% and 46%−156% of the ionising photon budget (at 7 < z < 8 and 8 < z < 9), when we assume a constant escape fraction of Lyman-continuum photon (f_esc = 0.14). The shape of the LF further shows the negligible contribution of faint galaxies to the Ṅ_ion. Additionally, we derived the cosmic star formation rate density (SFRD), finding results consistent with previous estimates. However, the sensitivity of GLIMPSE to lower SFRs reinforces the conclusion that very faint galaxies contribute very little to Ṅ_ion and the SFRD. Our results suggests that GLIMPSE has detected the bulk of the total [O III]+Hβ emission from star-forming galaxies, and that galaxies below our detection limits are likely minor contributors to cosmic re-ionisation.

Key words: galaxies: dwarf / galaxies: high-redshift / galaxies: luminosity function / mass function / galaxies: photometry / galaxies: star formation / dark ages / reionization / first stars

© The Authors 2026

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

During its first billion years, the Universe underwent significant transformations. After a gradual cooling, the Universe entered the so-called dark ages, a period in which photons decoupled from matter, rendering the Universe opaque to our telescopes (Ferrara & Pandolfi 2014; Adamo et al. 2025). This epoch ended around z ∼ 30 with the so-called epoch of re-ionisation (EoR) with the formation of massive and metal-free Population III stars (Ferrara 1998). These newly born objects emitted enough ionising photons to ionise the surrounding gas while dominating the recombination of hydrogen (Ferrara 1998; Bromm 2013). They emitted sufficient ionising photons to counteract recombination and ionise the surrounding gas in a patchy manner. Each bubble of ionised gas expanded and leaked more ionising photons, which ultimately re-ionised the whole Universe (e.g. Zaroubi 2013; D’Aloisio et al. 2015; Bosman et al. 2022; Robertson 2022; Korber et al. 2023; Meyer et al. 2025. The EoR is expected to have ended around z ∼ 5 − 6 according to measurements of the Gunn-Peterson effect in quasar spectra (e.g. Becker et al. 2001; Fan et al. 2006; Keating et al. 2020; Bosman et al. 2022) and the decreasing fraction of Lyman α emitting galaxies detected above z > 6 (Stark et al. 2011; Schenker et al. 2012). While galaxies are often considered as primary drivers of re-ionisation, the relative importance of fainter or brighter galaxies remains debated. A first scenario gives a prime role to faint galaxies because despite their small individual contribution to the ionising photon budget, their high number density might produce enough to re-ionise the Universe (e.g. Oesch et al. 2009; Finkelstein et al. 2019; Yeh et al. 2023; Atek et al. 2024; Simmonds et al. 2024b). Another scenario gives more weight to bright galaxies, which have a large individual contribution despite their rarity (Naidu et al. 2020). However, some studies also suggested that active galactic nuclei (AGNs) might also play a significant role in re-ionisation (Dayal et al. 2020; Madau et al. 2024; Maiolino et al. 2024; Grazian et al. 2024; Singha et al. 2025).

Until recently, probing galaxies during the EoR was hindered by technical limitations. The Hubble Space Telescope (HST) and Spitzer Space Telescope enabled us to probe the end of the EoR (z ≳ 6) in great detail (e.g. Ellis et al. 2013; Oesch et al. 2018; Atek et al. 2018; Bouwens et al. 2019) owing to very deep surveys such as the Hubble frontier field (HFF; Lotz et al. 2017). However, studies at higher redshift remained limited to the observations of a handful of very bright galaxies (z > 8 objects; (e.g. Zitrin et al. 2015b; Oesch et al. 2016) due to the wavelength coverage and limited sensitivity of HST, the resolution of Spitzer, and the transmission of IR light for ground-based astronomy. Since 2022, the James Webb Space Telescope (JWST) enables us to probe farther and deeper into the EoR with regular detections of record redshift galaxies (e.g. Curtis-Lake et al. 2023; Wang et al. 2023; Fujimoto et al. 2023; Harikane et al. 2024; Napolitano et al. 2025; Naidu et al. 2026) and studies of the faintest galaxies (e.g. Eisenstein et al. 2023; Bezanson et al. 2024; Suess et al. 2024; Finkelstein et al. 2025).

Uncovering the behaviour of the faintest galaxies (i.e. M_UV ≥ −17) is paramount to understanding the EoR. Faint galaxies are very numerous in the Universe, but remain difficult to study. The star formation in these galaxies is highly unstable and characterised by periods of intense activity, followed by prolonged quiescence (e.g. Endsley et al. 2025). On the one hand, the large abundance of extreme emission line emitters (EELGs), that is, of galaxies with emission lines reaching a rest-frame equivalent width of ≳1000 Å, shows the signature of low-mass galaxies undergoing an intense star formation period (Rinaldi et al. 2023). While rather rare at low redshift (Matthee et al. 2023), they appear to be more common at high redshift (e.g. Atek et al. 2011; van der Wel et al. 2011; Smit et al. 2014; Tang et al. 2019; Izotov et al. 2020; Onodera et al. 2020; Berg et al. 2021; Davis et al. 2024; Llerena et al. 2024; Rinaldi et al. 2025; Boyett et al. 2024; Endsley et al. 2024; Begley et al. 2025; Daikuhara et al. 2025). On the other hand, the numerous detections of galaxies with low star formation rates (SFRs), variously referred to as (mini-)quenched, in a phase of SFR downturn, or dormant (e.g. Gelli et al. 2023; Strait et al. 2023; Dome et al. 2024; Endsley et al. 2025; Looser et al. 2024, 2025; Mintz et al. 2025; Trussler et al. 2025; Covelo-Paz et al. 2026), shows the complex and bursty nature of star formation in the EoR, with short intense star formation periods followed by a longer period of very low star formation. EELGs are more readily detectable as their strong emission lines boost the filters through which they are observed (Schaerer & de Barros 2009). Faint low-mass galaxies, with low star formation, lack this additional boost, which limits their detection. Their importance is further confirmed by recent studies, however, which showed evidence for missing low-mass galaxies in a statistical sample, with a low SFR, which skews the statistics (e.g. Endsley et al. 2025; Simmonds et al. 2024b).

The luminosity function (LF) characterises the number density of galaxies for a certain luminosity, providing a powerful description of the galaxy population demographics. Measuring the LF in different fields enables us to compare galaxy populations and examine variations in different regions. The UV LF was extensively studied in the EoR to understand the contribution of galaxies to the re-ionisation (e.g. Richard et al. 2008; Oesch et al. 2009, 2018; Atek et al. 2018; Bouwens et al. 2020; Moutard et al. 2020; Bowler et al. 2020; Bouwens et al. 2022; Willott et al. 2024; Harikane et al. 2025; Weibel et al. 2025; Chemerynska et al. 2026, Atek et al. in prep.). Knowing the exact number density distribution enables us to assess the contribution of each type of galaxy by measuring their ionising photon-production rate (e.g. Naidu et al. 2020; Atek et al. 2024; Simmonds et al. 2024b). However, the challenges in constraining the Lyman-continuum escape fraction (f_esc) at high redshift means that this question remains unresolved. Jecmen et al. (2026) recently studied f_esc for the JWST GLIMPSE survey using its relation to β-slope (e.g. Chisholm et al. 2022). Their results indicated a constant f_esc ∼ 14%, consistent with the f_esc ∼ 10% from earlier studies (Robertson et al. 2013, 2015; Giovinazzo et al. 2026; Mascia et al. 2025), but extends these constraints to significantly fainter galaxies. Emerging evidence, including from Jecmen et al. (2026), further suggests that f_esc might evolve with galaxy mass, thereby significantly affecting the relative contribution of faint and bright galaxies to re-ionisation (e.g. Begley et al. 2022; Saldana-Lopez et al. 2023; Jecmen et al. 2026. In addition, the ionising photon-production efficiency ξ_ion also requires attention. By overestimating it, the measured steep UV LF results in an overshooting of the ionising photon-production rate, and therefore, in fast re-ionisation (Robertson et al. 2013, 2015; Muñoz et al. 2024; Simmonds et al. 2024b; Bosman & Davies 2024).

The ionisation of the interstellar medium of galaxies produces a wealth of information through emission lines (see Kewley et al. 2019). Hα is a major tracer of the instantaneous star formation in galaxies (Kennicutt & Evans 2012), but its observation with JWST/NIRCam is limited to redshift z ≤ 6.5. Therefore, we used the second-strongest Balmer series line Hβ to study star formation. However, due to the proximity of this line with [O III]λλ4960, 5008 (hereafter [O III]), we were unable to distinguish them with broadband photometry, motivating the study of the combined [O III]+Hβ complex. However, this adds complexity to tracing the star formation because the [O III] doublet and R3 = [O III]λ5008-to-Hβ ratio are heavily affected by the ISM metallicity (e.g. Curti et al. 2017; Maiolino & Mannucci 2019; Curti et al. 2020; Sanders et al. 2024; Scholte et al. 2025).

The [O III] or [O III]+Hβ LF has been measured multiple times for different fields (e.g. Colbert et al. 2013; Khostovan et al. 2015; De Barros et al. 2019; Khostovan et al. 2020; Bowman et al. 2021; Matthee et al. 2023; Sun et al. 2023; Nagaraj et al. 2023; Meyer et al. 2024; Wold et al. 2025), but was limited at high redshift by the wavelength coverage of WFC3/HST (z ≲ 3), the resolution of Spitzer/IRAC, and the poor infrared transmission through atmosphere for ground-based astronomy. The early results from (De Barros et al. 2019) were able to push the redshift limits to z ∼ 8 using a combination of Spitzer 3.6 μm and 4.5 μm, the deep fields of HST, and an assumed relation between the UV LF and the [O III] LF. Meyer et al. (2024) measured a non-biased [O III] LF using the JWST/NIRCam/WFSS FRESCO (Oesch et al. 2023) survey for relatively bright (M_UV ≲ −18) galaxies between 6.8 < z < 9.0. They found a number density lower than previously estimated and observed a rapid decline in the [O III] number density at z ≳ 7. Matthee et al. (2023) measured the [O III] LF at 5 < z < 7 using JWST/NIRCam/WFSS EIGER (Kashino et al. 2023) images and found little evolution of the [O III] LF when compared to lower-redshift ranges (z ≲ 5; Khostovan et al. 2015). However, until recently, the faint end of the LF was either indirectly constrained or simply fixed because only a few sources were detected. This is significant because the faint-end slope critically affects our understanding of the evolution of re-ionisation and the role of low-mass galaxies during this era (Naidu et al. 2020; Atek et al. 2024; Muñoz et al. 2024). Wold et al. (2025) was the first study to push the faint-end limits of the [O III] LF back at high redshift. To achieve this, they used photometric data from the strongly lensed Abell 2744 HFF field, observed by the JWST ultradeep NIRSpec and NIRCam observations before the epoch of re-ionisation team (UNCOVER, Bezanson et al. 2024), which made use of magnification to obtain deeper observations. So far, Wold et al. (2025) alone pushed the faint-end limits of the [O III] LF back at high redshift (L_{[O III]λ5008} ≳ 10⁴¹ erg s⁻¹ for Wold et al. (2025), whereas L_{[O III]λ5008} ≳ 10^{41.75 to 42} erg s⁻¹ for Matthee et al. (2023), Meyer et al. (2024).

To be able to reach the faintest galaxies and constrain the faint end of the [O III]+Hβ LF, we used the deepest JWST survey to date, GLIMPSE, on the strongly lensed HFF galaxy cluster Abell S1063. The unprecedented depth of the survey coupled to the strong magnification of some galaxies enabled us to reach the faintest galaxies ever observed at high redshift. We used a Lyman-break selection (Sect. 3.1) and constrained the strong-lensing (SL) model to deduce magnification, multiplets, and effective volumes (Sects. 3.2 and 3.3). We then made use of SED fits to measure their emission line fluxes (Sect. 3.4), we estimated the completeness of our sample (Sect. 3.6), and we constructed the [O III]+Hβ LFs and analysed their behaviour (Sect. 4). Finally, we discuss the implication of the resulting LFs with respect to the ionising photon budget required to re-ionise the Universe and the measured cosmic star formation rate density (SFRD; Sect. 5). We considered a flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7. All magnitudes are expressed in the AB system (Oke & Gunn 1983), and we adopted the (Chabrier 2003) IMF.

2. Data

2.1. Observations

The GLIMPSE survey is a cycle 2 JWST/NIRCam large program (GO-3293; PI Atek & Chisholm; Atek et al. 2025) that performed the deepest observations of the HFF (Lotz et al. 2017) galaxy cluster Abell S1063 (z = 0.348). The observations consist of one pointing with two modules: the first one centred on the lensed field and the second one on a nearby region. The field was observed with 7 JWST/NIRCam broadband filters (F090W, F115W, F150W, F200W, F277W, F356W and F444W) reaching 5σ depths of ∼30.9 mag and 2 JWST/NIRCam medium bands (F410M, F480M) reaching 5σ depth of ∼30.1 mag (Atek et al. 2025). In addition to GLIMPSE data, we also used shallow cycle 1 JWST/NIRCam observations with the F250M and F300M medium bands (Hashimoto et al. 2023), as these bands could provide some constraint on the SED fitting. In addition to the JWST/NIRCam observations, we also included the legacy data from the HST from the HFF survey (Lotz et al. 2017) and the beyond ultra-deep frontier fields and legacy observations survey (BUFFALO, Steinhardt et al. 2020). These observation provide deep imaging from HST/ACS (F606W, F814W) and HST/WFC3 (F105W, F125W, F140W, F160W), which constrain the rest-frame UV and blue-optical of the targeted galaxies in this study.

All the data were reduced following the procedure from (Endsley et al. 2024). In summary, the Point Spread Function (PSF) was built from the stars observed in the field. The JWST and HST images were then PSF matched to the reddest filter (F480M). The foreground bright cluster galaxies were subtracted following the method described in (Shipley et al. 2018; Weaver et al. 2024; Suess et al. 2024). Then, we performed source extraction using Sextractor (Bertin & Arnouts 1996) on an inverse-variance weighted stack of F277W+F356W+F444W, the broadband long wavelength filters of NIRCam. We extracted the photometry with an aperture size of D = 0.2″ and calculated the uncertainty by measuring the standard deviation from 2000 random photometric measurement at the same aperture in nearby empty regions. The small aperture size is motived by two key factors. First, the extreme depth of the GLIMPSE survey on the lensed field Abell S1063 results in a highly clustered field. Secondly, our primary science goal is the study of faint high-redshift galaxies, which are compact. Therefore, a smaller aperture ensures that we maximise the signal-to-noise for these galaxies, by reducing nearby contamination. We corrected the aperture using empirical measurements of the PSF curve of growth, enabling us to measure the fraction of flux missing. The total flux is corrected from aperture by a factor 0.412. The full detail of the reduction is described in the GLIMPSE survey paper (Atek et al. 2025).

2.2. Simulations

In addition to observational data, we also used simulated fields to assess our results with simulations. We used the JWST extragalactic mock catalog (JAGUAR, Williams et al. 2018), which is a catalogue of pre-JWST simulated galaxies with the aim of preparing for JWST science. These galaxies span a wide redshift range 0.2 < z < 15 for galaxies of masses M_* ≥ 10⁶ M_⊙. This simulation was mainly used to prepare for observations from the JWST advanced deep extragalactic survey (JADES, Rieke et al. 2023; Eisenstein et al. 2023). This simulation includes two branches: one containing star-forming galaxies and the other containing quiescent galaxies. In this paper, we only focus on the first realisation (R1) of the star-forming galaxies catalogue. To produce these catalogues, JAGUAR assumes an UV LF extrapolation from HST observations for its galaxy count, which can then be used to extrapolate the stellar mass of those galaxies. Galaxy models of BEAGLE (Chevallard & Charlot 2016) were then used to generate realistic SEDs of the galaxies. These models span many parameters, which are then validated against observations of that time, such as 3D-HST (Skelton et al. 2014). As observations back then were limited in specific regions of the parameter space, such as low stellar mass galaxies (M_* ≤ 10⁸ M_⊙), they relied on the theoretical knowledge of the time to forward-model galaxies. In practice, they enforce a sharp cut at M_* ≥ 10⁶ M_⊙, which will affect our very faint galaxies comparison (Williams et al. 2018). However, galaxies with M_UV ≤ −16 should not be affected. While these simulations do not reproduce exactly current observations, they still enable us to test and calibrate our methods (e.g. Meyer et al. 2024; Sun et al. 2023).

3. Method and results

3.1. Source selection

The selection of the present sample is based on a combination of the Lyman-break technique and photometric redshift estimates. The same approach is used by Chemerynska et al. (2026) at redshifts z > 9. We first apply a colour-colour selection that aims at identifying the flux dropout in the bands blueward of rest-frame Lyman-α, and at the same time minimise the contamination of low-redshift red sources by using the colours from the redder bands. For galaxies in the redshift range 7 < z < 9, we relied on a dropout in the F090W filter by adopting the following colour-colour selection criteria:

$\begin{matrix} \begin{matrix} m_{090} - m_{115} > 0.8 \\ \land m_{115} - m_{150} > 1.6 + (m_{150} - m_{200}) \\ \land m_{115} - m_{150} < 0.5 \lor m_{150} - m_{200} < 0.5 \\ \land {S / N}_{115} > 5 \lor {S / N}_{200} > 5 \end{matrix}, \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \begin{array}{l} m_{090}-m_{115}>0.8\\ \wedge \ m_{115}-m_{150}>1.6+(m_{150}-m_{200})\\ \wedge \ m_{115}-m_{150}<0.5 \vee m_{150}-m_{200} < 0.5 \\ \wedge \ \mathrm{S/N}_{115}>5 \vee \mathrm{S/N}_{200}>5 \end{array}, \end{aligned} $$$ (1)

where we also required galaxies to be detected with high significance, at a 5σ level, in at least one band. In addition, we verify that all galaxies remain undetected at a 2σ level in all of the HST bands. Furthermore, we restricted our search area to regions with at least three-dither overlap to mitigate contamination by spurious sources. Accordingly, the survey area presented in Sect. 3.3 has been corrected for these restrictions of the search area and we end up with an effective survey area of 4.30 to 4.34 arcmin².

The second step of this selection is based on photometric redshifts computed with Eazy (Brammer et al. 2008). We used the following parameters during the SED fitting procedure: the blue_sfhz_13 templates, dust attenuation following (Calzetti et al. 2000) and IGM attenuation following (Inoue et al. 2014). The Eazy redshift grid spanned the full range, from z = 0.01 to z = 30. We ensure that all dropout-selected galaxies have a best-fit solution at redshifts higher than z > 6. The final sample contains 173 galaxies between 7 < z < 9, and their M_UV distribution is shown in Fig. 1 ( $M_{UV} = - 16 . 57_{- 0.97}^{+ 1.63}$ $Mathematical equation: $ M_{\mathrm{UV}}=-16.57_{-0.97}^{+1.63} $$ ).

Fig. 1.

Relation between the de-lensed absolute UV magnitude M_UV of the GLIMPSE F090W dropouts selected in this work and their SL magnification μ. Beyond M_UV = −16, we need magnification to detect these objects. This is shown by the grey area, which is the minimum magnification required for a 5σ detection at z = 7 to 9 with F115W, our deepest band (see Atek et al. 2025). The top and right histograms show the normalised distribution of M_UV and μ for each redshift bins.

3.2. Gravitational lensing

In order to account for the strong gravitational lensing (SL) effect of Abell S1063, we use a new parametric SL mass model of the cluster, constructed with the Zitrin et al. (2015a) parametric method (see Furtak et al. 2023), revised to be fully analytic, i.e. not limited to a grid resolution. The lens model will be presented in detail in Furtak et al. (in prep.) and we refer the reader to that work for more details.

In short, the model comprises two smooth dark matter (DM) halos parametrised as pseudo-isothermal elliptical mass distribution (PIEMD; Kassiola & Kovner 1993), one centred on the brightest cluster galaxy (BCG), and the other on a group of galaxies in the north-east of the cluster as found by previous works (Bergamini et al. 2019; Beauchesne et al. 2024). In addition, the model includes 303 cluster member galaxies, parametrised as dual pseudo-isothermal ellipsoids (dPIEs; Elíasdóttir et al. 2007). The model is constrained with 75 multiple images of 28 sources, 24 of which have spectroscopic redshifts. Optimised using a Monte Carlo Markov chain (MCMC) analysis, the lens model achieves a final lens-plane image reproduction RMS of Δ_RMS = 0.54″. This model has previously been used various projects (e.g. Topping et al. 2024; Kokorev et al. 2025; Fujimoto et al. 2025; Chemerynska et al. 2026).

We use the lens model to compute the gravitational magnification at the coordinates and redshift of each objects. All sources, even in the off-cluster module, are magnified by at least a factor μ = 1.19. In Fig. 1, we show the explored M_UV against magnification μ. We see that for M_UV ≥ −16, stronger magnifications are required to even detect the objects, reaching max(μ)∼219 for our strongest-lensed galaxy.

We also need to account for multiple imaging in our object counts when deriving the LF in order to avoid double- or even triple-counting galaxies. This is done in an iterative process where image positions are injected into the lens model to predict eventual counter image positions. We then search the data for corresponding sources in the 3 arcseconds area surrounding the predicted counter-images location and carefully matched photometries and SEDs to make sure they are the same object. In the event of multiple imaging, we keep the counter-image with the highest F444W signal-to-noise. From our sample of 173 galaxies, 68 expect one or more counter-images. 9 of them can be found within our Lyman-Break Galaxy LBG sample at 7 < z < 9, which reduces our total sample size to 164 galaxies.

3.3. Estimating the effective survey volume

Estimating the effective volume of each source relies on our SL model of the field. As described in Atek et al. (2018), the effective volume is estimated as the integral of the co-moving volume that comes from the cosmology to which we apply completeness and magnification. It corresponds to the source plane volume in which we effectively observe a galaxy of observed magnitude m. Typically, a faint source can only be detected in highly magnified regions, which reduces the effective volume for these sources. The effective volume can be analytically computed and is defined as

$\begin{matrix} V_{eff} = \int_{0}^{\infty} d z \int_{μ > μ_{\min}} d Ω (μ, z) \frac{V_{com}}{d z} C (z, m, μ) \approx V \times C, \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \text{ V}_{\rm eff} = \int _0^\infty \mathrm{d}z\int _{\mu >\mu _{\rm min}}\mathrm{d}\Omega (\mu , z) \frac{\text{ V}_{\rm com}}{\mathrm{d}z} C(z, m, \mu ) \approx V \times C , \end{aligned} $$$ (2)

where V_eff is the effective volume, V_com is the co-moving cosmological volume, C(z, m, μ) the detection completeness for sources at redshift z, with apparent magnitude m and magnification μ, and dΩ(μ, z) is the surface element that depends on magnification and redshift (Atek et al. 2018). The volume uncertainties are obtained by propagating the magnification-dependent systematic error in the cumulated source plane area, the error on the completeness described in Sect. 3.6, and the variation of the effective volume in the M_UV bin. In this work, we extracted the completeness out of the integral by considering discrete bins of completeness and volumes instead of continuous functions. We refer to V as the average volume and to C as the average completeness over a given M_UV bin.

Figure 2 shows the evolution of the effective volume with respect to the galaxy UV magnitude M_UV. For brighter sources (M_UV < −17), a plateau is reached that corresponds to sources that should be detected in the entire image regardless of magnification, but for fainter sources (M_UV > −17), the magnification becomes important for the detection, which then quickly reduces the surveyed effective volume.

Fig. 2.

Top panel: Evolution of the effective volume probed for each UV magnitude M_UV for our two redshift bins 7 < z < 8 (blue) and 8 < z < 9 (orange) with 1σ uncertainties in the shaded area. Bottom panel: Completeness estimate for the same M_UV and redshift bins.

In summary, we integrated the source planes from the SL model for sources at different redshifts, with different apparent magnitudes and located in different magnification regions to obtain the effective volumes for each galaxy.

3.4. Line flux measurement of [O III]+Hβ

We measured line fluxes by fitting the SEDs of our dropout-selected galaxies. For this, we use the CIGALE (Boquien et al. 2019) SED fitting software with the fixed redshift obtained by EASY. We used a delayed-τ (sfhdelayed) star formation history with an e-folding time and an age interval between 1 Myr and 13.5 Gyr. Galaxies are bound to the redshift, and therefore, no galaxies exceed 500 Myr. We used the Bruzual&Charlot (bc03) (Bruzual & Charlot 2003) simple stellar populations assuming a Chabrier initial mass function (IMF; Chabrier 2003) and a fixed gas and stellar metallicity of 0.004 (the mass fraction of atoms heavier than helium, where Z_⊙ = 0.02). We include nebular templates with the ionisation parameter varying between log U = −1 and −4 and an electron density of 100 cm⁻³. We also account for dust attenuation, using a module based on a modification by Leitherer et al. (2002) of the law from Calzetti et al. (2000). We vary the E(B − V)_l colour excess of the nebular lines between 0 and 1.5 and a power law modifying the attenuation curve between –0.5 and 0. In addition, we also adopt the Milky Way extinction curve from Cardelli et al. (1989) updated by O’Donnell (1994) to better match high resolution extinction curves from Bastiaansen (1992). More details on the different fitting modules are described in Boquien et al. (2019). We measure the line flux using the internal logic from CIGALE, which uses a Bayesian estimation for each measured parameter by weighting templates by their likelihood. In CIGALE, the gas and stellar metallicity are considered as two separate quantities. Because of their degeneracies with other parameters, they end up being poorly constrained and might yield unphysical solutions with large differences between stellar and gas metallicities. Therefore, we decided to keep the metallicity fixed to Z = 0.004 in the fitting. As a test, we also varied the metallicity, but no significant differences in the LFs were observed. Therefore, for the performance improvement as well as coherence of the gas and stellar metallicity, we decided to keep the SED fitting metallicities parameters constant. Three representative galaxies with stamps cutouts and SED can be found in Appendix B.

To validate our measurement method, we assessed it with the JAGUAR simulations (Williams et al. 2018), which showed a good correlation (correlation coefficients of 0.93 and 0.79) between the absolute JAGUAR value of [O III]+Hβ flux and EW and the retrieved SED fitting estimation of the same parameters. In addition, we also validated our method through a direct measurement of the [O III]+Hβ flux and EW based on excess in the available F356W, F410M, F444W and F480M bands, which again showed a reasonably good correlation coefficient of 0.67 for the flux, but 0.20 for the EW. We see some divergence on the faint end for the EW, but this most probably comes from greater uncertainties in the direct measurement due to a lack of medium band availability. The details of the validation can be found in Appendix A.

3.5. Dust attenuation

Unless otherwise stated, the results reported in this paper correspond to observed quantities, which are not corrected for dust attenuation. To correct for attenuation, we use SED fitting, as described in Sect. 3.4, to infer the line colour excess E(B − V)_l which is defined in Calzetti et al. (2000) as the continuum E(B − V) divided by an empirical factor 0.44 ± 0.03. The inferred emission line fluxes (or luminosities) are then corrected for dust attenuation, adopting the extinction law from (Cardelli et al. 1989). Figure 3 shows the E(B − V)_l distribution of our sample. The median measurement lies at E(B − V) = 0.013, and the majority of the galaxies are compatible with a low attenuation of E(B − V)_l < 0.1. We do not find significant correlations between dust attenuation and parameters such as M_UV. Therefore, these low E(B − V)_l values imply that dust is unlikely to play an outsized role in the inferred [O III]+Hβ flux values.

Fig. 3.

Histogram of the dust attenuation parameter E(B − V)_l sampled from the SED fitting of our sample of 7 < z < 9 galaxies. N = 10 000 realisations of the measurement are performed for each galaxy, following the distributions of the SED fitting, and the histogram is normalised by N.

3.6. [O III]+Hβ completeness estimation

Given the non-trivial nature of oxygen lines, we used an indirect method to estimate the completeness of our sample via M_UV. We used the completeness estimation made by the GLIMPSE collaboration (e.g. Chemerynska et al. 2026), which follows the procedures described in Chemerynska et al. (2024b), Atek et al. (2018). In summary, synthetic galaxies were added into the source plane, mapped into the lens plane using the SL model (see Sect. 3.2), added to the images, and then selected using the same process as described earlier. This M_UV-based completeness does not exactly correspond to our [O III] completeness. However, M_UV and the [O III]+Hβ line flux correlate (Matthee et al. 2023). Figure 4 shows the correlation between the [O III]+Hβ line flux and the M_UV for the GLIMPSE sample and the sample from Meyer et al. (2024). The sample shows a good correlation between the two quantities, but shows a bimodal distribution of L_{[O III]+Hβ} for M_UV > −17. This bimodality emerges from the poorer constrains on the faintest observed objects. At this magnitude, the galaxies are near the observed detection limit of GLIMPSE and their signal-to-noise is limited. Therefore, degeneracies between fitting parameters increase and the flux estimation becomes uncertain as many model fit the data. These galaxies are then propagated to lower magnitude M_UV due to magnification, which gives this second mode in the [O III]+Hβ-to-UV relation. However, while their flux measurements are more uncertain, the values inferred from the model are physically motivated and the posterior distributions will be sampled to account for their variety. Therefore, despite the scatter, the quantities correlate well, and the completeness for M_UV approximately translates to [O III]+Hβ line flux.

Fig. 4.

Correlation between the [O III]+Hβ line flux against M_UV for our GLIMPSE sample (circles) and measurements from Meyer et al. (2024, shown as crosses). To match our observations, we only kept sources within the same 7 < z < 9 redshift range. Bayesian error estimates using CIGALE are also shown for the GLIMPSE measurements.

To measure uncertainties for the completeness, we assumed that it follows a binomial distribution B(n, p) where n is the total number of synthetic sources and p is the completeness for this type of source (depending on the M_UV and magnification μ of the source). This yields the estimation of the standard deviation for the completeness σ_C, with C the completeness, as given by

$\begin{matrix} δ C = \sqrt{\frac{C (1 - C)}{n}} . \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \delta C = \sqrt{\frac{C(1-C)}{n}} . \end{aligned} $$$ (3)

4. The [O III]+Hβ luminosity function at z ∼ 7 − 9

We now describe our standard approach to measuring the [O III]+Hβ LF. We computed the LF using the following equation:

$\begin{matrix} Φ (L) d {log}_{10} (L) = \sum_{i} \frac{1}{V_{i} C_{i}}, \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \Phi (L)\mathrm{d}\log _{10}(L) = \sum _i \frac{1}{V_i C_i} , \end{aligned} $$$ (4)

where L is the [O III]+Hβ line luminosity, V_i is the effective survey volume for a given source i, C_i is its associated completeness at given redshift and L.

To construct the LF, we first varied the line flux [O III]+Hβ (and M_UV for the completeness correction) from the measurements of the fitted SED templates. For this, we extracted the measurements of both quantities from CIGALE, as well as the associated template χ². By weighting each measurement by the likelihood ℒ = exp(−0.5χ²), we obtain a 2D histogram (N-dimensions when dust attenuation or other parameters are sampled) from which we sampled n = 10 000 realisations of the parameters. To avoid discretisation effects caused by the sampling of the histogram, we uniformly varied the measurements within their bins. The [O III]+Hβ line flux was treated on a logarithmic scale. As the SED was fitted to the observed photometry, the magnification was also sampled following a normal distribution, and was applied to M_UV and [O III]+Hβ flux.

From there, we measure the LF and its associated uncertainties from each realisation. While the LF is measured following Eq. (4), the uncertainties are more complex and come from multiple sources. First, we measured the intrinsic error of each luminosity bin coming from the uncertainties of the effective survey volume (which includes the completeness). By propagating uncertainties, we end up with the following:

$\begin{matrix} δ_{I} {log}_{10} Φ (L) = \frac{\sqrt{\sum_{i} {(\frac{δ V_{i}}{V_{i}^{2} C_{i}})}^{2} + {(\frac{δ C_{i}}{V_{i} C_{i}^{2}})}^{2}}}{ln (10) d {log}_{10} (L) \sum_{i} \frac{1}{V_{i} C_{i}}}, \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \delta _{\rm I}\log _{10}\Phi (L) = \frac{\sqrt{\sum _i \left(\frac{\delta V_i}{V_i^2C_i}\right)^2 + \left(\frac{\delta C_i}{V_iC_i^2}\right)^2 }}{\ln (10) d\log _{10}(L) \sum _i \frac{1}{V_i C_i}} , \end{aligned} $$$ (5)

where δ_Ilog₁₀Φ(L) is the intrinsic uncertainty on the log-LF, δV_i is the uncertainty on the volume and δC_i is the uncertainty on the completeness.

In addition to the intrinsic errors, we also account for Poissonian errors given by the number of sources in each bin of the LF. For a bin uncertainty given by $\sqrt{N}$ $Mathematical equation: $ \sqrt{N} $$ , with N the number of sources in each luminosity bin, we obtain the following:

$\begin{matrix} δ_{P} {log}_{10} Φ (L) = \frac{1}{\sqrt{N} ln (10) d {log}_{10} (L)} . \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \delta _{\rm P}\log _{10}\Phi (L) = \frac{1}{\sqrt{N}\ln (10)\mathrm{d}\log _{10}(L)} .\end{aligned} $$$ (6)

The intrinsic error and the Poissonian error are then summed in quadrature to obtain the final measurement of the LF for each realisation. We therefore have n = 10 000 realisations of the same LF, where each luminosity bin has a measurement with uncertainties. To combine them all into one final LF, we sample 100 subsample of each LF realisation, giving us 10 000 × 100 LF for which we measure mean and standard deviation.

In order to reduce biases from very low statistics due to non-detection in the faint end and small volume in the bright end, we applied two cuts to the LFs. First, each luminosity bin needs to have at least one source on average (⟨N⟩ ≥ 1), meaning that out of all realisations, on average, at least one object should lie in the luminosity bin. To avoid cases with large uncertainties on the [O III]+Hβ line flux, we measure the median value for each galaxy and we further enforced that each [O III]+Hβ bin needs at least one median object (N ≥ 1).

While we use a Bayesian approach for our statistics, we would like to estimate the quality of our results using some signal-to-noise ratio (S/N) estimation. For this, we assume a log-normal distribution of our luminosity bins and use the general definition of the expected value (E) and the variance (Var) to measure our signal-to-noise. Then, using μ and σ as the bins mean and standard deviation, we estimate the following equations:

$\begin{matrix} E & = exp (μ + \frac{σ^{2}}{2}) \end{matrix}$ $Mathematical equation: $$ \begin{aligned} E&= \exp (\mu + \frac{\sigma ^2}{2})\end{aligned} $$$ (7)

$\begin{matrix} Var & = (e^{σ^{2}} - 1) exp (2 μ + σ^{2}) \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \mathrm{Var}&= (e^{\sigma ^2} - 1) \exp (2\mu + \sigma ^2)\end{aligned} $$$ (8)

$\begin{matrix} S / N & = \frac{E}{\sqrt{Var}} . \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \mathrm{S/N} = \frac{E}{{\sqrt {\mathrm {Var}}}}. \end{aligned} $$$ (9)

The S/N of each bin is reported in their associated tables.

4.1. The [O III]+Hβ luminosity function at z ∼ 7 − 9

In Fig. 5 we show the [O III]+Hβ LF for our two redshift ranges, and their comparison to the values from literature. We separated our LF in luminosity bins of 0.5 dex, resulting in 9 bins for the z = 7 − 8 (⟨z⟩ ∼ 7.48) redshift range and 6 bins for the z = 8 − 9 (⟨z⟩ ∼ 8.33) redshift range. We observe large uncertainties on the bright and faint-ends of the LF. On the bright end, the low number statistics drives the Poissonian errors, but on the faint-end, the addition of the uncertainties caused by degeneracies between the different solutions in the SED fitting models and the small effective volume increases the uncertainties. Nevertheless, we have a relatively low uncertainty between 10⁴⁰erg s⁻¹ and 10⁴² erg s⁻¹, with S/N > 3. The strength of GLIMPSE is the depth reached. Indeed, by combining a photometry approach on a very deep lensed galaxy field, we are able to reach unprecedented [O III]+Hβ line fluxes. The deepest comparison from Wold et al. (2025) allowed us to securely probe the [O III]+Hβ LF down to ∼10^41.5 erg s⁻¹, while GLIMPSE goes down to ∼10³⁹ erg s⁻¹, which constrains the high redshift faint end of the [O III]+Hβ LF to unprecedented depths. The full detail of the [O III]+Hβ LF is provided in Table C.1.

Fig. 5.

Luminosity function of Hβ+[O III]λλ4960, 5008 for our two redshift ranges. The GLIMPSE data points give the mean and standard deviation of the LF, while the fits provide the median value as a solid line, and the 16–84 percentile in the shaded area. We added previous studies of the LF using JWST/NIRCam GRISM slitless spectroscopy (Meyer et al. 2024; Matthee et al. 2023; Sun et al. 2023), a JWST/NIRCam medium band survey by Wold et al. (2025) and a former Spitzer study by De Barros et al. (2019). All the JWST surveys specifically study the [O III]λ5007Å, so we converted them to [O III]+Hβ by considering their respective measured R3 factor as well as [O III]λ5008/[O III]λ4960 = 2.98 (Storey & Zeippen 2000). The data can be found in Table C.2 and the parametrisation can be found in Table 1.

Because most of the previous studies only computed the [O III]λ5008 emission line, we needed to convert these previous measurements into Hβ+[O III]λλ4960, 5008. To do this, we assumed an oxygen line ratio of [O III]λ5008/[O III]λ4960 = 2.98 (Storey & Zeippen 2000) and the R3 = [O III]λ5008/Hβ ratio measured or assumed in each study (R3 = 6.72 in Wold et al. (2025), 6.38 in Meyer et al. (2024), 6.30 in Matthee et al. (2023), 6.3 in Sun et al. (2023), and De Barros et al. (2019) report directly [O III]+Hβ). While this approach simplifies the conversion, since a rigorous treatment requires refitting each LF with individual R3 measurements, we adopt it to constrain the bright end of our LF, using the results from Meyer et al. (2024) as reference. The agreement between our LF and previous literature results in Fig. 5 suggests that this simplification does not significantly bias our results.

To parametrise our LFs, we use the common Schechter function (Schechter 1976) as given by

$\begin{matrix} ϕ (L) d L = ϕ^{*} {(\frac{L}{L^{*}})}^{α} exp (- \frac{L}{L^{*}}) d (\frac{L}{L^{*}}) . \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \phi (L)\mathrm{d}L = \phi ^*\left(\frac{L}{L^*}\right)^\alpha \exp \left(-\frac{L}{L^*}\right)\mathrm{d}\left(\frac{L}{L^*}\right) . \end{aligned} $$$ (10)

This function is parametrised by three variables: the normalisation ϕ^*, the characteristic luminosity L^* and the faint-end slope α. While we can probe very deep into the [O III]+Hβ LF, the field area (one NIRCam pointing in two modules) strongly limits the ability to probe brighter objects, which are less common than their faint counterparts. Indeed, we do not detect objects above 10⁴³ erg s⁻¹ and only a handful of objects above 10⁴² erg s⁻¹. Given that the characteristic luminosity is typically located in the vicinity of these luminosities, GLIMPSE does not enable us to strongly constrain the bright end of the [O III]+Hβ LF. Therefore, we combine our results with those of Meyer et al. (2024), who studied the bright end of the [O III]λ5008 LF using an unbiased sample of GRISM spectra. They measured the LF down to 10^41.75 erg s⁻¹ for the high-redshift (⟨z⟩ ∼ 7.9) sample, which is at a redshift comparable to our lower-redshift sample. Their [O III]λ5008 LF is converted to [O III]+Hβ thanks to the median R3 = [O III]λ5008/Hβ = 6.38 ± 0.85 measurement from their stacked spectroscopic spectra, enabling us to securely constrain the bright end of the LF, while leaving the faint-end constraints to GLIMPSE.

We used a MCMC approach to fit the LF. The likelihood ℒ that we use is defined as

$\begin{matrix} {log}_{10} L & = - \frac{1}{2} [\sum_{L}^{GLIMPSE} {(\frac{ϕ (L) - Φ (L, ϕ^{*}, L^{*}, α)}{σ_{ϕ} (L)})}^{2} \\ + & \sum_{L}^{FRESCO} {(\frac{ϕ (L) - Φ (L, ϕ^{*}, L^{*}, α)}{σ_{ϕ} (L)})}^{2}], \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \log _{10}\mathcal{L}&= -\frac{1}{2} \Bigg [\sum _{L}^\mathrm{GLIMPSE} \left(\frac{\phi (L) - \Phi (L,\phi ^*, L^*, \alpha )}{\sigma _\phi (L)}\right)^2 \nonumber \\ +&\sum _{L}^\mathrm{FRESCO} \left(\frac{\phi (L) - \Phi (L,\phi ^*, L^*, \alpha )}{\sigma _\phi (L)}\right)^2 \Bigg ], \end{aligned} $$$ (11)

where ϕ(L) is the mean observed density with uncertainty σ_ϕ(L) for a given luminosity L, Φ(L, ϕ^*, L^*, α) is the estimation of the parametrisation of the Schechter function at a given luminosity L. This parametrisation Φ corresponds to the convolution of the Schechter function with a Gaussian kernel, allowing us to account for the Eddington bias (Eddington 1913), which is a selection bias in which rare categories of objects are more often contaminated with more common categories of objects than the reverse because of their abundance. Finally, the first summation spans the luminosity range from GLIMPSE, and the second one from FRESCO (Meyer et al. 2024) at ⟨z⟩ ∼ 7.9. To estimate the uncertainties on the Schechter parameters, we sampled n = 10 000 realisations of the LF and fitted them with the above procedure. We then measured the median and 16–84% percentiles. The resulting fit parameters are listed in Table 1.

Table 1.

Schechter parametrisation of the LFs.

4.2. The characteristic luminosity and normalisation parameters

At the bright end of the LF, we obtain a characteristic luminosity of $L^{*} = 43 . 00_{- 0.14}^{+ 0.25}$ $Mathematical equation: $ L^* = 43.00_{-0.14}^{+0.25} $$ and $L^{*} = 42 . 85_{- 0.11}^{+ 0.14}$ $Mathematical equation: $ L^* = 42.85_{-0.11}^{+0.14} $$ for the ⟨z⟩ ∼ 7.48 and ⟨z⟩ ∼ 8.33 LFs respectfully. It agrees well with the literature results at similar redshifts (e.g. Meyer et al. 2024; Wold et al. 2025; De Barros et al. 2019) but differs from lower redshift studies, which indicates some evolution with redshift (e.g. Khostovan et al. 2015; Matthee et al. 2023). However, it is difficult to assess the impact of GLIMPSE in constraining the characteristic luminosity L^* because of the very limited sample. The uncertainties of GLIMPSE are large in the bright end, and the result is mostly driven by Meyer et al. (2024).

Regarding the normalisation, we obtain $ϕ^{*} = - 3 . 77_{- 0.33}^{+ 0.22}$ $Mathematical equation: $ \phi^* = -3.77_{-0.33}^{+0.22} $$ and $ϕ^{*} = - 3 . 51_{- 0.21}^{+ 0.17}$ $Mathematical equation: $ \phi^* = -3.51_{-0.21}^{+0.17} $$ for the ⟨z⟩ ∼ 7.48 and ⟨z⟩ ∼ 8.33 LFs respectfully. It agrees within uncertainties with studies at similar redshift (e.g. Meyer et al. 2024; Wold et al. 2025), except for De Barros et al. (2019) which quickly diverge from our observations. However, this last study did not measure [O III]+Hβ directly, but converted it from UV using a UV-to-[O III]+Hβ calibration with Spitzer and HST data, which can explain the visible difference. For lower redshifts, there seems to be an evolution between z = 0 − 3 with a decreasing normalisation with increasing redshift (Colbert et al. 2013; Khostovan et al. 2015, 2020; Bowman et al. 2021; Nagaraj et al. 2023), but little evolution between z = 3 − 9 (e.g. Khostovan et al. 2015; Matthee et al. 2023). GLIMPSE further confirms that trend with no statistically significant evolution between redshift 7 to 9. The evolution of all three fitting parameters is shown in Fig. 6.

Fig. 6.

Redshift evolution of the fitting parameters from Fig. 5.

4.3. Faint end of the [O III]+Hβ luminosity function

The unique intrinsic depth of GLIMPSE in combination with the strong lensing magnification enables us to derive for the first time the faint-end slope of the [O III]+Hβ LF (i.e. ≤10⁴¹ erg s⁻¹). We do find a significant differences between α at redshift 7 < z < 8 and 8 < z < 9 as we obtain $α = - 1 . 78_{- 0.06}^{+ 0.06}$ $Mathematical equation: $ \alpha = -1.78_{-0.06}^{+0.06} $$ and $α = - 1 . 55_{- 0.11}^{+ 0.11}$ $Mathematical equation: $ \alpha = -1.55_{-0.11}^{+0.11} $$ respectively, which is slightly above the 1σ threshold (see Table 1). This could be caused by the increasing difficulty of observing the faintest galaxies at higher-redshift, but this effect should be accounted for by completeness correction. However, an evolution of the metallicity could also be an explanation, with higher-redshift galaxies having less time in their history to synthesise metals.

We can compare this result to (Wold et al. 2025), who constrained the faint end of the [O III]λ5008 LF at z ∼ 7 from photometric data of the HFF A2744 strongly lensed field, assuming the R3 = 6.72 ratio measured by (Sun et al. 2023) to infer [O III]λ5008 from [O III]+Hβ. They obtain a slope $α = - 2 . 07_{- 0.23}^{+ 0.22}$ $Mathematical equation: $ \alpha=-2.07_{-0.23}^{+0.22} $$ , which is steeper than our result. We investigate possible reasons for this difference here below.

Firstly, the parametrisation is different as they use a double power law with a fixed value of log₁₀L^* = 42 erg s⁻¹, which could cause some different constraints on the fit. However, by refitting their observed LF with a Schechter function fixed at our log₁₀L^* ∼ 43.0 ± 0.3 erg s⁻¹ and ⟨z⟩ ∼ 7.9 data from Meyer et al. (2024) at the bright end, we obtain a similar $α = - 2 . 18_{- 0.13}^{+ 0.13}$ $Mathematical equation: $ \alpha = -2.18_{-0.13}^{+0.13} $$ . The faint-end slope is compatible between the two parametrisations, so the use of double power law with fixed characteristic luminosity does not explain the difference.

Another difference could be the ranges of faint-end luminosities probed by Wold et al. (2025) and GLIMPSE. The former is only constrained by a few data points, with limited detections, down to L ∼ 10⁴¹ erg s⁻¹, while GLIMPSE reaches to L ∼ 10³⁹ erg s⁻¹. By fitting our GLIMPSE LF again with a faint-end luminosity limit of L ≥ 10⁴¹ erg s⁻¹, which corresponds to the Wold et al. (2025) limit, we obtain $α = - 2 . 19_{- 0.12}^{+ 0.16}$ $Mathematical equation: $ \alpha = -2.19_{-0.12}^{+0.16} $$ for ⟨z⟩ ∼ 7.48 and $α = - 1 . 92_{- 0.25}^{+ 0.24}$ $Mathematical equation: $ \alpha = -1.92_{-0.25}^{+0.24} $$ for ⟨z⟩ ∼ 8.33, which is significantly steeper than our main result, and compatible with Wold et al. (2025). Further details on this can be found in Appendix C.

In addition to the bias induced by the differences in depth, the galaxies selected in Wold et al. (2025) do not include faint [O III]+Hβ galaxies. While we do not set limits on the [O III]+Hβ equivalent width, Wold et al. (2025) uses a colour-colour selection comparable to a cut of EW([O III]λ5008) > 500 Å (equivalent to EW([O III+Hβ) > 742 Å assuming their R3 = 6.72 and [O III]λ5008/[O III]λ4960 = 2.98). Figure 7 shows the distribution of [O III]+Hβ equivalent widths for the different redshift bins in GLIMPSE. A significant number of GLIMPSE galaxies (∼64%) have log(EW_{[O III]+Hβ}) ≤ 2.87, indicating that a significant number of galaxies reside below the detection threshold of Wold et al. (2025), which might skew their statistics (Endsley et al. 2025, 2024). Therefore, the differences between the selection method and survey depth explains the strong difference between the faint end of the GLIMPSE and Wold et al. (2025) [O III]+Hβ LF.

Fig. 7.

Distribution of the [O III]+Hβ equivalent width for each redshift bin. This figure is not completeness-corrected. The colours represent the two redshift ranges and the lines are KDE estimate of the distribution for clarity.

4.4. The faint-end slope of the nebular emission lines and UV luminosity functions

In Fig. 8 we report the evolution of the faint-end slope α for multiple LFs from past studies at various redshifts. We considered two types of studies: Nebular emission line LFs ([O III], Hα, Hβ) and UV LFs. First, we note a clear difference between their redshift evolution. While the faint-end slope of the nebular emission line LF, α_neb, stays approximately constant at α_neb ∼ −1.62 ± 0.21 (within large uncertainties) from redshift 0 to 9, the faint-end slope of the UV LFs, α_UV, steepens from α_UV ∼ −1.4 at z = 0 to α_UV ∼ −2.2 at redshift 9. This means that at high redshift, UV-faint galaxies outnumber UV-bright galaxies, while the proportion of nebular bright to faint galaxies stays approximately constant over time. This result is independent of the exact parametrisation of the LF – i.e. the Schechter function (Schechter 1976) or double power law (Dunlop & Peacock 1990) – as also shown in this figure.

Fig. 8.

Redshift evolution of the faint-end slope (α) in the Schechter (diamond) and double power-law (circle) fits for the UV continuum (grey symbols) and nebular emission line (coloured points) LFs. The orange squares with a black outline correspond to this study. The compilation of the data, including GLIMPSE, can be found in Table D.1.

Several processes or physical causes could explain the redshift evolution difference between the UV and the nebular faint-end slopes: (1) bursty or variable star formation histories, leading to an incoherence between M_UV and nebular lines due to their different timescales (2) a strong luminosity dependence of the metallicity, causing lower [O III]+Hβ emission at fainter luminosities, or (3) an intrinsic turnover of the UV LF at the faint-end, but below the range currently observed, possibly flattening the faint end of the [O III]+Hβ LF. We now discuss their effects on the LF one by one.

4.4.1. The impact of bursty star formation histories on the [O III]+Hβ LF

The different timescales of UV continuum or nebular line emission (≲100Myr and ≲10Myr, Kennicutt & Evans 2012) powered by star formation might cause a decoupling of the slope of the UV and [O III]+Hβ LF, as we would expect fainter galaxies to vary on shorter timescales (e.g. Endsley et al. 2024). Indeed, with variable or bursty star formation histories, these two quantities are expected to evolve differently, resulting in a large scatter between the relative [O III]+Hβ and UV emission. If L([[O III]+Hβ)/L_UV decreases systematically towards fainter galaxies due to an increased fraction of galaxies with decreasing star formation histories, this could thus in principle explain the relative flattening of the nebular LF with respect to the UV one.

In Fig. 9, we show the [O III]+Hβ-to-UV luminosity ratio over M_UV for our two redshift bins. We find that the ratio slowly evolves with M_UV, with roughly one order of magnitude change over the whole M_UV range. We checked the correlation of the data using a spearman test. For both redshifts bins, the correlation is > 3σ-significant (p-values of 8.14 × 10⁻⁹ and 7.07 × 10⁻⁵ for 7 < z < 8 and 8 < z < 9). Then, to test the effect that such an evolution would have on the LF, and following observations from Fig. 9, we consider a log-linear evolution from log₁₀(L_{[O III]+Hβ}/L_UV) = −2 to −3 between M_UV = −22 to − 12 (R = −0.1M_UV − 4.2). By sampling 10 million galaxies from a UV LF with an α = −2 and converting them to a [O III]+Hβ LF using the assumed relation, we obtain a [O III]+Hβ LF with α ∼ −1.89 (Fig. 10). This shows that the evolution of the [O III]+Hβ-to-UV ratio has the effect of flattening the [O III]+Hβ LF, which could explain the difference between our results and some of the observed UV LF. However, most of the UV LFs are steeper than α = −2 in Fig. 8 and require more than a bursty star formation.

Fig. 9.

Evolution of the [O III]+Hβ-to-UV ratio with M_UV for the 7 < z < 8 sample (left panel) and the 8 < z < 9 sample (right panel). In black is the FRESCO data from Meyer et al. (2024) with depth limit in dashed line, to cover the bright end. The dotted line shows the assumed evolution of the ratio following R = −0.1 × M_UV − 4.2, which enables us to flatten an nebular emission line LF from an intrinsically steep LF, as shown in Fig. 10. This law is arbitrary and chosen to follow the data in a simple manner. However, it does follow closely the linear fit of the data.

4.4.2. The impact of metallicity on the [O III]+Hβ LF

Since Hβ and other H recombination lines are directly related to the number of ionising photons absorbed by the neutral hydrogen (therefore, to the SFR), but forbidden metal lines such as [O III] also strongly depend on other physical parameters (metallicity and ionisation parameter in particular), it is important to separate their relative contribution to the measurement of [O III]+Hβ LF reported here. Ideally, this would allow us to infer the pure Hβ LF, providing a cleaner quantity to compare with the UV LF, which, like Hβ also samples massive stars (although on longer timescales). To achieve this, we now examine the impact of systematic variations of [O III]λ5008/Hβ as a function of galaxy luminosity.

Indeed, it is well known that R3 = [O III]λ5008/Hβ is on average a bimodal function of metallicity which peaks at 12 + log(O/H)∼8, with a maximum R3 ∼ 6 (or ([O III]λλ4960,5008)/Hβ ∼ 8), and decreases at higher and lower metallicities (e.g. Curti et al. 2017; Maiolino & Mannucci 2019; Nakajima et al. 2022; Sanders et al. 2024; Scholte et al. 2025). The relation between M_UV and metallicity is complex, but emerges from the mass-metallicity relation, which states that lower-mass galaxies have lower metallicity (e.g. Chemerynska et al. 2024a; Curti et al. 2020, 2024; Nakajima et al. 2024). We therefore expect the metallicity to decrease to fainter M_UV (Tremonti et al. 2004; Laseter et al. 2025). Unfortunately, direct measurements for galaxies as faint as those in our sample require deep medium-resolution and high signal-to-noise spectroscopy. Chemerynska et al. (2024a) studied the dependency of R3 on stellar mass for very faint galaxies at redshift z ∼ 6 − 7, by measuring the metallicities of 8 galaxies to be 12 + log(O/H)∼6.70 − 7.56 with R3 ∼ 1.4 − 5.5 at M_UV ∼ −15.34 to −17.17. The evolution of 12 + log(O/H) follows the expected decrease in metallicity with M_UV and stellar mass. The Chemerynska et al. (2024a) sample required ultra deep spectroscopy and is already very challenging to observe. However, GLIMPSE goes deeper than M_UV > −15 with photometry, which is vastly unexplored at this redshift and might yield lower R3, as shown by detections or tentative detections of galaxies with R3 < 1 and M_UV ranging from −11 to −16 (Vanzella et al. 2023; Fujimoto et al. 2025; Morishita et al. 2025; Hsiao et al. 2025). Such galaxies require higher magnifications and difficult to reach exposure time, making it challenging to statistically analyse them. Therefore, deep photometric studies such as GLIMPSE have their limitation, but enables indirect probes of the effect of metallicity on very faint (M_UV > −15) galaxies. Regarding the brighter end of [O III]+Hβ fluxes, studies such as Meyer et al. (2024) measured an extreme R3 = 6.38 ± 0.85 line ratio and M_UV ∼ −19.65 for their median stack, close to the maximum observed by Curti et al. (2017), Maiolino & Mannucci (2019). Such values were also confirmed by other JWST surveys (Nakajima et al. 2022; Sanders et al. 2024; Scholte et al. 2025).

To account for systematic variations of the [O III]-to-Hβ ratio with M_UV, we define the following simple function:

$\begin{matrix} R 3 (M_{UV}) = {\begin{matrix} 6.34 & if M_{UV} < - 19 \\ \frac{- 1.8 \times M_{UV} - 25.7}{1.34} & if M_{UV} \in [- 19, - 16.5] \\ \frac{- 0.5 \times M_{UV} - 4.25}{1.34} & if M_{UV} \in [- 16.5, - 12.5] \\ 1.49 & if M_{UV} \geq - 12.5 \end{matrix} \end{matrix}$ $Mathematical equation: $$ \begin{aligned} R3(M_{\rm UV}) = {\left\{ \begin{array}{ll} 6.34&\text{ if} M_{\rm UV} < -19\\ \frac{-1.8 \times M_{\rm UV} - 25.7}{1.34}&\text{ if} M_{\rm UV}\in [-19, -16.5]\\ \frac{-0.5 \times M_{\rm UV} - 4.25}{1.34}&\text{ if} M_{\rm UV}\in [-16.5, -12.5] \\ 1.49&\text{ if} M_{\rm UV} \ge -12.5 \end{array}\right.} \end{aligned} $$$ (12)

where R3(M_UV) varies between 1.49 and 6.34 for M_UV = −12.5 to −19. We used the measurements of Chemerynska et al. (2024a) and Meyer et al. (2024) as the foundation for this function, and we extrapolated for fainter galaxies to account for our measurements. We converted our ratio into the more common ratio following ([O III]λλ4960,5008)/Hβ = 1.34 × [O III]λ5008/Hβ. Therefore, at the bright end (M_UV ∼ −19), the high R3 increases the contribution of [O III]λ5008 (∼66%) to the total [O III]+Hβ relative to Hβ (∼10%). However, at the faint end (M_UV ∼ −12), the R3 ratio becomes so small that Hβ (∼33%) contributes comparably to [O III]+Hβ as [O III]λ5008 (∼50%). This evolving R3 ratio with M_UV predicts a steeper Hβ LF and flatter [O III] LF once Hβ and [O III]λ5008 are separated. By applying this separation to our LF from Fig. 5, we can split the contribution of Hβ and [O III]λ5008 as shown in Appendix C.1. These LFs confirm our assumptions and are parametrised by $α_{H β} = - 1 . 95_{- 0.08}^{+ 0.08}$ $Mathematical equation: $ \alpha_{\mathrm{H}\beta} = -1.95_{-0.08}^{+0.08} $$ and $- 1 . 68_{- 0.14}^{+ 0.13}$ $Mathematical equation: $ -1.68_{-0.14}^{+0.13} $$ , and $α_{[O I I I] λ 5008} = - 1 . 66_{- 0.05}^{+ 0.05}$ $Mathematical equation: $ \alpha_{[\mathrm{O}{\small { {\text{ III}}}}]\lambda5008} = -1.66_{-0.05}^{+0.05} $$ and $- 1 . 45_{- 0.10}^{+ 0.09}$ $Mathematical equation: $ -1.45_{-0.10}^{+0.09} $$ for redshift bins ⟨z⟩ ∼ 7.56 and ⟨z⟩ ∼ 8.45, as listed in Table 1. In short, we now obtain an Hβ LF whose faint-end slope is closer to, but still flatter than, that of the UV LFs at the same redshift (α_UV ∼ −2 to − 2.2 at z ∼ 8).

4.4.3. The early signs of an UV LF turnover

Finally, we explore the third possibility, namely that the flatter [O III]+Hβ LF (compared to the UV LF) could result from an intrinsic turnover of the galaxy UV LF at some faint limit M_UV ≥ −15.5 (Bouwens et al. 2022).

Previous studies investigated the possibility of finding turnovers in LFs to constrain the formation of the lowest mass haloes. On the theory side, a faint-end turnover of the UV LF is expected as a direct consequence of the photoevaporation of small dark matter haloes during the re-ionisation (Shapiro et al. 2004). Simulations observed such turnovers (e.g. Kuhlen et al. 2013; Ocvirk et al. 2016) by parametrising LFs with an additional turnover parameter (e.g. Jaacks et al. 2013) and tentative observations of that turnover were later attempted (Bouwens et al. 2017, 2022) on the HFF. This tentative detection of a turnover ended up ruling out a turnover in the UV LF above M_UV ≤ −15.5 for their z = 2 − 9 sample. In more details, at z ∼ 6, they rule out a turnover above M_UV ≤ −14.3. At higher redshift (z ∼ 9), the small number of detected sources back then limited the depth of the LF. Nevertheless, they argued that a turnover above M_UV ≤ −16 can be ruled out, but they would be consistent with a turnover at M_UV ∼ −15 (Bouwens et al. 2022; Atek et al. 2018). According to Fig. 4, this magnitude would roughly translate into a line flux of L_{[O III]+Hβ} ∼ 10⁴⁰ erg s⁻¹, and therefore, as GLIMPSE observes fainter galaxies, a turnover might be observed. This is roughly compatible with our LFs in Fig. 5 flattening, but does not prove the existence of this turnover. However, we can simulate the effect of a sharp UV LF faint-end turnover on the [O III]+Hβ LF. We define this sharp turnover by adding a simple component to Eq. 10. The final LF is given by

$\begin{matrix} ϕ (L) d L = ϕ^{*} {(\frac{L}{L^{*}})}^{α} exp (- \frac{L}{L^{*}}) exp (- {(\frac{L_{T}}{L})}^{δ}) d (\frac{L}{L^{*}}), \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \phi (L)\mathrm{d}L = \phi ^*\left(\frac{L}{L^*}\right)^\alpha \exp \left(-\frac{L}{L^*}\right)\exp \left(-\left(\frac{L_T}{L}\right)^\delta \right)\mathrm{d}\left(\frac{L}{L^*}\right), \end{aligned} $$$ (13)

where L_T is the turnover luminosity and δ turnover parameter for which negative values indicate the presence of a turnover. To simulate the scattering effect between the [O III]+Hβ and M_UV (as seen in Fig. 4), we converted L_UV to L_{[O III]+Hβ} using the following equation:

$\begin{matrix} {log}_{10} L_{[O III] +H β} = {log}_{10} L_{UV} - 2 + N (0, 0.4), \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \log _{10} L_{{[\text O{\small { {\text{ III}}}}]} \text{]+H}\beta} = \log _{10} L_{\rm UV} - 2 + \mathcal{N} (0, 0.4) , \end{aligned} $$$ (14)

where 𝒩(0.4) is a random number sampled from a normal distribution 𝒩(σ), with an additional a factor log₁₀(0.01) = − 2 to convert from UV to [O III]+Hβ, following the approximate ratio from Eq. (10). The choice of σ corresponds to the full-width half-maximum (FWHM) of ∼1 dex from Fig. 4. To show the impact of standard scattering from uncertainties, we sampled a UV LF with arbitrary parameters (as defined in Eq. (13)) and sampled the [O III]+Hβ emission associated with every UV measurement from Eq. (14). We observe that, while the UV shows a sharp faint-end turnover, the [O III]+Hβ LF instead shows a flattening compared to the UV, which can also mimic the observed flattening in Fig. 5.

Fig. 10.

Demonstration of the effects that arise from allowing the nebular emission lines to be fainter than the continuum in the faintest sources. This enables us to flatten an intrinsically steep UV LF (α ∼ −2, in grey) into a flatter nebular emission line LF (α ∼ −1.89, in black). In red is the Schechter least square fit (α = −1.89, log₁₀L^* = 42.09, log₁₀ϕ_* = −3.63) of the simulated [O III]+Hβ LF. The bright end is dashed, as the effect applied yields non-Schechter-like profile on the bright end.

We also note that because of our very simple scattering assumptions, the bright end is also strongly affected. In reality, the uncertainties are greater for faint galaxies than for brighter ones, and we therefore do not expect to see a major difference in the bright end. However, with the depth of GLIMPSE, we do not observe such a sharp turnover. This could be biased by the completeness estimation, as the very faint end is very incomplete, which may dominate the LF estimation. This will be discussed in a subsequent GLIMPSE publication (Atek et al. in prep.). Finally, UV LFs from JWST surveys did not suggest any turnover in UV at these redshifts or higher, and neither at these relatively similar M_UV (Chemerynska et al. 2026, Atek et al. in prep.). We therefore conclude that the data is not compatible with a turnover in the faint end of the LF.

4.5. Possible missing sources with our method

In this last part, we discuss the validity of the completeness correction, in particular the use of the M_UV completeness correction for our sample. As shown by Fig. 4, at a fixed M_UV, galaxies show variations up to ∼1 dex in [O III]+Hβ flux, meaning that in theory, we might be missing sources by only correcting for completeness on M_UV. In essence, this question requires the exploration of the whole range of [O III]+Hβ equivalent widths.

We first focus on weak [O III]+Hβ equivalent widths (e.g. ≲200 − 300 Å). Such emission lines are typically difficult to constrain with broadband photometry, as their effect is hard to distinguish from the continuum and the noise, hence these galaxies are solely detected via their continuum. Therefore, with our LBG selection, weak equivalent width galaxies are bounded to the detection of the continuum and should be correctly accounted for by our completeness estimations.

For stronger [O III]+Hβ equivalent widths (e.g. ≳1000 Å), their effect is visible and significant in the photometry (e.g. Smit et al. 2014). This increases the signal-to-noise of the affected filters and, in theory, increases the detection rate of galaxies. With our LBG selection, these galaxies would be detected, provided that their continuum is strong enough in other filters (see Eq. (1)). However, there might be situations where the continuum is not well enough detected (S/N_F115W < 5 and S/N_F200W < 5; see Eq. (1)), but the filters affected by [O III]+Hβ are still detected thanks to their strong emissions. In that case, the LBG selection does not find these sources, but are in fact not accounted for in our final completeness estimation. Correcting for this effect would require additional exploration of the completeness, by injecting galaxies with realistic [O III]+Hβ at a given M_UV.

The current statistics of z = 7 − 9 galaxies with ultra-faint M_UV is relatively low, making galaxy templates exploratory in this regime. But in addition to that, such UV-faint and strong [O III]+Hβ galaxies might simply be rare in the Universe. Endsley et al. (2024) observed a decrease in the average [O III]+Hβ equivalent width with fainter galaxies, with almost one dex difference between M_UV ∼ −20.1 and − 17.6 at z = 7 − 9. With our sample going down to M_UV ∼ −11, the decline might continue and making these galaxies even rarer. Therefore, the impact of these galaxies on the final LF should be limited.

5. Implications and discussion

We now present and discuss the main implications of our nebular LFs. We first use them to estimate the ionising photon budget produced by our galaxies, and then we infer the cosmic star formation rate at high redshift (z ∼ 7 − 9), where measurements of Hα are not feasible anymore with NIRCam. All of the following will assume the Eq. 12 function for the R3 ratio.

5.1. The ionising photon budget at z ∼ 7 − 9

For the Universe to re-ionise at its measured rate (Bouwens et al. 2015; Mason et al. 2019), galaxies needs to reach some ionising photon-production budget to overcome for recombination. The source of ionising photons comes from multiple objects (Star-forming galaxies and AGNs Robertson et al. 2015, 2013). Focusing on star formation, Balmer series emission lines enable us to estimate the ionising photon-production rate from star-forming galaxies. From the [O III]+Hβ LF, we can infer the Hβ (or Hα LF), which can then be summed to obtain the instantaneous ionising photon-production rate of galaxies, and thus, for a known escape fraction of ionising photons f_esc, the ionising photon emissivity Ṅ_ion from galaxies. It can be measured using the following equation:

$\begin{matrix} {\dot{N}}_{ion} & = \int_{L_{\min}}^{\infty} f_{esc} Q_{ion} (L) ϕ (L) d L \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \dot{N}_{\rm ion}&= \int _{L_{\rm min}}^\infty f_{\rm esc} Q_{\rm ion}(L)\phi (L) \mathrm{d}L\end{aligned} $$$ (15)

$\begin{matrix} = \int_{L_{\min}}^{\infty} f_{esc} \frac{L}{c_{α} (1 - f_{esc})} ϕ (L) d L, \end{matrix}$ $Mathematical equation: $$ \begin{aligned}&= \int _{L_{\rm min}}^\infty f_{\rm esc}\frac{L}{c_\alpha (1-f_{\rm esc})}\phi (L) \mathrm{d}L, \end{aligned} $$$ (16)

where Q_ion(L) is the ionising photon-production rate, which is related to the luminosity L = Q_ionc_α(1 − f_esc), c_α is the ratio between the line emissivity and the total recombination rate for which we adopt c_α = 1.37 × 10⁻¹² erg for an electron temperature of 10⁴ K (Schaerer 2003), f_esc is the escape fraction of ionising photons and ϕ(L) is the LF. In practice, we convert the measured [O III]+Hβ flux to Hβ using the ratio given by Eq. 12, we correct Hβ for dust attenuation (see Sect. 3.5) and then we multiply by 2.86 (Osterbrock & Ferland 2006) to obtain the dust-corrected Hα LF shown in Fig. C.4. To constrain the bright end, we converted the data from Meyer et al. (2024) in the same manner. However, as they reported a negligible dust attenuation in their stacked spectra, we directly converted their [O III]λ5008 data using their stacked ratio R3 ∼ 6.38 ± 0.85. As expected, the correction for dust is overall rather small, and the resulting LFs are therefore similar (see Figs. C.2 and C.4), as also seen from Table 1.

Figure 12 shows the resulting ionising photon emissivity at z ∼ 7 − 9 from GLIMPSE, as a function of the lower integration limit $L_{H α}^{\min}$ $Mathematical equation: $ L^{\mathrm{min}}_{\mathrm{H}\alpha} $$ , and for a standard fixed f_esc = 0.14 scenario (Jecmen et al. 2026). Probably the most striking result is that, since our faint-end slope is α > −2, the cumulative emissivity becomes already fairly flat over the range of Hα luminosities probed by GLIMPSE. For example, for an integration lower limit L_Hα = 10³⁹ erg s⁻¹, the ionising photon emissivity reaches ${log}_{10} {\dot{N}}_{ion} = 50 . 50_{- 0.07}^{+ 0.07}$ $Mathematical equation: $ \log_{10}\dot{N}_{\mathrm{ion}} = 50.50_{-0.07}^{+0.07} $$ and ${log}_{10} {\dot{N}}_{ion} = 50 . 26_{- 0.08}^{+ 0.10}$ $Mathematical equation: $ \log_{10}\dot{N}_{\mathrm{ion}} = 50.26_{-0.08}^{+0.10} $$ for our two redshift bins (7 < z < 8 and 8 < z < 9), which corresponds to 44%–58% and 36%–46% of the total re-ionisation ionising photon budget measured by Bouwens et al. (2015). However, for the more recent estimations of Ṅ_ion by Mason et al. (2019), our galaxies account for 31%−90% and 46%−156% of the budget for our two redshift bins (7 < z < 8 and 8 < z < 9), implying that star-forming galaxies provide sufficient ionising photons to re-ionise the Universe. Adopting a lower integration limit, e.g. L_Hα ∼ 10³⁸ erg s⁻¹, increases the total contribution of ionising photons only by about 1% – 10%. This means that, for a scenario with constant f_esc, the population of galaxies with very low star formation rates SFR≲(0.005 − 0.001) M_⊙ yr⁻¹, below the detection limit of GLIMPSE, contributes only a small fraction of the total production of ionising photons. The inferred values of Ṅ_ion exhibit variations depending on the assumed constant escape fraction, f_esc, since Ṅ_ion scales with $\frac{f_{esc}}{1 - f_{esc}}$ $Mathematical equation: $ \frac{f_{\mathrm{esc}}}{1-f_{\mathrm{esc}}} $$ . Specifically, Ṅ_ion is 0.5 dex lower for f_esc = 5% or 0.2 dex higher for f_esc = 20%.

According to our earlier scenario, the flattening of the [O III]+Hβ LF is driven by a population of faint galaxies characterised by bursty star formation history and lower metallicities. In this picture, many of these faint galaxies are not actively forming stars, thus diminishing their contribution to the ionising photon budget for re-ionisation. Consequently, the inferred Ṅ_ion stabilises at a fixed value.

To increase the contribution of fainter galaxies, scenarios involving varying escape fraction f_esc enables more ionising photons to escape and thus enhancing their role in cosmic re-ionisation. Such scenarios could emerge from the variation of f_esc with galaxy mass (e.g. Naidu et al. 2022; Begley et al. 2022; Flury et al. 2022; Saldana-Lopez et al. 2023; Pahl et al. 2023), metallicity (e.g. for dwarf irregulars Ramambason et al. 2022; Hunter et al. 2024), SFR (e.g. Giovinazzo et al. 2026), β-slope (e.g. Chisholm et al. 2022; Giovinazzo et al. 2026) or else. We refer to the GLIMPSE paper by Jecmen et al. (2026), which studies the f_esc in more detail.

For comparison, we also display in Fig. 12 the ionising photon emissivity inferred from the Hα LF, adopting α ∼ −2 to − 2.2 (And the same ϕ^* and L^* as the Hα LF). This slope corresponds to the range measured for some UV LF (see Sect. 4.3). Due to the shape of the Schechter function, when the faint-end slope α is smaller than −2, the ionising photon emissivity does not converge; instead, it continues to increase. This leads to a rapid overshooting of the ionising photon budget due to the enhanced contribution of fainter galaxies in the cosmic re-ionisation in that scenario. This problem has earlier been referred to as the ionising photon crisis, which arise by measuring the ionising photon-production rate from the UV LF. This requires the knowledge of the ionising photon-production efficiency ξ_ion. By assuming a high value for ξ_ion and an increase in f_esc for UV-faint galaxies, the Universe re-ionises much too rapidly (see Muñoz et al. 2024).

Our finding of a flattening of the LF therefore argues in favour of their third solution, where lower-mass galaxies have a diminishing impact on the ionising photon-production rate, in contrast to scenarios of low-mass galaxy driven re-ionisation proposed by other studies (e.g. Simmonds et al. 2024a; Atek et al. 2024). By computing ξ_ion = Q_ion/L_UV from our Hα and M_UV measurements, we obtain ${log}_{10} ξ_{ion} = 25 . 31_{- 0.62}^{+ 0.45}$ $Mathematical equation: $ \log_{10}\xi_{\mathrm{ion}} = 25.31_{-0.62}^{+0.45} $$ Hz erg⁻¹ for ⟨z⟩ ∼ 7.48 and $25 . 31_{- 0.59}^{+ 0.45}$ $Mathematical equation: $ 25.31_{-0.59}^{+0.45} $$ Hz erg⁻¹ for ⟨z⟩ ∼ 8.33, using f_esc = 0.14. For comparison with other paper, we also provide ${log}_{10} ξ_{ion, 0} = {log}_{10} ξ_{ion} + {log}_{10} (1 - f_{esc}) = 25 . 24_{- 0.61}^{+ 0.45}$ $Mathematical equation: $ \log_{10}\xi_{\mathrm{ion, 0}} = \log_{10}\xi_{\mathrm{ion}} + \log_{10}(1-f_{\mathrm{esc}}) = 25.24_{-0.61}^{+0.45} $$ Hz erg⁻¹ and $25 . 31_{- 0.59}^{+ 0.45}$ $Mathematical equation: $ 25.31_{-0.59}^{+0.45} $$ Hz erg⁻¹. The measurement does not evolve between our two redshift bins, and shows a decrease with increasing M_UV (see Fig. 9). We refer the reader to the GLIMPSE paper by Chisholm et al. (in prep) for a detailed analysis of ξ_ion at z ∼ 6.

This ξ_ion value is lower than measurements from the early-JWST ξ_ion (e.g. Simmonds et al. 2024a; Atek et al. 2024), which measured higher ionising photons production efficiency for fainter galaxies. However, later results from Simmonds et al. (2024b), which included lower-mass galaxies with a JWST photometric analysis, lead to a decrease in the ionising photon-production efficiency ξ_ion, enabling them to reduce the amount of ionising photons produced and not overestimate the ionising photon emissivity Ṅ_ion from Muñoz et al. (2024). Our findings agrees with this later result and does not overshoot the ionising photon emissivity required to re-ionise the Universe.

5.2. The cosmic star formation rate density

The intensity of Balmer series lines such as Hα are closely related to the SFR of the galaxy (Kennicutt & Evans 2012). We obtain the instantaneous SFR using the Hα line flux as follows:

$\begin{matrix} SFR = \frac{L_{H α}}{C_{X}} [M_{⊙} {yr}^{- 1}], \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \text{ SFR} = \frac{L_{H\alpha }}{C_X}\ \text{[}M_\odot \text{ yr}^{-1}\text{]} , \end{aligned} $$$ (17)

with the constant of proportionality log₁₀(C_X) = 41.27 calibrated by Kennicutt & Evans (2012). The calibration was performed with the Kroupa IMF (Kroupa & Weidner 2003), but because of the similarities between Kroupa and Chabrier IMF, the difference is minimal and therefore, we neglect conversion factors. The final SFR LF has the exact same shape as Fig. C.4, with a simple horizontal shift due to conversion to SFR. We included the SFR conversion on top of Fig. C.4.

Finally, to deduce the cosmic star formation rate density (SFRD), we simply integrate the Hα LF between a minimum luminosity value and theoretically infinity using

$\begin{matrix} SFRD = \int_{{SFRD}_{\min}}^{+ \infty} SFR \times Φ (log SFR) d log SFR . \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \text{ SFRD} = \int _{\mathrm{SFRD}_{\rm min}}^{+\infty } \text{ SFR}\times \Phi (\log \text{ SFR})\mathrm{d}\log \text{ SFR} .\end{aligned} $$$ (18)

We numerically integrate the SFR LF between a given lower limit and 10⁵ M_⊙ yr⁻¹ to obtain the SFRD following Eq. (18). Because GLIMPSE goes deeper than previous studies, who set a lower integration limit at SFRD_min ∼ 0.24 − 0.30 M_⊙ yr⁻¹ (hereafter, standard). This corresponds to log₁₀(L_Hα/erg s⁻¹) ∼ 40.75, which only covers the brighter galaxies of GLIMPSE and discard the faint-end. Therefore, we will compare our results with the same integration limit, as well as the deeper limit of 0.005 M_⊙ yr⁻¹ (hereafter: deep, and approximately corresponding to the faint end of the GLIMPSE detections of L_Hα ∼ 10³⁹ erg s⁻¹). In addition, while in principle Eq. (18) must be integrated to infinity, the uncertainties due to unexplored parameter space and numerical errors would artificially increase the measurement uncertainties at the bright end. Therefore, to avoid it, we set a reasonable integration upper limit of 1000 M_⊙ yr⁻¹ (i.e. L_Hα ∼ 10⁴⁴ erg s⁻¹).

In Fig. 13 we observe the redshift evolution of the cosmic star formation rate density (SFRD). We showed multiple literature points, coming from both UV measurements (Oesch et al. 2018; Bouwens et al. 2020) and Hα measurements (Bollo et al. 2023; Covelo-Paz et al. 2025; Fu et al. 2025), as well as two sets of GLIMPSE data points with different integration limits. The bright orange data point corresponds to standard integration limit, while the light orange to a deep integration limit of 0.005 M_⊙ yr⁻¹. GLIMPSE reaches deeper line fluxes, and shows flattening that suggest that fainter galaxies do not significantly contribute to the star formation rate density. GLIMPSE reaches deeper line flux limits and reveals a flattening in the LFs, suggesting that the contribution of fainter galaxies to the star formation rate density is not significant. Therefore, we decided to add a secondary measurement at lower integration limit. Table E.1 gives the SFRD for both integration limits and both redshift ranges, as well as the compilation of literature data points shown in Fig. 13.

Our results with the standard integration limits seem to be inline with an extrapolation of the Hα literature measurement, but it differs from UV measurements. However, both UV measurement comes from pre-JWST surveys, which had limited statistics in their faintest bins. All the nebular emission line measured SFRD from Hα come from JWST results, which is more sensitive to fainter galaxies, allowing us to retrieve the total SFRD better. In addition, UV and nebular emission lines do not trace the same star formation (Kennicutt & Evans 2012), which indicates a burstier star formation in the early Universe. Therefore, GLIMPSE does show a slightly higher total SFRD compared to previous UV-focus SFRD surveys, but inline for nebular emission lines and the compilation of Madau & Dickinson (2014).

Between the two integration limits, we observe a mild increase in the SFRD for the deeper integration limit, of 0.3dex (∼ × 1.99) and 0.1dex (∼ × 1.26) for ⟨z⟩ ∼ 7.48 and ⟨z⟩ ∼ 8.33 respectfully. This increase is only significant for the ⟨z⟩ ∼ 7.48 redshift bin, which shows an increasing importance of lower star-forming galaxies to the total SFRD with lowering redshift. This nevertheless agrees with the result showed in the previous section, with the fainter galaxies having a more limited impact on the total SFRD, with an addition of only a fraction of dex of SFRD despite the almost 2 orders of magnitude deeper data.

6. Conclusions

We measured the [O III]+Hβ LF for the strongly lensed field Abell S1063 with very deep JWST GLIMPSE observations. We selected our sample of 164 unique galaxies between redshift 7 < z < 9 using an LBG selection in combination with photometric redshifts obtained via SED fitting. We constrained the strongly lensed field using very deep HST and JWST images, removed the multiple counter images of some galaxies, and measured the completeness of our sample. We used an SED fitting to measure the [O III]+Hβ flux for each galaxy, the dust attenuation, and other physical quantities, and we then constructed the [O III]+Hβ LF for that sample. Combined with existing measurements at the bright end from Meyer et al. (2024), we thus determined the [O III]+Hβ LF down to unprecedented luminosities, reaching L_{[O III]+Hβ} ≥ 10³⁹ erg s⁻¹ (equivalent to galaxies as faint as M_UV ∼ −12; the average completeness at these magnitudes is low, but the bins have multiple detections due to magnification). The main results of our study are summarised below.

The [O III]+Hβ LFs show a relatively flat faint-end slope α ∼ −1.55 to − 1.78 at z ∼ 7 − 9, which contrasts with previous LFs. Studies such as Wold et al. (2025) already paved the way into fainter LFs by combining photometry of Abell 2744 and its lensing power, which enabled them to reach L_{[O III]+Hβ} ≥ 10⁴¹ erg s⁻¹. They measured a faint-end slope $α = - 2 . 07_{- 0.23}^{+ 0.22}$ $Mathematical equation: $ \alpha=-2.07_{-0.23}^{+0.22} $$ , which is significantly steeper than that of GLIMPSE. By fitting the GLIMPSE LF with a similar luminosity limit, we measured faint-end slopes α comparable to their results. On the other hand, spectroscopy-based LFs have more limited access to the faintest galaxies, with the best studies reaching L_{[O III]+Hβ} ≥ 10^41.75 erg s⁻¹. (e.g. Meyer et al. 2024; Matthee et al. 2023; Sun et al. 2023; De Barros et al. 2019). Our best comparison is with Meyer et al. (2024), but their sample typically probed the faint end of the LF only in a limited way, and they fixed the faint-end slope to steep values of α < −2.
This result contrasts with the steeper slope observed in UV LFs (α ≲ −2) generally found at z ≳ 7. To explain these differences between LFs of nebular emission lines and the UV continuum, we examined three scenarios: (i) A decrease in the [O III]+Hβ-to-UV ratio towards fainter galaxies (Fig. 4), which might be caused by more UV-faint (i.e. lower-mass) galaxies seen in a downturn phase (Endsley et al. 2025). (ii) A decrease in R3 = [O III]λ5008/Hβ towards fainter galaxies due to the decreasing metallicity, as observed by recent studies (e.g. Chemerynska et al. 2024a; Meyer et al. 2024), which implies that the Hβ and the [O III]λ5008 LFs have different faint-end slopes (Fig. 10). (iii) A possible sharp turnover of the UV-LF at the faint end below the current detection limits, which, due to natural scattering, would be observed as a simple flattening of the LF (Fig. 11). As no flatteing of the UV LF has so far been found, but burstiness and the evolution of R3=[O III]λ5008/Hβ have both been observed, we favour the first two solutions (i and ii) to explain the flatter nebular LF compared to the UV LF.
Assuming an average relation between the R3 ratio and M_UV (Eq. 12), we separated the contribution of Hβ and [O III]λ5008 to the observed [O III]+Hβ, leading to a steeper Hβ LF (α ∼ −1.68 to − 1.95; Fig. C.2) and a flatter [O III]λ5008 LF (α ∼ −1.45 to − 1.66; Fig. C.3). Therefore, the Hβ LF approaches but remains flatter than the UV LF (α ≤ −2), while the [O III]λ5008 LF quickly flattens out.
Correcting for dust attenuation (measured through an SED fitting), we converted the Hβ LF into a dust-corrected Hα LF (Fig. C.4), which allowed us to compute the total ionising photon emissivity and the cosmic SFR density of galaxies at z ∼ 7 − 9 (Fig. 12). Since the slope α_nebular > −2, we found that the total photon emissivity quickly saturates and does not significantly increase when it is integrated to objects fainter than those observed with GLIMPSE. Assuming a standard constant f_esc = 0.14, we reached 31%−90% and 46%−156% of the ionising photon budget required to drive re-ionisation at z ∼ 7.5 and z ∼ 8.3, respectively, as measured by Mason et al. (2019).
The SFRD (Fig. 13) showed comparable results to previous studies for comparable lower integration limits (∼ 0.3 M_⊙ yr⁻¹). Integrating down to our observed limit of SFR ∼ 0.005 M_⊙ yr⁻¹ (corresponding to L_Hα ∼ 10³⁹ erg s⁻¹), we obtained an SFRD higher by ∼0.1 − 0.3 dex, where the increase is only significant at ⟨z⟩ ∼ 7.48. The inclusion of fainter, lower SFR galaxies will not strongly increase the cosmic SFRD at these redshifts because the observed nebular LF, which provides the most robust measure of the instantaneous SFRD, is significantly flatter than the UV LF.

Fig. 11.

Simulation of the effect of scattering on a sharp faint-end turnover in a Schechter LF. We used the (B22, Bouwens et al. 2022) Schechter parametrisation at z ∼ 8, with a sharper turnover δ: L^* = 10^45.26 erg s⁻¹(M_UV = −20.9), α = −2.2, δ = 2.0, L_T = 10^42.98 erg s⁻¹(M_UV = −15.2). We include scatter that correspond to the natural scatter between [O III]+Hβ and M_UV (as seen in Fig. 4) by sampling from a normal distribution with parameters 𝒩(σ = 0.4).

Fig. 12.

Ionising photon emissivity as a function of the minimum Hα luminosity $L_{H α}^{\min}$ $Mathematical equation: $ L_{\mathrm{H}\alpha}^{\mathrm{min}} $$ as extrapolated from our observed [O III]+Hβ LF (see Sect. 5.1) in the two redshift bins we investigated. The ionising photon emissivity required to re-ionise the Universe at its measured rate, taken from Bouwens et al. (2015) is shown by the shaded areas and from Mason et al. (2019) in the hatched areas. All of our measurements assume f_esc = 0.14 (Jecmen et al. 2026), the solid lines are extrapolated from this work and the dotted lines assumes α = −2 to − 2.2, which is similar to the measured faint-end slope of the typical UV LF at this redshift. The black line shows the GLIMPSE depth following the data from Fig. C.4.

Fig. 13.

Evolution over redshift of the star formation rate density. We used the same colour scheme as Fig. 8, where grey corresponds to measurements based on UV, blue to Hα and orange to [O III] and [O III]+Hβ. The black line corresponds to the results of the review Madau & Dickinson (2014), which compiles the state-of-the-art measurements of the time. All the results have an integration lower limit of 0.24 − 0.30 M_⊙ yr⁻¹ Mpc⁻³ and are IMF-corrected to Chabrier (2003). The z ∼ 7 − 9 orange points correspond to GLIMPSE measurements, where the bright orange uses the standard 0.3 M_⊙ yr⁻¹ Mpc⁻³ lower integration limit and the light orange uses the deep 0.005 M_⊙ yr⁻¹ lower integration limit (approximately corresponding to L_Hα ∼ 10³⁹ erg s⁻¹). The non-exhaustive compilation of UV, Hα and [O III]+Hβ based measurements, including this work, are listed in Table E.1.

In short, the determination of the LF of nebular emission ([O III]+Hβ) from star-forming galaxies down to very faint fluxes yields flatter slopes (α ∼ −1.55 to − 1.78) of the LF than previously thought and measured for the UV LF at high-z (z ≳ 7). This indicates that GLIMPSE, combining ultra-deep JWST observations and strong gravitational lensing, has allowed us to reach the bulk of the star-forming galaxies at z ∼ 7 − 9, and thus, to determine their total ionising photon emissivity and the total cosmic SFRD of these objects. Our results suggest that faint galaxies contribute less to cosmic re-ionisation than previously thought because their number density flattens too rapidly to maintain a highly ionising photon-production rate. With conventional or plausible assumptions on the escape fraction of ionising photons, the observed galaxies are capable of driving cosmic re-ionisation at z ∼ 7 − 9.

Acknowledgments

DK thanks Emma Giovinazzo, Andrea Weibel, Callum Witten, Max Briel, Mengyuan Xiao for their useful discussions and suggestions throughout this project. 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 program #03293. Support for program #03293 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. IC acknowledges funding support from the Initiative Physique des Infinis (IPI), a research training program of the Idex SUPER at Sorbonne Université. LJF and AZ acknowledge support by Grant No. 2020750 from the United States-Israel Binational Science Foundation (BSF) and Grant No. 2109066 from the United States National Science Foundation (NSF); by the Ministry of Science & Technology, Israel; and by the Israel Science Foundation Grant No. 864/23. HA and IC acknowledge support from CNES, focused on the JWST mission and the Programme National Cosmology and Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. HA is supported by the French National Research Agency (ANR) under grant ANR-21-CE31-0838. RAM acknowledges support from the Swiss National Science Foundation (SNSF) through project grant 200020_207349. ASL acknowledges support from Knut and Alice Wallenberg Foundation. AA acknowledges support by the Swedish research council Vetenskapsrådet (VR 2021-05559, and VR consolidator grant 2024-02061). JBM acknowledges support from NSF Grants AST-2307354 and AST-2408637. This work has received funding from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00072, as well as from the Swiss National Science Foundation (SNSF) through project grant 200020_207349. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. Softwares: Emcee (Foreman-Mackey et al. 2013), Astropy (Collaboration et al. 2022), Numpy (Harris et al. 2020), Scipy (Virtanen et al. 2020), Matplotlib (Hunter 2007), Seaborn (Waskom 2021), Photutils (Bradley et al. 2024), LMFIT (Newville et al. 2025).

References

Adamo, A., Atek, H., Bagley, M. B., et al. 2025, Nat. Astron., 9, 1134 [Google Scholar]
Atek, H., Siana, B., Scarlata, C., et al. 2011, ApJ, 743, 121 [NASA ADS] [CrossRef] [Google Scholar]
Atek, H., Richard, J., Kneib, J.-P., & Schaerer, D. 2018, MNRAS, 479, 5184 [Google Scholar]
Atek, H., Labbé, I., Furtak, L. J., et al. 2024, Nature, 626, 975 [NASA ADS] [CrossRef] [Google Scholar]
Atek, H., Chisholm, J., Kokorev, V., et al. 2025 ArXiv e-prints [arXiv:2511.07542] [Google Scholar]
Bastiaansen, P. A. 1992, A&AS, 93, 449 [NASA ADS] [Google Scholar]
Beauchesne, B., Clément, B., Hibon, P., et al. 2024, MNRAS, 527, 3246 [Google Scholar]
Becker, R. H., Fan, X., White, R. L., et al. 2001, AJ, 122, 2850 [NASA ADS] [CrossRef] [Google Scholar]
Begley, R., Cullen, F., McLure, R. J., et al. 2022, MNRAS, 513, 3510 [NASA ADS] [CrossRef] [Google Scholar]
Begley, R., McLure, R. J., Cullen, F., et al. 2025, MNRAS, 537, 3245 [NASA ADS] [CrossRef] [Google Scholar]
Berg, D. A., Chisholm, J., Erb, D. K., et al. 2021, ApJ, 922, 170 [CrossRef] [Google Scholar]
Bergamini, P., Rosati, P., Mercurio, A., et al. 2019, A&A, 631, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bezanson, R., Labbe, I., Whitaker, K. E., et al. 2024, ApJ, 974, 92 [NASA ADS] [CrossRef] [Google Scholar]
Bollo, V., González, V., Stefanon, M., et al. 2023, ApJ, 946, 117 [NASA ADS] [CrossRef] [Google Scholar]
Bongiovanni, Á., Ramón-Pérez, M., Pérez García, A. M., et al. 2020, A&A, 635, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bosman, S. E. I., & Davies, F. B. 2024, A&A, 690, A391 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bosman, S. E. I., Davies, F. B., Becker, G. D., et al. 2022, MNRAS, 514, 55 [NASA ADS] [CrossRef] [Google Scholar]
Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 811, 140 [NASA ADS] [CrossRef] [Google Scholar]
Bouwens, R. J., Oesch, P. A., Illingworth, G. D., Ellis, R. S., & Stefanon, M. 2017, ApJ, 843, 129 [Google Scholar]
Bouwens, R. J., Stefanon, M., Oesch, P. A., et al. 2019, ApJ, 880, 25 [NASA ADS] [CrossRef] [Google Scholar]
Bouwens, R., González-López, J., Aravena, M., et al. 2020, ApJ, 902, 112 [NASA ADS] [CrossRef] [Google Scholar]
Bouwens, R. J., Illingworth, G., Ellis, R. S., Oesch, P., & Stefanon, M. 2022, ApJ, 940, 55 [NASA ADS] [CrossRef] [Google Scholar]
Bowler, R. A. A., Jarvis, M. J., Dunlop, J. S., et al. 2020, MNRAS, 493, 2059 [Google Scholar]
Bowman, W. P., Ciardullo, R., Zeimann, G. R., et al. 2021, ApJ, 920, 78 [Google Scholar]
Boyett, K., Bunker, A. J., Curtis-Lake, E., et al. 2024, MNRAS, 535, 1796 [NASA ADS] [CrossRef] [Google Scholar]
Bradley, L., Sipőcz, B., Robitaille, T., et al. 2024, https://doi.org/10.5281/zenodo.13989456 [Google Scholar]
Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503 [Google Scholar]
Bromm, V. 2013, Rep. Prog. Phys., 76, 112901 [NASA ADS] [CrossRef] [Google Scholar]
Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682 [NASA ADS] [CrossRef] [Google Scholar]
Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245 [Google Scholar]
Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
Chemerynska, I., Atek, H., Dayal, P., et al. 2024a, ApJ, 976, L15 [NASA ADS] [CrossRef] [Google Scholar]
Chemerynska, I., Atek, H., Furtak, L. J., et al. 2024b, MNRAS, 531, 2615 [NASA ADS] [CrossRef] [Google Scholar]
Chemerynska, I., Atek, H., Furtak, L. J., et al. 2026, MNRAS, 546, staf2267 [Google Scholar]
Chevallard, J., & Charlot, S. 2016, MNRAS, 462, 1415 [NASA ADS] [CrossRef] [Google Scholar]
Chisholm, J., Saldana-Lopez, A., Flury, S., et al. 2022, MNRAS, 517, 5104 [CrossRef] [Google Scholar]
Colbert, J. W., Teplitz, H., Atek, H., et al. 2013, ApJ, 779, 34 [NASA ADS] [CrossRef] [Google Scholar]
Collaboration, A., Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
Covelo-Paz, A., Giovinazzo, E., Oesch, P. A., et al. 2025, A&A, 694, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Covelo-Paz, A., Meuwly, C., Oesch, P. A., et al. 2026, A&A, 705, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Curti, M., Cresci, G., Mannucci, F., et al. 2017, MNRAS, 465, 1384 [Google Scholar]
Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944 [Google Scholar]
Curti, M., Maiolino, R., Curtis-Lake, E., et al. 2024, A&A, 684, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nat. Astron., 7, 622 [NASA ADS] [CrossRef] [Google Scholar]
Daikuhara, K., Morishita, T., Kodama, T., et al. 2025, ApJ, 989, 71 [Google Scholar]
D’Aloisio, A., McQuinn, M., & Trac, H. 2015, ApJ, 813, L38 [CrossRef] [Google Scholar]
Davis, K., Trump, J. R., Simons, R. C., et al. 2024, ApJ, 974, 42 [NASA ADS] [CrossRef] [Google Scholar]
Dayal, P., Volonteri, M., Choudhury, T. R., et al. 2020, MNRAS, 495, 3065 [Google Scholar]
De Barros, S., Oesch, P. A., Labbé, I., et al. 2019, MNRAS, 489, 2355 [NASA ADS] [CrossRef] [Google Scholar]
Dome, T., Tacchella, S., Fialkov, A., et al. 2024, MNRAS, 527, 2139 [Google Scholar]
Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011 [Google Scholar]
Donnan, C. T., McLure, R. J., Dunlop, J. S., et al. 2024, MNRAS, 533, 3222 [NASA ADS] [CrossRef] [Google Scholar]
Dunlop, J. S., & Peacock, J. A. 1990, MNRAS, 247, 19 [NASA ADS] [Google Scholar]
Eddington, A. S. 1913, MNRAS, 73, 359 [Google Scholar]
Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, ArXiv e-prints [arXiv:2306.02465] [Google Scholar]
Elíasdóttir, Á., Limousin, M., Richard, J., et al. 2007, ArXiv e-prints [arXiv:0710.5636] [Google Scholar]
Ellis, R. S., McLure, R. J., Dunlop, J. S., et al. 2013, ApJ, 763, L7 [NASA ADS] [CrossRef] [Google Scholar]
Endsley, R., Stark, D. P., Whitler, L., et al. 2024, MNRAS, 533, 1111 [NASA ADS] [CrossRef] [Google Scholar]
Endsley, R., Chisholm, J., Stark, D. P., Topping, M. W., & Whitler, L. 2025, ApJ, 987, 189 [Google Scholar]
Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, AJ, 132, 117 [NASA ADS] [CrossRef] [Google Scholar]
Ferrara, A. 1998, ApJ, 499, L17 [Google Scholar]
Ferrara, A., & Pandolfi, S. 2014, ArXiv e-prints [arXiv:1409.4946] [Google Scholar]
Finkelstein, S. L., D’Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36 [Google Scholar]
Finkelstein, S. L., Leung, G. C. K., Bagley, M. B., et al. 2024, ApJ, 969, L2 [NASA ADS] [CrossRef] [Google Scholar]
Finkelstein, S. L., Bagley, M. B., Arrabal Haro, P., et al. 2025, ApJ, 983, L4 [Google Scholar]
Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022, ApJ, 930, 126 [NASA ADS] [CrossRef] [Google Scholar]
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
Fu, S., Sun, F., Jiang, L., et al. 2025, ApJ, 987, 186 [Google Scholar]
Fujimoto, S., Arrabal Haro, P., Dickinson, M., et al. 2023, ApJ, 949, L25 [NASA ADS] [CrossRef] [Google Scholar]
Fujimoto, S., Naidu, R. P., Chisholm, J., et al. 2025, ApJ, 989, 46 [Google Scholar]
Furtak, L. J., Zitrin, A., Weaver, J. R., et al. 2023, MNRAS, 523, 4568 [NASA ADS] [CrossRef] [Google Scholar]
Gelli, V., Salvadori, S., Ferrara, A., Pallottini, A., & Carniani, S. 2023, ApJ, 954, L11 [NASA ADS] [CrossRef] [Google Scholar]
Giovinazzo, E., Oesch, P. A., Weibel, A., et al. 2026, A&A, 707, A352 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Grazian, A., Giallongo, E., Boutsia, K., et al. 2024, ApJ, 974, 84 [Google Scholar]
Harikane, Y., Nakajima, K., Ouchi, M., et al. 2024, ApJ, 960, 56 [NASA ADS] [CrossRef] [Google Scholar]
Harikane, Y., Inoue, A. K., Ellis, R. S., et al. 2025, ApJ, 980, 138 [NASA ADS] [Google Scholar]
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
Hashimoto, T., Álvarez-Márquez, J., Fudamoto, Y., et al. 2023, ApJ, 955, L2 [NASA ADS] [CrossRef] [Google Scholar]
Hayashi, M., Tanaka, M., Shimakawa, R., et al. 2018, PASJ, 70, S17 [NASA ADS] [Google Scholar]
Hsiao, T. Y. Y., Sun, F., Lin, X., et al. 2025, ArXiv e-prints [arXiv:2505.03873] [Google Scholar]
Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
Hunter, D. A., Elmegreen, B. G., & Madden, S. C. 2024, ARA&A, 62, 113 [NASA ADS] [CrossRef] [Google Scholar]
Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805 [NASA ADS] [CrossRef] [Google Scholar]
Izotov, Y. I., Schaerer, D., Worseck, G., et al. 2020, MNRAS, 491, 468 [NASA ADS] [CrossRef] [Google Scholar]
Jaacks, J., Thompson, R., & Nagamine, K. 2013, ApJ, 766, 94 [NASA ADS] [CrossRef] [Google Scholar]
Jecmen, M. C., Chisholm, J., Atek, H., et al. 2026, ArXiv e-prints [arXiv:2601.19995] [Google Scholar]
Kashino, D., Lilly, S. J., Matthee, J., et al. 2023, ApJ, 950, 66 [CrossRef] [Google Scholar]
Kassiola, A., & Kovner, I. 1993, ApJ, 417, 450 [Google Scholar]
Keating, L. C., Weinberger, L. H., Kulkarni, G., et al. 2020, MNRAS, 491, 1736 [NASA ADS] [CrossRef] [Google Scholar]
Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531 [NASA ADS] [CrossRef] [Google Scholar]
Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511 [Google Scholar]
Khostovan, A. A., Sobral, D., Mobasher, B., et al. 2015, MNRAS, 452, 3948 [Google Scholar]
Khostovan, A. A., Malhotra, S., Rhoads, J. E., et al. 2020, MNRAS, 493, 3966 [NASA ADS] [CrossRef] [Google Scholar]
Kokorev, V., Atek, H., Chisholm, J., et al. 2025, ApJ, 983, L22 [Google Scholar]
Korber, D., Bianco, M., Tolley, E., & Kneib, J.-P. 2023, MNRAS, 521, 902 [NASA ADS] [CrossRef] [Google Scholar]
Kroupa, P., & Weidner, C. 2003, ApJ, 598, 1076 [NASA ADS] [CrossRef] [Google Scholar]
Kuhlen, M., Madau, P., & Krumholz, M. R. 2013, ApJ, 776, 34 [NASA ADS] [CrossRef] [Google Scholar]
Laseter, I. H., Maseda, M. V., Simmonds, C., et al. 2025, ApJ, 988, 73 [Google Scholar]
Leitherer, C., Li, I.-H., Calzetti, D., & Heckman, T. M. 2002, ApJS, 140, 303 [Google Scholar]
Llerena, M., Amorín, R., Pentericci, L., et al. 2024, A&A, 691, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2024, Nature, 629, 53 [Google Scholar]
Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2025, A&A, 697, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97 [Google Scholar]
Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
Madau, P., Giallongo, E., Grazian, A., & Haardt, F. 2024, ApJ, 971, 75 [NASA ADS] [CrossRef] [Google Scholar]
Maiolino, R., & Mannucci, F. 2019, A&ARv, 27, 3 [Google Scholar]
Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2024, A&A, 691, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Mascia, S., Pentericci, L., Llerena, M., et al. 2025, A&A, 701, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Mason, C. A., Naidu, R. P., Tacchella, S., & Leja, J. 2019, MNRAS, 489, 2669 [Google Scholar]
Matthee, J., Mackenzie, R., Simcoe, R. A., et al. 2023, ApJ, 950, 67 [NASA ADS] [CrossRef] [Google Scholar]
Meyer, R. A., Oesch, P. A., Giovinazzo, E., et al. 2024, MNRAS, 535, 1067 [CrossRef] [Google Scholar]
Meyer, R. A., Roberts-Borsani, G., Oesch, P. A., & Ellis, R. S. 2025, MNRAS, 542, 1952 [Google Scholar]
Mintz, A., Setton, D. J., Greene, J. E., et al. 2025, ArXiv e-prints [arXiv:2506.16510] [Google Scholar]
Morishita, T., Liu, Z., Stiavelli, M., et al. 2025, ArXiv e-prints [arXiv:2507.10521] [Google Scholar]
Moutard, T., Sawicki, M., Arnouts, S., et al. 2020, MNRAS, 494, 1894 [NASA ADS] [CrossRef] [Google Scholar]
Muñoz, J. B., Mirocha, J., Chisholm, J., Furlanetto, S. R., & Mason, C. 2024, MNRAS, 535, L37 [CrossRef] [Google Scholar]
Nagaraj, G., Ciardullo, R., Bowman, W. P., Lawson, A., & Gronwall, C. 2023, ApJ, 943, 5 [Google Scholar]
Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109 [NASA ADS] [CrossRef] [Google Scholar]
Naidu, R. P., Matthee, J., Oesch, P. A., et al. 2022, MNRAS, 510, 4582 [CrossRef] [Google Scholar]
Naidu, R. P., Oesch, P. A., Brammer, G., et al. 2026, Open J. Astrophys., 9, 56033 [Google Scholar]
Nakajima, K., Ouchi, M., Xu, Y., et al. 2022, ApJS, 262, 3 [CrossRef] [Google Scholar]
Nakajima, K., Ouchi, M., Isobe, Y., et al. 2024, ArXiv e-prints [arXiv:2412.04541] [Google Scholar]
Napolitano, L., Castellano, M., Pentericci, L., et al. 2025, A&A, 693, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Newville, M., Otten, R., Nelson, A., et al. 2025, https://doi.org/10.5281/zenodo.598352 [Google Scholar]
Ocvirk, P., Gillet, N., Shapiro, P. R., et al. 2016, MNRAS, 463, 1462 [Google Scholar]
O’Donnell, J. E. 1994, ApJ, 422, 158 [Google Scholar]
Oesch, P. A., Carollo, C. M., Stiavelli, M., et al. 2009, ApJ, 690, 1350 [NASA ADS] [CrossRef] [Google Scholar]
Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129 [NASA ADS] [CrossRef] [Google Scholar]
Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105 [Google Scholar]
Oesch, P. A., Brammer, G., Naidu, R. P., et al. 2023, MNRAS, 525, 2864 [NASA ADS] [CrossRef] [Google Scholar]
Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713 [NASA ADS] [CrossRef] [Google Scholar]
Onodera, M., Shimakawa, R., Suzuki, T. L., et al. 2020, ApJ, 904, 180 [NASA ADS] [CrossRef] [Google Scholar]
Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei [Google Scholar]
Pahl, A. J., Shapley, A., Steidel, C. C., et al. 2023, MNRAS, 521, 3247 [NASA ADS] [CrossRef] [Google Scholar]
Parsa, S., Dunlop, J. S., McLure, R. J., & Mortlock, A. 2016, MNRAS, 456, 3194 [Google Scholar]
Ramambason, L., Lebouteiller, V., Bik, A., et al. 2022, A&A, 667, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Reddy, N. A., & Steidel, C. C. 2009, ApJ, 692, 778 [Google Scholar]
Richard, J., Stark, D. P., Ellis, R. S., et al. 2008, ApJ, 685, 705 [NASA ADS] [CrossRef] [Google Scholar]
Rieke, M. J., Robertson, B., Tacchella, S., et al. 2023, ApJS, 269, 16 [NASA ADS] [CrossRef] [Google Scholar]
Rinaldi, P., Caputi, K. I., Costantin, L., et al. 2023, ApJ, 952, 143 [NASA ADS] [CrossRef] [Google Scholar]
Rinaldi, P., Navarro-Carrera, R., Caputi, K. I., et al. 2025, ApJ, 981, 161 [Google Scholar]
Robertson, B. E. 2022, ARA&A, 60, 121 [NASA ADS] [CrossRef] [Google Scholar]
Robertson, B. E., Furlanetto, S. R., Schneider, E., et al. 2013, ApJ, 768, 71 [Google Scholar]
Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJ, 802, L19 [Google Scholar]
Saldana-Lopez, A., Schaerer, D., Chisholm, J., et al. 2023, MNRAS, 522, 6295 [NASA ADS] [CrossRef] [Google Scholar]
Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2024, ApJ, 962, 24 [NASA ADS] [CrossRef] [Google Scholar]
Schaerer, D. 2003, A&A, 397, 527 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schaerer, D., & de Barros, S. 2009, A&A, 502, 423 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schechter, P. 1976, ApJ, 203, 297 [Google Scholar]
Schenker, M. A., Stark, D. P., Ellis, R. S., et al. 2012, ApJ, 744, 179 [Google Scholar]
Scholte, D., Cullen, F., Carnall, A. C., et al. 2025, MNRAS, 540, 1800 [Google Scholar]
Shapiro, P. R., Iliev, I. T., & Raga, A. C. 2004, MNRAS, 348, 753 [NASA ADS] [CrossRef] [Google Scholar]
Shipley, H. V., Lange-Vagle, D., Marchesini, D., et al. 2018, ApJS, 235, 14 [Google Scholar]
Simmonds, C., Tacchella, S., Hainline, K., et al. 2024a, MNRAS, 527, 6139 [Google Scholar]
Simmonds, C., Tacchella, S., Hainline, K., et al. 2024b, MNRAS, 535, 2998 [NASA ADS] [CrossRef] [Google Scholar]
Singha, M., Malhotra, S., & Ely Rhoads, J. 2025, ArXiv e-prints [arXiv:2510.26990] [Google Scholar]
Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24 [Google Scholar]
Smit, R., Bouwens, R. J., Labbé, I., et al. 2014, ApJ, 784, 58 [NASA ADS] [CrossRef] [Google Scholar]
Sobral, D., Smail, I., Best, P. N., et al. 2013, MNRAS, 428, 1128 [NASA ADS] [CrossRef] [Google Scholar]
Stark, D. P., Ellis, R. S., & Ouchi, M. 2011, ApJ, 728, L2 [Google Scholar]
Steinhardt, C. L., Jauzac, M., Acebron, A., et al. 2020, ApJS, 247, 64 [Google Scholar]
Storey, P. J., & Zeippen, C. J. 2000, MNRAS, 312, 813 [NASA ADS] [CrossRef] [Google Scholar]
Strait, V., Brammer, G., Muzzin, A., et al. 2023, ApJ, 949, L23 [NASA ADS] [CrossRef] [Google Scholar]
Suess, K. A., Weaver, J. R., Price, S. H., et al. 2024, ApJ, 976, 101 [Google Scholar]
Sun, F., Egami, E., Pirzkal, N., et al. 2023, ApJ, 953, 53 [NASA ADS] [CrossRef] [Google Scholar]
Sun, L., Wang, X., Teplitz, H. I., et al. 2024, ApJ, 972, 8 [Google Scholar]
Tang, M., Stark, D. P., Chevallard, J., & Charlot, S. 2019, MNRAS, 489, 2572 [NASA ADS] [CrossRef] [Google Scholar]
Topping, M. W., Stark, D. P., Senchyna, P., et al. 2024, MNRAS, 529, 3301 [NASA ADS] [CrossRef] [Google Scholar]
Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898 [Google Scholar]
Trussler, J. A. A., Conselice, C. J., Adams, N., et al. 2025, MNRAS, 537, 3662 [Google Scholar]
van der Wel, A., Straughn, A. N., Rix, H.-W., et al. 2011, ApJ, 742, 111 [NASA ADS] [CrossRef] [Google Scholar]
Vanzella, E., Loiacono, F., Bergamini, P., et al. 2023, A&A, 678, A173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Med., 17, 261 [Google Scholar]
Wang, B., Fujimoto, S., Labbé, I., et al. 2023, ApJ, 957, L34 [NASA ADS] [CrossRef] [Google Scholar]
Waskom, M. L. 2021, J. Open Source Software, 6, 3021 [NASA ADS] [CrossRef] [Google Scholar]
Weaver, J. R., Cutler, S. E., Pan, R., et al. 2024, ApJS, 270, 7 [NASA ADS] [CrossRef] [Google Scholar]
Weibel, A., Oesch, P. A., Williams, C. C., et al. 2025, ArXiv e-prints [arXiv:2507.06292] [Google Scholar]
Williams, C. C., Curtis-Lake, E., Hainline, K. N., et al. 2018, ApJS, 236, 33 [CrossRef] [Google Scholar]
Willott, C. J., Desprez, G., Asada, Y., et al. 2024, ApJ, 966, 74 [NASA ADS] [CrossRef] [Google Scholar]
Wold, I. G. B., Malhotra, S., Rhoads, J. E., Weaver, J. R., & Wang, B. 2025, ApJ, 980, 200 [Google Scholar]
Yeh, J. Y.-C., Smith, A., Kannan, R., et al. 2023, MNRAS, 520, 2757 [NASA ADS] [CrossRef] [Google Scholar]
Zaroubi, S. 2013, Astrophys. Space Sci. Libr., 396, 45 [NASA ADS] [CrossRef] [Google Scholar]
Zitrin, A., Fabris, A., Merten, J., et al. 2015a, ApJ, 801, 44 [Google Scholar]
Zitrin, A., Labbé, I., Belli, S., et al. 2015b, ApJ, 810, L12 [Google Scholar]

Appendix A: Validating the SED fitting flux and EW measurement using flux excess

To assess the quality of CIGALE measurements on the GLIMPSE dataset, we measure an estimation of the [O III]+Hβ line flux and the equivalent width according to the excesses measured using the F444W, F410M and F480M filters. We assumed the line flux of [O III]λλ4960, 5008 to be separate from Hβ, giving us two unknown quantities, for which we need two equations. The reason for this separation is to account for cases where Hβ is in a filter, and [O III] is in another. In addition to that, we consider the continuum flux as another unknown quantity, requiring an additional equation. We then considered three equations of the same form, given by Eq. (A.1):

$\begin{matrix} W_{eff}^{A} F_{λ}^{A} = W_{eff}^{A} F_{λ, cont} + N_{H β}^{A} \times F_{H β} + N_{[O III]}^{A} \times F_{[O III]}, \end{matrix}$ $Mathematical equation: $$ \begin{aligned} W_{\rm eff}^{\text{A}}F_\lambda ^{\text{A}} = W_{\text{eff}}^{\text{A}}F_{\lambda , {\text{ cont}}} + N_{\mathrm{H}\beta }^\mathrm{A} \times F_{\mathrm{H}\beta } + N_{{[\mathrm{O}{\small { {\text{ III}}}} ]}}^A \times F_{{[\mathrm{O}{\small { {\text{ III}}}} ]}} ,\end{aligned} $$$ (A.1)

where A is the filter of interest, W_eff is the effective width of the filter in Å, F_λ is the spectral flux density in erg/s/cm²/Å, N_{Hβ;[O III]} is the coverage coefficient (which describes the proportion of the line covered by a filter compared to peak) and F_{Hβ;[O III]} is the monochromatic flux density in erg/s/cm².

With three filters, we can rewrite this system as Eq. (A.2) with A, B, and C being the three filters of interest. Because these equations mix monochromatic flux density and spectral flux densities, the vectors and the matrix are not uniform in units.

$\begin{matrix} (\begin{matrix} W_{eff}^{A} F_{λ}^{A} \\ W_{eff}^{B} F_{λ}^{B} \\ W_{eff}^{C} F_{λ}^{C} \end{matrix}) = (\begin{matrix} N_{H β}^{A} & N_{[O III]}^{A} & W_{eff}^{A} \\ N_{H β}^{B} & N_{[O III]}^{B} & W_{eff}^{B} \\ N_{H β}^{C} & N_{[O III]}^{C} & W_{eff}^{C} \end{matrix}) (\begin{matrix} F_{H β} \\ F_{[O III]} \\ F_{λ, cont} \end{matrix}) \end{matrix}$ $Mathematical equation: $$ \begin{aligned} \begin{pmatrix} W_{\rm eff}^\mathrm{A}F_\lambda ^\mathrm{A}\\ W_{\rm eff}^\mathrm{B}F_\lambda ^\mathrm{B}\\ W_{\rm eff}^\mathrm{C}F_\lambda ^\mathrm{C} \end{pmatrix} = \begin{pmatrix} N_{\mathrm{H}\beta }^{A}&N_{{[\mathrm{O}{\small { {\text{ III}}}} ]}}^{A}&W_{\rm eff}^\mathrm{A}\\ N_{\mathrm{H}\beta }^{B}&N_{{[\mathrm{O}{\small { {\text{ III}}}} ]}}^{B}&W_{\rm eff}^\mathrm{B}\\ N_{\mathrm{H}\beta }^{C}&N_{{[\mathrm{O}{\small { {\text{ III}}}} ]}}^{C}&W_{\rm eff}^\mathrm{C} \end{pmatrix} \begin{pmatrix} F_{\mathrm{H}\beta }\\ F_{{[\mathrm{O}{\small { {\text{ III}}}} ]}}\\ F_{\lambda , \text{ cont}} \end{pmatrix} \end{aligned} $$$ (A.2)

To compute the coverage coefficient N, which corresponds to the amount of line flux that we expect to retrieve in each filter for a given redshift, we used the convolution of the filter with a synthetic emission line. We built the synthetic emission lines separately for Hβ+[O III]λλ4960, 5008. In the former, we used a centred Gaussian of standard deviation 40Å, and for the latter two Gaussians separated by physical distances and a standard deviation of 40Å. We assumed a ratio of 2.98 between [O III]λ5007Å and [O III]λ4960Å (Storey & Zeippen 2000). Next, we convolved these synthetic emission lines with the JWST filters of interest and obtained a value of convolution depending on the wavelength. By normalising the convolution of the line at the redshift of interest by the maximum convolution, we obtain a coefficient N ∈ [0, 1], which we can measure for each line and filter. In the cases where the lines of interest lie outside medium filters or when medium filters are not detected or available, this method does not hold. When outside of the medium filters, an entire row of the coefficient matrix becomes zeros, leading to a singular matrix. In addition, when the medium filter is not detected, we lack one of the two equations. In these cases, we still try to estimate the [O III]+Hβ line flux by combining [O III]+Hβ, which reduces the number of equations needed.

In Fig. A.1 we show the comparison for the measurement of the equivalent width and line flux of [O III]+Hβ using SED fitting for GLIMPSE and JAGUAR data. The GLIMPSE data compares against the empirical measurements described above, and the JAGUAR data against the truth values obtained by the simulation. While we observe some scatter in equivalent width measurements, the line fluxes match well for both methods. All the measurements correlate, except for the empirically reconstructed equivalent width. For the equivalent width, the correlation coefficients are 0.79 for JAGUAR and 0.20 for empirical, and for line fluxes, 0.93 for JAGUAR and 0.67 for empirical. The comparison to GLIMPSE observations is more scattered than the comparison to the simulation. This might come from the lower amount of information used by the empirical method, which reduces the constraints on the measurement. In general, we observe an average good agreement between CIGALE and the two methods, validating our use of CIGALE for the final measurements.

Fig. A.1.

Comparison of measurement of the equivalent width and line flux of [O III]+Hβ using SED fitting (CIGALE) against the empirical method (blue dots) and the JAGUAR simulations (grey dots). Only sources with a valid empirical measurement were kept. For JAGUAR, we kept sources with redshift measurements within 1σ to ensure a fairer comparison with GLIMPSE. We did not display uncertainties for readability reasons. For the JAGUAR simulation, we added a linear fit on the equivalent widths and line fluxes (dashed grey lines).

Appendix B: Example of galaxies

In Fig. B.1 we display a few example SED fitting and galaxies cutouts. The three galaxies selected cover most situations: ID = 3644 is a relatively bright galaxy with limited magnification, ID = 6067 is a faint galaxy detected with limited magnification and finally ID = 56004 is very-faint but is strongly magnified by the lensing field.

Fig. B.1.

Galaxy stamp cutout and SED for three representative galaxies in the sample. Left: 3x3 panels show the cutouts for the nine GLIMPSE JWST filters. Right: the observed photometry for that galaxy (black), the photometry obtained using Bayesian inference of CIGALE (red) and the best fit SED (grey). The three galaxies were selected for being bright, faint without strong magnification, and very-faint with strong magnification.

Appendix C: Luminosity functions and associated data

In this appendix we report the additional LFs discussed in the paper that did not fit in it, as well as the data of all the LFs in this paper. All the tables in this section include the measurement of log₁₀Φ(L) for each luminosity bin L, for the two redshift bins. In addition, we report the average number of sources per bin, as well as the average completeness of these sources, including uncertainties.

C.1. [O III]+Hβ luminosity function

In Table C.1 we report the [O III]+Hβ LF data from Fig. 5. Only for this table, we report the average R3 used to transform our [O III]+Hβ LF to [O III]λ5008 and Hβ.

C.2. [O III]+Hβ luminosity function with flux cut to L ≥ 10⁴¹ erg s⁻¹

GLIMPSE probes [O III]+Hβ emitters fainter than any previous study, making it difficult to compare directly. Before GLIMPSE, Wold et al. (2025) was the deepest study of the [O III]+Hβ LF (reaching L_{[O III]+Hβ} ∼ 10⁴¹erg s⁻¹). As their faint-end slope is much steeper than GLIMPSE ( $α = - 2 . 07_{- 0.23}^{+ 0.22}$ $Mathematical equation: $ \alpha = -2.07_{-0.23}^{+0.22} $$ ), we truncate our [O III]+Hβ LF at 10⁴¹erg s⁻¹ to compare both studies and examine the behaviour of the fainter galaxies. Figure shows the [O III]+Hβ LF fixed at the limited depth. We observe a similar faint-end slope, which shows that the fainter galaxy population flattens the LF. The associated data are provided in Table C.2.

C.3. Separating [O III]λ5008 and Hβ from the [O III]+Hβ luminosity function

Differences in metallicity varies the R3 ratio (e.g. Maiolino & Mannucci 2019), which can explain differences between the [O III]+Hβ LF and the UV LF. By assuming an evolving R3 ratio following Eq. 12, we can separate the contribution of Hβ and [O III] in the [O III]+Hβ LF (see Sect. 4.4.2 for more detail and Table C.1 for the average R3 measurement in each luminosity bin). Figures and show the LF for Hβ and [O III]λ5008. Their associated data can be found in Tables C.3 and C.4.

Table C.1.

LF measured in Fig. 5.

Table C.2.

LF measured in Fig. C.1.

Fig. C.1.

LF of Hβ+[O III]λλ4960, 5008 for the different redshift ranges considered. For the fit of the LF, we limited the flux bins to the depth of (Wold et al. 2025), and discarded bins are displayed with a cross. We added previous studies of the LF using JWST/NIRCam GRISM instrument (Meyer et al. 2024; Matthee et al. 2023), a JWST/NIRCam wide field slitless spectroscopy study by (Sun et al. 2023), a JWST/NIRCam medium band survey by (Wold et al. 2025) and a former Spitzer study by De Barros et al. (2019). All the JWST surveys specifically study the [O III]λ5007Å, so we boosted their luminosity using their respective R3 ratio and [O III]λ5008/[O III]λ4960 = 2.98 (Storey & Zeippen 2000), according to our assumed line of ratios. This approximation matches the stacked median and flux measurements from these papers. The data can be found in Table C.3 and C.4 the parametrisation in Table 1.

Fig. C.2.

Hβ LF separated from the main [O III]+Hβ LF using the R3 ratio from Eq. 12. The Meyer et al. (2024) values are converted using their median R3 value. The data can be found in Table C.3 and the parametrisation can be found in Table 1.

Fig. C.3.

[O III]λ5008 LF separated from the main [O III]+Hβ LF using the R3 ratio from Eq. 12. The data can be found in Table C.4 and the parametrisation can be found in Table 1.

Fig. C.4.

Dust-corrected Hα LF obtained by converting [O III]+Hβ using the R3 ratio given by Eq. 12. The data from (Meyer et al. 2024) was converted using their median R3, assuming no dust attenuation as it was measured as negligible. Both used the usual Hα/Hβ = 2.86 (Osterbrock & Ferland 2006). We added the (Kennicutt & Evans 2012) conversion of Hα to SFR on top for a better understanding of the reached fluxes and implications. The data can be found in Table C.5 and the parametrisation in Table 1.

Table C.3.

Hβ luminosity function from Fig. C.2.

C.4. Dust correction to obtain the Hα luminosity function

To measure the ionising photon-production rate and the CSFRD, we need to correct the separated Hβ LF (Fig. C.2) for dust attenuation. Indeed, [O III] and Hβ are very close, and therefore, the dust attenuation is similar. However, Hα is affected differently from Hβ, which requires a correction. The detail of the correction can be found in Sect. 5.1. Figure shows the dust-corrected Hα LF with its associated data in Table C.5.

Table C.4.

[O III]λ5008 LF from Fig. C.3.

Table C.5.

Dust-corrected Hα LF from Fig. C.4.

Appendix D: Faint-end slope compilation

We report in Table D.1 the non-exhaustive list of measurement of the faint-end slope α for several nebular tracers, as well as UV. The measurement of GLIMPSE are already reported in Table 1.

Table D.1.

Non-exhaustive compilation of α measurements for [O III], Hβ, Hα and UV LFs.

Appendix E: Cosmic star formation rate density

We report in Table E.1 the non-exhaustive list of measurement of SFRD from the literature, as well as the associated measurement from GLIMPSE. We limited the compilation to a few papers with lower integration limit nearby 0.3M_⊙ yr⁻¹, as the average evolution of SFRD is well known. Every measurement was converted to the Chabrier (2003) IMF. The data is shown in Fig. 13.

Table E.1.

Non-exhaustive list of SFRD measurements shown in Fig. 13 from UV, Hα and [O III]+Hβ between z ∼ 0 − 10.

All Tables

Table 1.

Schechter parametrisation of the LFs.

	Fig. 2. Top panel: Evolution of the effective volume probed for each UV magnitude M_UV for our two redshift bins 7 < z < 8 (blue) and 8 < z < 9 (orange) with 1σ uncertainties in the shaded area. Bottom panel: Completeness estimate for the same M_UV and redshift bins.
In the text

	Fig. 3. Histogram of the dust attenuation parameter E(B − V)_l sampled from the SED fitting of our sample of 7 < z < 9 galaxies. N = 10 000 realisations of the measurement are performed for each galaxy, following the distributions of the SED fitting, and the histogram is normalised by N.
In the text

	Fig. 4. Correlation between the [O III]+Hβ line flux against M_UV for our GLIMPSE sample (circles) and measurements from Meyer et al. (2024, shown as crosses). To match our observations, we only kept sources within the same 7 < z < 9 redshift range. Bayesian error estimates using CIGALE are also shown for the GLIMPSE measurements.
In the text

	Fig. 6. Redshift evolution of the fitting parameters from Fig. 5.
In the text

	Fig. 7. Distribution of the [O III]+Hβ equivalent width for each redshift bin. This figure is not completeness-corrected. The colours represent the two redshift ranges and the lines are KDE estimate of the distribution for clarity.
In the text

	Fig. 8. Redshift evolution of the faint-end slope (α) in the Schechter (diamond) and double power-law (circle) fits for the UV continuum (grey symbols) and nebular emission line (coloured points) LFs. The orange squares with a black outline correspond to this study. The compilation of the data, including GLIMPSE, can be found in Table D.1.
In the text

	Fig. C.2. Hβ LF separated from the main [O III]+Hβ LF using the R3 ratio from Eq. 12. The Meyer et al. (2024) values are converted using their median R3 value. The data can be found in Table C.3 and the parametrisation can be found in Table 1.
In the text

	Fig. C.3. [O III]λ5008 LF separated from the main [O III]+Hβ LF using the R3 ratio from Eq. 12. The data can be found in Table C.4 and the parametrisation can be found in Table 1.
In the text

A GLIMPSE into the very faint end of the Hβ+[O III]λλ4960,5008 luminosity function at z ∼ 7 – 9 behind Abell S1063

1. Introduction

2. Data

2.1. Observations

2.2. Simulations

3. Method and results

3.1. Source selection

3.2. Gravitational lensing

3.3. Estimating the effective survey volume

3.4. Line flux measurement of [O III]+Hβ

3.5. Dust attenuation

3.6. [O III]+Hβ completeness estimation

4. The [O III]+Hβ luminosity function at z ∼ 7 − 9

4.1. The [O III]+Hβ luminosity function at z ∼ 7 − 9

4.2. The characteristic luminosity and normalisation parameters

4.3. Faint end of the [O III]+Hβ luminosity function

4.4. The faint-end slope of the nebular emission lines and UV luminosity functions

4.4.1. The impact of bursty star formation histories on the [O III]+Hβ LF

4.4.2. The impact of metallicity on the [O III]+Hβ LF

4.4.3. The early signs of an UV LF turnover

4.5. Possible missing sources with our method

5. Implications and discussion

5.1. The ionising photon budget at z ∼ 7 − 9

5.2. The cosmic star formation rate density

6. Conclusions

Acknowledgments

References

Appendix A: Validating the SED fitting flux and EW measurement using flux excess

Appendix B: Example of galaxies

Appendix C: Luminosity functions and associated data

C.1. [O III]+Hβ luminosity function

C.2. [O III]+Hβ luminosity function with flux cut to L ≥ 1041 erg s−1

C.3. Separating [O III]λ5008 and Hβ from the [O III]+Hβ luminosity function

C.4. Dust correction to obtain the Hα luminosity function

Appendix D: Faint-end slope compilation

Appendix E: Cosmic star formation rate density

All Tables

All Figures

A GLIMPSE into the very faint end of the Hβ+[O III]λλ4960,5008 luminosity function at z ∼ 7 – 9 behind Abell S1063

3.4. Line flux measurement of [O III]+Hβ

3.6. [O III]+Hβ completeness estimation

4. The [O III]+Hβ luminosity function at z ∼ 7 − 9

4.1. The [O III]+Hβ luminosity function at z ∼ 7 − 9

4.3. Faint end of the [O III]+Hβ luminosity function

4.4.1. The impact of bursty star formation histories on the [O III]+Hβ LF

4.4.2. The impact of metallicity on the [O III]+Hβ LF

C.1. [O III]+Hβ luminosity function

C.2. [O III]+Hβ luminosity function with flux cut to L ≥ 10⁴¹ erg s⁻¹

C.3. Separating [O III]λ5008 and Hβ from the [O III]+Hβ luminosity function