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

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). ¹⁷ 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.

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 dynesty (Speagle 2020) with a speed up in model generation enabled by a neural net emulator dubbed parrot (Mathews et al. 2023).

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 ≥9.0
  • (b)  
    0.1th percentile of the redshift posteriors ≥3.0
• 2.  
  Maximum-likelihood redshift ≥8.0
• 3.  
  ${\chi }_{\mathrm{red}}^{2}={\chi }^{2}/{n}_{\mathrm{bands}}\lt 5$
• 4.  
  Number of NIRCam filters >6
• 5.  
  S/N in photometry

  • (a)  
    F444W ≥2.0
  • (b)  
    F410M ≥2.0 if available
  • (c)  
    F277W ≥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 }_{\mathrm{red}}^{2}\lt 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.

Table 1.  z_phot > 9 Candidates Studied in This Paper

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

Notes.

^a Posterior moments of photometric redshifts and stellar masses from Prospector fits. ^b 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). ^c 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). ^d Possible data quality issue; see the Appendix for details.

Download table as:  ASCIITypeset image

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 ≳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.

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

$$\begin{eqnarray}&&\xi \propto {m}^{-x}.\end{eqnarray} \tag{ 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_⊙ 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_⊙ (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 χ² 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.

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 (Behroozi 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.

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).

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 ≳0.3 dex; increasing m_l from 0.08 M_⊙ to 0.5 or 1 M_⊙ further enhances the systematic shifts by ≳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 = 1 M_⊙, 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).

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

The bursty SFH prior of this paper causes negligible systematic changes in the inferred stellar masses. However, it tends to increase the uncertainty in the inferred SFR averaged over the most recent 10 Myr, SFR₁₀, 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₁₀₀, 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.

Zoom In Zoom Out Reset image size

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 ∼ 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.

The more flexible nebular emission modeling, enabled by the Cloudy emulator, Cue, causes differences in the inferred masses and SFRs that are within 1σ 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 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.

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 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.

Zoom In Zoom Out Reset image size

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_⊙, 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_⊙, 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_⊙ (1 M_⊙).

Zoom In Zoom Out Reset image size

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 (∼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 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.

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.