# CEERS: Increasing Scatter along the Star-Forming Main Sequence Indicates Early Galaxies Form in Bursts

Justin Cole, Casey Papovich, Steven Finkelstein, Micaela Bagley, Mark

Dickinson, Kartheik Iyer, L. Yung, Laure Ciesla, Ricardo Amorin, Pablo

Arrabal Haro, et al.

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

► To cite this version:

Justin Cole, Casey Papovich, Steven Finkelstein, Micaela Bagley, Mark Dickinson, et al. CEERS: Increasing Scatter along the Star-Forming Main Sequence Indicates Early Galaxies Form in Bursts. Astronomy & Astrophysics - A&A, 2025, 979 (2), pp.A205. 10.3847/1538-4357/ad9a6a. hal-04777890

HAL Id: hal-04777890

https://hal.science/hal-04777890v1

Submitted on 25 Mar 2025

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

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

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

# CEERS: Increasing Scatter along the Star-forming Main Sequence Indicates Early Galaxies Form in Bursts

Justin W. Cole¹,²,<fnref n="29" /> ID, Casey Papovich¹,² ID, Steven L. Finkelstein³ ID, Micaela B. Bagley³ ID, Mark Dickinson⁴ ID, Kartheik G. Iyer⁵ ID, L. Y. Aaron Yung⁶,⁷,<fnref n="30" /> ID, Laure Ciesla⁸ ID, Ricardo O. Amorín⁹,¹⁰ ID, Pablo Arrabal Haro⁴ ID, Rachana Bhatawdekar¹¹ ID, Antonello Calabrò¹² ID, Nikko J. Cleri¹,² ID, Alexander de la Vega¹³ ID, Avishai Dekel¹⁴ ID, Ryan Endsley³ ID, Eric Gawiser¹⁵ ID, Mauro Giavalisco¹⁶ ID, Nimish P. Hathi⁷ ID, Michaela Hirschmann¹⁷,¹⁸ ID,

Benne W. Holwerda¹⁹ ID, Jeyhan S. Kartaltepe²⁰ ID, Anton M. Koekemoer⁷ ID, Ray A. Lucas⁷ ID, Sara Mascia¹² ID, Bahram Mobasher¹³ ID, Pablo G. Pérez-González²¹ ID, Giulia Rodighiero²²,²³ ID, Kaila Ronayne¹,² ID, Sandro Tacchella²⁴,²⁵ ID, Benjamin J. Weiner²⁶ ID, and Stephen M. Wilkins²⁷,²⁸ ID

¹Department of Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242, USA; jwc68@tamu.edu

²George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA
³Department of Astronomy, The University of Texas at Austin, Austin, TX 78712-1205, USA

⁴NSF's National Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry Ave., Tucson, AZ 85719, USA

⁵Dunlap Institute for Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada

⁶Astrophysics Science Division, NASA Goddard Space Flight Center, 8800 Greenbelt Rd, Greenbelt, MD 20771, USA

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

⁸Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

⁹ ARAID Foundation, Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Unidad Asociada al CSIC, Plaza San Juan 1, E-44001 Teruel, Spain

¹⁰Departamento de Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La Serena 1720236, Chile

¹¹European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain

¹²INAF Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Rome, Italy

¹³Department of Physics and Astronomy, University of California, 900 University Ave, Riverside, CA 92521, USA

¹⁴Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem 91904, Israel

¹⁵Department of Physics and Astronomy, Rutgers, the State University of New Jersey, Piscataway, NJ 08854, USA

¹⁶University of Massachusetts Amherst, 710 North Pleasant Street, Amherst, MA 01003-9305, USA

¹⁷Institute of Physics, Laboratory for galaxy evolution, EPFL, Observatory of Sauverny, Chemin Pegasi 51, 1290 Versoix, Switzerland

¹⁸INAF, Astronomical Observatory of Trieste, Via G.P. Tiepolo 11, 34134 Trieste, Italy

¹⁹Physics & Astronomy Department, University of Louisville, Louisville, KY 40292, USA

²⁰Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA

²¹Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain

²²Department of Physics and Astronomy, Università degli Studi di Padova, Vicolo dell'Osservatorio 3, I-35122, Padova, Italy

²³INAF - Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122, Padova, Italy

²⁴Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK

²⁵Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK

²⁶MMT/Steward Observatory, University of Arizona, 933 N. Cherry St, Tucson, AZ 85721, USA

²⁷Astronomy Centre, University of Sussex, Falmer, Brighton, BN1 9QH, UK

²⁸Institute of Space Sciences and Astronomy, University of Malta, Msida MSD 2080, Malta

Received 2023 December 19; revised 2024 November 8; accepted 2024 December 2; published 2025 January 27

## Abstract

We present the star formation rate–stellar mass (SFR–M*) relation for galaxies in the Cosmic Evolution Early Release Science survey at 4.5 ≤ z ≤ 12. We model the JWST and Hubble Space Telescope rest-UV and rest-optical photometry of galaxies with flexible star formation histories (SFHs) using BAGPIPES. We consider SFRs averaged from the SFHs over 10 Myr (SFR$_{10}$) and 100 Myr (SFR$_{100}$), where the photometry probes SFRs on these timescales, effectively tracing nebular emission lines in the rest-optical (on ~10 Myr timescales) and the UV/optical continuum (on ~100 Myr timescales). We measure the slope, normalization and intrinsic scatter of the SFR–M* relation, taking into account the uncertainty and the covariance of galaxy SFRs and M*. From z ~ 5 to 9 there is larger scatter in the SFR$_{10}$–M* relation, with σ(log SFR$_{100}$) = 0.4 dex, compared to the SFR$_{100}$–M* relation, with σ(log SFR$_{10}$) = 0.1 dex. This scatter increases with redshift and increasing stellar mass, at least out to z ~ 7. These results can be explained if galaxies at higher redshift experience an increase in star formation variability and form primarily in short, active periods, followed by a lull in star formation (i.e., “napping” phases). We see a significant trend in the ratio R$_{SFR}$ = SFR$_{10}$/SFR$_{100}$ in which, on average, R$_{SFR}$ decreases with increasing stellar mass and increasing redshift. This yields a star formation “duty cycle” of ~40% for galaxies with log M*/M$_{⊙}$ ≥ 9.3 at z ~ 5, declining to ~20% at z ~ 9. Galaxies also experience longer lulls in star formation at higher redshift and at higher stellar mass, such that galaxies transition from periods of higher SFR variability at z ≳ 6 to smoother SFR evolution at z ≲ 4.5.

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

<footnote n="29">²⁹NASA FINESST Investigator.</footnote>

<footnote n="30">³⁰NASA Postdoctoral Fellow.</footnote>
*Unified Astronomy Thesaurus concepts: High-redshift galaxies (734); Galaxy evolution (594); Galaxy formation (595); Star formation (1569)*

# 1. Introduction

The beginning of galaxy formation and the early buildup of stellar mass in galaxies remains an unsolved area of astrophysical research. These are relevant both to understand the conditions that led to the onset of star formation in galaxies, and to understand when and how these galaxies reionized the intergalactic medium. Understanding early galaxy formation and stellar mass buildup has become even more pertinent in the era of JWST, where its data have identified via photometry, and confirmed via spectroscopy, galaxies at $z > 7$ and as high as $z \sim 12$ (see, e.g., S. L. Finkelstein et al. 2022; P. Arrabal Haro et al. 2023; E. Curtis-Lake et al. 2023; I. Labbé et al. 2023). In particular, JWST has progressed the study of these early galaxies to the point where we can characterize their properties from large, statistical samples with exquisite data.

From redshifts $z \sim 1$ to 6, the star formation history (SFH) of galaxies appears relatively smooth. The cosmic star formation rate (SFR) density grows from the earliest redshifts ($z \gtrsim 6$), peaks around $z \sim 2$, and is followed by a slow decline to the present, $z \sim 0$ (P. Madau & M. Dickinson 2014). This evolution is further evidenced by the evolution of the relation between the SFR and stellar mass of galaxies (the so-called “star-forming–main sequence,” or SFMS; K. G. Noeske et al. 2007; K. E. Whitaker et al. 2012; J. S. Speagle et al. 2014; B. Salmon et al. 2015; K. Iyer et al. 2018; E. Curtis-Lake et al. 2021; P. Popesso et al. 2023). Similarly, the interpretation of observational data has argued that the SFHs of individual galaxies follow this trend, albeit with a dependence on mass, such that more massive galaxies start and peak at earlier cosmic times (C. Papovich et al. 2011; C. Papovich et al. 2015; J. Pforr et al. 2012; N. A. Reddy et al. 2012; J. S. Speagle et al. 2014; C. Pacifici et al. 2016; K. G. Iyer et al. 2019). Therefore, from “cosmic noon” (the peak of the cosmic SFR density at $z \sim 2$) to the highest redshifts we have explored, galaxy formation has appeared relatively smooth and orderly for galaxies, at least as a population.

However, simulations of galaxy formation find that, while on average the SFHs of galaxies are smooth (e.g., K. Finlator et al. 2011; S. M. Wilkins et al. 2017, 2023; S. Bird et al. 2022), they are much more variable on shorter timescales, with larger variations in the SFR (K. Finlator et al. 2006; R. Davé et al. 2016; X. Ma et al. 2016; S. Tacchella et al. 2016; M. Donnari et al. 2019; H. Katz et al. 2023; G. Sun et al. 2023). G. Sun et al. (2023) recently illustrated this with simulations showing that galaxies at higher redshifts experience greater variability in their specific SFRs (sSFR $\equiv$ SFR/M$_*$). In these simulations, this stems from the fact that galaxies are more compact, with higher gas and stellar-mass surface densities. This leads to higher feedback (per unit mass) from supernovae and winds from massive stars, which more readily expel gas and disrupt gas accretion. These feedback mechanisms subsequently create a lull in the star formation. This naturally predicts more galaxies in “burst” modes of star formation (with elevated sSFR) and galaxies with temporary cessation of star formation (so-called “napping,” or “mini-quenched,” galaxies; T. J. Looser et al. 2023, 2024; T. Dome et al. 2024). The timescales of this feedback relative to the star formation and its

amplitude (i.e., strength per unit mass) have a strong impact on the SFR and sSFR (e.g., S. R. Furlanetto & J. Mirocha 2022).

With JWST, there has been renewed interest in the SFR variability in high-redshift galaxies. L. Whitler et al. (2023) and R. Endsley et al. (2024) argued that higher variability in the SFRs of galaxies is needed to explain the increase in the frequency of galaxies with evidence for strong emission-line equivalent widths (EWs). Other recent work from JWST imaging of low-mass ($10^8 M_\odot$) galaxies at $z \sim 5-6$ shows that $\sim 60\%$ of them have $H\alpha$/UV-flux ratios which indicate they are experiencing bursts and variability in their SFHs (S. Tacchella et al. 2023; Y. Asada et al. 2024). Yet other studies have invoked a higher variability in the SFRs of high-redshift galaxies to explain the higher-than-expected number density of galaxies measured at the bright end of the rest-UV luminosity functions (i.e., S. L. Finkelstein et al. 2023; C. A. Mason et al. 2023). This is also seen in updated simulations of galaxy formation (e.g., X. Shen et al. 2023; S. M. Wilkins et al. 2023), and inferred from studies of the evolution of galaxies along the SFR–$M_*$ relation (the SFR “gradient”; e.g., L. Ciesla et al. 2024).

Motivated by these observations and apparent problems, here we study the evolution of the SFR–$M_*$ plane using JWST data from the Cosmic Evolution Early Release Science (CEERS) survey. CEERS is a JWST survey of the Extended Groth Strip (EGS; O. Ilbert et al. 2013), providing seven photometric bands with NIRCam covering $\sim 1-5 \mu m$. This field also has deep legacy imaging from the Hubble Space Telescope’s (HST) Advanced Camera for Surveys (ACS) and Wide Field Camera 3 (WFC3; $\sim 0.6-1.6 \mu m$) and other ground-based and space telescopes (see N. A. Grogin et al. 2011; A. M. Koekemoer et al. 2011; M. Stefanon et al. 2017). We model the JWST/NIRCam and HST multiwavelength photometry with stellar population models using flexible SFHs to derive constraints on the stellar masses, SFRs, and redshifts. We then use the full posteriors from these fits to study the evolution of the slope, normalization, and intrinsic scatter of the SFR–$M_*$ relation (the SFMS), and the scatter of the SFR as a function of galaxy mass and cosmic time (redshift). We invoke a novel method to measure the evolution of the SFMS that accounts for the uncertainties and covariance from the posteriors of all the variables.

The outline for the rest of this paper is as follows. Section 2 presents the CEERS data used for this analysis and discusses our selection methods. Section 3 provides details of our analysis methods used to model the galaxy spectral energy distributions (SEDs) with stellar population models and flexible SFHs. We also discuss the ability of the modeling to distinguish between the SFRs on timescales of 10 and 100 Myr. In Section 4, we present measurements of the SFMS and its evolution. We provide a discussion of the implications of our results and compare the results to theoretical predictions in Section 5. In Section 6, we present our conclusions. We also include an Appendix A with tests of our ability to fit and recover the scatter of the SFMS using our methods. Throughout this paper, we assume a flat Lambda cold dark matter ($\Lambda$CDM) cosmology with values from Planck ($H_0 = 67.36$, $\Omega_M = 0.3158$). All magnitudes are presented in the AB system (J. B. Oke & J. E. Gunn 1983; M. Fukugita et al. 1996). We
assume a G. Chabrier (2003) initial mass function (IMF) for derived quantities such as the SFR and stellar mass.

## 2. Observations and Data Selection

We leverage the CEERS data set (S. L. Finkelstein et al. 2022; M. B. Bagley et al. 2023), providing JWST observations in the CANDELS (A. M. Koekemoer et al. 2011; N. A. Grogin et al. 2011) EGS field. The data set includes NIRCam imaging from JWST, and ACS and WFC3 imaging from HST. The detailed description of the NIRCam data reduction and methods for extracting the photometry can be found in M. B. Bagley et al. (2023).

Our work utilizes the 10 NIRCam pointings obtained by CEERS, with imaging in the F115W, F150W, F200W, F277W, F356W, F410M, and F444W filters, spanning a wavelength range from approximately 1–5 µm. The HST data includes imaging from the ACS F606W and F814W bands, and the WFC3 F125W, F140W, and F160W bands (with some data in F105W also available), covering 0.6–1.6 µm. The NIRCam data are aligned astrometrically to the HST imaging (see S. L. Finkelstein et al. 2023).

We adopt the photometric catalog from S. L. Finkelstein et al. (2024). Galaxies are detected using the F277W imaging, and photometry is measured using matched apertures applied to the NIRCam images and HST ACS and WFC3 images. The catalog then includes photometric data points from the 12 bandpasses spanning a wavelength range of 0.6–5 µm for ≈ 100,000 sources. In addition, we also make use of HST/WFC3 F105W band imaging that covers approximately 20% of the CEERS/NIRCam field. This provides an additional photometric point at ~1 µm for ~20% of the galaxies in the catalog. Below, we discuss the selection criteria applied to construct the sample used in this study (Sections 2.1 and 2.2) with additional explanations for sample selection provided in Section 3.

### 2.1. The CEERS 4.5 ≤ z ≤ 12 Sample

We selected galaxies with $m_{F277W} \le 28.2$, which is 1 magnitude brighter than the 5σ, point-source detection limit ($\approx 29.2$ mag; M. B. Bagley et al. 2023; S. L. Finkelstein et al. 2023). This ensures our sample has high fidelity through the benefit of higher signal-to-noise ratio (SNR) in multiple photometric bands. We also require galaxies to have photometric redshifts between $z_{peak} \sim 4.5$ and $z_{peak} \sim 12$, using redshifts ($z_{peak}$) derived from EAZY (G. B. Brammer et al. 2008). A description of the EAZY-derived redshifts can be found in S. L. Finkelstein et al. (2024). The lower redshift bound, $z > 4.5$, ensures that the photometric catalog covers the wavelength region around the Lyman limit/Lyα breaks for all galaxies, allowing for more accurate photometric redshifts. At the upper redshift bound, $z \le 12$, the photometric catalog covers wavelengths around the Balmer limit break (3646 Å) and the 4000 Å break (from line blanketing of metals in stellar photospheres). Therefore, our selection of $4.5 < z_{peak} \le 12$ provides rest-frame far-UV (FUV) to rest-optical photometry for the full sample, and this enables higher fidelity in derived quantities, such as the stellar mass and SFR estimates from SED fitting (J. Pforr et al. 2012). Applying the magnitude limit and redshift range selection to the catalog yields a parent sample of 5669 sources. We then inspected this sample visually and identified several sources of potential contamination, which

we screened through the application of additional selection criteria, described in Section 2.2 below.

The choice of JWST/NIRCam F277W as our selection bandpass is motivated by two primary factors. First, JWST/NIRCam F277W is the most sensitive bandpass. We demonstrate in Appendix B that selecting in the reddest band, JWST/NIRCam F444W, reduces the number of galaxies in our study by a factor of ~2, primarily affecting the uncertainties on inferred parameters presented in Section 4. However, the results are consistent within $1\sigma$, so this increase in uncertainties on these parameters does not change the qualitative conclusions drawn from this study. Second, JWST/NIRCam F277W is relatively unaffected by the presence of strong emission lines except for at redshifts $4.5 < z < 5.0$. We address these potential effects further in Section 4.1, where we estimate the effect of strong optical emission lines on the stellar-mass completeness function used in this study.

### 2.2. Matching Hubble Space Telescope Photometry, and Outlier Removal

We use three additional selection criteria to flag and exclude sources that stem from photometric contamination effects, and we also remove sources that have contamination from possible active galactic nuclei (AGN). First, because galaxies are matched to the legacy HST imaging, we require that the photometric measurements between the HST data and JWST data be consistent. Specifically, we require

$$ \frac{|m(\text{F150W}) - m(\text{F160W})|}{\sqrt{\sigma_{\text{F150W}}^2 + \sigma_{\text{F160W}}^2}} \le 2.5, \quad (1) $$

where $m(\text{F150W})$ ($m(\text{F160W})$) and $\sigma_{\text{F150W}}$ ($\sigma_{\text{F160W}}$) are the measured photometry and uncertainty in the JWST/NIRCam F150W (HST/WFC3 F160W) bandpass. This requirement forces the HST/F160W photometry to be within $2.5\sigma$ of the JWST/F150W photometry. In practice, we found that this selection requirement flags and removes obvious contamination between foreground sources, often because of differences in the point-spread function (PSF) between HST and JWST (e.g., objects with optical diffraction spikes). This also flagged and removed several spurious imaging artifacts and objects near the edges of the images. This selection removed 1231 galaxies (21% of 5669) from the initial parent sample.

Second, to ensure a robust, high-quality sample of galaxies, we added a requirement that the galaxies in our sample have well-defined photometric redshifts. Specifically, we exclude objects with multiple “peaks” in the photometric redshift probability density functions (PDFs, $P(z)$) or that have a very broad $P(z)$. In such situations, the uncertainty in the redshift $P(z)$ translates to very large uncertainties in the derived physical parameters like masses and SFRs (see, e.g., B. Mobasher et al. 2007; B. Darvish et al. 2015; B. Mobasher et al. 2015). As the goal of our study is to understand the latter, we exclude objects where the photometric-redshift accuracy would be the limiting factor. To make this selection, we defined a redshift, $P(z)$ “WIDTH_FLAG” using the EAZY redshift $P(z)$ functions for each source in the catalog. We start by using a method similar to the highest-density interval (e.g., C. A. L. Bailer-Jones et al. 2018). We start at the peak probability of the PDF and move away from the peak by stepping through redshift space until the area under the $P(z)$ curve integrates to 0.7 (roughly analogous
to the inner 16th to 84th percentiles). This then defines a $z_{low}$ and $z_{high}$ that encompass 70% of the probability distribution. We require that the difference between $z_{low}$ and $z_{high}$ meets the criteria

$$ \frac{z_{\text{high}} - z_{\text{low}}}{(1 + z_{\text{peak}})} \le 0.5. \qquad (2) $$

This ensures that we are only using galaxies that have relatively well-defined redshift probability spaces because 70% of the integrated $P(z)$ lies within $0.5 \times (1 + z_{peak})$.

This step removes 46% of the initial sample; however, it leaves galaxies with much more robust photometric redshifts, which translates to a higher-quality sample of objects. The galaxies removed by Equation 2 have poor constraints on $P(z)$, which leads to large uncertainties on the stellar mass and SFR constraints. We include the full probability density for each galaxy in the modeling below, which causes these galaxies to contribute very little signal to the SFR-$M_*$ relation. We therefore opt to remove these galaxies as they contribute more noise than signal. Nevertheless, we have also verified that our results and conclusions are unaffected by this decision. If we include the objects rejected by Equation 2, all of the median values for the parameters on the SFMS change negligibly, and the estimated uncertainties change by $\lesssim 15\%$ level. This supports our assertion that the objects with poor $P(z)$ constraints also have poor constraints on quantities such as stellar mass and SFR (similar findings are published in, B. Mobasher et al. 2007; B. Darvish et al. 2015; B. Mobasher et al. 2015). For this reason we present results using the sample that excludes objects that do not satisfy Equation 2, as the latter add more noise than signal, but our results are otherwise robust to this selection.

Lastly, we also reject galaxies that may contain an AGN. We used a Chandra-based X-ray source catalog (K. Nandra et al. 2015) of the EGS field to determine if any galaxies in our sample have high X-ray luminosities ($L_X/\text{erg s}^{-1} \ge 3 \times 10^{42}$, Y. Q. Xue et al. 2016; B. Luo et al. 2017). We matched these sources to our catalog using a search radius of $1''$, which yield only 3 sources with $L_X \ge 3 \times 10^{42} \text{ erg s}^{-1}$, which we excluded from the sample. Although this does not account for dust-obscured AGN, or faint AGN, by rejecting X-ray-bright AGN, this improves the likelihood that the galaxies in our sample have rest-UV-optical SEDs that are dominated by their stellar light and/or nebular emission from star-forming regions.

Following all the selection criteria, our final sample includes 1863 galaxies. We then proceed to model the galaxies' SEDs using stellar populations fits with the BAGPIPES code.

### 3. Spectral Energy Distribution Fitting

#### 3.1. BAGPIPES: Priors and Parameter Spaces

We fit the SEDs of all sources in our sample using BAGPIPES (v1.0.1; A. C. Carnall et al. 2018), a Bayesian SED fitting code designed to fit photometry and/or spectra of galaxies with a broad range of parameter spaces. BAGPIPES utilizes a $\chi^2$-likelihood function, assuming that photometric errors follow a normal distribution, and allows for a wide range in the prior parameter spaces. BAGPIPES then samples the posteriors for model parameters using the MULTINEST nested sampling algorithm (see F. Feroz et al. 2009; A. C. Carnall et al. 2018). Table 1 provides a summary of the parameters, their range, prior, and other information.

<table><caption>**Table 1**
Parameters, Parameter Spaces, and Priors Used in Our Custom Version of BAGPIPES</caption><thead><tr><th>Parameter</th><th>Parameter Space (Prior)<sup>a</sup></th></tr></thead><tbody><tr><td>SFH: Gaussian process</td><td>Four time bins: t<sub>X</sub></td></tr><tr><td>Dense Basis</td><td>(Dirichlet continuity)</td></tr><tr><td>(K. G. Iyer et al. 2019)</td><td>log(SFR) ∈ [-3, 4]<br>(Uniform in log<sub>10</sub>)</td></tr><tr><td>Stellar mass:</td><td>log(M<sub>*</sub>/M<sub>☉</sub>) ∈ [6, 15]<br>(Uniform)</td></tr><tr><td>Redshift:</td><td>z ∈ [0, 20]<br>(EAZY PDF)</td></tr><tr><td>Ionization parameter:</td><td>log(U) ∈ [-4, -1]<br>(Uniform)</td></tr><tr><td>Metallicity:</td><td>Z/Z<sub>☉</sub> ∈ [0.001, 1], Z<sub>stellar</sub> = Z<sub>neb</sub><br>(Uniform)</td></tr><tr><td>D. Calzetti et al. (2000) dust law:</td><td>A<sub>V</sub>/mag ∈ [0, 6]<br>(Uniform)</td></tr></tbody></table>

**Note.**

With BAGPIPES, we fit all galaxies assuming a G. Chabrier (2003) IMF and we allow BAGPIPES to fit for $log(M_*/M_☉) \in [6, 15]$ with a uniform prior, providing a wide prior space for stellar-mass measurements. We also assume a D. Calzetti et al. (2000) dust law, which has been suggested to be representative of high-redshift galaxies (R. A. A. Bowler et al. 2024). We employ Binary Population and Spectral Synthesis (BPASS, v2.2.1; J. J. Eldridge et al. 2017) stellar templates, and fit over a range of metallicity ($Z/Z_☉$) ∈ [0.001, 1]. We adopt the BPASS models that include the effects of mass transfer between binary stars, with an upper limit of 300 $M_☉$ to better reproduce the emission-line strength and ratios seen in some galaxies (e.g., C. C. Steidel et al. 2018; A. L. Strom et al. 2017; G. M. Olivier et al. 2022; R. L. Larson et al. 2023). We further allow for nebular emission (modeled with CLOUDY; G. J. Ferland et al. 2017), where the metallicity of the nebular gas is equal to that of the stellar population ($Z_{stars} = Z_{neb}$). The ionization parameter ranges over $log(U) \in [-4, -1]$ and we use a log-uniform prior.

For the SFH parameterization in BAGPIPES, we adopt the Gaussian process (GP) model from Dense Basis (K. G. Iyer et al. 2019). For a more detailed description of the GP implementation, see A. C. Carnall et al. (2019). In this case, the SFH is defined by $N$ parameters, $\{t_X\}$, where $X = 1, \dots, N$. Each $\{t_X\}$ is the lookback time at which a fraction, $1/N$, of the total stellar mass is formed. The value for each $t_X$ is allowed to vary in the fit to force this condition. K. G. Iyer et al. (2019) has demonstrated the robustness of this methodology in reconstructed bursty SFHs, which may be more indicative of star formation in the early Universe (T. Kimm et al. 2015; D. Ceverino et al. 2018; K. S. S. Barrow 2020; S. R. Furlanetto & J. Mirocha 2022; R. Kannan et al. 2022; G. Sun et al. 2023; R. Endsley et al. 2024).

For this work, we use $N = 4$ and fit for $\{t_X\}$ such that each lookback time bin contains a quarter (25%) of the total stellar mass formed, and we set a Dirichlet continuity prior (see J. Leja et al. 2017, 2019; K. G. Iyer et al. 2019, for a more detailed description) on the SFH, forcing continuity across the time bins.

<footnote><sup>a</sup> Parameter spaces are shown in brackets, and the prior imposed on each parameter space is shown in parentheses below each parameter space.</footnote>
Regarding the choice of $N = 4$, as we discuss in Section 3.2, the photometric data are sensitive to star formation on roughly three different timescales, which sets $N \ge 3$: (1) the data measuring the strength of nebular emission (e.g., [O III], $H\beta$, $H\alpha$) probe < 10 Myr timescales (the lifetime of O-type stars); (2) the data measuring the rest-UV probe ~100 Myr timescales (e.g., R. C. Kennicutt 1998); (3) the data measuring the region of the 4000 Å/Balmer break probe timescales of >500 Myr (the lifetimes of A-type and later-type stars). We set $N = 4$ to include these three time bins, and to allow for earlier episodes of star formation lost in the glare of the current stars (e.g., the “outshining” effect; C. Conroy 2013; C. Papovich et al. 2023; S. Tacchella et al. 2023; D. Narayanan et al. 2024). We have tested our methodology using $N = 3$ and $N = 6$ with negligible differences between inferred parameters, and therefore this choice does not impact our results. For the individual time bins, we employ a log-uniform prior for $log(SFR/M_\odot/yr^{-1}) \in [-3, 4]$, which allows for a wide range of SFRs for each time bin.

Modifying the method presented in K. Chworowsky et al. (2024), we have adapted BAGPIPES to use the $P(z)$ distributions as a prior on the redshift for each galaxy (this method was also implemented by C. Papovich et al. 2023). As mentioned in Section 2.2, we use the $P(z)$ for each source measured with EAZY. We opt to use this methodology as it provides a more physically motivated, and therefore more meaningful, distribution from which to draw redshift probabilities.

We performed the SED fitting using BAGPIPES on the Texas A&M High Performance Research Computing (HPRC) Grace cluster,<fnref n="31" /> which has 800 Intel 6248R 3.0 GHz 24-core processors with (at least) 384 Gb RAM. On a single core, the clock time to fit each galaxy ranged from 4 to 10 minutes. The total clock time to fit the full sample was ~520 core-hours (though the galaxies were fit in parallel so the wall-clock time was substantially shorter).

Figure 1 shows the results from the fits to one of our galaxies, CEERS_4386, which has a stellar mass representative of the mean stellar mass of the sample and high SFR compared to galaxies at $z \sim 5$. The corner plot (panel (a)) shows the posterior space for the stellar mass (accounting for mass loss due to stellar evolution), SFR, dust attenuation ($A_V$), and redshift using two-dimensional kernel density estimates (KDEs, where darker colored layers show more densely populated regions of the posterior spaces). The inset panels of Figure 1 show the fitted model SED (panel (b)) and the constraints on the SFH (panel (c)). Importantly, the thick blue and red lines in the latter show the SFR for the galaxy derived by averaging the SFH over the previous 10 and 100 Myr timescales, respectively. This will become important as the SED fitting is sensitive to features that vary on these timescales: The strength of the nebular emission responds to changes in the SFR over timescales of ~5 Myr (the lifetime of a low-metallicity, 25 $M_\odot$ O star), while the strength of the rest-frame UV continuum responds to changes in the SFR over timescales of ~20–100 Myr (the lifetime of stars more massive than 10 $M_\odot$ B stars; J. J. Eldridge & E. R. Stanway 2022).

The length ($\Delta t$) of the red and blue boxes demonstrates the timescale over which each SFR is averaged (100 and 10 Myr, respectively). The location of the boxes along the ordinate (SFR) shows the mean SFR when measured over each of these

10 and 100 Myr timescales. One of our conclusions is that the SED fits are able to differentiate the SFR on these timescales, and thereby place constraints on them. This impacts the interpretation of the SFR-$M_*$ relation. Before exploring the SFR-$M_*$ relation for the entire population of galaxies, we discuss how the photometric SEDs enable us to differentiate between galaxies with SFHs that are rising or declining over the past ~10 Myr compared to ~100 Myr timescales.

## 3.2. Star Formation Rates and Inference of Spectral Features

From the SED fitting, we are able to measure differences in the SFH when it is averaged over the previous 10 Myr ($SFR_{10}$) from those averaged over the previous 100 Myr ($SFR_{100}$), using the definitions below. There are several reasons for this. First, we have many photometric bands covering the galaxies’ SEDs from the redshifted $Ly\alpha$ break to the Balmer/4000 Å break. These are sensitive to the relative number of massive stars (B type and earlier, whose continua dominate the rest-UV) compared to less massive stars (A type and later), which contribute to the Balmer/4000 Å break. Second, these bandpasses are sensitive to the effects of emission lines in the rest-optical. This includes primarily [O III] $\lambda\lambda4960, 5008 + H\beta$, and $H\alpha$. These lines have the high EWs in star-forming galaxies (e.g., R. Smit et al. 2016; D. Sobral et al. 2018; C. C. Steidel et al. 2018; M. V. Maseda et al. 2023), and their strength appears to increase at redshifts $z \sim 2-7$ (see, e.g., R. Endsley et al. 2021; M. Stefanon et al. 2022; K. N. K. Boyett et al. 2022; R. Endsley et al. 2024). Because the EWs in the observed frame increases, the emission lines impact the fluxes in the medium- and broadband filters to a higher degree (e.g., C. Papovich et al. 2001). Importantly, because stellar lifetimes increase with decreasing mass, all of these features (nebular emission lines, the UV continuum, Balmer/4000 Å breaks) are sensitive to the SFH on different timescales (~3–10 Myr, ~20–100 Myr, and $\gtrsim 500$ Myr, respectively). This allows us to use the multiband data set to constrain the SFH.

As an example, these features are seen in the SED of CEERS_4386 in Figure 1(b), where both the JWST/NIRCam F277W and F356W photometric bands have elevated fluxes. Our interpretation from the SED modeling is that this is due to the presence of the [O III] + $H\beta$ and $H\alpha$ emission lines in these bands, respectively. For galaxies at $4.5 \lesssim z \lesssim 6$, the long-wavelength NIRCam bands contain $H\alpha$, $H\beta$, and [O III], while the shorter-wavelength bands contain the rest-frame UV coverage and the 4000 Å/Balmer breaks. For these galaxies we can place the tightest constraints on the SFH averaged over 10 and 100 Myr. For galaxies at $6 \lesssim z \lesssim 9$, all of these features remain accessible to NIRCam, except for $H\alpha$, which shifts to higher wavelengths. At $z \gtrsim 10$, we lose multiwavelength coverage of the rest-frame optical emission lines and the 4000 Å region in the NIRCam coverage, but we retain multiwavelength coverage of the rest-UV with the HST/ACS imaging. Motivated by these effects, we consider what features in the observed galaxy SEDs drive differences in the SFRs derived from the SFH when averaged over 10 and 100 Myr timescales.

Throughout, we define two quantities, $SFR_{10}$ and $SFR_{100}$, which are the SFR derived by averaging the SFH over 10 Myr and 100 Myr, respectively. We then derive these for every

<footnote n="31">³¹ See https://hprc.tamu.edu/kb/User-Guides/Grace, where the namesake of the cluster is https://en.wikipedia.org/wiki/Grace_Hopper.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** An example of the posterior space and best-fit models for a representative galaxy in our sample (CEERS_4386). The panels in (a) are the corner plot, which shows the two-dimensional KDEs of the multidimensional posterior space for stellar mass, the SFR (averaged over the previous 10 Myr, SFR<sub>10</sub>), dust attenuation (A<sub>V</sub>), and redshift estimated from the model fits. Panel (b) shows the median, best-fit SED as a solid black line with the inner-68th percentile as the gray region. The red squares and error bars show the measured photometry and uncertainties. The median redshift, stellar mass, and SFR<sub>10</sub> for this galaxy are listed in the panel inset. Panel (c) shows the inferred SFH, with the median SFH as a solid black line and the inner-68th percentile as the gray-shaded region. The blue and red rectangles show the SFR averaged over the previous 10 and 100 Myr, respectively (SFR<sub>10</sub> and SFR<sub>100</sub>, respectively).</figcaption></figure>

galaxy in the sample using the constraints on the SFH using

$$SFR_{10} = \frac{1}{10 \text{ Myr}} \int_{t_{\text{obs}}}^{t_{\text{obs}+10 \text{ Myr}}} SFR(t') \, dt', \quad (3)$$

$$SFR_{100} = \frac{1}{100 \text{ Myr}} \int_{t_{\text{obs}}}^{t_{\text{obs}+100 \text{ Myr}}} SFR(t') \, dt', \quad (4)$$

where in each case the integral spans a time from the current age of the Universe at the redshift of the galaxy (t<sub>obs</sub>) to the previous 10 or 100 Myr, respectively.

In Figures 2–4, we show examples of galaxies that have high, “normal,” and low SFR<sub>10</sub> relative to the median SFR for the galaxy sample at a given redshift. We note that for all galaxies there is less of a constraint at earlier ages, where all the SFH constraints allow for a possible, low level of star
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Example SEDs for galaxies with $log(M_*/M_\odot)$ ~9 and $log(SFR_{10}) \gtrsim 1.0$, such that they lie above the median SFR at their stellar mass compared to the SFMS. Each row shows one galaxy selected from each of the five redshift bins (increasing redshift from top to bottom). In the left panel of each row, the red squares and error bars show the observed photometry and uncertainties. Photometric points with SNR <1 are shown as black, downward-facing triangles to indicate upper limits. The median best-fit model SED for each galaxy is shown as the black solid line and the inner-68th percentile range as the gray-shaded region. The middle panel in each row shows the model SFH with the median (black line) and inner-68th percentile range for each galaxy. In the right panel of each row, we show the location of each galaxy (large square) in the SFR-$M_*$ plane compared to the other galaxies (small symbols) in the same redshift bin. In all cases (especially for $z < 9$), the galaxy SEDs show indications of strong emission lines, which causes an increase in the SFH in the previous ~10 Myr compared to the past average.</figcaption></figure>

formation in the distant past. This star formation is certainly not excluded by the modeling as the light from the stars would be lost in the glare of the current generation of stars (this is the well-known “outshining” effect; C. Papovich et al. 2006; C. Conroy 2013). C. Papovich et al. (2023) showed that even with JWST data one cannot exclude these early-formed stars (and certainly *some* earlier generations of stars must have

existed given the propensity of emission lines from metals in the galaxies). As we are interested in the current SFR averaged over the previous 10 and 100 Myr, possible earlier generations of stars in our galaxies do not change our conclusions here.

Figure 2 shows examples of galaxies that have high $SFR_{10}$. These galaxies all have $SFR_{10}$ values $\approx 0.3$ dex higher than the median $SFR_{10}$ of galaxies in the same redshift bin with
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Same as in Figure 2, but here we show example galaxies that have SFRs within 0.1 dex of the median SFR at their redshift. In these cases, the SFR over the past 10 Myr (i.e., SFR<sub>10</sub>) is consistent with the past average over the previous 100 Myr (SFR<sub>100</sub>).</figcaption></figure>

log(M*/M⊙) ∈ [8.7, 9.4]. Each galaxy shows evidence in its photometric SED for the presence of emission lines boosting the observed flux, specifically where bands contain the redshifted [O III] + Hβ and the Hα emission lines. The SED fits respond to these features, where the models that reproduce the data favor a recent “upturn” in the SFR that leads to a higher SFR<sub>10</sub>. This “upturn” in the SFR is not always sufficient to increase SFR<sub>100</sub> when the SFH is averaged over the longer timescale. The presence of these strong emission features allows BAGPIPES with the Dense Basis SFHs to infer that there is a recent increase in the SFR in these galaxies. Quantitatively, this increases the likelihood

for higher SFR<sub>10</sub> values compared to SFR<sub>100</sub>. In the examples in Figure 2, the SFR in the past 10 Myr is higher by log(SFR<sub>10</sub>) − log(SFR<sub>100</sub>) ≥ +0.3 dex, 1σ above the mean SFR at fixed stellar mass and their respective redshifts, with the exception of the highest-redshift bin.

Figure 3 shows examples of galaxies that have SFR<sub>10</sub> within ≈ 0.3 dex of the median SFR<sub>10</sub> at their redshift and with log(M*/M⊙) ∈ [8.7, 9.4]. In this sample of galaxies there is evidence for elevated flux in the photometry due to the presence of emission lines, and the blue slopes of the SEDs in the UV indicate ongoing, or very recent (≈ 100 Myr), star
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Same as in Figure 2, but for galaxies with $log(SFR_{10}/M_{\odot} yr^{-1}) \lesssim 0$, below the median SFR of the SFR–$M_{*}$ plane. In these cases, the lack of evidence for strong emission lines in the SEDs implies the SFR has declined over the previous ~10 Myr compared to the past average.</figcaption></figure>

formation. However, the emission-line strength is generally weaker, with lower EWs compared to the galaxies in Figure 2. There are also cases where the galaxy SEDs show evidence for Balmer breaks, implying the fits favor models with a larger fraction of evolved stellar populations. Taken together, the SED fits favor models where there is little difference between $SFR_{10}$ and $SFR_{100}$. On average, we find that $log(SFR_{10}) – log(SFR_{100}) \simeq 0.0$ dex for these galaxies within their $1\sigma$ uncertainties.

Lastly, Figure 4 shows examples of galaxies that have $SFR_{10} \approx 0.3$ dex lower than the median SFR of galaxies with

similar redshifts and $log(M_{*}/M_{\odot}) \in [8.7, 9.4]$. Among this sample, we observe little or no evidence for the presence of emission lines in the photometry of these galaxies. This implies little recent star formation in the past ~10 Myr, which is favored by the model fits as shown in the constraints on the SFHs in the figure. However, the SEDs *do* indicate evidence for Balmer/4000 Å breaks, which leads to higher stellar-mass-to-light ratios and lower SFRs. Taken together, the SED fits therefore favor SFHs where the $SFR_{10}$ is lower than $SFR_{100}$. For these galaxies we find $log(SFR_{10}) – log(SFR_{100}) \lesssim -0.3$ dex, $1\sigma$ below the mean $SFR_{10}$ at fixed mass and redshift.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** The SFR–$M_*$ relation illustrated as the star-forming main sequence (SFMS) for the full sample of galaxies above our estimated stellar-mass completeness limits. Data points show these galaxies, color-coded by redshift. Faded points show galaxies below the 90% completeness limits at their respective redshifts, which are indicated by the colored vertical lines at the bottom of each panel. The error bars indicate the 68% uncertainty on each point. The uncertainties on mass and SFR are covariant, and the gray-shaded ellipses in the left panel show representative $1\sigma$ and $2\sigma$ posterior ellipses in SFR–$M_*$ space for a single galaxy (the size of the ellipse may change, but the orientation of the covariance is similar for all galaxies in both SFR$_{10}$ and SFR$_{100}$). In the left panel, we show the SFR–$M_*$ relation using the SFRs measured over 100 Myr (SFR$_{100}$), and in the right panel, we show the same but with SFRs measured over 10 Myr (SFR$_{10}$). We observe little evolution in either the slope or normalization with redshift in the SFMS, but the scatter in the SFMS increases when using the SFR$_{10}$ values.</figcaption></figure>

Therefore, the SED fitting is able to provide meaningful constraints on the SFH over the past 10 and 100 Myr. This information is important for interpreting the evolution of the SFMS at these redshifts. This is primarily true for galaxies at $4.5 < z < 9$, where (as mentioned above) the NIRCam photometry samples the full SED from the rest-frame UV through the rest-frame optical, capturing the UV continuum, important nebular emission lines in the rest-optical, and the region around the Balmer/4000 Å break. In contrast, at $z > 9$, we lose the coverage of these redder features, which negatively impacts the constraints we can measure on the SFHs. We therefore primarily focus on the evolution of the SFRs and stellar masses at $4.5 \le z \le 9$. We retain the galaxies in the $9 < z < 12$ sample in our analysis as this provides hints for this population, but we acknowledge that additional data (e.g., MIRI imaging) are needed to cover the same rest-frame optical features as the $z \le 9$ sample.

## 4. Results

Using the SFR and stellar-mass constraints, we now focus on the SFR–$M_*$ plane to study the SFMS for the galaxies at $4.5 < z < 12$. Figure 5 shows the SFR–$M_*$ plane for all galaxies in our sample. The panels show the relation using the SFRs averaged over the previous 10 Myr (left panel, SFR$_{10}$) and 100 Myr (right panel, SFR$_{100}$) as defined in Equations (3) and (4), respectively. One immediate conclusion from these figures is that the scatter in the SFR–$M_*$ (in log(SFR)) increases dramatically when the SFRs are averaged over the shorter 10 Myr timescale. However, there is little evolution in either the slope or normalization with redshift in the SFMS.

It is important to note that we make no selection in sSFR in our analysis of the SFR–$M_*$ relation. However, all galaxies in our sample would be considered star-forming when using SFR$_{100}$, as none would qualify as “quiescent” galaxies.<fnref n="32" /> The lowest sSFR$_{100}$ in the sample is $log(sSFR_{100} \text{ yr}^{-1}) \simeq -9.5$, above the canonical cutoff for quiescent galaxies at $z \le 6$. For example, S. Tacchella et al. (2022) advocate for a sSFR selection of quiescent galaxies that can be as high as $(sSFR \text{ yr}^{-1}) \simeq -10$. Similarly, other studies (see, e.g., V. Estrada-Carpenter et al. 2023, and references therein) have shown that the sSFR separation between star-forming and quiescent galaxies slowly increases to higher sSFRs at higher redshifts, but even extrapolating these trends to the highest redshifts in our sample implies that all the galaxies in our sample are “star-forming” when considering sSFR$_{100}$. We therefore make no distinction between galaxies that are star-forming or quiescent in our sample.

### 4.1. The SFMS in the Epoch of Reionization

We parameterize the SFMS as a power-law function between the SFR and stellar mass, $M_*$, where $log \ SFR \sim \alpha \ log \ M_*$. We divided our sample into five bins of redshift (listed in Table 2).

<footnote n="32">³² There is at least one spectroscopically confirmed, quiescent galaxy in our sample, recently published by A. de Graaff et al. (2024b). In our analysis, the *P* (z) based on the broadband photometry favors solutions with some amount of ongoing star formation, so we would not classify it as “quiescent.” However, when we fix the redshift to that reported by A. de Graaff et al. (2024b) we recover the quiescent solution. We therefore may be biased against such objects with ambiguous nature. These galaxies should be rare (K. Chworowsky et al. 2024), and quantifying their impact will be possible with future JWST spectroscopic programs in these fields (e.g., A. de Graaff et al. 2024a; M. Dickinson et al. 2024)</footnote>
<table><caption>**Table 2**
Results from the MCMC Analysis for Both SFR<sub>10</sub> (top) and SFR<sub>100</sub> (Bottom)</caption><thead><tr><th>Redshift</th><th>Sample Size</th><th>Stellar Mass 90% Completeness</th><th>α</th><th>Median log(M<sub>*</sub>)</th><th>β(z)</th><th>σ(z)</th></tr></thead><tbody><tr><td colspan="7">MCMC Results Measured With SFR<sub>100</sub></td></tr><tr><td>All galaxies</td><td>1863</td><td></td><td>0.79 ± 0.04</td><td>8.74</td><td></td><td></td></tr><tr><td>4.5 &lt; z &le; 5</td><td>475</td><td>8.03</td><td></td><td></td><td><b>0.51</b> ± 0.05</td><td><b>0.10</b> ± 0.03</td></tr><tr><td>5 &lt; z &le; 6</td><td>921</td><td>8.15</td><td></td><td></td><td><b>0.52</b> ± 0.02</td><td><b>0.09</b> ± 0.02</td></tr><tr><td>6 &lt; z &le; 7</td><td>291</td><td>8.51</td><td></td><td></td><td>0.61 ± 0.06</td><td><b>0.09</b> ± 0.05</td></tr><tr><td>7 &lt; z &le; 9</td><td>147</td><td>8.63</td><td></td><td></td><td><b>0.64</b> ± 0.07</td><td><b>0.10</b> ± 0.07</td></tr><tr><td>9 &lt; z &le; 12</td><td>29</td><td>8.63</td><td></td><td></td><td><b>0.89</b> ± 0.05</td><td><b>0.17</b> ± 0.15</td></tr><tr><td colspan="7">MCMC Results Measured With SFR<sub>10</sub></td></tr><tr><td>All galaxies</td><td>1863</td><td></td><td>0.80 ± 0.01</td><td>8.74</td><td></td><td></td></tr><tr><td>4.5 &lt; z &le; 5</td><td>475</td><td>8.03</td><td></td><td></td><td><b>0.40</b> ± 0.04</td><td><b>0.34</b> ± 0.03</td></tr><tr><td>5 &lt; z &le; 6</td><td>921</td><td>8.15</td><td></td><td></td><td><b>0.27</b> ± 0.02</td><td><b>0.35</b> ± 0.02</td></tr><tr><td>6 &lt; z &le; 7</td><td>291</td><td>8.51</td><td></td><td></td><td><b>0.44</b> ± 0.06</td><td><b>0.42</b> ± 0.05</td></tr><tr><td>7 &lt; z &le; 9</td><td>147</td><td>8.63</td><td></td><td></td><td><b>0.37</b> ± 0.07</td><td><b>0.44</b> ± 0.06</td></tr><tr><td>9 &lt; z &le; 12</td><td>29</td><td>8.63</td><td></td><td></td><td>0.23 ± 0.05</td><td><b>0.41</b> ± 0.08</td></tr></tbody></table>

Note. Fitted parameters for Equation (5) are shown in bold (β(z) and σ²(z)), and we list α and the median stellar mass (M<sub>*</sub>) for clarity.

Two of these bins lie at $z < 6$ (nominally after the end of the Epoch of Reionization, EoR; G. D. Becker et al. 2021), and two bins at $6 < z < 9$ (during reionization). We also include a bin with $9 < z < 12$, but we note that at these redshifts the measurements of the SFR<sub>10</sub> for galaxies are less well constrained as the rest-frame optical nebular emission lines shift out of the NIRCam wavelength coverage.

To account for mass completeness, we include only galaxies in each redshift bin that have stellar masses above an estimated 90% completeness limit. We derived the latter following the method described in D. Marchesini et al. (2009). We take all galaxies with $m_{F277W} \le 27.7$ mag, 0.5 mag above our selection limit. We dim these galaxies (in magnitude) and lower their stellar mass by a corresponding amount in small increments until we reach the magnitude ($m_{F277W}$) at which 90% of the galaxies are still brighter than our selection criteria at a given redshift. We then take the stellar mass corresponding to the brightest galaxy that would fail our magnitude limit ($m_{F277W} \le 28.2$ mag) as the mass-completeness limit. We do this for each redshift bin as the stellar-mass threshold evolves with redshift. We list the stellar-mass completeness limits in Table 2.

Finally, there is a potential bias from strong emission lines on the galaxy selection. This would occur for our sample at $4.5 < z < 5.0$, where the JWST/NIRCam F277W contains the strong, rest-optical $H\beta+[OIII]$ emission lines. This could impact our selection function, as galaxies nearest the mass limit could be biased toward galaxies with stronger emission-line EWs. We tested if this impact is significant by estimating the average impact of these lines on the stellar-mass completeness limit, and refitting the SFMS relations. C. Simmonds et al. (2024) report an average EW($H\beta+[OIII]$) $\simeq 600$ Å with an upper bound (maximum) of EW $\simeq 2000$ Å for photometrically selected samples from JWST at $2 < z < 9$. We estimated the impact of these emission lines on the inferred magnitudes for galaxies at $4.5 < z < 5.0$ following C. Papovich et al. (2001). This shifts the magnitudes $-0.45$ and $-1.0$ mag, respectively, which corresponds to a change in stellar mass of 0.18 dex and 0.43 dex, respectively. We then repeated our fit of the SFMS for the SFR<sub>10</sub> values with these higher stellar-mass limits,

where the sample is now more mass complete as it includes both galaxies with these high emission-line EWs and those without. However, we find that changing the mass limit has very little effect on the inferred parameters of the SFMS. For an average EW([OIII]) $\sim 600$ Å (2000 Å), we get $\alpha = 0.81 \pm 0.04$ ($0.81 \pm 0.07$), $\beta(4.5 < z < 5) = 0.38 \pm 0.03$ ($0.37 \pm 0.03$), and $\sigma_{int}(4.5 < z < 5) = 0.31 \pm 0.03$ ($0.33 \pm -0.04$) In both cases, the resulting slope ($\alpha$), normalization ($\beta(z)$), and intrinsic scatter ($\sigma(z)$) change by less than the $1\sigma$ uncertainties listed in Table 2. We therefore conclude that galaxies with strong emission-line EWs do not overly bias the results, even at $4.5 < z < 5$, where the effect is potentially the most severe.

To model the SFMS, we define a power-law relation between the SFR and stellar mass as

$$ \log(\Psi_{*}) = \alpha \log\left(\frac{M_{*}}{10^{8.74} M_{\odot}}\right) + \beta(z) + \epsilon(z), \quad (5) $$

where $\Psi_{*}$ is the SFR in solar masses per year, and $\alpha$ is the slope of the SFMS. There are three parameters, $\alpha$, $\beta(z)$, and $\epsilon(z)$. Here, we assume a single value for the slope, $\alpha$, that does not evolve with redshift over $4.5 < z < 12$. This is an assumption supported by the fact that we observe no evidence in a change in slope from Figure 5, particularly in the SFR<sub>100</sub>-M<sub>*</sub> relation. However, this may be a result of the limited dynamic range we have in stellar mass in our sample.33 Larger-area and deeper data sets will be needed to provide better constraints on the evolution of this slope than what we can achieve here. We then allow for evolution in the normalization of the relation, $\beta(z)$, here defined at a mass of $\log(M_{*}/M_{\odot}) = 8.74$, the median stellar mass of our sample.

The variable $\epsilon(z)$ allows for the intrinsic scatter along the relation. We assume $\epsilon(z)$ is drawn from a normal distribution with mean, $\mu = 0$, and variance $= \sigma_{int}^2$; that is,

<footnote n="33">33 We verified this statement by fitting Equation (5) where the slope is allowed to evolve with redshift, $\hat{\alpha}(z)$. In this case, we find $\hat{\alpha}$ scatters around the value when assuming no evolution in redshift, $\alpha$, with $\hat{\alpha}(z) = 0.75 - 0.85$, with a typical uncertainty of 0.1 in each redshift bin from $4.5 < z < 9$. As we are most interested in the scatter about the SFMS, we choose to fit $\alpha$ with no redshift dependence as this provides tighter constraints on $\sigma$.</footnote>
$\epsilon(z) \sim N(0, \sigma_{\text{int}}(z)^2)$. The value of $\sigma_{\text{int}}(z)$ then represents the intrinsic scatter in the SFMS at each redshift (assuming here no dependence on stellar mass, but we consider such a dependence below).

We use a Monte Carlo sampler to determine the best fit and uncertainties. We use the Markov Chain Monte Carlo (MCMC) routine implemented with the Python package EMCEE (D. Foreman-Mackey et al. 2013) for this purpose. In the MCMC fit, for every walker in each iteration we fit a single $\alpha$ for all galaxies (in all redshift bins), where there is a separate $\beta$ ($z$) in each redshift bin. In each iteration of the fit, for each galaxy we resample from the covariance between the variables in the fit (i.e., $P(M_*, SFR)$). That is, we draw a new SFR and stellar mass from the two-dimensional posterior space for each galaxy, and in this way galaxies are allowed to move within their respective posterior spaces in each iteration of the MCMC fit. This is important as the SFR and stellar mass are strongly correlated (see gray ellipses in Figures 5–6).

The details of the likelihood used in the fit are provided in Appendix A. In summary, we split the likelihood into two parts. First, we model the relation

$$ \log(\Psi_*) = \alpha \log\left(\frac{M_*}{10^{8.72} M_{\odot}}\right) + \beta(z), \quad (6) $$

which becomes the first part of the likelihood function. Second, we subtract this relation from the observed data to model the scatter about the relation, $\epsilon(z)$. We then fit this in each redshift bin, taking the absolute value of the scatter, and use a $\chi^2$-squared statistic to determine the likelihood, assuming Gaussian scatter about the SFMS, to compare to the absolute value of a normally distributed sample with the walker's associated variance. The $\chi^2$-squared likelihood from this step is the second part of our likelihood function. We provide tests and further descriptions of this methodology in Appendix A and show that we are able to recover accurately the slope, normalization, and scatter in the relation.

We fit the SFMS separately for the SFR$_{10}$ and SFR$_{100}$ values. Figures 6 and 7 show the results for the SFMS from the MCMC fits. Table 2 lists the measured values for each parameter. We discuss the results in the next subsection.

## 4.2. The Slope and Normalization of the SFMS Relation

Figure 6 shows the results of the SFMS fits with SFR$_{100}$. The first five panels (left to right, top to bottom) of Figure 6 show the results for each of the five redshift bins ($4.5 < z < 12$), as labeled. The bottom-right panel shows the results for $\sigma_{\text{int}}(z)$, or the intrinsic scatter along the SFMS.

The SFMS we derive using the SFR$_{100}$ values is consistent with those from previous measurements at these redshifts (all measured using SFRs sensitive to the SFH over 100 Myr timescales; J. S. Speagle et al. 2014; B. Salmon et al. 2015; E. Curtis-Lake et al. 2021). The J. S. Speagle et al. (2014) relation is measured at lower redshifts. However, if we extrapolate this to the redshifts and values for our sample, the slope and normalization of the SFMS we derive agree. Similarly, our results also agree with those from B. Salmon et al. (2015), although they find evidence that the slope of the SFMS flattens with increasing redshift. Here, we assume a constant slope, where the evolution is then in the

normalization. The combination of these yields a similar SFMS relation (see Figure 6).

Figure 7 shows the SFMS fits using SFRs over 10 Myr timescales. The results from these fits are listed in the bottom rows of Table 2. The bottom-right panel of the figure again shows the intrinsic scatter, $\sigma_{\text{int}}(z)$, using SFR$_{10}$ (color-coded points and error bars), and the faded points in this panel show the $\sigma_{\text{int}}(z)$ values derived using the SFR$_{100}$ shown in Figure 6. Using SFR$_{10}$, we measure a slope for the SFMS of $\alpha = 0.80 \pm 0.01$, consistent with our measurement of the SFMS when using SFR$_{100}$. However, the normalization values for $\beta(z)$ are lower when using SFR$_{10}$ compared to SFR$_{100}$. This implies that on average the SFR$_{10}$ values are lower than SFR$_{100}$ (at fixed stellar mass). This is consistent with our measurement of the duty cycle for star formation (discussed below in Section 5.2).

We again find that the SFMS we derive using SFR$_{10}$ is consistent with values from the literature. The most relevant comparison is to the work of R. Smit et al. (2016), who used broadband photometry from Spitzer/IRAC to estimate the $H\alpha$ fluxes (and thus SFRs) for a spectroscopic sample of massive star-forming galaxies at $\langle z \rangle = 4.25$. As the $H\alpha$ flux traces the most massive stars it also will respond to variations on 10 Myr timescales, and thus similar to our SFR$_{10}$ values (R. C. Kennicutt 1998). R. Smit et al. (2016) derive a slope $\alpha = 0.78 \pm 0.23$, well matched to ours from the SFMS using SFR$_{10}$, $\alpha = 0.79 \pm 0.04$. Similarly, they derive a normalization for galaxies at $\log M_*/M_{\odot} = 10$ of $\beta = \log \text{SFR}/(M_{\odot} \text{ yr}^{-1} = 1.51 \pm 0.07)$. If we extrapolate our normalization to this mass we obtain $\beta(z = 4.5) = 1.40 \pm 0.04$, consistent at the $1\sigma$ level. This suggests some consensus in the shape of the SFMS at $z \sim 4-5$ when using galaxy SFRs averaged over 10 Myr. However, the scatter in the relation has yet to reach consensus (which we discuss in the next subsection).

The evolution in the normalization of the SFMS depends on the choice of SFR averaged over the different timescales. Using SFR$_{100}$ we measure a slow increase in normalization with increasing redshift, going from $\beta \simeq 0.5$ at $z = 4-6$ to $\beta \simeq 0.6$ at $z = 6-9$. This agrees with previous studies at $1 \le z \le 6$ (e.g., K. E. Whitaker et al. 2012; J. S. Speagle et al. 2014; B. Salmon et al. 2015; P. Popesso et al. 2023). However, this trend is markedly different for the normalization of the SFMS using the SFR$_{10}$ values. These are consistently lower than the normalization values derived using SFR$_{100}$, and show no evidence of evolution from $z = 4.5$ to 9. We do not consider the trend to higher redshifts, $z > 9$, to be significant as at these redshifts the strong optical emission lines tracing the SFR on 10 Myr timescales shift beyond the coverage of the JWST/NIRCam data. This has implications for the duty cycle (because on average the SFR$_{10}$ values are smaller than SFR$_{100}$), and we consider this in Section 5.2 below.

We see possible evidence that the normalization of the SFMS using SFR$_{10}$ is lower for galaxies in the $5 < z < 6$ redshift bin. This is notable as there is a purported galaxy overdensity in CEERS in this bin, at $z \simeq 5.1$ (K. Chworowsky et al. 2023), and it may be that this modifies the normalization of the SFMS relation through environmental processes, similar to other results at lower redshifts, $z \sim 0$ to $z \sim 1$ (L. J. Old et al. 2020; A. Boselli et al. 2023). It will be important to compare the results here with similar analyses in other fields to test if effects of the large-scale structure are biasing our conclusions.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** The SFR-$M_*$ relation for the full sample of galaxies using the SFRs averaged over the past 100 Myr (SFR$_{100}$). Each panel show the relation for a different redshift bin, as labeled. The data points show median values for each galaxy where the error bars show the inner-68th percentile of the posteriors on SFR$_{100}$ and $M_*$. Dark-filled data points lie above the 90% mass completeness limit (faint gray symbols lie below this limit and are excluded from the fit). The solid colored lines show the median best-fit relation from the MCMC fit (Equation 6), where the shaded region shows the $1\sigma$ error on this relation. For comparison, we show measurements from the literature at $z \lesssim 7$ (E. Curtis-Lake et al. 2021; B. Salmon et al. 2015) and we extrapolate the relation from J. S. Speagle et al. (2014) assuming that their relation continues its evolution out to $z = 12$. The solid green line shows the measurements at $z \gtrsim 6$ (L. Ciesla et al. 2024). The bottom-rightmost panel shows the inferred intrinsic scatter, $\sigma_{\text{int}}$, as a function of redshift with asymmetric errors estimated from the MCMC posterior.</figcaption></figure>

### 4.3. Scatter in the SFMS Relation

As illustrated in Figure 7 (lower-right panel), when using SFR$_{10}$ there is an apparent increase in the intrinsic scatter in the SFMS, $\sigma_{\text{int}}(z)$, with increasing redshift from $z \sim 4.5$ to $z \sim 9$. Fitting to the full samples (over all stellar masses), we see the intrinsic scatter increases from $\sigma_{\text{int}} = 0.34 \pm 0.03$ dex at $4.5 < z < 5$ to $\sigma_{\text{int}} = 0.44 \pm 0.06$ dex at $7 < z < 9$. In contrast, there is no evidence for evolution in $\sigma_{\text{int}}$ when using the SFRs averaged over the longer 100 Myr timescales, SFR$_{100}$. The

interpretation is that there is an increase in the variability on the shorter timescale of SFR$_{10}$ for galaxies at higher redshifts. However, these values for $\sigma_{\text{int}}$ are measured for the entire galaxy population at each redshift.

We also measured the evolution in the intrinsic scatter along the SFMS as a function of stellar mass in each redshift bin. We perform this measurement using both SFR$_{10}$ ($\sigma_{10}$) and SFR$_{100}$ ($\sigma_{100}$). For each stellar mass and redshift bin, we recalculated the scatter as half the distance from the 16th to 84th percentiles (i.e., the inner-68th percentile) of the distribution in SFR. We
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Same as Figure 6, but the SFR averaged over the previous 10 Myr (SFR<sub>10</sub>). In the bottom-rightmost panel, we show the results for σ<sub>int</sub> for SFR<sub>10</sub> as darker points. For comparison, the lighter points show the σ<sub>int</sub> measurements from SFR<sub>100</sub> (from Figure 6). The scatter derived from the SFR<sub>10</sub> measurements are substantially higher than those derived from SFR<sub>100</sub>.</figcaption></figure>

again assume that SFRs are distributed normally, and therefore
half of the inner-68th percentile is defined as σ and we can
estimate the uncertainty on the observed scatter as

$$
\delta\sigma = \frac{\sigma}{\sqrt{2(N-1)}}, \qquad (7)
$$

where *N* is the number of galaxies in the estimation of *σ*. We
show these measurements of the scatter and their uncertainties
in Figure 8. We include as a comparison measurements from
the literature (I. Shivaei et al. 2015) and predictions from
models and simulations (S. M. Wilkins et al. 2024;
L. Y. A. Yung et al. 2024).

There are several interesting trends in Figure 8. First, σ<sub>100</sub> is
roughly constant in both stellar mass and redshift over the

entire range of our sample, with σ<sub>100</sub> ≈ 0.2–0.3 dex. This is
consistent with previous results from lower redshifts, within the
quoted uncertainties. I. Shivaei et al. (2015) find that the scatter
in SFR derived for galaxies at z ~ 2 from rest-FUV data is
σ<sub>FUV</sub> ≃ 0.3 dex. C. Schreiber et al. (2015) obtain similar
results, finding σ ≃ 0.3 dex using far-IR data from Herschel for
massive star-forming galaxies at 2 < z < 4 with no apparent
evolution in redshift nor stellar mass. As these observables
probe star formation on ~100 Myr timescales, this is consistent
with our observed scatter in SFR<sub>100</sub>. L. Ciesla et al. (2024)
recently used JWST/NIRCam data to measure the scatter in the
UV absolute magnitudes of galaxies at z > 6, finding
σ<sub>UV</sub> ≃ 1.2 mag. Converting this to a scatter on the SFR yields
σ<sub>UV</sub> ≃ 0.5 dex, but this includes both the scatter in the SFR and
the dispersion in dust attenuation. As the latter is typically
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Observed scatter in the SFMS as a function of stellar mass and redshift. Each panel shows a different redshift, where the colors denote different redshift bins as in Figure 7 and as labeled. Square markers connected by solid lines show the observed $1\sigma$ scatter for SFR$_{10}$ (the SFR averaged over 10 Myr), and triangles connected by dotted lines show the observed $1\sigma$ scatter in SFR$_{100}$ (the SFR averaged over 100 Myr). Error bars on $\sigma_{\text{log SFR}}$ are estimated using Equation (7). In the panel for each redshift we show theoretical predictions for $\sigma_{\text{log(SFR)}}$ from the SC-SAM (black; L. Y. A. Yung et al. 2024) and FLARES (green; C. C. Lovell et al. 2021; S. M. Wilkins et al. 2024), where solid (dotted) lines show these predictions for SFR$_{10}$ (SFR$_{100}$) and the shaded regions indicate the 16th-to-84th percentiles. In the top-left panel (and bottom-rightmost panel), the “×” and “+” show the measured observed scatter along the SFMS at $z=2$ from I. Shivaei et al. (2015) for FUV SFRs and H$\alpha$ SFRs (as labeled). The bottom-rightmost panel compares the results for all redshift and stellar mass bins (error bars are removed, and points at the highest stellar masses are removed for clarity).</figcaption></figure>

$\sigma(A_{\text{FUV}}) = 0.5–0.6$ mag (J. Álvarez-Márquez et al. 2016), this may account for at least half of the scatter in $\sigma_{\text{UV}}$. We therefore consider these measurements consistent with our value for $\sigma_{100}$ here.

Second, at all redshifts $\sigma_{10}$ is higher than $\sigma_{100}$, with typical values of $\sigma_{10} \simeq 0.4–0.6$ dex. This is illustrated in the bottom-right panel of Figure 8. At lower redshifts, $4.5 < z < 6$, $\sigma_{10}$ shows no dependence on stellar mass, with a roughly constant value of $\sigma_{10} \simeq 0.5$ dex. At higher redshifts, $6 < z < 9$, there is a positive correlation between stellar mass and $\sigma_{10}$: the scatter increases with stellar mass, from $\sigma_{10} \simeq 0.4$ dex at $\log M_*/M_{\odot} \simeq 8.75$ to $\simeq 0.6$ dex at $\log M_*/M_{\odot} = 10$.

Comparing to measurements of $\sigma_{\text{SFR}}$ using H$\alpha$ data at $z \sim 2$ from I. Shivaei et al. (2015), the observed scatter in SFR$_{10}$ we obtain is higher by a factor of $\sim 2$. This may imply that there is an increase in the scatter when averaged over $\sim 10$ Myr timescales from $z \sim 2$ to $z > 4$ such that SFHs contain higher variability at higher redshifts. It is also possible that studies which require detection of H$\alpha$ could lead to a bias toward lower scatter (a form of Malmquist bias). I. Shivaei et al. (2015) included only galaxies with H$\alpha$ detections, where our sample includes all galaxies, including objects that lack evidence for strong nebular emission in their photometry (see Figure 4). The presence of these lower-SFR galaxies increases the scatter in
our measurement. Confirming the evolution in SFR<sub>10</sub> could be tested with spectroscopy around the regions of the strong rest-optical emission lines (e.g., Balmer lines, [O III], etc.) from future JWST spectroscopy for galaxies at z > 4.5.

At the higher redshifts, $6 < z < 9$, the scatter in SFR<sub>10</sub> appears to increase with increasing stellar mass ($\log M_*/M_\odot \gtrsim 9.5$). However, we consider this may be a systematic or selection effect. It is arguable that our samples are less complete in SFR at these stellar masses in the highest-redshift bins. We have estimated this effect to ensure the data are complete in stellar mass at 90% (see Section 4.1). However, this estimate assumes a slowly changing population in both luminosity and stellar mass. If there is substantially more scatter in the SFR at lower stellar masses, then the SFR distribution may be truncated and biased away from galaxies with lower SFRs. This would make our estimate incomplete and artificially tighten the distribution (another instance of Malmquist bias). This can be tested with deeper imaging data from JWST. However, the stellar-mass completeness will not impact the higher-mass end of the SFMS, $\log M_*/M_\odot \gtrsim 9.5$, which is where we observe a significant increase in the observed scatter from $4.5 < z < 5$ to $7 < z < 9$. We can identify no systematic effect that would impact this trend, which we therefore consider to be robust.

In Figure 8, we also compare to predictions for the scatter in the SFMS from the Santa Cruz semi-analytic model (SC-SAM; R. S. Somerville et al. 2015; L. Y. A. Yung et al. 2019) and First Light And Reionisation Epoch Simulations (FLARES; C. C. Lovell et al. 2021). We will return discuss the comparison between our results and these predictions in Section 5.3.

One potential source of systematic that we considered was the fact that at $z > 7.6$, $H\alpha$ shifts beyond the NIRCam coverage. The data remain sensitive to the $H\beta+[O III]$ emission for $9 < z < 12$. However, as $H\beta+[O III]$ lie in the reddest NIRCam band, there is a potential for increased uncertainties in the modeling. This fact leads to larger uncertainties on quantities like $\sigma(\log SFR_{10})$. This is apparent in Figure 7 (bottom-right panel) and Figure 8, which shows that the uncertainty on the scatter of the SFR<sub>10</sub> values increases from $6 < z < 7$ (where the data include $H\alpha$ emission) and $7 < z < 9$ (where the data do not). To improve the understanding of the accuracy of the uncertainties on the scatter of the SFMS likely will require improved measurements of the contribution of $H\beta+[O III]$ to the NIRCam broadband flux densities. This will be possible for large samples from forthcoming JWST spectroscopic programs (e.g., M. Dickinson et al. 2024).

### 4.4. The Distribution of Specific SFR

Given the apparent differences in the observed scatter in galaxy SFRs with stellar mass, we test whether these trends are also apparent in the distributions of sSFR. Figure 9 shows the sSFR<sub>10</sub> (the specific SFR using the SFRs averaged over 10 Myr timescales) as a funciton of redshift in two bins of stellar mass, $\log M_*/M_\odot \in [8.7, 9.3]$ and $\in [9.3, 10.3]$. To study the trends, we fit a parabolic function to the data in each bin of stellar mass allowing for evolution in redshift, z, following the form

$$ \log(\text{sSFR}(z)) = A \times z^2 + B. \qquad (8) $$

We measure the parameters using a simple least-squares regression. The fitted models for each mass bin are shown in both panels of Figure 9 for comparison, where we show the

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Specific SFRs (sSFRs) measured from SFR<sub>10</sub> as a function of redshift for the galaxies in our sample in bins of stellar mass. The top panel shows lower-mass galaxies with $10^{8.7} \le M_*/M_\odot < 10^{9.3}$. The bottom panel shows higher-mass galaxies with $10^{9.3} \le M_*/M_\odot < 10^{10.3}$. The gray circles in both panels show the mean $\log(\text{sSFR}_{10})$ for individual galaxies in each sample, where the error bars represent the inner-68th percentiles (these includes the contribution from the covariance between stellar mass and SFR; see Section 4.1). Vertical dotted lines show even-lookback time steps of 200 Myr to demonstrate the trends in both redshift and cosmic time. The blue and red lines and shaded areas show the best-fit and $1\sigma$ range from parabolic fits to the data top (lower-mass galaxies) and bottom panels (higher-mass galaxies), respectively. The curves for each bin of stellar mass are shown in both panels to contrast the results (but are dimmed for clarity), as labeled. At all redshifts the higher-mass galaxies have lower SFR<sub>10</sub>, where the offset is consistently 0.3 dex. This is significant at the $> 1.5\sigma$ level at all the redshifts of our sample.</figcaption></figure>

median best fits as the solid lines and $1\sigma$ range as the lighter shaded regions. At all redshifts the median sSFR<sub>10</sub> for higher-mass galaxies is lower by $\Delta\log(\text{sSFR}_{10}) \approx 0.3$ dex than those for lower-mass galaxies. This difference is significant at the $> 1.5\sigma$ level at all redshifts of the galaxies in our sample.

In addition, the scatter in the distribution of sSFR<sub>10</sub> is larger for the higher-mass galaxies. For the higher-mass galaxies, we measure $\sigma(\log \text{sSFR}_{10}) = 0.35$ dex at $z \sim 5$ and $0.45$ dex at $z \sim 8$. For the lower-mass galaxies, we find $\sigma(\log \text{sSFR}_{10}) = 0.36$ and $0.35$ at $z \sim 5$ and $8$, respectively. In all cases this is consistent with our observation above (Figure 8) that the scatter in the SFR<sub>10</sub> is increasing with increasing mass and redshift. Here, we find this is also the case when using the sSFRs.

We explore this further in Figure 10, where we fit KDEs to the sSFR<sub>10</sub> distributions in bins of redshift for the lower-mass (left column, solid black lines) and higher-mass galaxies (right column, solid black lines). We do not include galaxies at $9 < z < 12$ for the reason above that we lose the ability to accurately constrain SFR<sub>10</sub> as the rest-frame optical nebular emission lines shift beyond the wavelength coverage of NIRCam.

The sSFR<sub>10</sub> distributions for the lower-mass galaxies in each bin of redshift appear consistent with a unimodal distribution. We therefore fit the sSFR<sub>10</sub> distributions for the lower-mass galaxies in each redshift bin with a single-Gaussian model (illustrated in each panel) to the distributions for the lower-mass galaxies in each redshift bin. In each case these are consistent with the KDEs.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Specific SFR<sub>10</sub> (sSFR<sub>10</sub>) distributions for our galaxy sample in bins of mass and redshift. In each panel of the left column we show a KDE for the lower-mass galaxies with log(M<sub>*</sub>/M<sub>☉</sub>) < 9.3 as the black solid line. A best-fit Gaussian model fit to the distribution is shown as the blue dashed line. In each panel of the right column we show a KDE for the higher-mass galaxies, log(M<sub>*</sub>/M<sub>☉</sub>) ∈ [9.3, 10.3], as the black solid line. A best-fit double-Gaussian model is shown as the red dashed line. Each row shows the results for different bins of redshift (as labeled). The information inset in the panels gives the number of galaxies, *N*, along with the mean of the Gaussian fit(s), μ (μ<sub>1</sub> and μ<sub>2</sub> for the double-Gaussian fits), and Δμ = μ<sub>1</sub> - μ<sub>2</sub>, which shows the difference between the two means from the double-Gaussian fits (in units of dex). The short vertical lines show the means of the Gaussian fits, red lines show the fits to the higher-mass samples, and blue lines show the fits to the lower-mass samples (the latter is shown in both panels of each row for comparison). There is evidence for a bimodal distribution of galaxy sSFR<sub>10</sub> for the higher-mass sample of galaxies at all redshifts in our sample. Comparing the location of the means from the bimodal distribution for the higher-mass sample of galaxies with the peak from the unimodal distribution of the lower-mass sample of galaxies shows this appears to be the result of an increase in the number of galaxies with lower sSFR<sub>10</sub>.</figcaption></figure>

The sSFR<sub>10</sub> distributions for the higher-mass galaxies are more extended, and appear inconsistent with a unimodal distribution. We fit a single-Gaussian model to the each of distributions, and a Kolmogorov–Smirnov (KS) test fails to reject the hypothesis that these are drawn from the same parent distribution, with an average of the *p*-values of each redshift bin of *p* = 0.60. However, the distributions of sSFR<sub>10</sub> for the higher-mass galaxies appear bimodal. We therefore fit a double-Gaussian model to each of the sSFR<sub>10</sub> distributions

for the higher-mass galaxies, illustrated in each panel of Figure 10. In each panel of the right column, we list the means of the double-Gaussian models, and we list the difference between the means Δμ = μ₁ − μ₂. Here, we find that the KS test still does not reject the hypothesis that the distributions are drawn from the same parent distribution, with an average of the *p*-values of each redshift bin of *p* = 0.73, a ~22% increase in confidence over the single-Gaussian model.

The fact that the distributions of sSFR<sub>10</sub> for the higher-mass galaxies appear bimodal means they need an additional component, and this is not the case for the lower-mass galaxies. We can study this by comparing the location of the means (“peaks”) in the distributions of sSFR<sub>10</sub> for the lower-mass and higher-mass samples as a function of redshift. In each panel of Figure 10, we indicate the peak of the single component for the lower-mass sample (left columns, blue vertical lines). We also indicate the peaks of both the components for the higher-mass sample (right columns, red vertical lines). Comparing the peak values, we see that the peak of the sSFR<sub>10</sub> distribution of the lower-mass galaxies is always comparable, i.e., within ≲0.3 dex, to the peak of the higher-sSFR Gaussian for the higher-mass galaxies. This means that the higher-mass galaxies do not show evidence for higher sSFR<sub>10</sub> compared to the lower-mass galaxies. Rather, the higher-mass galaxies show a second component (modeled by a Gaussian here) with a peak at lower sSFR<sub>10</sub>, shifted by 0.6–1 dex compared to the high-sSFR<sub>10</sub> peak. This means that the higher-mass galaxies show evidence for a population of galaxies with lower sSFR<sub>10</sub>. This appears to be driving the increase in the total scatter in SFR<sub>10</sub> for galaxies at higher mass. We therefore conclude that the increase in scatter is a result of there being more galaxies with lower sSFR<sub>10</sub> (e.g., “napping” galaxies) rather than an increase in the relative number of galaxies with higher sSFR<sub>10</sub> (e.g., there is not an excess of “bursting” galaxies). We discuss the implications of this below.

## 5. Discussion

The results presented here favor a scenario where the mode of star formation in galaxies in the early Universe is more variable than that observed in galaxies at later cosmic times. The variable nature of star formation in early galaxies has been predicted in some simulations (e.g., G. Sun et al. 2023; S. M. Wilkins et al. 2024; L. Y. A. Yung et al. 2024), and has been proposed as one of the solutions to account for the higher number density of high-redshift galaxies at the bright end of the UV luminosity function (UVLF; e.g., C. A. Mason et al. 2023; S. L. Finkelstein et al. 2022; K. Chworowsky et al. 2023). Other studies based on SED modeling and measurements of galaxy EW distributions also favor “burstier” star formation in high-redshift galaxies (R. Endsley et al. 2022, 2024; L. Ciesla et al. 2024). This runs counter to studies at relatively lower redshifts that favor more slowly evolving SFHs (e.g., C. Papovich et al. 2011; N. A. Reddy et al. 2012). The fact that we find galaxy SFRs have more scatter on short timescales, ≲100 Myr, implies that the mechanism driving this increased variability also acts on these shorter timescales. However, “burstier” may not be an accurate description for the SFHs of these galaxies. The galaxies do not have excessively high sSFR<sub>10</sub> values, but rather show an increase in the frequency of galaxies with low sSFR<sub>10</sub> values. Therefore, the higher variability in galaxies in the EoR is driven by one or more
physical mechanisms that interrupt otherwise “normal” star
formation. Here, we discuss the implications on galaxy growth
at $z > 6$ (in the EoR) as a result of the increased variability in
the galaxy SFRs. We also discuss how our results compare to
theoretical predictions and to the findings from other recent
studies of galaxies at these epochs.

## 5.1. An Increase in the Scatter in the SFRs of Galaxies at the Epoch of Reionization

We see an increase in intrinsic scatter ($\sigma_{\text{int}}$) along the SFMS measured with SFRs on 10 Myr timescales at increasing redshift, and we see that the scatter in SFR$_{10}$ ($\sigma_{10}$) also increases with increasing stellar mass. This result is similar to the findings from many observational studies of galaxies at these epochs. One finding is that star formation must become increasingly more variable as we look to higher redshifts, based on the strength of galaxy nebular emission lines. R. Endsley et al. (2024) find that emission-line EW distributions are broader, ranging to larger values, for galaxies at $6 < z < 9$ in the JWST Advanced Deep Extragalactic Survey (JADES) field (M. J. Rieke et al. 2023; C. C. Williams et al. 2024). An explanation for this is that the galaxies contain a higher relative number of OB stars. R. Endsley et al. (2024) also find that a considerable subset of the SEDs in the CEERS sample show strong evidence for Balmer/4000 Å breaks, which indicates there is a substantial fraction of galaxies observed in CEERS which have SEDs dominated by older ($\gtrsim 100$ Myr) stellar populations. R. Endsley et al. (2024) concluded that these observations favor an increase in the variability of star formation at this epoch. This agrees with our findings. Moreover, R. Endsley et al. (2024) provide an independent analysis as they fit the galaxy SEDs using BEAGLE with a constant SFH (J. Chevallard & S. Charlot 2016). Therefore, this is further support for an increase in the scatter of SFRs of galaxies in the EoR.

Studies of the UVLF and its evolution at $z \gtrsim 7$ have invoked a high variability in star formation to explain the excess number density of galaxies at the bright end of the UVLF (e.g., S. L. Finkelstein et al. 2023; C. A. Mason et al. 2023). X. Shen et al. (2023) conclude that an increase in the scatter of UV absolute magnitudes would account for the differences of the UVLF, where they require $\sigma_{M_{UV}} \approx 0.75$ mag at the bright end and $\sigma_{M_{UV}} \approx 1.5$ mag at the faint end of the UVLF. In the JADES fields, L. Ciesla et al. (2024) measure $\sigma_{M_{UV}} \approx 1.2$ mag at $z > 9$ for galaxies with $-16.6 > M_{UV} > -18.4$, consistent with this idea. Comparing to our work here, we calculate a scatter in the SFR for galaxies at $7 < z < 9$ of $\sigma_{10} \approx 0.6$ dex for high-mass galaxies. Assuming $\sigma(\log \text{SFR}) \sim 0.4 \cdot \sigma(M_{UV})$, this corresponds to $\sigma_{M_{UV}} \sim 1.5$ mag. Therefore, the scatter we derive appears sufficient to account for the higher number density of objects at the bright end of the UVLF. We defer a detailed comparison of how the increase in the scatter of the SFRs impacts the UVLF and the stellar mass function to a future study.

At the high-mass end of the SFMS, we observe a bimodality
in the distribution of sSFR₁₀ (Figure 10), implying that galaxies
with log(M*/M☉) ≳ 9.3 are driving the intrinsic scatter along
the SFMS on 10 Myr timescales. Furthermore, our results show
this is apparently a result of an increase in the frequency of
galaxies with low sSFRs. There is growing evidence from other
studies that support this conclusion. L. Ciesla et al. (2024)
presented evidence for two populations of star-forming

galaxies at $7 < z < 9$ based on the SFR “gradient.” One
population has a SFR gradient moving them along the SFMS,
while the other population has a gradient moving them below
the SFMS (toward lower sSFRs). Qualitatively, the first
population is consistent with the galaxies we show in
Figure 3 that have SFRs which place them on the SFMS.
The second population includes the kinds of galaxies we show
in Figure 4, where we are observing them during a recent
increase in star formation, or in a declining phase (so-called
“napping” galaxies, which have experienced a cessation in star
formation). The latter population includes galaxies similar to
those seen in recent studies. For example, T. J. Looser et al.
(2024) report on a galaxy, JADES-GS-z7-01-QU, at $z = 7.3$
with no nebular emission lines, which they conclude ceased
forming stars 10–20 Myr earlier after a short, substantive
episode of star formation. T. J. Looser et al. (2023) present
additional evidence for galaxies at $5 < z < 11$ in this “napping”
stage (what they call “mini-quenched” galaxies). V. Gelli et al.
(2023) argue that ~30% of galaxies with stellar masses
log M*/M☉ ≲ 9.5 at these redshifts have SEDs indicating they
are in this “napping” mode, experiencing lower sSFRs. There is
therefore growing evidence that galaxies in the EoR experience
active periods of star formation, followed by periods of
cessation.

## 5.2. The Duty Cycle of Star Formation in the Epoch of Reionization

The scatter in the SFMS implies that galaxies spend some
period of time in starburst phases *and* in “napping” phases. The
fraction of the galaxy population in this phase places an
estimate on this time. We then define a “duty cycle” of star
formation as the relative amount of time galaxies spent in the
star-forming phase (when they lie on the SFMS).

We estimate the fraction of galaxies in active star-forming
stages using the ratio of the SFRs averaged over 10 Myr
compared to the SFRs averaged over 100 Myr. We define this
ratio as $R_{SFR} \equiv \text{SFR}_{10}/\text{SFR}_{100}$, similar to the quantity studied
by T. J. Looser et al. (2023). We then use this ratio to
quantitatively identify galaxies in a declining star-forming
stage, with $\log R_{SFR} < 0$, and galaxies that are in active star-
forming stages, with $\log R_{SFR} \ge 0$. Figure 11 shows the ratios
of $\log R_{SFR}$ as a function of stellar mass for galaxies in the
different redshift bins. Figure 11 illustrates that lower-mass
galaxies include galaxies with high and low values of $\log R_{SFR}$
in nearly equal amounts, with a distribution centered near unity.
In contrast, higher-mass galaxies appear to have a higher
fraction of objects with low ratios, $\log R_{SFR} < 0$. Quantifying
this trend, at redshifts $4.5 < z < 9$, >60% of higher-mass
galaxies with $\log M*/M☉ > 9.3$ have low ratios, $\log R_{SFR} < 0$.
That is, at these redshifts and stellar masses, galaxies are more
likely to be in a phase of lower SFR than in an active star-
forming phase. Generally, this implies that *on average* lower-
mass galaxies have SFHs that are constant on 10–100 Myr
timescales, while higher-mass galaxies have SFHs that have
declined over these timescales.

Figure 12 illustrates these trends by showing the fraction of galaxies with high $R_{SFR}$, that is, $f(\log R_{SFR} > 0)$ (top panel), in bins of stellar mass and redshift. Lower-mass galaxies with $\log(M*/M_{\odot}) < 9.3$ have near equal fractions above and below the unity line, with $f \approx 0.5$. It is interesting to see if these trends continue or change at lower masses, where there is some evidence that the fraction of galaxies with $\log R_{SFR} > 0$
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Ratio of the SFRs averaged over 10 Myr (SFR<sub>10</sub>) to the SFRs averaged over 100 Myr (SFR<sub>100</sub>) as a function of stellar mass and redshift, log R<sub>SFR</sub> = logSFR<sub>10</sub> - logSFR<sub>100</sub>. Each panel shows galaxies in a different redshift bin (as labeled; colored as in Figures 7 and 8). We divided the sample into a lower-mass sample (log(M<sub>*</sub>/M<sub>☉</sub>) ∈ [8.7, 9.3]) and a higher-mass sample (log(M<sub>*</sub>/M<sub>☉</sub>) ∈ [9.3, 10.3]), delineated by the dotted vertical lines. The thick horizontal solid line shows where the SFRs are equal, SFR<sub>10</sub> = SFR<sub>100</sub>. The gray-shaded region in each panel shows where our sample drops below 90% completeness in stellar mass.</figcaption></figure>

increases, indicating a higher fraction of bursts (e.g., T. J. Looser et al. 2023). Higher-mass galaxies with log(M*/M☉) > 9.3 have a lower fraction of galaxies with high R<sub>SFR</sub> values with f(log R<sub>SFR</sub> ≥ 0) of <40%. This fraction decreases with increasing redshift, falling to ≈ 15%–30% at 7 < z < 9.

We take the fraction of galaxies above the unity relation, log R<sub>SFR</sub> > 0, which yields an estimate of the duty cycle in a redshift and mass range. We further can estimate the average amount of time galaxies spend in the active star-forming stage by multiplying the duty cycle with the time period (Δt) spanned by the redshift range in each bin:

$$t(\log R > 0) = f(\log R > 0) \times \Delta t. \qquad (9)$$

We list the duty cycles and the time periods of active star formation for lower-mass and higher-mass galaxies in bins of redshift over the range 4.5 < z < 9 in Table 3.

Regarding the duty-cycle measurements, we see several interesting trends. First, higher-mass galaxies spend a shorter period in an active star-forming mode than lower-mass galaxies. This is more significant at higher redshifts, z > 5,

where the time spent in active stages for the higher-mass galaxies lies more than > 2σ away from the value for the lower-mass galaxies. The active star-forming stages for higher-mass galaxies have an average duration of ≈50–70 Myr. Lower-mass galaxies spend ≈80–130 Myr in this active mode.

These timescales in the active star-forming stage for high-redshift galaxies are similar to those predicted from recent predictions for (so-called) “feedback-free starbursts” (FFBs; see A. Dekel et al. 2023; Z. Li et al. 2024). FFBs are predicted to occur when the gas densities of galaxies lead to sufficiently short freefall times such that galaxies can form bursts of stars on timescales shorter than the lifetime of massive O stars. This FFB model also predicts that galaxies at z > 10 experience episodes of star formation that span ~80 Myr timescales. A. Dekel et al. (2023) demonstrate that feedback from supernova and stellar winds typically does not begin affecting star formation in neighboring gas clouds for ~2–3 Myr after the initial burst, and subsequent bursts are fed by continuous gas accretion over a ~80 Myr timescale (A. Dekel et al. 2023). Our measurement of the duty cycle suggests that if the FFB mechanism is the limiting physical mechanism for star
<table><caption>**Table 3**
Duty Cycle of Star Formation for Galaxies</caption><thead><tr><th rowspan="3">Redshift range</th><th rowspan="3">Δt (Myr)</th><th colspan="6">Galaxies with Stellar Mass</th></tr><tr><th colspan="3">9.3 ≤ log M*/M⊙ ≤ 10.3</th><th colspan="3">8.7 ≤ log M*/M⊙ < 9.3</th></tr><tr><th>N</th><th>f (log R<sub>SFR</sub> > 0)</th><th>t = f × Δt (Myr)</th><th>N</th><th>f (log R<sub>SFR</sub> > 0)</th><th>t = f × Δt (Myr)</th></tr></thead><tbody><tr><td>4.5 &lt; z ≤ 5</td><td>172</td><td>19</td><td>39 ± 4%</td><td>67 ± 7</td><td>137</td><td>47 ± 6%</td><td>81 ± 11</td></tr><tr><td>5 &lt; z ≤ 6</td><td>273</td><td>27</td><td>23 ± 8%</td><td>61 ± 21</td><td>262</td><td>39 ± 5%</td><td>107 ± 13</td></tr><tr><td>6 &lt; z ≤ 7</td><td>210</td><td>14</td><td>33 ± 8%</td><td>70 ± 16</td><td>122</td><td>57 ± 2%</td><td>120 ± 5</td></tr><tr><td>7 &lt; z ≤ 9</td><td>307</td><td>7</td><td>21 ± 13%</td><td>45 ± 27</td><td>51</td><td>63 ± 7%</td><td>132 ± 15</td></tr></tbody></table>

Note. Δt is the cosmic time spanned by each redshift range. f is the duty cycle when galaxies are in the “star-forming” mode and t (log R > 0) is the average duration of time galaxies spend in active periods of star formation. We show all timescales in megayears.

<figure><img src="image_16.png" /><figcaption>**Figure 12.** The fraction of galaxies that have ratios $R_{SFR} \equiv SFR_{10}/SFR_{100}$ greater than unity, $f(\log R_{SFR} > 0)$, which is an estimate of the “duty cycle” of star formation in galaxies. The markers show galaxies in different redshift bins, colored as in Figure 11, going from 4 < z < 5 (yellow) to 7 < z < 9 (purple), which are labeled for clarity. The errors shown are estimated using a bootstrapping method in which we resample each galaxy (at 4.5 ≤ z ≤ 9) and remeasure the fractions to account for uncertainties in redshift, stellar mass, and SFR. At all redshifts, the higher-mass samples have a higher fraction of galaxies with $\log R_{SFR} < 0$, but the change from the lower-mass to higher-mass galaxies is the most dramatic for galaxies at higher redshifts.</figcaption></figure>

formation at these epochs, then it should continue throughout the EoR, at least down to z ~ 6. One possible way to test these predictions is to measure accurate gas densities and gas fractions of galaxies at these redshifts.

The duty cycle for active star formation, $\log R_{SFR} > 0$, estimated in this study appears to be independent of redshift over the range 4.5 < z < 9. The average time spent in a period of active star formation is ≈ 65 Myr for higher-mass galaxies and ≈ 110 Myr for lower-mass galaxies. As noted above, this is similar to the predictions for the FFB scenario (though alternative theories exist; see Section 5.3 below). The FFB scenario also predicts that galaxies should spend time in a “lull” between active star-forming stages, where we predict these lulls should last from ~100 to 250 Myr based on the duty-cycle values we derive (see Table 3).

One important finding from our work is that the active periods of star formation are shorter, and that these lull periods are longer for the higher-mass galaxies (9.3 ≤ log M*/M⊙ ≤ 10.3) than lower-mass galaxies. However, we also find that these lull periods become shorter with decreasing redshift as the duty cycles increase (see Table 3). This naturally leads to more smooth star formation if the time period of the “lulls” continue to shorten as cosmic time increases (i.e., moving toward lower redshifts) such that the duty cycle trends toward unity. This accounts for more stochastic star formation at high redshifts, followed by more smooth star formation at lower redshifts. This is similar to the conclusion of L. Ciesla et al. (2024), but we find the transition from bursts to smooth star formation is more gradual and continues throughout the redshift range here, 4.5 < z < 9. Any physical mechanism or mechanisms regarding gas accretion, star formation, and feedback in these galaxies must account for these observations.

### 5.3. Physical Implications from Simulations

Studies of cosmological models and simulations have identified multiple pathways that fuel and regulate star formation during and after the EoR. The main physical picture involves the growth of dark-matter halos, gas accretion and cooling, star formation, and feedback. Current theoretical studies have yet to reach convergence on the physical processes and the strength of various effects. Regarding feedback, many studies rely on AGN, which can effectively disrupt gas accretion and drive cold gas from galaxies, causing a cessation of star formation (e.g., P. F. Hopkins et al. 2010). However, S. M. Wilkins et al. (2023) demonstrated that AGN feedback is less effective in reproducing observations toward the end of the EoR (4.5 < z < 6). This is supported by observational studies (e.g., R. Smit et al. 2016) that demonstrate the need for rising SFHs and AGN at z < 5 to reproduce observations. However, at higher redshifts these mechanisms are unable to reproduce the abundance of galaxies observed to have recently shut down their star formation (see, e.g., T. M. Heckman & P. N. Best 2014; S. M. Wilkins et al. 2023). Therefore, AGN may be insufficient to regulate star formation in these galaxies. Moreover, to explain the scatter in the SFRs we find, any AGN feedback would need to be episodic on timescales of ~100–200 Myr.

We can also use predictions from simulations and semi-analytical models to understand the global impact of physical effects. Figure 8 compares predictions from the SC-SAM model (L. Y. A. Yung et al. 2024) and from the FLARES simulation (S. M. Wilkins et al. 2024). Both the SC-SAM and
FLARES simulation include prescriptions for gas accretion and cooling, star formation, and stellar feedback. These all contribute to the scatter of galaxy SFRs on the SFMS. These appear to recover the scatter we measure in the SFRs on 100 Myr timescales ($\sigma_{100}$). As illustrated in Figure 8, we measure $\sigma_{100} \approx 0.2-0.3$ dex, which agrees with both the SC-SAM and FLARES predictions for the scatter on these timescales. However, both the SC-SAM and FLARES predictions fall short of the scatter we measure on 10 Myr timescales, where we measure $\sigma_{10} \approx 0.4-0.6$ dex while the predictions for the scatter on 10 Myr timescales is very similar to the predictions for this on 100 Myr timescales. This suggests that the timescales of physical effects pertaining to galaxy star formation are too long in these models and need to vary more quickly to account for our observations, and may require additional gas cloud physics (S. Tacchella et al. 2020). One possibility is the effect of FFBs, which are predicted to act on shorter timescales under the physical conditions expected for these high-redshift galaxies (A. Dekel et al. 2023 and Z. Li et al. 2024, and discussion in Section 5.2).

Other studies have demonstrated that multichannel stellar feedback processes are sufficient in regulating star formation and impact the sSFRs of galaxies. For example, C. A. Mason et al. (2023) demonstrate these processes are able to reproduce the evolution of the UVLF at high redshifts (see also X. Shen et al. 2023). Additionally, S. M. Wilkins et al. (2023) find that the effects of stellar feedback cause more massive star-forming galaxies at these redshifts to have lower sSFRs than lower-mass galaxies. These predictions are consistent with our findings (Figure 9).

## 6. Conclusions

In this paper, we study the relation of galaxy stellar mass and SFR at high redshifts. We focus on trends at $4.5 < z < 9$, and extending this in some cases to $9 < z < 12$. We study the SEDs of $>1800$ galaxies using multiwavelength data from JWST/NIRCam from CEERS, combined with legacy HST/ACS and WFC3 imaging from CANDELS. We model the SEDs with stellar population models that span a range of parameters, and include nebular emission and flexible SFHs. We show that the multiwavelength data set is sensitive to changes in the nebular emission lines of galaxies, which track the SFRs on timescales of 10 Myr (the lifetime of O-type stars), and to changes to the relative amount of rest-UV emission, which tracks the SFRs on timescales of 100 Myr. This allows us to compare the SFRs on short, 10 Myr timescales compared to longer, 100 Myr timescales.

We then use these results to measure and to analyze the evolution of the SFMS out to the redshifts. We implement an MCMC routine utilizing a two-part likelihood function that simultaneously measures the slope, normalization, and intrinsic scatter of the SFMS as a function of redshift. We study the distribution of galaxies in SFRs using $SFR_{10}$ and $SFR_{100}$, and in the $sSFR_{10}$ spaces. We use the ratio of $R_{SFR} = SFR_{10}/SFR_{100}$ to estimate evolution of galaxy SFHs and to estimate a duty cycle of star formation, all as functions of redshift and stellar mass. A summary of our primary findings is as follows.

The JWST/NIRCam multiwavelength photometry is well suited for measuring SFRs over 10 Myr timescales. This is primarily true for the sample of galaxies with $z \lesssim 9$, where the photometry covers both the rest-UV and Balmer regions of the observed SEDs, and samples the strength of emission lines in

the rest-optical ([O III], H$\beta$, etc). This allows better constraints on stellar mass, the SFR, and the SFH in galaxies. Figures 2–4 show examples of these fits and how the SFH constraints depend on the features of the model fits to the multi-wavelength SEDs.

There is minimal evidence that the evolution in the normalization of the SFMS changes from $4.5 < z < 12$ when using SFRs averaged over 10 or 100 Myr timescales (Figure 5).

The scatter of the SFMS is larger when using SFRs averaged on 10 Myr timescales ($\sigma_{10}$) compared to the SFRs averaged on 100 Myr timescales ($\sigma_{100}$). We find $\sigma_{10} \approx 0.3-0.5$ dex compared to $\sigma_{100} \approx 0.1-0.2$ dex, with no significant evolution in redshift for $4.5 < z < 9$ (Figures 6 and 7). The values of $\sigma_{100}$ do not appear to depend on stellar mass. However, while the scatter in the SFRs on 10 Myr timescales is also apparently constant as a function of stellar mass from $4.5 < z < 6$, it shows evidence that it increases with increasing stellar mass at $6 < z < 9$ (Figure 8).

The distributions of specific $SFR_{10}$ in galaxies evolve with both mass and redshift, but the $sSFR_{10}$ distributions for higher-mass galaxies are lower than those for lower-mass galaxies by $\approx 0.3$ dex at all redshifts, $4.5 < z < 9$ (Figure 9).

The distribution of $\log(sSFR_{10})$ for lower-mass galaxies, $8.7 \le \log M_*/M_{\odot} \le 9.3$, is consistent with a single Gaussian with a mean $\log(sSFR)$ that increases slightly with redshift from $\log(sSFR_{10})/\text{yr}^{-1}) = -8.4$ at $z \approx 5$ to $-8.0$ at $z \approx 8$. Higher-mass galaxies, $9.3 < \log M_*/M_{\odot} \le 10.3$, have a larger scatter in their $sSFR_{10}$. This is because their log sSFR distributions are apparently bimodal, which we model as a double Gaussian. Interestingly, the double-Gaussian models have a high sSFR peak (mean) that is consistent with that of the (single-) Gaussian models of lower-mass galaxies, while the second Gaussian model for the higher-mass galaxies is shifted to lower sSFRs (Figure 10). Therefore, the increase in the SFR scatter is not because of an excess of “bursting” galaxies, but rather because of an increase in the fraction of galaxies experiencing a lull in their SFRs at these redshifts.

We use the ratio of $SFR_{10}$ to $SFR_{100}$ to describe trends in galaxy SFHs. We define $\log R_{SFR} = \log SFR_{10} - \log SFR_{100}$ as a ratio of the current SFR (on 10 Myr timescales) to the past average (on 100 Myr timescales). The distribution of $R_{SFR}$ is centered near unity for lower-mass galaxies, while for higher-mass galaxies the distribution shifts to $\log R_{SFR} < 0$ (Figure 11). Generally, this implies that on average lower-mass galaxies have SFHs that are constant on 10–100 Myr timescales, while higher-mass galaxies have SFHs that have declined over these timescales.

We measure a duty cycle for the phase of active star formation in galaxies as the fraction of galaxies that have $\log R_{SFR} > 0$ as a function of stellar mass and redshift. For higher-stellar-mass galaxies, $9.3 < \log M_*/M_{\odot} \le 10.3$, the duty cycle is $\approx 30\%$, and increases from 10%–30% at $7 < z < 9$ to 35%–43% at $4.5 < z < 5$. For lower-mass galaxies, $8.7 \le \log(M_*/M_{\odot}) < 9.3$, the duty cycles are higher, averaging 40%–60% with little change in redshift (Figure 12 and Table 3).

Multiplying the time spanned by the redshift range and the duty cycle leads to an estimate for the duration of time galaxies spend in active star-forming stages (Table 3). For higher-stellar-mass galaxies, $9.3 < \log M_*/M_{\odot} \le 10.3$, this period of active star formation averages $\approx 60$ Myr, with no evidence for evolution over $4.5 < z < 9$. For lower-mass galaxies, the
periods of active star formation are longer, averaging ≈ 110 Myr over this redshift range. At least for the higher-mass galaxies, the fact that the timescales of active star formation are roughly constant but the duty cycle of star formation increases implies the periods of inactivity, when galaxies experience a lull in star formation, must decrease with cosmic time (decreasing redshift). This naturally implies that galaxies transition from a period of higher variability in their star formation at higher redshift to more smooth star formation at lower redshifts when the inactive, lull periods become short. We argue we are witnessing this transition from z ~ 9 to z ~ 4.5.

Future studies of the SFHs and properties of high-redshift galaxies with JWST and other facilities will test some of the findings in this work. Specifically, while current theoretical predictions account for some of our findings (Section 5), they are unable to account for the higher variability in star formation on shorter, 10 Myr timescales. To explain these findings requires new physics, possibly with novel ideas such as FFBs. These theories make predictions for the gas fractions and gas densities of galaxies and their dependence on mass, SFR, and redshift. All of these can be tested with future observations.

## Acknowledgments

We wish to thank all our colleagues in the CEERS collaboration for their hard work and valuable contributions on this project. We also thank all the scientists, engineers, and administrators that made JWST the incredible time machine that it is. J.W.C. acknowledges support from the NASA Headquarters under the Future Investigators in NASA Earth and Space Science and Technology (FINESST) award 22-ASTRO22-0264. J.W.C. would also like to thank Addy Evans, Jacqueline Antwi-Danso and all members of the Extragalactic Astronomy group at Texas A&M University for the countless hours of discussion, soundboarding and support they have provided. C.P. thanks Marsha and Ralph Schilling for generous support of this research. This material is based in part upon work supported by the National Science Foundation under grant AST-2009632. Portions of this research were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing (HPRC; http://hprc.tamu.edu). This work benefited from support from the George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University. This work acknowledges support from the NASA/ ESA/CSA James Webb Space Telescope through the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-03127. Support for program No. JWST- ERS01345 was provided through a grant from the STScI under NASA contract NAS5-03127.

## Appendix A

### Impact of SFR-Stellar Mass Covariance on the Scatter of the SFMS

The covariance between the SFR and $M_*$ for typical galaxies in our sample is correlated with a positive covariance (see Figure 5 and Section 4.1). For star-forming galaxies this covariance arises from the fact that the redder “color” of a galaxy is correlated with both quantities. That is, a redder color corresponds to an increase in the stellar-mass-to-light ratio (and thus increases the stellar mass) and a redder color typically

indicates higher levels of dust attenuation, which increases the SFR. We therefore test how this covariance affects the measurement of scatter along the SFMS. First, we test the fidelity of our two-function likelihood method described in Section 4.1 by applying the methodology with scatter only in the dependent variable, i.e., the SFR. Second, we test how well we are able to recover the scatter when including the effects of the covariance between the dependent and independent variables, i.e., SFR and stellar mass.

We use a model where galaxies have a linear relation between the variables such that

$$y_i = mx_i + b + \epsilon_i, \quad (A1)$$

where the index $i = \{1..N\}$ represents the number of data points (e.g., N galaxies), and $x_i$ is the independent variable and $y_i$ is the dependent variable. Here, $m$ and $b$ are the slope and normalization of the relation for the sample, and $\epsilon_i \sim N(0, \sigma^2)$ is the amount of scatter experienced by datum $i$ given a sample scatter of $\sigma$. We then run a MCMC, implemented with EMCEE (D. Foreman-Mackey et al. 2013), to infer $\hat{\theta} = \{\hat{m}, \hat{b}, \hat{\sigma}\}$, the slope, normalization, and intrinsic scatter, respectively. Note that we adopt the notation here where the “hat” on variables, i.e., $\hat{y}$, distinguishes variables that are estimates compared to true quantities. The symbol “~” means “is distributed as,” i.e., $y_i' \sim N(0, \hat{\sigma}^2)$ reads “$y_i'$ is distributed as a normal distribution.”

Here, we want to separate the slope and normalization from the intrinsic scatter in order to accurately describe the absolute intrinsic scatter, or the intrinsic scatter that reproduces the distribution $|N(0, \sigma^2)|$, while simultaneously inferring the full posterior space for $\hat{\theta}$. Because we are assuming that the variables are distributed along a linear relation and the scatter about that relation is distributed normally, this requires both a linear likelihood function and a $\chi^2$-likelihood function to accurately describe the distribution of the data about the linear relation. This is easily achieved by separating the two distributions and calculating their likelihoods separately.

Therefore, for each iteration of the MCMC fit, we split Equation (A1) to first derive the slope and normalization such that

$$y = \hat{m}x + \hat{b}, \quad (A2)$$

where $x$ and $y$ in our case contains the mass and SFR values for the full set of galaxies. We then calculate $\hat{y} = \hat{m}x + \hat{b}$, which is the estimate of the SFRs given the parameters $\hat{m}$ and $\hat{b}$. In the second step, we calculate

$$y - \hat{y} = y' - N(0, \hat{\sigma}^2). \quad (A3)$$

This is important as the measured estimates are used to infer the scatter in Equation (A3). Equation (A3) then describes the estimated measurement of the intrinsic scatter ($\hat{\sigma}$) along the linear relation. Subtracting off the linear relation leaves only points scattered about $y = 0$, and we model as a normal distribution.

The likelihood equations for these two functions then become

$$\mathcal{L}_1(\hat{m}, \hat{b} | x) = \prod_{i=1}^{N} \hat{m}x_i + \hat{b}, \quad (A4)$$

$$\mathcal{L}_2(\hat{\sigma} | y') \propto \prod_{i=1}^{N} \frac{-(y_i' - \epsilon_i)^2}{\hat{\sigma}^2}, \quad (A5)$$
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Simulated SFR–$M_*$ relations assuming different values of the intrinsic scatter. Blue circles show the simulated data and the solid black line shows the median best fit from the MCMC routine, with 1 and $2\sigma$ shown as the gray-shaded regions. The blue ellipses show the 1 and $2\sigma$ covariance matrix from which the simulated data are drawn.</figcaption></figure>

where $\epsilon_i \sim N(0, \hat{\sigma}^2)$, and we assume a $\chi^2$-likelihood under the assumption that the scatter about the linear function follows a normal distribution. Taking the natural logarithm ($\ln \mathcal{L} = \ell$) of the two functions provides the log-likelihood function used in our MCMC routine. They have the form

$$ \ell_1(\hat{m}, \hat{b}|\mathbf{x}) = \ln[\prod_{i=1}^N \hat{m}x_i + \hat{b}] = \sum_{i=1}^N \ln(\hat{m}x_i + \hat{b}), \quad (A6) $$

$$ \begin{aligned} \ell_2(\hat{\sigma}|\mathbf{y}') &= \ln\left[\frac{1}{\hat{\sigma}^{2N}} \prod_{i=1}^N (y'_i - \epsilon_i)^2\right] \\ &= -2N \ln \hat{\sigma} - \sum_{i=1}^N \ln[(y'_i - \epsilon_i)^2], \end{aligned} \quad (A7) $$

The final form of the likelihood function must then be the sum of Equations (A6) and (A7),

$$ \ell(\theta|\mathbf{x}, \mathbf{y}') = \ell_1 + \ell_2 = \sum_{i=1}^{N} \ln(\hat{m}x_i + \hat{b}) - \sum_{i=1}^{N} \ln[(y'_i - \epsilon_i)^2) - 2N \ln \hat{\sigma}], \quad (A8) $$

and is the total likelihood function that we implement in our MCMC routine.

To test the ability of this routine to measure the SFR–$M_*$ relation, we start by generating $N = 20,000$ random data points for $\tilde{x}'$ from a uniform distribution with $8 \le x \le 11$. We then compute $y_i = \alpha(x_i - M_0) + \beta$, where $\alpha = 1$, $\beta = 0$, and $M_0 = 8.25$, to allow the relation to span similar values as our data. This provides the set of $x_i, y_i$ for the $i = 1, \dots, N$ points. To introduce uncertainties on “stellar mass” and “SFR,” we generate a correlation matrix using the $P(M_*/M_\odot, \text{SFR})$ joint posterior space of a galaxy with log(SFR) ~ 0 and log($M_*/M_\odot$) ~ 9. We use this to generate 20,000 perturbations with means $\{\mu_x, \mu_y\} = (0,0)$ and $\{\sigma_x, \sigma_y\} = (0.2,0.2)$, $\{dx_i, dy_i\}$. We also generate 20,000 values for $\epsilon$ from a normal distribution, $\epsilon \sim N(0, 0.1)$, to simulate intrinsic scatter along the linear relation. Adding this to our simulated set of data, we get $\{\tilde{x}_i, \tilde{y}_i\} = \{x_i + dx_i, y_i + dy_i + \epsilon_i\}$. Finally, in order to ensure that we can accurately recover multiple values for intrinsic scatter along a linear relation, we scale the perturbations ($\{\epsilon_i\}$) to create two additional spaces with $\sigma = 0.2$ and $\sigma = 0.3$. We test our methodology on all three test cases, $\sigma = 0.1, 0.2$, and $0.3$, and the “data” are shown in Figure 13.

For each iteration of the MCMC simulation, we draw 2000 samples from the set of 20,000 simulated “galaxies” to simulate

<table><caption>**Table 4**

Results of Dual-likelihood MCMC Testing</caption><thead><tr><th>σ<sub>True</sub> = 0.1</th><th>σ<sub>True</sub> = 0.2</th><th>σ<sub>True</sub> = 0.3</th></tr></thead><tbody><tr><td>α̂ = 1.01 ± 0.02</td><td>α̂ = 1.00 ± 0.02</td><td>α̂ = 1.01 ± 0.02</td></tr><tr><td>β̂ = -0.01 ± 0.02</td><td>β̂ = 0.00 ± 0.03</td><td>β̂ = -0.01 ± 0.03</td></tr><tr><td>σ̂ = 0.07 ± 0.01</td><td>σ̂ = 0.14 ± 0.01</td><td>σ̂ = 0.21 ± 0.02</td></tr></tbody></table>

sampling the SFR–$M_*$ posterior spaces of a real sample. Additionally, any simulated data with log(SFR) ≤ –0.4 are excluded from the calculations. We do this to simulate an apparent selection bias against galaxies with low SFRs. We then compute the likelihood function for each walker, as described above. Figure 13 shows the results of our fitting method, with a sample of 2000 simulated galaxies to demonstrate the scatter produced by the covariance matrix and intrinsic scatter for $\sigma = 0.1, 0.2$, and 0.3 (left to right, respectively). The MCMC routine appears to trace the two-dimensional distribution of simulated galaxies very well, and we list the results in Table 4.

Figure 14 shows corner plots for each of the three test cases. As can be seen from both Figure 14 and Table 4, we underestimate $\sigma$ by ~30% in all cases. Additionally, we are able to recover the $\alpha$ and $\beta$ in all cases. Although the intrinsic scatter along the relation is underestimated, the constant offset, $\sigma_{MCMC} = 0.7\sigma_{\text{true}}$, will not affect the results of our analysis. We could take the $\epsilon(z)$ values measured in our analysis and increase them by a factor of $1/0.7 \approx 1.4$ without affecting the results of this paper. However, we only note that our results could be underestimated by this factor of ~1.4 and do not apply the correction to any of the results presented in our analysis. In order to address this constant offset between the measured values and the “true” values, we also tested this method using $\alpha_{\text{True}} = 0.8$, the value we measure for the SFMS in this analysis. We find that the constant offset remains, but decreases to $1/0.8 = 1.2$. We conclude that this constant offset is a result of an additional term that is introduced in the $\{dx_i, dy_i + \eta_i\}$ terms, causing an artificial tightening along the axis of the posterior spaces for galaxies. In the case of $\alpha_{\text{True}} = 1$, this artificial tightening effect is greatest because the posterior space for the galaxy we use as our “model” posterior space is parallel
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Corner plots from the MCMC routine showing the relationships between the posteriors of the slope (α), normalization (β), and the intrinsic scatter (σ). From left to right, we show the results assuming σ = 0.1, σ = 0.2, and σ = 0.3.</figcaption></figure>

to the simulated SFMS. Therefore, given the more randomly distributed orientation of the galaxy posterior spaces in our sample, this effect will be minimized through the random sampling of those posterior spaces.

**Appendix B**

**Impact of Selection Band on SFR-Stellar Mass Relation**

For this work, we used the NIRCam F277W photometry to select our fiducial sample of galaxies. This is motivated by the fact that the F277W filter is the most sensitive of the NIRCam bandpasses available (M. B. Bagley et al. 2023; S. L. Finkelstein et al. 2023). It also offers the best option among the redder NIRCam filters as it is mostly devoid of strong nebular emission lines over the redshift range of interest (see discussion below).

To study the impact of the selection band on our results, we repeated our analysis selecting galaxies in the F444W band. This band offers selection at redder wavelengths, potentially more sensitive to “redder” galaxies, but is shallower such that it is less sensitive to bluer galaxies, which constitute the majority of our sample.

We selected galaxies by requiring sources be detected 1 magnitude brighter than the 5σ detection limit in the F444W band, $m_{F444W} \le 28.6$ mag (S. L. Finkelstein et al. 2023). We then repeated the selection steps as outlined in Section 2.2. Overall, the number of sources detected drops by a factor of ~2, as we are using a shallower detection limit. However, we find that the distributions of $m_{F277W}$, $m_{F356W}$, and $m_{F444W}$ for the F444W-selected sample are statistically similar to the fiducial, F277W-selected sample, and the two methods are

consistent with drawing from the same parent sample of galaxies.

We then recalculated the SFR-$M_*$ relations using the F444W-selected sample. We repeated the MCMC analysis and inferred constraints on the slope (α'), normalizations (β(z)), and the intrinsic scatters (σint(z)) for all galaxies in the F444W selection. Table 5, comparable to Table 2, shows the inferred parameters for the F444W-selected SFR-$M_*$ relation. The results are qualitatively unchanged from the F277W-selected sample. The primary difference is that the inferred slope (α) for the F444W-selected sample is lower, i.e., the relation is flatter, than the slope derived from the F277W-selected sample (see Table 2). However, the results are consistent within their 1σ uncertainties and follow the same qualitative trends for both the SFR$_{10}$ and SFR$_{100}$ analyses. We therefore conclude that the choice of bandpass does not significantly alter our findings.

We also find higher scatter, σint, when using the F444W selection for both SFR$_{10}$ and SFR$_{100}$. We conclude the main reason for this is that the F444W selection keeps redder objects (in $m_{F277W} - m_{F444W}$), but misses bluer objects that are detected in the $m_{F277W}$ selection. The latter are “main-sequence” galaxies on the SFR-$M_*$ relation and constitute the majority of the F277W-selected sample. As a result, the scatter increases for the F444W-selected sample because we miss selecting the “typical” (average) star-forming galaxies. Nevertheless, the main trends presented here do not change in respect to the scatter in SFR$_{10}$ being higher than that in SFR$_{100}$. Furthermore, using the F444W selection, we see that the scatter in SFR$_{10}$ increases with redshift from 4.5 < z < 5 to 7 < z < 9, and the scatter in SFR$_{100}$ does not, qualitatively consistent with the trends found with the F277W-selected sample.
<table><caption>**Table 5**
Results of the MCMC Analysis for SFR<sub>10</sub> (top) and SFR<sub>100</sub> (bottom) for the F444W-selected Sample of Galaxies</caption><thead><tr><th>Redshift</th><th>Sample Size</th><th>α</th><th>Median log(M<sub>*</sub>)</th><th>β(z)</th><th>σ(z)</th></tr></thead><tbody><tr><td colspan="6">MCMC Results Measured With SFR<sub>100</sub></td></tr><tr><td>All galaxies</td><td>1348</td><td>0.73 ± 0.28</td><td>8.74</td><td></td><td></td></tr><tr><td>4.5 &lt; z &le; 5</td><td>253</td><td></td><td></td><td><strong>0.27</strong> ± 0.08</td><td><strong>0.37</strong> ± 0.06</td></tr><tr><td>5 &lt; z &le; 6</td><td>587</td><td></td><td></td><td><strong>0.39</strong> ± 0.04</td><td><strong>0.34</strong> ± 0.04</td></tr><tr><td>6 &lt; z &le; 7</td><td>190</td><td></td><td></td><td><strong>0.47</strong> ± 0.08</td><td><strong>0.28</strong> ± 0.06</td></tr><tr><td>7 &lt; z &le; 9</td><td>113</td><td></td><td></td><td><strong>0.49</strong> ± 0.12</td><td><strong>0.29</strong> ± 0.08</td></tr><tr><td>9 &lt; z &le; 12</td><td>35</td><td></td><td></td><td>0.72 ± 0.21</td><td><strong>0.33</strong> ± 0.14</td></tr><tr><td colspan="6">MCMC Results Measured With SFR<sub>10</sub></td></tr><tr><td>All galaxies</td><td>1348</td><td>0.65 ± 0.14</td><td>8.74</td><td></td><td></td></tr><tr><td>4.5 &lt; z &le; 5</td><td>253</td><td></td><td></td><td>0.71 ± 0.08</td><td><strong>0.45</strong> ± 0.07</td></tr><tr><td>5 &lt; z &le; 6</td><td>587</td><td></td><td></td><td><strong>0.70</strong> ± 0.06</td><td><strong>0.53</strong> ± 0.04</td></tr><tr><td>6 &lt; z &le; 7</td><td>190</td><td></td><td></td><td><strong>0.95</strong> ± 0.12</td><td><strong>0.56</strong> ± 0.11</td></tr><tr><td>7 &lt; z &le; 9</td><td>113</td><td></td><td></td><td><strong>0.27</strong> ± 0.14</td><td><strong>0.84</strong> ± 0.17</td></tr><tr><td>9 &lt; z &le; 12</td><td>35</td><td></td><td></td><td>– <strong>1.38</strong> ± 0.22</td><td><strong>0.50</strong> ± 0.12</td></tr></tbody></table>

Note. Fitted parameters for Equation (5) are shown in bold (β(z) and σ²(z)), and we list α and the median stellar mass (M<sub>*</sub>) for clarity.

## ORCID iDs

Justin W. Cole https://orcid.org/0000-0002-6348-1900

Casey Papovich https://orcid.org/0000-0001-7503-8482

Steven L. Finkelstein https://orcid.org/0000-0001-8519-1130

Micaela B. Bagley https://orcid.org/0000-0002-9921-9218

Mark Dickinson https://orcid.org/0000-0001-5414-5131

Kartheik G. Iyer https://orcid.org/0000-0001-9298-3523

L. Y. Aaron Yung https://orcid.org/0000-0003-3466-035X

Laure Ciesla https://orcid.org/0000-0003-0541-2891

Ricardo O. Amorín https://orcid.org/0000-0001-5758-1000

Pablo Arrabal Haro https://orcid.org/0000-0002-7959-8783

Rachana Bhatawdekar https://orcid.org/0000-0003-0883-2226

Antonello Calabrò https://orcid.org/0000-0003-2536-1614

Nikko J. Cleri https://orcid.org/0000-0001-7151-009X

Alexander de la Vega https://orcid.org/0000-0002-6219-5558

Avishai Dekel https://orcid.org/0000-0003-4174-0374

Ryan Endsley https://orcid.org/0000-0003-4564-2771

Eric Gawiser https://orcid.org/0000-0003-1530-8713

Mauro Giavalisco https://orcid.org/0000-0002-7831-8751

Nimish P. Hathi https://orcid.org/0000-0001-6145-5090

Michaela Hirschmann https://orcid.org/0000-0002-3301-3321

Benne W. Holwerda https://orcid.org/0000-0002-4884-6756

Jeyhan S. Kartaltepe https://orcid.org/0000-0001-9187-3605

Anton M. Koekemoer https://orcid.org/0000-0002-6610-2048

Ray A. Lucas https://orcid.org/0000-0003-1581-7825

Sara Mascia https://orcid.org/0000-0002-9572-7813

Bahram Mobasher https://orcid.org/0000-0001-5846-4404

Pablo G. Pérez-González https://orcid.org/0000-0003-4528-5639

Giulia Rodighiero https://orcid.org/0000-0002-9415-2296

Kaila Ronayne https://orcid.org/0000-0001-5749-5452

Sandro Tacchella https://orcid.org/0000-0002-8224-4505

Benjamin J. Weiner https://orcid.org/0000-0001-6065-7483

Stephen M. Wilkins https://orcid.org/0000-0003-3903-6935

## References

Álvarez-Márquez, J., Burgarella, D., Heinis, S., et al. 2016, A&A, 587, A122

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

Asada, Y., Sawicki, M., Abraham, R., et al. 2024, MNRAS, 527, 11372

Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJL, 946, L12

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58

Barrow, K. S. S. 2020, MNRAS, 491, 4509

Becker, G. D., D'Aloisio, A., Christenson, H. M., et al. 2021, MNRAS, 508, 1853

Bird, S., Ni, Y., Di Matteo, T., et al. 2022, MNRAS, 512, 3703

Boselli, A., Fossati, M., Roediger, J., et al. 2023, A&A, 669, A73

Bowler, R. A. A., Inami, H., Sommovigo, L., et al. 2024, MNRAS, 527, 5808

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

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019, ApJ, 873, 44

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

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

Chabrier, G. 2003, PASP, 115, 763

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

Chworowsky, K., Finkelstein, S. L., Spilker, J. S., et al. 2023, ApJ, 951, 49

Chworowsky, K., Finkelstein, S. L., Boylan-Kolchin, M., et al. 2024, AJ, 168, 113

Ciesla, L., Elbaz, D., Ilbert, O., et al. 2024, A&A, 686, A128

Conroy, C. 2013, ARA&A, 51, 393

Curtis-Lake, E., Chevallard, J., Charlot, S., & Sandles, L. 2021, MNRAS, 503, 4855

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

Darvish, B., Mobasher, B., Sobral, D., Scoville, N., & Aragon-Calvo, M. 2015, ApJ, 805, 121

Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265

de Graaff, A., Brammer, G., Weibel, A., et al. 2024a, arXiv:2409.05948

de Graaff, A., Setton, D. J., Brammer, G., et al. 2024b, NatAs, in press

Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, 523, 3201

Dickinson, M., Amorin, R., Arrabal Haro, P., et al. 2024, The CANDELS-Area Prism Epoch of Reionization Survey (CAPERS), JWST Proposal. Cycle 3 ID. #6368

Dome, T., Tacchella, S., Fialkov, A., et al. 2024, MNRAS, 527, 2139
Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485, 4817

Eldridge, J. J., & Stanway, E. R. 2022, ARA&A, 60, 455

Eldridge, J. J., Stanway, E. R., Xiao, L., et al. ed. J. J. Eldridge et al. 2017, The Lives and Death-Throes of Massive Stars, 329, 396

Endsley, R., Stark, D. P., Chevillard, J., & Charlot, S. 2021, MNRAS, 5229

Endsley, R., Stark, D. P., Bouwens, R. J., et al. 2022, MNRAS, 517, 5642

Endsley, R., Stark, D. P., Whitler, L., et al. 2024, MNRAS, 533, 1111

Estrada-Carpenter, V., Papovich, C., Momcheva, I., et al. 2023, ApJ, 951, 115

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, RMXAA, 53, 385

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

Finkelstein, S. L., Bagley, M. B., Arrabal Haro, P., et al. 2022, ApJL, 940, L55

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

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

Finlator, K., Davé, R., Papovich, C., & Hernquist, L. 2006, ApJ, 639, 672

Finlator, K., Oppenheimer, B. D., & Davé, R. 2011, MNRAS, 410, 1703

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fukugita, M., Ichikawa, T., Gunn, J. E., et al. 1996, AJ, 111, 1748

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

Gelli, V., Salvadori, S., Ferrara, A., Pallottini, A., & Carniani, S. 2023, ApJL, 954, L11

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589

Hopkins, P. F., Younger, J. D., Hayward, C. C., Narayanan, D., & Hernquist, L. 2010, MNRAS, 402, 1693

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Iyer, K., Gawiser, E., Davé, R., et al. 2018, ApJ, 866, 120

Iyer, K. G., Gawiser, E., Faber, S. M., et al. 2019, ApJ, 879, 116

Kannan, R., Smith, A., Garaldi, E., et al. 2022, MNRAS, 514, 3857

Katz, H., Rosdahl, J., Kimm, T., et al. 2023, OJAp, 6, 44

Kennicutt, R. C. 1998, ApJ, 498, 541

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

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Labbé, I., van Dokkum, P., Nelson, E., et al. 2023, Natur, 616, 266

Larson, R. L., Hutchison, T. A., Bagley, M., et al. 2023, ApJ, 958, 141

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

Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170

Li, Z., Dekel, A., Sarkar, K. C., et al. 2024, A&A, 690, A108

Looser, T. J., D'Eugenio, F., Maiolino, R., et al. 2023, arXiv:2306.02470

Looser, T. J., D'Eugenio, F., Maiolino, R., et al. 2024, Natur, 629, 53

Lovell, C. C., Vijayan, A. P., Thomas, P. A., et al. 2021, MNRAS, 500, 2127

Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2

Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Marchesini, D., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2009, ApJ, 701, 1765

Maseda, M. V., Lewis, Z., Matthee, J., et al. 2023, ApJ, 956, 11

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

Mobasher, B., Capak, P., Scoville, N. Z., et al. 2007, ApJS, 172, 117

Mobasher, B., Dahlen, T., Ferguson, H. C., et al. 2015, ApJ, 808, 101

Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, ApJS, 220, 10

Narayanan, D., Lower, S., Torrey, P., et al. 2024, ApJ, 961, 73

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJL, 660, L43

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

Old, L. J., Balogh, M. L., van der Burg, R. F. J., et al. 2020, MNRAS, 493, 5987

Olivier, G. M., Berg, D. A., Chisholm, J., et al. 2022, ApJ, 938, 16

Pacífici, C., Kassin, S. A., Weiner, B. J., et al. 2016, ApJ, 832, 79

Papovich, C., Dickinson, M., & Ferguson, H. C. 2001, ApJ, 559, 620

Papovich, C., Finkelstein, S. L., Ferguson, H. C., Lotz, J. M., & Giavalisco, M. 2011, MNRAS, 412, 1123

Papovich, C., Moustakas, L. A., Dickinson, M., et al. 2006, ApJ, 640, 92

Papovich, C., Labbé, I., Quadri, R., et al. 2015, AJ, 803, 26

Papovich, C., Cole, J. W., Yang, G., et al. 2023, ApJL, 949, L18

Pforr, J., Maraston, C., & Tonini, C. 2012, MNRAS, 422, 3285

Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526

Reddy, N. A., Pettini, M., Steidel, C. C., et al. 2012, ApJ, 754, 25

Rieke, M. J., Robertson, B., Tacchella, S., et al. 2023, ApJS, 269, 16

Salmon, B., Papovich, C., Finkelstein, S. L., et al. 2015, ApJ, 799, 183

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Shen, X., Vogelsberger, M., Boylan-Kolchin, M., Tacchella, S., & Kannan, R. 2023, MNRAS, 525, 3254

Shivaei, I., Reddy, N. A., Shapley, A. E., et al. 2015, ApJ, 815, 98

Simmonds, C., Tacchella, S., Hainline, K., et al. 2024, MNRAS, 527, 6139

Smit, R., Bouwens, R. J., Labbé, I., et al. 2016, ApJ, 833, 254

Sobral, D., Matthee, J., Darvish, B., et al. 2018, MNRAS, 477, 2817

Somerville, R. S., Popping, G., & Trager, S. C. 2015, MNRAS, 453, 4337

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Stefanon, M., Bouwens, R. J., Illingworth, G. D., et al. 2022, ApJ, 935, 94

Stefanon, M., Yan, H., Mobasher, B., et al. 2017, ApJS, 229, 32

Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123

Strom, A. L., Steidel, C. C., Rudie, G. C., et al. 2017, ApJ, 836, 164

Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., et al. 2023, ApJL, 955, L35

Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 457, 2790

Tacchella, S., Forbes, J. C., & Caplar, N. 2020, MNRAS, 497, 698

Tacchella, S., Conroy, C., Faber, S. M., et al. 2022, ApJ, 926, 134

Tacchella, S., Johnson, B. D., Robertson, B. E., et al. 2023, MNRAS, 522, 6236

Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJL, 754, L29

Whitler, L., Stark, D. P., Endsley, R., et al. 2023, MNRAS, 519, 5859

Wilkins, S. M., Feng, Y., Di Matteo, T., et al. 2017, MNRAS, 469, 2517

Wilkins, S. M., Vijayan, A. P., Lovell, C. C., et al. 2023, MNRAS, 518, 3935

Wilkins, S. M., Lovell, C. C., Irodotou, D., et al. 2024, MNRAS, 527, 7965

Williams, C. C., Alberts, S., Ji, Z., et al. 2024, ApJ, 968, 34

Xue, Y. Q., Luo, B., Brandt, W. N., et al. 2016, ApJS, 224, 15

Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Popping, G., & Davé, R. 2019, MNRAS, 483, 2983

Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Wilkins, S. M., & Gardner, J. P. 2024, MNRAS, 527, 5929