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

HAL
open science

# The star-forming and ionizing properties of dwarf z 6-9 galaxies in JADES: insights on bursty star formation and ionized bubble growth

Ryan Endsley, Daniel P. Stark, Lily Whitler, Michael W. Topping, Benjamin

D. Johnson, Brant Robertson, Sandro Tacchella, Stacey Alberts, William M.

Baker, Rachana Bhatawdekar, et al.

► To cite this version:

Ryan Endsley, Daniel P. Stark, Lily Whitler, Michael W. Topping, Benjamin D. Johnson, et al. The star-forming and ionizing properties of dwarf z 6-9 galaxies in JADES: insights on bursty star formation and ionized bubble growth. Monthly Notices of the Royal Astronomical Society, 2024, 533, pp.1111-1142. 10.1093/mnras/stae1857. insu-04851527

HAL Id: insu-04851527

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

Submitted on 20 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The star-forming and ionizing properties of dwarf $z \sim 6-9$ galaxies in JADES: insights on bursty star formation and ionized bubble growth

Ryan Endsley<sup>1,★</sup> Daniel P. Stark,<sup>2</sup> Lily Whitler<sup>1,2</sup> Michael W. Topping<sup>1,2</sup> Benjamin D. Johnson,<sup>3</sup>

Brant Robertson,<sup>4</sup> Sandro Tacchella<sup>1,5,6</sup> Stacey Alberts<sup>1,2</sup> William M. Baker<sup>1,5,6</sup>

Rachana Bhatawdekar<sup>1,7</sup> Kristan Boyett<sup>1,8,9</sup> Andrew J. Bunker,<sup>10</sup> Alex J. Cameron<sup>1,10</sup>

Stefano Carniani<sup>1,11</sup> Stephane Charlot<sup>1,12</sup> Zuyi Chen<sup>1,2</sup> Jacopo Chevallard<sup>1,10</sup> Emma Curtis-Lake<sup>1,13</sup>

A. Lola Danhaive,<sup>5,6</sup> Eiichi Egami,<sup>2</sup> Daniel J. Eisenstein<sup>1,3</sup> Kevin Hainline,<sup>2</sup> Jakob M. Helton,<sup>2</sup>

Zhiyuan Ji,<sup>2</sup> Tobias J. Looser,<sup>5,6</sup> Roberto Maiolino<sup>1,5,6,14</sup> Erica Nelson,<sup>15</sup> Dávid Puskás,<sup>5,6</sup>

George Rieke,<sup>2</sup> Marcia Rieke,<sup>2</sup> Hans-Walter Rix,<sup>16</sup> Lester Sandles,<sup>5,6</sup> Aayush Saxena<sup>1,10,17</sup>

Charlotte Simmonds<sup>1,5,6</sup> Renske Smit,<sup>18</sup> Fengwu Sun,<sup>2</sup> Christina C. Williams,<sup>19</sup>

Christopher N. A. Willmer,<sup>2</sup> Chris Willott<sup>20</sup> and Joris Witstok<sup>1,5,6</sup>

*Affiliations are listed at the end of the paper*

Accepted 2024 July 29. Received 2024 July 26; in original form 2023 June 21

## ABSTRACT

Reionization is thought to be driven by faint star-forming galaxies, but characterizing this population has long remained very challenging. Here, we utilize deep nine-band JADES (JWST Advanced Deep Extragalactic Survey)/NIRCam (Near-Infrared Camera) imaging to study the star-forming and ionizing properties of 756 $z \sim 6-9$ galaxies, including hundreds of very ultraviolet (UV)-faint objects ($M_{UV} > -18$). The faintest ($m \sim 30$) galaxies in our sample typically have stellar masses of $M_* \sim (1-3) \times 10^7 M_\odot$ and young light-weighted ages ($\sim 50$ Myr), though some show strong Balmer breaks implying much older ages ($\sim 500$ Myr). We find no evidence for extremely massive galaxies ($> 3 \times 10^{10} M_\odot$) in our sample. We infer a strong (factor $>2$) decline in the typical [O III]+H$\beta$ equivalent widths (EWs) towards very faint $z \sim 6-9$ galaxies, yet a weak UV luminosity dependence on the H$\alpha$ EWs at $z \sim 6$. We demonstrate that these EW trends can be explained if fainter galaxies have systematically lower metallicities as well as more recently declining star formation histories relative to the most UV-luminous galaxies. Our data provide evidence that the brightest galaxies are frequently experiencing a recent strong upturn in star formation rate. We also discuss how the EW trends may be influenced by a strong correlation between $M_{UV}$ and Lyman continuum escape fraction. This alternative explanation has dramatically different implications for the contribution of galaxies along the luminosity function to cosmic reionization. Finally, we quantify the photometric overdensities around two $z > 7$ strong Ly$\alpha$ emitters. One Ly$\alpha$ emitter lies close to a strong photometric overdensity, while the other shows no significant nearby overdensity, perhaps implying that not all strong $z > 7$ Ly$\alpha$ emitters reside in large ionized bubbles.

**Key words:** galaxies: evolution – galaxies: high-redshift – dark ages, reionization, first stars.

## 1 INTRODUCTION

The formation and assembly of galaxies within the first billion years of cosmic history directly influenced the large-scale ionization state of the Universe via the process of hydrogen reionization (e.g. Stark 2016; Dayal & Ferrara 2018; Robertson 2022). Recent Ly$\alpha$ forest measurements from a statistical sample of high-redshift quasars indicate that essentially all hydrogen atoms in the intergalactic medium (IGM) had been reionized by $z = 5.3$ (Bosman et al. 2022), approximately 1.1 Gyr after the big bang. Additional quasar and galaxy spectra at $z > 6$, along with measurements of the cosmic

microwave background, imply that reionization was about halfway complete $\sim 400$ Myr earlier at $z \sim 7-8$ (e.g. Davies et al. 2018; Mason et al. 2018; Planck Collaboration VI 2020; Yang et al. 2020; Goto et al. 2021). The primary agents of hydrogen reionization appear most likely to be star-forming galaxies given that constraints on the $z \gtrsim 6$ quasar luminosity function imply that active supermassive black holes were very rare at early times (e.g. Shen et al. 2020; Jiang et al. 2022; Matsuoka et al. 2023).

Prior to JWST, over 1000 Lyman-break galaxies at $z > 6$ had been identified from deep Hubble Space Telescope (HST) imaging (e.g. McLure et al. 2013; Atek et al. 2015; Bouwens et al. 2015, 2022b; Finkelstein et al. 2015, 2022; McLeod, McLure & Dunlop 2016; Livermore, Finkelstein & Lotz 2017; Ishigaki et al. 2018; Oesch et al. 2018; Salmon et al. 2020) as well as wide-area ground-based

<footnote>* E-mail: ryan.endsley@austin.utexas.edu</footnote>
imaging (e.g. Bowler et al. 2014, 2020; Stefanon et al. 2017, 2019; Ono et al. 2018; Endsley et al. 2021a; Harikane et al. 2022; Donnan et al. 2023). These studies demonstrated that the $z \sim 6-10$ galaxy ultraviolet (UV) continuum luminosity functions had very steep faint-end slopes ($\alpha \sim -2$; e.g. Atek et al. 2015; Livermore et al. 2017; Ishigaki et al. 2018; Bouwens et al. 2022b), clearly indicating that very UV-faint sources ($M_{UV} > -18$) greatly dominated the galaxy number counts during reionization.

The very high relative abundance of the faintest ($M_{UV} > -18$) $z \gtrsim 6$ galaxies has long motivated efforts to characterize their physical properties. Deep HST observations revealed that these systems often exhibit very blue rest-UV continuum slopes ($-2.5 \lesssim \beta \lesssim -2$ where $f_{\lambda} \propto \lambda^{\beta}$; e.g. Dunlop et al. 2012; Finkelstein et al. 2012; Rogers, McLure & Dunlop 2013; Bouwens et al. 2014; Bhatawdekar & Conselice 2021) as well as very compact rest-UV morphologies (effective radii $r_e \lesssim 200$ pc; Shibuya, Ouchi & Harikane 2015; Neufeld et al. 2022; Bouwens et al. 2022a). Based on results from local ($z \sim 0.3$) samples, such properties hint that very UV-faint reionization-era galaxies might be efficient at leaking ionizing photons into the IGM (Chisholm et al. 2022; Flury et al. 2022). However, detailed statistical studies of the star-forming and ionizing properties of this faint, abundant $z \gtrsim 6$ population have been precluded by the lack of deep data probing the rest-frame optical portion of their spectral energy distributions (SEDs).

For much of the past 20 yr, the only instrument capable of delivering any constraints on the rest-optical SEDs of high-redshift galaxies was the Infrared Array Camera (IRAC) on board the *Spitzer Space Telescope*. Even so, IRAC could only reach a $5\sigma$ imaging sensitivity of $m \sim 26.5$ ($M_{UV} \lesssim -20.5$ at $z \sim 7$) in the deepest pointings (120 h; Labbé et al. 2013, 2015; Oesch et al. 2013; Stefanon et al. 2021a). While dedicated observations in lensing fields were able to push this sensitivity to intrinsically fainter UV luminosities (Smit et al. 2014; Huang et al. 2016; Strait et al. 2020), statistical constraints on the rest-optical SEDs of the faintest reionization-era galaxies remained very limited. Consequently, in the lead-up to *JWST*, much attention was concentrated on understanding the physical properties of relatively bright ($M_{UV} \lesssim -20$) Lyman-break $z \gtrsim 6$ galaxies (e.g. Laporte et al. 2021; Endsley et al. 2021a; Stefanon et al. 2021b; Castellano et al. 2022; Schouws et al. 2022; Sommovigo et al. 2022; Tacchella et al. 2022; Witstok et al. 2022; Topping et al. 2022a; Whitler et al. 2023b), though even these analyses suffered from IRAC's poor angular resolution (FWHM~2 arcsec) and SED sampling (only two imaging filters at 3–5 µm, both with broad bandpasses).

Within the past year alone, data from *JWST* have greatly advanced our understanding of reionization-era galaxies. Rest-optical spectra from the Near-Infrared Spectrograph (NIRSpec, Jakobsen et al. 2022; Böker et al. 2023) have generally revealed very strong nebular line emission (e.g. [O III], $H\alpha$, and $H\beta$) coupled with signatures of hard radiation fields ([O III]/[O II] > 10) and low gas-phase metallicities ($\lesssim 0.1 Z_{\odot}$) among $z > 6$ galaxies (e.g. Arrabal Haro et al. 2023; Bunker et al. 2023; Cameron et al. 2023; Fujimoto et al. 2023; Larson et al. 2023; Sanders et al. 2023; Saxena et al. 2023; Tang et al. 2023; Tacchella et al. 2023b; Curti et al. 2024). This is consistent with expectations from findings prior to *JWST* (e.g. Labbé et al. 2013; Smit et al. 2014, 2015; Stark et al. 2015a, b, 2017; Roberts-Borsani et al. 2016; Laporte et al. 2017; Mainali et al. 2017, 2018; Schmidt et al. 2017; De Barros et al. 2019; Hutchison et al. 2019; Tang et al. 2019, 2021; Endsley et al. 2021a, b; Stefanon et al. 2022b). Less expected subpopulations of $z > 5$ galaxies are also emerging from early NIRSpec observations, including sources with (sometimes tentative) contributions from active galactic nuclei (AGNs, Harikane

et al. 2023; Kocevski et al. 2023; Larson et al. 2023; Onoue et al. 2023; Übler et al. 2023; Maiolino et al. 2024) as well as objects that are in a very inactive stage of star formation (Strait et al. 2023; Looser et al. 2024).

The Near-Infrared Camera (NIRCam; Rieke et al. 2005, 2023) on-board *JWST* is also proving to be an invaluable tool for studying the demographics of reionization-era galaxies. This is not only because of NIRCam's dramatic improvement in sensitivity, angular resolution, and SED sampling over IRAC, but also because Lyman-break $z \gtrsim 6$ selections with imaging are very efficient at yielding highly complete samples of objects with faint continua. Much of the early NIRCam studies on $z \gtrsim 6$ galaxies focused on data from the Cosmic Evolution Early Release Science (CEERS; Finkelstein et al. 2023) and GLASS (Treu et al. 2022) surveys which immediately pushed down to extremely deep *Hubble* imaging depths ($m_{AB} \sim 29$ mag at $5\sigma$) across 1–5 µm with just ~3 h of observations per photometric band. These data (among that from other early *JWST*/NIRCam surveys) have delivered a wealth of insight into faint ($M_{UV} \lesssim -19$) $z \gtrsim 6$ galaxies including their ages, stellar masses, nebular line strengths, UV slopes, dust attenuation strengths, and morphologies (e.g. Topping et al. 2022b; Barrufet et al. 2023; Chen et al. 2023; Endsley et al. 2023; Furtak et al. 2023b; Hsiao et al. 2023; Labbé et al. 2023; Laporte et al. 2023; Leethochawalit et al. 2023; Rodighiero et al. 2023; Treu et al. 2023; Dressler et al. 2023; Whitler et al. 2023a). While these studies have clearly advanced our understanding of reionization-era systems far beyond that possible prior to *JWST*, a detailed statistical analysis of the properties of very faint ($M_{UV} > -18$) $z \gtrsim 6$ galaxies has yet to be undertaken. Such an endeavour is clearly warranted as these very faint objects are often though to contribute substantially to cosmic reionization, and also are the likely progenitors of more typical galaxies at lower redshifts.

Here, we take steps to address this shortcoming by utilizing NIRCam imaging taken as part of the *JWST* Advanced Deep Extragalactic Survey (JADES; Eisenstein et al. 2023). JADES is a collaborative effort of the NIRCam, NIRSpec, and U.S. Mid-Infrared Instrument (MIRI; Wright et al. 2023) teams utilizing coordinated parallels to maximize science outcomes across the ≈770 h of allocated observing time. All data in JADES are being taken over the Great Observatories Origins Deep Survey (GOODS) fields in the Northern (GOODS-N) and Southern (GOODS-S) hemispheres (Giavalisco et al. 2004) which contain some of the deepest *Hubble* imaging ever obtained (e.g. Ellis et al. 2013; Illingworth et al. 2013). By adding exceptionally deep ($m_{AB} \sim 30-31$ mag $5\sigma$) NIRCam imaging in several bands from ~1–5 µm (often including two medium-band filters) across >200 arcmin² of the GOODS fields, JADES is opening a completely new window on the high-redshift Universe (e.g. Baker et al. 2023; Curtis-Lake et al. 2023; Robertson et al. 2023; Tacchella et al. 2023a; Dressler et al. 2024; Hainline et al. 2024; Helton et al. 2024).

The primary goal of this paper is to utilize the deep NIRCam imaging from JADES to measure the rest-UV + optical SEDs among a large sample of very faint ($M_{UV} > -18$) $z \sim 6-9$ galaxies and statistically characterize their physical properties in detail, comparing with more luminous systems. One of the key conclusions of this paper relates to the star formation histories (SFHs) of reionization-era galaxies, complementing the JADES/NIRCam investigation in Dressler et al. (2024) which utilized broad-band SEDs to infer the presence of multiple bursts by considering non-parametric SFHs. Here, we utilize the two long-wavelength NIRCam medium bands from JADES to focus on insight from nebular line emission about variations in star formation activity on short (~3 Myr) time-scales.
The structure of this paper is as follows. We begin by describing the imaging data, source extraction, photometric calculations, sample selection, and photoionization modelling in Section 2. Next, we present and discuss the UV luminosities, photometric redshifts, stellar masses, and light-weighted ages among our sample (Section 3). We then proceed to utilize photometric constraints on strong rest-optical nebular lines ([O III] + H β and H α; Section 4) to statistically analyse the star-forming and ionizing properties of relatively bright and very faint reionization-era galaxies (Section 5). Finally, we quantify the photometric overdensities surrounding two strong Ly α emitters (LAEs) at $z > 7$ within the JADES footprints to improve our understanding of the connection between strong Ly α and ionized bubbles deep in the epoch of reionization (Section 6). Our main conclusions are summarized in Section 7.

Throughout this paper, we quote all magnitudes in the AB system, report EWs in the rest frame, assume a Chabrier (2003) initial mass function (IMF) with limits of 0.1–300 $M_\odot$, and adopt a flat lambda-cold dark matter cosmology with parameters $h = 0.7$, $\Omega_M = 0.3$, and $\Omega_\Lambda = 0.7$.

## 2 IMAGING DATA AND SAMPLE SELECTION

In this section, we first describe the JADES/NIRCam imaging over the GOODS-S and GOODS-N fields as well as the complementary Hubble Advanced Camera for Surveys (ACS) imaging that assists with our Lyman-break dropout selection (Section 2.1). We next detail our procedure for source extraction and photometric measurements in Section 2.2, then describe our Lyman-break selections of $z \sim 6$ and $\sim 7-9$ galaxies (Section 2.3). Lastly, we describe the photoionization SED models used throughout this work to infer the physical properties of the galaxies within our sample (Section 2.4).

### 2.1 JWST and HST Imaging

We consider all JADES/NIRCam imaging data taken prior to 2023 February 10 (see Rieke & the JADES Collaboration 2023). These data include imaging in the broad-band F090W, F115W, F150W, F200W, F277W, F356W, and F444W filters,<fnref n="1" /> as well as the F410M medium-band filter across the full footprints in both the GOODS-S and GOODS-N fields. Additional medium-band F335M imaging is available over the majority of the imaging area taken up to 2023 February 10. Because the medium-band filters are particularly useful for determining whether galaxies show strong nebular line emission or continuum discontinuities (e.g. Balmer breaks), we here restrict our analysis to the ≈90 arcmin<fnref n="2" /> of JADES/NIRCam imaging that include data in all nine of the above filters. This imaging is split approximately evenly into the GOODS-S (≈44 arcmin<fnref n="2" />) and GOODS-N fields (≈46 arcmin<fnref n="2" />).

The NIRCam imaging reduction algorithm applied here largely follows that described in previous works (e.g. Bezanson et al. 2022; Merlin et al. 2022; Donnan et al. 2023; Endsley et al. 2023; Rieke & the JADES Collaboration 2023; Robertson et al. 2023; Tacchella et al. 2023a; Williams et al. 2023; Bagley et al. 2023) and is detailed below. We begin by processing the uncalibrated (\*rate\_fits) NIRCam exposures through the stage 1 step of the JWST Science Calibration Pipeline<fnref n="2" /> (v1.9.4) using the JWST Calibration Reference Data System (CRDS) context map jwst.1045.pmap. During this step, we implement a custom snowball (Rigby et al. 2023) masking algorithm largely based on the methods described in Bagley et al. (2023). From the output \*uncal\_fits files, we subtract prominent artefacts (e.g. wisps; Rigby et al. 2023) from the short-wavelength NIRCam imaging data using custom-built templates. These custom templates were created by combining all F090W, F115W, F150W, and F200W imaging data (as of 2023 February 1) from JADES as well as the deep, blank-field public programs CEERS (Finkelstein et al. 2023), Public Release IMaging for Extragalactic Research (PRIMER; PI: J. Dunlop), and Next Generation Deep Extragalactic Exploratory Public (NGDEEP) survey (Bagley et al. 2024). All snowball-masked and artefact-subtracted \*uncal\_fits files are visually inspected; for a very small fraction of imaging exposures, we implement additional masking of remaining prominent artefacts or remove the exposures from the reduction entirely if the overall data quality appears poor. From the resulting files, we build custom sky flats in all JADES imaging bands considered here, folding in additional data from CEERS, PRIMER, NGDEEP, and First Reionization Epoch Spectroscopically Complete Observations (FRESCO; Oesch et al. 2023).

The \*uncal\_fits files are then processed through the stage 2 step of the JWST Science Calibration Pipeline with our custom sky flats and the photometric zero-points from Boyer et al. (2022) that remain the most recent values implemented in CRDS. From the resulting calibrated \*_calember\_fits, we then subtract-off 1/f noise (Schlawin et al. 2020) using the sigma-clipped median values along given rows (on an amplifier-by-amplifier basis) and columns. Next, we subtract-off the 2D background using the SEP package (Barbary 2016) and follow the methods described in Bagley et al. (2023). The final \*_calember\_fits files are then processed through the stage 3 step of the Calibration Pipeline to create mosaics with a pixel scale of 30 mas pixel⁻<fnref n="1" />. Using the TWEAKREG software<fnref n="3" />, these mosaics are astrometrically matched to the HST imaging over the GOODS fields that have been registered to the Gaia frame (see below). With the primary intention of improving the astrometric alignment, we fold in F444W data from FRESCO (which significantly extends the covered area in each GOODS field) into our reduction and first align the F444W mosaics to the HST Wide Field Camera 3 (WFC3) F160W images, and then subsequently align all other NIRCam mosaics to that of F444W in the respective field.

For our dropout selection of $z \sim 6-9$ galaxies, we also utilize imaging from HST taken over the GOODS fields. The HST mosaics adopted here come from the Hubble Legacy Field archive (HLF; see Illingworth et al. 2013; Whitaker et al. 2019, and references therein) and are registered to the Gaia frame using the astrometry from the CHARGE images (Brammer, private communication) as described in Williams et al. (2023). All HLF mosaics used here (v2.0) have the same pixel scale as the final NIRCam mosaics (30 mas pixel⁻<fnref n="1" />). Here, we focus on using the optical ACS imaging for our dropout selections, which includes data in the F435W, F606W, F775W, F814W, and F850LP bands.

To obtain reliable photometric colours across the 0.4–5 µm SEDs, we must account for the fact that the angular resolution of the ACS and NIRCam bands considered here differ by up to a factor of ≈4. We therefore convolve all mosaics to the point spread function (PSF) of F444W which has the poorest angular resolution among these bands (FWHM ≈0.15 arcsec). We refer the interested reader to Endsley et al. (2023) for details on our procedure for constructing the

<footnote n="1">¹F070W imaging is also available over a single pointing as of 2023 February 10. Because this represents a small fraction of the total JADES/NIRCam imaging at the time, we here ignore the F070W data.</footnote>

<footnote n="2">²https://jwst-pipeline.readthedocs.io/en/latest/index.html</footnote>

<footnote n="3">³https://drizzlepac.readthedocs.io/en/latest/tweakreg.html</footnote>
empirical PSFs (one per band and field) and creating the convolution kernels using PHOTUTILS (Bradley et al. 2022).

## 2.2 Source extraction and photometry

To identify objects across the JADES footprints, we run SOURCE EXTRACTOR (Bertin & Arnouts 1996) on an inverse-variance weighted stack of the PSF-convolved F200W, F277W, F335M, F356W, F410M, and F444W mosaics.<fnref n="4" /> The photometry is computed in Kron (1980) apertures following commonly used procedures for high-redshift galaxies (e.g. Oesch et al. 2013; Bouwens et al. 2015; Finkelstein et al. 2015; Endsley et al. 2021a, 2023). First, we measure the photometry on the PSF-matched mosaics of each band within elliptical apertures with a Kron (1980) factor of $k = 1.2$, as these apertures have been shown to yield optimal signal-to-noise ratio (S/N, Finkelstein et al. 2022). The photometric errors are determined separately for each object and each band to take into account both the size and shape of the aperture, as well as the local background level from the varying exposure times across the mosaics. Specifically, the photometric errors are measured as the standard deviation of flux values computed within randomly placed elliptical apertures of the size/shape of interest) in nearby empty regions of the mosaic (as determined using SEP) with similar background pixel flux variance as the source of interest.

The photometry in the $k = 1.2$ apertures (and the associated error values) are corrected to total flux in two stages. The fluxes and their associated errors are first multiplied by the ratio of the flux measured in $k = 2.5$ apertures divided by that of the flux in the $k = 1.2$ apertures, with these measurements performed on the inverse-variance weighted stack used for the SOURCE EXTRACTOR detection. The second stage is correcting for flux outside the $k = 2.5$ apertures by assuming a point-source light profile outside these apertures and adopting our empirical F444W PSFs. We have verified that our measured photometry for $z \sim 6-9$ galaxies located in the current public release region of JADES (see Rieke & the JADES Collaboration 2023) is broadly consistent with the photometry released in that associated catalogue.

Because of NIRCam's sensitivity, a significant number of our identified Lyman-break dropout $z \sim 6-9$ galaxies have $k = 2.5$ apertures that contain one or more different objects, artificially boosting their recovered total fluxes. We therefore utilize the neighbour subtraction algorithm described in Endsley et al. (2023) to subtract-off any neighbouring object within the $k = 2.5$ aperture (as determined by the SOURCE EXTRACTOR segmentation map) and recompute the photometry on those images. Moreover, because the morphologies of $z \gtrsim 6$ galaxies are sometimes clumpy (e.g. Chen et al. 2023; Treu et al. 2023; Hainline et al. 2024), the galaxies we intend to study can be identified as multiple separate objects by SOURCE EXTRACTOR. If two or more nearby objects satisfy one of our Lyman-break dropout selections ($z \sim 6$ or $\sim 7-9$; see below) and have overlapping $k = 2.5$ apertures, we determine the smallest elliptical aperture that contains all pixels from their combined segmentation maps and recompute the photometry in that aperture, subtracting off other nearby objects when appropriate.

In Table 1, we summarize the range of $5\sigma$ depths in each HST and NIRCam band among our final Lyman-break galaxy sample

<table><caption>**Table 1.** Summary of the $5\sigma$ depths achieved in each band for our final Lyman-break sample of F775W and F090W dropouts. Because the depth for a given source depends not only on its position within the inhomogeneous GOODS HST and NIRCam imaging, but also the source size, we report the 10th, 50th, and 90th percentile $5\sigma$ depths in each band to show the range relevant to our sample. All photometry is computed on mosaics convolved to the F444W PSF (FWHM = 0.15 arcsec).</caption><thead><tr><th rowspan="2">Band</th><th colspan="3">5σ depths of z ~ 6 − 9 sample [AB mag]</th></tr><tr><th>10th percentile</th><th>50th percentile</th><th>90th percentile</th></tr></thead><tbody><tr><td>F435W</td><td>28.0</td><td>28.6</td><td>30.0</td></tr><tr><td>F606W</td><td>28.1</td><td>28.7</td><td>30.3</td></tr><tr><td>F775W</td><td>27.6</td><td>28.4</td><td>30.2</td></tr><tr><td>F814W</td><td>28.0</td><td>28.5</td><td>29.0</td></tr><tr><td>F850LP</td><td>27.2</td><td>27.8</td><td>29.3</td></tr><tr><td>F090W</td><td>28.1</td><td>29.2</td><td>29.9</td></tr><tr><td>F115W</td><td>28.5</td><td>29.5</td><td>30.2</td></tr><tr><td>F150W</td><td>28.4</td><td>29.4</td><td>30.1</td></tr><tr><td>F200W</td><td>28.5</td><td>29.5</td><td>30.1</td></tr><tr><td>F277W</td><td>28.7</td><td>29.7</td><td>30.4</td></tr><tr><td>F335M</td><td>28.1</td><td>29.1</td><td>29.6</td></tr><tr><td>F356W</td><td>28.8</td><td>29.7</td><td>30.2</td></tr><tr><td>F410M</td><td>28.0</td><td>29.0</td><td>29.7</td></tr><tr><td>F444W</td><td>28.2</td><td>28.9</td><td>29.5</td></tr></tbody></table>

of F775W and F090W dropouts (see Section 2.3). Given the inhomogeneous exposure maps of the GOODS imaging as well as our Kron aperture photometry, the achieved depth for a given galaxy depends on both its on-sky location as well as its angular size. We therefore report the 10th, 50th, and 90th percentile $5\sigma$ depths in each band to show the range relevant to our sample, noting that this is the depth achieved using the PSF-matched mosaics. The median $5\sigma$ depth of our final $z \sim 6-9$ galaxy sample is $m_{AB} \approx 28.5$ in F435W, F606W, F775W, and F814W, though $\approx 0.7$ mag shallower in F850LP ($m_{AB} = 27.8$). In NIRCam, the typical $5\sigma$ depths range from $m_{AB} \approx 29.0$ in F410M and F444W to $m_{AB} \approx 29.7$ in F277W and F356W. The deepest imaging in GOODS-S reaches $5\sigma$ depths of $m_{AB} \approx 30.2$ in several ACS and NIRCam bands for compact objects.

## 2.3 Selection of Lyman-break $z \sim 6-9$ galaxies

Our primary goal is to statistically characterize the physical properties of very UV-faint reionization-era galaxies for the first time, and moreover compare their properties to that of the more UV-luminous population. We therefore restrict our analysis to galaxies at $z \sim 6-9$ where the lower bound reflects when reionization remains significantly incomplete (e.g. Nasir & D’Aloisio 2020; Bosman et al. 2022; Zhu et al. 2022) while the upper bound of $z \sim 9$ is chosen such that the NIRCam SEDs remain sensitive to the [O III] + H$\beta$ emission line set, a key photometric probe of high-redshift galaxy properties (e.g. Roberts-Borsani et al. 2016; Stark et al. 2017; Tang et al. 2019; Endsley et al. 2021a).

Following the approach of many previous studies (e.g. Bunker et al. 2004; Bouwens et al. 2015; Ono et al. 2018; Stefanon et al. 2019; Endsley et al. 2023; Leethochawalit et al. 2023), we identify $z \sim 6-9$ galaxies via a Lyman-break colour selection. These Lyman-break selections rely on the appearance of a sharp spectral discontinuity at $\lambda_{rest} = 1216$ Å that is imposed by very strong absorption blueward of the Ly$\alpha$ line from H$\alpha$ in the intervening IGM (e.g. Inoue et al. 2014). Other studies have also selected high-redshift galaxies via photometric redshifts, including the study of Hainline et al. (2024) which aimed at establishing a census of $z > 8$ galaxies within the JADES/NIRCam data set and providing a first

<footnote n="4">⁴We do not expect that adopting such a red-sided detection image is strongly biasing our selection against very blue high-redshift galaxies. Topping et al. (2024) adopted the same detection image and none the less found a large sample of high-redshift galaxies with robust very blue UV slopes ($\beta < -2.8$).</footnote>
discussion of their colours and morphologies. Such photometric redshift selections have the advantage of folding in all photometric data points in determining the probability that a given object lies in the redshift interval of interest. However, here we explicitly choose not to select on photometric redshifts as doing so may bias our sample preferentially towards objects with strong rest-optical lines (which imprint unique long-wavelength NIRCam colour patterns) and the goal in this work is to characterize reionization-era galaxy physical properties. In Appendix A, we compute the photometric redshifts for all objects in our sample excluding the rest-optical photometry during the fits. The large majority of sources have a high probability (>90 per cent) of lying at high redshift ($z > 4$; see Fig. A1), and we have verified that accounting for photometric redshifts has no significant impact on our main results regarding the statistical ionizing and star-forming properties of our sample (see Appendix A).

Because the $Ly\alpha$ break spans a considerable range in observed wavelength across our targeted $z \sim 6-9$ interval, we divide our selection according to two separate sets of criteria, one for $z \sim 6$ and one for $z \sim 7-9$, which we discuss in sequence below. At $z \sim 6$, the $Ly\alpha$ break is located at $\approx 0.85 \ \mu m$ and thus galaxies at this redshift should appear as strong ACS/F775W dropouts (e.g. Stanway, Bunker & McMahon 2003; Bouwens et al. 2004; Bunker et al. 2004). We therefore begin our $z \sim 6$ galaxy selection with the following colour cuts:

(i) $F775W - F090W > 1.2$;

(ii) $F090W - F150W < 1.0$;

(iii) $F775W - F090W > F090W - F150W + 1.2$.

These criteria enforce the presence of a strong break between the F775W and F090W bands (at least 1.2 mag) as well as a much flatter colour between F090W and F150W. The second colour cut above sets an upper redshift limit of $z \approx 6.5$ where the $Ly\alpha$ break redshifts well into F090W. In these cuts, the F775W flux is set to its $1\sigma$ upper limit in cases where the S/N<1 in that band following past Lyman-break colour cut selections in the literature (e.g. Bouwens et al. 2015; Endsley et al. 2021a).

Additional analytic cuts are imposed to verify the presence of a strong $Ly\alpha$ break. First, we enforce S/N<2 in ACS/F435W since this band lies fully blueward of the Lyman-continuum (LyC) break at $z > 5.5$ and all LyC photons at these redshifts face extremely strong IGM attenuation. Second, we ensure that a strong dropout is also seen in ACS/F606W with either $F606W - F090W > 2.7$, or $F606W - F090W > 1.8$ if the S/N($F606W$)<2, again setting the $F606W$ flux to its $1\sigma$ upper limit when S/N<1. However, we ignore the cut in F435W S/N as well as the F606W dropout criteria for objects with extremely strong $Ly\alpha$ breaks ($F775W - F090W > 2.5$) which can confidently be identified as $z \sim 6$ galaxies even if photometric scatter impacts the F435W or $F606W$ fluxes. Finally, to ensure each source is real, we enforce the criteria that every $z \sim 6$ galaxy in our sample be detected at S/N>5 in at least one NIRCam band, as well as detected at S/N>3 in at least three NIRCam bands in addition to either ACS/F814W or ACS/F850LP.

Our $z \sim 7-9$ Lyman-break selection largely follows the same logic of the $z \sim 6$ criteria described above. Specifically, we begin with the colour cuts

(i) $F090W - F115W > 1.5$;

(ii) $F115W - F200W < 1.2$;

(iii) $F090W - F115W > F115W - F200W + 1.5$.

These cuts are satisfied by galaxies at $z \approx 7.0-9.0$ with blue rest-UV colours ($-2.5 \lesssim \beta \lesssim -2.0$), as well as extremely red ($-1.0 \lesssim \beta \lesssim 0.0$) objects at $z \approx 7-8$ where the Lyman-alpha break has not

yet shifted into F115W. As with the $z \sim 6$ criteria, we also impose the condition that S/N<2 in F435W. For the $z \sim 7-9$ selection, we additionally utilize the $\chi^2_{opt}$ parameter defined in Bouwens et al. (2015) as $\chi^2_{opt} = \sum \text{abs}(f)/f (f/\sigma)^2$ where $f$ and $\sigma$ represent, respectively, the measured flux density and its error in a given band, while $\text{abs}(f)$ is the absolute value of the flux density. After summing over the ACS F435W, F606W, and F775W bands, we enforce $\chi^2_{opt} < 5$. However, we ignore the S/N cut in F435W as well as the $\chi^2_{opt}$ cut for objects with extremely strong $Ly\alpha$ breaks ($F090W - F115W > 2.5$). Every selected $z \sim 7-9$ galaxy must be detected at S/N>5 in at least one NIRCam band redward of F090W, as well as detected at S/N>3 in at least three such bands.

We enforce a final set of cuts to both the $z \sim 6$ and $\sim 7-9$ samples to ensure that the NIRCam data are sufficiently sensitive for rest-optical colour measurements (tracing e.g. Balmer breaks or nebular line EWs). Following the IRAC-based approach of Endsley et al. (2021a), we require that $f(\text{FUV})/e(X) > 3$ where $f(\text{FUV})$ is the far-UV flux density corresponding to the inferred $M_{\text{UV}}$ value from the constant SFH SED fits (see Section 2.4), while $e(X)$ is the $1\sigma$ uncertainty in the flux density of band X. Here, this cut is enforced with the four reddest bands in the JADES filter set (i.e. $X \in \{F335M, F356W, F410M, F444W\}$).

Every selected $z \sim 6-9$ object was visually inspected in all HST and NIRCam mosaics (both at original and PSF-matched resolution) to remove spurious sources due to artifacts (mostly diffraction spikes) or diffuse emission from large, low-redshift objects. We also remove bright point-source objects that show colours consistent with brown dwarfs, which can mimic a strong F090W dropout. The measured colours are compared to the empirical SPEX $0.8-2.5 \ \mu m$ brown dwarf spectral library (Burgasser 2014), as well as the Sonora model brown dwarf spectral templates (Marley et al. 2021). After performing the final steps of neighbour subtraction and combined aperture photometry for the appropriate set of objects (see Section 2.2), we end up with a final sample of 280 F775W dropouts ($z \sim 6$) and 480 F090W dropouts ($z \sim 7-9$) across the $\approx 90 \ \text{arcmin}^2$ area with coverage in ACS F435W, F606W, F775W, F814W, and F850LP as well as NIRCam F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, and F444W.

## 2.4 Photoionization SED modelling

To infer the physical properties of each Lyman-break $z \sim 6-9$ galaxy in our JADES sample, we fit their 14-band ACS+NIRCam photometry with star-forming photoionization models. One of the primary goals of this work is to quantify the stellar masses implied by the rest-UV + optical SEDs. As has been discussed previously (e.g. Carnall et al. 2019; Leja et al. 2019; Lower et al. 2020; Johnson et al. 2021b; Tacchella et al. 2022; Topping et al. 2022a; Endsley et al. 2023; Whitler et al. 2023a, b), stellar mass estimates can change significantly depending on the assumed SFH. We therefore fit every galaxy in our sample with four different sets of models to quantify the systemic stellar mass uncertainties. Each of these four models is described in turn below, though we first note the model assumptions applied in every case.

For all four sets of SED model fits described below, we adopt a Chabrier (2003) stellar IMF with bounds 0.1–300 $M_{\odot}$ as well as an Small Magellanic Cloud (SMC) dust attenuation curve (Pei 1992). Moreover, in each case, we fit with log-uniform priors on stellar mass in the range $5 \le \log(M_*/M_{\odot}) \le 12$, on V-band dust optical depth in the range $-3 \le \log \tau_v \le 0.7$, on ionization parameter in the range $-4 \le \log U \le -1$, and on metallicity in
the range $-2.2 \le \log(Z/Z_\odot) \le -0.3$. The upper limit of $\approx 50$ per cent $Z_\odot$ is set to avoid unphysically high metallicities for the faint ($M_{UV} \gtrsim -20$) reionization-era galaxies that vastly dominate our sample. Moreover, the stellar and ISM metallicities are equivalent in all models considered here. We leave implementations of, for example, a stellar mass to metallicity relation prior or $\alpha$-enhanced metallicity models for future work.

All F775W and F090W dropouts are fit in the redshift range $z = 4-8$ and $6-10$, respectively, with uniform priors adopting the IGM attenuation model of Inoue et al. (2014). We intentionally choose not to make sample cuts based on photometric redshifts as this would likely bias our sample towards $z \sim 6-9$ objects with strong nebular lines (and hence young light-weighted ages) given the unique NIRCam colour patterns caused by [O III] + H$\beta$ and H$\alpha$ lines at high redshifts.

**BEAGLE constant SFH models.** We first consider models adopting a constant SFH (hereafter CSFH) using the BayEsian Analysis of GaLaxy sEds (BEAGLE) SED-fitting code (Chevallard & Charlot 2016). BEAGLE fits photometry against the suite of star-forming photoionization SED models described in Gutkin et al. (2016) which utilize isochrones computed by the PAdova and TRieste Stellar Evolution Code (PARSEC, Bressan et al. 2012; Chen et al. 2015). The Bayesian MULTINEST algorithm (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009) is implemented in BEAGLE to determine the posterior probability distribution for each free physical parameter in the fit. In the BEAGLE CSFH fits, the galaxy age is fit between 1 Myr and the age of the Universe at the sampled redshift with a log-uniform prior, where no star formation is assumed to have occurred prior to the fitted age in these models. Throughout this work, we refer to these BEAGLE CSFH ages as the light-weighted ages of the galaxies.

**BEAGLE two-component SFH models.** We allow for more flexible SFH models within context of BEAGLE following the approach of previous works (e.g. Stark et al. 2017; Endsley et al. 2021a; Castellano et al. 2022). Specifically, we adopt a two-component SFH (hereafter TcSFH) composed of (1) a delayed $\tau$ model component (SFR $\propto t e^{-t/\tau_{SF}}$ where t is the time since the onset of star formation) and (2) a CSFH component. The onset of the delayed $\tau$ model component is assumed to have been between $10^{1.35}$ ($\approx 22$) Myr age up to the age of the Universe at the fitted redshift, while the CSFH component defines the SFH over the most recent 1–20 Myr (log-uniform prior for both components). The SFR from the delayed component is not added to the CSFH component over its fitted recent time interval, allowing for both strong recent increases in the SFR, as well as strong recent drops in the SFR (and anything between). For the delayed $\tau$ component, the $\tau_{SF}$ parameter is fit with a log-uniform prior in the range 1 Myr to 30 Gyr. The SFR of the CSFH component is determined via the specific star formation rate (sSFR) over the associated time interval, where we adopt a log-uniform prior of $-14 \le \log(\text{sSFR/yr}^{-1}) \le -6$.

**PROSPECTOR continuity SFH prior models.** We also consider non-parametric SFH models using the SED-fitting code PROSPECTOR (Leja et al. 2019; Johnson et al. 2021b). PROSPECTOR adopts the star-forming photoionization models from FSPS (Conroy, Gunn & White 2009; Conroy & Gunn 2010) which by default utilizes products from the MESA Isochrones and Stellar Tracks project (MIST; Choi et al. 2016). Posterior probability distributions on fitted parameters are determined using the DYNESTY sampling package (Speagle 2020). We first consider the ‘continuity’ prior implemented in the non-parametric SFH models in PROSPECTOR. This continuity SFH prior weights against strong changes in the SFR over adjacent time bins (see Leja et al. 2019, for details), thereby preferentially yielding more extended SFHs for galaxies with young light-weighted ages.

The implementation of the non-parametric SFH models in this work largely follows that described in Whitler et al. (2023b). To summarize, the SFHs are composed of eight time bins where the SFR in each is a constant value, and the ratios of the SFR in adjacent time bins are fit by PROSPECTOR. The earliest time bin extends to a fitted formation redshift in the range $z_{\text{form}} = 10-30$ (uniform prior), and the two most recent time bins are fixed to 0–3 and 3–10 Myr while the remaining six time bins are divided evenly in logarithmic space to the fitted formation redshift.

**PROSPECTOR bursty SFH prior models.** Finally, we consider the ‘bursty’ SFH prior (Tacchella et al. 2022), a slight modification of the continuity prior described above. These bursty priors more easily allow for strong deviations in the SFR in adjacent time bins while still permitting very extended SFHs like the continuity prior.

For each of the four models described above, fiducial values on inferred physical properties for a given galaxy are taken as the median of the posterior probability distribution. Similarly, the associated $\pm 1\sigma$ error as the inner 68 per cent credible interval from the posterior. Reported absolute UV magnitudes, $M_{UV}$, are computed here from the continuum flux density at 1500 Å rest frame using the output redshift and SED posteriors.

Prior to running the SED fits, we add a 5 per cent systematic error to all photometric measurements, largely with the intention of being conservative about the precision of current photometric zero-points. In general, we find that these star-forming photoionization models yield acceptable fits to the $z \sim 6-9$ Lyman-break galaxies in our sample. For each of the four models sets described above, the median best-fitting $\chi^2$ value over the full sample ranges between 10–14 (14 photometric data points are fit in every case). Moreover, 90 per cent–98 per cent of the SED fits yield best-fitting $\chi^2 < 28$.

However, there are four objects in our sample with measured photometry that are very poorly fit by these star-forming only models. These four objects are those with best-fitting $\chi^2 > 100$ with the BEAGLE CSFH models and best-fitting $\chi^2 > 50$ for the other three models. Two of these objects appear to have steeply rising SEDs between F277W and F444W yet relatively blue colours between F115W and F277W, reminiscent of the ‘Little Red Dot’ (Matthee et al. 2024) population that has garnered much debate within the literature as to whether their light detected in NIRCam is primarily due to emission from stars or AGN (e.g. Akins et al. 2023; Endsley et al. 2023; Greene et al. 2024; Williams et al. 2024; Pérez-González et al. 2024). Both of these systems have been explored in Williams et al. (2024, IDs 121710 and 132229) and Pérez-González et al. (2024), where these works consider both AGN and star-forming SED model fits once adding in deep multiband MIRI imaging. A third object selected in our Lyman-break sample yet very poorly fit with star-forming only models shows very broad H$\alpha$ emission and is a candidate for hosting a black hole merger (ID 101013704 in Maiolino et al. 2023). Because the goal in this paper is to study reionization-era galaxies in context of star-forming models, we remove these four fit objects from our sample. Two of these objects were selected as F775W dropouts, while the other two were selected as F090W dropouts. Therefore, our final sample of Lyman-break $z \sim 6-9$ galaxies analysed below consists of 278 F775W dropouts ($z \sim 6$) and 478 F090W dropouts ($z \sim 7-9$).

## 3 THE PROPERTIES OF GALAXIES AT REDSHIFTS 6–9

In this section, we begin by discussing the overall properties of the sample, beginning with two basic parameters – photometric redshifts and absolute UV magnitudes – and then shift our focus to the inferred
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Distribution of photometric redshifts (left) and absolute UV magnitudes (right) among the sample of 278 F775W dropouts (top) and 478 F090W dropouts (bottom) considered in this work. Here, the adopted value for each galaxy is taken as the median value from the posterior of the BEAGLE CSFH fit outputs. We quote the mean (also shown with a vertical dashed line) and standard deviation of values for each dropout sample in their respective panel. Even though the photometric redshifts of our F090W dropouts extend to $z \approx 8.6$, our selection criteria is sensitive to blue galaxies at redshifts as high as $z \approx 9.0$ where the [O III] + H β emission lines remain in F444W (see Section 2.3).</figcaption></figure>

stellar masses (Section 3.1) and light-weighted ages (Section 3.2). We include a detailed discussion of the most massive objects in our sample, as well as those with the youngest and oldest light-weighted ages.

The F775W and F090W dropout subsets are largely comprised of galaxies with photometric redshifts in the range $5.5 < z_{phot} < 6.5$ and $6.5 < z_{phot} < 8.5$, respectively (see Fig. 1), consistent with expectations given our Lyman-break selection criteria. The mean photometric redshift is $\langle z_{phot} \rangle = 5.9$ and 7.3 for the F775W and F090W dropout subsets, respectively. Considering all 756 of these $z \sim 6-9$ JADES/NIRCam galaxies, our sample spans a very broad range in absolute UV magnitude from $M_{UV} = -22.0$ at the bright end to $M_{UV} = -16.4$ at the faintest end (see Fig. 1). Using a characteristic UV luminosity of $L_{UV}^* = -20.5$ (e.g. Bowler et al. 2017; Harikane et al. 2022), the $M_{UV}$ range of our sample corresponds to $\approx 0.02-4$ $L_{UV}^*$ thus covering a factor of $\approx 200$ in UV luminosity. The average absolute UV magnitude of the F775W dropout ($z \sim 6$) subset is $\langle M_{UV} \rangle = -18.6$, which is slightly brighter than that of the F090W dropout subset ($\langle M_{UV} \rangle = -18.3$) given that the sensitivity of our F775W selection is limited more by the HST/ACS depth. Both dropout subsets have very similar standard deviation in $M_{UV}$ (0.9–1.0 mag) and each include a handful of galaxies at very bright UV luminosities ($M_{UV} < -21.25$ or $L_{UV} > 2 L_{UV}^*$) as well as $\approx 100-200$ galaxies at the very faint end ($M_{UV} > -18$ or $L_{UV} < 0.1 L_{UV}^*$; see Fig. 1).

## 3.1 The stellar masses of Lyman-break $z \sim 6-9$ galaxies

The 756 Lyman-break $z \sim 6-9$ galaxies comprising our JADES sample possess a wide range of inferred stellar masses. From the BEAGLE CSFH models, we infer stellar masses spanning approximately

$3.5$ orders of magnitude from $M_* \approx 2 \times 10^6$ up to $\approx 7 \times 10^9 M_\odot$ (see Figs 2a and b). However, we find that the inferred stellar masses of individual $z \sim 6-9$ galaxies can change substantially depending on the assumed SFH, consistent with the results of previous works (e.g. Tacchella et al. 2022; Topping et al. 2022a; Endsley et al. 2023; Whitler et al. 2023a, b). In our sample, the BEAGLE CSFH stellar masses are, on average, approximately 0.2, 0.1, and 0.5 dex smaller than (i.e. $0.6 \times, 0.8 \times$, and $0.3 \times$) the stellar masses inferred from the BEAGLE CSFH fits, the PROSPECTOR bursty SFH prior fits, and the PROSPECTOR continuity SFH prior fits, respectively (see Figs 2e–g). Therefore, the BEAGLE CSFH fits tend to yield the lowest stellar masses among our sample while the PROSPECTOR continuity SFH prior fits typically yield the largest mass estimates. This is because the bulk of galaxies in our sample have young inferred ages in context of a CSFH ($\sim 50$ Myr; see Section 3.2) and in these models no stars are assumed to have formed prior to that inferred age. On the other hand, the PROSPECTOR continuity prior imposes a preference for significant star formation extending back to $z = 10-30$ ($\Delta t \sim 300-650$ Myr for a galaxy at $z \sim 7$) by weighting against strong time variability in the SFH.

For galaxies with the youngest light-weighted ages in our sample ($\text{age}_{\text{CSFH}} < 10$ Myr), the stellar masses inferred with the PROSPECTOR continuity SFH prior are typically 0.8 dex (i.e. $6 \times$) higher than that inferred with the BEAGLE CSFH setup, though this difference can rise to factors of $\sim 30-100$ in the most extreme cases (see Fig. 2g) as found in Whitler et al. (2023b) using a sample of very bright ($M_{UV} \lesssim -21$) $z \sim 7$ galaxies with IRAC coverage. The offset between constant and continuity prior SFH stellar masses is generally more moderate among galaxies with relatively old light-weighted ages ($\text{age}_{\text{CSFH}} > 100$ Myr), with a median and maximum offset of approximately 0.3 dex ($2 \times$) and 0.7 dex ($5 \times$), respectively (see Fig. 2g), again consistent with the IRAC-based findings of Whitler et al. (2023b). Consequently, we find that the range of stellar masses for our sample shifts upwards to $M_* \approx 7 \times 10^6 - 2 \times 10^{10} M_\odot$ with the PROSPECTOR continuity SFH prior fits (see Figs 2c and d), though we continue to conclude that the faintest ($M_{UV} \sim -17$) and youngest ($\text{age}_{\text{CSFH}} \lesssim 30$ Myr) galaxies tend to have the lowest stellar masses of $\sim 10^7 M_\odot$.

We quantify the relationship between UV luminosity and stellar mass among our sample by fitting the equation $\log(M_*/M_\odot) = a(M_{UV} + 19) + b$. This relation is derived at $z \sim 6$ and $\sim 7-9$ by fitting the F775W and F090W dropout subsets separately, and we moreover perform the fits using the BEAGLE CSFH and PROSPECTOR continuity SFH prior models to bracket the systematic uncertainties associated from different SFH assumptions. In these $M_*-M_{UV}$ fits, we only include objects with relatively bright UV luminosities ($M_{UV} \le -18$) as we aim to mitigate potentially significant biases at the faintest end, though acknowledge that future work directly accounting for incompleteness will be required. The fitted parameters ($a$ and $b$) and their uncertainties are taken as the median and standard deviation of values obtained from 1000 realizations of randomly sampling stellar masses and UV luminosities of each galaxy from the posteriors of their BEAGLE CSFH outputs, again keeping only those with $M_{UV} \le -18$. With this simple approach, we derive the $M_*-M_{UV}$ relations shown in panels Figs 2(a)–(d).

When adopting the BEAGLE CSFH models, we infer a typical stellar mass of $M_* \approx 2.1 \times 10^8$ and $\approx 3.4 \times 10^7 M_\odot$ for $z \sim 6$ galaxies with $M_{UV} = -20$ and $-18$, respectively. At $z \sim 7-9$, the typical inferred CSFH stellar masses are lower by $\approx 0.15-0.2$ dex at these $M_{UV}$ values, resulting in values of $M_* \approx 1.5 \times 10^8$ and $\approx 2.1 \times 10^7 M_\odot$, respectively. This redshift trend is consistent with previous
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: the relationship between absolute UV luminosity and stellar mass for our sample of F775W dropouts (top) and F090W dropouts (bottom). We show this relationship adopting the stellar masses from the BEAGLE CSFH (left panels) and PROSPECTOR continuity SFH prior fits (see Section 2.4). These two mass estimates generally bracket the systematic uncertainties associated with the assumed SFH and we show the fitted relation (dashed line) for each adopting the form log(M*/M⊙) = a(MUV+19) + b and restricting the fit to MUV ≤ −18 (the dashed line shows the relation extrapolated to fainter UV luminosities). Right: illustration of systematic differences in stellar mass estimates from adopting different assumptions on the SFH. We compare stellar mass estimates from the TcSFH BEAGLE fits (top), the bursty SFH prior PROSPECTOR fits (middle), and the continuity SFH prior PROSPECTOR fits (bottom), all relative to the BEAGLE CSFH stellar masses as a function of CSFH age.</figcaption></figure>

findings that, at fixed M<sub>UV</sub>, typical inferred stellar masses decrease slightly with increasing redshift at z ∼ 6 − 9 implying systematically higher sSFRs at earlier epochs (e.g. Stark et al. 2013; Duncan et al. 2014; Song et al. 2016; Bhatawdekar et al. 2019; Stefanon et al. 2021b).

The normalization of our derived M*−M<sub>UV</sub> relation is ≈0.5
dex (i.e. ≈ 3×) higher when adopting the stellar masses from
the PROSPECTOR continuity SFH prior fits, as expected from our
discussion above. The slope of the M*−M<sub>UV</sub> relation we derive in
each case (−0.48 ≤ a ≤ −0.40) is consistent with results of several
previous works (e.g. Stark et al. 2013; Duncan et al. 2014; González
et al. 2014; Song et al. 2016; Bhatawdekar et al. 2019; Kikuchihara
et al. 2020; Stefanon et al. 2021b). However, we are generally
inferring considerably (∼0.3–0.5 dex) lower stellar masses at fixed
M<sub>UV</sub> than found previously with HST and Spitzer data, even when
accounting for differences in the assumed SFHs throughout various
studies. We defer a more detailed comparison with literature results
to a future work, though quickly note that the z ∼ 7 − 9 CSFH M*−
M<sub>UV</sub> relation derived here is consistent with the early JWST results of
Endsley et al. (2023), perhaps reflecting the much richer information
gained from JWST on reionization-era galaxy rest-optical SEDs.

3.1.1 *The most massive Lyman-break z ~ 6 - 9 galaxies in JADES*

Shortly after the first JWST/NIRCam images were released, there
were reports of very high-redshift (z ∼ 7 − 11) candidates with
extremely large stellar masses of M* ∼ 10¹¹M⊙ identified over small
areas (≈40 arcmin²; Labbé et al. 2023) which clearly challenged

models of galaxy formation (e.g. Boylan-Kolchin 2023; Lovell et al.
2023). Since these initial reports,<fnref n="5" /> stellar mass estimates for nearly
all the candidates have been lowered to < 3 × 10¹⁰M⊙, though one
candidate M* = 10¹⁰.⁹M⊙ object remains in the Labbé et al. (2023)
sample. However, due the limited filter set available in these earliest
JWST data, solutions with 10–200× lower stellar masses can be
obtained for the most massive candidate in Labbé et al. (2023)
by adopting different model assumptions that yield similar (if not
considerably better) χ² (see Endsley et al. 2023). Here, we utilize
the deep nine-band NIRCam imaging of JADES (including two long-
wavelength medium bands) to build upon these initial studies by
investigating the potential abundance of extremely massive z > ∼ 6
galaxies.

There are no galaxies in our sample where the photometric data
clearly imply stellar masses of > 3 × 10¹⁰M⊙. Only a small number
of galaxies have inferred stellar masses around ~ 10¹⁰M⊙ when
applying the continuity SFH prior with PROSPECTOR which, as
discussed in Section 3.1, generally yields a maximum stellar mass
estimate for a given galaxy (see Fig. 2). Here, we discuss these most
massive z ~ 6 – 9 candidates in our JADES sample.

There are 13 galaxies with inferred stellar masses of ≥ 3 × 10⁹M⊙ from the PROSPECTOR continuity SFH prior fits among our JADES sample. Eleven and two of these 13 galaxies are in the F775W and F090W dropout subsets, respectively, and all of their SEDs are shown in Figs 3 and 4. Unsurprisingly, many of these objects show

<footnote n="5">⁵https://arxiv.org/abs/2207.12446v2</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The ACS + NIRCam SEDs for 12 of the 13 galaxies with the largest inferred stellar masses in our sample ($M_* \ge 10^{9.5} M_\odot$ from the PROSPECTOR continuity SFH prior fits; the last object is shown in Fig. 4). The measured photometry is shown with diamonds, while the model photometry from the PROSPECTOR continuity SFH prior SEDs is shown with circles. The black lines and corresponding shaded regions show the median posterior from the SEDs and their 68 per cent inner credible intervals, respectively, fixed at the median redshift from the posterior of each galaxy. The PROSPECTOR continuity SFH prior fits generally yield the largest stellar masses (see Fig. 2) and we find a maximum stellar mass of $10^{10.0} M_\odot$ across our full sample with these fits (see Fig. 4).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** The SED fit of the most massive $z \sim 6-9$ Lyman-break galaxy in our sample, JADES-GS + 53.03755 − 27.87491. We show the fit both excluding and including the MIRI F770W photometry for this object, with both adopting the PROSPECTOR continuity SFH prior. The inclusion of the MIRI data results in a consistent (yet slightly lower) stellar mass estimate for this galaxy, yielding a maximum stellar mass of $10^{10.0} M_\odot$ across our full sample (see also JADES-GN+189.13794 + 62.23601 in Fig. 3).</figcaption></figure>

NIRCam colours consistent with strong Balmer breaks implying that a relatively old ($\gtrsim 300$ Myr) stellar population is contributing significantly to the emergent light. However, a few of these galaxies show colours consistent with relatively weak Balmer breaks and are simply so luminous ($-22 \lesssim M_{UV} \lesssim -21.5$) that a substantial population of relatively old stars can be hidden within the observed SED if outshined by more recently formed stars.

The galaxy with the highest inferred stellar mass in the sample (JADES-GS + 53.03755 − 27.87491; $M_* = 10^{10.2\pm0.1} M_\odot$) lies at $z_{phot} \approx 6.1$, is very luminous in the rest UV ($M_{UV} = -21.6$), and shows a significant Balmer break as well as moderate nebular line emission from [O III] + H β and H α (Fig. 4). Fortuitously, this galaxy falls within the deep MIRI F770W parallel imaging of JADES (see Eisenstein et al. 2023) and is clearly detected at high S/N, granting improved constraints on the total amount of stellar mass allowed by the SED by extending measurements into the rest-frame NIR (e.g. Tang, Stark & Ellis 2022; Papovich et al. 2023). With the MIRI measurement included (Alberts et al. 2024), we infer a consistent though slightly ($\approx 0.2$ dex) lower stellar mass of $M_* = 10^{10.0\pm0.1} M_\odot$ (see Fig. 4). Because the F770W flux density (largely probing the rest-NIR continuum) is significantly higher than that in F410M (probing the rest-optical continuum), the model posteriors push to higher metallicity solutions ($\approx 0.3 Z_\odot$) over that without the MIRI data ($\approx 0.05-0.1 Z_\odot$). Such higher metallicity solutions alter the light-to-mass ratio across the rest-frame $\approx 0.4-1$ µm SED, resulting in a slightly lower stellar mass. A more systematic assessment of
how including MIRI data impacts stellar mass inferences among the high-redshift JADES sample will be presented in upcoming works (Florian et al., in preparation; Helton et al., in preparation; Ji et al., in preparation).

Aside from JADES-GS + 53.03755 – 27.87491, there is only one other galaxy in the $z \sim 6-9$ sample considered here with a comparable stellar mass: JADES-GN+189.13794 + 62.23601 with an estimated $M_* = 10^{10\pm0.1} M_\odot$ from the ACS+NIRCam data. JADES-GN+189.13794 + 62.23601 does not possess any MIRI constraints in the rest-frame NIR so we adopt this estimated mass as fiducial (in context of an extended SFH), though we note that previous work has shown that the addition of MIRI data generally lowers estimated stellar masses of high-redshift galaxies (Papovich et al. 2023). Overall, we find that the maximum inferred stellar mass (from the continuity PROSPECT fits) among all galaxies in our sample is $1.0 \times 10^{10} M_\odot$.

We remind the reader that, in this work, we are only considering Lyman-break selected $z \sim 6-9$ galaxies with photometry that can be reasonably explained by star-forming photoionization models. There are four Lyman-break $z \sim 6-9$ objects in JADES with photometry that is very poorly matched by star-forming only models and hence have been excluded from the sample considered here (see Section 2.4). Notably, these include compact objects with NIRCam colours similar to the most massive candidate reported in Labbé et al. (2023) and others subsequently identified in more recent data sets (Furtak et al. 2023a; Akins et al. 2023; Barro et al. 2024), including within JADES (Williams et al. 2024; Pérez-González et al. 2024). The contribution of AGN emission to the observed light in these systems remains heavily debated, and additional data is crucial to better estimate their stellar masses.

## 3.2 The constant SFH ages of Lyman-break $z \sim 6-9$ galaxies

In this subsection, we discuss the light-weighted ages (defined here as CSFH ages) inferred among our sample of 756 Lyman-break galaxies. A wide variety of rest UV + optical SED shapes are clearly found across our sample, indicating a large diversity in light-weighted ages among reionization-era galaxies. The typical galaxy in our sample shows an NIRCam SED consistent with a roughly fixed power-law continuum extending from the rest-UV to rest-optical, implying a relatively weak Balmer break and hence a young light-weighted age ($age_{CSFH} \sim 50$ Myr) consistent with previous findings from early JWST data (e.g. Endsley et al. 2023; Furtak et al. 2023b; Laporte et al. 2023; Leethochawalit et al. 2023; Morishita & Stiavelli 2023; Whitler et al. 2023a). The SEDs for a subset of these typical, young galaxies are shown in Fig. 5 including a variety of redshifts and UV luminosities. Many of the young galaxies show a significant photometric excess in at least one long-wavelength band consistent with strong [O III] + Hβ and/or Hα emission (EW~700–1000 Å; see e.g. IDs JADES-GN+189.19124+62.19952 and JADES-GS + 53.14154 – 27.82320 in Fig. 5) implying a substantial contribution of hot, massive stars to the SED. However, a subset of the young ($age_{CSFH} \sim 50$ Myr) galaxies have measured photometry consistent with surprisingly weak nebular line emission (e.g. IDs JADES-GN+189.19376+62.29430 and JADES-GS + 53.15167 – 27.80925 in Fig. 5); this subpopulation is discussed further in Section 5.2.

In addition to the typical young ($age_{CSFH} \sim 50$ Myr) galaxies, several sources in our sample show long-wavelength photometric excess patterns indicating extremely strong nebular line emission ([O III] + Hβ EW ~ 2000 – 5000 Å or Hα EW ~ 1000 – 2500 Å),

consistent with exceptionally young light-weighted ages ($age_{CSFH} \sim 3$ Myr). Other galaxies show strong Balmer breaks with signatures of weak-to-no nebular line emission ([O III] + Hβ EW~300 Å), consistent with relatively old light-weighted ages ($age_{CSFH} \sim 300–1000$ Myr). These subpopulations representing the extreme ends of the light-weighted age distribution are analysed in greater detail below.

### 3.2.1 Galaxies with the strongest Balmer breaks

The task of identifying reionization-era galaxies with prominent Balmer breaks has long been hindered by the sparse rest-optical SED sampling afforded by Spitzer/IRAC photometry, leading to strong degeneracy with SED solutions of extremely strong nebular line emission (e.g. Schaerer & de Barros 2010; Stark et al. 2013; Roberts-Borsani, Ellis & Laporte 2020). This degeneracy could be alleviated for galaxies situated in specific redshift intervals, leading to the early identification of a small number of $z \sim 7-9$ galaxies showing colours consistent with strong Balmer breaks (e.g. Hashimoto et al. 2018; Strait et al. 2020; Laporte et al. 2021; Endsley et al. 2021a; Tacchella et al. 2022). The list of $z \sim 6-9$ galaxies reported to exhibit prominent Balmer breaks has now steadily grown since the delivery of the first JWST data (e.g. Endsley et al. 2023; Labbé et al. 2023; Laporte et al. 2023), though in some cases the available NIRCam photometry could not confidently rule out solutions with extremely strong nebular line emission. Now equipped with deep imaging in two non-overlapping medium bands as well as three broad bands at $\approx 3-5 \mu$m, we identify and characterize the $z \sim 6-9$ Lyman-break galaxies in our JADES data that confidently show significant Balmer breaks.

For ease of comparison with existing literature (Binggeli et al. 2019), we quantify the strength of the Balmer break as $F_v(4200\AA)/F_v(3500\AA)$ where $F_v$ is the (continuum) flux density at the specified rest-frame wavelength. These flux densities are taken from the posterior SEDs from the BEAGLE fits and are corrected for dust using the inferred $A_V$ value associated with each step of the nested sampling chain. Here, we focus on the subset of galaxies with inferred $F_v(4200\AA)/F_v(3500\AA) > 1.25$ at $> 84$ per cent probability from the posteriors of both the CSFH and TcSFH BEAGLE fits. These objects are classified as those with confident strong Balmer breaks throughout this work.

We identify 13 and 9 galaxies with confident strong Balmer breaks [$F_v(4200\AA)/F_v(3500\AA) > 1.25$] in the F775W and F090W dropout samples, respectively, a subset of which is shown in Fig. 6. Notably, four of these 22 galaxies fall in the very UV-faint regime ($-18 < M_{UV} < -17$; see e.g. JADES-GS + 53.15472 – 27.81561 in Fig. 6) indicating that Balmer breaks do exist among reionization-era galaxies that were previously only identifiable in very deep HST imaging. The full sample of 22 galaxies with confident Balmer breaks spans absolute UV magnitudes of $-20.2 \le M_{UV} \le -17.1$ and redshifts of $z_{phot} = 5.1 - 7.9$. While we do not identify any objects with confident strong Balmer breaks at the brightest UV luminosities ($-22 \le M_{UV} \le -21$), we emphasize that our sample size is very limited at this end of the luminosity function given the search area considered in this work ($\approx 90$ arcmin$^2$). IRAC studies covering $> \text{deg}^2$ fields have revealed $z \sim 7-8$ galaxy candidates in this luminosity regime with potential strong Balmer breaks (e.g. Stefanon et al. 2017, 2019; Topping et al. 2022a; Whitler et al. 2023b).

The JADES galaxies with confident strong Balmer breaks are among those with the oldest light-weighted ages in the sample with
<figure><img src="image_8.png" /><figcaption>**Figure 5.** ACS + NIRCam SEDs for a subset of galaxies with fairly typical light-weighted ages among our sample (age<sub>CSFH</sub> ~ 50 Myr). We show galaxies at bright (M<sub>UV</sub> ≤ −19.5; upper panels) and faint (M<sub>UV</sub> > −18; lower panels) UV luminosities as well as those selected as F775W dropouts (z ~ 6; left panels) and F090W dropouts (z ~ 7 − 9; right panels). The model photometry from the BEAGLE CSFH fits are shown with circles; all quoted properties in this figure are inferred from the BEAGLE CSFH models.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** SEDs for a subset of galaxies that confidently show strong Balmer breaks (F<sub>v</sub>(4200 Å) / F<sub>v</sub>(3500 Å) > 1.25). All quoted properties are inferred from the BEAGLE CSFH fits, though the inferred Balmer breaks remain strong at the >84 per cent credible interval with the TcSFH BEAGLE fits. The figure format is very similar to Fig. 5.</figcaption></figure>

CSFH ages spanning age<sub>CSFH</sub> ≈ 250 − 1000 Myr. In context of these CSFH models, we may expect to see significant rest-optical emission line signatures in the photometry from O stars that have been produced over the past ~10 Myr. Indeed, a subset of galaxies with confident strong Balmer breaks show long-wavelength photometric

excesses implying [O III] + H β EWs ≈ 300–400 Å (e.g. IDs JADES-GN+189.16254+62.25824 and JADES-GS + 53.17954 − 27.77444 in Fig. 6). These emission line signatures imply moderate sSFRs over the most recent 10 Myr of sSFR<sub>10 Myr</sub> ≈ 1 − 10 Gyr<sup>−1</sup> with all four SED fitting procedures described in Section 2.4. However,
<figure><img src="image_10.png" /><figcaption>**Figure 7.** SEDs for a subset of galaxies that confidently exhibit extreme rest-optical line emission ([O III] + H β EW>1500 Å and/or H α EW>800 Å). All quoted properties are inferred from the BEAGLE CSFH fits, though the optical line EWs remain extreme at the >84 per cent credible interval with the TcSFH BEAGLE fits. The figure format is very similar to Fig. 5 and inferred H α EWs are only quoted for objects at z<sub>phot</sub> < 6.5 where NIRCam is sensitive to this line.</figcaption></figure>

we also identify a number of galaxies with strong Balmer breaks showing colours consistent with very weak to no rest-optical line emission ([O III] + Hβ and Hα EWs < 100 Å; e.g. JADES-GN+189.27602+62.19577 and JADES-GS + 53.15472 – 27.81561 in Fig. 6). For these galaxies, the SED fits allowing for more flexible SFHs tend to yield very low recent specific star formation rates (sSFR<sub>10 Myr</sub> < 0.3 Gyr<sup>-1</sup>) due to a declining SFH. The CSFH fits are forced to sSFR<sub>10 Myr</sub> > 1 Gyr<sup>-1</sup> solutions given the age of the Universe at z > 6, and thus always have some modest line emission in the models (minimum [O III] + Hβ EW~200 Å). Due to the very minor effect such weak emission lines can have on the observed broad-band photometry, we cannot reliably conclude whether any of these objects are in fact experiencing a declining SFH. However, NIRSpec spectra has proven to be highly effective at confirming the existence of such systems in the early Universe, with two post-starburst/micro-quenched galaxies now known at z = 5.2 – 7.3 (Strait et al. 2023; Looser et al. 2024). Dedicated follow-up of the candidates identified in this work would help characterize and clarify the abundance of these relatively inactive galaxies during the reionization era.

### 3.2.2 The most extreme line emitters

For about a decade, it has been known that $z > 4$ galaxies typically exhibit high EW (≥500 Å) rest-optical nebular emission lines (i.e. [O III] + Hβ or Hα) from their impact on broadband Spitzer/IRAC colours (e.g. Schaerer & de Barros 2009; Ono et al. 2012; Labbé et al. 2013; Stark et al. 2013; Rasappu et al. 2016). As more observations were dedicated to deep IRAC imaging, it became clear that a surprisingly large number of bright ($M_{UV} \lesssim -20$) $z \sim 7-8$ Lyman-break galaxies exhibited very high EW [O III] + Hβ emission (≥1500 Å, Smit et al. 2014, 2015; Roberts-Borsani et al. 2016; De

Barros et al. 2019; Endsley et al. 2021a; Stefanon et al. 2022a), implying efficient ionizing photon production and very young light-weighted ages. In the past year, JWST imaging has proven to be highly efficient at identifying many more of these $z \gtrsim 6$ extreme line emitters at fainter continuum luminosities, particularly when the nebular lines are situated in a medium band (e.g. Bouwens et al. 2023; Endsley et al. 2023; Rinaldi et al. 2023; Withers et al. 2023). In this subsection, we build upon these previous works by utilizing the deep nine-band NIRCam imaging from JADES (including F335M and F410M) to identify and characterize a statistical sample of $z \sim 6-9$ galaxies which confidently exhibit extreme rest-optical line emission.

Across our full sample of 756 $z \sim 6-9$ Lyman-break GOODS + JADES galaxies, we identify several systems with long-wavelength NIRCam colour patterns indicating extremely high EW rest-optical emission lines. For the purpose of explicitly investigating this population, we restrict our attention to galaxies where the BEAGLE posteriors yield >84 per cent probability of extremely high EWs from both the CSFH and TcSFH fits. By imposing this probability cut on the results from two SFH models we better ensure that the inferred presence of extreme line emission is likely not being confused with a very strong Balmer break.

There are 20 and 29 galaxies that confidently exhibit extreme [O III] + Hβ emission (EW>1500 Å) within our F775W and F090W dropout samples, respectively. We show the SEDs of a subset of these systems in Fig. 7, where objects were chosen to reflect a range of UV magnitudes and redshifts. The confident extreme [O III] + Hβ emitters in GOODS + JADES have absolute UV magnitudes spanning $-21.0 \le M_{UV} \le -16.9$ and thus encompass nearly the entire $M_{UV}$ range of our full $z \sim 6-9$ sample. Remarkably, this includes galaxies in the very UV-faint regime ($M_{UV} > -18$; see e.g. JADES-GS + 53.17835 – 27.77879 and JADES-GS + 53.19590 – 27.79240 in Fig. 7). The vast majority of these extreme [O III] + Hβ emitters
show colour patterns indicating they lie at redshifts consistent with the targeted windows of our dropout selections ($z_{phot} \approx 5.5-8.4$). Only two outlier sources show strong excesses in F277W suggesting $z_{phot} \sim 5$ which presumably entered our $z \sim 6$ Lyman-break selection due to photometric noise in F775W and/or F090W.

The confident extreme [O III] + Hβ emitters all show SEDs consistent with exceptionally young light-weighted ages (age<sub>CSFH</sub> ≈ 3 Myr) where the emergent light is heavily dominated by both stellar and nebular emission powered by recently formed O stars. Several of these galaxies exhibit a significant drop in flux density between the short-wavelength NIRCam bands and at least one of the long-wavelength medium bands implying a significant Balmer jump between the rest-UV and optical continuum (see e.g. IDs JADES-GS + 53.16900 - 27.80079 and JADES-GN+189.33577 + 62.18652 in Fig. 7). Such Balmer jumps are caused by strong nebular continuum emission which is most prominent at very young ages and low metallicities (e.g. Byler et al. 2017; Topping et al. 2022b). It is only for these objects that the inferred [O III] + Hβ EWs confidently reach values of >3000 Å (up to ≈5000 Å) as the data give direct evidence for weak rest-optical continuum emission yet extremely strong line emission from excesses in other bands.

The very young light-weighted ages indicate that these extreme [O III] + Hβ emitters have recently experienced a dramatic rise in SFR (e.g. Tang et al. 2019, 2022; Endsley et al. 2021a; Carnall et al. 2023; Tacchella et al. 2023b). In context of a CSFH, we infer an extremely large typical sSFR<sub>10 Myr</sub> ~ 300 Gyr<sup>-1</sup> for these extreme [O III] + Hβ emitters, which is simply the inverse of their typical light-weighted age (age<sub>CSFH</sub> ~ 3 Myr). But even with the PROSPECTOR continuity SFH prior fits (which generally yield the lowest sSFR estimates; see Section 3), we continue to infer a very high median sSFR averaged over the past 10 Myr for these systems (26 Gyr<sup>-1</sup>) with values up to 72 Gyr<sup>-1</sup>.

At $z \sim 6$, the JADES/NIRCam photometry provides constraints on not only the EW of [O III] + Hβ emission, but also Hα. We identify 31 galaxies within our F775W dropout subsample that have a >84 per cent posterior probability of Hα EW>800 Å and $z < 6.5$ from both the CSFH and two-component BEAGLE SFH fits. Here, we have added the $z < 6.5$ criteria as the Hα line begins to redshift out of F444W at higher redshifts. This sample of confident extreme Hα emitters overlaps significantly with the confident extreme [O III] + Hβ emitters (20 in both subsets) and thus has similar characteristics. The confident Hα emitters span absolute UV magnitudes of $-21.0 \le M_{UV} \le -17.4$ and light-weighted ages of age<sub>CSFH</sub> ≈ 1–10 Myr. A subset show colours implying prominent Balmer jumps in their continuum and have inferred Hα EWs confidently reaching values up to ≈2000–3000 Å (see e.g. JADES-GS + 53.16900 - 27.80079 in Fig. 7). Their typical specific star formation rates range from sSFR<sub>10 Myr</sub> = 150 Gyr<sup>-1</sup> with the BEAGLE CSFH fits down to sSFR<sub>10 Myr</sub> = 14 Gyr<sup>-1</sup> with the PROSPECTOR continuity prior SFH models. We place both the extreme [O III] + Hβ and extreme Hα emitter populations in context of the broader sample in the following section.

## 4 THE DEPENDENCE OF [O III] + Hβ AND Hα EQUIVALENT WIDTH ON UV LUMINOSITY

The EWs of rest-optical nebular lines have long provided key insight into the physical properties of reionization-era galaxies (e.g. Schaerer & de Barros 2009; Stark et al. 2013; Smit et al. 2014). Previous studies have found that relatively bright ($M_{UV} \lesssim -20$) $z \sim 7$ systems commonly show far higher [O III] + Hβ EWs than

galaxies even at $z \sim 2$ (e.g. De Barros et al. 2019; Endsley et al. 2021a; Boyett et al. 2022), implying generally less chemically enriched gas and more vigorous star formation at earlier times. However, the nature of very UV-faint ($M_{UV} \gtrsim -18$; $L_{UV} \lesssim 0.1 L_{UV}^*$) $z \gtrsim 6$ galaxies has remained far less clear given the difficulties in obtaining sensitive rest-UV + optical SED measurements among a statistical sample of such systems. Photometric constraints on the nebular emission line EWs of this relatively numerous galaxy population would provide a powerful first glimpse into their physical properties, as these lines are sensitive to several physical parameters including metallicity, SFH, and ionizing photon production and escape.

With our JADES/NIRCam sample of 756 $z \sim 6-9$ Lyman-break galaxies, we now investigate how the reionization-era [O III] + Hβ and Hα EW distributions correlate with UV luminosity over a large dynamic range ($-22 \lesssim M_{UV} \lesssim -16.5$). We focus on investigating how the EWs trend with UV luminosity given that the galaxies with the lowest stellar masses in our sample are biased towards those with the youngest light-weighted ages (see Fig. 2), and thus will be weighted towards higher EWs. We first describe the methodology for the EW distribution inference as a function of $M_{UV}$ in Section 4.1, following with the results on the [O III] + Hβ (Section 4.2) and Hα (Section 4.3) EW distributions.

### 4.1 Methodology

Here, we describe how we quantify the UV luminosity dependence on the [O III] + Hβ and Hα EW distributions in the reionization era. For this analysis, we divide our JADES/NIRCam sample into three different subsets separated by absolute UV magnitude: a bright ($M_{UV} \le -19.5$) subset, a faint ($-19.5 < M_{UV} \le -18$) subset, and a very faint ($M_{UV} > -18$) subset. We also aim to test whether there is any significant redshift evolution in the [O III] + Hβ EW distribution during the reionization era and thus infer this distribution for the F775W dropouts ($z \sim 6$) and F090W dropouts ($z \sim 7-9$), separately. Because the Hα line redshifts out of the reddest JADES/NIRCam band (F444W) at $z > 6.6$, we restrict our analysis of the Hα EW distribution to the F775W dropout sample.

To infer the EW distribution among a given subpopulation of our sample, we follow the Bayesian formalism described in Schenker et al. (2014, see also Endsley et al. 2021a; Boyett et al. 2022):

$$P(\theta) \propto \prod_i \int P_i(\text{EW}) P(\text{EW}|\theta) \, \text{dEW}. \quad (1)$$

In this equation, $\theta$ encapsulates the parameters describing the assumed functional form of the distribution (e.g. the median and variance for a lognormal distribution) while $P_i(\text{EW})$ is the probability distribution function (PDF) on EW for the $i$th galaxy in the subpopulation of interest. This formalism has two highly desired features. First, the expression in the integral of equation (1) is small when the assumed EW distribution ($P(\text{EW}|\theta)$) is inconsistent with the PDF of galaxy $i$, while the expression is large when the assumed distribution and PDF strongly overlap. Second, galaxies with narrow (i.e. precise) PDFs will have much stronger constraining power on the inferred $P(\theta)$ than galaxies with very broad (i.e. largely unconstrained) PDFs.

The PDFs for each galaxy are taken from the posterior of the SED fits, but are explicitly corrected for the prior on EW imposed by the SED fitting approach of interest (i.e. $P_i(\text{EW}) \equiv \text{posterior}_i(\text{EW}) / \text{prior}(\text{EW})$). In doing so, we ensure that our resulting inferred EW distributions are being driven by evidence from the photometry rather than the choice of priors in the SED
modelling. These model priors have a stronger influence on the output posterior PDFs for galaxies with lower S/N data, and hence would systematically impact the fainter populations more if we did not divide out the priors when computing $P_i$ (EW). To obtain the EW priors for each SED modelling setup, we fit a single photometric data point in the rest-UV thereby allowing the sampling algorithms to completely explore the allowed parameter space on [O III] + H β and H α EW.

In practice, we treat the integral in equation (1) as a Riemann sum to avoid having to approximate $P_i$ (EW) as some functional form via, for example, spline interpolation. Therefore, we rewrite equation (1) as

$$P(\theta) \propto \prod_i \left[ \sum_j P_{i,j}(\text{EW}) P_j(\text{EW}|\theta) \right] \quad (2)$$

where the index j represents a bin in EW. $P_{i,j}$ (EW) is the probability that object i has an EW in bin j, while $P_j$ (EW|θ) is the integrated probability of the assumed functional EW distribution within bin j. Therefore, both $P_{i,j}$ (EW) and $P_j$ (EW|θ) implicitly account for the width of bin j.

For the [O III] + H β EW inference, we divide the EW bins by values of log(EW/Å) = 2.5-3.7 with spacing of 0.05 dex, yielding 26 total bins (i.e. log(EW/Å) = [≤2.50, 2.50-2.55, 2.55-2.60,..., 3.65-3.70, ≥3.70]). We chose a lower bound of log(EW/Å) = 2.5 to reflect the fact that the photometric data can only place an effective upper limit on the inferred EWs of ≈300 Å for galaxies with the lowest S/N of ~3-5 in the long-wavelength NIRCam bands. The upper bound of log(EW/Å) = 3.7 was chosen given that the maximum inferred [O III] + H β EW in our sample is ≈5000 Å (see Section 3.2.2). Similarly, the H α EW bins applied in equation (2) are divided by values of log(EW/Å) = 2.5-3.5 with spacing of 0.05 dex, where the upper bound of log(EW/Å) = 3.5 reflects the maximum inferred H α EW of ≈3000 Å in the sample. We note that, with this approach, the inferred EW distribution values at ≲300 Å are effectively constrained by the high-EW tail of the distribution given that we assume a log-normal functional form.

For each galaxy subpopulation, we infer the EW distribution assuming a lognormal functional form defined by two parameters: the median of the EW distribution, $\mu_{EW4pt}$, and its standard deviation, $\sigma_{EW4pt}$. This assumption of a lognormal EW distribution is motivated by spectroscopic results at lower redshifts (Lee et al. 2007, 2012; Ly et al. 2011). Finally, we infer the EW distributions using both the BEAGLE CSFH and TcSFH models. Due to our findings in Appendix B and Section 5.2 (see also Fig. 14 and Fig. B1), we adopt the results from the TcSFH models as fiducial.

## 4.2 Results on [O III] + H β EW distributions

We begin by quantifying the inferred [O III] + H β EW distribution among the brightest galaxies in our sample, comparing to existing results from the literature. We then proceed to describe our findings for the very faint population that is now possible with the depth and area of JADES/NIRCam imaging. From our fiducial TcSFH BEAGLE SED fits, we infer a median [O III] + H β EW of $\mu_{EW4pt} = 740_{-80}^{+90}$ Å among the bright ($\langle M_{UV} \rangle = -20.1$) subset at $z \sim 7-9$. This is consistent with median [O III] + H β EWs reported in the literature for similarly UV-luminous HST and NIRCam selected galaxy samples at $z \sim 7-8$ (Labbé et al. 2013; De Barros et al. 2019; Endsley et al. 2021a, 2023; Stefanon et al. 2022a), as well as samples of very UV-bright ($-22.5 \lesssim M_{UV} \lesssim -21.5$) $z \sim 7$ galaxies selected from ground-based imaging (Endsley et al. 2021a, 2023). The inferred

standard deviation of the EW distribution among our bright $z \sim 7-9$ JADES/NIRCam subset ($0.32_{-0.04}^{+0.04}$ dex) is also consistent with values reported in the literature at similar luminosities (Endsley et al. 2021a, 2023).

Due to the strong degeneracy between Balmer breaks and nebular line emission on IRAC colours at $z \sim 6$, the [O III] + H β EW distribution has never been quantified in the literature at this epoch. With the rich (five-band) 3-5 µm SED sampling now afforded by JADES/NIRCam, we here infer that the median [O III] + H β EW among bright ($\langle M_{UV} \rangle = -20.0$) $z \sim 6$ galaxies is $\mu_{EW4pt} = 890_{-80}^{+90}$ Å, ≈0.08 dex ($\sim 1.2 \times$) higher than that inferred from our bright $z \sim 7-9$ subsample. The standard deviation of [O III] + H β EWs among the bright $z \sim 6$ subset is inferred to be very similar to that of the bright $z \sim 7-9$ subset (see Table 6).

Now equipped with a sample of over 600 Lyman-break $z \sim 6-9$ galaxies with ≈2-30× fainter UV luminosities ($-19.5 \le M_{UV} \lesssim -16.5$), we investigate whether the [O III] + H β EW distribution correlates significantly with UV luminosity in the reionization era. Using the fiducial TcSFH BEAGLE SED fits, we infer a median [O III] + H β EW of $\mu_{EW4pt} = 590_{-50}^{+40}$ and $380_{-60}^{+80}$ Å among the faint ($\langle M_{UV} \rangle = -18.7$) and very faint ($\langle M_{UV} \rangle = -17.5$) $z \sim 6$ subsets, respectively. This indicates a smooth and strong (factor of ≈2.3) decline in the typical [O III] + H β EW of $z \sim 6$ galaxies between the bright and very faint subsets, which are separated by an order of magnitude in average UV luminosity (see Fig. 8). Similarly, at $z \sim 7-9$, we infer a median [O III] + H β EW of $\mu_{EW4pt} = 500_{-40}^{+50}$ and $300_{-40}^{+40}$ Å among the faint ($\langle M_{UV} \rangle = -18.6$) and very faint ($\langle M_{UV} \rangle = -17.6$) subsets, respectively, indicating a ≈2.5× decline in typical [O III] + H β EW over ≈1 dex in $L_{UV}$ (Fig. 8). Therefore, in both the $z \sim 6$ and $\sim 7-9$ samples, we infer a strong UV luminosity dependence on the median [O III] + H β EW at $z \sim 6-9$ with high confidence ($\gtrsim 5\sigma$).

To help visualize these trends, we plot the fiducial [O III] + H β EWs (the median of the posterior probability distribution) for each individual galaxy as a function of $M_{UV}$ in Figs 9(a) and (b). In the UV-bright ($M_{UV} < -19.5$) bins, only a small fraction ($\approx 10$ per cent–15 per cent) of $z \sim 6-9$ LBGs have inferred [O III] + H β EWs of $\lesssim 300$ Å, while this fraction rises to ≈50 per cent–65 per cent in the very UV-faint ($-18 \le M_{UV} \lesssim -16.5$) bins. Therefore, even when just considering the fiducial EWs of each galaxy (i.e. not accounting for the priors or uncertainties on EW), we find that the median [O III] + H β EW is ≲300 Å at very low UV luminosities. With this very simple approach, we obtain much higher median [O III] + H β EWs of 880 and 700 Å in the UV-bright bins at $z \sim 6$ and $\sim 7-9$, respectively. None the less, we emphasize that the EW distribution inference approach (Section 4.1) should yield more robust constraints on the median EWs given that this method accounts for the SED-fitting priors and EW uncertainties for each object.

To understand why the NIRCam SEDs are driving the models to yield a strong decline in [O III] + H β EW with decreasing UV luminosity, we analyse the long-wavelength NIRCam colours sensitive to these lines. We begin by considering the F775W dropout ($z \sim 6$) sample where the large majority of galaxies lie at $5.5 < z < 6.5$ (see Fig. 1). For this sample, the F410M band provides a relatively clean probe of the rest-optical continuum given that it is free of contamination from the strongest rest-optical lines ([O III], H β, and H α) at $z \approx 5.6-6.7$. Meanwhile, the F356W and F335M bands are contaminated by [O III] + H β at redshifts associated with the majority of the $z \sim 6$ subset ($z \approx 5.4-6.9$ and ≈5.5–6.1, respectively). Therefore, we investigate whether the measured photometry of the $z \sim 6$ sample supports significant flattening in the F410M–F335M and F410M–F356W colours as
<table><caption>**Table 2.** Table summarizing the inferred [O III] + H β and H α EW distributions among our sample of 756 Lyman-break $z \sim 6-9$ galaxies. The EW distributions are assumed to be lognormal and are parametrized by the median EW, $\mu_{EW4pt}$, and standard deviation, $\sigma_{EW4pt}$. We divide our sample into six subpopulations divided by F775W dropouts ($z \sim 6$) and F090W dropouts ($z \sim 7-9$), as well as into bright ($M_{UV} < -19.5$), faint ($-19.5 < M_{UV} < -18$), and very faint ($M_{UV} > -18$) subsets. For each of the six subpopulations, we infer the EW distribution from the posteriors of the BEAGLE CSFH and TcSFH models. We limit the H α EW distribution results to the $z \sim 6$ sample given that H α redshifts out of NIRCam at $z > 6.5$. We recommend using the EW distributions derived from the TcSFH model outputs given our findings in Section 5.2 (see Fig. 14).</caption><thead><tr><th>Redshift</th><th>M<sub>UV</sub> bin</th><th>SED model</th><th>⟨M<sub>UV</sub>⟩</th><th>N<sub>gal</sub></th><th>μ<sub>EW4pt</sub> ([O III] + H β) [Å]</th><th>σ<sub>EW4pt</sub> ([O III] + H β) [dex]</th><th>μ<sub>EW4pt</sub> (Hα) [Å]</th><th>σ<sub>EW4pt</sub> (Hα) [dex]</th></tr></thead><tbody><tr><td rowspan="6">z ~ 6</td><td rowspan="2">Bright</td><td>BEAGLE CSFH</td><td>-20.0</td><td>64</td><td>930<sup>+70</sup><sub>-80</sub></td><td>0.26<sup>+0.03</sup><sub>-0.03</sub></td><td>760<sup>+50</sup><sub>-30</sub></td><td>0.17<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>BEAGLE TcSFH</td><td>-20.0</td><td>65</td><td>890<sup>+90</sup><sub>-80</sub></td><td>0.31<sup>+0.03</sup><sub>-0.03</sub></td><td>710<sup>+50</sup><sub>-60</sub></td><td>0.24<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td rowspan="2">Faint</td><td>BEAGLE CSFH</td><td>-18.7</td><td>138</td><td>680<sup>+50</sup><sub>-50</sub></td><td>0.33<sup>+0.03</sup><sub>-0.03</sub></td><td>740<sup>+30</sup><sub>-30</sub></td><td>0.18<sup>+0.01</sup><sub>-0.02</sub></td></tr><tr><td>BEAGLE TcSFH</td><td>-18.7</td><td>138</td><td>590<sup>+40</sup><sub>-50</sub></td><td>0.37<sup>+0.04</sup><sub>-0.03</sub></td><td>580<sup>+40</sup><sub>-30</sub></td><td>0.26<sup>+0.01</sup><sub>-0.02</sub></td></tr><tr><td rowspan="2">Very faint</td><td>BEAGLE CSFH</td><td>-17.4</td><td>76</td><td>460<sup>+80</sup><sub>-80</sub></td><td>0.49<sup>+0.09</sup><sub>-0.07</sub></td><td>850<sup>+40</sup><sub>-40</sub></td><td>0.14<sup>+0.03</sup><sub>-0.04</sub></td></tr><tr><td>BEAGLE TcSFH</td><td>-17.5</td><td>75</td><td>380<sup>+80</sup><sub>-60</sub></td><td>0.51<sup>+0.08</sup><sub>-0.07</sub></td><td>630<sup>+60</sup><sub>-60</sub></td><td>0.31<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td rowspan="6">z ~ 7-9</td><td rowspan="2">Bright</td><td>BEAGLE CSFH</td><td>-20.1</td><td>52</td><td>780<sup>+70</sup><sub>-70</sub></td><td>0.26<sup>+0.04</sup><sub>-0.03</sub></td><td>—</td><td>—</td></tr><tr><td>BEAGLE TcSFH</td><td>-20.1</td><td>55</td><td>740<sup>+90</sup><sub>-80</sub></td><td>0.32<sup>+0.04</sup><sub>-0.04</sub></td><td>—</td><td>—</td></tr><tr><td rowspan="2">Faint</td><td>BEAGLE CSFH</td><td>-18.6</td><td>218</td><td>580<sup>+60</sup><sub>-40</sub></td><td>0.42<sup>+0.04</sup><sub>-0.03</sub></td><td>—</td><td>—</td></tr><tr><td>BEAGLE TcSFH</td><td>-18.6</td><td>223</td><td>500<sup>+50</sup><sub>-40</sub></td><td>0.45<sup>+0.04</sup><sub>-0.03</sub></td><td>—</td><td>—</td></tr><tr><td rowspan="2">Very faint</td><td>BEAGLE CSFH</td><td>-17.5</td><td>208</td><td>330<sup>+50</sup><sub>-40</sub></td><td>0.43<sup>+0.07</sup><sub>-0.05</sub></td><td>—</td><td>—</td></tr><tr><td>BEAGLE TcSFH</td><td>-17.6</td><td>200</td><td>300<sup>+40</sup><sub>-40</sub></td><td>0.48<sup>+0.06</sup><sub>-0.05</sub></td><td>—</td><td>—</td></tr></tbody></table>

we move from the bright to very faint subsets. In the bright $z \sim 6$ subset, the median F410M–F335M and F410M–F356W colours are $0.98^{+0.08}_{-0.10}$ and $0.63^{+0.02}_{-0.07}$ mag, respectively, consistent with typical contamination from strong [O III] + H β emission. These colours are typically flatter among the very faint subpopulation ($0.67^{+0.10}_{-0.05}$ and $0.54^{+0.06}_{-0.13}$ mag, respectively), consistent with generally weaker [O III] + H β EWs at lower UV luminosities. The uncertainties on median colours are derived using bootstrap resampling.

We conduct a similar colour analysis with the F090W dropout sample. Across nearly the entire redshift range of our F090W dropout selection, the rest-optical continuum is probed free of strong lines by either F335M ($z \approx 6.3-8.2$) or F356W ($z \approx 7.2-8.8$), while the [O III] + H β emission is captured by either F410M ($z \approx 7.0-7.6$) or F444W ($z \approx 7.0-9.0$). We therefore investigate whether the measured photometry supports significant flattening in the F335M–F410M, F335M–F444W, F356W–F410M, and F356W–F444W colours as we move from the bright to very faint subsets of our $z \sim 7-9$ sample. In the bright subset, the median values for all four above colours are noticeably red ($0.49^{+0.06}_{-0.12}$, $0.34^{+0.03}_{-0.07}$, $0.46^{+0.03}_{-0.09}$, and $0.27^{+0.08}_{-0.02}$, respectively), consistent with typical contamination of strong [O III] + H β emission. On the other hand, in the very faint subset we measure much flatter median colours ($0.20^{+0.04}_{-0.06}$, $0.00^{+0.05}_{-0.02}$, $0.23^{+0.04}_{-0.06}$, and $-0.01^{+0.08}_{-0.03}$, respectively) suggestive of typically much weaker [O III] + H β EWs.

Because the sample size of each $M_{UV}$ bin is fairly large ($N_{gal} = 55-200$; see Table 2), we do not expect these median colours to be strongly impacted by noise. None the less, we employ a Monte Carlo simulation to quantify the likelihood of measuring the median colours of the very faint subsets (accounting for noise) under the assumption that they indeed have true colours represented by their respective bright subset of galaxies. For this test, we take the measured photometry of all galaxies in the bright subset for a given dropout sample (randomly sampled with replacement), uniformly renormalize the photometry of each bright galaxy such that it has an $M_{UV}$ value equivalent to a galaxy in the very faint subset in the

same dropout sample (again randomly sampled with replacement), and add Gaussian noise to their photometry in each band using the uncertainties from the associated very faint galaxy. This Monte Carlo simulation is performed 5000 times and, in each iteration, we compute the median colours of interest.

The distributions of median colours from the Monte Carlo simulations are shown in Fig. 10, along with the actual measured median colours of the very faint subsets plotted as vertical dashed lines in each panel. For nearly every colour probing [O III] + H β emission, we find that very few ($\le 0.1$ per cent) of the 5000 Monte Carlo iterations yield median recovered colours as flat (or flatter) than that measured among the very faint subset. While the median F410M–F356W Monte Carlo colours in the F775W dropout subset are consistent (at the $\approx 1\sigma$ level) with that measured among the very faint subset, F356W is less sensitive to strong [O III] + H β emission relative to F335M given its broader bandwidth and thus the significance in colour difference is expected to be considerably weaker. Moreover, the smaller size of the very faint F775W dropout sample makes its median colour uncertainties $\approx 2 \times$ higher relative to those of the F090W dropout sample (see Fig. 10). This higher noise likely contributes to the lack of a strong difference in the median F410M–F356W Monte Carlo colours relative to that measured for the very faint F775W dropout sample.

Having now validated that the inferred trend in median [O III] + H β EW with UV luminosity is consistent with expectations given the measured long-wavelength colours, we consider whether the data imply any significant changes in the width of the EW distribution towards fainter UV luminosities. Even though the median EW of the very faint population declines substantially relative to the brightest $M_{UV}$ bin, we continue to identify a number of very faint galaxies with prominent long-wavelength colours implying extremely high [O III] + H β EWs of $>1500$ Å (see e.g. ID JADES-GS + 53.17835 – 27.77879 and JADES-GS + 53.19590 – 27.79240 in Fig. 7). Accordingly, our assumed lognormal EW distributions are inferred to broaden considerably towards fainter UV luminosities at
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The inferred distributions of [O III] + H β EW (top) and H α EW (bottom) among our sample using the TcSFH BEAGLE SED fit outputs. The [O III] + H β EW distributions are inferred for both the *z* ~ 6 and ~ 7 − 9 subsets, while the H α EW distribution is inferred only at *z* ~ 6 given this line redshifts out of NIRCam at *z* > 6.5. The different curves show the EW distribution inferred among three different UV luminosity bins (see Table 2). We find a strong, highly significant decline in the typical [O III] + H β EW with UV luminosity, yet a very weak UV luminosity dependence on H α EW (at least at *z* ~ 6).</figcaption></figure>

both *z* ~ 6 and ~ 7, with σ<sub>EW4pt</sub> increasing from ≈ 0.30±0.03 dex in the bright subsets to ≈ 0.50±0.06 dex in the very faint subsets (see Table 2). From these distributions, we infer that the fraction of galaxies with extremely high [O III] + H β EWs (>1500 Å) declines significantly towards fainter UV luminosities, going from 24<sup>+5</sup><sub>-4</sub> per cent (18<sup>+5</sup><sub>-4</sub> per cent) in the bright *z* ~ 6 (*z* ~ 7 − 9) subset to 13<sup>+4</sup><sub>-3</sub> per cent (8<sup>+2</sup><sub>-2</sub> per cent) in the respective very faint subset. Conversely, the fraction of galaxies with relatively very low [O III] + H β EWs (<300 Å) rises dramatically towards the faint-end of the luminosity function, increasing from 6<sup>+3</sup><sub>-2</sub> per cent (11<sup>+5</sup><sub>-4</sub> per cent) in the bright *z* ~ 6 (*z* ~ 7 − 9) subset to 42<sup>+7</sup><sub>-7</sub> per cent (50<sup>+5</sup><sub>-5</sub> per cent) in the respective very faint subset.

4.3 Results on Hα EW distributions

Due to the poor rest-optical SED sampling of IRAC photometry, previous efforts to quantify H α EWs at high redshifts have primarily been limited to the redshift range z ≈ 4.0 − 5.0 where several works have inferred typical H α EWs of ∼400 Å among HST-selected galaxies (Stark et al. 2013; Rasappu et al. 2016; Smit et al. 2016; Faisst et al. 2019; Lam et al. 2019). Most recently, Stefanon et al. (2022b) stacked the IRAC 5.8 µm data of ≈100 z ∼ 8 galaxies and recovered a considerably larger typical H α EW (≥1000 Å), perhaps implying a significant redshift evolution in the median H α EW at high redshifts. With the JADES NIRCam data, we are now able to infer the H α EW distribution at z ∼ 6 and moreover investigate whether it correlates strongly with UV luminosity as may be expected given our results on the [O III] + H β EW distribution above.

From the BEAGLE TcSFH SED fits, we infer very similar median H α EWs among all three UV luminosity bins at z ~ 6 (≈600–700 Å in each; see Table 2). Combined with the literature results discussed above, this implies that the typical H α EWs increase by ≈0.2 dex (≈ 1.6×) between z ~ 4.5 and ~ 6. But perhaps the most striking result here is that we find no evidence for a strong change in the typical z ~ 6 H α EW across 1 dex in UV luminosity, in marked contrast with our findings on the [O III] + H β EW distribution in Section 4.2. Even when just using the fiducial H α EWs for each object (Fig. 9c), we find that the median H α EW lies in a narrow (0.1 dex) range of 530–670 Å for every M$_{UV}$ bin.

We again verify that the measured long-wavelength NIRCam colours are consistent with little-to-no change in the median H α EW between the bright and very faint z ~ 6 subsets following the same approach as described in Section 4.2. At z ~ 6, the F410M band provides a relatively clean probe of the rest-optical continuum (see Section 4.2), while the F444W band contains H α . The median F410M − F444W colour of the bright z ~ 6 subset (0.35<sup>+0.05</sup><sub>-0.04</sub> mag) is nearly identical to that of the very faint subset (0.33<sup>+0.07</sup><sub>-0.14</sub> mag), consistent with a weak trend in typical H α EW with UV luminosity. The Monte Carlo simulation test described in Section 4.2 also reveals that even accounting for S/N differences in the bright versus very faint z ~ 6 subsets often yields nearly identical median F410M − F444W colours between the two populations (see Fig. 10). We thus conclude that the photometry is consistent with little change in typical z ~ 6 H α EWs over the UV luminosity range spanned by our sample (−22 ≤ M$_{UV}$ ≤ −16.4).

There is moderate evidence (≈ 2σ) for a slight (≈0.05 dex) increase in the width of the z ~ 6 H α EW distribution towards very faint UV luminosities (see Table 2). This is consistent with the slightly larger fraction of very faint z ~ 6 galaxies with fiducial H α EWs <300 Å (see Fig. 9c) even though the median EW remains effec-
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The inferred fiducial EWs of each individual galaxy in our $z \sim 6-9$ LBG sample plotted as a function of $M_{UV}$. These fiducial EWs correspond to the median of the posterior probability distribution from the BEAGLE TcSFH fits. Those with $EW<300$ Å are plotted as upper limits at that value (see Section 4.1), and we quote the fraction of such objects in each $M_{UV}$ bin. The trends obtained from these fiducial EWs alone (i.e. without correcting for the SED-fitting priors or directly accounting for object uncertainties) are broadly consistent with the results of the more robust EW distribution inference method (Fig. 8).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Illustration of results from Monte Carlo simulation test validating that the substantial decline in median [O III] + Hβ EW at fainter UV luminosities is very unlikely to be caused by photometric noise. The shaded histograms show the distribution of recovered median colours (one colour combination per panel) from 5000 iterations of degrading the photometry of objects in our brightest ($M_{UV} < -19.5$) bin to the S/N of galaxies in the faintest bin ($M_{UV} > -18$). For nearly all colours sensitive to [O III] + Hβ contamination, very few ($\le 0.1$ per cent) of the 5000 Monte Carlo iterations yield median recovered colours as flat (or flatter) as that measured among the very faint subset (vertical dashed lines). The median F410M–F444W colour of the very faint $z \sim 6$ subset are however consistent with that recovered from the bright subset (after S/N degradation), implying little to no UV luminosity dependence on the median Hα EW.</figcaption></figure>

tively unchanged. None the less, overall the inferred Hα EW distributions for the three UV luminosity bins are nearly indistinguishable (see Fig. 8). If we assume that the Hα EW distribution can be described as a single lognormal form across the entire UV luminosity range of our $z \sim 6$ sample ($-22 \le M_{UV} \le -16.4$), we obtain $\mu_{EW4pt} = 600_{-10}^{+40}$ Å and $\sigma_{EW4pt} = 0.27 \pm 0.01$ dex. Given the strong UV luminosity dependence on [O III] + Hβ EW inferred in Section 4.2, it is no surprise that we then infer a strong UV luminosity dependence on the [O III]λ5007/Hα ratio (hereafter [O III]/Hα) as well. In the bright $z \sim 6$ subset, we infer a median [O III]/Hα ≈ 1.6 while in the faint and very faint subsets we infer median [O III]/Hα ratios of ≈1.3 and ≈0.9, respectively.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Illustration of how the EWs of [O III] + H β (blue solid lines) and H α emission correlate with metallicity (panel a), recent SFH, quantified by the ratio of SFR averaged over the past 3 Myr divided by the SFR averaged over the past 50 Myr (panel b), and LyC photon escape fraction assuming ionization-bounded conditions (i.e. 'picket-fence' models; panel c). In both panels, the EWs are normalized by their values adopting a 50 Myr CSFH model with 0.2 Z⊙, log U = -2.5, fesc = 0, and τV = 0 with the Gutkin, Charlot & Bruzual (2016) models implemented in BEAGLE.</figcaption></figure>

# 5 IMPLICATIONS FOR STAR FORMATION AND IOINIZING EFFICIENCY IN THE REIONIZATION ERA

In the previous section, we demonstrated that the JADES/NIRCam data indicate a strong decline in [O III] + H β EW with UV luminosity at z ~ 6 - 9, yet an approximately constant H α EW distribution across an order of magnitude in L$_{UV}$ at z ~ 6 (Fig. 8). Here, we begin by discussing what physically may be driving the different L$_{UV}$ trends for [O III] + H β and H α EWs then explore the results implied from the SED fits, with the ultimate goal of considering potential implications for the nature of very UV-faint reionization-era galaxies.

## 5.1 The role of metallicity in explaining the EW trends

One factor likely influencing the UV luminosity trends in nebular EWs is metallicity. Due to the deficit of oxygen atoms at very low metallicity, the strength of nebular [O III] + H β emission declines substantially at ≲0.2 Z⊙ (see Fig. 11a). Because fainter, lower mass galaxies are generally expected to be less chemically enriched, this likely contributes to much weaker [O III] + H β EWs in our faintest subset. We consider how much lower the typical metallicity of very faint (⟨M$_{UV}$}⟩ ≈ -17.5) z ~ 6 - 9 galaxies would have to be relative to the bright subset (⟨M$_{UV}$}⟩ ≈ -20) to completely explain the ≈ 2.4× decline in median [O III] + H β EW inferred in Section 4.2. For this test, we utilize the Gutkin et al. (2016) models employed in BEAGLE which we note assume equal stellar and gas-phase metallicities (i.e. no α-enhancement). From the Gutkin et al. (2016) models, we find that the metallicity would need to drop by ≈1 dex to solely explain the substantial decline in [O III] + H β EW (see Fig. 11a).

Such a dramatic (≈1 dex) decline in metallicity would have a significant impact on the H α EWs. Due to changes in stellar opacity, models with lower metallicity lead to more efficient ionizing photon production and hence larger H α EWs. From the Gutkin et al. (2016) star-forming models employed in our BEAGLE fits, we expect that lowering the metallicity by 1 dex from 0.2 to 0.02 Z⊙ would result in a ≈50 per cent increase in H α EWs. Therefore, while lower metallicities at fainter UV luminosities may partially explain the decreasing [O III] + H β EWs, something else must also be changing

between the bright and very faint populations to explain why the H α EW distribution remains nearly constant at z ~ 6.

## 5.2 Evidence for bursty star formation histories at z ≳ 6

One potential way to offset the expected increase in H α EW at fainter L$_{UV}$ due to lower metallicity is by invoking differences in the recent SFHs of bright versus very faint reionization-era galaxies (we discuss an alternative solution with differing ionizing photon escape fractions below<fnref n="6" />). Because the Balmer and [O III] nebular lines are primarily powered by hot, massive O stars that have formed in the most recent ≈3 Myr, a strong drop in the SFR over this time-scale results in weaker emission lines, decreasing the EWs of H α and [O III] + H β (see Fig. 11b). If UV-faint galaxies are more often caught in a recent downturn of star formation (relative to the UV-bright population), this may explain why their H α EWs remain comparable to that of the more UV luminous population. Such a change in recent SFHs would also contribute to the reduced [O III] + H β EWs in the very UV-faint population (Fig. 8), working in tandem with the shift toward lower metallicities.

To better assess whether the photometric data support systematic changes in metallicity and recent SFHs across the wide UV luminosity range spanned by our sample (-22 ≲ M$_{UV}$ ≲ -16.5), we utilize the results of our TcSFH BEAGLE fits. For each galaxy, we compute the inferred ratio of the SFR averaged over the past 3 Myr (SFR$_{3 Myr}$) to that over the past 50 Myr (SFR$_{50 Myr}$) as a proxy for whether their SFH has recently declined or risen on time-scales most relevant for nebular line emission. The 50 Myr time-scale adopted in the denominator of this ratio was chosen to reflect the typical light-weighted age of the sample, though we note that our findings below are qualitatively unchanged if we instead adopt a 30 or 100 Myr time-scale. We quantify the distributions of SFR$_{3 Myr}$/SFR$_{50 Myr}$ and metallicity in the same three UV luminosity bins as considered for the EW distribution analysis (see Section 4.1), again adopting a lognormal distribution for each M$_{UV}$ bin and applying equation (2). Here, we focus our analysis solely on the z ~ 6 sample given the lack of H α information for the large majority of the z ~ 7 - 9 sample.

<footnote n="6">⁶Other mechanisms may also be responsible for decreasing the [O III] and Balmer line EWs at fainter UV luminosities, including changes in IMF, contributions from low-luminosity AGN, or α-enhanced metallicities.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Illustration of how the strong decline in [O III] + H β EW with UV luminosity yet approximately constant H α EW distributions at z ~ 6 (see Fig. 8) can be explained by either the combination of lower metallicity and more recently declining SFHs in the very UV-faint population (top) or lower metallicities and higher LyC escape fractions in the very UV-faint population (bottom). These two scenarios have different implications for contributions of galaxies along the UV luminosity function to cosmic reionization, as also illustrated in the inferred ionizing photon production efficiency distributions in the right panels. Deep spectra will be required to better distinguish the mechanisms dominating the UV-luminosity-dependent EW trends found here.</figcaption></figure>

Because the SFR$_3$ Myr/SFR$_{50}$ Myr ratio is by definition restricted to values $\le \frac{50}{3}$ (this maximum value occurs when no stars are assumed to form between 3–50 Myr ago), we truncate the lognormal distributions at this upper limit (i.e. $P(\text{SFR}_3\text{Myr}/\text{SFR}_{50\text{Myr}} > \frac{50}{3} \| \theta) = 0$). Similarly, the lognormal distributions for metallicity are truncated outside the fitted parameter space in our SED fits ($-2.2 \le \log(Z/Z_\odot) \le -0.3$). When quoting median inferred values for metallicity and SFR$_3$ Myr/SFR$_{50}$ Myr below we account for this truncation, though the standard deviations are quoted directly from the assumed lognormal functional form.

The TcSFH BEAGLE outputs suggest a significant anticorrelation between UV luminosity and metallicity at $z \sim 6$ (Fig. 12a), with a median inferred metallicity of $0.27^{+0.13}_{-0.08}$, $0.11^{+0.01}_{-0.01}$, and $0.06^{+0.01}_{-0.01} Z_\odot$ in the bright ($\langle M_{UV} \rangle = -20.0$), faint ($\langle M_{UV} \rangle = -18.7$), and very faint ($\langle M_{UV} \rangle = -17.5$) subsets, respectively. Therefore, in context of these models, changes in metallicity contribute significantly to the different UV luminosity trends with [O III] + H β and H α EWs. We check whether these inferred metallicities are consistent with predictions from simulations under very simple assumptions on the connection between $M_{UV}$ and stellar mass and assuming that our inferred metallicities largely reflect the gas-phase metallicity.

In our $z \sim 6$ sample, the typical stellar mass at $M_{UV} \approx -20$ is $M_* \sim 4 \times 10^8 M_\odot$ while it is $\sim 2 \times 10^7 M_\odot$ at $M_{UV} \approx -17.5$ (see Fig. 2). The inferred metallicities above from the TcSFH models are reasonably consistent (factor of $\lesssim 2$ difference) with the predicted $z \sim 6$ gas-phase metallicities at the respective stellar masses from the ILLUSTRIS TNG (Torrey et al. 2019), FIRSTLIGHT (Langan, Ceverino & Finlator 2020), and ASTREAUS (Ucci et al. 2023) simulations, though systematically $\approx 5 \times$ higher than that predicted from the FIRE simulations (Ma et al. 2016).

We now investigate whether the TcSFH BEAGLE fit outputs also imply differences in the recent SFHs of bright versus very faint $z \sim 6$ galaxies. Using equation (2), we infer a median SFR$_3$ Myr/SFR$_{50}$ Myr ratio of $1.5^{+0.3}_{-0.2}$, $0.8^{+0.1}_{-0.1}$, and $0.6^{+0.9}_{-0.2}$ in the bright, faint, and very faint subsets, respectively. Therefore, under context of the adopted SED models, we do find evidence that brighter $z \sim 6$ galaxies tend to have more very-recently rising SFHs than the faintest sources (Fig. 12b). Notably, not all bright (very faint) galaxies are inferred to have very-recently rising (declining) SFHs as we recover substantial scatter in the SFR$_3$ Myr/SFR$_{50}$ Myr ratio at fixed UV luminosity. The standard deviation of the assumed lognormal distribution is inferred to be $0.52^{+0.09}_{-0.07}$, $0.62^{+0.08}_{-0.06}$, and $1.21^{+0.66}_{-0.33}$ dex in the bright, faint, and
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Example inferred SFHs of two $z \sim 6$ extreme emission line galaxies (top; SEDs shown in Fig. 7) and two $z \sim 6$ galaxies showing SEDs consistent with young light-weighted ages yet weak line emission (bottom; SEDs shown in Fig. 14). The median SFHs (solid lines) along with their inner 68 per cent credible intervals (grey regions) are taken from the BEAGLE TcSFH models. The extreme emission line galaxies are experiencing a rapid upturn in SFR, while the young, weak line emitters are inferred to have recently experienced a strong downturn in SFR, consistent with bursty SFHs among reionization-era galaxies.</figcaption></figure>

very faint subsets, respectively, implying generally more variation in recent SFHs at fainter luminosities (and hence lower stellar masses).

There are two particular subsets of our sample that provide evidence in favour of bursty SFHs among reionization-era galaxies. The first subset consists of galaxies showing extremely strong nebular line emission. The very high [O III] + Hβ EWs (~2000–5000 Å) confidently inferred among ≈50 galaxies in our sample (see Section 3.2.2) imply that a considerable fraction of Lyman-break $z \sim 6–9$ galaxies have recently experienced a rapid strong upturn in SFR, yielding SEDs that are completely dominated by hot, massive O stars. We plot the inferred SFHs of two of the confident $z \sim 6$ extreme emission line galaxies in the top panels of Fig. 13, where both are inferred to have substantially (≥5×) higher SFRs over the most recent 3 Myr relative to that over the past 10–50 Myr (see also e.g. Tacchella et al. 2023b). Such a dramatic recent increase in SFR is consistent with a scenario in which these extreme emission line galaxies are undergoing a burst of star formation. We note that a variety of recent papers are building empirical evidence for bursty SFHs among $z \ge 6$ galaxies using different methods (Strait et al. 2023; Dressler et al. 2023, 2024; Looser et al. 2023, 2024; Dome et al. 2024). Bursty star formation in the early Universe is also predicted in several simulations and models of galactic evolution (e.g. Kimm et al. 2015; Ceverino, Klessen & Glover 2018; Faucher-Giguère 2018; Ma et al. 2018; Furlanetto & Mirocha 2022).

Another subset of galaxies in our JADES sample may be caught during a relatively inactive phase of star formation. We identify $z \sim 6$

galaxies exhibiting NIRCam SEDs consistent with relatively weak Balmer breaks (thus implying young light-weighted ages) yet also weak [O III] + Hβ and Hα emission. The SEDs of two of these galaxies are shown in Fig. 14 for illustrative purposes (we discuss the abundance of similar objects within our full $z \sim 6$ JADES sample below). In the left panels of Fig. 14, we show the fitted SEDs from the BEAGLE CSFH models which clearly struggle to reproduce the observed photometry. For both galaxies, the CSFH models are forced to extremely low-metallicity solutions (≈0.01 Z⊙) to simultaneously explain the weak Balmer breaks[^7] ($F_v$(4200 Å) / $F_v$(3500 Å) ≈ 1 as implied by the fairly flat F200W–F277W colours) as well as the low [O III] + Hβ EWs (≈300–400 Å as implied by the weak photometric excesses in F335M and F356W). The extremely low-metallicity solutions in turn cause the models to drastically overpredict the Hα EWs relative to that implied by the observed F410M–F444W colours, suggesting that the SEDs of these galaxies cannot be well explained by our fiducial CSFH models.

When allowing for more flexible SFHs, we obtain SED models that yield much better fits to the photometry of these two illustrative $z \sim 6$ galaxies. As shown in the middle panels of Fig. 14, the BEAGLE TcSFH models can reproduce the NIRCam photometry in all nine bands with moderately low-metallicity solutions (≈0.05–0.1 Z⊙) for these two $z \sim 6$ galaxies. In these TcSFH solutions, the two galaxies are currently experiencing a strong decline in SFR that has followed a burst of star formation which occurred ~5–50 Myr ago (see bottom panels of Fig. 13). Such SFHs yield relatively few O stars that are still alive to power the [O III] and Balmer lines, while the B and A stars that formed ≥10 Myr ago continue to dominate the rest UV + optical SEDs. It is thus possible that these galaxies were in an extreme emission line phase just ~10–30 Myr earlier, and are now caught in a relatively passive phase of star formation yielding weak nebular lines yet also young SEDs from the continua. We note that a significant population of reionization-era galaxies with young light-weighted ages and weak [O III] + Hβ emission was first tentatively identified in Endsley et al. (2023) from early JWST data, though the lack of Hα constraints for the $z \sim 6.5–8$ sample considered therein left much uncertainty on the physical origin of the weak line emission (see also Prieto-Lyon et al. 2023; Simmonds et al. 2023). With the improved imaging sensitivity and SED sampling of JADES and GOODS, we have now extended this investigation to a $z \sim 6$ sample and find that our fiducial CSFH models cannot adequately explain the photometry of at least some of these very early galaxies.

In addition to the two illustrative objects shown in Fig. 14, we identify a few other $z \sim 6$ galaxies that confidently show SEDs consistent with young light-weighted ages and weak line emission, perhaps indicating that they have recently experienced a strong downturn in SFR. To formally quantify the number of such candidate objects in our JADES sample, we begin by selecting $z \sim 6$ galaxies with >84 per cent probability of having [O III] + Hβ and Hα EWs of <800 Å from the TcSFH posterior outputs. We then apply the colour-colour cut F200W–F277W<F150W–F200W+0.3 to ensure that the selected galaxies do not show a strong Balmer break. While we could apply cuts on the inferred Balmer break strength

<footnote>[^7]: While the CSFH model fits to these two galaxies yield fairly high light-weighted ages (~100 Myr), the resulting Balmer breaks in the models remain weak due to strong nebular continuum emission. As has been discussed elsewhere (e.g. Byler et al. 2017; Topping et al. 2022b), nebular continuum emission is stronger at lower metallicities and the fractional contribution of nebular emission to the total (stellar + nebular) continuum is larger immediately blueward of the Balmer break than redward of this break.</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Illustration of two $z_{phot} \approx 6.0$ galaxies (one per row) showing SEDs consistent with young light-weighted ages as well as weak [O III] + H β and H α emission. We find that our fiducial CSFH models (with $f_{esc} = 0$; left panels) struggle to match the data while models allowing for strong recent drops in SFR (middle panels; see Fig. 13 for inferred SFHs) or very high LyC escape fractions (right panels) provide considerably better fits to the data.</figcaption></figure>

from the TcSFH fits directly, we find that such an approach is not ideal given that the Balmer break strength can increase dramatically on short (~10 Myr) time-scales when allowing for strong drops in SFR. Finally, because our goal is to identify galaxies that are relatively poorly fit with the CSFH models yet considerably better fit by the TcSFH models, we only include those with a best-fitting $\chi^2_{CSFH} > 14$ (we are fitting 14 data photometric data points) and $\chi^2_{CSFH} - \chi^2_{TcSFH} > 10$. This rather strict selection results in a total of seven $z \sim 6$ galaxies, including the two illustrative examples shown in Fig. 14. The majority of these seven galaxies have inferred $SFR_3_{Myr}/SFR_{50_{Myr}}$ ratios of < 0.2 with >84 per cent confidence from the TcSFH posterior outputs.

We emphasize that the seven objects described above are not an exhaustive census of $z \sim 6$ galaxies in our sample showing potential evidence of a recent strong downturn in SFR. They are simply the most confident cases and hence are limited to galaxies with the highest signal-to-noise photometry [$S/N(F277W) \gtrsim 30$ before adding the 5 per cent systematic error]. The inferred $SFR_3_{Myr}/SFR_{50_{Myr}}$ ratio distributions shown in Fig. 12 demonstrate the collective evidence in our sample for a significant population of $z \sim 6$ galaxies (particularly at the UV-faint end) with recently declining SFHs. This statistical, photometric evidence for strongly declining SFHs among some reionization-era galaxies is further supported by the spectroscopic confirmation of a $z = 7.3$ ('mini-)quenched' galaxy from recent JWST observations (Dome et al. 2024; Looser et al. 2024; for similar examples at $z \sim 5$ see also Strait et al. 2023; Looser et al. 2023).

### 5.2.1 The different nature of UV-bright and UV-faint reionization-era galaxies

The TcSFH models potentially provide new insight into the nature of reionization-era galaxies along the UV luminosity function. In

context of these models, galaxies at the bright end are weighted more towards systems that have very recently experienced an increase in star formation rate ($SFR_3_{Myr}/SFR_{50_{Myr}} > 1$). Objects with such SFHs will have a greater fraction of emergent light arising from hot and exceptionally luminous O stars, resulting in a temporary boost to their light-to-mass ratio. Therefore, UV-luminous reionization-era galaxies may more often consist of relatively low-mass objects that have been upscattered in luminosity due to their SFH, particularly when considering the exponential decline of the mass function at high masses. However, given the large systematic uncertainties in stellar mass estimates among $z \sim 6-9$ galaxies (see Section 3), drawing such a conclusion directly from the data is quite challenging. If early galaxies do indeed have very bursty SFHs, this would help explain the apparent abundance of bright $z > 10$ galaxies recently identified by JWST as one would then expect a higher fraction of low-mass galaxies to be upscattered above the detection threshold (e.g. Mason, Trenti & Treu 2023; Mirocha & Furlanetto 2023; Sun et al. 2023).

In contrast, galaxies much further down the UV luminosity function would be less weighted towards systems with recent upturns in their SFR in context of the TcSFH models. None the less, the substantial scatter in the $SFR_3_{Myr}/SFR_{50_{Myr}}$ ratio inferred among our very UV-faint $z \sim 6$ subset ($\approx 1$ dex) does imply considerable diversity in the recent SFHs of this population. This is directly demonstrated by the presence of galaxies with extremely strong [O III] + H β EWs ($\gtrsim 1500$ Å) within the very UV-faint bin (e.g. JADES-GS + 53.17835 – 27.77879 in Fig. 7 and the $M_{UV} = -16.7$ galaxy at $z_{spec} = 7.28$ described in Saxena et al. 2023) indicating a strong rise in their SFHs on ~3 Myr time-scales. It may be the case that galaxies considerably fainter than the exponential cut-off of the UV luminosity (and mass) function are more equally likely to be a relatively low-mass object experiencing a recent strong rise in SFR (resulting in high $L_{UV}/M_*$) as they are to be
a relatively high-mass object experiencing a recent strong decline in SFR (yielding relatively low $L_{UV}/M_*$), thereby yielding the broader SFR$_3$ Myr/SFR$_{50}$ Myr distribution at fainter $M_{UV}$.

## 5.3 An alternative explanation: efficient LyC escape

Above, we have demonstrated how bursty SFHs (when coupled with lower metallicities at fainter $M_{UV}$) can explain the $z \sim 6$ [O III] + H β and H α EW trends with UV luminosity found in Sections 4.2 and 4.3. In this subsection, we discuss how these EW trends can alternatively be explained by allowing for a large fraction of ionizing photons to escape the reionization-era galaxies in our sample. Because the [O III] and Balmer nebular lines are powered by ionizing photons interacting with dense gas in galaxies, removing a substantial fraction of these ionizing photons will result in fewer recombinations and hence weaker lines. We estimate the LyC ($\lambda_{rest} < 912$ Å) photon escape fractions required to explain the photometric data (and resulting [O III] + H β and H α EW trends) under the assumption of a CSFH by using the 'picket-fence' model within BEAGLE. In these picket fence models (e.g. Conselice et al. 2000; Heckman et al. 2011; Rivera-Thorsen et al. 2017), a given fraction (hereafter $f_{esc}$) of LyC photons are assumed to escape the H II regions without ionizing the gas while still assuming ionization-bounded conditions. Therefore, all nebular line and continuum emission (independent of wavelength) are reduced by a factor of $1 - f_{esc}$ relative to the base CSFH BEAGLE models (see Fig. 11c).

In the picket-fence BEAGLE SED fits, we adopt a log-uniform prior on $f_{esc}$ between 0.001 and 0.8 with all other priors equivalent to the nominal CSFH fits. We then infer lognormal distributions for both metallicity and $f_{esc}$ using equation (2) in the three different UV luminosity bins for our $z \sim 6$ sample. While density-bounded photoionization models may be more appropriate at the highest $f_{esc}$ values ($\gtrsim 0.5$), our goal here is to simply gain an initial sense of the escape fractions required to explain the photometric data (and the implied EW distributions) under the assumption of CSFHs. Adopting density-bounded models (Plat et al. 2019) would generally result in a slower decline of the [O III] EWs with increasing $f_{esc}$ relative to the Balmer line EWs given that the high-ionization zones are the last to be significantly impacted by ionizing photon escape in this context. None the less, the difference in EWs between the ionization-bounded and density-bounded BEAGLE models is $\lesssim 10$ per cent at $f_{esc} < 0.8$.

The picket-fence CSFH BEAGLE fits also yield an anticorrelation between UV luminosity and metallicity at $z \sim 6$ (see Fig. 12d) as physically expected, with median values of $0.35^{+0.09}_{-0.09}$, $0.38^{+0.08}_{-0.10}$, and $0.05^{+0.03}_{-0.02}Z_{\odot}$ among the bright ($\langle M_{UV} \rangle = -20.0$), faint ($\langle M_{UV} \rangle = -18.7$), and very faint ($\langle M_{UV} \rangle = -17.5$) $z \sim 6$ subsets, respectively. However, the inferred distribution of metallicities in a given UV luminosity bin is considerably different in shape from that inferred from the TcSFH BEAGLE fits adopting $f_{esc} = 0$ (cf. panel (a) in Fig. 12). With the picket-fence models, both the bright and faint subsets have inferred metallicity distributions that peak at fairly high values ($\gtrsim 0.3 Z_{\odot}$) while the metallicity distribution for the very faint bin is exceptionally broad with an inferred standard deviation of $>0.8$ dex resulting in a fairly uniform distribution across the $-2.2 \le \log(Z/Z_{\odot}) \le -0.3$ parameter space allowed in our SED fits. It is unclear whether such a distribution in $z \sim 6$ galaxy metallicity is physical at such faint UV luminosities and low stellar masses (maximum $M_* = 10^9 M_{\odot}$; see Fig. 2). None the less, we proceed with investigating what correlation between UV luminosity and $f_{esc}$ is required in context of these models to explain the photometry of the $z \sim 6$ sample.

Under the assumption of CSFHs, the trends in [O III] + H β and H α EW found in Sections 4.2 and 4.3 imply a strong increase in $f_{esc}$ towards very faint UV luminosities at $z \sim 6$ (Fig. 12d). With the picket-fence BEAGLE SED models, we infer median $f_{esc}$ values of $0.07^{+0.05}_{-0.03}$, $0.21^{+0.04}_{-0.04}$, and $0.35^{+0.05}_{-0.05}$ in the bright, faint, and very faint subsets, respectively, when adopting lognormal distributions for $f_{esc}$ in each luminosity bin. Moreover, the fraction of galaxies with substantial ionizing photon leakage ($f_{esc} > 20$ per cent) is inferred to be $28^{+8}_{-6}$ per cent, $52^{+6}_{-6}$ per cent, $76^{+9}_{-10}$ per cent in each of the respective $M_{UV}$ bins. As can be seen in Fig. 12, the inferred distribution of $f_{esc}$ simultaneously narrows and peaks closer to the allowed upper bound of $f_{esc} = 0.8$ towards fainter UV luminosity. The generally bluer UV slopes, smaller rest-UV sizes, and lower inferred dust attenuation among fainter $z \gtrsim 6$ galaxies (e.g. Bouwens et al. 2014; Shibuya et al. 2015; Bhatawdekar & Conselice 2021) do imply that these galaxies would typically be more efficient leakers of LyC photons based on local constraints (e.g. Chisholm et al. 2022; Flury et al. 2022). However, it is unclear whether a dramatic (factor of $\sim 5$) change in median $f_{esc}$ between $M_{UV} \approx -20$ and $\approx -17.5$ at $z \sim 6$ would be expected. No evidence yet exists among local samples that the LyC escape fraction strongly correlates with far-UV luminosity nor stellar mass (Izotov et al. 2021; Flury et al. 2022).

We consider whether the picket-fence CSFH models can explain the photometry of the seven $z \sim 6$ galaxies confidently showing SEDs consistent with young light-weighted ages as well as relatively weak nebular lines discussed in Section 5.2. Five of these seven galaxies are considerably better fit with the picket-fence models relative to standard CSFH models ($\Delta\chi^2 > 10$) and, in each of these five cases, the inferred LyC escape fraction is $f_{esc} > 50$ per cent with $>84$ per cent probability from the output posterior. We show the picket-fence CSFH model SEDs for the two illustrative example galaxies in the rightmost panels of Fig. 14. For one of these objects, the best-fitting picket-fence CSFH model has $\chi^2 = 19.5$ across the 14 fitted data points, implying a somewhat poorer fit relative to the best-fitting TcSFH model ($\chi^2 = 9.8$). However, in four of the seven galaxies under consideration here, there is little-to-no preference in the data for the TcSFH models over the picket-fence CSFH models ($\Delta\chi^2 < 5$). As discussed in the following subsection, these two different scenarios (bursty SFH or substantial LyC escape) have very different implications for reionization, and deep spectroscopic observations will ultimately be required to better determine the physical origin of very early galaxies showing young light-weighted ages and relatively weak line emission.

## 5.4 Potential implications for reionization

Above, we have demonstrated that the $z \sim 6$ [O III] + H β and H α EW trends found in Sections 4.2 and 4.3 can be explained by a combination of decreasing metallicity with UV luminosity, in addition to either (1) systematically more very-recently rising SFHs among brighter systems with small LyC leakage at all luminosities or (2) systematically much greater LyC leakage among fainter systems with CSFHs at all luminosities (Fig. 12). To better understand the potential implications for the ionizing efficiency of bright vs. very faint galaxies in these two limiting cases of model assumptions (i.e. bursty SFH with $f_{esc} \approx 0$ or CSFH with $0 < f_{esc} < 0.8$), we infer the distribution of ionizing photon production efficiencies ($\xi_{ion}$) as a function of UV magnitude in each case. Here, we are considering the $\xi_{ion}$ defined as the production rate of LyC photons divided by the observed luminosity at rest-frame 1500 Å (i.e. no correction is made to $L_{1500}$ for dust attenuation or nebular continuum emission).
With our TcSFH BEAGLE models that adopt $f_{\text{esc}} = 0$, we infer median ionizing photon production efficiencies of $\log[\xi_{\text{ion}}/(\text{erg}^{-1} \text{ Hz})] = 25.65^{+0.03}_{-0.03}$, $25.48^{+0.03}_{-0.02}$, and $25.35^{+0.05}_{-0.05}$ in the bright ($\langle M_{\text{UV}} \rangle = -20.0$), faint ($\langle M_{\text{UV}} \rangle = -18.7$), and very faint ($\langle M_{\text{UV}} \rangle = -17.5$) $z \sim 6$ subsets, respectively (Fig. 12c). This is consistent with expectations given that the very faint population is inferred to generally have more very-recently declining SFHs with these models (Fig. 12b), resulting in a lower fraction of O stars producing the bulk of ionizing photons. We note that the typical $\xi_{\text{ion}}$ we infer among our faintest subset is in excellent agreement with the value reported in Atek et al. (2024) for a sample of four $z_{\text{spec}} = 6.2-6.8$ galaxies with $-17.2 \le M_{\text{UV}} \le -16.5$ ($\log[\xi_{\text{ion}}/(\text{erg}^{-1} \text{ Hz})] = 25.4 \pm 0.2$). The standard deviation of the assumed lognormal distribution in $\xi_{\text{ion}}$ is inferred to moderately increase from $\approx 0.20$ dex in the bright and faint bins to $\approx 0.3$ dex in the very faint bin, also consistent with the inferred larger variation in recent SFHs among the very faint population. This broadening reflects that, while the median $\xi_{\text{ion}}$ decreases with UV luminosity, there remains a significant subset of a very UV-faint $z \sim 6$ galaxies with very efficient production of hydrogen ionizing photons ($\xi_{\text{ion}} \gtrsim 10^{25.7} \text{ erg Hz}^{-1}$; see Fig. 7). Similar sources at extremely low UV luminosities ($M_{\text{UV}} \sim -16$) have now been spectroscopically confirmed (Atek et al. 2024).

The UV luminosity dependence on $\xi_{\text{ion}}$ is inferred to be much weaker with the CSFH models that allow for significant LyC leakage (Fig. 12f). From these fits, the median inferred $\log[\xi_{\text{ion}}/(\text{erg}^{-1} \text{ Hz})]$ is $25.75^{+0.03}_{-0.02}$, $25.64^{+0.02}_{-0.01}$, and $25.60^{+0.03}_{-0.02}$ in the respective UV luminosity bins, with the standard deviation in the range $\approx 0.15-0.18$ dex in all bins. Such a small UV luminosity dependence on $\xi_{\text{ion}}$ yet large dependence on $[\text{O III}] + \text{H}\beta$ is accommodated in the models by pushing $f_{\text{esc}}$ to the very large values ($\gtrsim 0.5$) in the very faint population. We also note that the CSFH models allow for less dynamic range in the inferred $\xi_{\text{ion}}$ value of an individual galaxy given that a significant population of O stars is always assumed to be present (given the allowed age of a galaxy at $z \gtrsim 6$), in contrast to scenarios with the TcSFH models where the SFR can sharply decline in the most recent $\sim 3$ Myr (see Fig. 13). These differences in model assumptions contribute to the systematically higher and narrower inferred $\xi_{\text{ion}}$ distributions with the CSFH models.

The outputs of the TcSFH and picket-fence CSFH models therefore have very different implications for the contribution of galaxies to reionization along the UV luminosity function. The picket-fence CSFH models yield a substantial, systematic increase in $f_{\text{esc}}$ towards very faint UV luminosities yet a weak UV-luminosity dependence on $\xi_{\text{ion}}$, implying that the faintest galaxies vastly dominate the ionizing photon budget for reionization. However, the TcSFH models assume that all galaxies regardless of $M_{\text{UV}}$ have small $f_{\text{esc}}$, but that the most UV-luminous galaxies have systematically larger $\xi_{\text{ion}}$. While the assumption that all reionization-era galaxies have small $f_{\text{esc}}$ may very well be invalid, so too may be the assumption that all such galaxies have a CSFH given the abundance of extreme line emitters. Our primary intent here is to demonstrate that existing uncertainties on what is physically driving the EW trends with $M_{\text{UV}}$ directly result in large uncertainties on the relative contribution of galaxies along the UV luminosity function to reionization.

In summary, the UV-luminosity dependence on $[\text{O III}] + \text{H}\beta$ and $\text{H}\alpha$ EWs implied by the NIRCam data provide empirical evidence on a statistical basis for distinct differences in the nature of UV-bright and UV-faint $z \sim 6-9$ galaxies. We have demonstrated that these trends can be explained by either lower metallicity and more recently declining SFHs in the very UV-faint population, or lower metallicities and higher LyC escape fractions in the very UV-faint population. The degeneracy of bursty SFHs and LyC leakage on the nebular

lines cannot be broken with photometry alone and it may very well be the case that both factors play a significant, perhaps connected, role (e.g. Fletcher et al. 2019; Tang et al. 2019; Barrow et al. 2020; Ma et al. 2020). As noted previously in this sub-section, additional mechanisms such as different IMFs, low-luminosity AGN, or $\alpha$-enhanced metallicities could also play a significant role in driving the EW distribution trends found in Sections 4.2 and 4.3. These various scenarios will have different implications for the contributions of galaxies along the UV luminosity function to cosmic reionization, and thus deep spectroscopic efforts are critical to better determine the physical mechanisms responsible for the $[\text{O III}] + \text{H}\beta$ and $\text{H}\alpha$ EW trends.

# 6 OVERDENSITIES AROUND LY$\alpha$ EMITTERS AT $z > 7$

It has long been proposed that galaxies exhibiting strong Ly$\alpha$ emission (EW $> 25$ Å) at $z > 7$ often trace strong galaxy overdensities that are capable of generating large ionized bubbles prior to the completion of reionization (e.g. Wyithe & Loeb 2005; Castellano et al. 2016; Hutter et al. 2017; Stark et al. 2017; Weinberger et al. 2018; Jung et al. 2020; Tilvi et al. 2020; Endsley & Stark 2022; Larson et al. 2022; Leonova et al. 2022; Lu et al. 2024). In this scenario, strong $z > 7$ LAEs would commonly act as signposts of the earliest sites of structure formation. Here, we use the improved photometric redshifts enabled by the JADES/NIRCam GOODS-N and GOODS-S imaging to better test the connection between overdensities and strong $z > 7$ Ly$\alpha$ emission.

We search the literature for $z > 7$ galaxies with known high-confidence Ly$\alpha$ detections in the GOODS fields (i.e. S/N $> 10$ in Ly$\alpha$ alone or S/N $> 10$ of a systemic line plus a $> 7\sigma$ Ly$\alpha$ detection), resulting in three objects: JADES-GS-z7-LA at $z_{\text{Ly}\alpha} = 7.281$ (Saxena et al. 2023), z8_GND_5296 at $z_{\text{Ly}\alpha} = 7.508$ (Finkelstein et al. 2013; Tilvi et al. 2016; Jung et al. 2019, 2020), and z7_GND_16863 at $z_{\text{Ly}\alpha} = 7.599$ (Jung et al. 2019, 2020). Two of these systems fall within the JADES/NIRCam footprint (z8_GND_5296 and JADES-GS-z7-LA) and we investigate whether either show a significant photometric overdensity. The galaxy z8_GND_5296 is a UV-bright ($M_{\text{UV}} = -21.5$) system in the GOODS-N field with a most recent Ly$\alpha$ EW measurement of $33 \pm 4$ Å from Jung et al. (2020). JADES-GS-z7-LA is, on the other hand, a very UV-faint ($M_{\text{UV}} = -16.7$) galaxy in the GOODS-S field with a near-unity Ly$\alpha$ escape fraction yielding an extremely large EW of $400 \pm 90$ Å implying that it lies within a very large ionized bubble (Saxena et al. 2023). Both of these strong LAEs exhibit very high $[\text{O III}] + \text{H}\beta$ EWs of $\approx 1500$ Å placing them in the top $\approx 17$ per cent and $\approx 7$ per cent of the $[\text{O III}] + \text{H}\beta$ EW distribution, respectively, given their UV luminosities (see Section 4.2).

Here, we focus on quantifying the photometric overdensities (following e.g. Castellano et al. 2016; Endsley et al. 2021b; Leonova et al. 2022; Whitler et al. 2024) of the two $z > 7$ LAEs in JADES noted above given that deep spectroscopic follow-up of $z \sim 7-8$ candidates surrounding these rare sources remains incomplete. We thus aim to supplement our base $z \sim 7-9$ selection criteria with additional cuts to better identify neighbouring galaxies lying relatively close in redshift space to the two strong $z > 7$ LAEs in JADES. Fortunately, both of these LAEs lie at $z = 7.0-7.6$ where relatively high-precision photometric redshifts can be obtained by exploiting the medium-band photometry, thereby improving our photometric overdensity estimates. At $z \approx 7.0-7.6$, the $[\text{O III}]\lambda\lambda4959,5007$ doublet falls in the F410M filter yet outside of the F356W bandpass and thus galaxies in this redshift range with large $[\text{O III}] + \text{H}\beta$
EWs will show remarkably red F356W–F410M colours. From the inferred [O III] + H β EW distributions in the JADES dataset (see Table 2 and Fig. 8), we expect the majority of $M_{UV} < -18$ galaxies at $z \approx 7.3$ will exhibit [O III] + H β EWs of >400 Å which, assuming a flat rest-optical continuum (in $F_v$), will result in F356W–F410M colours of $\gtrsim 0.6$ mag at $z \approx 7.0–7.6$. Such prominent colours should be relatively easy to identify for all but the faintest objects in our deep JADES imaging. Moreover, at $z \approx 7.0–7.6$, the F335M band is entirely redward of the Balmer break at 3648 Å and thus we would generally expect flat F335M–F356W colours given the typically little dust attenuation inferred for the sample ($A_V \lesssim 0.01$ mag).

Guided by the above expectations for the colours of $z \approx 7.0–7.6$ galaxies, we identify potential neighbours of the very UV-faint ($M_{UV} = -16.7$) $z = 7.28$ LAE JADES-GS-z7-LA and the UV-bright ($M_{UV} = -21.5$) $z = 7.51$ LAE z8_GND_5296 by supplementing our F090W dropout selection criteria (Section 2.3) with the below additional cuts:

(i) F356W–F410M > 0.6;

(ii) F335M – F356W < 0.3;

(iii) F150W < 29.0.

Here, we have added the condition of $F150W < 29$ (corresponding to $M_{UV} \lesssim -18.0$) to ensure that our selection of candidate neighbours is fairly complete at $z \approx 7.0–7.6$ given the different depth tiers of the JADES imaging (Eisenstein et al. 2023). The on-sky positions of the candidate $z \approx 7.0–7.6$ galaxies in the JADES GOODS-S and GOODS-N footprints are shown in Fig. 15, along with the positions of the two strong $z > 7$ LAEs covered by JADES. In the GOODS-S field containing the very UV-faint LAE, we identify 43 potential $z \approx 7.0–7.6$ galaxies with $F150W < 29$ across the 42.2 arcmin² JADES area with F335M coverage. Across the 47.2 arcmin² JADES/GOODS-N area with F335M coverage, we find 42 $z \approx 7.0–7.6$ galaxies which includes the UV-bright ($M_{UV} = -21.5$) LAE that is the source of interest in this field (z8_GND_5296). The very similar numbers of $z \approx 7.0–7.6$ candidates in each field despite significantly shallower NIRCam exposure times over much of GOODS-N (Eisenstein et al. 2023) is consistent with past findings that the (completeness-corrected) number density of $z \sim 7$ Lyman-break galaxies is considerably higher over GOODS-N than GOODS-S (e.g. Bouwens et al. 2015; Finkelstein et al. 2015). Below, we first discuss the surface density of potential $z \approx 7.0–7.6$ neighbours surrounding the LAEs selected from JADES, then proceed to calculate the surrounding photometric overdensities folding in the expected completeness of our selection.

It is immediately clear from Fig. 15 that the surface density of $z \approx 7.0–7.6$ galaxies is unusually high near the GOODS-S LAE JADES-GS-z7-LA. In particular, about 1 arcmin east of this very UV-faint LAE lies a complex of six distinct $z \approx 7.0–7.6$ galaxies with a maximum separation of 23 arcsec from one another. The surface density of $z \approx 7.0–7.6$ galaxies in this very small (0.073 arcmin²) region of GOODS-S is 87.3 galaxies/arcmin⁻², a factor of ≈70× higher than that on average over the deep imaging region covering JADES-GS-z7-LA. The six objects comprising this dense concentration of $z \approx 7.0–7.6$ galaxies have absolute UV magnitudes in the range -19.3 ≤ $M_{UV}$ ≤ -18.2 and [O III] + H β EWs ranging between ≈400–2000 Å. Moreover, there are two UV-bright (-20.4 ≤ $M_{UV}$ ≤ -20.2) galaxies located approximately 1 arcmin to the south of the UV-faint LAE JADES-GS-z7-LA.

We also identify five candidate $z \approx 7.0–7.6$ galaxies lying near (< 1.5 arcmin separation) the UV-bright LAE in GOODS-N (z8_GND_5296), all located roughly to the southeast (see Fig. 15).

<figure><img src="image_18.png" /><figcaption>**Figure 15.** On-sky distributions of candidate $z \approx 7.0–7.6$ galaxies in the two GOODS fields (circles) with the known high-confidence $z > 7$ LAE in each field shown as a star. JADES-GS-z7-LA is a very UV-faint ($M_{UV} = -16.7$) galaxy at $z = 7.28$ with near-unity Ly α escape (EW = 400 ± 90 Å; Saxena et al. 2023) while z8_GND_5296 is a UV-bright ($M_{UV} = -21.5$) galaxy at $z = 7.51$ with a less extreme Ly α EW of 33 ± 4 Å (Finkelstein et al. 2013; Jung et al. 2020). We identify a dense group of six faint (-19.3 ≤ $M_{UV}$ ≤ -18.2) neighbouring $z \approx 7.0–7.6$ galaxies approximately 1 arcmin to the east of JADES-GS-z7-LA, as well as two bright ($M_{UV} \approx -20.3$) such galaxies ≈1 arcmin to the south.</figcaption></figure>
Considering the smallest rectangular area that encloses all six of these objects, we obtain a surface density of 3.6 galaxies/arcmin⁻² which is approximately 4× higher than the average over the full JADES/GOODS-N footprint with F335M. These five relatively nearby objects have inferred [O III] + Hβ EWs of ≈500–1500 Å and UV magnitudes of −20.0 ≤ M<sub>UV</sub> ≤ −18.4, indicating that they are all substantially (≥ 4×) fainter than z8_GND_5296.

We now quantify the photometric overdensities around the two strong $z = 7.3 - 7.5$ LAEs in JADES by comparing the nearby $z \approx 7.0 - 7.6$ galaxy counts to that expected from literature luminosity functions at this epoch. Our procedure largely follows that of Whitler et al. (2024) which we outline below. As a first step towards comparing with the expected cosmic mean density, we must first correct our measured $z \approx 7.0 - 7.6$ galaxy counts for selection completeness. To this end, we analytically estimate our selection completeness as a function of true redshift and absolute UV magnitude (as in Endsley et al. 2021a; Whitler et al. 2024) by generating mock SEDs each with a flat continuum (in $F_v$) and IGM attenuation incorporated using the analytic model from Inoue et al. (2014). Given that our selection for this overdensity analysis depends on [O III] + Hβ emission via the F356W–F410M and F335M–F356W colour cuts above, we also add [O III] + Hβ emission into the mock SEDs following our results on the $z \sim 7 - 9$ EW distribution in Section 4.2, considering additional constraints on the $z \sim 7$ EW distribution among the very bright ($-22.5 \lesssim M_{UV} \lesssim -21$) population from Endsley et al. (2023). Specifically, for a given $M_{UV}$, we assume a lognormal [O III] + Hβ EW distribution defined by median EW log($\mu_{EW4pt}/\AA$) = $\max([2.87, -0.16(M_{UV} + 20) + 2.87])$ and standard deviation $\sigma_{EW4pt} = \min([0.25, 0.08(M_{UV} + 20) + 0.3])$ dex and adopt a fixed [O III]$\lambda5007/H\beta$ ratio of 6 (e.g. Tang et al. 2019). In doing so, we account for our strong incompleteness to objects with relatively weak [O III] + Hβ emission (EW $\lesssim 400$ Å) for this overdensity analysis, though we note that such objects are expected to be in the minority at $F150W < 29$ ($M_{UV} \lesssim -18$) given the derived EW distributions.

In the selection completeness simulations, for each grid point in redshift ($z = 6.5 - 8.5$, 0.05 unit spacing) and absolute UV magnitude ($-22 \le M_{UV} \le -17$, 0.25 mag spacing), we generate 1000 mock SEDs each with a different [O III] + Hβ EW value pulled randomly from the parametrized distribution above. To account for photometric scatter in our selection, we add Gaussian noise to the ‘true’ photometry derived from each mock SED. The photometric noise is computed as a function of $M_{UV}$ separately for each of the four different imaging regions (medium, deep, and very deep in GOODS-S as well as medium in GOODS-N) by fitting a linear relationship between $M_{UV}$ and the logarithm of the photometric error in each band using data from all $z \sim 7 - 9$ galaxies selected in each imaging region. We find that such a linear relationship adequately captures the general trend between photometric uncertainty and $M_{UV}$ imposed by typically larger galaxy sizes (and hence Kron apertures) at brighter luminosities.

The resulting simulated completeness of our $z \approx 7.0 - 7.6$ selection is shown in Fig. 16 for each of the four different JADES/NIRCam imaging regions. Consistent with our targeted redshift interval, the estimated completeness is $\gtrsim 50$ per cent at $z \approx 7.05 - 7.60$ at the brightest luminosities and very low at redshifts $z < 6.9$ and $> 7.7$ for all four regions. We also find that the completeness is $\gtrsim 50$ per cent at absolute UV magnitudes of $M_{UV} \lesssim -18.5$ in the very deep and deep imaging regions in GOODS-S, while such relatively high completeness is achieved at a slightly brighter magnitude threshold of $M_{UV} \lesssim -19$ in the medium-depth regions of GOODS-S and GOODS-N. These selection completeness grids are convolved with the UV luminosity function of Bouwens et al. (2022b) to estimate

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Simulated completeness of the $z \approx 7.0 - 7.6$ selection criteria where we supplement the base $z \sim 7 - 9$ selection criteria (Section 2.3) with additional colour cuts that utilize the long-wavelength medium-band filters to identify [O III] + Hβ emitters within this redshift range (see Section 6). We show the simulated completeness as a function of $M_{UV}$ and redshift in the four different JADES/NIRCam imaging regions where the cut-offs at $M_{UV} \approx -18$ are due to the requirement of $F150W < 29$ for this selection.</figcaption></figure>

the expected surface density counts of $z \approx 7.0 - 7.6$ galaxies in each imaging region under the assumption of cosmic mean density. Following the methods of Whitler et al. (2024), we compute the amplitude of the photometric overdensities around JADES-GS-z7-LA and z8_GND_5296 as a function of angular separation by dividing the GOODS-S and GOODS-N footprints into rings of varying concentric circular radii centered on each $z > 7$ LAE. In each ring, the amplitude of the photometric overdensity at that angular separation is computed as the surface density of objects identified by our $z \approx 7.0 - 7.6$ selection divided by the cosmic mean.

Given its near-unity Ly$\alpha$ escape fraction and relatively small Ly$\alpha$ velocity offset ($\approx 120$ km s$^{-1}$), JADES-GS-z7-LA is expected to reside in a very large ($R \gtrsim 3$ physical Mpc) ionized bubble (Saxena et al. 2023). As discussed in Saxena et al. (2023), this object cannot have created such a large bubble on its own given its very low luminosity ($M_{UV} = -16.7$), implying that the local ionizing photon budget must be dominated by neighbouring galaxies. With the JADES/GOODS-S imaging, we have identified a very dense group of $z \approx 7.0 - 7.6$ galaxies located $\approx 50$ arcsec ($\approx 250$ kpc in projection) east of JADES-GS-z7-LA (Fig. 15) where the surrounding photometric overdensity is very large ($\approx 25$) on small scales (see Fig. 17). It is conceivable that these objects all lie at very similar redshifts as JADES-GS-z7-LA and are assisting in powering a large ionized bubble in the vicinity (e.g. Lu et al. 2024). The two UV-bright ($M_{UV} \approx -20.3$) systems located at a very similar projected distance to the south of JADES-GS-z7-LA (Fig. 15) may be significantly contributing as well. While the surface density of $z \approx 7.0 - 7.6$ galaxies measured in the $\lesssim 2$ arcmin radius ($\approx 600$ kpc projected) around JADES-GS-z7-LA itself is consistent with the cosmic mean
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Measurements of photometric overdensity as a function of concentric angular radii around the two strong $z > 7$ LAEs covered by the JADES/NIRCam imaging. Neither object shows a remarkably strong photometric overdensity across the $z \approx 7.0-7.6$ redshift range probed by our photometric selection when centring the concentric rings on the LAEs themselves. However, there is a dense group of $z \approx 7.0-7.6$ galaxies near JADES-GS-z7-LA (see Fig. 15) where the photometric overdensity is very high ($\approx$25) on small scales ($<0.2$ arcmin). The underdensity in GOODS-S on large scales is consistent with previous findings (Bouwens et al. 2015; Finkelstein et al. 2015).</figcaption></figure>

(see Fig. 17), this can easily rise by a factor of $\sim$2–3 if the nearby galaxies are largely located near JADES-GS-z7-LA along the line of sight (i.e. $\Delta z \lesssim 0.1$).

In contrast to the UV-faint LAE described above, we do not identify any significant photometric overdensity near the UV-bright ($M_{UV} = -21.5$) $z = 7.51$ LAE z8_GND_5296 (Fig. 17). The surface density of $z \approx 7.0-7.6$ objects on both small ($<1.5$ arcmin; $<450$ kpc projected) and large ($6-10$ arcmin; $\sim 2-3$ Mpc projected) separations from z8_GND_5296 are consistent with the cosmic mean. This may indicate that the strong Ly$\alpha$ emission from z8_GND_5296 is due more to its physical properties. With an estimated [O III] + H$\beta$ EW of $1540_{-160}^{+220}$ Å, z8_GND_5296 lies in the high-end tail of the EW distribution for UV-bright $z \sim 7-9$ galaxies (see Fig. 8). Accordingly, this source has a large inferred ionizing photon production efficiency of $\log[\xi_{ion}/(\text{erg}^{-1} \text{Hz})] = 25.73_{-0.05}^{+0.06}$ (from the TcSFH BEAGLE SED fits) indicating that it has a relatively large intrinsic Ly$\alpha$ EW (e.g. Tang et al. 2019) thus requiring stronger attenuation (from either the host galaxy or IGM) to make it appear as a weak LAE.

The detection of the systemic [C III]λ1908 line from z8_GND_5296 may also indicate that its Ly$\alpha$ photons are redshifted far into the damping wing (velocity offset $\approx$400 km s$^{-1}$; Hutchison et al. 2019) as commonly found among similarly bright $z \sim 7$ galaxies (see Endsley et al. 2022 and references therein), resulting in relatively low Ly$\alpha$ opacity through the IGM (e.g. Stark et al. 2017; Mason et al. 2018). However, this is speculative as it is unclear which component of the [C III] doublet is currently detected given skyline contamination in the ground-based spectrum (Hutchison et al. 2019). Detailed spectroscopic measurements (e.g. Meyer et al. 2024; Witstok et al. 2024) are required to better assess the relative role of internal properties to the strong Ly$\alpha$ emission seen from this UV-bright LAE, as well as measure the redshifts of the neighbouring $z \approx 7.0-7.6$

candidates to better constrain its surrounding overdensity (see also Jung et al. 2020).

## 7 SUMMARY

We have used deep nine-band NIRCam imaging taken as part of the JADES program to characterize the star-forming and ionizing properties of a very large ($N = 756$) sample of Lyman-break $z \sim 6-9$ galaxies. This sample includes hundreds of reionization-era galaxies in the very UV-faint regime ($-18 \le M_{UV} \lesssim -16.5$; see Fig. 1) where it was effectively impossible to measure the rest-optical SEDs (at least among a statistical sample) with Spitzer/IRAC. Our main conclusions are summarized below.

(i) The faintest galaxies in our sample ($M_{UV} \sim -17$) tend to have inferred stellar masses of $(1-3) \times 10^7 \text{M}_\odot$ while those at the brightest end of our sample ($-22 \lesssim M_{UV} \lesssim -21$) have $\sim 100 \times$ higher stellar masses (Fig. 2). Masses inferred from models with priors weighted towards extended SFHs (i.e. ‘continuity’ priors) are systematically $\approx 0.5$ dex higher than those assuming CSFH, though can rise to factors of $10-100 \times$ higher in the most extreme cases for galaxies with the youngest light-weighted ages ($\text{age}_{\text{CSFH}} \sim 3$ Myr; see Fig. 2g).

(ii) There are no galaxies in our sample where the photometric data clearly imply very large stellar masses ($> 3 \times 10^{10} \text{M}_\odot$). There are only 13 galaxies in our sample with inferred stellar masses $> 3 \times 10^9 \text{M}_\odot$ when adopting the continuity prior that favours extended SFHs (Fig. 3). Some of these objects are relatively faint ($M_{UV} \sim -19$) and show strong Balmer breaks consistent with old stellar populations. Others are simply so luminous ($M_{UV} \sim -22$) such that a massive evolved stellar population can be hidden under the light of more recently formed stars. The most massive candidate in our sample has an inferred mass of $1 \times 10^{10} \text{M}_\odot$ when including rest-frame NIR SED constraints from parallel MIRI imaging (Fig. 4).

(iii) The typical galaxy in our sample shows a NIRCam SED consistent with a young light-weighted age of $\text{age}_{\text{CSFH}} \sim 50$ Myr (Fig. 5). However, we confidently identify several individual galaxies showing strong Balmer breaks implying much older ages ($\text{age}_{\text{CSFH}} \sim 250-1000$ Myr; Fig. 6). Some of the objects with confident strong Balmer breaks show signs of recent star formation from significant emission line signatures while others show no indications of emission lines in the photometry. Follow-up spectroscopy is necessary to determine if any galaxies in this latter subpopulation may be (temporarily) quenched.

(iv) We find a strong, highly significant decline in the typical [O III] + H$\beta$ EWs of $z \sim 6-9$ galaxies towards lower UV luminosities (median EW $\approx$800 Å at $M_{UV} = -20$ yet $\approx 350$ Å at $M_{UV} = -17.5$; see Fig. 8 and Table 2). We verify that this [O III] + H$\beta$ EW trend with $M_{UV}$ is reflected in the typical NIRCam colours of our sample (Fig. 10). The [O III] + H$\beta$ EW distribution is also found to broaden considerably at lower UV luminosities (standard deviation $\approx 0.3$ dex at $M_{UV} = -20$ yet $\approx 0.5$ dex at $M_{UV} = -17.5$). We infer a slight ($\approx 0.05-0.1$ dex) decline in the typical [O III] + H$\beta$ EW between the $z \sim 6$ sample and the $z \sim 7-9$ sample at fixed $M_{UV}$.

(v) In remarkable contrast to the strong decline in [O III] + H$\beta$ EW with $M_{UV}$, we find that the H$\alpha$ EW distribution changes very weakly (if at all) across a factor of $\approx 10$ in UV luminosity among our $z \sim 6$ sample (where H$\alpha$ falls in F444W; Fig. 8). If we assume that the $z \sim 6$ H$\alpha$ EW distribution does not change with $M_{UV}$ over the range probed by our $z \sim 6$ sample, we infer a lognormal distribution with median EW = $630_{-40}^{+10}$ Å and a standard deviation of $0.26 \pm 0.01$ dex.
(vi) We demonstrate that the [O III] + Hβ and Hα EW trends with M<sub>UV</sub> can be explained by a combination of lower metallicity and systematically more recently declining SFHs at lower UV luminosities (top panels of Fig. 12). In this interpretation, the brightest z ~ 6 galaxies in our sample (⟨M<sub>UV</sub>⟩ = −20.0) have a median inferred metallicity ≈ 3× higher than that of the faintest z ~ 6 galaxies (⟨M<sub>UV</sub>⟩ = −20.0; Z ≈ 0.06Z<sub>☉</sub>). Moreover, the median inferred ratio of SFR averaged over the past 3 Myr to that averaged over the past 50 Myr (SFR<sub>3 Myr</sub>/SFR<sub>50 Myr</sub>) is found to be ≈ 3.5× higher among the brightest z ~ 6 objects relative to the faintest galaxies. The brightest galaxies are frequently inferred to be experiencing a recent strong upturn in SFR (i.e. (SFR<sub>3 Myr</sub>/SFR<sub>50 Myr</sub> > 1), while the faintest galaxies are inferred to be a more even mixture of objects experiencing a recent strong rise in SFR as a recent strong downturn in SFR.

(vii) There are two particular subsets of our sample which provide evidence in favour of bursty SFHs. A substantial fraction of galaxies in our sample (≈10 per cent–20 per cent depending on M<sub>UV</sub>) are inferred to exhibit extremely high [O III] + Hβ EWs (>1500 Å; Fig. 7) implying that they have recently experienced a dramatic (factor ≳ 5) increase in SFR over the past 3 Myr (upper panels of Fig. 13). Another subset of galaxies in our JADES sample show relatively weak Balmer breaks yet also weak nebular line signatures (Fig. 14) implying that they may be experiencing a lull in star formation activity that followed a burst of star formation that occurred ~5–50 Myr ago (lower panels of Fig. 13). These two subpopulations may therefore be caught during the peaks and trough of bursty SFHs among reionization-era galaxies. We note that other recent studies are building empirical evidence for bursty SFHs among early galaxies (e.g. Strait et al. 2023; Looser et al. 2023, 2024; Dressler et al. 2024; Dome et al. 2024).

(viii) We also discuss how the [O III] + Hβ and Hα EW trends with M<sub>UV</sub> may be influenced by a strong correlation between UV luminosity and LyC escape fraction (bottom panels of Fig. 12). When allowing for substantial LyC escape yet enforce CSFHs, we infer that the faintest z ~ 6 galaxies in our sample (⟨M<sub>UV</sub>⟩ = −17.5) are typically very efficient at leaking LyC photons into the IGM (median f<sub>esc</sub>≈0.5), while the brightest (⟨M<sub>UV</sub>⟩ = −20.0) objects leak only a moderate fraction of their LyC photons (median f<sub>esc</sub>≈0.08). We discuss how this scenario has very different implications for the contribution of galaxies along the luminosity function to cosmic reionization compared to the interpretation of bursty SFHs. This highlights the need for deep spectroscopic follow-up to better determine the physical origin of the [O III] + Hβ and Hα EW trends with UV luminosity found in this analysis.

(ix) Finally, we quantify the photometric overdensities around two strong LAEs at z > 7 in the JADES footprint. The very UV-faint (M<sub>UV</sub> = −16.7) z = 7.28 LAE lies close to a very dense concentration of z ≈ 7.0–7.6 galaxies (Fig. 15) that may have helped generate a large ionized bubble leading to efficient Lyα transmission through the largely neutral IGM. However, the UV-bright (M<sub>UV</sub> = −21.5) LAE shows no significant nearby overdensity (Fig. 17) perhaps suggesting that efficient ionizing photon production and a large Lyα velocity offset strongly contributed to its Lyα detection, and that not all strong LAEs in the reionization era necessarily occupy large ionized bubbles.

The JADES results presented here provide a valuable observational baseline against which to compare predictions of reionization-era galaxy properties from models (e.g. Mutch et al. 2016; Rosdahl et al. 2018; Wilkins et al. 2020, 2023; Wu et al. 2020; Ceverino et al. 2021; Hutter et al. 2021; Kannan et al. 2021; Lewis et al. 2022; Hirschmann

et al. 2023; Seeyave et al. 2023). Such future work will help build insight into the nature and assembly of the faintest (and brightest) z ~ 6 – 9 galaxies, and their role in cosmic reionization.

ACKNOWLEDGEMENTS

The authors thank the anonymous referees for constructive comments which improved the quality of this paper. RE thanks John Chisholm, Charlotte Mason, and Adele Plat for enlightening discussions that improved this work. EE, DJE, BDJ, BR, GR, MR, FS, DPS, LW, and CNAW acknowledge support from JWST/NIRCam contract to the University of Arizona/NAS5-02015. DPS acknowledges additional support from the National Science Foundation through the grant AST-2109066. LW acknowledges additional support from the National Science Foundation Graduate Research Fellowship under grant no. DGE-2137419. WB acknowledges support by the Science and Technology Facilities Council (STFC), ERC advanced grant 695671 ‘QUENCH’. KB acknowledges support by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. AJB, AJC, JC, IEBW, AS, and GCJ acknowledge funding from the ‘First-Galaxies’ Advanced Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 789056). SC acknowledges support by European Union’s HE ERC Starting grant no. 101040227 – WINGS. ECL acknowledges support of an STFC Webb Fellowship (ST/W001438/1) ALD thanks the University of Cambridge Harding Distinguished Postgraduate Scholars Programme and Technology Facilities Council (STFC) Center for Doctoral Training (CDT) in Data intensive science at the University of Cambridge (STFC grant no. 2742605) for a PhD studentship. DJE is also supported as a Simons Investigator. TJL and RM acknowledge support by the Science and Technology Facilities Council (STFC) and by the ERC through advanced grant 695671 ‘QUENCH’. RM also acknowledges funding from a research professorship from the Royal Society. DP acknowledges support by the Huo Family Foundation through a P.C. Ho PhD Studentship. RS acknowledges support from a STFC Ernest Rutherford Fellowship (ST/S004831/1). The research of CCW is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. JW acknowledges support from the ERC advanced grant 695671, ‘QUENCH’, and the Fondation MERAC.

This material is based in part upon High Performance Computing (HPC) resources supported by the University of Arizona TRIF, UITS, and Research, Innovation, and Impact (RII) and maintained by the UArizona Research Technologies department. The authors acknowledge use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST-1828315.

This work made use of the following software: NUMPY (Harris et al. 2020); MATPLOTLIB (Hunter 2007); SCIPY (Virtanen et al. 2020); ASTROPY,<fnref n="8" /> a community-developed core PYTHON package for Astronomy (Astropy Collaboration 2013; Price-Whelan et al. 2018); SOURCE EXTRACTOR (Bertin & Arnouts 1996) via SEP (Barbary 2016); PHOTUTILS (Bradley et al. 2022); BEAGLE (Chevallard & Charlot 2016); MULTINEST (Feroz & Hobson 2008; Feroz et al. 2009); PROSPECTOR (Johnson et al. 2021b); DYNESTY (Speagle 2020); SEDPY (Johnson 2021); and FSPS (Conroy et al. 2009; Conroy & Gunn 2010) via PYTHON-FSPS (Johnson et al. 2021a).

<footnote n="8">⁸https://www.astropy.org/</footnote>
DATA AVAILABILITY

The HST ACS data used in this work are available from the HLF
archive (https://archive.stsci.edu/prepds/hlf/). All raw JWST JADES
data utilized in this work is now publicly available via the Mikulski
Archive for Space Telescopes (https://mast.stsci.edu/portal/Mashup
/Clients/Mast/Portal.html), with a subset released as part of high-
level science products (https://archive.stsci.edu/hlsp/jades).

REFERENCES

Akins H. B. et al., 2023, ApJ, 956, 61

Alberts S. et al., 2024, preprint (arXiv:2405.15972)

Arrabal Haro P. et al., 2023, ApJ, 951, L22

Astropy Collaboration, 2013, A&A, 558, A33

Atek H. et al., 2015, ApJ, 800, 18

Atek H. et al., 2024, Nature, 626, 975

Bagley M. B. et al., 2024, ApJ, 965, L6

Bagley M. B. et al., 2023, ApJ, 946, L12

Baker W. M. et al., 2023, preprint (arXiv:2306.02472)

Barbary K., 2016, J. Open Source Softw., 1, 58

Barro G. et al., 2024, ApJ, 963, 128

Barrow K. S. S., Robertson B. E., Ellis R. S., Nakajima K., Saxena A., Stark D. P., Tang M., 2020, ApJ, 902, L39

Barrufet L. et al., 2023, MNRAS, 522, 449

Bertin E., Arnouts S., 1996, A&AS, 117, 393

Bezanson R. et al., 2022, preprint (arXiv:2212.04026)

Bhatawdekar R., Conselice C. J., 2021, ApJ, 909, 144

Bhatawdekar R., Conselice C. J., Margalef-Bentabol B., Duncan K., 2019, MNRAS, 486, 3805

Binggeli C. et al., 2019, MNRAS, 489, 3827

Böker T. et al., 2023, PASP, 135, 038001

Bosman S. E. I. et al., 2022, MNRAS, 514, 55

Bouwens R. J. et al., 2004, ApJ, 616, L79

Bouwens R. J. et al., 2014, ApJ, 793, 115

Bouwens R. J. et al., 2015, ApJ, 803, 34

Bouwens R. J. et al., 2023, MNRAS, 523, 1036

Bouwens R. J., Illingworth G. D., van Dokkum P. G., Oesch P. A., Stefanon M., Ribeiro B., 2022a, ApJ, 927, 81

Bouwens R. J., Illingworth G., Ellis R. S., Oesch P., Stefanon M., 2022b, ApJ, 940, 55

Bowler R. A. A. et al., 2014, MNRAS, 440, 2810

Bowler R. A. A., Dunlop J. S., McLure R. J., McLeod D. J., 2017, MNRAS, 466, 3612

Bowler R. A. A., Jarvis M. J., Dunlop J. S., McLure R. J., McLeod D. J., Adams N. J., Milvang-Jensen B., McCracken H. J., 2020, MNRAS, 493, 2059

Boyer M. L. et al., 2022, Res. Notes Am. Astron. Soc., 6, 191

Boyett K. N. K., Stark D. P., Bunker A. J., Tang M., Maseda M. V., 2022, MNRAS, 513, 4451

Boylan-Kolchin M., 2023, Nat. Astron., 7, 731

Bradley L. et al., 2022, astropy/photutils: 1.5.0, available at: https://doi.org/10.5281/zenodo.6825092

Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127

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

Bunker A. J., Stanway E. R., Ellis R. S., McMahon R. G., 2004, MNRAS, 355, 374

Burgasser A. J., 2014, in Singh H. P., Prugniel P., Vauglin I., eds, ASI Conf. Ser. 11, International Workshop on Stellar Spectral Libraries, p. 7, preprint (arXiv:1406.4887)

Byler N., Dalcanton J. J., Conroy C., Johnson B. D., 2017, ApJ, 840, 44

Cameron A. J. et al., 2023, A&A, 677, A115

Carnall A. C. et al., 2023, MNRAS, 518, L45

Carnall A. C., Leja J., Johnson B. D., McLure R. J., Dunlop J. S., Conroy C., 2019, ApJ, 873, 44

Castellano M. et al., 2016, ApJ, 818, L3

Castellano M. et al., 2022, A&A, 662, A115

Ceverino D., Hirschmann M., Klessen R. S., Glover S. C. O., Charlot S., Feltre A., 2021, MNRAS, 504, 4472

Ceverino D., Klessen R. S., Glover S. C. O., 2018, MNRAS, 480, 4842

Chabrier G., 2003, PASP, 115, 763

Chen Y., Bressan A., Girardi L., Marigo P., Kong X., Lanza A., 2015, MNRAS, 452, 1068

Chen Z., Stark D. P., Endsley R., Topping M., Whitler L., Charlot S., 2023, MNRAS, 518, 5607

Chevallard J., Charlot S., 2016, MNRAS, 462, 1415

Chisholm J. et al., 2022, MNRAS, 517, 5104

Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102

Conroy C., Gunn J. E., 2010, ApJ, 712, 833

Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486

Conselice C. J., Gallagher J. S., Calzetti D., Homeier N., Kinney A., 2000, AJ, 119, 79

Curti M. et al., 2024, A&A, 684, A75

Curtis-Lake E. et al., 2023, Nat. Astron., 7, 622

Davies F. B. et al., 2018, ApJ, 864, 142

Dayal P., Ferrara A., 2018, Phys. Rep., 780-782, 1

De Barros S., Oesch P. A., Labbé I., Stefanon M., González V., Smit R., Bouwens R. J., Illingworth G. D., 2019, MNRAS, 489, 2355

Dome T., Tacchella S., Fialkov A., Dekel A., Ginzburg O., Lapiner S., Looser T. J., 2024, MNRAS, 527, 2139

Donnan C. T. et al., 2023, MNRAS, 518, 6011

Dressler A. et al., 2023, ApJ, 947, L27

Dressler A. et al., 2024, ApJ, 964, 150

Duncan K. et al., 2014, MNRAS, 444, 2960

Dunlop J. S., McLure R. J., Robertson B. E., Ellis R. S., Stark D. P., Cirasuolo M., de Ravel L., 2012, MNRAS, 420, 901

Eisenstein D. J. et al., 2023, preprint (arXiv:2306.02465)

Ellis R. S. et al., 2013, ApJ, 763, L7

Endsley R. et al., 2022, MNRAS, 517, 5642

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

Endsley R., Stark D. P., Charlot S., Chevallard J., Robertson B., Bouwens R. J., Stefanon M., 2021b, MNRAS, 502, 6044

Endsley R., Stark D. P., Chevallard J., Charlot S., 2021a, MNRAS, 500, 5229

Endsley R., Stark D. P., Whitler L., Topping M. W., Chen Z., Plat A., Chisholm J., Charlot S., 2023, MNRAS, 524, 2312

Faisst A. L., Capak P. L., Emami N., Tacchella S., Larson K. L., 2019, ApJ, 884, 133

Faucher-Giguère C.-A., 2018, MNRAS, 473, 3717

Feroz F., Hobson M. P., 2008, MNRAS, 384, 449

Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601

Finkelstein S. L. et al., 2012, ApJ, 756, 164

Finkelstein S. L. et al., 2013, Nature, 502, 524

Finkelstein S. L. et al., 2015, ApJ, 810, 71

Finkelstein S. L. et al., 2022, ApJ, 928, 52

Finkelstein S. L. et al., 2023, ApJ, 946, L13

Fletcher T. J., Tang M., Robertson B. E., Nakajima K., Ellis R. S., Stark D. P., Inoue A., 2019, ApJ, 878, 87

Flury S. R. et al., 2022, ApJ, 930, 126

Fujimoto S. et al., 2023, ApJ, 949, L25

Furlanetto S. R., Mirocha J., 2022, MNRAS, 511, 3895

Furtak L. J. et al., 2023a, ApJ, 952, 142

Furtak L. J., Shuntov M., Atek H., Zitrin A., Richard J., Lehnert M. D., Chevallard J., 2023b, MNRAS, 519, 3064

Giavalisco M. et al., 2004, ApJ, 600, L93

González V., Bouwens R., Illingworth G., Labbé I., Oesch P., Franx M., Magee D., 2014, ApJ, 781, 34

Goto H. et al., 2021, ApJ, 923, 229

Greene J. E. et al., 2024, ApJ, 964, 39

Gutkin J., Charlot S., Bruzual G., 2016, MNRAS, 462, 1757

Hainline K. N. et al., 2024, ApJ, 964, 71

Harikane Y. et al., 2022, ApJS, 259, 20

Harikane Y. et al., 2023, ApJ, 959, 39
Harris C. R. et al., 2020, Nature, 585, 357

Hashimoto T. et al., 2018, Nature, 557, 392

Heckman T. M. et al., 2011, ApJ, 730, 5

Helton J. M. et al., 2024, ApJ, 962, 124

Hirschmann M. et al., 2023, MNRAS, 526, 3610

Hsiao T. Y.-Y. et al., 2023, ApJ, 949, L34

Huang K.-H. et al., 2016, ApJ, 817, 11

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Hutchison T. A. et al., 2019, ApJ, 879, 70

Hutter A., Dayal P., Müller V., Trott C. M., 2017, ApJ, 836, 176

Hutter A., Dayal P., Yepes G., Gottlöber S., Legrand L., Ucci G., 2021, MNRAS, 503, 3698

Illingworth G. D. et al., 2013, ApJS, 209, 6

Inoue A. K., Shimizu I., Iwata I., Tanaka M., 2014, MNRAS, 442, 1805

Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2018, ApJ, 854, 73

Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Chisholm J., Thuan T. X., Fricke K. J., Verhamme A., 2021, MNRAS, 503, 1734

Jakobsen P. et al., 2022, A&A, 661, A80

Jiang L. et al., 2022, Nat. Astron., 6, 850

Johnson B. D., 2021, bd-j/sedpy: sedpy v0.2.0, https://doi.org/10.5281/zenodo.4582723

Johnson B. D., Leja J., Conroy C., Speagle J. S., 2021b, ApJS, 254, 22

Johnson B. et al., 2021a, dmf/python-fsps: python-fsps v0.4.1rc1, Zenodo, doi:10.5281/zenodo.4737461

Jung I. et al., 2019, ApJ, 877, 146

Jung I. et al., 2020, ApJ, 904, 144

Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M., Hernquist L., 2022, MNRAS, 511, 4005

Kikuchihara S. et al., 2020, ApJ, 893, 60

Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900

Kocevski D. D. et al., 2023, ApJ, 954, L4

Kron R. G., 1980, ApJS, 43, 305

Labbé I. et al., 2013, ApJ, 777, L19

Labbé I. et al., 2015, ApJS, 221, 23

Labbé I. et al., 2023, Nature, 616, 266

Lam D. et al., 2019, A&A, 627, A164

Langan I., Ceverino D., Finlator K., 2020, MNRAS, 494, 1988

Laporte N., Ellis R. S., Witten C. E. C., Roberts-Borsani G., 2023, MNRAS, 523, 3018

Laporte N., Meyer R. A., Ellis R. S., Robertson B. E., Chisholm J., Roberts-Borsani G. W., 2021, MNRAS, 505, 3336

Laporte N., Nakajima K., Ellis R. S., Zitrin A., Stark D. P., Mainali R., Roberts-Borsani G. W., 2017, ApJ, 851, 40

Larson R. L. et al., 2022, ApJ, 930, 104

Larson R. L. et al., 2023, ApJ, 953, L29

Lee J. C. et al., 2012, PASP, 124, 782

Lee J. C., Kennicutt R. C., Funes S. J. J. G., Sakai S., Akiyama S., 2007, ApJ, 671, L113

Leethochawalit N. et al., 2023, ApJ, 942, L26

Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019, ApJ, 876, 3

Leonova E. et al., 2022, MNRAS, 515, 5790

Lewis J. S. W. et al., 2022, MNRAS, 516, 3389

Livermore R. C., Finkelstein S. L., Lotz J. M., 2017, ApJ, 835, 113

Looser T. J. et al., 2023, preprint (arXiv:2306.02470)

Looser T. J. et al., 2024, Nature, 629, 53

Lovell C. C., Harrison I., Harikane Y., Tacchella S., Wilkins S. M., 2023, MNRAS, 518, 2511

Lower S., Narayanan D., Leja J., Johnson B. D., Conroy C., Davé R., 2020, ApJ, 904, 33

Lu T.-Y., Mason C., Hutter A., Mesinger A., Qin Y., Stark D. P., Endsley R., 2024, MNRAS, 528, 4872

Ly C., Lee J. C., Dale D. A., Momcheva I., Salim S., Staudaher S., Moore C. A., Finn R., 2011, ApJ, 726, 109

Ma X. et al., 2018, MNRAS, 478, 1694

Ma X., Hopkins P. F., Faucher-Giguère C.-A., Zolman N., Muratov A. L., Kereš D., Quataert E., 2016, MNRAS, 456, 2140

Ma X., Quataert E., Wetzel A., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2020, MNRAS, 498, 2001

Mainali R. et al., 2018, MNRAS, 479, 1180

Mainali R., Kollmeier J. A., Stark D. P., Simcoe R. A., Walth G., Newman A. B., Miller D. R., 2017, ApJ, 836, L14

Maiolino R. et al., 2023, preprint (arXiv:2308.01230)

Maiolino R. et al., 2024, Nature, 627, 59

Marley M. S. et al., 2021, ApJ, 920, 85

Mason C. A., Trenti M., Treu T., 2023, MNRAS, 521, 497

Mason C. A., Treu T., Dijkstra M., Mesinger A., Trenti M., Pentericci L., de Barros S., Vanzella E., 2018, ApJ, 856, 2

Matsuoka Y. et al., 2023, ApJ, 949, L42

Matthee J. et al., 2024, ApJ, 963, 129

McLeod D. J., McLure R. J., Dunlop J. S., 2016, MNRAS, 459, 3812

McLure R. J. et al., 2013, MNRAS, 432, 2696

Merlin E. et al., 2022, ApJ, 938, L14

Meyer R. A. et al., 2024, preprint (arXiv:2405.05111)

Mirocha J., Furlanetto S. R., 2023, MNRAS, 519, 843

Morishita T., Stiavelli M., 2023, ApJ, 946, L35

Mutch S. J., Geil P. M., Poole G. B., Angel P. W., Duffy A. R., Mesinger A., Wyithe J. S. B., 2016, MNRAS, 462, 250

Nasir F., D'Aloisio A., 2020, MNRAS, 494, 3080

Neufeld C. et al., 2022, MNRAS, 516, 2162

Oesch P. A. et al., 2013, ApJ, 773, 75

Oesch P. A. et al., 2023, MNRAS, 525, 2864

Oesch P. A., Bouwens R. J., Illingworth G. D., Labbé I., Stefanon M., 2018, ApJ, 855, 105

Ono Y. et al., 2012, ApJ, 744, 83

Ono Y. et al., 2018, PASJ, 70, S10

Onoue M. et al., 2023, ApJ, 942, L17

Papovich C. et al., 2023, ApJ, 949, L18

Pei Y. C., 1992, ApJ, 395, 130

Pérez-González P. G. et al., 2024, ApJ, 968, 4

Planck Collaboration VI, 2020, A&A, 641, A6

Plat A., Charlot S., Bruzual G., Feltre A., Vidal-García A., Morisset C., Chevallard J., Todt H., 2019, MNRAS, 490, 978

Price-Whelan A. M. et al., 2018, AJ, 156, 123

Prieto-Lyon G. et al., 2023, A&A, 672, A186

Rasappu N., Smit R., Labbé I., Bouwens R. J., Stark D. P., Ellis R. S., Oesch P. A., 2016, MNRAS, 461, 3886

Rieke M. J. et al., 2023, PASP, 135, 028001

Rieke M., Kelly D., Horner S., NIRCam Team, 2005, American Astronomical Society Meeting Abstracts, 207, 115.09

Rieke M., the JADES Collaboration, 2023, ApJS, 269, 16

Rigby J. et al., 2023, PASP, 135, 048001

Rinaldi P. et al., 2023, ApJ, 952, 143

Rivera-Thorsen T. E. et al., 2017, A&A, 608, L4

Roberts-Borsani G. W. et al., 2016, ApJ, 823, 143

Roberts-Borsani G. W., Ellis R. S., Laporte N., 2020, MNRAS, 497, 3440

Robertson B. E. et al., 2023, Nat. Astron., 7, 611

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

Rodighiero G., Bisigello L., Iani E., Marasco A., Grazian A., Sinigaglia F., Cassata P., Gruppioni C., 2023, MNRAS, 518, L19

Rogers A. B., McLure R. J., Dunlop J. S., 2013, MNRAS, 429, 2456

Rosdahl J. et al., 2018, MNRAS, 479, 994

Salmon B. et al., 2020, ApJ, 889, 189

Sanders R. L., Shapley A. E., Topping M. W., Reddy N. A., Brammer G. B., 2023, ApJ, 955, 54

Saxena A. et al., 2023, A&A, 678, A68

Schaerer D., de Barros S., 2009, A&A, 502, 423

Schaerer D., de Barros S., 2010, A&A, 515, A73

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

Schlawin E., Leisenring J., Misselt K., Greene T. P., McElwain M. W., Beatty T., Rieke M., 2020, AJ, 160, 231

Schmidt K. B. et al., 2017, ApJ, 839, 17

Schouws S. et al., 2022, ApJ, 928, 31
Seeyave L. T. C. et al., 2023, MNRAS, 525, 2422

Shen X., Hopkins P. F., Faucher-Giguère C.-A., Alexander D. M., Richards
G. T., Ross N. P., Hickox R. C., 2020, MNRAS, 495, 3252

Shibuya T., Ouchi M., Harikane Y., 2015, ApJS, 219, 15

Simmonds C. et al., 2023, MNRAS, 523, 5468

Smit R. et al., 2014, ApJ, 784, 58

Smit R. et al., 2015, ApJ, 801, 122

Smit R., Bouwens R. J., Labbé I., Franx M., Wilkins S. M., Oesch P. A., 2016, ApJ, 833, 254

Sommovigo L. et al., 2022, MNRAS, 513, 3122

Song M. et al., 2016, ApJ, 825, 5

Speagle J. S., 2020, MNRAS, 493, 3132

Stanway E. R., Bunker A. J., McMahon R. G., 2003, MNRAS, 342, 439

Stark D. P. et al., 2015a, MNRAS, 450, 1846

Stark D. P. et al., 2015b, MNRAS, 454, 1393

Stark D. P. et al., 2017, MNRAS, 464, 469

Stark D. P., 2016, ARA&A, 54, 761

Stark D. P., Schenker M. A., Ellis R., Robertson B., McLure R., Dunlop J., 2013, ApJ, 763, 129

Stefanon M. et al., 2017, ApJ, 851, 43

Stefanon M. et al., 2019, ApJ, 883, 99

Stefanon M. et al., 2021a, ApJS, 257, 68

Stefanon M., Bouwens R. J., Illingworth G. D., Labbé I., Oesch P. A., Gonzalez V., 2022b, ApJ, 935, 94

Stefanon M., Bouwens R. J., Labbé I., Illingworth G. D., Gonzalez V., Oesch P. A., 2021b, ApJ, 922, 29

Stefanon M., Bouwens R. J., Labbé I., Illingworth G. D., Oesch P. A., van Dokkum P., Gonzalez V., 2022a, ApJ, 927, 48

Strait V. et al., 2020, ApJ, 888, 124

Strait V. et al., 2023, ApJ, 949, L23

Sun G., Faucher-Giguère C.-A., Hayward C. C., Shen X., 2023, MNRAS,
526, 2665

Tacchella S. et al., 2022, ApJ, 927, 170

Tacchella S. et al., 2023a, ApJ, 952, 74

Tacchella S. et al., 2023b, MNRAS, 522, 6236

Tang M. et al., 2023, MNRAS, 526, 1657

Tang M., Stark D. P., Chevallard J., Charlot S., 2019, MNRAS, 489, 2572

Tang M., Stark D. P., Chevallard J., Charlot S., Endsley R., Congiu E., 2021, MNRAS, 501, 3238

Tang M., Stark D. P., Ellis R. S., 2022, MNRAS, 513, 5211

Tilvi V. et al., 2016, ApJ, 827, L14

Tilvi V. et al., 2020, ApJ, 891, L10

Topping M. W. et al., 2022a, MNRAS, 516, 975

Topping M. W. et al., 2024, MNRAS, 529, 4087

Topping M. W., Stark D. P., Endsley R., Plat A., Whitler L., Chen Z., Charlot S., 2022b, ApJ, 941, 153

Torrey P. et al., 2019, MNRAS, 484, 5587

Treu T. et al., 2022, ApJ, 935, 110

Treu T. et al., 2023, ApJ, 942, L28

Übler H. et al., 2023, A&A, 677, A145

Ucci G. et al., 2023, MNRAS, 518, 3557

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Weinberger L. H., Kulkarni G., Haehnelt M. G., Choudhury T. R., Puchwein E., 2018, MNRAS, 479, 2564

Whitaker K. E. et al., 2019, ApJS, 244, 16

Whitler L., Endsley R., Stark D. P., Topping M., Chen Z., Charlot S., 2023a,
MNRAS, 519, 157

Whitler L., Stark D. P., Endsley R., Chen Z., Mason C., Topping M. W., Charlot S., 2024, MNRAS, 529, 855

Whitler L., Stark D. P., Endsley R., Leja J., Charlot S., Chevallard J., 2023b,
MNRAS, 519, 5859

Wilkins S. M. et al., 2020, MNRAS, 493, 6079

Wilkins S. M. et al., 2023, MNRAS, 522, 4014

Williams C. C. et al., 2023, ApJS, 268, 64

Williams C. C. et al., 2024, ApJ, 968, 34

Withers S. et al., 2023, ApJ, 958, L14

Witstok J. et al., 2022, MNRAS, 515, 1751

Witstok J. et al., 2024, A&A, 682, A40

Wright G. S. et al., 2023, PASP, 135, 048003

Wu X., Davé R., Tacchella S., Lotz J., 2020, MNRAS, 494, 5636

Wyithe J. S. B., Loeb A., 2005, ApJ, 625, 1

Yang J. et al., 2020, ApJ, 897, L14

Zhu Y. et al., 2022, ApJ, 932, 76

APPENDIX A: CONSIDERING A POTENTIAL
POPULATION OF LOW-REDSHIFT
CONTAMINANTS

As discussed in Section 2.3, we choose to select galaxies based on Lyman-break colour criteria rather than photometric redshifts. Utilizing photometric redshifts can bias our sample preferentially towards objects with strong rest-optical lines since such lines imprint unique long-wavelength NIRCam colour patterns not easily reproduced with lower redshift solutions. None the less, we aim to quantify the potential low-redshift interloper fraction of our sample and assess whether the presence of such interlopers might impact our main conclusions. To this end, we calculate the photometric redshifts of all galaxies in our final Lyman-break cut sample but excluding bands which fall in the rest-optical regime. That is, we only fit using bands blueward of (and including) F150W and F200W for the F775W and F090W dropout samples, respectively. We use the BEAGLE TcSFH setup described in Section 2.4 and allow redshifts in the range $z = 0 - 10$ with a uniform prior.

For nearly every sample subset considered in Section 4, ≳90
per cent of galaxies have a high probability (>90 per cent) of
lying at z > 4 (see Fig. A1). The only exception is the very faint
($M_{UV} < -18$) subset of the F090W dropout sample, in which ≈80
per cent of the galaxies have a large probability (>80 per cent) of
lying at high redshifts (Fig. A1).

In Section 4.2, we found that the median [O III] + H β EW of reionization-era galaxies is substantially lower at fainter UV luminosities. Here, we verify that this result is not likely due to a potential moderate population of low-redshift contaminants in the very faint subsets. We account for the high-redshift probabilities of

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

**Figure A1.** Cumulative distribution of the fraction of galaxies in our final Lyman-break sample (Section 2.3) as a function of the probability that their photometric redshift is *z* > 4. Here, we exclude photometric bands covering the rest-optical regime of each dropout sample (see Appendix A). We plot each of the six galaxy subsamples used in Section 4 separately.
each object when inferring the [O III] + H β EW distributions by replacing equation (2) with the following:

$$P(\theta) \propto \prod_i \left[ \sum_j P_{i,j}(\text{EW}) P_j(\text{EW}|\theta) \right]^{P_i(z>4)} . \quad (A1)$$

By raising the term in brackets to the power of $P_i(z > 4)$ (i.e. the probability that galaxy $i$ lies at $z > 4$), we weight the contribution of each source to the EW distribution constraint based on the likelihood that it lies at high redshift.

With equation (A1), we infer median [O III] + H β EWs of 400$^{+70}_{-70}$ Å among the very faint F775W and F090W dropout subsets, respectively, using the TcSFH model outputs. These median EWs are entirely consistent with our fiducial values reported in Table 2. This indicates that a potential moderate (≈10 per cent) population of low-redshift interlopers in our faintest bins cannot explain the strong UV luminosity dependence on the median [O III] + H β found in Section 4.2. We have also verified that adopting equation (A1) does not change our conclusions on the lack of a strong UV luminosity dependence on the H α EW distribution at $z \sim 6$ (Section 4.3).

## APPENDIX B: ASSESSMENT OF THE GOODNESS OF SED FITS

Here, we explore the goodness of our BEAGLE SED fits which we utilize in our EW distribution analysis (Section 4). We show the best-fitting $\chi^2$ values of the F775W dropout ($z \sim 6$) and F090W dropout ($z \sim 7-9$) samples in the top and bottom panels, respectively, of Fig. B1. In the F775W dropout sample, it is clear that a significant fraction of objects are considerably better fit with the BEAGLE TcSFH models than the CSFH models, motivating us to adopt the EW distribution results from the TcSFH fits as fiducial. The reason why the TcSFH models provide substantially better fits to a subset of $z \sim 6$ galaxies is discussed in detail in Section 5.2. Briefly, these models allow for solutions where the continuum is consistent with a young stellar population (age $\lesssim 30$ Myr) while still yielding relatively weak [O III] + H β and H α emission from a recent downturn in SFR. Such a contrast in fit quality is not seen in the F090W dropout sample given that the NIRCam data are not sensitive to H α emission at $z \gtrsim 7$, which is necessary to determine if weak [O III] emission can be explained by extremely low metallicity (see Section 5.2).

We also check that the BEAGLE TcSFH models generally deliver acceptable fits to the rest-optical photometry alone. In Fig. B2, we show the distributions of best-fitting $\chi^2$ values limited to the rest-

<figure><img src="image_22.png" /><figcaption>Figure B1. A comparison of the best-fitting $\chi^2$ values (using all HST and NIRCam bands) between the BEAGLE TcSFH model fits (y-axis) and the BEAGLE CSFH fits (x-axis). We show the comparison for the F775W dropout ($z \sim 6$) sample and the F090W dropout ($z \sim 7-9$) sample in the top and bottom panels, respectively.</figcaption></figure>
<figure><img src="image_23.png" /></figure>

**Figure B2.** The distribution of best-fitting $\chi^2$ values from the BEAGLE TcSFH fits for each sample subset used in our [O III] + H$\beta$ and H$\alpha$ EW distribution analyses (Section 4). Here, we are only considering the bands that probe the rest-optical, so F200W (F277W) and redder for the F775W (F090W) dropout samples, resulting in six (five) fitted bands. The rest-optical photometry is well fitted by these models for the large majority of objects in each redshift and UV luminosity bin, with no strong differences in fit quality in the faintest bins.

optical bands, considering each redshift and UV luminosity subset adopted in Section 4 separately. We find that the best-fitting rest-optical $\chi^2$ values are less than the number of fitted rest-optical bands (six for F775W dropouts, 5 = five for F090W dropouts) for the large majority of objects (≥85 per cent) in every subset. Moreover, we find no strong differences in rest-optical fit quality among the faintest bins relative to the brighter bins.

This paper has been typeset from a TeX/LATEX file prepared by the author.

<footnote n="1">¹Department of Astronomy, University of Texas, Austin, TX 78712, USA</footnote>

<footnote n="2">²Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ 85721 USA</footnote>

<footnote n="3">³Center for Astrophysics | Harvard and Smithsonian, 60 Garden St., Cambridge, MA 02138, USA</footnote>

<footnote n="4">⁴Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 96054, USA</footnote>

<footnote n="5">⁵Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="6">⁶Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="7">⁷European Space Agency, ESA/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, NL</footnote>

<footnote n="8">⁸School of Physics, University of Melbourne, Parkville, VIC 3010, Australia</footnote>

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

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

<footnote n="11">¹¹Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

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

<footnote n="13">¹³Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield AL10 9AB, UK</footnote>

<footnote n="14">¹⁴Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="15">¹⁵Department for Astrophysical and Planetary Science, University of Colorado, Boulder, CO 80309, USA</footnote>

<footnote n="16">¹⁶Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany</footnote>

<footnote n="17">¹⁷Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="18">¹⁸Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK</footnote>

<footnote n="19">¹⁹NSF's National Optical-Infrared Astronomy Research Laboratory, 950 North Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="20">²⁰NRC Herzberg, 5071 West Saanich Rd, Victoria, BC V9E 2E7, Canada</footnote>