Reionization is thought to be driven by faint starforming galaxies but characterizing this population has long remained very challenging. Here we utilize deep nineband JADES JWST Advanced Deep Extragalactic SurveyNIRCam NearInfrared Camera imaging to study the starforming and ionizing properties of 756 zsim 69 galaxies including hundreds of very ultraviolet UVfaint objects MmathrmUVgt 18. The faintest msim 30 galaxies in our sample typically have stellar masses of Mast sim 13times 107 mathrm Modot and young lightweighted ages sim50 Myr though some show strong Balmer breaks implying much older ages sim500 Myr. We find no evidence for extremely massive galaxies gt 3times 1010 mathrm Modot in our sample. We infer a strong factor gt 2 decline in the typical O IIIH beta equivalent widths EWs towards very faint zsim 69 galaxies yet a weak UV luminosity dependence on the H alpha EWs at zsim 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 UVluminous 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 MmathrmUV 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 zgt7 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 zgt 7 Ly alpha emitters reside in large ionized bubbles.
['Astrophysics - Astrophysics of Galaxies']
The starforming and ionizing properties of dwarf z 69 galaxies in JADES insights on bursty star formation and ionized bubble growth
{'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 1 ★ , Daniel P. Stark 2 , Lily Whitler 2 , Michael W. Topping 2 , Benjamin D. Johnson 3 , Brant Robertson 4 , Sandro Tacchella 5 , 6 , Stacey Alberts 2 , William M. Baker 5 , 6 , Rachana Bhatawdekar 7 , Kristan Boyett 8 , 9 , Andrew J. Bunker 10 , Alex J. Cameron 10 , Stefano Carniani 11 , Stephane Charlot 12 , Zuyi Chen 2 , Jacopo Chevallard 10 , Emma Curtis-Lake 13 , A. Lola Danhaive 5 , 6 , Eiichi Egami 2 , Daniel J. Eisenstein 3 Kevin Hainline 2 , Jakob M. Helton 2 , Zhiyuan Ji 2 , Tobias J. Looser 5 , 6 , Roberto Maiolino 5 , 6 , 14 , Erica Nelson 15 , Dávid Puskás 5 , 6 , George Rieke 2 , Marcia Rieke 2 , Hans-Walter Rix 16 , Lester Sandles 5 , 6 Aayush Saxena 10 , 17 , Charlotte Simmonds 5 , 6 , Renske Smit 18 , Fengwu Sun 2 , Christina C. Williams 19 , Christopher N. A. Willmer 2 , Chris Willott 20 , and Joris Witstok 5 , 6  \n1 Department of Astronomy, University of Texas, Austin, TX 78712, USA  \n- 2 Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ 85721 USA\n- 3 Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge MA 02138 USA\n- 4 Department of Astronomy and Astrophysics University of California, Santa Cruz, 1156 High Street, Santa Cruz CA 96054, USA\n- 5 Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK\n- 6 Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK\n- 7 European Space Agency, ESA/ESTEC, Keplerlaan 1, 2201 AZ Noordwĳk, NL\n- 8 School of Physics, University of Melbourne, Parkville 3010, VIC, Australia\n- 9 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia\n- 10 Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK\n- 11 Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy\n- 12 Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France\n- 13 Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield AL10 9AB, UK\n- 14 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK\n- 15 Department for Astrophysical and Planetary Science, University of Colorado, Boulder, CO 80309, USA\n- 16 Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany\n- 17 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK\n- 18 Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK\n- 19 NSF's National Optical-Infrared Astronomy Research Laboratory, 950 North Cherry Avenue, Tucson, AZ 85719, USA\n- 20 NRC Herzberg, 5071 West Saanich Rd, Victoria, BC V9E 2E7, Canada  \nAccepted XXX. Received YYY; in original form ZZZ", '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/NIRCam imaging to study the star-forming and ionizing properties of 756 𝑧 ∼ 6 -9 galaxies, including hundreds of very UV-faint objects ( 𝑀 UV > -18). The faintest ( 𝑚 ∼ 30) galaxies in our sample typically have stellar masses of 𝑀 ∗ ∼ ( 1 -3 ) × 10 7 𝑀 ⊙ and young light-weighted ages ( ∼ 50 Myr), though some show strong Balmer breaks implying much older ages ( ∼ 500 Myr). We find no evidence for extremely massive galaxies ( > 3 × 10 10 𝑀 ⊙ ) in our sample. We infer a strong (factor > 2) decline in the typical [OIII] + H 𝛽 EWs towards very faint 𝑧 ∼ 6 -9 galaxies, yet a weak UV luminosity dependence on the H 𝛼 EWs at 𝑧 ∼ 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 SFR. We also discuss how the EW trends may be influenced by a strong correlation between 𝑀 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 𝑧 > 7 strong Ly 𝛼 emitters. One Ly 𝛼 emitter lies close to a strong photometric overdensity while the other shows no significant nearby overdensity, perhaps implying that not all strong 𝑧 > 7 Ly 𝛼 emitters reside in large ionized bubbles.  \nKey words: galaxies: high-redshift - galaxies: evolution - dark ages, reionization, first stars  \n★  \nE-mail: ryan.endsley@austin.utexas.edu', '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 𝛼 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 𝑧 = 5 . 3 (Bosman et al. 2022), approximately 1.1 Gyr after the Big Bang. Additional quasar and galaxy spectra at 𝑧 > 6, along with measurements of the cosmic microwave background, imply that reionization was about halfway complete ∼ 400 Myr earlier at 𝑧 ∼ 7 -8 (e.g. Davies et al. 2018; Mason et al. 2018; Planck Collaboration et al. 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 𝑧 ≳ 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).  \nPrior to JWST , over 1000 Lyman-break galaxies at 𝑧 > 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 et al. 2016; Livermore et al. 2017; Ishigaki et al. 2018; Oesch et al. 2018; Salmon et al. 2020) as well as wide-area ground-based imaging (e.g. Bowler et al. 2014, 2020; Ono et al. 2018; Stefanon et al. 2017, 2019; Endsley et al. 2021a; Harikane et al. 2022; Donnan et al. 2023). These studies demonstrated that the 𝑧 ∼ 6 -10 galaxy UV continuum luminosity functions had very steep faint-end slopes ( 𝛼 ∼ -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.  \nThe very high relative abundance of the faintest (M UV > -18) 𝑧 ≳ 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 ≲ 𝛽 ≲ -2 where 𝑓 𝜆 ∝ 𝜆 𝛽 ; e.g. Dunlop et al. 2012; Finkelstein et al. 2012; Rogers et al. 2013; Bouwens et al. 2014; Bhatawdekar & Conselice 2021) as well as very compact rest-UV morphologies (effective radii 𝑟 𝑒 ≲ 200 pc; Shibuya et al. 2015; Bouwens et al. 2022a; Neufeld et al. 2022). Based on results from local ( 𝑧 ∼ 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 𝑧 ≳ 6 population have been precluded by the lack of deep data probing the rest-frame optical portion of their spectral energy distributions (SEDs).  \nFor much of the past 20 years, 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 𝜎 imaging sensitivity of 𝑚 ∼ 26 . 5 (M UV ≲ -20 . 5 at 𝑧 ∼ 7) in the deepest pointings (120 hours; 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 ≲ -20) Lymanbreak 𝑧 ≳ 6 galaxies (e.g. Laporte et al. 2021; Stefanon et al. 2021b; Endsley et al. 2021a; Castellano et al. 2022; Schouws et al. 2022;  \nSommovigo et al. 2022; Tacchella et al. 2022; Topping et al. 2022a; Witstok et al. 2022; Whitler et al. 2023c), 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).  \nWithin the past year alone, data from JWST have greatly advanced our understanding of reionization-era galaxies. Rest-optical spectra fromthe Near-Infrared Spectrograph (NIRSpec; Jakobsen et al. 2022; Böker et al. 2023) have generally revealed very strong nebular line emission (e.g. [OIII], H 𝛼 , H 𝛽 ) coupled with signatures of hard radiation fields ([OIII]/[OII] > 10) and low gas-phase metallicities ( ≲ 0.1 𝑍 ⊙ ) among 𝑧 > 6galaxies (e.g. Arrabal Haro et al. 2023; Bunker et al. 2023; Cameron et al. 2023; Curti et al. 2023; Fujimoto et al. 2023; Larson et al. 2023; Sanders et al. 2023; Saxena et al. 2023; Tacchella et al. 2023b; Tang et al. 2023). 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 𝑧 > 5 galaxies are also emerging from early NIRSpec observations, including sources with (sometimes tentative) contributions from active galactic nuclei (AGN; Harikane et al. 2023; Larson et al. 2023; Kocevski et al. 2023; Maiolino et al. 2023a; Onoue et al. 2023; Übler et al. 2023) as well as objects that are in a very inactive stage of star formation (Looser et al. 2023a; Strait et al. 2023).  \nThe 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 Lymanbreak 𝑧 ≳ 6 selections with imaging are very efficient at yielding highly complete samples of objects with faint continua. Much of the early NIRCam studies on 𝑧 ≳ 6 galaxies focused on data from the Early Release Science (ERS) CEERS (Finkelstein et al. 2023) and GLASS (Treu et al. 2022) surveys which immediately pushed down to extremely deep Hubble imaging depths ( 𝑚 AB ∼ 29 mag at 5 𝜎 ) across 1-5 𝜇 m with just ∼ 3 hours 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 ≲ -19) 𝑧 ≳ 6 galaxies including their ages, stellar masses, nebular line strengths, UV slopes, dust attenuation strengths, and morphologies (e.g. Endsley et al. 2023; Topping et al. 2022b; Barrufet et al. 2023; Chen et al. 2023; Dressler et al. 2023b; Furtak et al. 2023; Hsiao et al. 2023; Labbé et al. 2023; Laporte et al. 2023; Leethochawalit et al. 2023; Rodighiero et al. 2023; Treu et al. 2023; Whitler et al. 2023b). 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) 𝑧 ≳ 6 galaxies has yet to be undertaken. Such an endeavor 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.  \nHere we take steps to address this shortcoming by utilizing NIRCamimagingtakenaspartoftheJWSTAdvancedDeepExtragalactic 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 hours 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  \net 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 ( 𝑚 AB ∼ 30 -31 mag 5 𝜎 ) NIRCam imaging in several bands from ∼ 1-5 𝜇 m (often including two medium-band filters) across > 200 arcmin 2 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; Dressler et al. 2023a; Hainline et al. 2023; Helton et al. 2023; Robertson et al. 2023; Tacchella et al. 2023a).  \nThe 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) 𝑧 ∼ 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 of reionization-era galaxies, complementing the JADES/NIRCam investigation in Dressler et al. (2023a) which utilized broadband 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) timescales.  \nThe structure of this paper is as follows. We begin by describing the imaging data, source extraction, photometric calculations, sample selection, and photo-ionization modelling in §2. Next, we present and discuss the UV luminosities, photometric redshifts, stellar masses, and light-weighted ages among our sample (§3). We then proceed to utilize photometric constraints on strong rest-optical nebular lines ([OIII]+H 𝛽 and H 𝛼 ; §4) to statistically analyze the star-forming and ionizing properties of relatively bright and very faint reionizationera galaxies (§5). Finally, we quantify the photometric overdensities surrounding two strong Ly 𝛼 emitters at 𝑧 > 7 within the JADES footprints to improve our understanding of the connection between strong Ly 𝛼 and ionized bubbles deep in the epoch of reionization (§6). Our main conclusions are summarized in §7.  \nThroughout 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 𝑀 ⊙ , and adopt a flat Λ CDM cosmology with parameters ℎ = 0 . 7, Ω M = 0 . 3, and Ω Λ = 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 (§2.1). We next detail our procedure for source extraction and photometric measurements in §2.2, then describe our Lyman-break selections of 𝑧 ∼ 6 and 𝑧 ∼ 7 -9 galaxies (§2.3). Lastly, we describe the photo-ionization SED models used throughout this work to infer the physical properties of the galaxies within our sample (§2.4).', '2.1 JWST and HST Imaging': "We consider all JADES/NIRCam imaging data taken prior to February 10, 2023 (see Rieke & the JADES Collaboration 2023). These data include imaging in the broad-band F090W, F115W, F150W, F200W, F277W, F356W, and F444W filters 1 as well as the F410M  \nmedium-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 Feb. 10, 2023. 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 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 2 ) and GOODS-N fields ( ≈ 46 arcmin 2 ).  \nTheNIRCamimagingreductionalgorithmappliedherelargelyfollows that described in previous works (e.g. Endsley et al. 2023; Merlin et al. 2022; Bezanson et al. 2022; Bagley et al. 2023b; Donnan et al. 2023; Rieke & the JADES Collaboration 2023; Robertson et al. 2023; Tacchella et al. 2023a; Williams et al. 2023a) 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 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. (2023b). From the output *\\_uncal.fits files, we subtract prominent artifacts (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 Feb. 1 2023) from JADES as well as the deep, blank-field public programs CEERS (Finkelstein et al. 2023), PRIMER (PI J. Dunlop), and NGDEEP (Bagley et al. 2023a). All snowball-masked and artifact-subtracted *\\_uncal.fits files are visually inspected; for a very small fraction of imaging exposures, we implement additional masking of remaining prominent artifacts 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 FRESCO (Oesch et al. 2023).  \nThe *\\_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 zeropoints from Boyer et al. (2022) that remain the most recent values implemented in CRDS. From the resulting calibrated *\\_cal.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 /s.pc/e.pc/p.pc package (Barbary 2016) and follow the methods described in Bagley et al. (2023b). The final *\\_cal.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. Using the /t.pc/w.pc/e.pc/a.pc/k.pc/r.pc/e.pc/g.pc software 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 /WFC3 F160W images, and then subsequently align all other NIRCam mosaics to that of F444W in the respective field.  \nTable 1. Summary of the 5 𝜎 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 𝜎 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 '' ).  \nFor our dropout selection of 𝑧 ∼ 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 (G. Brammer private communication) as described in Williams et al. (2023a). All HLF mosaics used here (v2.0) have the same pixel scale as the final NIRCam mosaics (30 mas/pixel). 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.  \nTo obtain reliable photometric colors 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 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 empirical PSFs (one per band and field) and creating the convolution kernels using /p.pc/h.pc/o.pc/t.pc/u.pc/t.pc/i.pc/l.pc/s.pc (Bradley et al. 2022).", '2.2 Source Extraction and Photometry': "To identify objects across the JADES footprints, we run S/o.pc/u.pc/r.pc/c.pc/e.pc E/x.pc/t.pc/r.pc/a.pc/c.pc/t.pc/o.pc/r.pc (Bertin & Arnouts 1996) on an inverse-variance weighted stack of the PSF-convolved F200W, F277W, F335M, F356W, F410M, and F444W mosaics 4 . The photometry is computed in Kron (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  \nwithin elliptical apertures with a Kron (1980) factor of 𝑘 = 1 . 2, as these apertures have been shown to yield optimal 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 /s.pc/e.pc/p.pc) with similar background pixel flux variance as the source of interest.  \nThe photometry in the 𝑘 = 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 𝑘 = 2 . 5 apertures divided by that of the flux in the 𝑘 = 1 . 2 apertures, with these measurements performed on the inverse-variance weighted stack used for the /s.pc/o.pc/u.pc/r.pc/c.pc/e.pc /e.pc/x.pc/t.pc/r.pc/a.pc/c.pc/t.pc/o.pc/r.pc detection. The second stage is correcting for flux outside the 𝑘 = 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 𝑧 ∼ 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 catalog.  \nBecause of NIRCam's sensitivity, a significant number of our identified Lyman-break dropout 𝑧 ∼ 6 -9 galaxies have 𝑘 = 2 . 5 apertures that contain one or more different objects, artificially boosting their recovered total fluxes. We therefore utilize the neighbor subtraction algorithm described in Endsley et al. (2023) to subtract off any neighboring object within the 𝑘 = 2 . 5 aperture (as determined by the /s.pc/o.pc/u.pc/r.pc/c.pc/e.pc /e.pc/x.pc/t.pc/r.pc/a.pc/c.pc/t.pc/o.pc/r.pc segmentation map) and recompute the photometry on those images. Moreover, because the morphologies of 𝑧 ≳ 6 galaxies are sometimes clumpy (e.g. Chen et al. 2023; Hainline et al. 2023; Treu et al. 2023), the galaxies we intend to study can be identified as multiple separate objects by /s.pc/o.pc/u.pc/r.pc/c.pc/e.pc /e.pc/x.pc/t.pc/r.pc/a.pc/c.pc/t.pc/o.pc/r.pc. If two or more nearby objects satisfy one of our Lyman-break dropout selections ( 𝑧 ∼ 6 or 𝑧 ∼ 7 -9; see below) and have overlapping 𝑘 = 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.  \nIn Table 1, we summarize the range of 5 𝜎 depths in each HST and NIRCam band among our final Lyman-break galaxy sample of F775W and F090W dropouts (see §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 𝜎 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 𝜎 depth of our final 𝑧 ∼ 6 -9 galaxy sample is 𝑚 AB ≈ 28 . 5 in F435W, F606W, F775W, and F814W, though ≈ 0.7 mag shallower in F850LP ( 𝑚 AB = 27 . 8). In NIRCam, the typical 5 𝜎 depths range from 𝑚 AB ≈ 29 . 0 in F410M and F444W to 𝑚 AB ≈ 29 . 7 in F277W and F356W. The deepest imaging in GOODS-S reaches 5 𝜎 depths of 𝑚 AB ≈ 30 . 2 in several ACS and NIRCam bands for compact objects.", '2.3 Selection of Lyman-break z ∼ 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 𝑧 ∼ 6 -9  \nwhere 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 𝑧 ∼ 9 is chosen such that the NIRCam SEDs remain sensitive to the [OIII]+H 𝛽 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).  \nFollowing 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 𝑧 ∼ 6 -9 galaxies via a Lyman-break color selection. These Lyman-break selections rely on the appearance of a sharp spectral discontinuity at 𝜆 rest = 1216 Å that is imposed by very strong absorption blueward of the Ly 𝛼 line from HI 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. (2023) which aimed at establishing a census of 𝑧 > 8 galaxies within the JADES/NIRCam data set and providing a first discussion of their colors 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 color 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%) of lying at high redshift ( 𝑧 > 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).  \nBecause the Ly 𝛼 break spans a considerable range in observed wavelength across our targeted 𝑧 ∼ 6 -9 interval, we divide our selection according to two separate sets of criteria, one for 𝑧 ∼ 6 and one for 𝑧 ∼ 7 -9, which we discuss in sequence below. At 𝑧 ∼ 6, the Ly 𝛼 break is located at ≈ 0.85 𝜇 m and thus galaxies at this redshift should appear as strong ACS/F775W dropouts (e.g. Stanway et al. 2003; Bunker et al. 2004; Bouwens et al. 2004). We therefore begin our 𝑧 ∼ 6 galaxy selection with the following color cuts:  \n- (i) F775W -F090W > 1.2\n- (ii) F090W -F150W < 1.0\n- (iii) F775W -F090W > F090W -F150W + 1.2.  \nThese criteria enforce the presence of a strong break between the F775WandF090Wbands(at least 1.2 mags) as well as a much flatter color between F090W and F150W. The second color cut above sets an upper redshift limit of 𝑧 ≈ 6 . 5 where the Ly 𝛼 break redshifts well into F090W. In these cuts, the F775W flux is set to its 1 𝜎 upper limit in cases where the S/N < 1 in that band following past Lymanbreak color cut selections in the literature (e.g. Bouwens et al. 2015; Endsley et al. 2021a).  \nAdditional analytic cuts are imposed to verify the presence of a strong Ly 𝛼 break. First, we enforce S/N < 2 in ACS/F435W since this band lies fully blueward of the Lyman-continuum break at 𝑧 > 5 . 5 and all Lyman-continuum 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 𝜎 upper limit when S/N < 1. However, we ignore the cut in F435W S/N as well as the F606W dropout criteria  \nfor objects with extremely strong Ly 𝛼 breaks (F775W -F090W > 2.5) which can confidently be identified as 𝑧 ∼ 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 𝑧 ∼ 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.  \nOur 𝑧 ∼ 7 -9 Lyman-break selection largely follows the same logic of the 𝑧 ∼ 6 criteria described above. Specifically, we begin with the color cuts  \n- (i) F090W -F115W > 1.5\n- (ii) F115W -F200W < 1.2\n- (iii) F090W -F115W > F115W -F200W + 1.5.  \nThese cuts are satisfied by galaxies at 𝑧 ≈ 7 . 0 -9 . 0 with blue restUV colors ( -2 . 5 ≲ 𝛽 ≲ -2 . 0), as well as extremely red ( -1 . 0 ≲ 𝛽 ≲ 0 . 0) objects at 𝑧 ≈ 7 -8 where the Lyman-alpha break has not yet shifted into F115W. As with the 𝑧 ∼ 6 criteria, we also impose the condition that S/N < 2 in F435W. For the 𝑧 ∼ 7 -9 selection, we additionally utilize the 𝜒 2 opt parameter defined in Bouwens et al. (2015) as 𝜒 2 opt ≡ " abs ( 𝑓 )/ 𝑓 ( 𝑓 / 𝜎 ) 2 where 𝑓 and 𝜎 represent, respectively, the measured flux density and its error in a given band, while abs ( 𝑓 ) is the absolute value of the flux density. After summing over the ACS F435W, F606W, and F775W bands, we enforce 𝜒 2 opt < 5. However, we ignore the S/N cut in F435W as well as the 𝜒 2 opt cut for objects with extremely strong Ly 𝛼 breaks (F090W -F115W > 2.5). Every selected 𝑧 ∼ 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.  \nWe enforce a final set of cuts to both the 𝑧 ∼ 6 and 𝑧 ∼ 7 -9 samples to ensure that the NIRCam data are sufficiently sensitive for rest-optical color measurements (tracing e.g. Balmer breaks or nebular line EWs). Following the IRAC-based approach of Endsley et al. (2021a), we require that 𝑓 ( FUV )/ 𝑒 ( X ) > 3 where 𝑓 ( FUV ) is the far-UV flux density corresponding to the inferred M UV value from the constant star formation history SED fits (see §2.4) while 𝑒 ( X ) is the 1 𝜎 uncertainty in the flux density of band 𝑋 . Here, this cut is enforced with the four reddest bands in the JADES filter set (i.e. 𝑋 ∈ {F335M, F356W, F410M, F444W}).  \nEvery selected 𝑧 ∼ 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 colors consistent with brown dwarfs, which can mimic a strong F090W dropout. The measured colors are compared to the empirical SPEX 0.8-2.5 𝜇 mbrown 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 neighbor subtraction and combined aperture photometry for the appropriate set of objects (see §2.2), we end up with a final sample of 280 F775W dropouts ( 𝑧 ∼ 6) and 480 F090W dropouts ( 𝑧 ∼ 7 -9) across the ≈ 90 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 𝑧 ∼ 6 -9 galaxy in our JADES sample, we fit their 14-band ACS+NIRCam photometry with star-forming photoionization models. One of the  \nprimary 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; Endsley et al. 2023; Tacchella et al. 2022; Topping et al. 2022a; Whitler et al. 2023c,b), stellar mass estimates can change significantly depending on the assumed star formation history (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.  \nFor all four sets of SED model fits described below, we adopt a Chabrier (2003) stellar initial mass function (IMF) with bounds 0.1-300 𝑀 ⊙ as well as an SMC dust attenuation curve (Pei 1992). Moreover, in each case, we fit with log-uniform priors on stellar mass in the range 5 ≤ log ( 𝑀 ∗ / 𝑀 ⊙ ) ≤ 12, on V-band dust optical depth in the range -3 ≤ log 𝜏 V ≤ 0 . 7, on ionization parameter in the range -4 ≤ log 𝑈 ≤ -1, and on metallicity in the range -2 . 2 ≤ log( 𝑍 / 𝑍 ⊙ ) ≤ -0 . 3. The upper limit of ≈ 50% 𝑍 ⊙ is set to avoid unphysically high metallicities for the faint (M UV ≳ -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 e.g. a stellar mass to metallicity relation prior or 𝛼 -enhanced metallicity models for future work.  \nAll F775W and F090W dropouts are fit in the redshift range 𝑧 = 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 𝑧 ∼ 6 -9 objects with strong nebular lines (and hence young light-weighted ages) given the unique NIRCam color patterns caused by [OIII]+H 𝛽 and H 𝛼 lines at high redshifts.  \n/b.pc/e.pc/a.pc/g.pc/l.pc/e.pc Constant SFH models: We first consider models adopting a constant SFH (hereafter CSFH) using the BayEsian Analysis of GaLaxy sEds (/b.pc/e.pc/a.pc/g.pc/l.pc/e.pc) SED-fitting code (Chevallard & Charlot 2016). /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc 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 (/p.pc/a.pc/r.pc/s.pc/e.pc/c.pc; Bressan et al. 2012; Chen et al. 2015). The Bayesian /m.pc/u.pc/l.pc/t.pc/i.pc/n.pc/e.pc/s.pc/t.pc algorithm (Feroz & Hobson 2008; Feroz et al. 2009) is implemented in /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc to determine the posterior probability distribution for each free physical parameter in the fit. In the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc 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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH ages as the light-weighted ages of the galaxies.  \n/b.pc/e.pc/a.pc/g.pc/l.pc/e.pc Two-component SFH models: We allow for more flexible SFH models within context of /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc 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 𝜏 model component (SFR ∝ 𝑡 𝑒 -𝑡 / 𝜏 SF where 𝑡 is the time since the onset of star formation) and 2) a constant SFH component. The onset of the delayed 𝜏 model component is assumed to have been between 10 1 . 35 ( ≈ 22) Myr age up to the age of the Universe at the fitted redshift, while the constant SFH 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 constant SFH 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 𝜏 component, the 𝜏 SF parameter is fit with a log-uniform prior in the range 1 Myr to 30 Gyr. The SFR of the  \nconstant SFH component is determined via the specific star formation rate (sSFR) over the associated time interval, where we adopt a log-uniform prior of -14 ≤ log(sSFR/yr -1 ) ≤ -6.  \n/p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior models: We also consider non-parametric SFH models using the SED-fitting code /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc (Leja et al. 2019; Johnson et al. 2021b). /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc adopts the star-forming photoionization models from /f.pc/s.pc/p.pc/s.pc (Conroy et al. 2009; Conroy & Gunn 2010) which by default utilizes products from the MESA Isochrones and Stellar Tracks project (/m.pc/i.pc/s.pc/t.pc; Choi et al. 2016). Posterior probability distributions on fitted parameters are determined using the /d.pc/y.pc/n.pc/e.pc/s.pc/t.pc/y.pc sampling package (Speagle 2020). We first consider the 'continuity' prior implemented in the nonparametric SFH models in /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc. 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.  \nThe implementation of the non-parametric SFH models in this work largely follows that described in Whitler et al. (2023c). To summarize, the SFHs are composed of 8 time bins where the SFR in each is a constant value, and the ratios of the SFR in adjacent time bins are fit by /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc. The earliest time bin extends to a fitted formation redshift in the range 𝑧 form = 10 -30 (uniform prior), and the two most recent time bins are fixed to 0-3 Myr and 3-10 Myr while the remaining 6 time bins are divided evenly in logarithmic space to the fitted formation redshift.  \n/p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc 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 star formation histories like the continuity prior.  \nFor 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 ± 1 𝜎 error as the inner 68% 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.  \nPrior to running the SED fits, we add a 5% 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 𝑧 ∼ 6 -9 Lyman-break galaxies in our sample. For each of the four models sets described above, the median best-fitting 𝜒 2 value over the full sample ranges between 10-14 (14 photometric data points are fit in every case). Moreover, 90-98% of the SED fits yield best-fitting 𝜒 2 < 28.  \nHowever, 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 𝜒 2 > 100 with the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH models and best-fitting 𝜒 2 > 50 for the other three models. Two of these objects appear to have steeply rising SEDs between F277W and F444W yet relatively blue colors between F115W and F277W, reminiscent of the 'Little Red Dot' (Matthee et al. 2023) 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. 2023; Williams et al. 2023b; Pérez-González et al. 2024). Both of these systems have been explored in Williams et al. (2023b, 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 multi-band MIRI imaging. A third object selected in  \nFigure 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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc 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 𝑧 ≈ 8 . 6, our selection criteria is sensitive to blue galaxies at redshifts as high as 𝑧 ≈ 9 . 0 where the [OIII]+H 𝛽 emission lines remain in F444W (see §2.3).  \n<!-- image -->  \nour Lyman-break sample yet very poorly fit with star-forming only models shows very broad H 𝛼 emission and is a candidate for hosting a black hole merger (ID 10013704 in Maiolino et al. 2023b). Because the goal in this paper is to study reionization-era galaxies in context of star-forming models, we remove these four poorly-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 𝑧 ∼ 6 -9 galaxies analyzed below consists of 278 F775W dropouts ( 𝑧 ∼ 6) and 478 F090W dropouts ( 𝑧 ∼ 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 stellar masses (§3.1) and light-weighted ages (§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.  \nThe F775W and F090W dropout subsets are largely comprised of galaxies with photometric redshifts in the range 5 . 5 < 𝑧 phot < 6 . 5 and 6 . 5 < 𝑧 phot < 8 . 5, respectively (see Fig. 1), consistent with expectations given our Lyman-break selection criteria. The mean photometric redshift is ⟨ 𝑧 phot ⟩ = 5 . 9 and 7 . 3 for the F775W and F090W dropout subsets, respectively. Considering all 756 of these 𝑧 ∼ 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 𝐿 ∗ UV = -20 . 5 (e.g. Bowler et al. 2017; Harikane et al. 2022), the M UV range of our sample corresponds to ≈ 0.02 - 4 𝐿 ∗ UV thus covering a factor of ≈ 200 in UV luminosity. The  \naverage absolute UV magnitude of the F775W dropout ( 𝑧 ∼ 6) subset is ⟨ M UV ⟩ = -18 . 6, which is slightly brighter than that of the F090W dropout subset ( ⟨ M UV ⟩ = -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.91.0 mag) and each include a handful of galaxies at very bright UV luminosities (M UV < -21 . 25 or 𝐿 UV > 2 𝐿 ∗ UV ) as well as ≈ 100-200 galaxies at the very faint end (M UV > -18 or 𝐿 UV < 0 . 1 𝐿 ∗ UV ; see Fig. 1).', '3.1 The Stellar Masses of Lyman-break z ∼ 6 -9 Galaxies': 'The756Lyman-break 𝑧 ∼ 6 -9 galaxies comprising our JADES sample possess a wide range of inferred stellar masses. From the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH models, we infer stellar masses spanning approximately 3.5 orders of magnitude from 𝑀 ∗ ≈ 2 × 10 6 𝑀 ⊙ up to 𝑀 ∗ ≈ 7 × 10 9 𝑀 ⊙ (see Fig. 2a,b). However, we find that the inferred stellar masses of individual 𝑧 ∼ 6 -9 galaxies can change substantially depending on the assumedSFH,consistent with the results of previous works (e.g. Endsley et al. 2023; Tacchella et al. 2022; Topping et al. 2022a; Whitler et al. 2023c,b). In our sample the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH stellar masses are, on average, approximately 0.2, 0.1, and 0.5 dex smaller than (i.e. 0.6 × , 0.8 × , and 0.3 × ) the stellar masses inferred from the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc two-component SFH fits, the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc bursty SFH prior fits, and the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior fits, respectively (see Fig. 2eg). Therefore, the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH fits tend to yield the lowest stellar masses among our sample while the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc 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 ( ∼ 50 Myr; see §3.2) and in these models no stars are assumed to have formed prior to that inferred age. On the other hand, the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity prior imposes a preference for significant star formation extending back to 𝑧 = 10 -30 ( Δ 𝑡 ∼ 300 -650 Myr for a galaxy at 𝑧 ∼ 7) by weighting against strong time-variability in the SFH.  \nFor galaxies with the youngest light-weighted ages in our sample (age CSFH < 10Myr), the stellar masses inferred with the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior are typically 0.8 dex (i.e. 6 × ) higher than that inferred with the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH setup, though this difference can rise to factors of ∼ 30-100 in the most extreme cases (see Fig. 2g) as found in Whitler et al. (2023c) using a sample of very bright (M UV ≲ -21) 𝑧 ∼ 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 (age CSFH > 100 Myr), with a median and maximum offset of approximately 0.3 dex (2 × ) and 0.7 dex (5 × ), respectively (see Fig. 2g), again consistent with the IRAC-based findings of Whitler et al. (2023c). Consequently, we find that the range of stellar masses for our sample shifts upwards to 𝑀 ∗ ≈ 7 × 10 6 𝑀 ⊙ -2 × 10 10 𝑀 ⊙ with the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFHprior fits (see Fig. 2c,d), though we continue to conclude that the faintest (M UV ∼ -17) and youngest (age CSFH ≲ 30 Myr) galaxies tend to have the lowest stellar masses of ∼ 10 7 𝑀 ⊙ .  \nWe quantify the relationship between UV luminosity and stellar mass among our sample by fitting the equation log ( 𝑀 ∗ / 𝑀 ⊙ ) = 𝑎 GLYPH<16> M UV + 19 GLYPH<17> + 𝑏 . This relation is derived at 𝑧 ∼ 6 and 𝑧 ∼ 7 -9 by fitting the F775W and F090W dropout subsets separately, and we moreover perform the fits using the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH and /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior models to bracket the systematic uncertainties associated from different SFH assumptions. In these 𝑀 ∗ -M UV fits, we only include objects with relatively bright UV luminosities (M UV ≤ -18) as we aim to mitigate potentially significant biases  \n<!-- image -->  \nFigure 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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc constant SFH (left panels) and /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior fits (see §2.4). These two mass estimates generally bracket the systematic uncertainties associated with the assumed SFH and we show the fitted relation (black line) for each adopting the form log( 𝑀 ∗/ 𝑀 ⊙ ) = 𝑎 (M UV + 19) + 𝑏 and restricting the fit to M UV ≤ -18 (the dashed line shows the relation extrapolated to fainter UV luminosities). Each point is color-coded by the inferred CSFH age of the associated galaxy. Right: Illustration of systematic differences in stellar mass estimates from adopting different assumptions on the SFH. We compare stellar mass estimates from the two-component SFH (TcSFH) /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc fits (top), the bursty SFH prior /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc fits (middle), and the continuity SFH prior /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc fits (bottom), all relative to the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH stellar masses as a function of CSFH age.  \n<!-- image -->  \nat the faintest end, though acknowledge that future work directly accounting for incompleteness will be required. The fitted parameters ( 𝑎 and 𝑏 ) 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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH outputs, again keeping only those with M UV ≤ -18. With this simple approach, we derive the 𝑀 ∗ -M UV relations shown in panels Fig. 2a-d.  \nWhenadopting the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH models, we infer a typical stellar mass of 𝑀 ∗ ≈ 2 . 1 × 10 8 𝑀 ⊙ and ≈ 3 . 4 × 10 7 𝑀 ⊙ for 𝑧 ∼ 6 galaxies with M UV = -20 and M UV = -18, respectively. At 𝑧 ∼ 7 -9, the typical inferred CSFH stellar masses are lower by ≈ 0.15-0.2 dex at these M UV values, resulting in values of 𝑀 ∗ ≈ 1 . 5 × 10 8 𝑀 ⊙ and ≈ 2 . 1 × 10 7 𝑀 ⊙ , respectively. This redshift trend is consistent with previous findings that, at fixed M UV , typical inferred stellar masses decrease slightly with increasing redshift at 𝑧 ∼ 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).  \nThe normalization of our derived 𝑀 ∗ -M UV relation is ≈ 0.5 dex (i.e. ≈ 3 × ) higher when adopting the stellar masses from the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior fits, as expected from our discussion above. The slope of the 𝑀 ∗ -M UV relation we derive in each case ( -0 . 48 ≤ 𝑎 ≤ -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 UV than  \nfound 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 𝑧 ∼ 7 -9 CSFH 𝑀 ∗ -M UV 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 𝑧 ∼ 6 -9 Galaxies in JADES': 'Shortly after the first JWST /NIRCam images were released, there were reports of very high-redshift ( 𝑧 ∼ 7 -11) candidates with extremely large stellar masses of 𝑀 ∗ ∼ 10 11 𝑀 ⊙ identified over small areas ( ≈ 40 arcmin 2 ; Labbé et al. 2023) which clearly challenged models of galaxy formation (e.g. Boylan-Kolchin 2023; Lovell et al. 2023). Since these initial reports 5 , stellar mass estimates for nearly all the candidates have been lowered to < 3 × 10 10 𝑀 ⊙ , though one candidate 𝑀 ∗ = 10 10 . 9 𝑀 ⊙ 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) 𝜒 2 (see Endsley et al. 2023). Here, we utilize the deep 9-band NIRCam imaging of JADES (including two longwavelength medium bands) to build upon these initial studies by  \nGalaxies with the Largest Stellar Masses (see also Fig. 4)  \nFigure 3. The ACS+NIRCam SEDs for 12 of the 13 galaxies with the largest inferred stellar masses in our sample ( 𝑀 ∗ ≥ 10 9 . 5 𝑀 ⊙ from the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior fits; the last object is shown in Fig. 4). The measured photometry is shown with blue diamonds while the model photometry from the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior SEDs is shown with red circles. The black lines and corresponding shaded regions show the median posterior from the SEDs and their 68% inner credible intervals, respectively, fixed at the median redshift from the posterior of each galaxy. The /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior fits generally yield the largest stellar masses (see Fig. 2) and we find a maximum stellar mass of 10 10 . 0 𝑀 ⊙ across our full sample with these fits (see Fig. 4).  \n<!-- image -->  \nFigure 4. The SED fit of the most massive 𝑧 ∼ 6 -9 Lyman-break galaxy in our sample, JADES-GS+53.03755 -27.87491. We show the fit both excluding (green) and including (red) the MIRI F770W photometry for this object, with both adopting the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc 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 𝑀 ⊙ across our full sample (see also JADES-GN+189.13794+62.23601 in Fig. 3).  \n<!-- image -->  \ninvestigating the potential abundance of extremely massive 𝑧 ≳ 6 galaxies.  \nThere are no galaxies in our sample where the photometric data clearly imply stellar masses of > 3 × 10 10 𝑀 ⊙ . Only a small number of galaxies have inferred stellar masses around ∼ 10 10 𝑀 ⊙ when applying the continuity SFH prior with /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc which, as discussed in §3.1, generally yields a maximum stellar mass estimate for a given galaxy (see Fig. 2). Here, we discuss these most massive 𝑧 ∼ 6 -9 candidates in our JADES sample.  \nThere are 13 galaxies with inferred stellar masses of ≥ 3 × 10 9 𝑀 ⊙ from the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior fits among our JADES sample. Eleven and two of these thirteen 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 NIRCam colors consistent with strong Balmer breaks implying that a relatively old ( ≳ 300 Myr) stellar population is contributing significantly to the emergent light. However, a few of these galaxies show colors consistent with relatively weak Balmer breaks and are simply so luminous ( -22 ≲ M UV ≲ -21 . 5) that a substantial population of relatively old stars can be hidden within the observed SED if outshined by more recently-formed stars.  \nThe galaxy with the highest inferred stellar mass in the sample (JADES-GS+53.03755 -27.87491; 𝑀 ∗ = 10 10 . 2 ± 0 . 1 𝑀 ⊙ ) lies at 𝑧 phot ≈ 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 [OIII]+H 𝛽 and H 𝛼 (Fig. 4). Fortuitously, this galaxy falls within the deep MIRI F770W parallel imaging of JADES (see  \nFigure 5. ACS+NIRCam SEDs for a subset of galaxies with fairly typical light-weighted ages among our sample (age CSFH ∼ 50 Myr). We show galaxies at bright (M UV ≤ -19 . 5; upper panels) and faint (M UV > -18; lower panels) UV luminosities as well as those selected as F775W dropouts ( 𝑧 ∼ 6; left panels) and F090W dropouts ( 𝑧 ∼ 7 -9; right panels). The model photometry from the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc constant star formation history (CSFH) fits are shown with red circles; all quoted properties in this figure are inferred from the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH models.  \n<!-- image -->  \nGalaxies with Typical Light-Weighted Ages (  \n∼  \n50 Myr)  \nEisenstein 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 near infrared (e.g. Tang et al. 2022; Papovich et al. 2023). With the MIRI measurement included (Alberts et al. 2024), we infer a consistent though slightly ( ≈ 0.2 dex) lower stellar mass of 𝑀 ∗ = 10 10 . 0 ± 0 . 1 𝑀 ⊙ (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 ( ≈ 0.3 𝑍 ⊙ ) over that without the MIRI data ( ≈ 0.05-0.1 𝑍 ⊙ ). Such higher metallicity solutions alter the light-to-mass ratio across the rest-frame ≈ 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 prep.; Helton et al. in prep; Ji et al. in prep).  \nAside from JADES-GS+53.03755 -27.87491, there is only one other galaxy in the 𝑧 ∼ 6 -9 sample considered here with a comparable stellar mass: JADES-GN+189.13794+62.23601 with an estimated 𝑀 ∗ = 10 10 ± 0 . 1 𝑀 ⊙ from the ACS+NIRCam data. JADESGN+189.13794+62.23601 does not possess any MIRI constraints in the rest-frame near-infrared 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 /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc fits) among all galaxies in our sample is 1.0 × 10 10 𝑀 ⊙ .  \nWe remind the reader that, in this work, we are only considering Lyman-break selected 𝑧 ∼ 6 -9 galaxies with photometry that can be reasonably explained by star-forming photoionization models. There are four Lyman-break 𝑧 ∼ 6 -9 objects in JADES with photometry  \nthat is very poorly matched by star-forming only models and hence have been excluded from the sample considered here (see §2.4). Notably, these include compact objects with NIRCam colors similar to the most massive candidate reported in Labbé et al. (2023) and others subsequently identified in more recent data sets (Furtak et al. 2022; Akins et al. 2023; Barro et al. 2023), including within JADES (Williams et al. 2023b; 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 ∼ 6 -9 Galaxies': 'In this sub-section, we discuss the light-weighted ages (defined here as CSFHages)inferred amongoursampleof756Lyman-breakgalaxies. 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 a 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 ∼ 50 Myr) consistent with previous findings from early JWST data (e.g. Endsley et al. 2023; Furtak et al. 2023; Laporte et al. 2023; Leethochawalit et al. 2023; Morishita & Stiavelli 2023; Whitler et al. 2023b). 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 [OIII]+H 𝛽 and/or H 𝛼 emission (EW ∼ 7001000 Å; see e.g. IDs JADES-GN+189.19124+62.19952 and JADESGS+53.14154 -27.82320 in Fig. 5) implying a substantial contribution of hot, massive stars to the SED. However, a subset of the young  \n(age CSFH ∼ 50 Myr) galaxies have measured photometry consistent with surprisingly weak nebular line emission (e.g. IDs JADESGN+189.19376+62.29430 and JADES-GS+53.15167 -27.80925 in Fig. 5); this sub-population is discussed further in §5.2.  \nIn addition to the typical young (age CSFH ∼ 50 Myr) galaxies, several sources in our sample show long-wavelength photometric excess patterns indicating extremely strong nebular line emission ([OIII]+H 𝛽 EW ∼ 2000 -5000 Å or H 𝛼 EW ∼ 1000 -2500 Å), consistent with exceptionally young light-weighted ages (age CSFH ∼ 3 Myr). Other galaxies show strong Balmer breaks with signatures of weak-to-no nebular line emission ([OIII]+H 𝛽 EW ≲ 300 Å), consistent with relatively old light-weighted ages (age CSFH ∼ 300 -1000 Myr). These sub-populations representing the extreme ends of the light-weighted age distribution are analyzed 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 et al. 2020). This degeneracy could be alleviated for galaxies situated in specific redshift intervals, leading to the early identification of a small number of 𝑧 ∼ 7 -9 galaxies showing colors 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 𝑧 ∼ 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 ≈ 3-5 𝜇 m, we identify and characterize the 𝑧 ∼ 6 -9 Lyman-break galaxies in our JADES data that confidently show significant Balmer breaks.  \nFor ease of comparison with existing literature (Binggeli et al. 2019), we quantify the strength of the Balmer break as 𝐹 𝜈 (4200 Å) / 𝐹 𝜈 (3500 Å) where 𝐹 𝜈 is the (continuum) flux density at the specified rest-frame wavelength. These flux densities are taken from the posterior SEDs from the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc fits and are corrected for dust using the inferred 𝐴 𝑉 value associated with each step of the nested sampling chain. Here we focus on the subset of galaxies with inferred 𝐹 𝜈 (4200 Å) / 𝐹 𝜈 (3500 Å) > 1.25 at > 84% probability from the posteriors of both the CSFH and two-component SFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc fits. These objects are classified as those with confident strong Balmer breaks throughout this work.  \nWeidentify 13 and 9 galaxies with confident strong Balmer breaks ( 𝐹 𝜈 (4200 Å) / 𝐹 𝜈 (3500 Å) > 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 reionizationera 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 ≤ M UV ≤ -17 . 1 and redshifts of 𝑧 phot = 5 . 1 -7 . 9. While we do not identify any objects with confident strong Balmer breaks at the brightest UV luminosities ( -22 ≲ M UV ≲ -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 ( ≈ 90 arcmin 2 ). IRAC studies covering  \n> deg 2 fields have revealed 𝑧 ∼ 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. 2023c).  \nThe JADES galaxies with confident strong Balmer breaks are among those with the oldest light-weighted ages in the sample with CSFH ages spanning age CSFH ≈ 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 [OIII]+H 𝛽 EWs ≈ 300-400 Å (e.g. IDs JADESGN+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 10 Myr ≈ 1 -10 Gyr -1 with all four SED fitting procedures described in §2.4. However, we also identify a number of galaxies with strong Balmer breaks showing colors consistent with very weak to no rest-optical line emission ([OIII]+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 10 Myr ≲ 0 . 3 Gyr -1 ) due to a declining SFH. The CSFH fits are forced to sSFR 10 Myr > 1 Gyr -1 solutions given the age of the Universe at 𝑧 > 6, and thus always have some modest line emission in the models (minimum [OIII]+H 𝛽 EW ∼ 200 Å). Due to the very minor effect such weak emission lines can have on the observed broadband 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 𝑧 = 5 . 2 -7 . 3 (Looser et al. 2023a; Strait et al. 2023). 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 𝑧 > 4 galaxies typically exhibit high EW ( ≳ 500 Å) rest-optical nebular emission lines (i.e. [OIII]+H 𝛽 or H 𝛼 ) from their impact on broadband Spitzer /IRAC colors (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 ≲ -20) 𝑧 ∼ 7 -8 Lyman-break galaxies exhibited very high EW [OIII]+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 𝑧 ≳ 6 extreme line emitters at fainter continuum luminosities, particularly when the nebular lines are situated in a medium band (e.g. Endsley et al. 2023; Bouwens et al. 2023; Rinaldi et al. 2023; Withers et al. 2023). In this sub-section, we build upon these previous works by utilizing the deep 9-band NIRCam imaging from JADES (including F335M and F410M) to identify and characterize a statistical sample of 𝑧 ∼ 6 -9 galaxies which confidently exhibit extreme rest-optical line emission.  \nAcross our full sample of 756 𝑧 ∼ 6 -9 Lyman-break GOODS+JADES galaxies, we identify several systems with longwavelength NIRCam color patterns indicating extremely high EW rest-optical emission lines. For the purpose of explicitly investigat-  \nFigure 6. SEDs for a subset of galaxies that confidently show strong Balmer breaks ( 𝐹 𝜈 (4200 Å) / 𝐹 𝜈 (3500 Å) > 1.25). All quoted properties are inferred from the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH fits, though the inferred Balmer breaks remain strong at the > 84% credible interval with the two-component SFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc fits. The figure format is very similar to Fig. 5.  \n<!-- image -->  \nGalaxies with the Strongest Balmer Breaks  \nGalaxies with Extreme Rest-optical Line Emission  \nFigure 7. SEDs for a subset of galaxies that confidently exhibit extreme rest-optical line emission ([OIII]+H 𝛽 EW > 1500 Å and/or H 𝛼 EW > 800 Å). All quoted properties are inferred from the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH fits, though the optical line EWs remain extreme at the > 84% credible interval with the two-component SFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc fits. The figure format is very similar to Fig. 5 and inferred H 𝛼 EWs are only quoted for objects at 𝑧 phot < 6 . 5 where NIRCam is sensitive to this line.  \n<!-- image -->  \ning this population, we restrict our attention to galaxies where the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc posteriors yield > 84% probability of extremely high EWs from both the constant SFH and two-component SFH fits. By imposing this probability cut on the results from two SFH models we  \nbetter ensure that the inferred presence of extreme line emission is likely not being confused with a very strong Balmer break.  \nThere are 20 and 29 galaxies that confidently exhibit extreme [OIII]+H 𝛽 emission (EW > 1500 Å) within our F775W and F090W dropout samples, respectively. We show the SEDs of a subset of  \nthese systems in Fig. 7, where objects were chosen to reflect a range of UV magnitudes and redshifts. The confident extreme [OIII]+H 𝛽 emitters in GOODS+JADES have absolute UV magnitudes spanning -21 . 0 ≤ M UV ≤ -16 . 9 and thus encompass nearly the entire M UV range of our full 𝑧 ∼ 6 -9 sample. Remarkably, this includes galaxies in the very UV-faint regime (M UV > -18; see e.g. JADESGS+53.17835 -27.77879 and JADES-GS+53.19590 -27.79240 in Fig. 7). The vast majority of these extreme [OIII]+H 𝛽 emitters show color patterns indicating they lie at redshifts consistent with the targeted windows of our dropout selections ( 𝑧 phot ≈ 5 . 5 -8 . 4). Only two outlier sources show strong excesses in F277W suggesting 𝑧 phot ∼ 5 which presumably entered our 𝑧 ∼ 6 Lyman-break selection due to photometric noise in F775W and/or F090W.  \nThe confident extreme [OIII]+H 𝛽 emitters all show SEDs consistent with exceptionally young light-weighted ages (age CSFH ≈ 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 continua (see e.g. IDs JADESGS+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 [OIII]+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.  \nThe very young light-weighted ages indicate that these extreme [OIII]+H 𝛽 emitters have recently experienced a dramatic rise in star formation rate (e.g. Tang et al. 2019, 2022; Endsley et al. 2021a; Carnall et al. 2023; Tacchella et al. 2023b). In context of a constant star formation history, we infer an extremely large typical sSFR 10 Myr ∼ 300 Gyr -1 for these extreme [OIII]+H 𝛽 emitters, which is simply the inverse of their typical light-weighted age (age CSFH ∼ 3 Myr). But even with the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity SFH prior fits (which generally yield the lowest sSFR estimates; see §3), we continue to infer a very high median sSFR averaged over the past 10 Myr for these systems (26 Gyr -1 ) with values up to 72 Gyr -1 .  \nAt 𝑧 ∼ 6, the JADES/NIRCam photometry provides constraints on not only the EW of [OIII]+H 𝛽 emission, but also H 𝛼 . We identify 31 galaxies within our F775W dropout sub-sample that have a > 84% posterior probability of H 𝛼 EW > 800 Å and 𝑧 < 6 . 5 from both the CSFH and two-component /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc SFH fits. Here, we have added the 𝑧 < 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 [OIII]+H 𝛽 emitters (20 in both subsets) and thus has similar characteristics. The confident H 𝛼 emitters span absolute UV magnitudes of -21 . 0 ≤ M UV ≤ -17 . 4 and light-weighted ages of age CSFH ≈ 1-10 Myr. A subset show colors implying prominent Balmer jumps in their continuum and have inferred H 𝛼 EWs confidently reaching values up to ≈ 20003000 Å (see e.g. JADES-GS+53.16900 -27.80079 in Fig. 7). Their typical specific star formation rates range from sSFR 10 Myr = 150 Gyr -1 with the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH fits down to sSFR 10 Myr = 14 Gyr -1 with the /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc continuity prior SFH models. We place both the extreme [OIII]+H 𝛽 and extreme H 𝛼 emitter populations in context of the broader sample in the following section.', '4 THE DEPENDENCE OF [OIII]+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 ≲ -20) 𝑧 ∼ 7 systems commonly show far higher [OIII]+H 𝛽 EWs than galaxies even at 𝑧 ∼ 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 ≳ -18; 𝐿 UV ≲ 0 . 1 𝐿 ∗ UV ) 𝑧 ≳ 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, star formation history, and ionizing photon production and escape.  \nWith our JADES/NIRCam sample of 756 𝑧 ∼ 6 -9 Lyman-break galaxies, we now investigate how the reionization-era [OIII]+H 𝛽 and H 𝛼 EW distributions correlate with UV luminosity over a large dynamic range ( -22 ≲ M UV ≲ -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 §4.1, following with the results on the [OIII]+H 𝛽 (§4.2) and H 𝛼 (§4.3) EWdistributions.', '4.1 Methodology': 'Here, we describe how we quantify the UV luminosity dependence on the [OIII]+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 ≤ -19 . 5) subset, a faint ( -19 . 5 < M UV ≤ -18) subset, and a very faint (M UV > -18) subset. We also aim to test whether there is any significant redshift evolution in the [OIII]+H 𝛽 EW distribution during the reionization era and thus infer this distribution for the F775Wdropouts ( 𝑧 ∼ 6) and F090W dropouts ( 𝑧 ∼ 7 -9), separately. Because the H 𝛼 line redshifts out of the reddest JADES/NIRCam band (F444W) at 𝑧 > 6 . 6, we restrict our analysis of the H 𝛼 EW distribution to the F775W dropout sample.  \nTo infer the EW distribution among a given sub-population of our sample, we follow the Bayesian formalism described in Schenker et al. (2014) (see also Endsley et al. 2021a; Boyett et al. 2022):  \n𝑃 ( 𝜃 ) ∝ GLYPH<214> 𝑖 ∫ 𝑃 𝑖 ( EW ) 𝑃 ( EW | 𝜃 ) 𝑑 EW . (1)  \nIn this equation, 𝜃 encapsulates the parameters describing the assumed functional form of the distribution (e.g. the median and variance for a log-normal distribution) while 𝑃 𝑖 ( EW ) is the probability distribution function (PDF) on EW for the 𝑖 th galaxy in the subpopulation of interest. This formalism has two highly-desired features. First, the expression in the integral of Eq. 1 is small when the assumed EW distribution ( 𝑃 ( EW | 𝜃 ) ) is inconsistent with the PDF of galaxy 𝑖 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 in-  \nTable 2. Table summarizing the inferred [OIII]+H 𝛽 and H 𝛼 EWdistributions among our sample of 756 Lyman-break 𝑧 ∼ 6 -9 galaxies. The EW distributions are assumed to be log-normal and are parametrized by the median EW, 𝜇 EW , and standard deviation, 𝜎 EW . We divide our sample into six sub-populations divided by F775W dropouts ( 𝑧 ∼ 6) and F090W dropouts ( 𝑧 ∼ 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 sub-populations, we infer the EW distribution from the posteriors of the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc constant star formation history (CSFH) and two-component star formation history (TcSFH) models. We limit the H 𝛼 EW distribution results to the 𝑧 ∼ 6 sample given that H 𝛼 redshifts out of NIRCam at 𝑧 > 6 . 5. We recommend using the EW distributions derived from the TcSFH model outputs given our findings in §5.2 (see Fig. 14).  \nalaxies with very broad (i.e. largely unconstrained) PDFs.  \nThe 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. 𝑃 𝑖 ( EW ) ≡ posterior 𝑖 ( EW ) / prior ( 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 probability distribution functions for galaxies with lower signal-tonoise data, and hence would systematically impact the fainter populations more if we did not divide out the priors when computing 𝑃 𝑖 ( 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 [OIII]+H 𝛽 and H 𝛼 EW.  \nIn practice, we treat the integral in Eq. 1 as a Riemann sum to avoid having to approximate 𝑃 𝑖 ( EW ) as some functional form via, e.g., spline interpolation. Therefore, we rewrite Eq. 1 as  \n𝑃 ( 𝜃 ) ∝ GLYPH<214> 𝑖      ∑︁ 𝑗 𝑃 𝑖, 𝑗 ( EW ) 𝑃 𝑗 ( EW | 𝜃 )      (2)  \nwhere the index 𝑗 represents a bin in EW. 𝑃 𝑖, 𝑗 ( EW ) is the probability that object 𝑖 has an EW in bin 𝑗 , while 𝑃 𝑗 ( EW | 𝜃 ) is the integrated probability of the assumed functional EW distribution within bin 𝑗 . Therefore, both 𝑃 𝑖, 𝑗 ( EW ) and 𝑃 𝑗 ( EW | 𝜃 ) implicitly account for the width of bin 𝑗 .  \nFor the [OIII]+H 𝛽 EWinference, 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  \nbound of log(EW/Å) = 3.7 was chosen given that the maximum inferred [OIII]+H 𝛽 EWin our sample is ≈ 5000 Å (see §3.2.2). Similarly, the H 𝛼 EW bins applied in Eq. 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 𝛼 EWof ≈ 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.  \nFor each galaxy sub-population, we infer the EW distribution assuming a log-normal functional form defined by two parameters: the median of the EW distribution, 𝜇 EW , and its standard deviation, 𝜎 EW . This assumption of a log-normal 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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc constant star formation history (CSFH) and two-component star formation history (TcSFH) models. Due to our findings in Appendix B and §5.2 (see also Figs. 14 and B1), we adopt the results from the TcSFH models as fiducial.', '4.2 Results on [OIII]+H 𝛽 EWDistributions': 'We begin by quantifying the inferred [OIII]+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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc SED fits, we infer a median [OIII]+H 𝛽 EW of 𝜇 EW = 740 + 90 -80 Å among the bright ( ⟨ M UV ⟩ = -20 . 1) subset at 𝑧 ∼ 7 -9. This is consistent with median [OIII]+H 𝛽 EWs reported in the literature for similarly UV-luminous HST and NIRCam selected galaxy samples at 𝑧 ∼ 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 ≲ M UV ≲ -21 . 5) 𝑧 ∼ 7 galaxies selected from ground-based imaging (Endsley et al. 2021a, 2023). The inferred  \nFigure 8. The inferred distributions of [OIII]+H 𝛽 EW (top) and H 𝛼 EW (bottom) among our sample using the two-component SFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc SED fit outputs. The [OIII]+H 𝛽 EWdistributions are inferred for both the 𝑧 ∼ 6 and 𝑧 ∼ 7 -9 subsets while the H 𝛼 EW distribution is inferred only at 𝑧 ∼ 6 given this line redshifts out of NIRCam at 𝑧 > 6 . 5. The different colors show the EW distribution inferred among three different UV luminosity bins (see Table 2). We find a strong, highly-significant decline in the typical [OIII]+H 𝛽 EWwith UV luminosity, yet a very weak UV luminosity dependence on H 𝛼 EW(at least at 𝑧 ∼ 6).  \n<!-- image -->  \nstandard deviation of the EW distribution among our bright 𝑧 ∼ 7 -9 JADES/NIRCamsubset(0.32 + 0 . 04 -0 . 04 dex) is also consistent with values reported in the literature at similar luminosities (Endsley et al. 2021a, 2023).  \nDue to the strong degeneracy between Balmer breaks and nebular line emission on IRAC colors at 𝑧 ∼ 6, the [OIII]+H 𝛽 EWdistribution has never been quantified in the literature at this epoch. With the rich (five-band) 3-5 𝜇 mSEDsamplingnowaffordedbyJADES/NIRCam, wehere infer that the median [OIII]+H 𝛽 EWamongbright( ⟨ M UV ⟩ = -20 . 0) 𝑧 ∼ 6 galaxies is 𝜇 EW = 890 + 90 -80 Å, ≈ 0.08 dex ( ≈ 1.2 × ) higher than that inferred from our bright 𝑧 ∼ 7 -9 sub-sample. The standard deviation of [OIII]+H 𝛽 EWsamongthebright 𝑧 ∼ 6subset is inferred to be very similar to that of the bright 𝑧 ∼ 7 -9 subset (see Table 2).  \nNow equipped with a sample of over 600 Lyman-break 𝑧 ∼ 6 -9 galaxies with ≈ 2-30 × fainter UV luminosities ( -19 . 5 ≤ M UV ≲ -16 . 5), we investigate whether the [OIII]+H 𝛽 EW distribution correlates significantly with UV luminosity in the reionization era. Using the fiducial TcSFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc SED fits, we infer a median [OIII]+H 𝛽 EW of 𝜇 EW = 590 + 40 -50 Å and 380 + 80 -60 Å among the faint ( ⟨ M UV ⟩ = -18 . 7) and very faint ( ⟨ M UV ⟩ = -17 . 5) 𝑧 ∼ 6 subsets, respectively. This indicates a smooth and strong (factor of ≈ 2.3) decline in the typical [OIII]+H 𝛽 EW of 𝑧 ∼ 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 𝑧 ∼ 7 -9 we infer a median [OIII]+H 𝛽 EW of 𝜇 EW = 500 + 50 -40 Å and 300 + 40 -40 Å among the faint ( ⟨ M UV ⟩ = -18 . 6) and very faint ( ⟨ M UV ⟩ = -17 . 6) subsets, respectively, indicating a ≈ 2.5 × decline in typical [OIII]+H 𝛽 EW over ≈ 1 dex in 𝐿 UV (Fig. 8). Therefore, in both the 𝑧 ∼ 6 and 𝑧 ∼ 7 -9 samples, we infer a strong UV luminosity dependence on the median [OIII]+H 𝛽 EWat 𝑧 ∼ 6 -9 with high confidence ( ≳ 5 𝜎 ).  \nTo help visualize these trends, we plot the fiducial [OIII]+H 𝛽 EWs (the median of the posterior probability distribution) for each individual galaxy as a function of M UV in Fig. 9a,b. In the UV-bright (M UV < -19 . 5) bins, only a small fraction ( ≈ 10-15%) of 𝑧 ∼ 6 -9 LBGs have inferred [OIII]+H 𝛽 EWs of < 300 Å while this fraction rises to ≈ 50-65% in the very UV-faint ( -18 ≤ M UV ≲ -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 [OIII]+H 𝛽 EW is ≲ 300 Å at very low UV luminosities. With this very simple approach, we obtain much higher median [OIII]+H 𝛽 EWsof 880 Å and 700 Å in the UV-bright bins at 𝑧 ∼ 6 and 𝑧 ∼ 7 -9, respectively. Nonetheless, we emphasize that the EWdistribution inference approach (§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.  \nTo understand why the NIRCam SEDs are driving the models to yield a strong decline in [OIII]+H 𝛽 EW with decreasing UV luminosity, we analyze the long-wavelength NIRCam colors sensitive to these lines. We begin by considering the F775W dropout ( 𝑧 ∼ 6) sample where the large majority of galaxies lie at 5 . 5 < 𝑧 < 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 ([OIII], H 𝛽 , and H 𝛼 ) at 𝑧 ≈ 5 . 6 -6 . 7. Meanwhile. the F356W and F335M bands are contaminated by [OIII]+H 𝛽 at redshifts associated with the majority of the 𝑧 ∼ 6 subset ( 𝑧 ≈ 5 . 4 -6 . 9 and 𝑧 ≈ 5 . 5 -6 . 1, respectively). Therefore, we investigate whether the measured photometry of the 𝑧 ∼ 6 sample supports significant flattening in the F410M -F335M and F410M -F356W colors as we move from the bright to very faint subsets. In the bright 𝑧 ∼ 6 subset, the median F410M -F335M and F410M -F356W colors are 0 . 98 + 0 . 08 -0 . 10 and 0 . 63 + 0 . 02 -0 . 07 mag, respec-  \nFigure 9. The inferred fiducial EWs of each individual galaxy in our 𝑧 ∼ 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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc TcSFH fits. Those with EW < 300 Å are plotted as upper limits at that value (see §4.1 with red markers, 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).  \n<!-- image -->  \nstrong [OIII]+H 𝛽 emission. These colors 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 [OIII]+H 𝛽 EWs at lower UV luminosities. The uncertainties on median colors are derived using bootstrap resampling.  \nWe conduct a similar color 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 ( 𝑧 ≈ 6 . 3 -8 . 2) or F356W ( 𝑧 ≈ 7 . 2 -8 . 8) while the [OIII]+H 𝛽 emission is captured by either F410M ( 𝑧 ≈ 7 . 0 -7 . 6) or F444W ( 𝑧 ≈ 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 colors as we move from the bright to very faint subsets of our 𝑧 ∼ 7 -9 sample. In the bright subset, the median values for all four above colors 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 [OIII]+H 𝛽 emission. On the other hand, in the very faint subset we measure much flatter median colors (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 [OIII]+H 𝛽 EWs.  \nBecause the sample size of each M UV bin is fairly large ( 𝑁 gal = 55-200; see Table 2), we do not expect these median colors to be strongly impacted by noise. Nonetheless, we employ a Monte Carlo simulation to quantify the likelihood of measuring the median colors of the very faint subsets (accounting for noise) under the assumption that they indeed have true colors 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 re-normalize 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  \nvery faint galaxy. This Monte Carlo simulation is performed 5000 times and, in each iteration, we compute the median colors of interest.  \nThe distributions of median colors from the Monte Carlo simulations are shown in Fig. 10, along with the actual measured median colors of the very faint subsets plotted as vertical dashed lines in each panel. For nearly every color probing [OIII]+H 𝛽 emission, we find that very few ( ≤ 0.1%) of the 5,000 Monte Carlo iterations yield median recovered colors as flat (or flatter) than that measured among the very faint subset. While the median F410M -F356WMonteCarlo colors in the F775W dropout subset are consistent (at the ≈ 1 𝜎 level) with that measured among the very faint subset, F356W is less sensitive to strong [OIII]+H 𝛽 emission relative to F335M given its broader bandwidth and thus the significance in color difference is expected to be considerably weaker. Moreover, the smaller size of the very faint F775W dropout sample makes its median color uncertainties ≈ 2 × 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 colors relative to that measured for the very faint F775W dropout sample.  \nHaving now validated that the inferred trend in median [OIII]+H 𝛽 EW with UV luminosity is consistent with expectations given the measured long-wavelength colors, 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 colors implying extremely high [OIII]+H 𝛽 EWs of > 1500 Å (see, e.g., ID JADESGS+53.17835 -27.77879 and JADES-GS+53.19590 -27.79240 in Fig. 7). Accordingly, our assumed log-normal EW distributions are inferred to broaden considerably towards fainter UV luminosities at both 𝑧 ∼ 6 and 𝑧 ∼ 7, with 𝜎 EW 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 [OIII]+H 𝛽 EWs ( > 1500 Å) declines significantly towards fainter UV luminosities, going from 24 + 5 -4 % (18 + 5 -4 %) in the bright 𝑧 ∼ 6 ( 𝑧 ∼ 7 -9) subset to 13 + 4 -3 %(8 + 2 -2 %) in  \nFigure 10. Illustration of results from Monte Carlo simulation test validating that the substantial decline in median [OIII]+H 𝛽 EW at fainter UV luminosities is very unlikely to be caused by photometric noise. The shaded histograms show the distribution of recovered median colors (one color 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 colors sensitive to [OIII]+H 𝛽 contamination, very few ( ≤ 0.1%) of the 5,000 Monte Carlo iterations yield median recovered colors as flat (or flatter) as that measured among the very faint subset (vertical dashed lines). The median F410M -F444W color of the very faint 𝑧 ∼ 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.  \n<!-- image -->  \nthe respective very faint subset. Conversely, the fraction of galaxies with relatively very low [OIII]+H 𝛽 EWs ( < 300 Å) rises dramatically towards the faint-end of the luminosity function, increasing from 6 + 3 -2 %(11 + 5 -4 %) in the bright 𝑧 ∼ 6 ( 𝑧 ∼ 7 -9) subset to 42 + 7 -7 % (50 + 5 -5 %) in the respective very faint subset.', '4.3 Results on H 𝛼 EWDistributions': "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 𝑧 ≈ 4 . 0 -5 . 0 where several works have inferred typical H 𝛼 EWsof ∼ 400 Å among HST -selected galaxies (Stark et al. 2013; Smit et al. 2016; Rasappu 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 𝑧 ∼ 8 galaxies and recovered a considerably larger typical H 𝛼 EW( ≳ 1000 Å), perhaps implying a significant redshift evolution in the median H 𝛼 EWat high redshifts. With the JADES NIRCam data, we are now able to infer the H 𝛼 EW distribution at 𝑧 ∼ 6 and moreover investigate whether it correlates strongly with UV luminosity as may be expected given our results on the [OIII]+H 𝛽 EWdistribution above.  \nFrom the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc TcSFH SED fits, we infer very similar median H 𝛼 EWs among all three UV luminosity bins at 𝑧 ∼ 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 𝑧 ∼ 4 . 5 and 𝑧 ∼ 6. But perhaps the most striking result here is that we find no evidence for a strong change in the typical 𝑧 ∼ 6 H 𝛼 EW across 1 dex in UV luminosity, in marked contrast with our findings on the [OIII]+H 𝛽 EWdistribution in §4.2.  \nEven 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.  \nWe again verify that the measured long-wavelength NIRCam colors are consistent with little-to-no change in the median H 𝛼 EW between the bright and very faint 𝑧 ∼ 6 subsets following the same approach as described in §4.2. At 𝑧 ∼ 6, the F410M band provides a relatively clean probe of the rest-optical continuum (see §4.2) while the F444W band contains H 𝛼 . The median F410M -F444W color of the bright 𝑧 ∼ 6 subset (0 . 35 + 0 . 05 -0 . 04 mag) is nearly identical to that of the very faint subset (0 . 33 + 0 . 07 -0 . 14 mag), consistent with a weak trend in typical H 𝛼 EWwith UV luminosity. The Monte Carlo simulation test described in §4.2 also reveals that even accounting for signal-tonoise differences in the bright versus very faint 𝑧 ∼ 6 subsets often yields nearly identical median F410M -F444W colors between the two populations (see Fig. 10). We thus conclude that the photometry is consistent with little change in typical 𝑧 ∼ 6 H 𝛼 EWs over the UV luminosity range spanned by our sample ( -22 ≤ M UV ≤ -16 . 4).  \nThere is moderate evidence ( ≈ 2 𝜎 ) for a slight ( ≈ 0.05 dex) increase in the width of the 𝑧 ∼ 6 H 𝛼 EWdistribution towards very faint UV luminosities (see Table 2). This is consistent with the slightly larger fraction of very faint 𝑧 ∼ 6 galaxies with fiducial H 𝛼 EWs < 300 Å (see Fig. 9c) even though the median EW remains effectively unchanged. Nonetheless, 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 log-normal form across the entire UV luminosity range of our 𝑧 ∼ 6 sample ( -22 ≤ M UV ≤ -16 . 4), we obtain 𝜇 EW =600 + 40 -10 Åand 𝜎 EW = 0.27 ± 0.01 dex. Given the strong UV luminosity dependence on [OIII]+H 𝛽 EWinferred in §4.2, it is no surprise that we then infer  \n<!-- image -->  \n/circledot  \n<!-- image -->  \nFigure 11. Illustration of how the EWs of [OIII]+H 𝛽 (blue solid lines) and H 𝛼 emission correlate with metallicity (panel a), recent star formation history, 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 𝑍 ⊙ , log 𝑈 = -2 . 5, 𝑓 esc = 0, and 𝜏 𝑉 = 0 with the Gutkin et al. (2016) models implemented in /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc.  \n<!-- image -->  \na strong UV luminosity dependence on the [OIII] 𝜆 5007/H 𝛼 ratio (hereafter [OIII]/H 𝛼 ) as well. In the bright 𝑧 ∼ 6 subset, we infer a median [OIII]/H 𝛼 ≈ 1.6 while in the faint and very faint subsets we infer median [OIII]/H 𝛼 ratios of ≈ 1.3 and ≈ 0.9, respectively.", '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 [OIII]+H 𝛽 EWwith UV luminosity at 𝑧 ∼ 6 -9, yet an approximately constant H 𝛼 EWdistribution across an order of magnitude in L UV at 𝑧 ∼ 6 (Fig. 8). Here, we begin by discussing what physically may be driving the different L UV trends for [OIII]+H 𝛽 and H 𝛼 EWsthen explore the results implied from the SEDfits, 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 [OIII]+H 𝛽 emission declines substantially at ≲ 0.2 𝑍 ⊙ (see Fig. 11a). Because fainter, lower-mass galaxies are generally expected to be less chemically enriched, this likely contributes to much weaker [OIII]+H 𝛽 EWs in our faintest subset. We consider how much lower the typical metallicity of very faint ( ⟨ M UV ⟩ ≈ -17 . 5) 𝑧 ∼ 6 -9 galaxies would have to be relative to the bright subset ( ⟨ M UV ⟩ ≈ -20) to completely explain the ≈ 2.4 × decline in median [OIII]+H 𝛽 EW inferred in §4.2. For this test, we utilize the Gutkin et al. (2016) models employed in /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc 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 [OIII]+H 𝛽 EW(see Fig. 11a).  \nSuch 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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc fits, we expect that lowering the metallicity by 1 dex from 0.2 𝑍 ⊙ to 0.02 𝑍 ⊙ would result  \nin a ≈ 50% increase in H 𝛼 EWs. Therefore, while lower metallicities at fainter UV luminosities may partially explain the decreasing [OIII]+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 𝑧 ∼ 6.', '5.2 Evidence for Bursty Star Formation Histories at z ≳ 6': "One potential way to offset the expected increase in H 𝛼 EWatfainter 𝐿 UV due to lower metallicity is by invoking differences in the recent star formation histories of bright versus very faint reionization-era galaxies (we discuss an alternative solution with differing ionizing photon escape fractions below 6 ). Because the Balmer and [OIII] 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 timescale results in weaker emission lines, decreasing the EWs of H 𝛼 and [OIII]+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 [OIII]+H 𝛽 EWs in the very UV-faint population (Fig. 8), working in tandem with the shift toward lower metallicities.  \nTo better assess whether the photometric data support systematic changes in metallicity and recent star formation histories across the wide UV luminosity range spanned by our sample ( -22 ≲ M UV ≲ -16 . 5), we utilize the results of our two-component star formation history (TcSFH) /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc fits. For each galaxy, we compute the inferred ratio of the star formation rate 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 timescales most relevant for nebular line emission. The 50 Myr timescale adopted in the denominator of this ratio was chosen to reflect the typical lightweighted age of the sample, though we note that our findings below are qualitatively unchanged if we instead adopt a 30 Myr or 100 Myr timescale. We quantify the distributions of SFR 3 Myr / SFR 50 Myr  \nand metallicity in the same three UV luminosity bins as considered for the EW distribution analysis (see §4.1), again adopting a lognormal distribution for each M UV bin and applying Eq. 2. Here, we focus our analysis solely on the 𝑧 ∼ 6 sample given the lack of H 𝛼 information for the large majority of the 𝑧 ∼ 7 -9 sample. Because the SFR 3 Myr / SFR 50 Myr ratio is by definition restricted to values ≤ 50 3 (this maximum value occurs when no stars are assumed to form between 3-50 Myr ago), we truncate the log-normal distributions at this upper limit (i.e. 𝑃 (SFR 3 Myr / SFR 50 Myr > 50 3 | 𝜃 ) = 0). Similarly, the log-normal distributions for metallicity are truncated outside the fitted parameter space in our SED fits ( -2 . 2 ≤ log( 𝑍 / 𝑍 ⊙ ) ≤ -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 log-normal functional form.  \nThe TcSFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc outputs suggest a significant anti-correlation between UV luminosity and metallicity at 𝑧 ∼ 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 𝑍 ⊙ in the bright ( ⟨ M UV ⟩ = -20 . 0), faint ( ⟨ M UV ⟩ = -18 . 7), and very faint ( ⟨ M UV ⟩ = -17 . 5) subsets, respectively. Therefore, in context of these models, changes in metallicity contribute significantly to the different UV luminosity trends with [OIII]+H 𝛽 and H 𝛼 EWs.Wecheckwhethertheseinferredmetallicities 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 𝑧 ∼ 6 sample, the typical stellar mass at M UV ≈ -20 is 𝑀 ∗ ∼ 4 × 10 8 𝑀 ⊙ while it is ∼ 2 × 10 7 𝑀 ⊙ at M UV ≈ -17 . 5 (see Fig. 2). The inferred metallicities above from the TcSFH models are reasonably consistent (factor of ≲ 2 difference) with the predicted 𝑧 ∼ 6 gas-phase metallicities at the respective stellar masses from the /i.pc/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pc /t.pc/n.pc/g.pc (Torrey et al. 2019), F/i.pc/r.pc/s.pc/t.pcL/i.pc/g.pc/h.pc/t.pc (Langan et al. 2020), and /a.pc/s.pc/t.pc/r.pc/e.pc/a.pc/u.pc/s.pc (Ucci et al. 2023) simulations, though systematically ≈ 5 × higher than that predicted from the /f.pc/i.pc/r.pc/e.pc simulations (Ma et al. 2016).  \nWe now investigate whether the TcSFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc fit outputs also imply differences in the recent SFHs of bright vs. very faint 𝑧 ∼ 6 galaxies. Using Eq. 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 𝑧 ∼ 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 veryrecently 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 log-normal distribution is inferred to be 0 . 52 + 0 . 09 -0 . 07 dex, 0 . 62 + 0 . 08 -0 . 06 dex, 1 . 21 + 0 . 66 -0 . 33 dex in the bright, faint, and very faint subsets, respectively, implying generally more variation in recent SFHs at fainter luminosities (and hence lower stellar masses).  \nThere are two particular subsets of our sample that provide evidence in favor of bursty SFHs among reionization-era galaxies. The first subset consists of galaxies showing extremely strong nebular line emission. The very high [OIII]+H 𝛽 EWs( ∼ 2000-5000 Å) confidently inferred among ≈ 50 galaxies in our sample (see §3.2.2) imply that a considerable fraction of Lyman-break 𝑧 ∼ 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 star formation histories of two of the confident 𝑧 ∼ 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  \nwith 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 𝑧 ≳ 6 galaxies using different methods (Dome et al. 2023; Dressler et al. 2023b,a; Looser et al. 2023a,b; Strait et al. 2023). 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 et al. 2018; Faucher-Giguère 2018; Ma et al. 2018; Furlanetto & Mirocha 2022).  \nAnother subset of galaxies in our JADES sample may be caught during a relatively inactive phase of star formation. We identify 𝑧 ∼ 6 galaxies exhibiting NIRCam SEDs consistent with relatively weak Balmer breaks (thus implying young light-weighted ages) yet also weak [OIII]+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 𝑧 ∼ 6 JADES sample below). In the left panels of Fig. 14, we show the fitted SEDs from the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc 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 𝑍 ⊙ ) to simultaneously explain the weak Balmer breaks 7 ( 𝐹 𝜈 ( 4200 ) / 𝐹 𝜈 ( 3500 ) ≈ 1 as implied by the fairly flat F200W -F277W colors) as well as the low [OIII]+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 over-predict the H 𝛼 EWsrelative to that implied by the observed F410M -F444W colors, suggesting that the SEDs of these galaxies cannot be well explained by our fiducial CSFH models.  \nWhenallowingfor more flexible star formation histories, we obtain SEDmodelsthat yield much better fits to the photometry of these two illustrative 𝑧 ∼ 6 galaxies. As shown in the middle panels of Fig. 14, the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc TcSFH models can reproduce the NIRCam photometry in all nine bands with moderately low-metallicity solutions ( ≈ 0.050.1 𝑍 ⊙ ) for these two 𝑧 ∼ 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 star formation histories yield relatively few O stars that are still alive to power the [OIII] 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 ∼ 1030 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 reionizationera galaxies with young light-weighted ages and weak [OIII]+H 𝛽 emission was first tentatively identified in Endsley et al. (2023) from early JWST data, though the lack of H 𝛼 constraints for the 𝑧 ∼ 6 . 5 -8 sample considered therein left much uncertainty on the physical origin of the weak line emission (see also Simmonds et al. 2023; Prieto-Lyon et al. 2023). With the improved imaging sensitivity and SED sampling of JADES and GOODS, we have now extended this investigation to a 𝑧 ∼ 6 sample and find that our fiducial CSFH models cannot adequately explain the photometry of at least some of these very early galaxies.  \nFigure 12. Illustration of how the strong decline in [OIII]+H 𝛽 EW with UV luminosity yet approximately constant H 𝛼 EW distributions at 𝑧 ∼ 6 (see Fig. 8) can be explained by either the combination of lower metallicity and more recently declining star formation histories 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.  \n<!-- image -->  \nConstant SFH with 0  \n<  \nf  \nesc  \n<  \n0  \n.  \n8 ;  \nz  \n6  \n<!-- image -->  \n/circledot  \n/circledot  \nIn addition to the two illustrative objects shown in Fig. 14, we identify a few other 𝑧 ∼ 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 𝑧 ∼ 6 galaxies with > 84% probability of having [OIII]+H 𝛽 and H 𝛼 EWs of < 800 Å from the TcSFH posterior outputs. We then apply the color-color 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 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) timescales 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 𝜒 2 CSFH > 14 (we are fitting 14 data photometric data points) and 𝜒 2 CSFH -𝜒 2 TcSFH > 10. This rather strict selection results in a total of seven 𝑧 ∼ 6 galaxies, including the two illustrative examples shown in Fig. 14. The majority of these  \nseven galaxies have inferred SFR 3 Myr / SFR 50 Myr ratios of < 0 . 2 with > 84% confidence from the TcSFH posterior outputs.  \nWe emphasize that the seven objects described above are not an exhaustive census of 𝑧 ∼ 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) ≳ 30 before adding the 5% 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 𝑧 ∼ 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 𝑧 = 7 . 3 '(mini-)quenched' galaxy from recent JWST observations (Looser et al. 2023a; Dome et al. 2023; for similar examples at 𝑧 ∼ 5 see also Looser et al. 2023b; Strait et al. 2023).  \n∼  \nSFHs of Extreme Line Emitters  \nFigure 13. Example inferred star formation histories of two 𝑧 ∼ 6 extreme emission line galaxies (top; SEDs shown in Fig. 7) and two 𝑧 ∼ 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% credible intervals (gray regions) are taken from the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc 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.  \n<!-- image -->", '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 star formation histories 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 lowmass objects that have been up-scattered 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 𝑧 ∼ 6 -9 galaxies (see §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 𝑧 > 10 galaxies recently identified by JWST as one would then expect a higher fraction of low-mass galaxies to be up-scattered above the detection threshold (e.g. Mason et al. 2023; Mirocha & Furlanetto 2023; Sun et al. 2023).  \nIn 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. Nonetheless, the substantial scatter in the SFR 3 Myr / SFR 50 Myr ratio inferred among  \nour very UV-faint 𝑧 ∼ 6 subset ( ≈ 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 [OIII]+H 𝛽 EWs ( ≳ 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 𝑧 spec = 7 . 28 described in Saxena et al. 2023) indicating a strong rise in their SFHs on ∼ 3 Myr timescales. 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 𝐿 UV / 𝑀 ∗ ) as they are to be a relatively highmass object experiencing a recent strong decline in SFR (yielding relatively low 𝐿 UV / 𝑀 ∗ ), 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 𝑧 ∼ 6 [OIII]+H 𝛽 and H 𝛼 EW trends with UV luminosity found in §4.2-4.3. In this sub-section, 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 [OIII] 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 Lyman-continuum (LyC; 𝜆 rest < 912 Å) photon escape fractions required to explain the photometric data (and resulting [OIII]+H 𝛽 and H 𝛼 EW trends) under the assumption of a constant SFH by using the 'picket-fence' model within /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc. In these picket fence models (e.g. Conselice et al. 2000; Heckman et al. 2011; Rivera-Thorsen et al. 2017), a given fraction (hereafter 𝑓 esc ) of LyC photons are assumed to escape the HII 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𝑓 esc relative to the base CSFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc models (see Fig. 11c).  \nIn the picket-fence /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc SED fits, we adopt a log-uniform prior on 𝑓 esc between 0.001 and 0.8 with all other priors equivalent to the nominal CSFH fits. We then infer log-normal distributions for both metallicity and 𝑓 esc using Eq. 2 in the three different UV luminosity bins for our 𝑧 ∼ 6 sample. While density-bounded photoionization models may be more appropriate at the highest 𝑓 esc values ( ≳ 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 constant star formation histories. Adopting density-bounded models (Plat et al. 2019) would generally result in a slower decline of the [OIII] EWs with increasing 𝑓 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. Nonetheless, the difference in EWs between the ionization-bounded and density-bounded /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc models is ≲ 10% at 𝑓 esc < 0 . 8.  \nThe picket-fence CSFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc fits also yield an anti-correlation between UV luminosity and metallicity at 𝑧 ∼ 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 𝑍 ⊙ among the bright ( ⟨ M UV ⟩ = -20 . 0), faint ( ⟨ M UV ⟩ = -18 . 7), and very faint ( ⟨ M UV ⟩ = -17 . 5) 𝑧 ∼ 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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc fits adopting 𝑓 esc = 0 (c.f. panel  \nFigure 14. Illustration of two 𝑧 phot ≈ 6 . 0 galaxies (one per row) showing SEDs consistent with young light-weighted ages as well as weak [OIII]+H 𝛽 and H 𝛼 emission. We find that our fiducial CSFH models (with 𝑓 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.  \n<!-- image -->  \n(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 ( ≳ 0.3 𝑍 ⊙ ) 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 ≤ log( 𝑍 / 𝑍 ⊙ ) ≤ -0 . 3 parameter space allowed in our SED fits. It is unclear whether such a distribution in 𝑧 ∼ 6 galaxy metallicity is physical at such faint UV luminosities and low stellar masses (maximum 𝑀 ∗ = 10 9 𝑀 ⊙ ; see Fig. 2). Nonetheless, we proceed with investigating what correlation between UV luminosity and 𝑓 esc is required in context of these models to explain the photometry of the 𝑧 ∼ 6 sample.  \nUnder the assumption of constant star formation histories, the trends in [OIII]+H 𝛽 and H 𝛼 EW found in §4.2-4.3 imply a strong increase in 𝑓 esc towards very faint UV luminosities at 𝑧 ∼ 6 (Fig. 12d). With the picket-fence /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc SED models, we infer median 𝑓 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 log-normal distributions for 𝑓 esc in each luminosity bin. Moreover, the fraction of galaxies with substantial ionizing photon leakage ( 𝑓 esc > 20%) is inferred to be 28 + 8 -6 %, 52 + 6 -6 %, 76 + 9 -10 % in each of the respective M UV bins. As can be seen in Fig. 12, the inferred distribution of 𝑓 esc simultaneously narrows and peaks closer to the allowed upper bound of 𝑓 esc = 0.8 towards fainter UV luminosity. The generally bluer UV slopes, smaller rest-UV sizes, and lower inferred dust attenuation among fainter 𝑧 ≳ 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 ∼ 5) change in median 𝑓 esc between M UV ≈ -20 and M UV ≈ -17 . 5 at 𝑧 ∼ 6 would be expected. No evidence yet exists among local  \nsamples that the LyC escape fraction strongly correlates with far-UV luminosity nor stellar mass (Izotov et al. 2021; Flury et al. 2022).  \nWe consider whether the picket-fence CSFH models can explain the photometry of the seven 𝑧 ∼ 6 galaxies confidently showing SEDs consistent with young light-weighted ages as well as relatively weak nebular lines discussed in §5.2. Five of these seven galaxies are considerably better fit with the picket-fence models relative to standard CSFH models ( Δ 𝜒 2 > 10) and, in each of these five cases, the inferred LyC escape fraction is 𝑓 esc > 50% with > 84% 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 𝜒 2 = 19 . 5 across the 14 fitted data points, implying a somewhat poorer fit relative to the best-fitting TcSFH model ( 𝜒 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 ( Δ 𝜒 2 < 5). As discussed in the following sub-section, 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 𝑧 ∼ 6 [OIII]+H 𝛽 and H 𝛼 EW trends found in §4.2-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 muchgreater LyC leakage among fainter systems with constant SFHs at all luminosities (Fig. 12). To better understand the potential impli-  \nthe ionizing efficiency of bright vs. very faint galaxies in these two limiting cases of model assumptions (i.e. bursty SFH with 𝑓 esc ≈ 0 or CSFH with 0 < 𝑓 esc < 0 . 8), we infer the distribution of ionizing photon production efficiencies ( 𝜉 ion ) as a function of UV magnitude in each case. Here, we are considering the 𝜉 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 𝐿 1500 for dust attenuation or nebular continuum emission).  \nWith our TcSFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc models that adopt 𝑓 esc = 0, we infer median ionizing photon production efficiencies of log [ 𝜉 ion /( erg -1 Hz )] = 25 . 65 + 0 . 03 -0 . 03 , 25 . 48 + 0 . 03 -0 . 02 , and 25 . 35 + 0 . 05 -0 . 05 in the bright ( ⟨ M UV ⟩ = -20 . 0), faint ( ⟨ M UV ⟩ = -18 . 7), and very faint ( ⟨ M UV ⟩ = -17 . 5) 𝑧 ∼ 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 𝜉 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 𝑧 spec = 6 . 2 -6 . 8 galaxies with -17 . 2 ≤ M UV ≤ -16 . 5 (log [ 𝜉 ion /( erg -1 Hz )] = 25.4 ± 0.2). The standard deviation of the assumed log-normal distribution in 𝜉 ion is inferred to moderately increase from ≈ 0.20 dex in the bright and faint bins to ≈ 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 𝜉 ion decreases with UV luminosity, there remains a significant subset of a very UV-faint 𝑧 ∼ 6 galaxies with very efficient production of hydrogen ionizing photons ( 𝜉 ion ≳ 10 25 . 7 erg Hz -1 ; see Fig. 7). Similar sources at extremely low UV luminosities (M UV ∼ -16) have now been spectroscopically confirmed (Atek et al. 2024).  \nThe UV luminosity dependence on 𝜉 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 [ 𝜉 ion /( erg -1 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 ≈ 0.15-0.18 dex in all bins. Such a small UV luminosity dependence on 𝜉 ion yet large dependence on [OIII]+H 𝛽 is accommodated in the models by pushing 𝑓 esc to the very large values ( ≳ 0.5) in the very faint population. We also note that the CSFH models allow for less dynamic range in the inferred 𝜉 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 𝑧 ≳ 6), in contrast to scenarios with the TcSFH models where the SFR can sharply decline in the most recent ∼ 3 Myr (see Fig. 13). These differences in model assumptions contribute to the systematically higher and narrower inferred 𝜉 ion distributions with the CSFH models.  \nThe 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 𝑓 esc towards very faint UV luminosities yet a weak UV-luminosity dependence on 𝜉 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 UV have small 𝑓 esc , but that the most UV-luminous galaxies have systematically larger 𝜉 ion . While the assumption that all reionization-era galaxies have small 𝑓 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 UV directly result in  \nlarge uncertainties on the relative contribution of galaxies along the UV luminosity function to reionization.  \nIn summary, the UV-luminosity dependence on [OIII]+H 𝛽 and H 𝛼 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 𝑧 ∼ 6 -9 galaxies. We have demonstrated that these trends can be explained by either lower metallicity and more recently declining star formation histories in the very UV-faint population, or lower metallicities and higher LyC escape fractions in the very UVfaint 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 𝛼 -enhanced metallicities could also play a significant role in driving the EW distribution trends found in §4.2-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 [OIII]+H 𝛽 and H 𝛼 EWtrends.', '6 OVERDENSITIES AROUND LY 𝛼 EMITTERS AT Z > 7': "It has long been proposed that galaxies exhibiting strong Ly 𝛼 emission (EW > 25 Å) at 𝑧 > 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. 2023). In this scenario, strong 𝑧 > 7 Ly 𝛼 emitters (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 GOODSS imaging to better test the connection between overdensities and strong 𝑧 > 7 Ly 𝛼 emission.  \nWe search the literature for 𝑧 > 7 galaxies with known highconfidence Ly 𝛼 detections in the GOODS fields (i.e. S/N > 10 in Ly 𝛼 alone or S/N > 10 of a systemic line plus a > 7 𝜎 Ly 𝛼 detection), resulting in three objects: JADES-GS-z7-LA at 𝑧 Ly 𝛼 = 7 . 281 (Saxena et al. 2023), z8\\_GND\\_5296 at 𝑧 Ly 𝛼 = 7 . 508 (Finkelstein et al. 2013; Tilvi et al. 2016; Jung et al. 2019, 2020), and z7\\_GND\\_16863 at 𝑧 Ly 𝛼 = 7 . 599 (Jung et al. 2019, 2020). Two of these systems fall within the JADES/NIRCam footprint (z8\\_GND\\_5296 and JADESGS-z7-LA) and we investigate whether either show a significant photometric overdensity. The galaxy z8\\_GND\\_5296 is a UV-bright (M UV = -21 . 5) system in the GOODS-N field with a most recent Ly 𝛼 EW measurement of 33 ± 4 Å from Jung et al. (2020). JADESGS-z7-LA is, on the other hand, a very UV-faint (M UV = -16 . 7) galaxy in the GOODS-S field with a near-unity Ly 𝛼 escape fraction yielding an extremely large EW of 400 ± 90 Å implying that it lies within a very large ionized bubble (Saxena et al. 2023). Both of these strong LAEs exhibit very high [OIII]+H 𝛽 EWs of ≈ 1500 Å placing them in the top ≈ 17% and ≈ 7% of the [OIII]+H 𝛽 EW distribution, respectively, given their UV luminosities (see §4.2).  \nHere 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. 2023a) of the two 𝑧 > 7 LAEs in JADES noted above given that deep spectroscopic follow-up of 𝑧 ∼ 7 -8 candidates surrounding these rare sources remains incomplete. We thus aim to supplement our base 𝑧 ∼ 7 -9 selection criteria with ad-  \nto better identify neighboring galaxies lying relatively close in redshift space to the two strong 𝑧 > 7 LAEs in JADES. Fortunately, both of these LAEs lie at 𝑧 = 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 𝑧 ≈ 7 . 0 -7 . 6, the [OIII] 𝜆𝜆 4959,5007 doublet falls in the F410M filter yet outside of the F356W bandpass and thus galaxies in this redshift range with large [OIII]+H 𝛽 EWs will show remarkably red F356W -F410M colors. From the inferred [OIII]+H 𝛽 EW distributions in the JADES dataset (see Table 2 and Fig. 8), we expect the majority of M UV < -18 galaxies at 𝑧 ≈ 7 . 3 will exhibit [OIII]+H 𝛽 EWsof > 400 Å which, assuming a flat rest-optical continuum (in 𝐹 𝜈 ), will result in F356W -F410Mcolors of ≳ 0.6 mag at 𝑧 ≈ 7 . 0 -7 . 6. Such prominent colors should be relatively easy to identify for all but the faintest objects in our deep JADES imaging. Moreover, at 𝑧 ≈ 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 colors given the typically little dust attenuation inferred for the sample ( 𝐴 𝑉 ≲ 0 . 01 mag).  \nGuided by the above expectations for the colors of 𝑧 ≈ 7 . 0 -7 . 6 galaxies, we identify potential neighbors of the very UV-faint (M UV = -16 . 7) 𝑧 = 7 . 28 Ly 𝛼 emitter JADES-GS-z7-LA and the UVbright (M UV = -21 . 5) 𝑧 = 7 . 51 Ly 𝛼 emitter z8\\_GND\\_5296 by supplementing our F090W dropout selection criteria (§2.3) with the below additional cuts:  \n- (i) F356W - F410M > 0.6\n- (ii) F335M - F356W < 0.3\n- (iii) F150W < 29.0  \nHere, we have added the condition of F150W < 29 (corresponding to M UV ≲ -18 . 0) to ensure that our selection of candidate neighbors is fairly complete at 𝑧 ≈ 7 . 0 -7 . 6 given the different depth tiers of the JADES imaging (Eisenstein et al. 2023). The on-sky positions of the candidate 𝑧 ≈ 7 . 0 -7 . 6 galaxies in the JADES GOODSS and GOODS-N footprints are shown in Fig. 15, along with the positions of the two strong 𝑧 > 7 LAEs covered by JADES. In the GOODS-S field containing the very UV-faint LAE, we identify 43 potential 𝑧 ≈ 7 . 0 -7 . 6 galaxies with F150W < 29 across the 42.2 arcmin 2 JADES area with F335M coverage. Across the 47.2 arcmin 2 JADES/GOODS-N area with F335M coverage, we find 42 𝑧 ≈ 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 𝑧 ≈ 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 𝑧 ∼ 7 Lymanbreak galaxies is considerably higher over GOODS-N than GOODSS (e.g. Bouwens et al. 2015; Finkelstein et al. 2015). Below, we first discuss the surface density of potential 𝑧 ≈ 7 . 0 -7 . 6 neighbors surrounding the LAEs selected from JADES, then proceed to calculate the surrounding photometric overdensities folding in the expected completeness of our selection.  \nIt is immediately clear from Fig. 15 that the surface density of 𝑧 ≈ 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 𝑧 ≈ 7 . 0 -7 . 6 galaxies with a maximum separation of 23 arcsec from one another. The surface density of 𝑧 ≈ 7 . 0 -7 . 6 galaxies in this very small (0.073 arcmin 2 ) region of GOODS-S is 87.3 galaxies/arcmin -2 , 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 𝑧 ≈ 7 . 0 -7 . 6 galaxies have absolute UV magnitudes  \non  \ni  \nt  \nna  \ni  \nl  \nc  \ne  \nD  \nMUV  \nRight AscensionFigure 15. On-sky distributions of candidate 𝑧 ≈ 7 . 0 -7 . 6 galaxies in the two GOODS fields (colored circles) with the known high-confidence 𝑧 > 7 Ly 𝛼 emitter in each field shown as a blue star. JADES-GS-z7-LA is a very UV-faint (M UV = -16 . 7) galaxy at 𝑧 = 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 𝑧 = 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) neighboring 𝑧 ≈ 7 . 0 -7 . 6 galaxies approximately 1 arcmin to the East of JADES-GS-z7-LA, as well as two bright (M UV ≈ -20 . 3) such galaxies ≈ 1 arcmin to the South.  \n<!-- image -->  \n53.0  \nin the range -19 . 3 ≤ M UV ≤ -18 . 2 and [OIII]+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.  \nWe also identify five candidate 𝑧 ≈ 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). Considering the smallest rectangular area that encloses all six of these objects, we obtain a surface density of 3.6 galaxies/arcmin -2 which is approximately 4 × higher than the average over the full JADES/GOODS-N footprint with F335M. These five relatively nearby objects have inferred [OIII]+H 𝛽 EWs of ≈ 500-1500 Å and UV magnitudes of -20 . 0 ≤ M UV ≤ -18 . 4, indicating that they are all substantially ( ≥ 4 × ) fainter than z8\\_GND\\_5296.  \nWe now quantify the photometric overdensities around the two strong 𝑧 = 7 . 3 -7 . 5 LAEs in JADES by comparing the nearby 𝑧 ≈ 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. (2023a) which we outline below. As a first step towards comparing with the expected cosmic mean density, we must first correct our measured 𝑧 ≈ 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. 2023a) by generating mock SEDs each with a flat continuum (in 𝐹 𝜈 ) and IGM attenuation incorporated using the analytic model from Inoue et al. (2014). Given that our selection for this overdensity analysis depends on [OIII]+H 𝛽 emission via the F356W -F410M and F335M -F356W color cuts above, we also add [OIII]+H 𝛽 emission into the mock SEDs following our results on the 𝑧 ∼ 7 -9 EW distribution in §4.2, considering additional constraints on the 𝑧 ∼ 7 EW distribution among the very bright ( -22 . 5 ≲ M UV ≲ -21) population from Endsley et al. (2023). Specifically, for a given M UV , we assume a log-normal [OIII]+H 𝛽 EW distribution defined by median EW log( 𝜇 EW /Å) = max([2.87 , -0 . 16 ( M UV + 20 ) + 2 . 87]) and standard deviation 𝜎 EW = min([0.25 , 0 . 08 ( M UV + 20 ) + 0 . 3]) dex and adopt a fixed [OIII] 𝜆 5007/H 𝛽 ratio of 6 (e.g. Tang et al. 2019). In doing so, we account for our strong incompleteness to objects with relatively weak [OIII]+H 𝛽 emission (EW ≲ 400 Å) for this overdensity analysis, though we note that such objects are expected to be in the minority at 𝐹 150 𝑊 < 29 (M UV ≲ -18) given the derived EW distributions.  \nIn the selection completeness simulations, for each grid point in redshift ( 𝑧 = 6 . 5 -8 . 5, 0.05 unit spacing) and absolute UV magnitude ( -22 ≤ M UV ≤ -17, 0.25 mag spacing), we generate 1000 mock SEDs each with a different [OIII]+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 GOODSS 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 𝑧 ∼ 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.  \nThe resulting simulated completeness of our 𝑧 ≈ 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 ≳ 50% at 𝑧 ≈ 7 . 05 -7 . 60 at the brightest  \nFigure 16. Simulated completeness of the 𝑧 ≈ 7 . 0 -7 . 6 selection criteria where we supplement the base 𝑧 ∼ 7 -9 selection criteria (§2.3) with additional color cuts that utilize the long-wavelength medium band filters to identify [OIII]+H 𝛽 emitters within this redshift range (see §6). We show the simulated completeness as a function of M UV and redshift in the four different JADES/NIRCam imaging regions where the cutoffs at M UV ≈ -18 are due to the requirement of F150W < 29 for this selection.  \n<!-- image -->  \nluminosities and very low at redshifts 𝑧 < 6 . 9 and 𝑧 > 7 . 7 for all four regions. We also find that the completeness is ≳ 50% at absolute UV magnitudes of M UV ≲ -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 ≲ -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 the expected surface density counts of 𝑧 ≈ 7 . 0 -7 . 6 galaxies in each imaging region under the assumption of cosmic mean density. Following the methods of Whitler et al. (2023a), 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 𝑧 > 7 Ly 𝛼 emitter. In each ring, the amplitude of the photometric overdensity at that angular separation is computed as the surface density of objects identified by our 𝑧 ≈ 7 . 0 -7 . 6 selection divided by the cosmic mean.  \nGiven its near-unity Ly 𝛼 escape fraction and relatively small Ly 𝛼 velocity offset ( ≈ 120 km/s), JADES-GS-z7-LA is expected to reside in a very large ( 𝑅 ≳ 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 neighboring galaxies. With the JADES/GOODSS imaging, we have identified a very dense group of 𝑧 ≈ 7 . 0 -7 . 6 galaxies located ≈ 50 arcsec ( ≈ 250 kpc in projection) East of JADESGS-z7-LA (Fig. 15) where the surrounding photometric overdensity is very large ( ≈ 25) on small scales (see Fig. 17). It is conceivable  \nFigure 17. Measurements of photometric overdensity as a function of concentric angular radii around the two strong 𝑧 > 7 Ly 𝛼 emitters covered by the JADES/NIRCam imaging. Neither object shows a remarkably strong photometric overdensity across the 𝑧 ≈ 7 . 0 -7 . 6 redshift range probed by our photometric selection when centering the concentric rings on the LAEs themselves. However, there is a dense group of 𝑧 ≈ 7 . 0 -7 . 6 galaxies near JADES-GS-z7-LA (see Fig. 15) where the photometric overdensity is very high ( ≈ 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).  \n<!-- image -->  \nthat 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. 2023). The two UV-bright (M UV ≈ -20 . 3) systems located at a very similar projected distance to the South of JADESGS-z7-LA (Fig. 15) may be significantly contributing as well. While the surface density of 𝑧 ≈ 7 . 0 -7 . 6 galaxies measured in the ≲ 2 arcmin radius ( ≈ 600 kpc projected) around JADES-GS-z7-LA itself is consistent with the cosmic mean (see Fig. 17), this can easily rise by a factor of ∼ 2-3 if the nearby galaxies are largely located near JADES-GS-z7-LA along the line of sight (i.e. Δ 𝑧 ≲ 0 . 1).  \nIn contrast to the UV-faint Ly 𝛼 emitter described above, we do not identify any significant photometric overdensity near the UVbright (M UV = -21 . 5) 𝑧 = 7 . 51 Ly 𝛼 emitter z8\\_GND\\_5296 (Fig. 17). The surface density of 𝑧 ≈ 7 . 0 -7 . 6 objects on both small ( < 1.5 arcmin; < 450 kpc projected) and large (6-10 arcmin; ∼ 2 -3 Mpc projected) separations from z8\\_GND\\_5296 are consistent with the cosmic mean. This may indicate that the strong Ly 𝛼 emission from z8\\_GND\\_5296 is due more to its physical properties. With an estimated [OIII]+H 𝛽 EW of 1540 + 220 -160 Å, z8\\_GND\\_5296 lies in the high-end tail of the EW distribution for UV-bright 𝑧 ∼ 7 -9 galaxies (see Fig. 8). Accordingly, this source has a large inferred ionizing photon production efficiency of log [ 𝜉 ion /( erg -1 Hz )] = 25.73 + 0 . 06 -0 . 05 (from the TcSFH /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc SED fits) indicating that it has a relatively large intrinsic Ly 𝛼 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 Ly 𝛼 emitter.  \nThe detection of the systemic [CIII] 𝜆 1908 line from z8\\_GND\\_5296 may also indicate that its Ly 𝛼 photons are redshifted far into the damping wing (velocity offset ≈ 400 km/s; Hutchison et al. 2019) as commonly found among similarly bright 𝑧 ∼ 7 galaxies (see Endsley et al. 2022 and references therein), resulting in relatively low Ly 𝛼 opacity through the IGM (e.g. Stark et al. 2017; Mason et al.  \n2018). However, this is speculative as it is unclear which component of the [CIII] 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 𝛼 emission seen from this UV-bright Ly 𝛼 emitter, as well as measure the redshifts of the neighboring 𝑧 ≈ 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 𝑧 ∼ 6 -9 galaxies. This sample includes hundreds of reionization-era galaxies in the very UV-faint regime ( -18 ≤ M UV ≲ -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.  \n(i) The faintest galaxies in our sample (M UV ∼ -17) tend to have inferred stellar masses of (1-3) × 10 7 𝑀 ⊙ while those at the brightest end of our sample ( -22 ≲ M UV ≲ -21) have ∼ 100 × higher stellar masses (Fig. 2). Masses inferred from models with priors weighted towards extended SFHs (i.e. 'continuity' priors) are systematically ≈ 0.5 dex higher than those assuming CSFH, though can rise to factors of 10-100 × higher in the most extreme cases for galaxies with the youngest light-weighted ages (age CSFH ∼ 3 Myr; see Fig. 2g).  \n(ii) There are no galaxies in our sample where the photometric data clearly imply very large stellar masses ( > 3 × 10 10 𝑀 ⊙ ). There are only 13 galaxies in our sample with inferred stellar masses > 3 × 10 9 𝑀 ⊙ when adopting the continuity prior that favors extended SFHs (Fig. 3). Some of these objects are relatively faint (M UV ∼ -19) and show strong Balmer breaks consistent with old stellar populations. Others are simply so luminous (M UV ∼ -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 × 10 10 𝑀 ⊙ when including rest-frame nearinfrared SED constraints from parallel MIRI imaging (Fig. 4).  \n(iii) The typical galaxy in our sample shows a NIRCam SED consistent with a young light-weighted age of age CSFH ∼ 50 Myr (Fig. 5). However, we confidently identify several individual galaxies showing strong Balmer breaks implying much older ages (age CSFH ∼ 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 sub-population may be (temporarily) quenched.  \n(iv) We find a strong, highly-significant decline in the typical [OIII]+H 𝛽 EWsof 𝑧 ∼ 6 -9 galaxies towards lower UV luminosities (median EW ≈ 800 Å at M UV = -20 yet ≈ 350 Å at M UV = -17 . 5; see Fig. 8 and Table 2). We verify that this [OIII]+H 𝛽 EW trend with M UV is reflected in the typical NIRCam colors of our sample (Fig. 10). The [OIII]+H 𝛽 EW distribution is also found to broaden considerably at lower UV luminosities (standard deviation ≈ 0.3 dex at M UV = -20 yet ≈ 0.5 dex at M UV = -17 . 5). We infer a slight ( ≈ 0.05-0.1 dex) decline in the typical [OIII]+H 𝛽 EW between the 𝑧 ∼ 6 sample and the 𝑧 ∼ 7 -9 sample at fixed M UV .  \n(v) In remarkable contrast to the strong decline in [OIII]+H 𝛽 EW  \nwith M UV , we find that the H 𝛼 EWdistribution changes very weakly (if at all) across a factor of ≈ 10 in UV luminosity among our 𝑧 ∼ 6 sample (where H 𝛼 falls in F444W; Fig. 8). If we assume that the 𝑧 ∼ 6 H 𝛼 EWdistribution does not change with M UV over the range probed by our 𝑧 ∼ 6 sample, we infer a log-normal distribution with median EW=630 + 10 -40 Å and a standard deviation of 0.26 ± 0.01 dex.  \n(vi) We demonstrate that the [OIII]+H 𝛽 and H 𝛼 EW trends with M UV 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 𝑧 ∼ 6 galaxies in our sample ( ⟨ M UV ⟩ = -20 . 0) have a median inferred metallicity ≈ 3 × higher than that of the faintest 𝑧 ∼ 6 galaxies ( ⟨ M UV ⟩ = -20 . 0; 𝑍 ≈ 0 . 06 𝑍 ⊙ ). Moreover, the median inferred ratio of SFR averaged over the past 3 Myr to that averaged over the past 50 Myr (SFR 3 Myr / SFR 50 Myr ) is found to be ≈ 3.5 × higher among the brightest 𝑧 ∼ 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 3 Myr / SFR 50 Myr > 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.  \n(vii) There are two particular subsets of our sample which provide evidence in favor of bursty SFHs. A substantial fraction of galaxies in our sample ( ≈ 10-20% depending on M UV ) are inferred to exhibit extremely high [OIII]+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 sub-populations may therefore be caught during the peaks and trough of bursty star formation histories among reionization-era galaxies. We note that other recent studies are building empirical evidence for bursty SFHs among early galaxies (e.g. Dome et al. 2023; Dressler et al. 2023a; Looser et al. 2023a,b; Strait et al. 2023).  \n(viii) Wealso discuss how the [OIII]+H 𝛽 and H 𝛼 EWtrends with M UV 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 constant SFHs, we infer that the faintest 𝑧 ∼ 6 galaxies in our sample ( ⟨ M UV ⟩ = -17 . 5) are typically very efficient at leaking LyC photons into the IGM (median 𝑓 esc ≈ 0.5) while the brightest ( ⟨ M UV ⟩ = -20 . 0) objects leak only a moderate fraction of their LyC photons (median 𝑓 esc ≈ 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 [OIII]+H 𝛽 and H 𝛼 EW trends with UV luminosity found in this analysis.  \n(ix) Finally, we quantify the photometric overdensities around two strong Ly 𝛼 emitters at 𝑧 > 7 in the JADES footprint. The very UVfaint ( 𝑀 UV = -16 . 7) 𝑧 = 7 . 28 Ly 𝛼 emitter lies close to a very dense concentration of 𝑧 ≈ 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 ( 𝑀 UV = -21 . 5) Ly 𝛼 emitter 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 Ly 𝛼 emitters in the reionization era necessarily occupy large ionized bubbles.  \nTheJADESresultspresentedhereprovideavaluable 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; Hirschmann et al. 2022; Lewis et al. 2022; Seeyave et al. 2023). Such future work will help build insight into the nature and assembly of the faintest (and brightest) 𝑧 ∼ 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 & 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 & GCJ acknowledge funding from the "FirstGalaxies" Advanced Grant from the European Research Council (ERC) under the European Union\'s Horizon 2020 research and innovation programme (Grant agreement No. 789056) S.C acknowledges support by European Union\'s HE ERC Starting Grant No. 101040227 - WINGS. ECL acknowledges support of an STFC Webb Fellowship (ST/W001438/1) A.L.D. 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 number 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.  \nThis 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.  \nThis work made use of the following software: /n.pc/u.pc/m.pc/p.pc/y.pc (Harris et al. 2020); /m.pc/a.pc/t.pc/p.pc/l.pc/o.pc/t.pc/l.pc/i.pc/b.pc (Hunter 2007); /s.pc/c.pc/i.pc/p.pc/y.pc (Virtanen et al. 2020); /a.pc/s.pc/t.pc/r.pc/o.pc/p.pc/y.pc 8 , a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013; Price-Whelan et al.  \n2018); S/o.pc/u.pc/r.pc/c.pc/e.pc E/x.pc/t.pc/r.pc/a.pc/c.pc/t.pc/o.pc/r.pc (Bertin & Arnouts 1996) via /s.pc/e.pc/p.pc (Barbary 2016); /p.pc/h.pc/o.pc/t.pc/u.pc/t.pc/i.pc/l.pc/s.pc (Bradley et al. 2022); /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc (Chevallard & Charlot 2016); /m.pc/u.pc/l.pc/t.pc/i.pc/n.pc/e.pc/s.pc/t.pc (Feroz & Hobson 2008; Feroz et al. 2009); /p.pc/r.pc/o.pc/s.pc/p.pc/e.pc/c.pc/t.pc/o.pc/r.pc Johnson et al. (2021b); /d.pc/y.pc/n.pc/e.pc/s.pc/t.pc/y.pc (Speagle 2020); /s.pc/e.pc/d.pc/p.pc/y.pc (Johnson 2021); and /f.pc/s.pc/p.pc/s.pc (Conroy et al. 2009; Conroy & Gunn 2010) via /p.pc/y.pc/t.pc/h.pc/o.pc/n.pc-/f.pc/s.pc/p.pc/s.pc (Johnson et al. 2021a).', 'DATA AVAILABILITY': 'The HST ACS data used in this work are available from the Hubble Legacy Field 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': 'Figure A1. Cumulative distribution of the fraction of galaxies in our final Lyman-break sample (§2.3) as a function of the probability that their photometric redshift is 𝑧 > 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 sub-samples used in §4 separately.  \n<!-- image -->  \nphot  \n```\nWhitler L., Endsley R., Stark D. P., Topping M., Chen Z., Charlot S., 2023b, MNRAS, 519, 157 Whitler L., Stark D. P., Endsley R., Leja J., Charlot S., Chevallard J., 2023c, 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., 2023a, arXiv e-prints, p. arXiv:2301.09780 Williams C. C., et al., 2023b, arXiv e-prints, p. arXiv:2311.07483 Withers S., et al., 2023, arXiv e-prints, p. arXiv:2304.11181 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\n```', 'APPENDIX A: CONSIDERING A POTENTIAL POPULATION OF LOW-REDSHIFT CONTAMINANTS': 'As discussed in §2.3, we choose to select galaxies based on Lymanbreak color 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 longwavelength NIRCam color patterns not easily reproduced with lower redshift solutions. Nonetheless, we aim to quantify the potential lowredshift 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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc two-component SFH setup described in §2.4 and allow redshifts in the range 𝑧 = 0 -10 with a uniform prior.  \nFor nearly every sample subset considered in §4, ≳ 90%ofgalaxies  \nhave a high probability ( > 90%) of lying at 𝑧 > 4 (see Fig. A1). The only exception is the very faint (M UV > -18) subset of the F090W dropout sample, in which ≈ 80% of the galaxies have a large probability ( > 80%) of lying at high redshifts (Fig. A1).  \nIn §4.2, we found that the median [OIII]+H 𝛽 EWof reionizationera 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 each object when inferring the [OIII]+H 𝛽 EWdistributions by replacing Eq. 2 with the following:  \n𝑃 ( 𝜃 ) ∝ GLYPH<214> 𝑖      ∑︁ 𝑗 𝑃 𝑖, 𝑗 ( EW ) 𝑃 𝑗 ( EW | 𝜃 )      𝑃 𝑖 ( 𝑧> 4 ) (A1)  \nBy raising the term in brackets to the power of 𝑃 𝑖 ( 𝑧 > 4 ) (i.e. the probability that galaxy 𝑖 lies at 𝑧 > 4), we weight the contribution of each source to the EW distribution constraint based on the likelihood that it lies at high redshift.  \nWith Equation A1, we infer median [OIII]+H 𝛽 EWs of 400 + 70 -70 Å and 320 + 40 -50 Å among the very faint F775W and F090W dropout subsets, respectively, using the TcSFH model outputs. These median EWsare entirely consistent with our fiducial values reported in Table 2. This indicates that a potential moderate ( ≲ 10%) population of low-redshift interlopers in our faintest bins cannot explain the strong UV luminosity dependence on the median [OIII]+H 𝛽 found in §4.2. We have also verified that adopting Eq. A1 does not change our conclusions on the lack of a strong UV luminosity dependence on the H 𝛼 EWdistribution at 𝑧 ∼ 6 (§4.3).', 'APPENDIX B: ASSESSMENT OF THE GOODNESS OF SED FITS': 'Here we explore the goodness of our /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc SED fits which we utilize in our EW distribution analysis (§4). We show the best-fitting 𝜒 2 values of the F775W dropout ( 𝑧 ∼ 6) and F090W dropout ( 𝑧 ∼ 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 /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc 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 𝑧 ∼ 6 galaxies is discussed in detail in §5.2. Briefly, these models allow for solutions where the continuum is consistent with a young stellar population (age ≲ 30 Myr) while still yielding relatively weak [OIII]+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 is not sensitive to H 𝛼 emission at 𝑧 ≳ 7, which is necessary to determine if weak [OIII] emission can be explained by extremely low metallicity (see §5.2).  \nWe also check that the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc TcSFH models generally deliver acceptable fits to the rest-optical photometry alone. In Fig. B2, we show the distributions of best-fitting 𝜒 2 values limited to the restoptical bands, considering each redshift and UV luminosity subset adopted in §4 separately. We find that the best-fitting rest-optical 𝜒 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%) in every subset. Moreover, we find no strong differences in rest-optical fit quality among the faintest bins relative to the brighter bins.  \n<!-- image -->  \n50  \nFigure B1. A comparison of the best-fitting 𝜒 2 values (using all HST and NIRCam bands) between the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc TcSFH model fits (y-axis) and the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc CSFH fits (x-axis). We show the comparison for the F775W dropout ( 𝑧 ∼ 6) sample and the F090W dropout ( 𝑧 ∼ 7 -9) sample in the top and bottom panels, respectively.  \n<!-- image -->  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure B2. The distribution of best-fitting 𝜒 2 values from the /b.pc/e.pc/a.pc/g.pc/l.pc/e.pc TcSFH fits for each sample subset used in our [OIII]+H 𝛽 and H 𝛼 EW distribution analyses (§4). Here we are only considering the bands that probe the restoptical, so F200W (F277W) and redder for the F775W (F090W) dropout samples, resulting in six (five) fitted bands. The rest-optical photometry is well fit 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.  \n<!-- image -->'}