# Quantifying the Effects of Known Unknowns on Inferred High-redshift Galaxy Properties: Burstiness, IMF, and Nebular Physics

Bingjie Wang, Joel Leja, Hakim Atek, Ivo Labbé, Yijia Li, Rachel Bezanson, Gabriel Brammer, Sam E. Cutler, Pratika Dayal, Lukas J. Furtak, et al.

► To cite this version:

Bingjie Wang, Joel Leja, Hakim Atek, Ivo Labbé, Yijia Li, et al. Quantifying the Effects of Known Unknowns on Inferred High-redshift Galaxy Properties: Burstiness, IMF, and Nebular Physics. The Astrophysical Journal, 2024, 963, 10.3847/1538-4357/ad187c. insu-04851581

HAL Id: insu-04851581

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

Submitted on 20 Dec 2024

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

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

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

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

# Quantifying the Effects of Known Unknowns on Inferred High-redshift Galaxy Properties: Burstiness, IMF, and Nebular Physics

Bingjie Wang (王冰洁)<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Joel Leja<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Hakim Atek<fnref n="4" /> ID, Ivo Labbé<fnref n="5" /> ID, Yijia Li (李轶佳)<fnref n="1" />,<fnref n="3" /> ID, Rachel Bezanson<fnref n="6" /> ID, Gabriel Brammer<fnref n="7" /> ID, Sam E. Cutler<fnref n="8" /> ID, Pratika Dayal<fnref n="9" /> ID, Lukas J. Furtak<fnref n="10" /> ID, Jenny E. Greene<fnref n="11" /> ID, Vasily Kokorev<fnref n="9" /> ID, Richard Pan<fnref n="12" /> ID, Sedona H. Price<fnref n="6" /> ID, Katherine A. Suess<fnref n="13" />,<fnref n="16" /> ID, John R. Weaver<fnref n="8" /> ID, Katherine E. Whitaker<fnref n="7" />,<fnref n="8" /> ID, and Christina C. Williams<fnref n="14" />,<fnref n="15" /> ID

Received 2023 October 9; revised 2023 December 21; accepted 2023 December 22; published 2024 February 29

## Abstract

The era of the James Webb Space Telescope ushers stellar population models into uncharted territories, particularly at the high-redshift frontier. In a companion paper, we apply the Prospector Bayesian framework to jointly infer galaxy redshifts and stellar population properties from broadband photometry as part of the UNCOVER survey. Here we present a comprehensive error budget in spectral energy distribution (SED) modeling. Using a sample selected to have photometric redshifts higher than 9, we quantify the systematic shifts stemming from various model choices in inferred stellar mass, star formation rate (SFR), and age. These choices encompass different timescales for changes in the star formation history (SFH), nonuniversal stellar initial mass functions (IMF), and the inclusion of variable nebular abundances, gas density, and ionizing photon budget. We find that the IMF exerts the strongest influence on the inferred properties: the systematic uncertainties can be as much as 1 dex, 2–5 times larger than the formal reported uncertainties in mass and SFR, and importantly, exceed the scatter seen when using different SED fitting codes. Although the assumptions on the lower end of the IMF induce degeneracy, our findings suggest that a common practice in the literature of assessing uncertainties in SED-fitting processes by comparing multiple codes is substantively underestimating the true systematic uncertainty. Highly stochastic SFHs change the inferred SFH by much larger than the formal uncertainties, and introduce ~0.8 dex systematics in SFR averaged over a short timescale and ~0.3 dex systematics in average age. Finally, employing a flexible nebular emission model causes ~0.2 dex systematic increase in mass and SFR, comparable to the formal uncertainty. This paper constitutes an initial step toward a complete uncertainty estimate in SED modeling.

Unified Astronomy Thesaurus concepts: Galaxy formation (595); H II regions (694); High-redshift galaxies (734); Initial mass function (796); Spectral energy distribution (2129); Star formation (1569); Galaxy evolution (594)

## 1. Introduction

Since the initial data releases from the James Webb Space Telescope (JWST), rapid progress has been made in photometrically identifying $z > 9$ candidates, leading to a new census of the early stages of galaxy formation (e.g., Finkelstein et al. 2022; Naidu et al. 2022; Donnan et al. 2023; Harikane et al. 2023). Among the surprising findings is the apparent overabundance of high-redshift sources in comparison to theoretical

predictions, or extrapolation of lower-redshift luminosity functions (Bouwens et al. 2023; Castellano et al. 2023; Ferrara et al. 2023; Mason et al. 2023; Mauerofer & Dayal 2023; Pérez-González et al. 2023). Spectroscopic follow-ups of some of the candidates generally confirm their high-redshift nature, but also reveal lower-redshift interlopers (Wang et al. 2023a; Arrabal Haro et al. 2023; Curtis-Lake et al. 2023).

Two questions immediately emerge from the early searches for $z > 9$ candidates: the reliability of photometric redshifts and uncertainties in their key properties as inferred by stellar population synthesis (SPS) models. The flexibility in SPS modeling enabled by Bayesian spectral energy distribution (SED)-fitting codes (Chevallard & Charlot 2016; Carnall et al. 2018; Johnson et al. 2021) means that these systematics can be investigated in detail. In a companion paper (Wang et al. 2024), we release a SPS catalog, consisting of ~50,000 redshifts and

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

<footnote n="1">¹ Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA; bwang@psu.edu</footnote>

<footnote n="2">² Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="3">³ Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="4">⁴ Institut d'Astrophysique de Paris, CNRS, Sorbonne Université, 98 bis Boulevard Arago, F-75014, Paris, France</footnote>

<footnote n="5">⁵ Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Melbourne, VIC 3122, Australia</footnote>

<footnote n="6">⁶ Department of Physics & Astronomy and PITT PACC, University of Pittsburgh, Pittsburgh, PA 15260, USA</footnote>

<footnote n="7">⁷ Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen, Jagtvej 128, København N, DK-2200, Denmark</footnote>

<footnote n="8">⁸ Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA</footnote>

<footnote n="9">⁹ Kapteyn Astronomical Institute, University of Groningen, 9700 AV Groningen, The Netherlands</footnote>

<footnote n="10">¹⁰ Department of Physics, Ben-Gurion University of the Negev, P.O. Box 653, Be'er-Sheva 84105, Israel</footnote>

<footnote n="11">¹¹ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="12">¹² Department of Physics & Astronomy, Tufts University, MA 02155, USA</footnote>

<footnote n="13">¹³ Kavli Institute for Particle Astrophysics and Cosmology and Department of Physics, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="14">¹⁴ NSF's National Optical-Infrared Astronomy Research Laboratory, Tucson, AZ 85719, USA</footnote>

<footnote n="15">¹⁵ Steward Observatory, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="16">¹⁶ NHFP Hubble Fellow.</footnote>
stellar population properties for sources in the Abell 2744 galaxy cluster field, as part of the UNCOVER survey (Bezanson et al. 2022). The UNCOVER SPS catalog is a large-scale application of the Prospector Bayesian framework (Johnson et al. 2021), and is representative of the state of the art in SPS modeling. Importantly, it provides a baseline for evaluating alternative SPS model choices in an exotic early universe now made eminently accessible by JWST.

SPS models include multiple ingredients, such as an initial mass function (IMF), stellar isochrones, and stellar spectra for the construction of simple stellar populations, star formation history (SFH) models, and stellar metallicity models for the construction of composite stellar populations, a model for dust attenuation and emission, and nebular continuum and line emission (see Conroy 2013 for a review, and Section 5.2 of Wang et al. 2024 for a discussion in the context of the UNCOVER SPS catalog). Even with data of exquisite quality and wavelength coverage at $z \lesssim 3$, the dependence of inferred parameters on modeling assumptions is already visible (Pacifici et al. 2023). At $z \gtrsim 10$, SPS models have yet to undergo robust tests. This implies that the dominant uncertainty for high-redshift measurements is not usually measurement uncertainty, but instead systematic uncertainty, which can potentially be orders of magnitudes higher than the former. The implications of different stellar libraries and functional forms of SFHs for JWST observations have been studied at some length (Whitler et al. 2023). This work focuses on a different set of systematics, encompassing the following three model choices.

**First, a nonuniversal IMF.** The relative numbers of stars as a function of the stellar mass influence most of the observable properties of stellar populations. The strongest effect is expected on the inferred stellar mass: at rest-frame optical/ultraviolet (UV) wavelengths for rising SFHs, we only observe $M \gtrsim 20 M_\odot$ stars (e.g., Figure 5 in Conroy 2013), meaning that a considerable fraction of the reported stellar mass relies on an extrapolation from the assumed IMF shape. Despite such importance, a complete theoretical understanding of the origin of the shape of the IMF is still lacking (see Bastian et al. 2010; Smith 2020 for recent reviews). Introduced in Salpeter (1955), the functional form of the IMF above the solar mass scale remains remarkably consistent with a slope of $\sim -1.3$. Since then, follow-up studies of the Galactic regions suggest deviations from the nominal IMF (Kroupa 2001; Chabrier 2003), but it is typically assumed to be constant across time and environment. There is, however, mounting evidence in the nearby universe against a universal IMF (Treu et al. 2010; Conroy & van Dokkum 2012; Geha et al. 2013; McWilliam et al. 2013; van Dokkum et al. 2017). Recent studies have claimed that IMF variations may be particularly relevant in the context of early star-forming activities (Chon et al. 2022; Katz et al. 2022; Pacucci et al. 2022; Steinhardt et al. 2023). In this work, we aim to quantify the systematic effects of IMFs that cover a reasonable range as suggested by observational studies.

**Second, outshining.** It has long been known that bursty star formation causes light from recent star formation to outshine that of older stellar populations (Papovich et al. 2001). Here burstiness refers to the pattern of sporadic star-forming activities interspersed with less active or quiescent phases. Such stochasticity is expected at high redshift as the feedback timescale approaches or exceeds the dynamical time of the system, decreasing and delaying the regulatory effect on the

star formation rate (SFR; Tacchella et al. 2016; Faucher-Giguère 2018; Dome et al. 2024). Observations also hint at bursty SFHs (Caputi et al. 2007; Smit et al. 2015; Guo et al. 2016; Looser et al. 2023). A recent bursty star formation in the detected high-redshift galaxies is also strongly motivated by their extremely blue UV slopes, which are explained by the ejection of dust due to sudden enhancement of radiation pressure (Tsuna et al. 2023; Ziparo et al. 2023). The recent SFH significantly influences galaxy SEDs across the electromagnetic spectrum and particularly in the rest-frame UV, meaning that the inference of galaxy properties depends on the assumed level of burstiness (Furlanetto & Mirocha 2022; Endsley et al. 2023). Various previous studies have investigated the effect of varying SFHs and priors (Carnall et al. 2018; Leja et al. 2019a; Suess et al. 2022; Duan et al. 2023). We focus on assumed burstiness in the context of nonparametric SFHs, and test it by varying the prior on the change in SFRs in adjacent age bins. Our bursty prior represents a more extreme case than what is implied by the bursty continuity prior in Tacchella et al. (2022). This is an intentional choice as the goal is to perform a sensitivity test rather than propose a new physical model.

**Third, nebular physics.** Line and continuum emission can make up a significant fraction of the total flux for stars at low metallicity and at young ages, up to 50% of the flux in broadband filters for highly star-forming galaxies (Pacifici et al. 2015), implying increased importance at high redshift (Charlot & Longhetti 2001; Anders & Fritze-v., 2003; Reines et al. 2010; Smit et al. 2014). The most sophisticated approach to incorporating nebular physics into SED fitting thus far relies on a model grid precomputed using photoionization simulations (Gutkin et al. 2016; Byler et al. 2017). While such an approach better accounts for variations in the recent SFH and nonsolar metallicity than analytic prescriptions (Fioc & Rocca-Volmerange 1999; Leitherer et al. 1999; Mollá et al. 2009), the accessible parameter space is still limited. The diverse chemical environment that may be present at early epochs (e.g., Bunker et al. 2023; Katz et al. 2023), and the possibility of strong optical nebular emission lines at $z < 6$ mimicking JWST/NIRCam photometry of high-redshift galaxies (McKinney et al. 2023; Zavala et al. 2023) renew the interest in the modeling of nebular emission with increased flexibility in elemental abundance ratios and separate parameterizations of gas-phase number density and ionizing photon density, while retaining consistency with the ionizing flux emitted by the model stellar populations. We do so by integrating a neural net emulator for the spectral synthetic code Cloudy (Chatzikos et al. 2023), dubbed Cue (Y. Li et al. 2024, in preparation), into Prospector.

The structure of this paper is as follows. Section 2 summarizes the JWST observations, including our sample of $z > 9$ candidates. Section 3 details our fiducial as well as new modeling approaches regarding the IMF, SFH prior, and nebular physics. Section 4 presents the effects of new models on the inferred parameters. We discuss the error budget as well as possible ways forward in Section 5, and conclude in Section 6.

Where applicable, we adopt the best-fit cosmological parameters from the Wilkinson Microwave Anisotropy Probe 9 yr results: $H_0 = 69.32$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_M = 0.2865$, and $\Omega_\Lambda = 0.7135$ (Hinshaw et al. 2013). Unless otherwise
mentioned, we report the median of the posterior, and 1σ error
bars are in the 16th and 84th percentiles.

## 2. Data

### 2.1. Photometry

This paper utilizes the latest photometric catalog (Weaver et al. 2024), made available as part of the second data release from the UNCOVER survey (Bezanson et al. 2022).<fnref n="17" /> The catalog includes all publicly available imaging data of Abell 2744 obtained from JWST/NIRCam, Hubble Space Telescope/Advanced Camera for Surveys (HST/ACS), and HST/WFC3. In particular, the JWST observations consist of three programs: UNCOVER (PIs Labbé & Bezanson, JWST-GO-2561; Bezanson et al. 2022), the Early Release Science program GLASS (PI: Treu, JWST-ERS-1324; Treu et al. 2022), and a Director’s Discretionary program (PI: Chen, JWST-DD-2756). These observations cover a wavelength range of ~1–5 µm in the observed frame, using eight filters: F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W. The HST data, obtained from the public archive, include HST-GO-11689 (PI: Dupke), HST-GO-13386 (PI: Rodney), HST-DD-13495 (PI: Lotz; Lotz et al. 2017), and HST-GO-15117 (PI: Steinhardt; Steinhardt et al. 2020). These additional observations cover a wavelength range of ~0.4–1.6 µm in the observed frame, using seven filters: F435W, F606W, F814W, F105W, F125W, F140W, and F160W.

### 2.2. Fiducial Redshifts and Stellar Population Properties

Fiducial redshifts and galaxy properties are taken from the UNCOVER SPS catalog (Wang et al. 2024). Here we only briefly reiterate the modeling components for completeness. All parameters, including redshift, are inferred jointly using the Prospector inference framework (Johnson et al. 2021), adopting the MIST stellar isochrones (Choi et al. 2016; Dotter 2016) and MILES stellar library (Sánchez-Blázquez et al. 2006) from FSPS (Conroy & Gunn 2010). A nonparametric form is assumed for the SFH, defined by mass formed in seven logarithmically spaced time bins (Prospector-α; Leja et al. 2017). Only the contribution of obscured active galactic nuclei (AGN) is considered following Nenkova et al. (2008a, 2008b), and parameterized by the normalization and dust optical depth in the mid-infrared (Conroy et al. 2009; Leja et al. 2018). A mass function prior, and a dynamic SFH (M, z) prior are included to optimize the photometric inference over the wide parameter space covered by deep JWST surveys (Wang et al. 2023b). Lensing magnification, computed from the public UNCOVER v1.1 lensing maps by Furtak et al. (2023), is performed on the fly to ensure consistency with the scale-dependent priors. The model consists of 18 free parameters, and sampling is performed using the dynamic nested sampler dynasty (Speagle 2020) with a speed up in model generation enabled by a neural net emulator dubbed parrot (Mathews et al. 2023).

### 2.3. High-redshift Candidates

We select $z_{phot} > 9$ candidates directly from the newly public SPS catalog, demonstrating the quality of the catalog as well as

the exciting parameter space now accessible via deep JWST surveys. Specifically, the criteria are

1. Redshift posteriors in the Prospector-β (Wang et al. 2023b) fits

(a) 50th percentile of the redshift posteriors $\ge 9.0$

(b) 0.1th percentile of the redshift posteriors $\ge 3.0$

2. Maximum-likelihood redshift $\ge 8.0$

3. $\chi^2_{red} = \chi^2/n_{bands} < 5$

4. Number of NIRCam filters $> 6$

5. S/N in photometry

(a) F444W $\ge 2.0$

(b) F410M $\ge 2.0$ if available

(c) F277W $\ge 6.0$

The first criterion defines our selection—criterion 1(a) imposes the redshift cut, whereas 1(b) decreases the probability of low-redshift interpolators given the multi-modality commonly seen in the posterior distribution. The rest of the criteria is to ensure the purity of the sample: $\chi^2_{red} < 5$ is applied to ensure a reasonable fit; that is, to remove data defects in the photometric catalog that are difficult to describe by SPS models. The signal-to-noise ratio (S/N) cut in the photometry is determined empirically after visual inspection of the 31 candidates given the first four criteria.

Our selection eventually yielded 18 candidates, which are listed in Table 1. In the same table, we crossmatch our sample to existing candidates in the literature, and include available spectroscopic redshifts as well. We note that the spectroscopic redshifts are systematically lower than the photometric redshifts by ~1.5σ–2σ, which may be a manifestation of the Eddington bias claimed in Serjeant & Bakx (2023). Further details on this sample are supplemented in the Appendix, which include the new candidates from this paper, a possible artifact, and discrepant photometric and spectroscopic redshifts.

In principle, the cutoff based on the maximum-likelihood redshifts is unnecessary if the posterior mass is located correctly and consistently in the full sample. However, as it is still in the early stage of photometry calibration for JWST, we impose this criterion for a cleaner sample. In the future, removing this would provide a more independent approach for selecting high-redshift candidates, complementing the traditional color selection technique (e.g., Castellano et al. 2022; Atek et al. 2023b).

We now proceed to investigate the effects of various modeling choices on the inferred redshifts and galaxy properties. It is worth noting that in order to reach a definitive conclusion regarding the signals of interest of this paper, relatively high S/N spectra with medium to high resolution would likely be necessary. For instance, rest-frame near-infrared spectra with S/N $\ge 100$ are needed for placing constraints on the IMF in the local universe (Conroy & van Dokkum 2012; van Dokkum & Conroy 2012). Given the difficulty in acquiring such data at high redshift, the scope of this paper is restricted to exploring the systematic effects on inferred properties when fitting broadband photometry.

## 3. Varying Known Unknowns in Stellar Population Inference

While the current inference framework is highly flexible, the known unknowns in the SPS modeling may induce systematic shifts in the inferred parameters comparable to or larger than the uncertainties in the fiducial model. The companion paper

<footnote n="17">¹⁷ All data products are accessible via https://jwst-uncover.github.io/#releases.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Model choices tested in this paper. (Left) The fiducial Chabrier IMF is plotted in black. A flat ($x = 1.8$) IMF is shown in orange, whereas a steep ($x = 2.8$) IMF is shown in green. $f$ is the ratio of the numbers of low-mass ($0.08 < M_{\odot} < 0.5$) to massive ($5 < M_{\odot} < 120$) stars. (Middle) Prior distributions of specific SFRs for a galaxy observed at $z = 10$ are plotted as functions of lookback time in gigayears. The fiducial rising SFH prior (Wang et al. 2023b) is illustrated by the black curves, while the extremely bursty prior is illustrated by the cyan curves. Shading indicates the 16th and 84th quantiles. (Right) The BPT diagram demonstrates the model coverage in the standard Cloudy grid (Byler et al. 2017) and in the more flexible model enabled by Cue. The gray dashed curve separating AGNs and starburst galaxies is taken from Kauffmann et al. (2003).</figcaption></figure>

<table><caption>**Table 1**
$z_{phot} > 9$ Candidates Studied in This Paper</caption><thead><tr><th>ID</th><th>R.A.<br>(J2000 deg)</th><th>decl.<br>(J2000 deg)</th><th>z<sub>phot</sub><sup>a</sup></th><th>M<sub>*</sub><sup>a</sup><br>(log M<sub>☉</sub>)</th><th>z<sub>spec</sub><sup>b</sup></th><th>Other IDs<sup>c</sup></th></tr></thead><tbody><tr><td>4745</td><td>3.61720</td><td>-30.42554</td><td>9.957<sup>+0.290</sup><sub>-0.301</sub></td><td>9.215<sup>+0.219</sup><sub>-0.150</sub></td><td>9.325 (3686; F23)</td><td>dhz1 / 2065</td></tr><tr><td>7720</td><td>3.60383</td><td>-30.41581</td><td>9.157<sup>+0.149</sup><sub>-0.250</sub></td><td>7.285<sup>+0.231</sup><sub>-0.193</sub></td><td></td><td></td></tr><tr><td>14019</td><td>3.57087</td><td>-30.40158</td><td>13.631<sup>+0.591</sup><sub>-0.466</sub></td><td>7.866<sup>+0.250</sup><sub>-0.170</sub></td><td>13.079 (UNCOVER-z13; W23)</td><td></td></tr><tr><td>14088</td><td>3.59250</td><td>-30.40146</td><td>10.104<sup>+0.097</sup><sub>-0.124</sub></td><td>7.526<sup>+0.295</sup><sub>-0.192</sub></td><td>9.880 (13151; F23)</td><td></td></tr><tr><td>27025</td><td>3.56707</td><td>-30.37786</td><td>10.975<sup>+0.226</sup><sub>-0.331</sub></td><td>7.834<sup>+0.235</sup><sub>-0.140</sub></td><td>10.071 (26185; F23)</td><td>uhz1 / 21623</td></tr><tr><td>28397</td><td>3.56152</td><td>-30.37494</td><td>10.098<sup>+0.201</sup><sub>-0.163</sub></td><td>7.059<sup>+0.149</sup><sub>-0.133</sub></td><td></td><td></td></tr><tr><td>30818</td><td>3.51193</td><td>-30.37186</td><td>9.868<sup>+0.750</sup><sub>-0.353</sub></td><td>8.881<sup>+0.366</sup><sub>-0.235</sub></td><td></td><td></td></tr><tr><td>38095</td><td>3.59011</td><td>-30.35974</td><td>10.928<sup>+0.250</sup><sub>-0.261</sub></td><td>8.163<sup>+0.172</sup><sub>-0.150</sub></td><td>10.255 (37126; F23)</td><td>ghz1 / 26928</td></tr><tr><td>39753</td><td>3.51356</td><td>-30.35680</td><td>13.474<sup>+0.816</sup><sub>-0.689</sub></td><td>8.047<sup>+0.267</sup><sub>-0.205</sub></td><td>12.393 (UNCOVER-z12; W23)</td><td>39074</td></tr><tr><td>42485</td><td>3.51374</td><td>-30.35157</td><td>10.272<sup>+0.521</sup><sub>-0.608</sub></td><td>8.391<sup>+0.224</sup><sub>-0.221</sub></td><td></td><td>ghz4</td></tr><tr><td>42658<sup>d</sup></td><td>3.52286</td><td>-30.35129</td><td>18.915<sup>+0.432</sup><sub>-0.445</sub></td><td>8.207<sup>+0.181</sup><sub>-0.158</sub></td><td></td><td>ghz9 / 52008</td></tr><tr><td>45758</td><td>3.47874</td><td>-30.34554</td><td>11.067<sup>+0.248</sup><sub>-0.257</sub></td><td>8.477<sup>+0.222</sup><sub>-0.185</sub></td><td></td><td></td></tr><tr><td>46843</td><td>3.50731</td><td>-30.34321</td><td>11.119<sup>+0.398</sup><sub>-0.450</sub></td><td>7.898<sup>+0.216</sup><sub>-0.174</sub></td><td></td><td></td></tr><tr><td>54328</td><td>3.49899</td><td>-30.32476</td><td>12.559<sup>+0.043</sup><sub>-0.054</sub></td><td>8.673<sup>+0.158</sup><sub>-0.161</sub></td><td>12.117 (GLASS-z12; B23)</td><td>ghz2 / 70846</td></tr><tr><td>55541</td><td>3.45142</td><td>-30.32181</td><td>11.341<sup>+0.253</sup><sub>-0.283</sub></td><td>8.412<sup>+0.187</sup><sub>-0.130</sub></td><td></td><td>ghz8 / 73667</td></tr><tr><td>55950</td><td>3.45136</td><td>-30.32072</td><td>10.802<sup>+0.310</sup><sub>-0.869</sub></td><td>8.307<sup>+0.219</sup><sub>-0.154</sub></td><td></td><td>ghz7 / 81198</td></tr><tr><td>57149</td><td>3.45471</td><td>-30.31689</td><td>9.985<sup>+0.462</sup><sub>-0.452</sub></td><td>9.043<sup>+0.321</sup><sub>-0.280</sub></td><td></td><td>83338</td></tr><tr><td>61062</td><td>3.49177</td><td>-30.29831</td><td>9.118<sup>+0.124</sup><sub>-0.149</sub></td><td>8.461<sup>+0.165</sup><sub>-0.160</sub></td><td></td><td></td></tr></tbody></table>

**Notes.**

<sup>a</sup> Posterior moments of photometric redshifts and stellar masses from Prospector fits.

<sup>b</sup> Spectroscopic redshifts. All are determined from NIRSpec/Prism data, except GLASS-z12, which is based on an [O III] detection from the Atacama Large Millimeter/submillimeter Array (Bakx et al. 2023). F23 refers to Fujimoto et al. (2023), and W23 refers to Wang et al. (2023a).

<sup>c</sup> Photometrically identified candidates in the literature. Alphabetical IDs are from Castellano et al. (2022), and numerical IDs are from Atek et al. (2023a); GLASS-z12 is also selected in Naidu et al. (2022).

<sup>d</sup> Possible data quality issue; see the Appendix for details.

(Wang et al. 2024) discusses some of the assumptions that are particularly relevant to the UNCOVER SPS catalog. This work further examines the influences of a steeper/flatter IMF, bursty star formation, and more flexible nebular emission modeling, as summarized in Figure 1. We elaborate on both fiducial and alternative models below.

### 3.1. IMF

The fiducial IMF is taken from Chabrier (2003). We consider two alternative IMFs described by a power law of slope $x$ over

the full mass range as

$$ \xi \propto m^{-x}. \qquad (1) $$

Two variations are examined, $x = 1.8$ and $2.8$, covering a reasonable range as suggested from observations (Cappellari et al. 2012; Martín-Navarro et al. 2015; van Dokkum et al. 2017). For the main analyses, we define all the IMFs over the interval $0.08 < M_{*} < 120 M_{\odot}$ for a controlled experiment. However, the minimum ($m_l$) and maximum stellar mass cutoffs remain uncertain. The low-mass end of the IMF is particularly
difficult to constrain since these stars are largely invisible
outside the Milky Way. We thus additionally consider two
cases where $m_l = 0.5$ and $1 M_\odot$ (Kroupa 2001; Chabrier 2003).

Since we intend to test the scenario in which the IMF varies
across cosmic time, and changing the IMF affects the scale of
the predicted spectrum, it is necessary to modify the mass prior
in Wang et al. (2023b) accordingly. We find the normalization
factor by minimizing the $\chi^2$ between spectra predicted by
different IMFs, and simply scale the mass prior by this factor.
Strictly speaking, IMF($z$) describes an instantaneous SFR, but
the stellar mass is determined by effective IMF from all stars
formed. We thus note that our modification on the mass prior
serves only as an approximation to mitigate the bias that may
be induced by the scale difference among different IMFs. The
fact that we assume the same IMF over the entire history
justifies this simplification.

Effectively, the same stellar mass function prior is applied to all
fits. Therefore, our choice of the prior does not influence the
systematic effects to be quantified as these are the relative shifts in
the inferred parameter with respect to the fiducial values. In
addition, the mass prior is constructed by using the observed mass
functions between $0.2 < z < 3$ from Leja et al. (2020); for $z > 3$, we
adopt the nearest-neighbor solution, i.e., the $z=3$ mass functions,
in the absence of reliable high-resolution rest-frame optical selected
mass functions at $z > 3$. As mentioned in Wang et al. (2023b), this
choice allots a conservatively high probability for yet-to-be-
discovered populations of high-mass, high-redshift galaxies.

## 3.2. Timescales for Changes in SFH

All SFHs are modeled non-parametrically via mass formed
in seven logarithmically spaced time bins (Leja et al. 2017). In
our fiducial setup, we use the dynamic SFH prior from Wang
et al. (2023b). The expectation values of the logarithmic SFR
ratios in adjacent time bins are based on the observed cosmic
SFR density in the empirical UniverseMachine model (Beh-
roozi et al. 2019), favoring rising SFHs in the early universe
and falling SFHs in the late universe, and adjusted as a function
of mass to reflect downsizing (Cowie et al. 1996; Thomas et al.
2005). The distribution of the logarithmic SFR ratios is
modeled as a student's t- distribution.

An additional bursty SFH is generated by replacing the
student's t-distribution with a uniform prior spanning from -5
to 5. This means that the SFR is effectively uncorrelated
between the time bins, and can jump up to 10 orders of
magnitude in a single time step. Our prior is similar in concept
to the bursty continuity prior in Tacchella et al. (2022), but
represents a more extreme bursty scenario.

## 3.3. Nebular Physics

Prospector models the nebular emission self-consistently
via a precomputed Cloudy grid (Byler et al. 2017). The
parameter space is limited for a practical reason: grid-based
approaches necessarily suffer from the curse of dimensionality.

We replace the grid with a neural net emulator, dubbed Cue
(Y. Li et al. 2024, in preparation), which emulates the
continuum and line emission from a single H II region based
on the spectral synthetic code Cloudy version 22.00
(Chatzikos et al. 2023). Cue is tailored to flexibly model
exotic chemistry and unusual ionizing properties of galaxies,
and is well suited for galaxies in the early universe. It models
the ionizing input spectra as a flexible four-part piecewise

continuous power law, along with freedom in gas density, total
ionizing photon budget, O/H, C/O, and N/O. While in
principle this power-law approximation allows us to infer the
ionizing spectrum separately from the model stellar spectra,
here we fit this piecewise power law to the stellar ionizing
spectra from the model for each spectrum generated. With the
exception of carbon and nitrogen, the element abundances are
scaled with the oxygen abundance, based on solar abundances
and dust depletion factor specified in Dopita et al. (2000).

The newly gained flexibility, however, does not a priori
forbid unlikely nebular abundances. Certain chemical ratios are
preferred on galactic scales due to the production rate of
various elements; for example, N/O as a function of oxygen
abundance has a well-defined shape following from the
secondary production of nitrogen. We thus place a Gaussian
prior on the N/O ratio against oxygen abundance based on the
observed trend. The Gaussian mean and standard deviation are
estimated from the data compiled by Groves et al. (2004). In
the low oxygen abundance regime where no data is available,
we take values from the model grid in Gutkin et al. (2016).

# 4. Results

The effect of these changes on the stellar mass, SFR, and
stellar age are summarized in Figure 2, where these key
parameters inferred from alternative models are plotted against
the fiducial values. Furthermore, the preference of a given
model by the data is quantified via the Bayesian evidence in
Figure 3. The results regarding each model are presented in
more detail below.

## 4.1. Alternative IMFs

Changing the slope of the IMF induces a systematic shift in
the inferred stellar mass and SFR by ~0.3–1 dex, as suggested
in Figure 2. It has a negligible effect on the inferred stellar age.
The effects of varying the minimum stellar mass cutoff are
shown in the same way in Figure 4. More specifically,
assuming a top-heavy ($x=1.8$) IMF systematically decreases
the inferred stellar mass and SFR by $\gtrsim$0.3 dex; increasing $m_l$
from 0.08 $M_\odot$ to 0.5 or $1 M_\odot$ further enhances the systematic
shifts by $\gtrsim$0.1 dex. Assuming a top-light ($x=2.8$) IMF
increases the inferred stellar mass and SFR by ~1 dex;
increasing $m_l$, however, significantly decreases the systematic
shifts. With $x=2.8$, $m_l \approx 1 M_\odot$, we recover the fiducial
Chabrier masses and SFRs. Our results are consistent with
Woodrum et al. (2023), in which different IMFs are used to fit a
subsample in the JADES photometric catalog (Eisenstein et al.
2023; Rieke et al. 2023).

Importantly, the Bayes factors (Trotta 2008), plotted in
Figure 3, indicate that the data are almost completely agnostic
to the different forms of IMFs. Further intuition can be gained
from examining a single fit. This is illustrated in Figure 5,
where the posterior median model photometry and spectra from
different IMF cases are plotted. All models can fit the observed
photometric data equally well, and the model spectra only
differ marginally from each other. The longer wavelength range
covered by JWST/MIRI may offer additional distinguishing
power. Similar results are found for the rest of the sample.

## 4.2. Bursty Star Formation Prior

The bursty SFH prior of this paper causes negligible
systematic changes in the inferred stellar masses. However, it
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Changes in inferred parameters due to different model choices. The stellar mass, SFR, and mass-weighted age (MWA) are plotted against the fiducial values in each panel. (Upper panel) The comparison between results assuming an IMF with $x = 1.8$ ($x = 2.8$) and the fiducial Chabrier IMF are shown in orange (green). Significant systematic shifts in stellar mass and SFR are found. (Middle panel) The fiducial model assumes a rising SFH prior (Wang et al. 2023b), while the alternative model assumes a bursty prior. The SFRs averaged over the most recent 10 Myr (SFR$_{10}$, unfilled circles) show significant scatter driven by the stochasticity in the SFH, but the SFRs averaged over the most recent 100 Myr (SFR$_{100}$, filled circles) are less affected. The bursty prior also causes a systematic decrease in stellar age. (Lower panel) The fiducial model computes nebular emission and continuum from a precomputed Cloudy grid (Byler et al. 2017), while the alternative model uses a Cloudy emulator dubbed Cue. The extra flexibility enabled by Cue systematically increases the inferred mass and SFR by ~0.1–0.2 dex, similar in size to the $1\sigma$ measurement uncertainties indicated by the error bars.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Model comparison. Logarithm of the Bayes factor comparing different model choices is calculated for each candidate. A value of $B'/B > (<)1$ represents an increase (decrease) of the support in favor of the alternative model given the observed data (Trotta 2008). The boundaries between different shading indicate weak evidence ($|B'/B| = 1$), moderate evidence ($|B'/B| = 2.5$), and strong evidence ($|B'/B| = 5$), respectively. In most cases, the data are agnostic of the model choices.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Same as Figure 2, but showing the changes in the inferred parameters with increasing minimum stellar mass cutoffs, $m_l$, of a given IMF.</figcaption></figure>

tends to increase the uncertainty in the inferred SFR averaged over the most recent 10 Myr, SFR$_{10}$, with a subset of the fits having these recent SFRs that significantly deviate from their fiducial values in Figure 2. In contrast, the inferred SFRs averaged over a longer timescale of the most recent 100 Myr, SFR$_{100}$, are much less affected. Additionally, the bursty prior leads to a systematic decrease in the inferred average stellar ages of all galaxies, of order 0.2–0.5 dex.

The left panel of Figure 6 shows the inferred SFHs for the object where the most dramatic change in SFR is observed. The increase in the stochasticity in the SFH due to the bursty prior is clearly visible, and at almost all lookback times the SFR is systematically changed by much larger than the size of the formal uncertainties.

Similar to the IMF case, diagnostics based on the Bayes factor, as well as the median model photometry and spectra show that the observed photometric data have little power to discriminate between bursty and smooth SFHs in general. This is in spite of the fact that here we contrast extremely opposite scenarios of an extremely bursty prior with a relatively smooth prior. We note that this null result differs from that of a recent study, which claims that stochastic SFHs are found at $z \sim 9$ (Ciesla et al. 2023) using the JADES photometric catalog (Eisenstein et al. 2023; Rieke et al. 2023). Given the relatively weak constraining power of broadband photometry coupled with the uncertain redshifts, more informative data, in particular spectroscopic observations, are likely required to reach a definitive conclusion.

### 4.3. Flexible Nebular Emission Modeling

The more flexible nebular emission modeling, enabled by the Cloudy emulator, Cue, causes differences in the inferred

masses and SFRs that are within $1\sigma$ of the formal uncertainties. However, consistent increases in the masses and SFRs are found of order ~0.2 dex. We note that this mass offset is more significant than an extrapolation from studies at lower redshift (e.g., the ~0.12 dex difference found in Li & Leja 2022), which is a natural consequence of the $z > 9$ population being more sensitive to nebular modeling than the regular galaxy population.

It is also a curious finding that the scatter in the model comparison plot of Figure 3 is noticeably larger than those of the above two cases—this suggests that, unlike the IMF and SFH timescales, nebular properties may well be able to be constrained from broadband photometry, and this constraining power would be substantially increased by medium-band photometry.

To gain more insights into the influence of the flexible model, we show the emission line ratios calculated from posterior samples for one of the candidates in Figure 6. The recent measurement of GN-z11 is included for reference (Bunker et al. 2023). Interestingly, the line ratios of GN-z11 seem to occupy a space that is better sampled by Cue. We note that this particular galaxy is claimed to host a bright AGN, which also explains its peculiar abundance (Maiolino et al. 2023; Scholtz et al. 2023). Line emission due to AGN is capable of being modeled with Cue; however, here we assume the ionizing spectra of our galaxies are dominated by stars, and use Cue such that the ionizing input follows the stellar ionizing input. We opt not to allow full freedom in the ionizing spectra since the photometric data lack the constraining power to distinguish between the stellar and AGN contributions. The ability to sample the parameter space occupied by GN-z11 thus indicates that the emulator is flexible enough to model a stellar
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Model photometry and spectra for one of the galaxies as an example. (a) Photometric data are shown as gray squares. Median model spectra assuming the different IMFs are plotted as functions of observed wavelengths in the left panel, and the shading indicates the 16th and 84th quantiles. Median model photometries are overplotted in the same colors. The right panel is the same, but compares the fiducial model with one that assumes a bursty prior. In all cases, the model photometric points fit the data equally well, and appear to be identical, whereas the model spectra exhibit marginal difference. (b) Same as the upper panel, but includes the longer wavelength range covered by JWST/MIRI. The lower panels here show the differences between the photometries of the alternative models and the fiducial model, normalized by the $1\sigma$ uncertainty in the fiducial model photometry.</figcaption></figure>

ionizing spectra mimicking an AGN spectrum. These results suggest that while more flexible nebular models have a relatively less strong effect on bulk properties inferred from SED fitting than key uncertainties like the IMF and SFH timescales, they will be very important in interpreting more detailed spectroscopic observations—and they can be observationally constrained. We discuss the various implications of these results in the following section.

# 5. Discussion

## 5.1. Systematic Uncertainties from Fixed Parameters in SED Modeling

Two immediate conclusions can be reached from Section 4: first, changing the IMF and the SFH prior produces substantial systematic effects on some of the key stellar population

parameters; second, the broadband photometric observations available from JWST cannot distinguish between the model choices.

It is thus important to ask, given the reasonable model choices and data, what the *real* uncertainties are for high-redshift galaxies. One way to address this question is by comparing the systematic uncertainties to the reported formal uncertainties. Here a formal uncertainty refers to the 16th–84th percentiles of posterior distributions from a single fit, where the different model choices cannot yet be marginalized. We compare the two kinds of uncertainties as follows.

For each fit, we subtract the medians from the posteriors, which are then stacked to provide an estimate of the average formal uncertainty. The resulting distribution for each of the inferred parameters is well approximated by a Gaussian function. We quantify the systematic uncertainty in a similar
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Inferred properties for one of the galaxies as an example. (Left) The inferred SFH assuming a rising SFH prior (Wang et al. 2023b) is plotted as a function of lookback time in gigayears in black, whereas that assuming a bursty prior is plotted in cyan. The bursty prior increases stochasticity as expected, and importantly, the SFR is systematically changed by much larger than the size of the formal uncertainties shown in gray shading at almost all lookback times. (Right) Emission line ratios calculated from posterior samples based on the standard Cloudy grid (Byler et al. 2017) are shown in black, whereas those from the Cloudy emulator, Cue, are shown in purple. The measurement of GN-z11 is included for reference (Bunker et al. 2023). Cue enables the exploration of a wider parameter space in nebular conditions.</figcaption></figure>

manner, except that we subtract the median calculated from the fiducial fit from the posteriors. A Gaussian function is fitted to the stacked posteriors, of which the mean corresponds to the average systematic shift, and the standard deviation indicates the range of possible systematic shifts.

The typical formal and systematic uncertainties of stellar mass, SFR, and stellar age are reported in Figure 7. For the case of a nonuniversal IMF, the usually unaccounted systematic uncertainty in stellar mass (0.4 ± 0.3 dex for a flatter, top-heavy IMF, and 1.0 ± 0.2 dex for a steeper, top-light IMF) can be two to five times larger than the formal uncertainty (~0.2 dex). The same applies to SFR, where the formal uncertainty is found to be ~0.2 dex, but the systematic uncertainties are 0.3 ± 0.3 and 1.0 ± 0.1 dex for a top-heavy and a top-light IMF, respectively. In addition, theoretical models have indicated even more exotic IMF shapes at high redshift (e.g., a flat IMF in Chon et al. 2022), or possible environmental dependences (Gunawardhana et al. 2011; Jeřábková et al. 2018; Yan et al. 2019). Here we only intend to cover the range of the IMF slopes suggested by observations.

Critically, a steep IMF, motivated on both theoretical (Chabrier et al. 2014) and observational grounds (Conroy et al. 2017; van Dokkum et al. 2024), leads to systematics that always dominates over the error budget from all other uncertainties normally considered. These values are also greater than the scatter caused by fitting a sample of galaxies with different SED-fitting codes, which are found to be ~0.1 and ~0.2 dex for stellar mass and SFR, respectively (Pacifici et al. 2023). This means that a common practice in the literature of assessing uncertainties in SED fitting processes by comparing multiple codes is substantively underestimating the true uncertainty, as these codes typically use similar assumptions for the IMF, the smoothness of the SFH, and nebular emission.

We note that the minimum and the maximum stellar mass cutoffs, which are held fixed in the above analyses, may further affect the systematic uncertainties. Of particular interest is the low-mass end of the power law, as for a steeper, top-light IMF, the stellar mass would be dominated by stars around the lower limit of the IMF. Thus far, we focus on the discussion on the varying IMF slopes because low-mass stars make up a

negligible fraction of the total light. We additionally include Figure 8 showing the systematic uncertainties with increasing lower mass limits to illustrate the degeneracy between $m_l$ and inferred stellar mass. For the top-light IMF, assuming $m_l = 0.5 M_\odot$, the systematic uncertainties of stellar mass and SFR decrease to 0.3 ± 0.2 dex, comparable to the formal uncertainties; assuming $m_l = 1 M_\odot$, the inferred parameters exhibit negligible differences compared to the fiducial Chabrier values. It is, however, also easy to imagine performing a similar analysis for IMFs with increasing hidden mass, either owing to lower $m_l$ or steeper slopes at the low-mass end. These alternatives are difficult to be ruled out as low-mass stars would be largely invisible outside the Milky Way, and would increase the systematic uncertainties. For the top-heavy IMF, the results are less sensitive to the lower limit. The stellar masses are systematically decreased by 0.5 ± 0.4 (0.6 ± 0.4) dex, and the SFRs are similarly decreased by 0.4 ± 0.4 dex, if $m_l = 0.5 M_\odot$ ($1 M_\odot$).

The same analysis is performed for the SFH priors, which assume different timescales for changes. Also plotted in Figure 7, we find that stochastic SFHs, which are expected at high redshift (Tacchella et al. 2016; Faucher-Giguère 2018; Legrand et al. 2022; Dome et al. 2024), produce significant systematics in SFR averaged over the most recent 10 Myr (0.8 ± 1.5 dex), and in stellar age (0.3 ± 0.3 dex). The formal uncertainties in these two parameters are ~0.2 and ~0.3 dex, respectively. It is worth noting that the systematics in the SFR can be mitigated by averaging the SFR over a longer timescale. Using the most recent 100 Myr, the systematic uncertainty in SFR decreases to 0.1 ± 0.5 dex. In addition, it may be surprising to see that the bursty prior does not lead to significant systematics in the inferred stellar mass (0.0 ± 0.4 dex). This is likely because of the construction of our fiducial SFH prior from Wang et al. (2023b), which has been shown to recover stellar mass reasonably well even in the presence of very bursty SFHs from simulations (Narayanan et al. 2024).

In contrast to the substantial systematic shifts found in the above two cases, a more flexible nebular emission model results in systematic uncertainties in mass and SFR (0.2 ± 0.3 dex) that are comparable to the formal uncertainty
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Formal reported uncertainties vs. systematics. (Upper panel) The black curve is a Gaussian fit to the stacked posteriors with the median subtracted for the Chabrier stellar masses. The fit suggests that the typical reported uncertainty, σ, is about 0.20 dex. The orange (green) curve is a Gaussian fit to the stacked posteriors of masses assuming x = 1.8 (x = 2.8) IMF, with the median of the Chabrier masses subtracted. The mean, -0.38 (0.97), describes the typical systematic shift after changing the IMF, whereas the standard deviation, 0.31 (0.16), indicates the spread in the systematics. Overplotted as the gray shade is the scatter driven by using different SED-fitting codes (Pacifici et al. 2023). The other two plots contrast the uncertainties and the systematics in the SFRs and stellar ages, respectively. Critically, the systematics in the inferred stellar mass and SFR due to different IMF choices can be at least comparable to and often larger than the standard uncertainties (posterior moments), and the scatter seen when fitting the object with different codes. (Middle panel) The same as the upper panel, but showing the systematics introduced by adopting a bursty prior instead of a rising SFH prior (Wang et al. 2023b). The stochasticity in the SFH leads to a large spread in the SFR posteriors. The typical systematic uncertainty in SFR average over the most recent 10 Myr is of order ~0.84 dex, and with a large spread of ±1.47 dex. However, SFR averaged over the most recent 100 Myr is less prone to the change in priors. The typical systematic uncertainty decreases to 0.05 dex, with a spread of ±0.50 dex. The typical systematic uncertainty in age is -0.31 ± 0.34 dex. (Lower panel) The typical systematic uncertainty in inferred mass introduced by modeling the nebular physics with more flexibility while retaining consistency with the model stellar populations is comparable (~0.16 dex) to the formal uncertainty. A similar systematic uncertainty in SFR (~0.15 dex) is also observed.</figcaption></figure>

~0.2 dex). The smaller systematic changes mean that even in the case of very high emission line equivalent widths expected in the first galaxies, the assumption of solar-scaled nebular abundance patterns, which are not at all expected to hold at early times, do not cause a factor of 2 or more systematics in the inferred masses, SFRs, and ages. Furthermore, the larger

scatter in the Bayesian evidence and the differences in the observed line ratios suggest that this is a solvable systematic with future observations—spectroscopy, or even photometry with relatively higher spectral resolution.

A counterintuitive finding, however, is that the estimates for stellar mass and SFR of one galaxy, ID 14088, using the model
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Same as Figure 7, but showing the uncertainties for the cases assuming different lower limits, $m_l$, of the IMF. For the top-heavy IMF, increasing $m_l$ leads to larger systematic uncertainties in inferred stellar mass and SFR. For the top-light IMF, increasing $m_l$ leads to smaller systematic shifts owing to less hidden mass.</figcaption></figure>

with increased flexibility are much stronger constrained than the original model. In addition, the Bayesian evidence strongly favors the original model for this particular galaxy. Upon examination, the fiducial fit has a bimodal posterior distribution for the stellar-phase metallicity, but one of the solutions is disfavored by Cue, thus decreasing the error bars. A possible reason is that while Cue finds a marginally better single mode with increased flexibility, the likelihood of this mode is not high enough to compensate for the added four parameters.

Finally, we note that the contribution from AGN can have a significant influence on the stellar mass estimates (D'Silva et al. 2023). In this work, we only model the mid-infrared emission from the AGN, and find that our sample has low AGN contribution. The mean AGN to galaxy bolometric luminosity is ~0.004. However, we caution that the wavelength coverage of our photometric data is not expected to process constraining power for the particular model considered. The AGN modeling in Prospector will be improved in future works.

## 5.2. Implications for Early Galaxy Formation and Possible Ways Forward

Given the significant systematics in parameters inferred by SED fitting combined with the lack of short-term prospects for observational constraints on said parameters, one could argue for including these systematic uncertainties in our error budget for all galaxies. Indeed, the IMF-induced systematics are often larger than the uncertainties and systematics usually considered in the stellar mass functions (Behroozi et al. 2010; Conroy 2013; Speagle et al. 2014; Grazian et al. 2015; Santini et al. 2015; Carnall et al. 2018; Leja et al. 2019b; Furtak et al. 2021).

More recently, various papers have challenged the validity of galaxy formation models or ΛCDM in light of newly discovered galaxy populations with JWST (e.g., Casey et al. 2023; Labbé et al. 2023; Xiao et al. 2023). A nonuniversal IMF can serve as a possible explanation for the overly massive galaxies, as pointed out in Boylan-Kolchin (2023) and Steinhardt et al. (2023). Taking into account the systematic uncertainties found in this paper can, for example, bring the star formation efficiency expected from the Labbé et al. (2023) sample into 1σ consistency with a reasonable value of ≲0.32 at $z>8$, and that expected from the Xiao et al. (2023) sample to ≲0.2 at $z>5$, thus removing the need for exotic cosmological models (e.g., Gong et al. 2023; Padmanabhan & Loeb 2023; Parashari & Laha 2023). This is intriguing given that the modifications in the power spectrum necessary for solving the tension with ΛCDM can induce conflicts at lower redshifts (Sabti et al. 2024).

Unfortunately, observationally constraining the IMF beyond the Milky Way is a known challenge. Thus far most of the meaningful constraints come from independent measurements of the mass-to-light ratio from gravitational lensing (Treu et al. 2010; van Dokkum et al. 2024) or dynamical modeling (Cappellari et al. 2013; Posacki et al. 2015), or from direct observational signatures of low-mass stars in high S/N, R ~ 3000 near-IR spectra (Conroy & van Dokkum 2012; van Dokkum et al. 2017). These approaches are difficult to generalize to a statistical sample at a cosmological distance.

Therefore, it seems necessary to fully incorporate the systematic uncertainties in the error budget. This is, however, nontrivial via conventional techniques. To start, including IMF
as a free parameter requires rebuilding the stellar library in every fit, which makes it infeasible to sample the likelihood surface within a reasonable amount of time. Possible solutions include machine-learning accelerated SED fitting of individual galaxies (Hahn & Melchior 2022; Wang et al. 2023c), or population-level inference (Alsing et al. 2023; Li et al. 2024). For the former, once a suitable training set is constructed, it becomes straightforward to marginalize over IMF variations. The latter approach directly makes the population distribution the inference objective, which drastically increases the computational efficiency.

The overabundance of high-redshift galaxies is another interesting discovery from early JWST observations (e.g., Yung et al. 2024). Bursty star formation has been invoked as a possible explanation (Pallottini & Ferrara 2023; Shen et al. 2023; Sun et al. 2023). While we find that adapting a bursty prior can fit the photometric observations equally well, the current data do not offer evidence for or against such expectations. Understanding the prevalence of bursty star formation, and quantifying the resulting selection effects (i.e., galaxies oscillating in and out of the observable sample due to SFR variability) likely requires at least spectroscopy in order to measure the different SFR and SFH indicators that can constrain burstiness. Perhaps the most convincing evidence would stem from exploring the population distributions of these spectral indicators.

As for nebular physics, the Cloudy emulator, Cue, allows for the exploration of a wide parameter space while retaining consistency with stellar populations; particularly relevant to this work is the Baldwin-Phillips-Terlevich (BPT) region occupied by GN-z11, which is sparsely sampled in the standard Cloudy grid (Byler et al. 2017). This indicates that, with more informative spectroscopic or medium-band photometric data, the flexible model is a promising approach to characterize the exotic nebular conditions in the early universe.

Having discussed the systematic effects on the inferred galaxies properties, we would like to add that none of the alternative models leads to dramatic changes in the inferred photometric redshifts. For the *z* > 9 range studied in this work, the photometric redshifts are almost entirely determined by the Lyman break. The Cue parameterization does not influence the broadband photometry substantially enough to change the results. This means that the uncertainties in the photometric redshifts are mainly driven by the modeling of Lyα and its interaction with the intergalactic medium, although nebular emission lines such as Hα have been shown to be able to resemble the Lyman break in photometric data by spectroscopic follow-ups (Arrabal Haro et al. 2023). In the near future, we plan to incorporate a phenomenological model to describe the radiative transfer process (e.g., damping wings; Curtis-Lake et al. 2023). With that being said, the transfer of resonant lines in gas is an intricate process, the accurate description of which requires numerical radiative transfer simulations (e.g., Gronke et al. 2017; Michel-Dansac et al. 2020). There is much room for improvement to fully capture the complexity of the propagation of Lyα.

Finally, we note that the systematics analyzed for the *z*phot > 9 sample of this paper are expected to be present across redshifts, albeit with varying degrees of importance. On the observational front, JWST has been revealing new galaxy populations and establishing statistical samples, reaching uncharted regions in the cosmic history. Sophisticated

modeling of galaxy populations that encapsulates the full uncertainties, analogous to the framework proposed for redshift distribution inference in Alsing et al. (2023), would be invaluable in forming a coherent narrative of galaxy evolution over the observed dynamic range.

## 6. Conclusions

In this paper, we quantify the systematic uncertainties from different model choices on inferred high-redshift galaxy properties. We conduct our experiment on a *z* > 9 sample selected from the publicly available UNCOVER SPS catalog (Wang et al. 2024). As a large-scale application of the Prospector Bayesian inference framework, where redshifts and stellar population parameters are inferred jointly, the UNCOVER SPS catalog provides a baseline against which the fundamental assumptions in SED modeling can be tested. In particular, we investigate the three model choices as follows.

First, we fit the photometry assuming a flatter/steeper IMF in addition to the fiducial Chabrier IMF. The photometric data are firmly unable to distinguish between particular forms of IMF. We find that a flatter IMF (*x* = 1.8) tends to systematically decrease the inferred stellar mass by ~0.4 dex, and SFR by ~0.3 dex. A steeper (*x* = 2.8) IMF typically increases the inferred stellar mass and SFR both by ~1.0 dex. These values are notably larger than the reported uncertainty of ~0.2 dex on mass and SFR from fits assuming a fixed IMF, and also the scatter found when fitting data using different SED codes (~0.1 dex in mass and ~0.2 dex in SFR; Pacifici et al. 2023). However, these results can be sensitive to the lower mass limits of the IMF, which increase or decrease the systematic uncertainties depending on the change in the amount of hidden mass. Observationally, this is hard to measure, since the total integrated light is dominated by massive stars. Nevertheless, taken together with the lack of IMF constraints at high redshift and the motivations for a nonuniversal IMF on observational as well as theoretical grounds, we infer that the usually reported uncertainties on stellar mass and SFR are likely to be underestimated.

Second, we adopt an extremely bursty star formation prior, in contrast to the fiducial rising SFH prior proposed in Wang et al. (2023b). Similar to the varying IMF case, the Bayes factor shows no particular preference given the photometric data. The bursty prior leads to >1 dex deviations from the fiducial SFRs averaged over the most recent 10 Myr in some cases, however, the large systematics can be mitigated by averaging the SFRs over longer timescales, e.g., 100 Myr. The typical systematic decrease is ~0.8 dex in inferred SFR₁₀, which is greater than the formal uncertainty of ~0.2 and ~0.1 dex in SFR₁₀₀. Therefore, care must be taken when attempting to infer burstiness by comparing SFRs averaged over different timescales. Additionally, the bursty prior results in marginally younger stellar ages for the full population, of an order of ~0.3 dex, comparable to the formal uncertainty. It is worth emphasizing that, while being able to fit the observations equally well, the inferred SFR assuming the bursty prior is systematically changed by much larger than the size of the formal uncertainties at almost all lookback times. This means that the inference on SFH from broadband photometry is prior dominated rather than likelihood dominated; in other words, choosing the appropriate prior is critical for SFH studies, as pointed out in Leja et al. (2019a), Suess et al. (2022), Tacchella et al. (2022).
Third, we implement a Cloudy emulator (Cue; Y. Li et al. 2024, in preparation.) into Prospector, which allows us to model the exotic nebular conditions in the early universe while retaining consistency with the assumed stellar populations. This results in systematic uncertainties in mass and SFR (~0.2 dex) that are comparable to the formal uncertainty. Interestingly, utilizing Cue also leads to better explored parameter space on the BPT diagram where GN-z11 resides. This implies that flexible nebular emission modeling would be valuable for both describing and interpreting chemical evolution given more informative (higher resolution) data.

To conclude, the new redshift frontier established by JWST pushes stellar population models into new, exciting, and largely uncalibrated regimes. This work represents a first step toward a comprehensive error budget in inferred galaxy properties from SED modeling, paving the way for an accurate accounting of the intricate processes governing galaxy formation and evolution.

## Acknowledgments

We thank the anonymous referee for the helpful comments. We are also grateful to Ben Johnson for his help in developing the Cloudy emulator. This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452. B.W. and J.L. acknowledge support from JWST-GO-02561.022-A. H.A. is supported by CNES, focused on the JWST mission. H.A. acknowledges support from the Program National Cosmology and Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. P.D. acknowledges support from the Dutch Research Council (NWO) through the award of the VIDI grant 016.VIDI.189.162 ("ODIN") and the European Commission's and University of Groningen's CO-FUND Rosalind Franklin program. L.J.F acknowledges support by grant No. 2020750 from the United States-Israel Binational Science Foundation (BSF) and grant No. 2109066 from the United States National Science Foundation (NSF). J.E.G. acknowledges NSF/AAG grant No. 1007094. The work of C. C.W. is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy, Inc. (AURA) under a cooperative agreement with the NSF.

This work is based in part on observations made with the NASA/ESA/CSA JWST. The data were obtained from the

Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by AURA, under NASA contract NAS 5-03127 for JWST. These observations are associated with JWST-GO-2561, JWST-ERS-1324, and JWST-DD-2756. Support for program JWST-GO-2561 was provided by NASA through a grant from the Space Telescope Science Institute under NASA contract NAS 5-26555. This research is also based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute under NASA contract NAS 526555. These observations are associated with programs HST-GO-11689, HST-GO-13386, HST-GO/DD-13495, HST-GO-13389, HST-GO-15117, and HST-GO/DD-17231. The specific observations analyzed can be accessed via 10.17909/nftp-e621, and the SPS catalog corresponds to the following copy deposited at zenodo: 10.5281/zenodo.10223792.

Computations for this research were performed on the Pennsylvania State University's Institute for Computational and Data Sciences' Roar supercomputer. This publication made use of the NASA Astrophysical Data System for bibliographic information.

*Facility:* HST (ACS, WFC3), JWST (NIRCam, NIRSpec).

*Software:* Astropy (Astropy Collaboration et al. 2013, 2018, 2022), Dynasty (Speagle 2020), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), Prospector (Johnson et al. 2021), SciPy (Virtanen et al. 2020)

## Appendix

We supplement further details on the $z > 9$ sample of this paper in this appendix. Figure 9 shows the new candidates, whereas Figure 10 shows a possible data artifact. Additionally, we list the three objects with $z_{spec} > 9$ that are not included in the sample of this paper in Table 2. ID 11701 only has data in six NIRCam filters, and ID 21347 is identified as a broad-line AGN (Kokorev et al. 2023). Both sources show uncommon spectral signatures, and it is thus understandable that our photometric redshifts are inaccurate. ID 31955 has a likely spurious detection in a dropout band, albeit with large uncertainty, which skews the posterior mass in redshift to lower values.

<table><caption>**Table 2**
Discrepant Photometric and Spectroscopic Redshifts</caption><thead><tr><th>ID</th><th>ID<sub>MSA</sub></th><th>R.A.</th><th>decl.</th><th>z<sub>spec</sub><sup>a</sup></th><th>z<sub>phot, ml</sub></th><th>z<sub>phot, med</sub></th><th>Notes</th></tr></thead><tbody><tr><td>11701</td><td>10646</td><td>3.63696</td><td>-30.40636</td><td>8.510</td><td>0.275</td><td>0.272<sup>+0.012</sup><sub>-0.012</sub></td><td>Only have data available in six JWST filter bands; will be presented in J. Weaver et al. (2024, in preparation).</td></tr><tr><td>21347</td><td>20466</td><td>3.64041</td><td>-30.38644</td><td>8.500</td><td>1.305</td><td>1.309<sup>+0.132</sup><sub>-0.098</sub></td><td>Broad-line AGN (Kokorev et al. 2023)</td></tr><tr><td>31955</td><td>31028</td><td>3.54417</td><td>-30.37032</td><td>9.740</td><td>0.522</td><td>0.567<sup>+0.132</sup><sub>-0.094</sub></td><td>Noisy photometric measurement in F814W</td></tr></tbody></table>

**Note.**

<footnote>a NIRSpec/Prism confirmation conducted during the second phase of the UNCOVER survey (Fujimoto et al. 2023; S. Price et al. 2024, in preparation).</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** New candidates in this paper. For each candidate, we show cutouts including two HST bands (F606W, and F814W) and three JWST bands (F150W, F277W, and F444W), observed photometry in black, best-fit model in red, and the marginalized probability distribution of redshift posteriors. The black circle on each cutout indicates an aperture size of 0.′′32, within which the photometry is extracted.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** A possible artifact. ID 42658 only appears in the first epoch, and falls in the detector gap in the second epoch.</figcaption></figure>

ORCID iDs

Bingjie Wang (王冰洁) https://orcid.org/0000-0001-9269-5046

Joel Leja https://orcid.org/0000-0001-6755-1315

Hakim Atek https://orcid.org/0000-0002-7570-0824

Ivo Labbé https://orcid.org/0000-0002-2057-5376

Yijia Li (李轶佳) https://orcid.org/0000-0002-0682-3310

Rachel Bezanson https://orcid.org/0000-0001-5063-8254

Gabriel Brammer https://orcid.org/0000-0003-2680-005X

Sam E. Cutler https://orcid.org/0000-0002-7031-2865

Pratika Dayal https://orcid.org/0000-0001-8460-1564

Lukas J. Furtak https://orcid.org/0000-0001-6278-032X

Jenny E. Greene https://orcid.org/0000-0002-5612-3427

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

Vasily Kokorev https://orcid.org/0000-0002-5588-9156

Richard Pan https://orcid.org/0000-0002-9651-5716

Sedona H. Price https://orcid.org/0000-0002-0108-4176

Katherine A. Suess https://orcid.org/0000-0002-1714-1905

John R. Weaver https://orcid.org/0000-0003-1614-196X

Katherine E. Whitaker https://orcid.org/0000-0001-7160-3632

Christina C. Williams https://orcid.org/0000-0003-2919-7495

References

Alsing, J., Peiris, H., Mortlock, D., Leja, J., & Leistedt, B. 2023, ApJS, 264, 29

Anders, P., & Fritze-v. Alvensleben, U. 2003, A&A, 401, 1063

Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, Natur, 622, 707

Astropy Collaboration, Price-Whelan, A. M., & Lim, P. L. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., & Tollerud, E. J. 2013, A&A, 558, A33

Atek, H., Chemerynska, I., Wang, B., et al. 2023a, MNRAS, 524, 5486

Atek, H., Shuntov, M., Furtak, L. J., et al. 2023b, MNRAS, 519, 1201

Bakx, T. J. L. C., Zavala, J. A., Mitsuhashi, I., et al. 2023, MNRAS, 519, 5076

Bastian, N., Covey, K. R., & Meyer, M. R. 2010, ARA&A, 48, 339

Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143

Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379

Bezanson, R., Labbe, I., Whitaker, K. E., et al. 2022, arXiv:2212.04026

Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036

Boylan-Kolchin, M. 2023, NatAs, 7, 731

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

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

Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2012, Natur, 484, 485

Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2013, MNRAS, 432, 1862

Caputi, K. I., Lagache, G., Yan, L., et al. 2007, ApJ, 660, 97

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

Casey, C. M., Akins, H. B., Shuntov, M., et al. 2023, arXiv:2308.10932

Castellano, M., Fontana, A., Treu, T., et al. 2022, ApJL, 938, L15

Castellano, M., Fontana, A., Treu, T., et al. 2023, ApJL, 948, L14
Chabrier, G. 2003, PASP, 115, 763

Chabrier, G., Hennebelle, P., & Charlot, S. 2014, ApJ, 796, 75

Charlot, S., & Longhetti, M. 2001, MNRAS, 323, 349

Chatzikos, M., Bianchi, S., Camilloni, F., et al. 2023, RMxAAG, 59, 327

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

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Chon, S., Ono, H., Omukai, K., & Schneider, R. 2022, MNRAS, 514, 4639

Ciesla, L., Elbaz, D., Ilbert, O., et al. 2023, arXiv:2309.15720

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

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

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

Conroy, C., & van Dokkum, P. G. 2012, ApJ, 760, 71

Conroy, C., van Dokkum, P. G., & Villaume, A. 2017, ApJ, 837, 166

Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839

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

Dome, T., Tacchella, S., Fialkov, A., et al. 2024, MNRAS, 527, 2139

Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011

Dopita, M. A., Kewley, L. J., Heisler, C. A., & Sutherland, R. S. 2000, ApJ, 542, 224

Dotter, A. 2016, ApJS, 222, 8

D'Silva, J. C. J., Driver, S. P., Lagos, C. D. P., et al. 2023, ApJL, 959, L18

Duan, Q., Conselice, C. J., Li, Q., et al. 2023, arXiv:2309.14961

Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv:2306.02465

Endsley, R., Stark, D. P., Whitler, L., et al. 2023, MNRAS, 524, 2312

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

Ferrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, 522, 3986

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

Fioc, M., & Rocca-Volmerange, B. 1999, arXiv:astro-ph/9912179

Fujimoto, S., Wang, B., Weaver, J., et al. 2023, arXiv:2308.11609

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

Furtak, L. J., Atek, H., Lehnert, M. D., Chevallard, J., & Charlot, S. 2021, MNRAS, 501, 1568

Furtak, L. J., Zitrin, A., Weaver, J. R., et al. 2023, MNRAS, 523, 4568

Geha, M., Brown, T. M., Tumlinson, J., et al. 2013, ApJ, 771, 29

Gong, Y., Yue, B., Cao, Y., & Chen, X. 2023, ApJ, 947, 28

Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96

Gronke, M., Dijkstra, M., McCourt, M., & Oh, S. P. 2017, A&A, 607, A71

Groves, B. A., Dopita, M. A., & Sutherland, R. S. 2004, ApJS, 153, 9

Gunawardhana, M. L. P., Hopkins, A. M., Sharp, R. G., et al. 2011, MNRAS, 415, 1647

Guo, Y., Rafelski, M., Faber, S. M., et al. 2016, ApJ, 833, 37

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

Hahn, C., & Melchior, P. 2022, ApJ, 938, 11

Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19

Hunter, J. D. 2007, CSE, 9, 90

Jeřábková, T., Hasani Zonoozi, A., Kroupa, P., et al. 2018, A&A, 620, A39

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

Katz, H., Rosdahl, J., Kimm, T., et al. 2022, MNRAS, 510, 5603

Katz, H., Saxena, A., Cameron, A. J., et al. 2023, MNRAS, 518, 592

Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055

Kokorev, V., Fujimoto, S., Labbe, I., et al. 2023, ApJL, 957, L7

Kroupa, P. 2001, MNRAS, 322, 231

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

Legrand, L., Hutter, A., Dayal, P., et al. 2022, MNRAS, 509, 595

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

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

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

Leja, J., Johnson, B. D., Conroy, C., & van Dokkum, P. 2018, ApJ, 854, 62

Leja, J., Johnson, B. D., Conroy, C., et al. 2019b, ApJ, 877, 140

Leja, J., Speagle, J. S., Johnson, B. D., et al. 2020, ApJ, 893, 111

Li, J., Melchior, P., Hahn, C., & Huang, S. 2024, AJ, 167, 16

Li, Y., & Leja, J. 2022, ApJ, 940, 88

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

Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97

Maiolino, R., Scholtz, J., Witstok, J., et al. 2023, arXiv:2305.12492

Martín-Navarro, I., La Barbera, F., Vazdekis, A., Falcón-Barroso, J., & Ferreras, I. 2015, MNRAS, 447, 1033

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

Mathews, E. P., Leja, J., Speagle, J. S., et al. 2023, ApJ, 954, 132

Mauerhofer, V., & Dayal, P. 2023, MNRAS, 526, 2196

McKinney, J., Finnerty, L., Casey, C. M., et al. 2023, ApJL, 946, L39

McWilliam, A., Wallerstein, G., & Mottini, M. 2013, ApJ, 778, 149

Michel-Dansac, L., Blaizot, J., Garel, T., et al. 2020, A&A, 635, A154

Mollá, M., García-Vargas, M. L., & Bressan, A. 2009, MNRAS, 398, 451

Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022, ApJL, 940, L14

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

Nenkova, M., Sirocky, M. M., Ivezić, Ž., & Elitzur, M. 2008a, ApJ, 685, 147

Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezić, Ž., & Elitzur, M. 2008b, ApJ, 685, 160

Pacifici, C., da Cunha, E., Charlot, S., et al. 2015, MNRAS, 447, 786

Pacifici, C., Iyer, K. G., Mobasher, B., et al. 2023, ApJ, 944, L14

Pacucci, F., Dayal, P., Harikane, Y., Inoue, A. K., & Loeb, A. 2022, MNRAS, 514, L6

Padmanabhan, H., & Loeb, A. 2023, ApJL, 953, L4

Pallottini, A., & Ferrara, A. 2023, A&A, 677, L4

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

Parashari, P., & Laha, R. 2023, MNRAS, 526, L63

Pérez-González, P. G., Costantin, L., Langeroodi, D., et al. 2023, ApJL, 951, L1

Posacki, S., Cappellari, M., Treu, T., Pellegrini, S., & Ciotti, L. 2015, MNRAS, 446, 493

Reines, A. E., Nidever, D. L., Whelan, D. G., & Johnson, K. E. 2010, ApJ, 708, 26

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

Sabti, N., Muñoz, J. B., & Kamionkowski, M. 2024, PhRvL, 132, 061002

Salpeter, E. E. 1955, ApJ, 121, 161

Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, MNRAS, 371, 703

Santini, P., Ferguson, H. C., Fontana, A., et al. 2015, ApJ, 801, 97

Scholtz, J., Witten, C., Laporte, N., et al. 2023, arXiv:2306.09142

Serjeant, S., & Bakx, T. J. L. C. 2023, NatAs, 7, 1143

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

Smit, R., Bouwens, R. J., Franx, M., et al. 2015, ApJ, 801, 122

Smit, R., Bouwens, R. J., Labbé, I., et al. 2014, ApJ, 784, 58

Smith, R. J. 2020, ARA&A, 58, 577

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

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

Steinhardt, C. L., Jauzac, M., Acebron, A., et al. 2020, ApJS, 247, 64

Steinhardt, C. L., Kokorev, V., Rusakov, V., Garcia, E., & Sneppen, A. 2023, ApJL, 951, L40

Suess, K. A., Leja, J., Johnson, B. D., et al. 2022, ApJ, 935, 146

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

Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 458, 242

Tacchella, S., Finkelstein, S. L., Bagley, M., et al. 2022, ApJ, 927, 170

Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673

Treu, T., Auger, M. W., Koopmans, L. V. E., et al. 2010, ApJ, 709, 1195

Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110

Trotta, R. 2008, ConPh, 49, 71

Tsuna, D., Nakazato, Y., & Hartwig, T. 2023, MNRAS, 526, 4801

van Dokkum, P., Brammer, G., Wang, B., Leja, J., & Conroy, C. 2024, NatAs, 8, 119

van Dokkum, P., Conroy, C., Villaume, A., Brodie, J., & Romanowsky, A. J. 2017, ApJ, 841, 68

van Dokkum, P. G., & Conroy, C. 2012, ApJ, 760, 70

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Wang, B., Fujimoto, S., Labbé, I., et al. 2023a, ApJL, 957, L34

Wang, B., Leja, J., Bezanson, R., et al. 2023b, ApJL, 944, L58

Wang, B., Leja, J., Labbé, I., et al. 2024, ApJS, 270, 12

Wang, B., Leja, J., Villar, V. A., & Speagle, J. S. 2023c, ApJL, 952, L10

Weaver, J. R., Cutler, S. E., Pan, R., et al. 2024, ApJS, 270, 7

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

Woodrum, C., Rieke, M., Ji, Z., et al. 2023, arXiv:2310.18464

Xiao, M., Oesch, P., Elbaz, D., et al. 2023, arXiv:2309.02492

Yan, Z., Jerabkova, T., Kroupa, P., & Vazdekis, A. 2019, A&A, 629, A93

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

Zavala, J. A., Buat, V., Casey, C. M., et al. 2023, ApJL, 943, L9

Ziparo, F., Ferrara, A., Sommovigo, L., & Kohandel, M. 2023, MNRAS, 520, 2445