# The Star-Forming Main Sequence in JADES and CEERS at $z>1.4$: Investigating
the Burstiness of Star Formation

Leonardo Clarke University of California, Los Angeles, 475 Portola Plaza, Los
Angeles, CA, 90095, USA  Alice E. Shapley University of California, Los
Angeles, 475 Portola Plaza, Los Angeles, CA, 90095, USA  Ryan L. Sanders NHFP
Hubble Fellow  University of Kentucky, 506 Library Drive, Lexington, KY,
40506, USA  Michael W. Topping Steward Observatory, University of Arizona, 933
N Cherry Avenue, Tucson, AZ 85721, USA  Gabriel B. Brammer Cosmic Dawn Center
(DAWN), Denmark  Niels Bohr Institute, University of Copenhagen, Jagtvej 128,
DK-2200 Copenhagen N, Denmark  Trinity Bento  University of California, Los
Angeles, 475 Portola Plaza, Los Angeles, CA, 90095, USA  Naveen A. Reddy
Department of Physics & Astronomy, University of California, Riverside, 900
University Avenue, Riverside, CA 92521, USA  Emily Kehoe  University of
California, Los Angeles, 475 Portola Plaza, Los Angeles, CA, 90095, USA

###### Abstract

We have used public JWST/NIRSpec and JWST/NIRCam observations from the CEERS
and JADES surveys in order to analyze the star-forming main sequence (SFMS)
over the redshift range $1.4\leq z<7$. We calculate the star-formation rates
(SFRs) of the galaxy sample using three approaches: Balmer line luminosity,
spectral energy distribution (SED) fitting, and UV luminosity. We find a
larger degree of scatter about the SFMS using the Balmer-based SFRs compared
to the UV-based SFRs. Because these SFR indicators are sensitive to star
formation on different time scales, the difference in scatter may be evidence
of bursty star-formation histories in the early universe. We additionally
compare the H$\alpha$-to-UV luminosity ratio (L(H$\alpha$)/$\nu$Lν,1600) for
individual galaxies in the sample and find that 29%$-$52% of the ratios across
the sample are poorly described by predictions from a smooth star-formation
history. Measuring the burstiness of star formation in the early universe has
multiple significant implications, such as deriving accurate physical
parameters from SED fitting, explaining the evolution of the UV luminosity
function, and providing constraints for sub-grid models of feedback in
simulations of galaxy formation and evolution.

††facilities: JWST, HST††software: Astropy (Astropy Collaboration et al.,
2022), emcee (Foreman-Mackey et al., 2013), jwst (Bushouse et al., 2023),
msaexp (Brammer, 2023), scipy (Virtanen et al., 2020)

##  1 Introduction

The galaxy star-forming main sequence (SFMS) is a relation that exists between
the stellar masses and the star-fomation rates (SFRs) of galaxies. This
relation was first noted in a large sample of galaxies from the Sloan Digital
Sky Survey (SDSS, York et al., 2000; Brinchmann et al., 2004), and its
evolution as a function of redshift has also been studied extensively (e.g.,
Noeske et al., 2007; Whitaker et al., 2012, 2014; Speagle et al., 2014;
Shivaei et al., 2015; Santini et al., 2017). Studies of the galaxy SFMS are
typically aimed at constraining galaxy star-formation histories (SFHs) over
time, and the tightness of the SFMS implies that galaxies of similar masses
generally have similar SFHs characterized by a smooth buildup of stellar mass
over time.

The intrinsic scatter in the SFMS is thought to encapsulate information about
the stocasticity or “burstiness” of galaxy SFHs. State-of-the-art
hydrodynamical simulations of galaxy formation such as FIRE-2 (Hopkins et al.,
2018), Illustris TNG (Springel et al., 2018), and EAGLE (Schaye et al., 2015)
make predictions of galaxy SFHs that depend on sub-grid models of feedback and
star formation. Comparing observational constraints on galaxy SFHs with
simulation predictions that are based on sub-grid models provides an
opportunity to test and improve upon our understanding of how galaxies build
up their stellar mass over time.

In addition to constraining simulations, bursty SFHs have been proposed as a
possible explanation for the excess of galaxies at the bright end of the UV
luminosity function at $z\gtrsim 10$ (e.g., Mason et al., 2023; Mirocha &
Furlanetto, 2023; Sun et al., 2023a, b; Shen et al., 2023). In short, the
burstiness of galaxies may upscatter galaxies of fixed stellar mass to high UV
luminosities, resulting in an excess of UV-bright galaxies at early times.
Characterizing the burstiness of SFHs at early times is also important for
SED-based determinations of stellar mass. It has been shown, for example, that
assuming a constant and/or smooth SFH during SED fitting can underestimate the
stellar masses of star-forming galaxies on the order of $\sim$0.1-0.8 dex
(Domínguez et al., 2015; Endsley et al., 2023b; Whitler et al., 2023; Wang et
al., 2024). Additionally, biases can be introduced into SED-fitting procedures
due to recent star-formation episodes outshining older stellar populations
(e.g., Shapley et al., 2005; Narayanan et al., 2024)

There are several observational signatures that give insight into the SFHs of
galaxies such as the intrinsic scatter in the SFMS, the H$\alpha$-to-UV ratio,
and the deviation in log(SFR) from the SFMS (e.g., Weisz et al., 2012; Atek et
al., 2022; Asada et al., 2024). A larger intrinsic scatter, for example, may
indicate that galaxies undergo significant excursions from the SFMS over the
course of their growth in the form of frequent bursts. Similarly, H$\alpha$
and UV emission are sensitive to SFRs on timescales of $\sim$5 Myr and
$\sim$100 Myr, respectively, and their ratio can constrain the rate of
fluctuations in the SFR. These observational signatures serve as metrics that
can be directly compared with theoretical predictions of SFHs, thereby placing
constraints on the proper implementation of sub-grid models of feedback and
star formation.

There is a large body of work that has compared observational signatures of
bursty SFHs with theoretical predictions for galaxies at various redshifts
(e.g., Weisz et al., 2012; Domínguez et al., 2015; Emami et al., 2019; Asada
et al., 2024). It has been found, for example, that there may be a larger
degree of intrinsic scatter in the SFMS for galaxies at low stellar masses
($\log({\rm M}_{*}/{\rm M}_{\odot})\lesssim 8$), though a number of
observational studies find no mass dependence in the intrinsic scatter (e.g.,
Whitaker et al., 2012; Salmon et al., 2015; Schreiber et al., 2015; Kurczynski
et al., 2016). However, many of these previous observational studies rely
solely on photometric data which probes star formation on $>$100 Myr time
scales, or measure H$\alpha$ using photometric excesses which can be less
certain when compared with direct spectroscopic measurement. With the
increased spectroscopic sensitivity and wavelength coverage of the Near-
Infrared Spectrograph (NIRSpec; Jakobsen et al., 2022) on the James Webb Space
Telescope (JWST), studies of the SFMS can now be done routinely for
representative samples of high-redshift galaxies (e.g., Looser et al., 2023;
Simmonds et al., 2024; Narayanan et al., 2024; Ciesla et al., 2023; Curti et
al., 2024).

In this study, we use public JWST data from the Near-Infrared Camera (NIRCam;
Rieke et al., 2023) and JWST/NIRSpec observations from the JWST Advanced Deep
Extragalactic Survey (JADES) and the Cosmic Evolution Early Release Science
(CEERS) survey in conjunction with 3DHST photometric observations to measure
the SFMS and its intrinsic scatter for galaxies at $1.4<z<7$. We measure SFRs
using three different indicators: H$\alpha$ luminosity, UV luminosity, and SED
fitting. The use of different indicators allows us to measure the intrinsic
scatter in the SFR on different time scales, and therefore gain insights into
the burstiness of galaxy star-formation histories in the early universe.

Throughout this study, we adopt $H_{0}=70\ {\rm km\ s^{-1}\ Mpc}$,
$\Omega_{m}=0.3$, and $\Omega_{\Lambda}=0.7$ as cosmological parameters, and
we assume a Chabrier (2003) initial mass function. Additionally, we use the
solar abundances reported by Asplund et al. (2009) with $12+\log({\rm
O/H})_{\odot}=8.69$ and the solar metallicity as $Z_{\odot}=0.014$.

##  2 Observations and Measurements

###  2.1 Observations

####  2.1.1 CEERS

We use public NIRSpec data from the CEERS program (Program ID: 1345
Finkelstein et al., 2023; Arrabal Haro et al., 2023). Our analysis is based on
6 NIRSpec Micro-Shutter Assembly (MSA) pointings in the AEGIS field, which
utilized the grating/filter combinations of G140M/F100LP, G235M/F170LP, and
G395M/F290LP, providing a spectral resolution of $R\sim 1000$ over the
wavelength range approximately $1-5\mu$m. Each pointing was observed for a
total of 3107 sec in each grating/filter combination. A 3-point nod pattern
was adopted for each observation along each 3-shutter MSA slit. In total, the
6 MSA pointings contained 321 slits and 318 distinct targets.

####  2.1.2 JADES

In addition to CEERS, we made use of publicly released data from the JADES
program (Program ID: 1210, Eisenstein et al., 2023) in the GOODS-S
extragalactic legacy field. The full data release consisted of NIRCam imaging
in multiple photometric bands as well as NIRSpec spectra taken with the
PRISM/CLEAR, G140M/F070LP, G235M/F170LP, G395M/F290LP, and G395H/F290LP
grating/filter configurations. Throughout this analysis, we utilized NIRSpec
observations in the G140M/F070LP, G235M/F170LP, and G395M/F290LP
grating/filter combinations.

Observations were conducted over the course of three different visits. In
order to reduce noise introduced by detector defects, the MSA configuration
was shifted to place targets on different slits for each visit. A total of 198
galaxies were observed with the aforementioned grating/filter configuration,
117 of which we were able to measure robust redshifts. Of this subset of 117
galaxies considered in this study, 37 were observed in all three visits, 25
were observed in two visits, and 55 were only observed during a single visit.
For each visit, the exposure time was 2.3 hr in each grating, resulting in a
total exposure time of 6.9 hr per grating for objects observed in all three
visits. A more detailed description of the JADES NIRSpec observations can be
found in Bunker et al. (2023).

###  2.2 Data Reduction

####  2.2.1 CEERS

The CEERS data reduction was performed using the STScI jwst
pipeline.111https://jwst-pipeline.readthedocs.io/en/latest/index.html The
details of the data reduction process can be found in Sanders et al. (2023).
In short, the data were reduced using the pipeline, and exposures were co-
added to produce a set of 310 individual 2D spectra. 1D spectra were then
extracted from the 2D spectra according to the procedure described in Section
2.2.2, resulting in a sample of 252 1D spectra. These spectra were then
corrected for wavelength-dependent slit losses (see Reddy et al., 2023) and
scaled so that the flux densities matched the broadband NIRCam photometry (or
the 3DHST photometry where NIRCam imaging wasn’t available).

####  2.2.2 JADES

We downloaded reduced, two-dimensional (2D) NIRSpec spectra from the DAWN JWST
Archive (DJA)222https://dawn-cph.github.io/dja/spectroscopy/nirspec/ whose
reduction closely resembles the reduction procedure applied to the CEERS
NIRSpec data. The reduction was performed using a custom data reduction
pipeline, MsaExp333https://github.com/gbrammer/msaexp, which takes the Stage 2
data products from the MAST database444https://archive.stsci.edu/hlsp/jades as
inputs and performs the wavelength calibrations, flat-fielding, photometric
calibrations, and a point-source slit-loss correction on each of the
individual exposures. Each of the 2D exposures was then co-added and combined
into a final 2D spectrum for each object. A more detailed explanation of the
reduction of the 2D spectra can be found in Heintz et al. (2023) and Heintz et
al. (2024).

The one-dimensional (1D) spectra were extracted from the 2D spectra by hand
using an optimal extraction aperture (Horne, 1986) consistent with the
procedure described in Sanders et al. (2023). During the extraction procedure,
strong emission lines were identified by eye in the 2D spectra and were fit
with a Gaussian profile to make initial estimates of their redshifts.

Because the exposures were taken over multiple visits, there were uncertain
variations in the position of the target on the slit during each visit. Since
the exposures in each of these visits were combined into a single 2D frame,
deriving a single solution for wavelength-dependent slit losses, accounting
for the extended nature of galaxy targets, became non-trivial. However, we
calculate that the slit-loss correction factors vary less than 5% between
$1-5\ \mu m$. Thus, we do not apply slit-loss corrections to the JADES spectra
beyond the default point-source correction. We instead account for extended-
source slit losses during the flux calibration stage.

In order to obtain a robust absolute flux calibration, we first ensured that
the relative flux calibration between gratings was accurate before scaling the
spectra in all gratings to match the available photometry. To ensure an
accurate relative flux calibration, we applied a multiplicative scaling factor
to the spectra in both the G140M and the G395M gratings, leaving the spectra
in the G235M grating fixed. The relative scaling was applied under the
assumption that all targets with the same exposure time should reach the same
relative sensitivity between gratings. In other words, the median value of the
error spectrum in the G140M grating ($med(\sigma_{G140M})$) relative to that
for the G235M grating ($med(\sigma_{G395M})$) should be the same for all
targets with the same exposure time. The same should be true for
$med(\sigma_{G395M})$ relative to $med(\sigma_{G235M})$. Thus, we scaled each
spectrum in the G140M grating such that the value of
$med(\sigma_{G140M})/med(\sigma_{G235M})$ for each target was equal to the
median value of $med(\sigma_{G140M})/med(\sigma_{G235M})$ across the whole
sample. The same procedure was used to scale the spectra in the G395M grating
relative to the G235M grating.

Subsequently, we produced synthetic photometric observations of all targets by
passing the spectra through six JWST filter curves (F090W, F115W, F200W,
F277W, F356W, and F444W). We then compared the synthetic photometry to the
observed photometry and calculated a scaling factor in each filter. For the
targets in which both the synthetic and the observed photometry were measured
with a signal-to-noise ratio $>5$, we applied the median value of the scaling
factors in all signicicantly detected filters to the spectra in all three
gratings. This was done to each individual target where there were available
photometric observations. Of the 117 spectra for which an initial redshift was
determined during extraction, 56 targets had no absolute flux scaling applied
due to lack of photometric observations or low signal-to-noise photometry.
However, for the 61 objects with significantly detected photometry, the median
absolute scaling factor is 1.405, and the standard deviation is 1.481 (0.17
dex), so we do not apply a multiplicative scaling factor to those 56 spectra.

###  2.3 SED and Emission-Line Fitting

The spectral energy distribution (SED) of each galaxy was fit using the full
suite of available photometry from JWST/NIRCam, HST/WFC3, Spitzer/IRAC, and
ground-based surveys to cover the rest-frame UV to the rest-frame infrared.

For both the CEERS and JADES samples, we used publicly available multi-
wavelength photometric catalogs constructed by G.
Brammer555https://s3.amazonaws.com/grizli-v2/JwstMosaics/v7/index.html. The
CEERS catalog includes 7 HST bands (F435W, F606W, F814W, F105W, F125W, F140W,
and F160W), and 7 JWST/NIRCam bands (F115W, F150W, F200W, F277W, F356W, F410M,
and F444W), while the JADES catalog includes 9 HST bands (F435W, F606W, F775W
F814W, F850LP, F105W, F125W, F140W, and F160W), and 14 JWST/NIRCam bands (both
the broadband filters F090W, F115W, F150W, F200W, F277W, F356W, and F444W, and
the medium-band filters F182M, F210M, F335M, F410M, F430M, F460M, and F480M).
For objects lacking NIRCam coverage, we used the SEDs cataloged by the 3D-HST
team (Skelton et al., 2014; Momcheva et al., 2016) in the AEGIS and GOODS-S
fields for CEERS and JADES, respectively. These SEDs include ground-based and
HST optical and near-IR photometry, and Spitzer/IRAC $3.6-8.0$ $\mu$m
measurements.

The SED fitting was done using the FAST algorithm described in Kriek et al.
(2015). Following the redshift and stellar-mass-dependent criteria described
in Section 2.2 of Du et al. (2018), we assume either a 1.4-times solar
metallicity and a Calzetti et al. (2000) dust extinction law (hereafter
$1.4Z_{\odot}$+Calz), or a sub-solar metallicity and a Small Magellanic Cloud
(SMC; Gordon et al., 2003) extinction curve (hereafter 0.28$Z_{\odot}$+SMC).
For galaxies in the range $1.4<z\leq 2.7$, we assumed 0.28$Z_{\odot}$+SMC for
galaxies below $10^{10.45}\ M_{\odot}$, and $1.4Z_{\odot}$+Calz otherwise. For
galaxies in the range $2.7<z\leq 3.4$, we assumed 0.28$Z_{\odot}$+SMC for
galaxies below $10^{10.66}\ M_{\odot}$, and $1.4Z_{\odot}$+Calz otherwise. For
all galaxies above $z=3.4$, we assumed 0.28$Z_{\odot}$+SMC. The best-fit SED
was then used to model the continuum for each galaxy during the initial line
fitting procedure. Using the line fluxes calculated during this initial line
fitting procedure, we corrected the photometry for contributions from emission
lines that were detected with $>$5$\sigma$ significance. The SEDs were then
re-fit using FAST on the emission-line-corrected photometry. The resulting
best-fit SEDs were used to model the continuum during the final run of the
line fitting.

The emission lines were fit using a custom Python procedure that utilizes
scipy.optimize.curve_fit() to perform a chi-squared minimization with a
Levenberg-Marquardt algorithm. Each emission line was modeled individually
with a single Gaussian profile; however, adjacent emission lines (e.g.,
H$\gamma$ and [O iii]$\lambda 4364$, or H$\alpha$ and [N ii]$\lambda\lambda
6550,6585$) were fit simultaneously. The flux ratio of the [N
ii]$\lambda\lambda 6550,6585$ line doublet was fixed with a ratio of 1:3
during the fitting procedure since the intrinsic line strengths of the doublet
members are fixed to the ratio of their Einstein A coefficients. The 68%
confidence interval (and associated 1$\sigma$ uncertainties) of the line
fluxes were obtained via Monte Carlo simulations, perturbing the flux density
spectrum according to the error spectrum 500 times and calculating the 16th
and 84th percentiles of the resulting flux distribution.

Figure 1: Redshift distribution of the combined CEERS and JADES sample that we
analyze consisting of 146 spectroscopically confirmed star-forming galaxies
with robust SED measurements. Of these, 106 galaxies come from CEERS and 40
galaxies come from JADES.

###  2.4 Sample Properties and SFR calculations

Of the 117 JADES galaxies with measured redshifts, in addition to ground-based
photometry, 110 have JWST/NIRCam and 3DHST photometry, 3 have only 3DHST
photometry, and 4 have no space-based photometric observations. Of the 231
galaxies with measured redshifts in the CEERS sample, 217 had a combination of
NIRCam and/or 3DHST measurements. We further restricted our sample to galaxies
with $>$5$\sigma$ detections of at least two Balmer lines for dust correction,
and we selected for star-forming galaxies by imposing a cutoff in specific SFR
of $\log{\rm(sSFR(SED)/yr^{-1})}>-11$, where sSFR(SED) is the sSFR calculated
from the SED fitting routine. Finally, we flagged and removed AGN-dominated
galaxies by identifying either an [N ii]/H$\alpha$ flux ratio $>$0.5, or
through the presence of broad Balmer emission. Overall, the combined CEERS and
JADES sample that we analyze consists of 146 galaxies in the redshift range
$1.4\leq z<7$.

Figure 2: H$\alpha$-based SFRs vs. stellar masses for the combined JADES and
CEERS sample binned by redshift. The darker and lighter shaded regions denote
the 1$\sigma_{\rm int}$ and 2$\sigma_{\rm int}$ intrinsic scatter intervals
about the SFMS, respectively. The literature curves from Shivaei et al.
(2015), Speagle et al. (2014), Topping et al. (2021), and Popesso et al.
(2023) were shifted down in SFR by 0.32, 0.32, 0.37, and 0.32 dex,
respectively to match our low-metallicity H$\alpha$ to SFR conversion. The
slightly different shift for Topping et al. (2021) reflects the use of the Hao
et al. (2011) conversion, whereas the other works adopt the conversion from
Kennicutt (1998). Figure 3: SED-based SFRs vs. stellar masses for the combined
JADES and CEERS sample binned by redshift. The darker and lighter shaded
regions denote the 1$\sigma_{\rm int}$ and 2$\sigma_{\rm int}$ intrinsic
scatter intervals about the SFMS, respectively. The lower normalization of our
sample compared to the literature is largely due to the fact that we assume a
SMC dust law, whereas most literature studies have assumed a Calzetti dust
law. Figure 4: UV-based SFRs vs. stellar masses for the combined JADES and
CEERS sample binned by redshift. The darker and lighter shaded regions denote
the 1$\sigma_{\rm int}$ and 2$\sigma_{\rm int}$ intrinsic scatter intervals
about the SFMS, respectively. The lower normalization of our sample compared
to the literature is largely due to the fact that we assume a SMC dust law,
whereas most literature studies have assumed a Calzetti dust law.

We show the redshift distribution of the combined JADES and CEERS NIRSpec
sample in Figure 1. The full sample extends between $z\sim 1-7$, providing
insights into the high-redshift universe enabled by the near-IR sensitivity
that JWST/NIRSpec provides. Prior to JWST, the properties and evolution of
star-forming galaxies from the local universe out to $z\sim 3$ had been
investigated extensively in the literature. Our sample, however, contains 94
galaxies at $z>2.7$ with $>$5$\sigma$ Balmer line detections, Balmer-
decrement-based E($B-V$) measurements, and SED information, enabling an
analysis of the properties of star-forming galaxies both during and after the
epoch of reionization. We calculated the SFRs of the galaxies in our sample
using the dust-corrected H$\alpha$ luminosity, the SED-based SFR from FAST,
and the dust-corrected UV luminosity. To correct the observed H$\alpha$
luminosity for dust attenuation, we assumed a Cardelli et al. (1989) dust law,
with E($B-V$) values derived from the H$\alpha$/H$\beta$ Balmer decrement. We
assumed an intrinsic ratio of H$\alpha$/H$\beta$ of 2.79, corresponding to
Case B recombination where $T_{e}=15,000\ \rm K$. In cases where H$\alpha$ was
not detected due to lack of wavelength coverage (15 galaxies), we used the
H$\gamma$/H$\beta$ ratio to dust-correct the emission lines, assuming an
intrinsic ratio of 0.475. Then, using the dust-corrected H$\beta$ flux, we
inferred the H$\alpha$ flux with the assumed intrinsic H$\alpha$/H$\beta$
ratio of 2.79. In cases where H$\beta$ was not detected due to lack of
wavelength coverage (13 galaxies), we used the H$\gamma$/H$\alpha$ ratio to
calculate E($B-V$), assuming an intrinsic ratio of 0.173. All Balmer line
fluxes have been measured accounting for Balmer absorption in the underlying
stellar continuum during the fitting procedure.

After correcting for dust attenuation, we converted the H$\alpha$ luminosities
to SFRs according to the following equation:

$\log\left({\rm{\frac{SFR}{M_{\odot}\
yr^{-1}}}}\right)=\log\left({\rm\frac{L_{H\alpha}}{erg\ s^{-1}}}\right)+C$ (1)

where $C$ is a metallicity-dependent conversion calculated from a set of BPASS
models (Stanway & Eldridge, 2018) with an upper stellar mass limit of 100 M⊙
(Reddy et al., 2022). For galaxies fit with the $1.4Z_{\odot}$+Calz
assumption, we used a conversion of $C=-41.37$, corresponding to a stellar
metallicity of $Z_{*}=0.02$. For galaxies fit with the 0.28$Z_{\odot}$+SMC
assumption, we used a conversion of $C=-41.59$, corresponding to a stellar
metallicity of $Z_{*}=0.004$. To calculate a final uncertainty on the measured
H$\alpha$-based SFR, we add the corresponding absolute flux calibration
uncertainty in quadrature with the measurement uncertainties, adding 0.17 dex
for those galaxies in JADES, and 0.23 dex for those in CEERS.

In order to correct the UV luminosity for dust attenuation, we calculated the
UV slope, $\beta$, by performing a linear regression to the observed rest-UV
photometry for each galaxy, which yields a median beta value of $-1.55$ with a
median uncertainty on beta of 0.36. We converted $\beta$ to the attenuation at
1600 Å, $A_{1600}$, using the following metallicity- and dust-law-dependent
conversion from Reddy et al. (2018):

$A_{1600}=\begin{cases}1.82\beta+4.43,&{\rm 1.4Z_{\odot}}\text{+Calz}\\\
0.93\beta+2.52,&0.28{\rm Z_{\odot}}\text{+SMC}\end{cases}$ (2)

We then converted the dust-corrected UV luminosity to a SFR using the Murphy
et al. (2011) and Hao et al. (2011) conversion:

$\log\left({\rm\frac{SFR}{M_{\odot}\ yr^{-1}}}\right)=\log\left({\rm\frac{\nu
L_{\nu,1600}}{erg\ s^{-1}}}\right)-43.35$ (3)

adjusted for a Chabrier (2003) IMF by adding $\log_{10}(0.77)$ to equation 3.
We plot the SFRs calculated using H$\alpha$, SED-fitting, and UV luminosity
vs. the stellar masses in Figures 2, 3, and 4, respectively.

In our analysis of the SFMS, we divide the galaxy sample into bins of redshift
in order to track changes in galaxy properties with cosmic time. We restrict
our analysis to objects with $z\geq 1.4$ for the purpose of comparison to
previous related work from the MOSDEF survey (e.g. Shivaei et al., 2015;
Topping et al., 2021). We define our lowest redshift bin, $1.4<z\leq 2.7$, to
capture galaxies up to $z=2.7$, beyond which the H$\alpha$ line is shifted out
of the wavelength range of ground-based observatories. The following redshift
bins ($2.7<z\leq 4$, $4<z\leq 5$, $5<z\leq 6$, and $6<z\leq 7$) represent the
new rest-optical redshift frontier enabled by JWST/NIRSpec.

##  3 Fitting The Star-forming Main Sequence

We performed a linear fit to the SFMS in Figures 2, 3, and 4 assuming an
equation of the form:

$\displaystyle\log({\rm SFR})$
$\displaystyle=m\times\log\left({\rm\frac{M_{*}}{10^{9.16}M_{\odot}}}\right)$
$\displaystyle+b+\mathcal{N}(0,\rm\sigma_{int}^{2})$ (4)

where $m$ is the slope, $b$ is the intercept, $10^{9.16}$ is the median
stellar mass of the sample, and $\rm\sigma_{int}$ is the intrinsic scatter
distributed normally ($\mathcal{N}$) about the best-fit relation. We
normalized the stellar masses by $10^{9.16}$, the median stellar mass, during
the fitting to minimize the covariance between the best-fit slope and the
intercept. Previous works have found evidence for a turnover in the SFMS at
high masses, leading some authors to use a different parametrization (e.g.,
Lee et al., 2015; Popesso et al., 2023). However, our galaxy sample
predominantly lies below the characteristic mass of $\sim$$10^{10.5}\
M_{\odot}$ where this turnover occurs, so we proceed using a linear fit.

We estimated the parameters $m$, $b$, and $\rm\sigma_{int}$ using the
following likelihood function $\mathcal{L}$ presented by Hogg et al. (2010):

$\displaystyle\ln{\mathcal{L}}$
$\displaystyle=-\sum_{i=1}^{N}\frac{1}{2}\ln[\Sigma_{i}^{2}+V]-\sum_{i=1}^{N}\frac{\Delta_{i}^{2}}{2[\Sigma_{i}^{2}+V]}$
(5) $\displaystyle\Delta_{i}$
$\displaystyle=\boldsymbol{\hat{\upsilon}}^{T}\begin{bmatrix}x_{i}\\\
y_{i}\end{bmatrix}-\frac{b}{\sqrt{1+m^{2}}}$ (6)
$\displaystyle\boldsymbol{\hat{\upsilon}}$
$\displaystyle=\frac{1}{\sqrt{1+m^{2}}}\begin{bmatrix}-m\\\ 1\end{bmatrix}$
(7)

where $x_{i}=\log{(\rm M_{*}/10^{9.16}\ M_{\odot})}_{i}$,
$y_{i}=\log{\rm(SFR/M_{\odot}\ yr^{-1})}_{i}$, $\Delta_{i}$ is the orthogonal
distance of each point $(x_{i},y_{i})$ to the best-fit line, $N$ is the number
of data points, $\boldsymbol{\hat{\upsilon}}$ is the unit vector orthogonal to
the best-fit relation, $\Sigma_{i}^{2}$ is the covariance matrix of each
measurement projected onto the best-fit line, defined as:

$\displaystyle\Sigma_{i}^{2}=\boldsymbol{\hat{\upsilon}}^{T}\begin{bmatrix}\sigma_{x,i}^{2}&\sigma_{xy,i}\\\
\sigma_{xy,i}&\sigma_{y,i}^{2}\end{bmatrix}\boldsymbol{\hat{\upsilon}}$ (8)

with $\sigma_{x,i}$ being the uncertainty on $x_{i}$, $\sigma_{y,i}$ being the
uncertainty on $y_{i}$, $\sigma_{xy,i}$ being the covariance between the two
measurements for a given galaxy, and $V$ is a matrix encapsulating the
intrinsic scatter, defined as:

$V=\boldsymbol{\hat{\upsilon}}^{T}\begin{bmatrix}0&0\\\
0&\rm\sigma_{int}^{2}\end{bmatrix}\boldsymbol{\hat{\upsilon}}$ (9)

Throughout this analysis, we assume that $\sigma_{xy,i}=0$ for all galaxies.
We emphasize that this is likely not the case when the SFR is estimated using
SED fitting, since these two quantities are not derived independently.
However, we do not quantify this covariance, and instead we proceed with this
caveat in mind when analyzing the SED-based SFMS.

We multiply our likelihood function $\mathcal{L}$ by a uniform prior
probability distribution, restricting each parameter to the following ranges:
$m=[0.5,1.5]$, $b=[-5,5]$, and $\sigma_{\rm int}=[0,1.5]$. To estimate the
best-fit parameters, we use the Python Markov Chain Monte Carlo (MCMC)
implementation emcee666https://emcee.readthedocs.io/en/stable/ developed by
Foreman-Mackey et al. (2013). We use 48 walkers and a maximum of 5000 steps to
sample the likelihood function, always discarding the first 100 “burn-in”
samples before evaluating the best-fit parameters. These hyperparameters
ensure that the number of steps is always greater than 50 times the
autocorrelation time. We evaluate the likelihood function within the bounds
defined by the priors, and we calculate the best-fit parameters for the SFMS
in each redshift bin for each SFR indicator. For each parameter, we report the
median value of the MCMC sample distribution, and we calculate the $1\sigma$
confidence interval by reporting the 16th and the 84th percentile values. The
results of this analysis are given in Table 1.

For each SFR indicator, we find a consistent level of intrinsic scatter at all
redshifts, with the inverse-variance-weighted average values being
$\langle\sigma_{{\rm int,H}\alpha}\rangle=0.48\pm 0.05$, $\langle\sigma_{{\rm
int,SED}}\rangle=0.42\pm 0.03$, and $\langle\sigma_{{\rm
int,UV}}\rangle=0.26\pm 0.02$. In Figure 5, we illustrate how our
$\rm\sigma_{int}$ measurements vary over each redshift bin for each SFR
indicator. We see that over the majority of the redshift range that we probe,
$\sigma_{\rm int,H\alpha}$ and $\sigma_{\rm int,SED}$ are greater than
$\sigma_{\rm int,UV}$.

The scatter in the SED-based SFMS is more consistent with the H$\alpha$-based
SFMS than with the UV-based SFMS, contrary to our expectations given that the
SED-based and UV-based SFRs are both measured using stellar continuum
emission. However, the SED-based SFMS scatter may be overestimated due to
anticorrelated errors between the stellar mass and SFR in the SED fitting.
Since the SFR and stellar mass are both sensitive to the normalization of the
SED model to the photometry, the two are highly correlated. If there is
positive covariance between the stellar mass and the SFR uncertainties from
SED fitting, then setting $\rm\sigma_{xy,i}=0$ and minimizing Equation 5 will
underestimate the intrinsic scatter. However, some studies (e.g., Salmon et
al., 2015; Kurczynski et al., 2016) have found in the analysis of their SED
fitting approach that the SFR vs. stellar mass covariance was negative, driven
largely by the uncertainties on the UV extinction. This negative covariance
would lead to an overestimate of the intrinsic scatter in the SED-based SFMS
using our methodology. Since we do not explicitly test for the covariance
between the stellar mass and the SED-based SFR, it is difficult to assess the
accuracy of our $\rm\sigma_{int,SED}$ calculations. However, since the
measurements of SFR(H$\alpha$) and SFR(UV) are made independently of the SED
fitting process, we proceed assuming that the uncertainties in the H$\alpha$\-
and UV-based SFRs are not correlated, or are only weakly correlated with the
stellar mass uncertainties.

In each redshift bin, we additionally calculate the intrinsic scatter below
the median stellar mass of the sample ($\rm M_{*}<10^{9.16}\ M_{\odot}$) and
above the median stellar mass of the sample ($\rm M_{*}>10^{9.16}\ M_{\odot}$)
to investigate whether we detect the mass-dependence of $\rm\sigma_{int}$
reported by other studies (e.g., Cole et al., 2023; Santini et al., 2017). To
calculate the mass-dependent $\rm\sigma_{int}$, we calculate the maximum
likelihood value of Equation 5, treating $V$ as a piecewise function:

$V=\begin{cases}V({\rm\sigma_{int,low}}),&{\rm M_{*}<10^{9.16}\ M_{\odot}}\\\
V({\rm\sigma_{int,high}}),&{\rm M_{*}>10^{9.16}\ M_{\odot}}\end{cases}$ (10)

We then performed 500 Monte Carlo simulations, perturbing each data point
$(x_{i},y_{i})$ by the measurement uncertainty $(\sigma_{x,i},\sigma_{y,i})$
and allowing galaxies to move between mass bins. The results of our
simulations are also tabulated in Table 1. We do not find strong evidence for
a mass dependence of the intrinsic scatter in the SFMS based on this analysis.

Figure 5: The intrinsic scatter in the SFMS $\rm\sigma_{int}$ as a function of
redshift. The points have an arbitrary offset in the x-axis to more easily see
overlapping points. The blue, orange, and green points show the intrinsic
scatter in the SFMS as determined from H$\alpha$, SED fitting, and UV
luminosity, respectively. Table 1: Star-forming Main Sequence best-fit
parameters as a function of redshift and SFR indicator.

SFR Indicator | Redshfit Bin | m | b | $\sigma_{\rm int}$ |  $\sigma_{\rm int,low}$aaThe intrinsic scatter in the SFMS as measured for galaxies with stellar masses less than the sample median of $10^{9.16}\ M_{\odot}$. |  $\sigma_{\rm int,high}$bbThe intrinsic scatter in the SFMS as measured for galaxies with stellar masses greater than the sample median of $10^{9.16}\ M_{\odot}$. | N
---|---|---|---|---|---|---|---
H$\alpha$ | $1.4<z\leq 7.0$ | $0.69^{+0.07}_{-0.08}$ | $0.56^{+0.05}_{-0.05}$ | $0.53^{+0.05}_{-0.04}$ | $0.57^{+0.05}_{-0.04}$ | $0.59^{+0.04}_{-0.04}$ | 146
H$\alpha$ | $1.4<z\leq 2.7$ | $0.95^{+0.15}_{-0.18}$ | $0.16^{+0.09}_{-0.09}$ | $0.56^{+0.10}_{-0.08}$ | $0.89^{+0.33}_{-0.28}$ | $0.58^{+0.06}_{-0.06}$ | 52
H$\alpha$ | $2.7<z\leq 4.0$ | $0.67^{+0.10}_{-0.12}$ | $0.59^{+0.09}_{-0.08}$ | $0.42^{+0.09}_{-0.07}$ | $0.45^{+0.07}_{-0.08}$ | $0.44^{+0.07}_{-0.07}$ | 32
H$\alpha$ | $4.0<z\leq 5.0$ | $0.74^{+0.15}_{-0.20}$ | $0.81^{+0.12}_{-0.12}$ | $0.54^{+0.15}_{-0.10}$ | $0.26^{+0.09}_{-0.10}$ | $0.74^{+0.11}_{-0.09}$ | 25
H$\alpha$ | $5.0<z\leq 6.0$ | $0.75^{+0.13}_{-0.16}$ | $0.83^{+0.12}_{-0.11}$ | $0.44^{+0.13}_{-0.10}$ | $0.51^{+0.08}_{-0.08}$ | $0.32^{+0.12}_{-0.13}$ | 23
H$\alpha$ | $6.0<z\leq 7.0$ | $0.69^{+0.13}_{-0.26}$ | $0.97^{+0.20}_{-0.15}$ | $0.42^{+0.21}_{-0.14}$ | $0.26^{+0.10}_{-0.11}$ | $0.62^{+0.21}_{-0.21}$ | 14
SED | $1.4<z\leq 7.0$ | $0.53^{+0.02}_{-0.04}$ | $0.51^{+0.04}_{-0.04}$ | $0.45^{+0.03}_{-0.03}$ | $0.44^{+0.02}_{-0.02}$ | $0.47^{+0.02}_{-0.02}$ | 146
SED | $1.4<z\leq 2.7$ | $0.62^{+0.07}_{-0.09}$ | $0.23^{+0.06}_{-0.06}$ | $0.39^{+0.05}_{-0.04}$ | $0.34^{+0.03}_{-0.03}$ | $0.39^{+0.02}_{-0.01}$ | 52
SED | $2.7<z\leq 4.0$ | $0.72^{+0.11}_{-0.14}$ | $0.61^{+0.09}_{-0.09}$ | $0.48^{+0.09}_{-0.07}$ | $0.29^{+0.03}_{-0.03}$ | $0.54^{+0.04}_{-0.04}$ | 32
SED | $4.0<z\leq 5.0$ | $0.55^{+0.04}_{-0.08}$ | $0.59^{+0.10}_{-0.09}$ | $0.44^{+0.10}_{-0.07}$ | $0.27^{+0.03}_{-0.03}$ | $0.55^{+0.07}_{-0.05}$ | 25
SED | $5.0<z\leq 6.0$ | $0.55^{+0.04}_{-0.09}$ | $0.81^{+0.10}_{-0.10}$ | $0.44^{+0.10}_{-0.07}$ | $0.40^{+0.04}_{-0.03}$ | $0.41^{+0.06}_{-0.05}$ | 23
SED | $6.0<z\leq 7.0$ | $0.63^{+0.09}_{-0.17}$ | $0.90^{+0.14}_{-0.12}$ | $0.36^{+0.15}_{-0.10}$ | $0.36^{+0.07}_{-0.06}$ | $0.08^{+0.11}_{-0.08}$ | 14
UV | $1.4<z\leq 7.0$ | $0.54^{+0.03}_{-0.04}$ | $0.58^{+0.03}_{-0.03}$ | $0.34^{+0.02}_{-0.02}$ | $0.35^{+0.02}_{-0.01}$ | $0.34^{+0.01}_{-0.01}$ | 145
UV | $1.4<z\leq 2.7$ | $0.59^{+0.05}_{-0.07}$ | $0.32^{+0.04}_{-0.05}$ | $0.31^{+0.04}_{-0.03}$ | $0.24^{+0.02}_{-0.02}$ | $0.33^{+0.01}_{-0.01}$ | 52
UV | $2.7<z\leq 4.0$ | $0.68^{+0.05}_{-0.06}$ | $0.63^{+0.05}_{-0.05}$ | $0.23^{+0.04}_{-0.03}$ | $0.20^{+0.02}_{-0.02}$ | $0.26^{+0.03}_{-0.03}$ | 31
UV | $4.0<z\leq 5.0$ | $0.56^{+0.04}_{-0.06}$ | $0.68^{+0.05}_{-0.05}$ | $0.23^{+0.05}_{-0.04}$ | $0.16^{+0.03}_{-0.02}$ | $0.30^{+0.04}_{-0.04}$ | 25
UV | $5.0<z\leq 6.0$ | $0.53^{+0.03}_{-0.05}$ | $0.88^{+0.08}_{-0.07}$ | $0.34^{+0.07}_{-0.06}$ | $0.32^{+0.06}_{-0.04}$ | $0.30^{+0.07}_{-0.07}$ | 23
UV | $6.0<z\leq 7.0$ | $0.64^{+0.08}_{-0.09}$ | $0.88^{+0.08}_{-0.07}$ | $0.18^{+0.08}_{-0.06}$ | $0.22^{+0.05}_{-0.06}$ | $0.08^{+0.04}_{-0.08}$ | 14

**footnotetext: We are unable to measure a UV-based SFR for one galaxy in the
$2.7<z\leq 4.0$ bin due to a lack of rest-UV photometric observations.

##  4 Discussion

###  4.1 SFR Intrinsic Scatter

We analyze both H$\alpha$\- and UV-based SFRs in order to characterize the
burstiness of galaxy SFHs in the early universe. We see in Figure 5 that the
H$\alpha$-based SFMS is characterized by a larger intrinsic scatter than the
UV-based SFMS, with an average difference of $\langle\sigma_{{\rm
int,H}\alpha}\rangle-\langle\sigma_{{\rm int,UV}}\rangle=0.21\pm 0.05$ dex. We
argue that this difference in $\rm\sigma_{int}$ is consistent with time
fluctuations in galaxy SFHs, though systematic uncertainties may also
contribute (see Section 4.3).

####  4.1.1 Comparison to theoretical works

Some previous theoretical works have predicted that galaxies should have
bursty SFHs at high redshift. Domínguez et al. (2015) analyzed hydrodynamic
simulations (gasoline; Wadsley et al., 2004) of 11 star-forming galaxies and
calculated the standard deviation in the H$\alpha$ luminosity values in 5 Myr
time steps from $z=2.74$ to $z=1.97$. They found that the standard deviation
in the H$\alpha$ luminosities ($L_{\rm H\alpha}$) of galaxies between $10^{8}$
and $10^{9}\ M_{\odot}$ was 0.45 dex, while that for galaxies with stellar
masses $\gtrsim$$10^{9}\ M_{\odot}$ was 0.22 dex. We compare these predictions
to our observations of the H$\alpha$-based SFMS in the $1.4<z\leq 2.7$ bin
assuming that our measurements of $\sigma_{\rm int,H\alpha}$ (and thus,
$L_{\rm H\alpha}$) capture the average deviations in SFR from a shared SFH.
The scatter in $L_{\rm H\alpha}$ that Domínguez et al. (2015) found in their
$10^{8}\text{-}10^{9}\ M_{\odot}$ galaxies (0.45 dex) is lower than our
SFR(H$\alpha$) intrinsic scatter of $0.56^{+0.10}_{-0.08}$ dex measured across
the full mass range of our sample. This discrepancy also holds with respect to
$\sigma_{\rm int,Ha}$ in our lower-mass bin ($<$$10^{9.16}\ M_{\odot}$), which
we determine to be $0.90^{+0.32}_{-0.28}$ dex. However, the uncertainty on
$\sigma_{\rm int,H\alpha}$ in the lower mass bin has large uncertainties owing
to the relatively few galaxies in this mass range and the measurement
uncertainties moving galaxies between mass bins. With respect to the higher-
mass ($\gtrsim$$10^{9}\ M_{\odot}$) sample of Domínguez et al. (2015), their
predicted standard deviation in $L_{H\alpha}$ of 0.22 dex is lower than what
we find for $\sigma_{\rm int,H\alpha}$ across the whole mass range in our
sample as well as in our high stellar mass bin ($0.56^{+0.10}_{-0.08}$ and
$0.58^{+0.06}_{-0.06}$, respectively).

When considering the UV-based SFMS, the standard deviation in $L_{UV}$ that
Domínguez et al. (2015) predicted also presents some differences when compared
with our measurements. The predicted SFR(UV) scatter from Domínguez et al. was
0.28 dex at $M_{*}=10^{8}\text{--}10^{9}\ M_{\odot}$, and 0.16 at
$M_{*}\gtrsim$10${}^{9}\ M_{\odot}$. The scatter in their
$10^{8}\text{--}10^{9}\ M_{\odot}$ sample is consistent with our measurement
of $\sigma_{\rm int,UV}=0.31^{+0.04}_{-0.03}$ across our full mass range
within the 1$\sigma$ confidence interval. However, the scatter in our low-mass
bin of $0.24^{+0.02}_{-0.02}$ dex is slightly smaller than the 0.28 dex
prediction. Whether this difference comes down to small sample numbers in our
low-mass bin or from differences in galaxy SFHs is still uncertain, and more
observations of galaxies at lower masses will improve this scatter
measurement. In their $\gtrsim$10${}^{9}\ M_{\odot}$ galaxies, the predicted
scatter of 0.16 is smaller than what we find in both our high-mass galaxies
($0.33^{+0.01}_{-0.01}$ dex) and across the whole mass range.

Our finding that the scatter in $L_{\rm H\alpha}$ predicted by Domínguez et
al. (2015) is smaller than our measured $\sigma_{\rm int,H\alpha}$ suggests
that systematic measurement uncertainties may be contributing to our inferred
scatter, or that the output from these simulations underpredicts the
stochasticity of SFRs in galaxies between $z=2.74$ and $z=1.97$. A larger,
mass-complete spectroscopic sample would help to distinguish between these
interpretations.

We also compare our data to the intrinsic scatter values measured in the FIRE
and IllustrisTNG simulations. In the FIRE simulations, Sparre et al. (2017)
reports the predicted scatter in the SFMS to be between 0.4–0.5 dex for
H$\alpha$-based SFRs and 0.3–0.4 dex for UV-based SFRs. For IllustrisTNG,
Donnari et al. (2019) reports that in the stellar mass range
$10^{9}\text{--}10^{10.5}\ M_{\odot}$, their 100 cMpc simulation predicts a
scatter of 0.2 dex for a 200-Myr-timescale SFR indicator (which we assume
closely matches the UV luminosity) and 0.25 dex for an instantaneous SFR
indicator (which we assume closely matches the H$\alpha$ luminosity). In the
case of FIRE, our measured $\sigma_{\rm int,H\alpha}$ and $\sigma_{\rm
int,UV}$ across the full mass range at $1.4<z\leq 2.7$ agree with the
predictions from FIRE at $z=2$. However, our measured scatter in both the
H$\alpha$-based and UV-based SFMS exceeds the predictions made by the
IllustrisTNG simulations at $z=2$.

The disagreement of our observations with IllustrisTNG may come down to the
simulation’s treatment of star formation and feedback that are on scales too
small to resolve. These processes are thought to contribute to the
stochasticity of star formation (e.g., Hopkins et al., 2018), and the
inability to resolve these processes may affect predictions regarding the
burstiness of SFHs. Considering this, our observation that $\sigma_{\rm
int,H\alpha}>\sigma_{\rm int,UV}$ out to $z\sim 7$ is consistent with the
interpretation that SFHs are bursty. The caveat to this, of course, is that
the systematic uncertainties associated with converting H$\alpha$ flux to SFR
are sub-dominant with respect to fluctuations in the galaxy SFHs. We discuss
this caveat further in 4.3.

Whether the scatter in the SFMS varies as a function of stellar mass at high
redshift does not have a clear consensus among theoretical works. Domínguez et
al. (2015) found that the scatter in H$\alpha$ luminosity varies as a function
of stellar mass, with a scatter of 1.0 dex at the low-mass end ($10^{7}\ \rm
M_{\odot}$) and 0.15 dex at the high-mass end ($10^{9.5}\ \rm M_{\odot}$).
Several other theoretical works also find a mass dependence in the SFMS
scatter (e.g., Ma et al., 2018; Ceverino et al., 2018; Katsianis et al., 2019;
Legrand et al., 2022), while some studies find no such trend (e.g., Donnari et
al., 2019; Pallottini & Ferrara, 2023). In our galaxy sample, we do not see
any consistent mass dependence on $\sigma_{\rm int}$ in the SFMS between the
low- and high-mass bins across redshift.

####  4.1.2 Comparison to observational works

In addition to theoretical predictions, previous observational works have
constrained $\sigma_{\rm int}$ in the SFMS. Studies that use SFR indicators
tracing star formation on 10-Myr timescales typically find intrinsic scatter
values of $\sigma_{\rm int}\sim 0.2\text{--}0.5$ dex (e.g., Kashino et al.,
2013; Shivaei et al., 2015; Cole et al., 2023). Many other studies find an
intrinsic scatter of 0.2–0.3 dex when analyzing SFR indicators that trace star
formation over longer timescales (e.g. Whitaker et al., 2012; Schreiber et
al., 2015; Tomczak et al., 2016).

Figure 6: The H$\alpha$-based star-forming main sequence as determined from
JWST observations in this study (purple) and MOSDEF observations from Sanders
et al. (2021) (gray). The MOSDEF data points have been shifted down by 0.32
dex to match the H$\alpha$ luminosity to SFR conversion that we use in this
paper.

In the case of our measurements of $\sigma_{\rm int}$ from both UV-based and
H$\alpha$-based SFRs, we find results that are consistent with previous
literature. In Figure 6, we compare our results to those found by Sanders et
al. (2021) in the MOSDEF survey, which uses an expanded version of the sample
from Shivaei et al. (2015). We note that our measured $\sigma_{\rm
int,H\alpha}$ exceeds that found from the MOSDEF survey in a similar redshift
range. However, our sample is more sensitive to lower H$\alpha$ luminosities
and lower stellar masses. Thus, our interpretation is that including galaxies
in this new parameter space increases the scatter when compared with previous
work (Shivaei et al., 2015).

As is the case in theoretical work, observational studies yield no consistent
answer regarding the mass-dependence of the scatter in the SFMS. Some studies
report an increase in scatter with increasing stellar mass (e.g., Guo et al.,
2013; Popesso et al., 2019a, b; Sherman et al., 2021; Cole et al., 2023), some
do not report any monotonic mass-dependence on the scatter (e.g., Noeske et
al., 2007; Schreiber et al., 2015; Shivaei et al., 2015; Davies et al., 2022),
and some report a decrease in scatter with increasing stellar mass (e.g.,
Willett et al., 2015).

When comparing our low-redshift sample with the MOSDEF (Kriek et al., 2015;
Shivaei et al., 2015) sample, however, we tentatively observe a mass
dependence. The median stellar mass of our JWST sample in the lowest redshift
bin is $10^{9.48}\ M_{\odot}$, while the median stellar mass of the MOSDEF
sample spanning $2.0<z<2.6$ is $10^{9.88}\ M_{\odot}$. The fact that these
samples span different ranges in stellar mass and exhibit different amounts of
scatter in the SFMS provide tentative evidence for mass-dependence of
$\sigma_{\rm int}$ when probing a large enough range in $M_{*}$. A larger,
mass-complete spectroscopic sample is needed, however, to make a definitive
statement about the presence or absence of mass dependence on the intrinsic
SFMS scatter.

We also note that we do not measure any strong redshift evolution in
$\sigma_{\rm int}$ for any of the SFR indicators, as illustrated in Figure 5.
There has been increased interest in determining whether an evolution in
bursty SFHs can explain the excess of UV-bright galaxies at $z>10$ (e.g., Sun
et al., 2023a, b; Shen et al., 2023; Mason et al., 2023). Though our sample
does not extend past $z\sim 7$, a larger, higher-redshift sample with smaller
measurement uncertainties on all of the SFR indicators would allow for a
better determination of the evolution of $\sigma_{\rm int}$ with cosmic time.

###  4.2 The H$\alpha$-to-UV luminosity ratio

Figure 7: The H$\alpha$-to-UV luminosity ratio vs. stellar mass for our galaxy sample. The gray hatched region shows the parameter space that we are not sensitive to assuming a limiting survey flux and a UV-luminosity to stellar mass ratio (which we derive from the parameters reported in Table 1). We take the limiting flux in our sample to be $2.0\times 10^{-19}\ \rm erg\ s^{-1}\ cm^{-2}$ since this is the lowest H$\alpha$ flux we observe at a significance greater than 5$\sigma$. The dashed gray lines denote the uncertainty on the region of our sensitivity limits due to the intrinsic scatter in SFR(UV) vs. stellar mass relation. The horizontal blue hatched region shows the equilibrium range of H$\alpha$-to-UV luminosity calculated by Mehta et al. (2023). Table 2: The fractions of galaxies in vs. out of equilibrium SFHs in each redshift bin. Redshift Bin | $f_{eq}^{\tablenotemark{\scriptsize a}}$ | $f_{above}^{\tablenotemark{\scriptsize b}}$ | $f_{below}^{\tablenotemark{\scriptsize c}}$ | $\sigma_{\rm H\alpha/UV}^{\tablenotemark{\scriptsize d}}$
---|---|---|---|---
$1.4<z\leq 2.7$ | 0.54 | 0.13 | 0.33 | $0.3^{+0.05}_{-0.05}$
$2.7<z\leq 4$ | 0.71 | 0.13 | 0.16 | $0.3^{+0.07}_{-0.06}$
$4<z\leq 5$ | 0.52 | 0.32 | 0.16 | $0.36^{+0.09}_{-0.08}$
$5<z\leq 6$ | 0.48 | 0.17 | 0.35 | $0.46^{+0.1}_{-0.09}$
$6<z\leq 7$ | 0.64 | 0.21 | 0.14 | $0.29^{+0.13}_{-0.1}$

Another common metric used to characterize the burstiness of SFHs is the
L(H$\alpha$)/$\nu\rm L_{\nu,1600}$ ratio. Since the H$\alpha$ and UV
luminosities respond to changes in the SFR on different timescales ($\sim$5
Myr and $\sim$100 Myr, respectively), measuring their ratio gives information
about the nature of SFHs for a galaxy population (see Flores Velázquez et al.,
2021; Rezaee et al., 2023, for more detailed discussions). Using stellar
population synthesis models and a constant SFH, Mehta et al. (2023) found that
the $\log({\rm L(H\alpha)/\nu\rm L_{\nu,1600}})$ ratio reaches an equilibrium
after $\sim$100-200 Myr, reaching values the range of $-1.93$ and $-1.78$ for
metallicities ($\log{(Z/Z_{\odot})}$) between $-2$ and 0. Using this
equilibrium range, Asada et al. (2024) analyzed the H$\alpha$-to-UV luminosity
ratio in their sample of $z\sim 4.7\text{--}6.5$ galaxies, finding that 60% of
their galaxies deviate from a constant SFH by greater than $1\sigma$.

We plot the L(H$\alpha$)/$\nu\rm L_{\nu,1600}$ ratio of our galaxy sample in
Figure 7, along with the equilibrium range calculated by Mehta et al. (2023)
indicated by a blue, hatched, horizontal band. We additionally shade in the
regions below which we do not expect to detect galaxies due to the sensitivity
limits of the NIRSpec spectra ($\sim$$2\times 10^{-19}\ \rm erg\ s^{-1}\
cm^{-2}$). We convert this limiting flux to a lower limit on the detectable
H$\alpha$ luminosity at the median redshift of the galaxy sample in each
redshift bin. We then calculate the predicted value of $\nu L_{\nu,1600}$ as a
function of stellar mass by translating our SFR(UV) vs. stellar mass linear
fit parameters in Table 1 to a relationship between $\nu L_{\nu,1600}$ and
stellar mass. This results in a curve that represents our predicted
sensitivity to $L_{\rm H\alpha}/\nu L_{\nu,1600}$ as a function of stellar
mass.

In Table 2, for each redshift bin, we show the fraction of galaxies whose
H$\alpha$-to-UV ratios are consistent with the equilibrium value to within
1$\sigma$ ($f_{eq}$) and the fraction of galaxies with H$\alpha$-to-UV ratios
above and below the equilibrium value ($f_{above}$ and $f_{below}$,
respectively). We also calculate the standard deviation in the H$\alpha$-to-UV
ratio, $\sigma_{\rm H\alpha/UV}$. Across redshift bins, we see that between
48%–71% of galaxies fall within the predicted equilibrium range of
H$\alpha$-to-UV ratios. Our measured $f_{eq}$ represents an upper limit
however, since it is not able to distinguish between galaxies that have been
steadily forming stars for 100-200 Myr and galaxies that are transitioning
between star-forming and quenching episodes. Our data reveal no monotonic
trend of $f_{eq}$ with redshift, and instead indicate that the galaxy
population is characterized by a mixture of smooth and bursty SFHs.

We also consider $f_{above}$ and $f_{below}$ in an attempt to constrain the
duty cycle of star formation in the early universe. We expect that a low duty
cycle of star formation will result in $f_{below}>f_{above}$ and vice versa.
Our data reveal no consistent trend with $f_{above}$ or $f_{below}$ vs.
redshift. In the lowest two redshift bins, $f_{below}>f_{above}$, suggesting a
low duty cycle of star formation in this redshift range. However,
$f_{above}>f_{below}$ in the $4<z\leq 5$ bin, which would suggest a high duty
cycle. Endsley et al. (2023b) suggested in their analysis of $z\sim 6$
galaxies that UV-faint (and thus low-mass) galaxies may be more likely to be
observed in recent downturns in their star formation. This would produce low
H$\alpha$ to UV ratios at low stellar masses in our sample in the high
redshift bins. Due to the relatively low numbers of galaxies in these highest
two redshift bins, however, it is difficult to ascertain how physically
meaningful $f_{above}$ and $f_{below}$ are. Additionally, the sensitivity
limits of the surveys begin to hinder our ability to measure $\log(\rm
L_{H\alpha}/\nu L_{\nu,1600})$ for low-mass galaxies at these redshifts. A
future deep spectroscopic analysis of $z\sim 6$ galaxies would be able to
provide a better view of the star-forming behavior of galaxies at this epoch.
Finally, we tabulate the measurement-subtracted scatter in the $\log(\rm
L_{H\alpha}/\nu L_{\nu,1600})$ distribution for each redshift bin. The scatter
ranges between 0.30–0.46 dex, suggesting that time-dependent fluctuations in
galaxy SFHs may result in differences in SFR that vary by a factor of
$\sim$2–3 when determined from H$\alpha$ vs. UV luminosity.

###  4.3 Important considerations

We have based our discussion thus far upon measurements of galaxy SFRs derived
from H$\alpha$ and UV luminosities. There are, however, systematic
uncertainties in these conversions that are difficult to quantify, and they
may not be negligible.

The conversion between H$\alpha$ luminosity to SFR, for example, introduces
some uncertainties. Deriving a SFR from hydrogen recombination line
luminosities assumes that the escape fraction of ionizing radiation from
galaxies ($f_{esc}$) is zero, which may not be the case in the early universe,
especially for faint galaxies which are thought to be key contributors to
cosmic Reionization (e.g., Endsley et al., 2023a). Direct measurements of
$f_{esc}$ at $z\sim 3$ by Pahl et al. (2021) suggest that galaxies typically
have escape fractions of 5–10%. Indirect studies of the SEDs of galaxies at
higher redshifts suggest a range of escape fractions, with typical values of
$\sim 13$% (Mascia et al., 2024), reaching up to 50% in the most extreme cases
(Topping et al., 2022; Endsley et al., 2023a). However, the highest inferred
escape fractions are characteristic of galaxies with the faintest UV
luminosities ($M_{UV}\sim-17.5$) and bluest UV slopes ($\beta\sim-3$),
properties not representative of the majority of our sample. Thus, it is
unlikely that uncertainties on $f_{esc}$ dominate the systematics in this
study.

The conversion factor that we choose to translate betwen H$\alpha$ luminosity
and SFR also introduces systematic uncertainties. We choose these conversion
factors self-consistently with the metallicity assumptions used to derive
stellar population properties from the SED fitting. However, calculating the
ionizing photon production based on the H$\alpha$ luminosity likely involves
more nuance than what is captured by choosing between two different
metallicity-dependent conversion factors. The difference between the two
conversion factors that we use is 0.22 dex, so the magnitude of the systematic
uncertainty on SFR(H$\alpha$) is likely on this order. In addition, the fact
that we use a single UV luminosity to SFR conversion factor rather than
accounting for metallicity dependence introduces some uncertainty.

There are also uncertainties involved in choosing which dust attenuation law
to apply to the stellar component of each galaxy in the sample. Several
studies analyzing the infrared excess vs. UV slope relation in high redshift
galaxies have determined that higher redshift, lower metallicity galaxies are
better described by an SMC attenuation law, while local, metal-rich galaxies
are better described by a Calzetti et al. (2000) dust law (e.g., Reddy et al.,
2018). However, though this trend may apply on average to galaxy populations,
the specific dust law that best applies to an individual galaxy is difficult
to determine, and assuming that a particular dust law applies to many galaxies
may lead to increased scatter in the inferred SFRs (e.g., Shivaei et al.,
2015, 2018). We also note the possibility of missing nebular emission coming
from optically-thick dust-enshrouded star-forming regions (Kennicutt & Evans,
2012). However, observations of Paschen lines by Reddy et al. (2023) in
galaxies at $z\sim 1\text{--}3$ suggest that this effect would boost the
SFR(H$\alpha$) by 25% at most. Additionally, we do not expect this effect to
be particularly pronounced since most of our sample lies at low stellar mass
($\log(\rm M_{*}/M_{\odot})<10$), where dust attenuation is less significant
(e.g., Whitaker et al., 2017; McLure et al., 2018).

Finally, variations in the stellar IMF may contribute significantly to derived
galaxy properties (e.g., Wang et al., 2024). Throughout this analysis, we
assume a Chabrier (2003) IMF, however, observations of galaxy properties in
the early universe have led some authors to entertain the possibility that
significant variations in the IMF exist in the early universe (e.g., Katz et
al., 2022; Pacucci et al., 2022). It is not feasible to determine the specific
form of the stellar IMF in each galaxy (see Hennebelle & Grudíc, 2024, for a
more recent review on the subject), but this may be a source of uncertainty in
our measurements if the IMF for early galaxies does vary significantly from
the Chabrier (2003) functional form.

##  5 Conclusions

In this study, we have analyzed the galaxy SFMS at $1.4<z<7$ using a combined
CEERS and JADES NIRSpec sample. We measure the SFRs using different
indicators: H$\alpha$, SED-fitting output, and UV luminosity. With these
various indicators, we are able to investigate the burstiness of galaxy SFHs
in the early universe. Our main conclusions are as follows:

  1. 1.

We find that the intrinsic scatter in the H$\alpha$-based SFMS is larger than
the UV-based SFMS by $0.21\pm 0.05$ dex on average. This difference is
consistent with the idea that galaxies with bursty SFHs will deviate from the
SFMS on short timescales, resulting in increased scatter for short timescale
indicators such as H$\alpha$ compared to longer timescale indicators such as
the FUV luminosity.

  2. 2.

We do not see evidence for mass-dependence in the intrinsic scatter about the
SFMS within our sample. We do, however, see tentative evidence for larger
scatter at lower masses when comparing our sample with the MOSDEF survey,
suggesting that a large dynamic range in stellar mass is required to probe any
mass dependence in $\sigma_{\rm int}$.

  3. 3.

We also find evidence for bursty SFHs in some galaxies when analyzing the
H$\alpha$-to-UV luminosity ratio. Depending on the redshift, 29–52% of
galaxies in our sample deviate from the ratio predicted for a stellar
population with a constant SFH over a 100 Myr timespan. The degree of scatter
about this equilibrium ratio also suggests the influence of burstiness on the
SFHs of the galaxies in our sample. In order to fully see how this ratio
behaves as a function of stellar mass, a deep, mass-complete spectroscopic
sample is required.

This study represents a preliminary analysis of the characteristics of the
SFHs of star-forming galaxies in the early universe with multiple SFR
indicators. Evidence for bursty star-formation histories has important
implications for several areas, namely in the determination of stellar
population properties through SED fitting and matching of theoretical
simulations of galaxy formation with observations. We anticipate the promise
of repeating this analysis with order-of-magnitude larger samples at similar
redshifts covered by multiple SFR indicators.

##  6 Acknowledgements

We acknowledge the CEERS and JADES teams for their effort to design, execute,
and make public their observational surveys. We would also like to thank Kate
Whitaker, Zoë Haggard, James Bullock, and his research group for useful
discussions. We also acknowledge support from NASA grant JWST-GO-01914. This
work is based on observations made with the NASA/ESA/CSA James Webb Space
Telescope. The data were obtained from the Mikulski Archive for Space
Telescopes at the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS5-03127 for JWST. Data were also obtained from the DAWN JWST
Archive maintained by the Cosmic Dawn Center. The specific observations
analyzed can be accessed via
https://archive.stsci.edu/doi/resolve/resolve.html?doi=10.17909/z7p0-8481
(catalog DOI: 10.17909/z7p0-8481) and
https://archive.stsci.edu/doi/resolve/resolve.html?doi=10.17909/8tdj-8n28
(catalog DOI: 10.17909/8tdj-8n28).

## References

  * Arrabal Haro et al. (2023) Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, ApJ, 951, L22, doi: 10.3847/2041-8213/acdd54
  * Asada et al. (2024) Asada, Y., Sawicki, M., Abraham, R., et al. 2024, MNRAS, 527, 11372, doi: 10.1093/mnras/stad3902
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Atek et al. (2022) Atek, H., Furtak, L. J., Oesch, P., et al. 2022, MNRAS, 511, 4464, doi: 10.1093/mnras/stac360
  * Brammer (2023) Brammer, G. 2023, msaexp: NIRSpec analyis tools, 0.6.17, Zenodo, doi: 10.5281/zenodo.7299500
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151, doi: 10.1111/j.1365-2966.2004.07881.x
  * Bunker et al. (2023) Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2306.02467, doi: 10.48550/arXiv.2306.02467
  * Bushouse et al. (2023) Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline, 1.9.4, Zenodo, doi: 10.5281/zenodo.7577320
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Ceverino et al. (2018) Ceverino, D., Klessen, R. S., & Glover, S. C. O. 2018, MNRAS, 480, 4842, doi: 10.1093/mnras/sty2124
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Ciesla et al. (2023) Ciesla, L., Elbaz, D., Ilbert, O., et al. 2023, arXiv e-prints, arXiv:2309.15720, doi: 10.48550/arXiv.2309.15720
  * Cole et al. (2023) Cole, J. W., Papovich, C., Finkelstein, S. L., et al. 2023, arXiv e-prints, arXiv:2312.10152, doi: 10.48550/arXiv.2312.10152
  * Curti et al. (2024) Curti, M., Maiolino, R., Curtis-Lake, E., et al. 2024, A&A, 684, A75, doi: 10.1051/0004-6361/202346698
  * Davies et al. (2022) Davies, L. J. M., Thorne, J. E., Bellstedt, S., et al. 2022, MNRAS, 509, 4392, doi: 10.1093/mnras/stab3145
  * Domínguez et al. (2015) Domínguez, A., Siana, B., Brooks, A. M., et al. 2015, MNRAS, 451, 839, doi: 10.1093/mnras/stv1001
  * Donnari et al. (2019) Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485, 4817, doi: 10.1093/mnras/stz712
  * Du et al. (2018) Du, X., Shapley, A. E., Reddy, N. A., et al. 2018, ApJ, 860, 75, doi: 10.3847/1538-4357/aabfcf
  * Eisenstein et al. (2023) Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv e-prints, arXiv:2306.02465, doi: 10.48550/arXiv.2306.02465
  * Emami et al. (2019) Emami, N., Siana, B., Weisz, D. R., et al. 2019, ApJ, 881, 71, doi: 10.3847/1538-4357/ab211a
  * Endsley et al. (2023a) Endsley, R., Stark, D. P., Whitler, L., et al. 2023a, MNRAS, 524, 2312, doi: 10.1093/mnras/stad1919
  * Endsley et al. (2023b) —. 2023b, arXiv e-prints, arXiv:2306.05295, doi: 10.48550/arXiv.2306.05295
  * Finkelstein et al. (2023) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13, doi: 10.3847/2041-8213/acade4
  * Flores Velázquez et al. (2021) Flores Velázquez, J. A., Gurvich, A. B., Faucher-Giguère, C.-A., et al. 2021, MNRAS, 501, 4812, doi: 10.1093/mnras/staa3893
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Gordon et al. (2003) Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279, doi: 10.1086/376774
  * Guo et al. (2013) Guo, K., Zheng, X. Z., & Fu, H. 2013, ApJ, 778, 23, doi: 10.1088/0004-637X/778/1/23
  * Hao et al. (2011) Hao, C.-N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124, doi: 10.1088/0004-637X/741/2/124
  * Heintz et al. (2023) Heintz, K. E., Watson, D., Brammer, G., et al. 2023, arXiv e-prints, arXiv:2306.00647, doi: 10.48550/arXiv.2306.00647
  * Heintz et al. (2024) Heintz, K. E., Brammer, G. B., Watson, D., et al. 2024, arXiv e-prints, arXiv:2404.02211, doi: 10.48550/arXiv.2404.02211
  * Hennebelle & Grudíc (2024) Hennebelle, P., & Grudíc, M. 2024, arXiv e-prints, arXiv:2404.07301, doi: 10.48550/arXiv.2404.07301
  * Hogg et al. (2010) Hogg, D. W., Bovy, J., & Lang, D. 2010, arXiv e-prints, arXiv:1008.4686, doi: 10.48550/arXiv.1008.4686
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Horne (1986) Horne, K. 1986, PASP, 98, 609, doi: 10.1086/131801
  * Jakobsen et al. (2022) Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80, doi: 10.1051/0004-6361/202142663
  * Kashino et al. (2013) Kashino, D., Silverman, J. D., Rodighiero, G., et al. 2013, ApJ, 777, L8, doi: 10.1088/2041-8205/777/1/L8
  * Katsianis et al. (2019) Katsianis, A., Zheng, X., Gonzalez, V., et al. 2019, ApJ, 879, 11, doi: 10.3847/1538-4357/ab1f8d
  * Katz et al. (2022) Katz, H., Rosdahl, J., Kimm, T., et al. 2022, MNRAS, 510, 5603, doi: 10.1093/mnras/stac028
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ARA&A, 36, 189, doi: 10.1146/annurev.astro.36.1.189
  * Kennicutt & Evans (2012) Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531, doi: 10.1146/annurev-astro-081811-125610
  * Kriek et al. (2015) Kriek, M., Shapley, A. E., Reddy, N. A., et al. 2015, ApJS, 218, 15, doi: 10.1088/0067-0049/218/2/15
  * Kurczynski et al. (2016) Kurczynski, P., Gawiser, E., Acquaviva, V., et al. 2016, ApJ, 820, L1, doi: 10.3847/2041-8205/820/1/L1
  * Lee et al. (2015) Lee, N., Sanders, D. B., Casey, C. M., et al. 2015, ApJ, 801, 80, doi: 10.1088/0004-637X/801/2/80
  * Legrand et al. (2022) Legrand, L., Hutter, A., Dayal, P., et al. 2022, MNRAS, 509, 595, doi: 10.1093/mnras/stab3034
  * Looser et al. (2023) Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2023, arXiv e-prints, arXiv:2306.02470, doi: 10.48550/arXiv.2306.02470
  * Ma et al. (2018) Ma, X., Hopkins, P. F., Garrison-Kimmel, S., et al. 2018, MNRAS, 478, 1694, doi: 10.1093/mnras/sty1024
  * Mascia et al. (2024) Mascia, S., Pentericci, L., Calabrò, A., et al. 2024, A&A, 685, A3, doi: 10.1051/0004-6361/202347884
  * Mason et al. (2023) Mason, C. A., Trenti, M., & Treu, T. 2023, MNRAS, 521, 497, doi: 10.1093/mnras/stad035
  * McLure et al. (2018) McLure, R. J., Dunlop, J. S., Cullen, F., et al. 2018, MNRAS, 476, 3991, doi: 10.1093/mnras/sty522
  * Mehta et al. (2023) Mehta, V., Teplitz, H. I., Scarlata, C., et al. 2023, ApJ, 952, 133, doi: 10.3847/1538-4357/acd9cf
  * Mirocha & Furlanetto (2023) Mirocha, J., & Furlanetto, S. R. 2023, MNRAS, 519, 843, doi: 10.1093/mnras/stac3578
  * Momcheva et al. (2016) Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27, doi: 10.3847/0067-0049/225/2/27
  * Murphy et al. (2011) Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67, doi: 10.1088/0004-637X/737/2/67
  * Narayanan et al. (2024) Narayanan, D., Lower, S., Torrey, P., et al. 2024, ApJ, 961, 73, doi: 10.3847/1538-4357/ad0966
  * Noeske et al. (2007) Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43, doi: 10.1086/517926
  * Pacucci et al. (2022) Pacucci, F., Dayal, P., Harikane, Y., Inoue, A. K., & Loeb, A. 2022, MNRAS, 514, L6, doi: 10.1093/mnrasl/slac035
  * Pahl et al. (2021) Pahl, A. J., Shapley, A., Steidel, C. C., Chen, Y., & Reddy, N. A. 2021, MNRAS, 505, 2447, doi: 10.1093/mnras/stab1374
  * Pallottini & Ferrara (2023) Pallottini, A., & Ferrara, A. 2023, A&A, 677, L4, doi: 10.1051/0004-6361/202347384
  * Popesso et al. (2019a) Popesso, P., Concas, A., Morselli, L., et al. 2019a, MNRAS, 483, 3213, doi: 10.1093/mnras/sty3210
  * Popesso et al. (2019b) Popesso, P., Morselli, L., Concas, A., et al. 2019b, MNRAS, 490, 5285, doi: 10.1093/mnras/stz2635
  * Popesso et al. (2023) Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526, doi: 10.1093/mnras/stac3214
  * Reddy et al. (2023) Reddy, N. A., Topping, M. W., Sanders, R. L., Shapley, A. E., & Brammer, G. 2023, ApJ, 948, 83, doi: 10.3847/1538-4357/acc869
  * Reddy et al. (2018) Reddy, N. A., Oesch, P. A., Bouwens, R. J., et al. 2018, ApJ, 853, 56, doi: 10.3847/1538-4357/aaa3e7
  * Reddy et al. (2022) Reddy, N. A., Topping, M. W., Shapley, A. E., et al. 2022, ApJ, 926, 31, doi: 10.3847/1538-4357/ac3b4c
  * Rezaee et al. (2023) Rezaee, S., Reddy, N. A., Topping, M. W., et al. 2023, MNRAS, 526, 1512, doi: 10.1093/mnras/stad2842
  * Rieke et al. (2023) Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023, PASP, 135, 028001, doi: 10.1088/1538-3873/acac53
  * Salmon et al. (2015) Salmon, B., Papovich, C., Finkelstein, S. L., et al. 2015, ApJ, 799, 183, doi: 10.1088/0004-637X/799/2/183
  * Sanders et al. (2023) Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2023, ApJ, 955, 54, doi: 10.3847/1538-4357/acedad
  * Sanders et al. (2021) Sanders, R. L., Shapley, A. E., Jones, T., et al. 2021, ApJ, 914, 19, doi: 10.3847/1538-4357/abf4c1
  * Santini et al. (2017) Santini, P., Fontana, A., Castellano, M., et al. 2017, ApJ, 847, 76, doi: 10.3847/1538-4357/aa8874
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521, doi: 10.1093/mnras/stu2058
  * Schreiber et al. (2015) Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74, doi: 10.1051/0004-6361/201425017
  * Shapley et al. (2005) Shapley, A. E., Steidel, C. C., Erb, D. K., et al. 2005, ApJ, 626, 698, doi: 10.1086/429990
  * Shen et al. (2023) Shen, X., Vogelsberger, M., Boylan-Kolchin, M., Tacchella, S., & Kannan, R. 2023, MNRAS, 525, 3254, doi: 10.1093/mnras/stad2508
  * Sherman et al. (2021) Sherman, S., Jogee, S., Florez, J., et al. 2021, MNRAS, 505, 947, doi: 10.1093/mnras/stab1350
  * Shivaei et al. (2015) Shivaei, I., Reddy, N. A., Shapley, A. E., et al. 2015, ApJ, 815, 98, doi: 10.1088/0004-637X/815/2/98
  * Shivaei et al. (2018) Shivaei, I., Reddy, N. A., Siana, B., et al. 2018, ApJ, 855, 42, doi: 10.3847/1538-4357/aaad62
  * Simmonds et al. (2024) Simmonds, C., Tacchella, S., Hainline, K., et al. 2024, MNRAS, 527, 6139, doi: 10.1093/mnras/stad3605
  * Skelton et al. (2014) Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24, doi: 10.1088/0067-0049/214/2/24
  * Sparre et al. (2017) Sparre, M., Hayward, C. C., Feldmann, R., et al. 2017, MNRAS, 466, 88, doi: 10.1093/mnras/stw3011
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15, doi: 10.1088/0067-0049/214/2/15
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676, doi: 10.1093/mnras/stx3304
  * Stanway & Eldridge (2018) Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 479, 75, doi: 10.1093/mnras/sty1353
  * Sun et al. (2023a) Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., & Shen, X. 2023a, MNRAS, 526, 2665, doi: 10.1093/mnras/stad2902
  * Sun et al. (2023b) Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., et al. 2023b, ApJ, 955, L35, doi: 10.3847/2041-8213/acf85a
  * Tomczak et al. (2016) Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2016, ApJ, 817, 118, doi: 10.3847/0004-637X/817/2/118
  * Topping et al. (2022) Topping, M. W., Stark, D. P., Endsley, R., et al. 2022, ApJ, 941, 153, doi: 10.3847/1538-4357/aca522
  * Topping et al. (2021) Topping, M. W., Shapley, A. E., Sanders, R. L., et al. 2021, MNRAS, 506, 1237, doi: 10.1093/mnras/stab1793
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wadsley et al. (2004) Wadsley, J. W., Stadel, J., & Quinn, T. 2004, New A, 9, 137, doi: 10.1016/j.newast.2003.08.004
  * Wang et al. (2024) Wang, B., Leja, J., Atek, H., et al. 2024, ApJ, 963, 74, doi: 10.3847/1538-4357/ad187c
  * Weisz et al. (2012) Weisz, D. R., Johnson, B. D., Johnson, L. C., et al. 2012, ApJ, 744, 44, doi: 10.1088/0004-637X/744/1/44
  * Whitaker et al. (2017) Whitaker, K. E., Pope, A., Cybulski, R., et al. 2017, ApJ, 850, 208, doi: 10.3847/1538-4357/aa94ce
  * Whitaker et al. (2012) Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29, doi: 10.1088/2041-8205/754/2/L29
  * Whitaker et al. (2014) Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104, doi: 10.1088/0004-637X/795/2/104
  * Whitler et al. (2023) Whitler, L., Stark, D. P., Endsley, R., et al. 2023, MNRAS, 519, 5859, doi: 10.1093/mnras/stad004
  * Willett et al. (2015) Willett, K. W., Schawinski, K., Simmons, B. D., et al. 2015, MNRAS, 449, 820, doi: 10.1093/mnras/stv307
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513