Inferred galaxy properties during Cosmic Dawn from early JWST photometry results

ABSTRACT

Early photometric results from JWST have revealed a number of galaxy candidates above redshift 10. The initial estimates of inferred stellar masses and the associated cosmic star formation rates are above most theoretical model predictions up to a factor of 20 in the most extreme cases, while this has been moderated after the recalibration of NIRCam and subsequent spectroscopic detections. Using these recent JWST observations, we use galaxy scaling relations from cosmological simulations to model the star formation history to very high redshifts, back to a starting halo mass of 10⁷ M_⊙, to infer the intrinsic properties of the JWST galaxies. Here, we explore the contribution of supermassive black holes, stellar binaries, and an excess of massive stars to the overall luminosity of high-redshift galaxies. Despite the addition of alternative components to the spectral energy distribution, we find stellar masses equal to or slightly higher than previous stellar mass estimates. Most galaxy spectra are dominated by the stellar component, and the exact choice for the stellar population model does not appear to make a major difference. We find that four of the 12 high-redshift galaxy candidates are best fit with a non-negligible active galactic nuclei component, but the evidence from the continuum alone is insufficient to confirm their existence. Upcoming spectroscopic observations of z > 10 galaxies will confirm the presence and nature of high-energy sources in the early Universe and will constrain their exact redshifts.

galaxies: formation, galaxies: high-redshift, galaxies: photometry, (galaxies:) quasars: supermassive black holes

1 INTRODUCTION

JWST has provided us with data from a time never seen before. We are now seeing galaxies, as they first form when the Universe was only a few hundred million years old. The conclusions that we draw as a community will continue to inform our models for many years to come, and thus, the early JWST results act as our initial conditions for this new era of astronomy. The first observations of these high-redshift galaxies are coming from two major surveys: the Cosmic Evolution Early Release Survey (CEERS; Finkelstein et al. 2017) and the Grism Lens-Amplified Survey from Space (GLASS; Treu et al. 2022). CEERS is a photometric and spectroscopic survey covering 100 square arcminutes of sky utilizing the NIRSpec, NIRCam, and MIRI instruments. One of CEERS’s specific goals is to discover galaxies between 9 < z < 13 and observe their spectra to constrain details on early galaxy formation. One of the first and more striking results coming from the CEERS survey is the detection of a galaxy at z = 11.8 (Finkelstein et al. 2022). The GLASS survey is aimed at answering questions related to reionization and baryon usage and recycling within faintly magnified galaxies by utilizing the NIRISS, NIRSpec, and the NIRCam instruments. Some of the first results in this survey showed lensed low-metallicity and star forming galaxies at z ≥ 7 (Roberts-Borsani et al. 2022) and higher redshift galaxies at z ∼ 9−15 (Castellano et al. 2022b), demonstrating the power of high-redshift observations that can be done by JWST.

These first galaxies began to form approximately 100 Myr after the big bang as small dark matter haloes (minihaloes) merged to form larger haloes, hierarchically building a home for future galaxies. While the term ‘first galaxy’ is not well defined, it is commonly agreed that the first galaxies were ones that were able to form stars and withstand their feedback (Bromm & Yoshida 2011), in order to continue to merge and grow. At early times, the first stars to form were metal-free [Population III (Pop III)] stars that are thought to have been very massive from the lack of efficient cooling (Yoshida et al. 2003). Simulations have shown that Pop III stars primarily formed in minihaloes with masses ∼10⁵⁻⁶ M_⊙ at z > 10 (Yoshida et al. 2003; Bromm & Yoshida 2011; Schauer et al. 2019; Skinner & Wise 2020). If we consider these dark matter minihaloes to be the hosts of the first galaxies, their initial formation would have occurred around z ∼ 30 and continued to grow and merge with one another until they reach the atomic cooling limit at M_halo ≃ 10⁸ M_⊙ (e.g. Wise & Abel 2007). These haloes are then able to continuously form stars in a less bursty fashion, and if massive enough, can withstand the feedback. These galaxies are ones that we may be more familiar with in the present day. This process of going from small dark matter haloes, to the first generation of bursty star formation, to atomically cooling haloes drove the early stages of reionization (Wise et al. 2014). Observing these first galaxies with JWST is a massive step forward in understanding the importance of these sources in relation to reionization.

Still, at such early times, far from the views of any telescopes that have looked to the cosmos so far, simulations remain the only tool informing our understanding of galaxy formation in the very early Universe (see Somerville & Davé 2015; Wechsler & Tinker 2018 for reviews). The Renaissance simulations (Xu, Wise & Norman 2013; Chen et al. 2014; Xu et al. 2014; O’Shea et al. 2015) are a suite of hydrodynamical adaptive mesh refinement zoom-in simulations focusing on different regions of the universe; namely overdense, normal, and low-density regions. These simulations were initially used to study Pop III star formation, but have since been used to study reionization, first galaxy formation, and black hole (BH) growth. In particular, Chen et al. (2014) used the Renaissance simulations to provide various scaling relations for galaxies forming at z ≥ 15. They found that haloes below the atomic cooling threshold show bursty star formation as opposed to haloes above the atomic cooling threshold, where the star formation becomes more efficient.

The SPHINX simulations are another suite of cosmological hydrodynamical adaptive mesh refinement simulations studying the effects of various astrophysical processes on reionization. In their first paper, Rosdahl et al. (2018) studied the effects of binary systems on reionization. They found that binary stars are necessary to include if the simulation volume is to be reionized by z ∼ 6. Including binary stars leads to higher escape fractions as compared to the single-star counterparts. The Obelisk simulation followed up on the SPHINX simulations, finding that stellar radiation predominately drove reionization (Trebitsch et al. 2021). More recently, Katz et al. (2022a) used the SPHINX simulations to study the C and O abundances in galaxies during reionization and found that the stellar population within galaxies at z > 6 may need to be under the influence of a top-heavy initial mass function (IMF) in order to reproduce the observed abundances. At these high redshifts, this is somewhat expected since the first generation of stars are expected to come from a more top-heavy IMF. FLARES (First Light and Reionization Epoch Simulations) is also a cosmological suite following galaxy formation in the early universe based on the EAGLE simulations, a smoothed particle hydrodynamic (SPH) simulation suite (Schaye et al. 2015). Vijayan et al. (2021) used FLARES to predict photometric data for high-redshift galaxies, in preparation for JWST. They presented a ultraviolet (UV) luminosity function (LF) that matched well with the observed data at that time. More recently, Wilkins et al. (2023) used FLARES to study the stellar histories of galaxies in the early universe, within the redshift range that JWST aims to cover. They found a stellar mass–metallicity relationship at very high redshifts, which can be probed and verified by results from JWST, and that the environment does not play a huge role in affecting the stellar masses, metallicities, star formation rates (SFRs), or ages. As more data comes from JWST, the results from such simulations will be further constrained and will help us understand galaxy formation at these early times.

As the first galaxies begin to grow into the large galaxies we see today, so do the BHs that reside in their hosts and over time, will migrate towards their centres. Beginning in the 1960s, quasars were first discovered as radio bright sources in the sky (Schmidt 1963), only later to be discovered and labelled as quasars, or supermassive black holes (SMBHs) growing rapidly from their surroundings. Since 2001, nearly 100 quasars have been discovered with masses of |$\lower0.3em\hbox{$\, \buildrel\gt \over \sim \, $}10^{9} { {\rm M_\odot }}$| at z > 6 (see Inayoshi et al. 2022, for a collection of these sources). The existence of these SMBHs at such early times requires an explanation for how they grew so quickly in such a short period of time. There are many possible origins of these BHs, including stellar mass remnants from Pop III stars, intermediate-mass BHs formed from stellar collisions, or direct collapse BHs formed from collapsing gas clouds (Valiante et al. 2017). From here, models of active galactic nuclei (AGNs) have helped answer some questions regarding how efficient accretion can be at growing these BHs to be very massive. Assuming that the galaxies that are observed by JWST are some of the brightest and thus most massive galaxies, they are likely to have an SMBH at their centres.

2 THE FIRST SIX MONTHS OF JWST RESULTS

As soon as data from JWST became publicly available in early July 2022, preprints started to be posted to the arXiv. Bright, high redshift galaxies were being reported (e.g. Finkelstein et al. 2022; Labbe et al. 2022; Castellano et al. 2022b) and along with it, a growing concern that these galaxies were not abiding by the laws of Lambda cold dark matter (ΛCDM) cosmology. But with new telescopes comes the natural growing pains of calibration and data re-analysis. Calibrating the instruments aboard JWST is a critical step to ensure that we are able to interpret the data as accurately as possible (e.g. Bagley et al. 2022; Rigby et al. 2022). At the time of writing, calibration is still ongoing which will require the reprocessing of already published data. An example of the effects of recalibration and reanalysis comes from Adams et al. (2023). Using post-launch calibrations, they search for galaxies at z > 9 and compare with results from other studies. They find that compared to the studies using the early July 2022 calibration, a large percentage of the reportedly high-redshift galaxies are actually at much lower redshifts (see their table 4 and references therein), and there is not a lot of overlap of sources between studies of the same fields. The reasons for these issues have not been fully explored, but they do find that some sources are affected by the post-flight calibration, in that a lower redshift is obtained when using the new calibrations. Boyer et al. (2022) found that there is a large flux offset in the NIRCam filters, resulting in magnitude differences of ∼0.01−0.2 mag, and reported new zero-points for the NIRCam filters. Calibration will take place throughout Cycle 1 of JWST, and will continue to be updated, requiring reanalysis of already published data. While calibration is an important issue, improved reanalysis of these objects as time goes on may lead to different results. For example, Finkelstein et al. (2022) initially reported finding a highly star-forming galaxy at z ∼ 14 with a high stellar mass. After improving their astrometry methods, certain filters were better aligned with each other and with the aperture, resulting in an overall brighter flux. This resulted in the galaxy’s redshift being decreased to z ∼ 11.8. While this is still a very high redshift galaxy, this is an example of how initial measurements of these early galaxies will likely change with the improved calibration and analysis of JWST data.

Alongside the rush of reports on these high-redshift galaxies, theoretical explorations on the abundance and limits of these galaxies had begun to be reported as well. Boylan-Kolchin (2022) presented a straightforward calculation of the comoving number density of haloes and the maximum amount of stellar mass contained within a halo as a function of redshift assuming maximum star formation efficiency (SFE). They found that some early measurements of JWST galaxies with high stellar masses are found in a volume much smaller than expected, and it appears that there are two galaxies whose stellar mass density lies above the comoving stellar mass density assuming maximum SFE, i.e. their stellar mass is larger than the amount of baryons available to their host haloes. The analysis presented in Boylan-Kolchin (2022) provided some of the first examples of the possible tension with ΛCDM cosmology. Relieving this tension, Mason, Trenti & Treu (2023) investigated the upper limit of the UV LF by modelling the SFR and including a simple stellar spectrum as a function of time. They found that when the SFE is at a maximum, the upper limit to the UV LF is about four orders of magnitude higher than what is currently detected, showing that these detections are still in the realm of ΛCDM cosmology. They also found that when the SFE is significantly decreased to a more realistic value and is calculated as a function of the halo mass, and if dust attenuation is negligible at high redshifts, the UV LF more closely matches currently observed values. Mason et al. (2023) go on to show that younger star-forming galaxies that have rapidly evolved are the main galaxies that are detectable due to their increased UV magnitude, implying that we are selectively only seeing the brightest galaxies, and they may not represent the entire population of galaxies at that time. Another possible explanation comes from Ferrara, Pallottini & Dayal (2022), who constructed a minimal physical model of the UV LF as a function of time. They modelled the SFR using the classical Schmidt-type expression and altered the SFE by including feedback from supernova (SN). They then fit their UV LF to the ALMA REBELS survey to calibrate for dust attenuation, and find that some high redshift bright detections at z > 11 are best described with a UV LF that does not include dust attenuation, supporting results from Mason et al. (2023). Inayoshi et al. (2022) also determined the upper bound of the LF at three different redshift ranges and compare to the LF with varying SFEs. They found that current detections fit well with SFEs that closely match SFEs in starburst galaxies. They also inspected the product of the SFE and the UV radiative efficiency, providing another constraint. Meeting this constraint could mean a combination of periods of very efficient metal-poor and relatively efficient metal-free star formation. They concluded that weakly accreting BHs and quasars could also contribute to the UV continuum. Results from Inayoshi et al. (2022) agree with results from Mason et al. (2023) in that these high-redshift massive galaxies do not violate ΛCDM cosmology, and they also provide us with constraints on their stellar population.

Groups have also been reanalysing the JWST data to provide more insight and tools for future JWST observations. Furtak et al. (2023) detailed their spectral energy distribution (SED) fitting process for a number of weakly lensed galaxies, and provided their galactic best-fitting parameters. They found that the galaxies at z ≥ 10 are young and highly star forming, which aligns with results from Mason et al. (2023). They calculated the mass-to-light ratios of each galaxy and found no redshift-dependent evolution of this relation. Importantly, the galaxies they analyse are all weakly magnified by gravitational lensing, and thus they are only looking at the brightest galaxies in the UV, echoing concerns about selection bias from Mason et al. (2023). Higher magnification is needed to observe galaxies on the fainter end of the UV spectrum. Furtak et al. (2023) concluded that spectra for these galaxies are needed in order to confirm their redshift, as they find that for some galaxies that do not have strong Balmer breaks, a low-redshift SED can also describe the data. Finally, they say that their data remain within the upper limits ΛCDM cosmology.

Some results are showing the importance of considering other sources contributing to the SED. Brinchmann (2022) analysed JWST NIRSpec data to compare high-redshift sources with local analogues. They found that a few of the high-redshift sources they analyse are likely affected by an AGN due to emission lines associated with AGN activity at lower redshifts. This introduces the important idea that the UV flux may be coming from more than just a stellar population, and additional sources should be included. Steinhardt et al. (2022) generated new photometric templates in order to more accurately fit the JWST photometry at high redshifts, rather than using templates based on IMFs calibrated to the local Universe. This is important due to a major environmental difference at early times, namely, the increased temperature of the cosmic microwave background (CMB), resulting in a bottom-light IMF (Smith et al. 2009). They find that a change in the IMF produces a 0.2–0.3 dex increase in stellar mass and that a change in the gas temperature produces a much stronger effect. Since their work is based on gas temperature being heated by the CMB, this is a lower bound on the gas temperature, and in reality, the gas temperature could be much higher, resulting in an even bottom-lighter IMF. They state that this could mean that fitting high redshifts galaxies with the local IMF could be overestimating the stellar mass by ≥1 dex. Steinhardt et al. (2022) assert that the ‘early galaxy problem’ is solved by using fits that account for the CMB heating gas at higher redshifts. Volonteri, Habouzit & Colpi (2022) found that the unobscured AGN emission in these JWST high-z galaxy candidates is generally fainter than the stellar component and will be difficult to detect in colour–colour selections; however, overmassive SMBHs may be detectable that could point towards it being seeded by a massive BH.

Given the theory discussions so far, the ‘early galaxy problem’ may not be as big a problem as initially thought. A top-heavy IMF may help explain the higher UV magnitudes and may be a more appropriate choice to model high-redshift galaxies, as explained by Steinhardt et al. (2022). Including an AGN and more luminous sources may be another missing piece to the puzzle, since there have been some emission lines detected that are indicative of an AGN in some galaxies so far (Brinchmann 2022). Adjusting our models for each of these pieces may provide explanations to the high stellar masses found in the high-redshift galaxies. To tackle this problem, we investigate how an AGN, binary stars, and top-heavy IMFs may affect the UV magnitude of a set of JWST galaxies by modelling the halo mass and the star formation history. From the star formation history, we determine the stellar spectrum given an IMF choice and add to that to an AGN spectrum. We vary a set of free parameters in order to find the set that best fit the observed rest-frame UV and optical magnitudes. We describe our model in Section 3, detailing the star formation history in Section 3.1, the stellar spectrum in Section 3.2 and the AGN spectrum in Section 3.3. Details about dust attenuation, computing the magnitudes, and tests are described in Sections 3.4, 3.5, and 3.6, respectively. We vary our model using a Monte Carlo Markov Chain (MCMC), described in Section 3.7. Our results are described and discussed in Sections 4 and 5, respectively. We conclude with Section 6.

3 METHODS

To explore how the JWST photometry is affected by different properties of a galaxy, we create a model made of three discrete parts, varying the star formation history, the stellar population, and the central SMBH. We then use emcee (Foreman-Mackey et al. 2013), an MCMC ensemble sampler, to vary certain free parameters in order to determine the most probable set of characteristics that would explain the observed photometry. We run our analysis pipeline on each of the JWST galaxies listed in Table 1, coming from Adams et al. (2023, hereafter A22) and Harikane et al. (2022, hereafter H22), where we only consider the candidates in the latter that are also selected by other works. The complete analysis pipeline can be seen in Fig. 1 and the model choices and free parameters are listed in Tables 2 and 3, respectively. Throughout our calculations, we use the cosmological parameters from Planck Collaboration VI (2020): Ω_M = 0.3153, Ω_Λ = 0.6847, Ω_b = 0.04923, σ₈ = 0.8111, and h = 0.6736, where the parameters are their usual definitions.

Figure 1.

Flowchart of the analysis pipeline. The red parameters indicate individual model choices, while the blue parameters indicate parameters that are continuously varied by emcee. To generate a model SED, model A or B is chosen to model the star formation history. The outputs are then fed to the stellar and AGN spectrum models. The stellar spectrum is controlled by choosing one IMF and one stellar population. The AGN spectrum is controlled by two parameters varied by emcee. Once the model SED (the grey line shows the total SED, the pink line the AGN component, and the blue line the stellar component) is generated and attenuation is applied, the AB magnitudes and χ² are calculated. emcee then runs the analysis again for a set of red parameters and varies the blue parameters.

Table 1.

Sample JWST galaxies.

ID	z	log10(M⋆/M⊙)	Other detections (alternative ID)
Adams et al. (2023, A22)
1514	9.85	9.8	–
1696	9.59	8.7	At23 (SMACS-z10c)
2462	9.50	9.5	At23 (SMACS-z10b)
2779	9.51	8.7	–
6115	10.94	8.4	–
6878	9.59	9.3	At23 (SMACS-z10a)
10234	11.42	7.8	–
Harikane et al. (2022, H22)
GL-z9 − 1	10.49	|$9.04^{+0.61}_{-0.74}$|	C22, C22B (GHZ1), D23 (17487), N22 (GLz11)
GL-z9 − 2	10.46	|$8.64^{+0.16}_{-1.04}$|	C22, C22B (GHZ4)
GL-z12 − 1	12.28	|$8.36^{+0.90}_{-0.28}$|	C22 (GHZ2), D23 (1698), N22 (GLz13)
CR2-z12 − 1	11.63	|$7.88^{+0.95}_{-0.26}$|	F22 (Maisie’s Galaxy)
CR2-z16 − 1	16.25	|$9.20^{+1.06}_{-0.72}$|	D23 (93316)

Note. Numeric and alphanumeric IDs are taken from A22 and H22, respectively. Photometric redshifts and stellar masses are shown from these two works with the exception of GL-z9 − 2, where the stellar mass is from C22.

Sources: At23 (Atek et al. 2023), C22 (Castellano et al. 2022a), C22B (Castellano et al. 2022b), D23 (Donnan et al. 2023), F22 (Finkelstein et al. 2022), N22 (Naidu et al. 2022).

Table 2.

Fixed model choices.

	Name	Explanation
Growth Model	A	Press-Schechter
	B	Exponential Growth
IMF	Chabrier	α1 = exp cut-off; α2 = −2.3
		M1 = 1.0 M⊙; Mmax = 100 M⊙
	Standard	α1 = −1.30; α2 = −2.35
		M1 = 0.5 M⊙; Mmax = 100 M⊙
	Top heavy	σ = 1 M⊙;  Mc = 10 M⊙
		Mmin = 1 M⊙; Mmax = 500 M⊙
Stellar Population	Binary	Binary & non-binary population spectra
	Non-binary	Only single stellar population spectra

Note. A single model consists of a choice of growth model, IMF, and stellar population. Note that the Chabrier and Standard IMF models come from BPASS (Eldridge et al. 2017) and the Top-Heavy IMF model comes from Yggdrasil (Zackrisson et al. 2011). emcee is run separately for each fixed model.

Table 3.

Free parameter choices.

	Definition	Range
MHalo	Mass of the halo in M⊙	9 < log10(M/M⊙) < 12
α	Controls halo growth rate (Wechsler et al. 2002)	0.4 < α < 0.8
δBH	Alters AGN slope (Yang et al. 2022)	−1 ≤ δBH ≤ 1
fBH	Controls BH mass in terms of stellar mass	−4 ≤ log10(fBH) ≤ 0
τV	Optical depth in the visual band	0 < τV < 2

Note. Parameter choices. These values are varied by emcee to map out parameter space. The α parameter is only required for Model B.

3.1 Star formation history

In order to produce a model SED that best represents the observed photometry for a given JWST galaxy, the history of the galaxy must be determined. This includes the halo mass, stellar mass, stellar ages, and stellar metallicity throughout time. These components are inputs to the spectrum models in Sections 3.2 and 3.3. As described in Section 3.7, the halo mass at the observed redshift will be varied to determine the most probable set of parameters to fit the observed luminosity. The stellar mass is computed from the fits in Behroozi, Wechsler & Conroy (2013, their equation 3), given the host halo mass. While this equation is fit only for haloes between 0 < z < 8, Chen et al. (2014) found that these relations fit well with their results for z ≤ 15, thus we use the z = 15 fit for redshifts above that value. This fit comparison can be seen in Fig. 2. In order to better model haloes at higher masses, we extrapolate the low-mass power law to higher masses. We also limit the stellar mass to a maximum SFE of 0.03 (Pillepich et al. 2018), as shown in Fig. 2. This relation is also consistent with Riaz, Hartwig & Latif (2022), who determined the stellar mass–halo mass (SMHM) relation from a semi-analytic model of halo masses at high redshifts, including the contribution of Pop III stars.

Figure 2.

The SMHM relation implemented in our model compared to different sources. The solid lines labelled by a redshift come from the fit from Behroozi et al. (2013). The similarly coloured dashed lines are our linear fits to the low-mass power laws at the noted redshifts. They converge to a stellar mass of 0.03M_Halo, since we limit the stellar mass to a maximum SFE, indicated by the opaque pink line. The opaque black line shows the fit from Chen et al. (2014). Notice how this relation matches well with the Behroozi et al. (2013) fit at z = 15. The opaque blue line shows the relation from Riaz et al. (2022).

The stellar-mass–halo-mass relation is not well studied for massive haloes at such high redshifts. Behroozi et al. (2013) only modelled this relation for haloes from redshift 0 to 8. At redshift 8, the largest halo mass that they modelled was 10¹¹ M_⊙, which is reasonable given that a halo of that mass is very rare at z = 8. Although comparison with Chen et al. (2014) shows that extending the SMHM relation to redshift 15 is consistent and matches simulations very well at lower masses, the maximum halo mass in Chen et al. (2014) was 10⁹ M_⊙. Equation 3 in Behroozi et al. (2013) is not meant to be used for halo masses larger than 10¹¹ M_⊙ at redshift 8. This means the plateau seen in Fig. 2 is unrealistic. The plateau is an artefact of applying their equation beyond the intended range of halo mass and redshift. The turnover describes the transition to elliptical galaxies which occurs at lower redshifts. Therefore, we cannot simply use the Behroozi et al. (2013) equation as is for large halo masses. The SMHM relation is often modelled as a double power law (Yang et al. 2012; Moster, Naab & White 2013) so in the absence of a better model, we argue that the most sensible choice is to transition to a high-mass power law with a physically motivated maximum slope, determined by the maximum SFE. Our SFE limiting causes the low-mass SMHM relation to transition to a shallower slope that more accurately describes the high-mass regime.

We explore two models to generate a halo growth history based on the halo mass of the JWST galaxy calculated previously at the observed redshift. The first model comes from the ellipsoidal variant of Press–Schechter formalism (Sheth, Mo & Tormen 2001). Given the halo mass and observed redshift, we directly calculate the rarity of the halo (ν) and the halo mass history beginning at the observed redshift. The second model is an exponential growth history from Wechsler et al. (2002): M(z) = M₀e^−αz, where M₀ is the present-day halo mass and α is a free parameter controlling the rate of halo growth. The higher the α value, the more rapid the halo growth. For both models, the halo growth history is calculated back to when M_Halo = 10⁷ M_⊙ in steps of 5 Myr. Given the halo mass growth history, we now directly calculate the stellar mass history as given by Behroozi et al. (2013) and shown in Fig. 2. This allows us to calculate the SFR, stellar ages, and metallicities as a function of time. The SFR is calculated simply by taking the difference in the stellar mass between two time steps: SFR = (M_{⋆, i} − M_{⋆, i − 1})/(t_i − t_{i − 1}). The stellar age is simply the lookback time from the photometric redshift estimate. The stellar metallicities are calculated from Torrey et al. (2019, T19), where they describe the mass–metallicity relation from the IllustrisTNG simulation suite. They present the Log(O/H) + 12 galactic metallicity value for a series of stellar mass bins – log₁₀M_⋆/M_⊙ = [8, 8.5, 9, 9.5, 10, 10.5] – from 0 < z < 10 (see their fig. 7). We take the y-intercept points and interpolate between these lines for our given stellar mass and redshift. We then convert this into a total metallicity (Z_Tot) in units of Z_⊙. To calculate the stellar metallicity (Z_⋆) of a single stellar population (SSP) at each time-step, we multiply the change in total metallicity by the mass fraction of new stellar mass created: Z_⋆ = (Z_Tot,i − Z_Tot,i-1) × (M_⋆,i/(M_⋆,i − M_{⋆,i − 1})). Below 10⁸ and above 10^10.5 M_⊙, we use the minimum and maximum of the T19 relationship at 10⁸ and 10^10.5 M_⊙, respectively. While this relation is fit only to a particular mass range and up to z ∼ 10, it generally matches well with simulations. The T19 metallicity differs slightly from metallicity relations from other simulations for the same stellar mass, and redshift. Here, we present comparisons with three other simulations: Wilkins et al. (2023, W22), Abe et al. (2021, Abe21), and Jeon et al. (2015, J15). We note that the MZR we use from T19 is a galactic gas-phase MZR, not specifically a stellar MZR. However, we use this as an approximation for the stellar MZR. J15 and Abe21 also use a gas-phase metallicity; however, W22 use a stellar metallicity from a smoothing kernel applied to the star particles in their simulation. In each of these comparisons, we measure the difference in UV magnitude of an SSP at a given redshift, age, and stellar mass, only differing in the metallicity.

Compared to W22, with M_⋆ = 10⁹ M_⊙, z = 10, Age = 45 Myr, Z_Wilkins = 0.48 Z_⊙, Z_Torrey = 0.16 Z_⊙, using the T19 metallicity produces brighter UV flux by 0.15 mag.

In the following comparisons, the simulated stellar mass quoted did not form from an SSP at one age, however, since the authors do not present a full mass–metallicity relation as a function of redshift, the best we can do is compare the spectra of SSPs at various ages.

We compared to the fiducial simulation from Abe21, with M_⋆ = 7.5 × 10⁴ M_⊙ h⁻¹, z = 9, Z_Abe = 7.4 × 10⁻³Z_⊙, Z_Torrey = 0.12 Z_⊙, for SSPs with ages ranging from 10 to 200 Myr. Using the T19 metallicity resulted in a dimmer UV magnitude than when using the metallicity from Abe21. The range of the difference in magnitudes was 0.14–0.35 mag. The difference monotonically increased with increasing stellar age.

We performed the same comparison using parameters from J15. Namely, M_⋆ = 3 × 10⁴ M_⊙, z = 10.5, Z_Jeon = 10⁻⁴Z_⊙, Z_Torrey = 0.095 Z_⊙, for SSPs with ages ranging from 10 to 200 Myr. Using the much larger T19 metallicity resulted in a much brighter UV magnitude than when using the metallicity from J15. The range of the difference in magnitudes was 0.60–1.93 mag. Again, the difference monotonically increased with increasing stellar age.

This shows that although the T19 MZR shows different metallicities than other simulations, the change in the UV magnitude compared to W22 and Abe21 is small. However, the difference in metallicity between T19 and J15 is large and results a large difference in flux, especially for old stellar populations. However, given that T19 has the most robust MZR for high-redshift with stellar masses spanning many orders of magnitude, it is an ideal model for generating spectra for our SED fitting pipeline.

From each growth history model, the stellar masses, ages, and metallicities are passed to the next two components of our pipeline in order to generate a stellar and AGN spectrum.

3.2 Stellar spectrum

Given the stellar mass, age, and metallicity histories for each galaxy from the star formation history model, we calculate a stellar spectrum at the observed redshift by summing the stellar spectra of SSPs forming at 5 Myr intervals. We obtain the stellar spectra of each SSP from Binary Population and Spectral Synthesis (BPASS; Eldridge et al. 2017) and Yggdrasil (Zackrisson et al. 2011).

As the BPASS name suggests, it is used for simulating stellar populations as well as hosting a library of spectra generated from simulations using a set of IMFs and stellar populations at various stellar ages and metallicities. BPASS specializes in synthesizing binary stellar populations, creating spectra with similar binary system distributions compared to the local Universe. The Yggdrasil spectra instead use SSPs and focus on Pop III and low-metallicity stars for a better comparison with the observations of the first galaxies. The Yggdrasil SED data base models Pop III (characteristic masses of 100 M_⊙), Pop III.2 (characteristic masses of 10 M_⊙), and metal-poor Pop II stars.

In our parameter exploration, we use the Chabrier (2003) IMF and a ‘standard’ IMF from BPASS and a ‘Pop III.2’ (top-heavy) IMF from Yggdrasil. The SMHM relation used to generate our star formation history is based on the Chabrier IMF, but we explore the effect of different IMFs on the stellar spectrum here. We include the top-heavy IMF in our analysis since the first galaxies formed, while the CMB temperature was higher and the interstellar medium was predominantly metal-poor, cooling less efficiently than solar-metallicity gas and possibly reducing fragmentation into low-mass stars, leading to a larger number of more massive stars than seen in the local Universe.

The IMF functions modelled in BPASS follow a broken power law:

$$\begin{eqnarray} N(M \lt M_{\mathrm{max}}) \propto \int _{0.1}^{M_1} \left(\frac{M}{{ {\rm M_\odot }}}\right)^{\alpha _ 1} \, {\rm d}M \\ + \ M_{1}^{\alpha _1} \int _{M_{1}}^{M_{\mathrm{max}}} \left(\frac{M}{{ {\rm M_\odot }}} \right)^{\alpha _2} \, {\rm d}M. \end{eqnarray}$$

(1)

Here, each IMF modelled in BPASS varies the α parameters as well as M₁ and M_max that are listed in Table 2. The Chabrier IMF in BPASS is modelled with an exponential cut-off in the low-mass regime to recover the smooth transition from low- to high-mass stars that is typical of the Chabrier IMF. Our chosen ‘standard’ IMF will have a more discontinuous drop in N(M) from low- to high-mass stellar mass ranges.

The IMF function for Yggdrasil adopts a lognormal IMF taking the form

$$\begin{eqnarray} \mathrm{ln} \left(\frac{{\rm d}N}{{\rm d}\mathrm{ ln}(M)} \right) = A-\frac{1}{2\sigma ^2} \left[\mathrm{ln} \left(\frac{M}{M_\mathrm{c}} \right) \right]^2, \end{eqnarray}$$

(2)

where σ is the width of the distribution, M_c is the characteristic mass, and A is an arbitrary normalization (Tumlinson 2006), whose values are also listed in Table 2.

Each SED available in BPASS is organized into a specific IMF, metallicity, stellar age, and stellar population – binary or non-binary. Available absolute stellar metallicities in BPASS are 0.00001, 0.0001, 0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.01, 0.014, 0.02, 0.03, and 0.04, where the stellar ages range from 10⁶⁻¹¹ yr in increments of 0.1 dex. Since the Yggdrasil spectra assume a fixed metallicity, they are organized only by stellar age. All BPASS spectra are calculated between 10⁻³  and  10 |$\mu$|m, and all Yggdrasil spectra are calculated between 0.0091 and 9 |$\mu$|m. The stellar spectra for each burst of star formation are chosen by taking the age and metallicity of each SSP and selecting the nearest available neighbour in the database for the appropriate IMF and stellar population model. The total stellar SED is the sum of the SEDs from the individual SSPs that form every 5 Myr, starting when M_Halo = 10⁷ M_⊙ and ending at the observed redshift.

3.3 AGN spectrum

To include the effects of an AGN on our model galaxy, we implement a disk continuum from the SKIRTOR model (Stalevski et al. 2016) with a modification from Yang et al. (2022). To better match observations, Yang et al. (2022) introduced a free parameter, δ_AGN, in order to let the optical spectral slope deviate from the intrinsic shape in the piecewise power law:

$$\begin{eqnarray} \lambda L(\lambda) \propto \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\lambda ^{1.2} & 0.001 \le \lambda \le 0.01\ \ (\mu {\rm m}) \\ \lambda ^{0} & 0.01 \lt \lambda \le 0.1\ \ (\mu {\rm m}) \\ \lambda ^{-0.5 + \delta _{\rm AGN}} & 0.1 \lt \lambda \le 5\ \ (\mu {\rm m}) \\ \lambda ^{-3} & 5 \lt \lambda \le 1000\ \ (\mu {\rm m}) \\ \end{array}\right.. \end{eqnarray}$$

(3)

We take the AGN luminosity to be Eddington-limited L_AGN = 1.26 × 10³⁸ f_Edd f_BH (M_⋆/M_⊙) erg s⁻¹. The above continuum represents the face-on accretion disc luminosity. Dust obscuration and inclination angle are not included in this spectrum. We include δ_AGN as a free parameter in our study. We normalize the AGN SED with the bolometric luminosity L_AGN using the BH mass as a percentage of the stellar mass (f_BH ≡ M_BH/M_⋆) determined from the star formation history at the observed redshift, where we vary f_BH in our MCMC analysis. Because of the degeneracy of the Eddington and BH–stellar mass ratios, we do not vary f_Edd, and the variations in f_BH will account for any variations in the product of these two parameters.

3.4 Dust attenuation

Once we compute both the stellar SED and the AGN SED, we sum them to obtain the total SED of the galaxy. This spectrum is still the intrinsic, unattenuated spectrum in the rest-frame of the galaxy. In reality, as UV photons travel through the (dusty) interstellar medium, they are absorbed and re-emitted in various emission lines and in the infrared. Following the same method used in H22, we apply two models of dust attenuation; however, we do not consider nebular emission. The effects of this will be discussed in Section 5. First, absorption from the neutral IGM from Madau (1995) and second, extinction via the starburst reddening curve from Calzetti et al. (2000). To attenuate the spectrum, we compute an effective optical depth, τ_eff, and scale the intrinsic luminosity by |${\rm e}^{-\tau _\mathrm{eff}}$| to compute the attenuated luminosity. There are three contributions to the IGM absorption from Madau (1995):

• Metal-line blanketing:

  $$\begin{eqnarray} \tau _\mathrm{metal}(\lambda) = A_\mathrm{metal}\bigg (\frac{\lambda _\mathrm{obs}}{\lambda _\alpha } \bigg)^{1.68}\ \ (\lambda _\mathrm{obs} \lt \lambda _\alpha), \end{eqnarray}$$

  (4)

  where λ_obs = λ(1 + z) is the observer frame wavelength of the SED, A_metal = 0.0017, and λ_α = 1215.67 Å.

• Lyman series line blanketing:

  $$\begin{eqnarray} \tau _\mathrm{LB}(\lambda) = \sum _{j=2,i} A_j \bigg (\frac{\lambda _\mathrm{obs}}{\lambda _j} \bigg)^{3.46} \ \ (\lambda _{i+1} \lt \lambda _\mathrm{obs}/(1+z) \lt \lambda _i), \end{eqnarray}$$

  (5)

  where λ_j are the wavelengths of photons emitted by neutral hydrogen from an electron transition from excited state j to the first excited state. We compute 17 terms for τ_LB for j = 2 to j = 18, given the coefficients A_j¹ and 1/λ_j = R_H(1 − 1/j²) with |$R_\mathrm{H} \approx 1.0968 \times 10^{-3}\, \mathrm{nm}^{-1}$|⁠.

• Photoelectric absorption by intervening systems:

  $$\begin{eqnarray} \begin{split} \tau _\mathrm{PA}(\lambda) \simeq \ & 0.25 x_\mathrm{c}^3 (x_\mathrm{em}^{0.46}-x_\mathrm{c}^{0.46}) + 9.4x_\mathrm{c}^{1.5} (x_\mathrm{em}^{0.18}-x_c^{0.18}) \\ - &0.7x_c^3 (x_\mathrm{c}^{-1.32}-x_\mathrm{em}^{-1.32})-0.023 (x_\mathrm{em}^{1.68}-x_\mathrm{c}^{1.68}), \end{split} \end{eqnarray}$$

  (6)

  where x_c ≡ λ_obs/λ_L, λ_L = 911.75 Å, and x_em ≡ 1 + z. We use the approximation² of equation 16 in Madau (1995).

We also apply dust extinction to the spectrum using the starburst reddening curve from Calzetti et al. (2000), adding another contribution to the total optical depth

$$\begin{eqnarray} \tau _\mathrm{red}(\lambda) = \frac{k^{\prime }(\lambda)}{R_V^{\prime }}\, \tau _\mathrm{V} \end{eqnarray}$$

(7)

with

$$\begin{eqnarray} k^{\prime }(\lambda) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}2.659(-2.156 + 1.509/\lambda \\ \ \quad - 0.198/\lambda ^2 + 0.011/\lambda ^3) + R^{\prime }_\mathrm{V} & [0.12,\, 0.63)\, \mu \mathrm{m} \\ 2.659(-1.857 + 1.040/\lambda) + R^{\prime }_\mathrm{v} & [0.63,\, 2.20]\, \mu \mathrm{m} \end{array}\right. \end{eqnarray}$$

(8)

and |$R^{\prime }_\mathrm{V} = 4.05$|⁠. τ_V is the V-band optical depth and is a free parameter in our model. Finally, we calculate the dust attenuated luminosity from the intrinsic luminosity as L_obs(λ) = L(λ)exp [ − τ_eff(λ)] with τ_eff(λ) = τ_metal + τ_LB + τ_PA + τ_red.

3.5 Computing magnitudes from the SED

We convert the rest-frame attenuated luminosity density to a flux density given the redshift z and luminosity distance D_L to the galaxy. Using the JWST NIRCam filters, we convolve the throughput with the flux density and integrate the spectrum to obtain the total flux in each filter. All of our apparent magnitudes are computed in the AB system using

$$\begin{eqnarray} m_{\mathrm{AB}} = -2.5 \log _{10} \Bigg (\frac{\int f_\nu (\nu) T(\nu) \frac{d\nu }{\nu } }{f_{\mathrm{AB}} \int T(\nu) \frac{d\nu }{\nu }} \Bigg), \end{eqnarray}$$

(9)

where

$$\begin{eqnarray} f_\nu (\nu) = \frac{(L_\lambda ^\mathrm{stellar} + L_\lambda ^\mathrm{AGN})/\nu }{4\pi D_L^2 (1+z)}. \end{eqnarray}$$

(10)

Here, f_ν is the flux density per unit frequency; T is the throughput for a given NIRCam filter; f_AB = 3630.78 Jy is the zero-point for the AB magnitude system, and L_λ is the luminosity density (per unit wavelength). For the absolute UV magnitude, we follow Donnan et al. (2023) and convolve the flux density with a top-hat filter centred on 1500 Å with a width of 100 Å.

3.6 Validating the SED Model

Before running a full MCMC analysis, we validate our model spectra against the photometry of GL-z10 from Naidu et al. (2022, hereafter N22), using the same SED fitting parameters as them. We chose GL-z10 because this is the same object analyzed by H22, named GL-z9-1, which we analyze and discuss in this paper. The reason we used the SED fitting parameters from N22 rather than H22 is because N22 include the extinction A_V and the median stellar age, while H22 do not. The goal with this test is to confirm that our dust attenuation model is similar to the N22 model. We used star formation model B, the chab100 IMF, and we did not include an AGN or stellar binaries, i.e. we used a stellar-only SED. Rather than choosing a halo mass and computing the stellar mass from the Behroozi et al. (2013) SMHM relation, we used the stellar mass of GL-z10 directly from N22. This eliminates three of the five free parameters from our model, M_halo, f_BH, and δ_AGN. This leaves two free parameters: the exponential growth rate α and the V-band optical depth τ_V. We chose α = 0.65 such that the median stellar age is 165 Myr, similar to the value of |$163^{+20}_{-133}$| Myr for the same object analysed by N22. We used the dust attenuation A_5500Å; = 0.3 which corresponds to an optical depth τ_V ≈ 0.28, where τ_V = 0.4A_5500Å;/log₁₀(e). We found that when using the SED parameters from N22, our SED was brighter than the photometry by approximately 1 mag and the slope was slightly flatter. However, we were able to match the slope and normalization when we decreased the stellar mass by a factor of 0.65 and increased the dust attenuation by a factor of 1.5 (see Fig. 3).

Figure 3.

Model stellar spectra for GL-z10 from our initial test assuming no AGN, no stellar binaries, a chab100 IMF, and exponential halo growth with α = 0.65. The blue line shows the spectra using the stellar mass and dust attenuation from N22: log₁₀(M_⋆/M_⊙) = 9.6 and A_V = 0.3. The grey line shows the spectrum when decreasing the stellar mass by a factor of 0.65, while increasing the dust attenuation by a factor of 1.5. The red points with error bars show the magnitudes in NIRCam filters F150W, F200W, F277W, F356W, and F444W. The flux in each filter was taken from Naidu et al. (2022).

Ideally, using the fitted parameters from N22 with our stellar-only spectra would match the photometry perfectly, but our model produces a higher luminosity for a given stellar mass. However, the stellar mass and dust attenuation needed to match the photometry are within the uncertainties in stellar mass and dust attenuation from N22. For reference, these values are |$\log _{10}(\textrm {M}_\star /{ {\rm M_\odot }}) = 9.6^{+0.2}_{-0.4}$| and |$A_\mathrm{V} = 0.3^{+0.4}_{-0.2}$| and our values are log₁₀(M_⋆/M_⊙) = 9.41 and A_V = 0.45. From these results, it was unclear whether this discrepancy was due to a problem without our model or differences between the models. It is not a perfect ‘apples-to-apples’ comparison because we use a different star formation history model and different SED libraries.

To investigate this further, we compared our model spectra with fitted spectra of galaxies from the Hubble Ultra Deep Field (HUDF; McLure et al. 2011). The goal of the HUDF analysis was to compare our model to theirs in the most similar way possible. To do this, we turned off the AGN in our model and we chose the standard chab_100 IMF. The models we compared with also use a Chabrier (2003) IMF and stellar-only spectra without binaries. We only compared with HUDF galaxies that were modelled using a single burst of star formation. This eliminates the star-formation history from our model, i.e. we use neither SFH model A or B, we just use a spectrum from an SSP given the stellar mass, age, metallicity, and extinction in the V band from McLure et al. (2011). This means that the only differences between our model and theirs are the SED templates and the treatment of dust. Since we trust that BPASS provides high-quality SED templates, the primary difference is only our treatment of dust. If there was a systematic error in our dust attenuation code, we would expect this to be apparent because our SEDs of the HUDF galaxies would look much different than the SEDs from McLure et al. (2011).

Using the redshift, stellar mass, stellar age, metallicity, and V-band extinction from McLure et al. (2011) as inputs to our model, we found that the model spectra matched the HUDF galaxies reasonably well (see Fig. 4 for an example of one model), particularly at lower wavelengths. To assess the quality of our spectra, we computed the difference in absolute UV magnitude between our model and theirs, ΔM_UV. For example, for object HUDF_1173 with τ_V = 0, ΔM_UV = 0.05, and for object HUDF_1016 with τ_V = 1.01, ΔM_UV = 0.22. We chose to show HUDF_1016 in Fig. 4 because although our spectra do not fit the photometry as well as some of the other examples, HUDF_1016 has the largest optical depth and the primary reason for this test was to verify that our dust attenuation was working properly.

Figure 4.

Model stellar spectra for HUDF_1016 from our test assuming no AGN, no stellar binaries, a chab100 IMF, and a single burst of star formation. This is only one of the eight galaxies we compared with and it has neither the best nor the worst fit.

Fig. 4 shows that the spectra does not fit the 4.5 |$\mu{\rm m}$| flux for HUDF_1016. We tested seven other HUDF/ERS galaxies and all but one of the others fit the 4.5 |$\mu {\rm m}$| flux better than this one. Three of the objects fit the 4.5 |$\mu {\rm m}$| flux well. For these three objects, the difference between the observed IRAC2 (4.5 |$\mu {\rm m}$|⁠) magnitudes and the model magnitudes were 0.05, 0.26, and 0.29 mag. Overall, HUDF_1016, is neither the best nor the worst example from the objects we tested. Four of the of the seven other objects match the UV flux better than this one.

Since we are able to match the photometry of HUDF_1016 which has a large optical depth and the other objects, this confirms that the discrepancy between our model of GL-z10 and the model from N22 is not due to a problem with our dust attenuation. These results, combined with the fact that we can match the photometry of GL-z10 using N22 fitting parameters (within their uncertainty) validates our model. We now describe how we determine the posterior distributions of our input parameters for each object in our sample.

3.7 Exploring the parameter space

To find the parameters that produce a spectrum with the photometry that best matches the observed JWST values, we use the MCMC ensemble sampler emcee (Foreman-Mackey et al. 2013). The MCMC algorithm uses Bayesian inference to determine the probability distributions of the model parameters which also allows us to quantify the most likely values and their uncertainties. For each model permutation, we use emcee to vary four or five continuous parameters: the halo mass M_halo, the BH mass to stellar mass ratio f_BH, the AGN spectrum power law slope offset δ_AGN, the V-band optical depth τ_V, and for models using exponential halo growth, we also vary the growth rate α. emcee works by sending out ‘walkers’ to explore the parameter space. As the walkers move, we compute the photometric magnitudes associated with the walker’s position in the parameter space and compute the posterior probability of the parameters given the JWST observations. The posterior probability is given by

$$\begin{eqnarray} p(\Theta |O_{\mathrm{JWST}}) = \frac{p(O_{\mathrm{JWST}}|\Theta) p(\Theta)}{ p(O_{\mathrm{JWST}}) } , \end{eqnarray}$$

(11)

where Θ is the set of model parameters {M_halo, f_BH, δ_AGN, τ_V, α}; O_JWST is the set of JWST observed magnitudes in each filter; p(O_JWST|Θ) is the likelihood of a JWST observation given our parameters; p(Θ) are the parameter priors, and p(O_JWST) is the probability of the JWST observations. In practice, we do not need to know p(O_JWST) because we are only interested in maximizing the likelihood of the model parameters and p(O_JWST) is a constant that does not affect this process. We assume flat priors in the ranges for the parameters described in Table 3.

$$\begin{eqnarray} p(\Theta) = \left\lbrace \begin{array}{ll}1 \ , & \mathrm{if}\ \theta \in [\theta _{\mathrm{min}},\, \theta _{\mathrm{max}}]\ \ \forall \ \ \theta \in \Theta \\ 0 \ , & \mathrm{otherwise} \end{array} \right.\ \end{eqnarray}$$

(12)

We define the log of our likelihood function as −χ²/8 between the model magnitudes and the observed magnitudes because emcee seeks to maximize the log likelihood function. In our case, this means minimizing χ²,

$$\begin{eqnarray} \mathrm{ln}\ p(\Theta |O_{\mathrm{\it JWST}}) = -\frac{1}{2}\sum _{i} \biggl (\frac{m_i^{\mathrm{model}}-m_i^{\mathrm{obs}}}{\sigma _i^{\mathrm{obs,+}}-\sigma _i^{\mathrm{obs,-}}}\biggl)^2, \end{eqnarray}$$

(13)

where we sum over each filter i. Here, m_i are the magnitudes in each filter, and |$\sigma _i^+$| and |$\sigma _i^-$| are the upper and lower uncertainties in the observed magnitudes, respectively. Therefore, the reduced χ² is

$$\begin{eqnarray} \chi ^2 = \sum _{i} \biggl (\frac{m_i^{\mathrm{model}}-m_i^{\mathrm{obs}}}{\sigma _i^{\mathrm{obs}}}\biggl)^2 = -8\, \mathrm{ln}\ p(\Theta |O_{\mathrm{\it JWST}}), \end{eqnarray}$$

(14)

where |$\sigma _i = 0.5(\sigma _i^+ - \sigma _i^-)$|⁠.

We ran emcee for each of the seven model choices for each of the 12 galaxies in our sample for a total of 84 different runs. For each run we used 288 chains of walkers with a maximum of 20 000 steps, though not all models made it to 20 000 steps in the allocated computing time. Each model was run in parallel on a cluster using 144 cores each for a maximum of 8 h. To verify that the Markov chains had converged to a steady state distribution, we checked that the chains were at least 50 times longer than the integrated autocorrelation time, |$\mathrm{\hat{\tau }}$|⁠, and that the autocorrelation time had not changed by more than 1 per cent between steps at the end of the run. The autocorrelation time is the number of steps it takes for the walkers to settle into the true distribution of the parameter space. Eight out of the 84 models did not meet these standards. These were the Yggdrasil models for object 10234, and six other models for object 6115. However, the best-fitting models had chains longer than 50|$\mathrm{\hat{\tau }}$|⁠, and all but three of the best models had chains longer than 100|$\mathrm{\hat{\tau }}$|⁠. The best model for each object – except 6115 and GL-z9-2 – also met the condition of |$\mathrm{\hat{\tau }}$| not changing by more than 1 per cent by the time the runs finished. For the two objects where this condition was not met, the fractional change in |$\mathrm{\hat{\tau }}$| was only marginally larger than 1 per cent. For details on how the autocorrelation time is calculated and its relevance to convergence, see Foreman-Mackey et al. (2013).

At each step in the chain, there is a certain chance that a walker will move to a new position in the parameter space. In general, a walker has a lower – although non-zero – chance to move to a new position if the posterior probability p(Θ|O_JWST) at the new position is very low. When the move is rejected, the walker will remain in place and the next step in the chain will be a duplicate of the previous step. This often occurs when a walker reaches an edge of the parameter space and the next step would take it outside the range of the priors. For this reason, we remove the duplicate samples and only compute statistics for the distributions using the unique samples. We also discard the initial 3|$\mathrm{\hat{\tau }}$| steps as ‘burn-in’ to ensure we only work with parameters that are sampled from the true distribution. We calculate the median and the uncertainty in the distributions of each parameter from the 16th, 50th, and 84th percentiles of the distributions.

4 RESULTS

For most galaxies, we find that our analysis pipeline determines a stellar mass this is consistent with A22/H22 within their uncertainties. For half of the objects, the best-fitting parameters are higher than their best-fitting values but still within their uncertainties. Fig. 5 shows the median stellar mass from each model for each galaxy. The black lines indicate the stellar mass from A22/H22, except for galaxy GL-z9-2 which comes from Castellano et al. (2022a). For most galaxies, the BPASS models lie above the A22/H22 stellar mass, whereas the Yggdrasil models lie well below the black line and tend to have much larger error bars. While we set out to see if we could determine stellar masses that are lower than what the literature predicts by including other luminous sources, we found the opposite. We discuss the exact reasons for this overestimate in Section 5.

Figure 5.

The emcee-predicted median stellar mass for each galaxy according to our seven fixed models is shown in with lettered markers. The error bars represent the 16th and 84th percentiles of the emcee distribution after discarding any data prior to three times the integrated autocorrelation time, which indicates when the walkers begin converging. The black line is the stellar mass determined by A22 and H22, and the black shaded region shows their error if provided. Note that the stellar mass and error shown in black for galaxy GL-z9-2 is predicted by Castellano et al. (2022a), not A22 or H22. It is evident that a significant number of our predictions overestimate the stellar mass as compared to those predicted by H22 and A22. The AGN-dominated Yggdrasil spectra, which are used in models F and G, recover significantly lower stellar masses as expected when there is contribution from a central BH.

A useful property about this approach is that we can estimate a typical star formation history of the galaxies, beginning at its onset, with any scatter being captured by the MCMC. The full star formation history for each galaxy and their best-fitting model can be seen in Figs 6–8. The best-fitting spectra for each galaxy is shown in Figs 9–11. We present the best fit and median parameters for the best-fitting model for all galaxies in Tables 4 and 5, respectively.

Figure 6.

The halo growth and star formation histories for each galaxy from A22 given the best-fitting parameters. Upper left: Halo mass history. Upper right: Stellar mass history. The scatter points indicate the predicted stellar mass from JWST. Lower left: SFR history averaged over the last 50 Myr. Lower right: SFE history.

Figure 7.

The halo growth and star formation histories for each galaxy from H22 given the best-fitting parameters. Upper left: Halo mass history. Upper right: Stellar mass history. The scatter points indicate the predicted stellar mass from JWST. Lower left: SFR history averaged over the last 50 Myr. Lower right: SFE history.

Figure 8.

The stellar population from A22 (top panels) and H22 (bottom panels), given the best-fitting parameters. Left-hand panels: The cumulative star formation history as a function of stellar age. Right-hand panels: The normalized distribution of stellar metallicities.

Figure 9.

Fitted SED given the best-fitting models for galaxies (a) 6115, (b) 1514, (c) 1696, and (d) 2462. Blue triangles indicate the measured photometry from our fitted spectra and the pink circles are the observed photometry values and errors. Annotated are the measured redshifts and stellar masses from A22, as well as our stellar mass given the best-fitting halo mass and our χ² value.

Figure 10.

Fitted SED given the best-fitting models for galaxies (a) 2779, (b) 6878, (c) 10234, and (d) GL-z9-1. Blue triangles indicate the measured photometry from our fitted spectra and the pink circles are the observed photometry values and errors. Annotated are the measured redshifts and stellar masses from A22 and H22, as well as our stellar mass given the best-fitting halo mass and our χ² value.

Figure 11.

Fitted SED given the best-fitting models for galaxies (a) GL-z9-2, (b) GL-z12-1, (c) CR2-z12-1, and (d) CR2-z16-1. Blue triangles indicate the measured photometry from our fitted spectra, and the pink circles are the observed photometry values and errors. Annotated are the measured redshifts and stellar masses from H22, as well as our stellar mass given the best-fitting halo mass and our χ² value.

Table 4.

Best-fitting model parameters.

Galaxy	Model	log10(Mhalo/M⊙)	log10(M⋆/M⊙)	log10(fBH)	δAGN	τV	α	χ2
10234	B_chab100_bin	10.09	8.44	–4.00	–0.94	0.00	0.80	78.75
1514	A_chab100_sin	11.10	9.58	–2.65	–0.97	1.07	–	6.35
1696	A_chab100_sin	10.68	9.16	–3.99	–0.99	0.24	–	5.35
2462	A_chab100_sin	11.47	9.95	–3.99	–0.98	0.76	–	9.55
2779	B_chab100_bin	10.37	8.67	–3.95	–0.80	0.19	0.40	7.04
6115	A_yggdrasil_sin	9.70	7.53	–1.04	–0.84	0.00	–	2.28
6878	A_chab100_sin	10.71	9.18	–2.79	–0.88	0.41	–	0.66
CR2-z12 − 1	B_chab100_bin	10.18	8.65	–3.82	0.73	0.02	0.41	6.07
CR2-z16 − 1	A_chab100_sin	11.04	9.52	–3.53	–0.77	0.36	–	4.12
GL-z12 − 1	B_yggdrasil_sin	10.48	8.96	–3.96	–0.93	0.00	0.44	4.05
GL-z9 − 1	A_chab100_bin	10.77	9.25	–2.71	–0.09	0.36	–	0.43
GL-z9 − 2	B_chab100_bin	10.34	8.72	–3.91	–0.81	0.01	0.42	4.47

Table 5.

Median of the parameter distributions and their uncertainties for the best-fitting models.

Galaxy	Model	log10(Mhalo/M⊙)	log10(M⋆/M⊙)	log10(fBH)	δAGN	τV	α
10234	B_chab100_bin	10.11|$^{+0.02}_{-0.02}$|	8.50|$^{+0.04}_{-0.05}$|	−3.65|$^{+0.22}_{-0.32}$|	−0.44|$^{+0.38}_{-0.66}$|	0.02|$^{+0.01}_{-0.02}$|	0.73|$^{+0.07}_{-0.05}$|
1514	A_chab100_sin	11.08|$^{+0.10}_{-0.09}$|	9.56|$^{+0.10}_{-0.09}$|	−3.28|$^{+0.42}_{-0.44}$|	−0.10|$^{+0.55}_{-0.63}$|	0.99|$^{+0.08}_{-0.09}$|	–
1696	A_chab100_sin	10.67|$^{+0.04}_{-0.04}$|	9.15|$^{+0.04}_{-0.04}$|	−3.56|$^{+0.27}_{-0.34}$|	−0.15|$^{+0.55}_{-0.69}$|	0.25|$^{+0.03}_{-0.04}$|	–
2462	A_chab100_sin	11.44|$^{+0.07}_{-0.07}$|	9.92|$^{+0.07}_{-0.07}$|	−3.59|$^{+0.26}_{-0.34}$|	−0.09|$^{+0.59}_{-0.69}$|	0.74|$^{+0.07}_{-0.07}$|	–
2779	B_chab100_bin	10.35|$^{+0.04}_{-0.04}$|	8.63|$^{+0.09}_{-0.07}$|	−3.43|$^{+0.34}_{-0.40}$|	−0.00|$^{+0.60}_{-0.62}$|	0.29|$^{+0.09}_{-0.09}$|	0.55|$^{+0.10}_{-0.13}$|
6115	A_yggdrasil_sin	10.12|$^{+0.05}_{-0.03}$|	8.39|$^{+0.10}_{-0.07}$|	−3.19|$^{+0.48}_{-0.81}$|	−0.30|$^{+0.53}_{-0.73}$|	0.07|$^{+0.04}_{-0.06}$|	–
6878	A_chab100_sin	10.70|$^{+0.05}_{-0.05}$|	9.17|$^{+0.05}_{-0.05}$|	−3.33|$^{+0.38}_{-0.40}$|	−0.10|$^{+0.56}_{-0.65}$|	0.37|$^{+0.04}_{-0.05}$|	–
CR2-z12 − 1	B_chab100_bin	10.21|$^{+0.10}_{-0.12}$|	8.69|$^{+0.14}_{-0.12}$|	−3.36|$^{+0.39}_{-0.50}$|	−0.16|$^{+0.53}_{-0.66}$|	0.21|$^{+0.11}_{-0.12}$|	0.62|$^{+0.12}_{-0.11}$|
CR2-z16 − 1	A_chab100_sin	10.92|$^{+0.23}_{-0.20}$|	9.40|$^{+0.23}_{-0.20}$|	−2.99|$^{+0.60}_{-0.72}$|	−0.12|$^{+0.56}_{-0.66}$|	0.29|$^{+0.16}_{-0.16}$|	–
GL-z12 − 1	B_yggdrasil_sin	10.44|$^{+0.10}_{-0.11}$|	8.91|$^{+0.10}_{-0.11}$|	−3.33|$^{+0.40}_{-0.51}$|	−0.22|$^{+0.50}_{-0.66}$|	0.09|$^{+0.06}_{-0.08}$|	0.60|$^{+0.12}_{-0.12}$|
GL-z9 − 1	A_chab100_bin	10.84|$^{+0.25}_{-0.17}$|	9.31|$^{+0.25}_{-0.17}$|	−2.99|$^{+0.59}_{-0.67}$|	−0.10|$^{+0.55}_{-0.63}$|	0.39|$^{+0.15}_{-0.14}$|	–
GL-z9 − 2	B_chab100_bin	10.27|$^{+0.07}_{-0.07}$|	8.57|$^{+0.14}_{-0.15}$|	−3.30|$^{+0.42}_{-0.52}$|	−0.17|$^{+0.53}_{-0.65}$|	0.18|$^{+0.10}_{-0.14}$|	0.64|$^{+0.11}_{-0.09}$|

Note. The upper and lower uncertainties are computed from the 84th and 15th percentiles of the distributions, respectively.

All galaxies follow a somewhat similar halo growth history, although for all haloes, the limit on the SFE can be seen limiting the stellar mass when there is an abrupt knee in the SFR, stellar mass, and SFE plots. The most massive progenitor of the A22 and H22 galaxies reaches the atomic cooling limit, where galaxy formation commences, at redshifts z = 20−30 when the Universe’s age is less than 200 Myr. Pop III star formation in the progenitors may occur as early as z ∼ 40 when their halo masses reach ∼10⁶ M_⊙ as these rare objects undergo an early and rapid assembly. The galaxies from H22 generally exist at higher redshifts than the A22 galaxies, resulting in shorter star formation sequences with the oldest SSPs having ages of 100–200 Myr and a lower, less diverse metallicity distribution as compared to the lower redshift galaxies from A22. The oldest stars in the A22 galaxies show an age range of 150–300 Myr, shown in Fig. 8. In general, all galaxies are metal poor, as expected given their stellar masses and high redshifts (Torrey et al. 2019).

4.1 Galaxy groupings

There are some important comparisons to make across all the galaxies we modelled. In this section, we group the galaxies based on different aspects of the fitting results.

Quality of the spectral fitting: The best measure of how well our models fit the photometry is χ². Objects GL-z9-1 and 6878 have the lowest χ² values, with their best-fitting models having χ² < 0.7. For these objects, the model magnitudes in every filter all lie within the error bars of the photometric magnitudes from A22 and H22 (see Table A1 in the Appendix). For the other objects, the model magnitudes in at least one filter lie outside the photometric standard error from A22 and H22 and their best-fitting χ² ranges from 2.28 to 9.55 except for 10234. These objects fit the photometry well in nearly every filter. 10 234 has the worst fit by far with χ² > 78 for all models. Our pipeline was not able to reproduce the photometry of this galaxy. This is most likely because F200W is quite bright and F350W is quite dim. This means the spectrum must have a steep negative slope in the 2–3.5 |$\mu$|m range. Our models were not able to find a spectrum with a steep enough slope, while also matching filters F150W and F444W.

Table 6.

Best parameters, median of the parameter distributions and their uncertainties for GL-z9 − 1. Similar tables for all objects are given in the online supplementary material.

GL-z9 − 1
Best parameters
Model	log10(Mhalo/M⊙)	log10(M⋆/M⊙)	log10(fBH)	δAGN	τV	α	χ2
A_chab100_bin	10.77	9.25	–2.71	–0.09	0.36	–	0.43
A_chab100_sin	10.60	9.08	–2.20	–0.55	0.30	–	0.45
A_135_100_bin	10.82	9.30	–2.74	–0.14	0.35	–	0.44
B_chab100_bin	10.80	9.28	–2.77	–0.16	0.39	0.74	0.43
B_135_100_bin	10.82	9.30	–2.83	0.01	0.36	0.80	0.44
B_yggdrasil_sin	9.48	6.97	–0.02	–0.27	0.02	0.44	3.61
A_yggdrasil_sin	9.53	7.06	–0.14	–0.25	0.01	–	3.62
Median parameters
Model	log10(Mhalo/M⊙)	log10(M⋆/M⊙)	log10(fBH)	δAGN	τV	α	χ2	χ2(θmed)
A_chab100_bin	10.84|$^{+0.25}_{-0.17}$|	9.31|$^{+0.25}_{-0.17}$|	−2.99|$^{+0.59}_{-0.67}$|	−0.10|$^{+0.55}_{-0.63}$|	0.39|$^{+0.15}_{-0.14}$|	–	5.14|$^{+3.38}_{-7.68}$|	1.58
A_chab100_sin	10.71|$^{+0.21}_{-0.16}$|	9.19|$^{+0.21}_{-0.16}$|	−2.88|$^{+0.64}_{-0.68}$|	−0.14|$^{+0.53}_{-0.64}$|	0.27|$^{+0.12}_{-0.13}$|	–	6.08|$^{+3.30}_{-7.53}$|	3.82
A_135_100_bin	10.88|$^{+0.27}_{-0.18}$|	9.36|$^{+0.27}_{-0.18}$|	−2.99|$^{+0.60}_{-0.67}$|	−0.11|$^{+0.54}_{-0.63}$|	0.38|$^{+0.15}_{-0.14}$|	–	5.19|$^{+3.40}_{-7.68}$|	1.70
B_chab100_bin	10.86|$^{+0.24}_{-0.17}$|	9.34|$^{+0.24}_{-0.17}$|	−3.01|$^{+0.57}_{-0.66}$|	−0.12|$^{+0.53}_{-0.62}$|	0.37|$^{+0.15}_{-0.15}$|	0.60|$^{+0.12}_{-0.12}$|	5.04|$^{+3.28}_{-7.46}$|	1.57
B_135_100_bin	10.91|$^{+0.26}_{-0.18}$|	9.38|$^{+0.26}_{-0.18}$|	−3.03|$^{+0.56}_{-0.66}$|	−0.12|$^{+0.53}_{-0.62}$|	0.36|$^{+0.15}_{-0.15}$|	0.60|$^{+0.12}_{-0.11}$|	5.10|$^{+3.31}_{-7.49}$|	1.55
B_yggdrasil_sin	10.59|$^{+0.51}_{-0.38}$|	9.07|$^{+0.86}_{-0.38}$|	−2.21|$^{+0.69}_{-1.22}$|	−0.28|$^{+0.43}_{-0.67}$|	0.35|$^{+0.19}_{-0.20}$|	0.59|$^{+0.11}_{-0.12}$|	12.69|$^{+5.24}_{-8.16}$|	21.28
A_yggdrasil_sin	10.56|$^{+0.49}_{-0.38}$|	9.04|$^{+0.86}_{-0.38}$|	−2.21|$^{+0.72}_{-1.23}$|	−0.25|$^{+0.45}_{-0.68}$|	0.37|$^{+0.19}_{-0.20}$|	–	13.41|$^{+5.47}_{-8.24}$|	19.80

Note. The upper and lower uncertainties are computed from the 84th and 16th percentiles of the distributions, respectively. χ²(θ_med) is χ² for the spectra computed from the median parameters in columns 3–8 rather than the median of the χ² distribution.

We find that for seven of the twelve objects (10234, 1514, 1696, 2462, 2779, 6878, and GL-z9-2) models A_yggdrasil_sin and B_yggdrasil_sin fit significantly worse than the models which use BPASS. This is typically because the Yggdrasil spectra do not have a pronounced Balmer break like the BPASS spectra do. This makes it difficult for the Yggdrasil models to fit the photometry of F444W at z ≤ 11. For the higher redshift galaxies, CR2-z12-1 and CR2-z16-1, the BPASS models are still best but χ² for the Yggdrasil models is not significantly larger. While the Yggdrasil models typically do not perform as well as the BPASS models, there are two objects (6115 and GL-z12-1) for which the Yggdrasil spectra is the best model; however, in 6115 the AGN component dominates over the stellar one so we cannot conclude that this suggests a top-heavy IMF in this object. Conversely, the stellar component dominates the spectrum of GL-z12-1 so this object may have a top heavy IMF. We will further discuss these two objects and the IMF variations in Sections 4.2 and 4.3.

Stellar mass: The best-fitting parameters from our MCMC analysis show that six of the objects have higher best-fitting stellar masses than predicted by A22 and H22. These are objects 10234, 1696, 2462, CR2-z12-1, CR2-z16-1, and GL-z12-1. Of these, CR2-z12-1 shows the largest difference in stellar mass than the value from the literature. The best-fitting stellar mass of CR2-z12-1 is log₁₀(M_⋆/M_⊙) = 8.65 compared to |$\textrm {log}_{10} (\mathrm{M}_{\star }/{ {\rm M_\odot }}) = 7.88_{-0.27}^{0.95}$| from H22. Although this difference is large, our value is still within the H22 error bars. This is true for all H22 objects. Even when our stellar masses differ greatly from the H22 values, the best-fitting stellar masses are all within the error bars from H22.

Four of the objects have comparable stellar masses to the values from A22 and H22. These are 1514, 2779, 6878, and GL-z9-2. All of these objects have |log₁₀(M_⋆,model/M_⋆,JWST)| < 0.22. Of these, 2779’s stellar mass is closest to the value from A22 with log₁₀(M_⋆,model/M_⋆,JWST) ≈ −0.03. For objects 1514, 2779, and 6878, their stellar masses, while comparable, are still lower than the values from A22 and H22. Objects 1514 and 6878 may have lower stellar masses because the spectrum is not completely dominated by the stellar component due to the presence of an AGN.

Object 6115 shows the largest difference in stellar mass. The stellar mass of 6115 is log₁₀(M_⋆/M_⊙) = 7.53 compared to log₁₀(M_⋆/M_⊙) = 8.40 from A22. This enormous difference is because the spectrum of 6115 is completely dominated by the AGN component. With a large f_BH, the AGN spectrum fits the photometry of 6115 well with χ² = 2.28. This object will be discussed in more detail in Section 4.2.

Component Spectra and Black Hole Fractions: We find that the spectra are dominated by the stellar component for most of the objects in our sample. This is the case for 10234, 1696, 2462, 2779, CR2-z12-1, CR2-z16-1, GL-z12-1, and GL-z9-2. For these objects, log₁₀(f_BH) is near the minimum value of –4, indicating they do not have a bright AGN, or if there is an AGN, it is obscured. There are three objects, 1514, 6878, and GL-z9-1, which have a noticeable contribution from an AGN. These models have larger log₁₀(f_BH) values ranging from –2.79 to –2.65. Object 6115 is much different from all other objects. It is the only one whose spectrum is completely dominated by a bright AGN with log₁₀(f_BH) = −1.04 because both BPASS and Yggdrasil stellar models could not provide a good fit to the photometry.

4.2 Individual objects

While many of the objects can be grouped by their similarities as shown previously, there are a few objects that stand out. Here we describe the most interesting objects in more detail.

6115: This object has a photometric redshift of |$z = 10.94_{-0.15}^{+0.12}$| and a stellar mass of log₁₀(M_⋆/M_⊙) = 8.4 (A22). The best-fits from all models predict lower stellar masses than A22, with Yggdrasil models predicting larger best-fitting stellar masses than the BPASS models by 1 dex, and smaller BH masses also by 1 dex. All χ² values are below 2.72 for this galaxy, indicating a good reproduction of photometry. Independent of the model choice in the best-fitting model, the photometry is reproduced with a massive BH with log₁₀(f_BH) ≳ −1; however, the median parameters suggest a more typical central BH with f_BH ≃ 10⁻³.

Best-fitting parameters: We find the spectra that best fits this object comes from model A_yggdrasil_sin, with a best-fitting stellar mass of log₁₀(M_⋆/M_⊙) = 7.53 (see Fig. 9a). This is the lowest stellar mass of all the objects in our sample. The best-fitting parameters for all models are in the supplementary data and the best-fitting spectra and photometry can be seen plotted in Fig. 9(a). Our parameters match the photometry in almost every filter except F277W, which is outside the error bounds of A22 by a magnitude of 0.05. The predicted optical depth in every best-fitting model was between 0 and 0.02, with four out of the seven models predicting an optical depth of 0.

Median parameters: We find a median stellar mass of |$\textrm {log}_{10}(\mathrm{{\it M}}_{\star }/{ {\rm M_\odot }}) = 8.50_{-0.16}^{+0.19}$| from the best model. This larger stellar mass compared to the best-fitting value is due to the lower median |$\mathrm{log}_{10}(f_\mathrm{BH}) = -3.19_{-1.03}^{+0.49}$|⁠. This shows the trade-off between a luminous BH and a bright stellar component. The reason for the difference in stellar mass and BH mass is the f_BH distribution is peaked at around log₁₀(f_BH) ≈ −3 but has a flat tail that extends to larger values where the best-fitting value is found. The median parameters can be seen in the supplementary data. All models have similar median fits and no significant trend is found. To see a posterior distribution of spectra for the galaxy along with the χ² distribution, refer to Fig. 12.

Figure 12.

Model spectra coloured by χ² using the best model from four galaxies: GL-z9-1 with model B_chab100_bin (top left), 6115 model A_yggdrasil_sin (top right), 6878 model A_chab100_sin (bottom left), 1514 model A_chab100_sin (bottom right). For each galaxy, the SEDs were generated by randomly sampling 400 sets of model parameters θ from within the 16th–84th percentile range of each posterior parameter distribution, i.e. within the median +/− standard errors shown in Table 5. The pink line shows the SED computed from the median parameters of the best model. The red markers with error bars show the JWST magnitudes and standard error in filters F150W, F200W, F277W, F365W, and F444W.

Unlike the other objects, the best-fitting corner plot³ for object 6115 shows a noticeable peak in δ_AGN while still exploring the full parameter space at the wings of the distribution. This is because when f_BH is large, the AGN spectrum becomes important and δ_AGN can be constrained to fit the photometry. Whereas for other objects with a small f_BH, varying δ_AGN has little to no effect on the spectrum.

GL-z9-1: This galaxy stands out from the rest in our sample because it is the only one that has comparable stellar and AGN luminosities for the non-Yggdrasil spectra. While the stellar component is still brighter, the AGN spectrum makes a significant contribution. This is shown in the rightmost panel of Fig. 13. It is also the object with the best fit to the photometry with χ² = 0.43. This object has a photometric redshift of |$z = 10.49_{-0.72}^{+0.53}$| and a stellar mass of |$\textrm {log}_{10} (\mathrm{\it M}_{\star }/{ {\rm M_\odot }}) = 9.04_{-0.74}^{+0.61}$| (H22). For all BPASS models, the best-fitting parameters and the median parameters give a slightly larger stellar mass than the value from H22. The best-fitting spectra match the JWST photometry in every filter. Our magnitudes are all within the error of the photometric magnitudes from H22. This is shown in Fig. 10(d) and Table A1.

Figure 13.

A comparison of SED models showing a stellar dominated SED, an AGN dominated SED, and an SED where the stellar and AGN components are both contributing to the total. The grey lines indicate the total SED, the blue lines indicate the stellar spectrum, and the pink lines indicate the AGN spectrum. Left-hand panel: SED generated for galaxy 1696 from model A_chab100_sin. This model shows a completely stellar dominated SED. Middle panel: SED generated for galaxy 1696 for model B_yggdrasil_sin. This model shows a completely AGN-dominated SED. Right-hand panel: SED generated for galaxy GL-z9-1 for the best-fitting model A_chab100_bin. This model shows comparable contributions from the stellar and AGN components to the total SED.

Best-fitting parameters: We find the spectra that best fits this object comes from model A_chab100_bin, with a best-fitting stellar mass of log₁₀(M_⋆/M_⊙) = 9.25 (see Fig. 10d). The best-fitting parameters for all models can be seen in the Table 6. All BPASS models fit this galaxy extremely well with χ² ranging from 0.43 to 0.45 so there is nothing particularly special about model A_chab100_bin. Models using the Yggdrasil spectra for a top heavy IMF have the largest χ² ≈ 3.6. Unlike the best models for most other galaxies, GL-z9-1 consistently has a relatively large values of log₁₀(f_BH) ranging from –2.83 to –2.20.

Median parameters: We find a median stellar mass of |$\textrm {log}_{10} (\mathrm{M}_{\star }/{ {\rm M_\odot }}) = 9.31_{-0.17}^{+0.25}$| from the best model. The median parameters can be seen in the Table 6. Compared to the best-fitting parameters, the median parameters give larger stellar masses and smaller BH fractions. As with the best-fitting parameters, the Yggdrasil models have the largest BHs. As expected, the spectra from the median parameters have larger χ² than the best-fitting spectra.

The corner plot for GL-z9-1 has a couple of interesting features. The M_halo versus f_BH plot shows a negative slope with a slight tail that extends to larger f_BH. This shows the trade-off between BH mass and halo mass. To match the photometry, if the AGN is brighter, the stellar component must be dimmer. This results in wider distributions for both f_BH and M_halo than when only one component dominates. The δ_AGN distribution is skewed toward negative values with a peak at –0.5. Since f_BH is larger, the AGN component of the spectra is non-negligible, and therefore δ_AGN becomes an important parameter. The peak near –0.5 is important because this is consistent with measurements of δ_AGN from observations (Yang et al. 2022).

The top left-hand panel of Fig. 12 shows the distribution of spectra within the range of the median parameters and their standard errors, coloured by χ², for the best model of GL-z9-1. Compared to the other panels, GL-z9-1’s spectra have a wider distribution in magnitudes at low χ². This is mainly because the observed NIRCam/F277W magnitude is relatively bright causing bump in the photometry, whereas our spectra are more monotonic in this bandwidth.

4.3 Variation on models

The corner plots for all models reveal how emcee constrains the parameters that best match the data. From this, we can see what parameters are most and least important to the model.

Host halo mass: By far and unsurprisingly, the most important parameter that controls the goodness of the fit to the JWST photometry is the halo mass M_halo. As will be discussed shortly, the AGN model cannot match the photometry very well on its own while the BPASS stellar spectrum can. Since the halo mass directly controls the stellar mass and the star formation history, M_halo is the parameter that controls the results the most. This can be seen in the corner plots where M_halo, and therefore M_⋆, typically have very peaked distributions. When the Yggdrasil spectrum is used instead of BPASS, emcee has a harder time fitting these spectra to the photometry, and the M_halo distribution can sometimes be doubly peaked. In these cases, the AGN spectra will play a more dominant role in fitting the photometry since the Yggdrasil spectra cannot fit the data very well, which we will discuss in Section 5.

Halo growth models (A or B): The exact history of the haloes does not have a significant impact. For most of the galaxies in our sample, χ² does not differ greatly between models using growth model A or B, indicating that there is not a strong preference for one halo growth model over the other. This can be seen in the distributions of α in our corner plots, where the α parameter tends to have a fairly broad distribution. The value of α that emcee lands on is typically around 0.6, which is the median of our range of α. This implies that the most recent star formation is what is most important, as the history of the halo does not play a major role.

IMF: Both Yggdrasil models tend to predict higher BH masses, controlled by the f_BH parameter. This results in smaller stellar masses since the AGN spectra tends to dominate the fitting, but the Yggdrasil models also tend to have higher χ² values than the BPASS models, indicating a systematic issue between matching Yggdrasil spectra to the observed photometry. Some potential causes of this systematic error are discussed in Section 5. There is no significant trend between BPASS models; both the Chabrier and Standard IMF models produce similar results regardless of which one is chosen as the best-fit model. Their χ² values are all consistent with each other and most of the BPASS models are able to fit the spectra with low χ², showing good fits to photometry.

Single versus binary stellar population: Only BPASS models both single and binary stellar populations, while Yggdrasil only models SSPs. Of the 10 galaxies that had best-fit parameters predicted by a BPASS IMF model, six of them have binary stellar population spectra while four have SSP spectra. For each galaxy, BPASS IMFs produce similar parameter values regardless of stellar population.

δ_AGN : The AGN spectrum does not fit the photometry very well by itself, and therefore, it is not typically used by emcee as the dominant fitting spectra. Because of this, the slope of the AGN spectrum, controlled by our parameter δ_AGN does not vary that much, and when f_BH is small, δ_AGN does not affect the SED at all so the best-fit value can vary greatly. This will be discussed further in Section 5.

5 DISCUSSION

There are a number of reasons why we may be obtaining higher best-fit stellar masses than those predicted by A22 and H22, given our theoretical model. These may include an insensitivity to the halo growth model, a dominant stellar spectra component, and ill-fitting AGN and top-heavy IMF stellar spectra, and the absence of nebular emission in our spectra. We will explore each of these in this section, and focus individually on the AGN spectra and the IMFs in Sections 5.1 and 5.2, respectively.

A common theme across many of our models and most of the objects analyzed is that the distributions of the exponential halo growth rate α spans the full range of the priors from 0.4 to 0.8. The median values of α are typically symmetric about 0.6 (see Tables in supplementary data). This indicates that varying α does not strongly affect the probability of our spectra matching the JWST photometry for a given object. Recall that our haloes grow exponentially starting from the redshift when they had a mass of 10⁷M_⊙ until the observed redshift, sampled every 5 Myr. The growth rate does not have a significant effect on the spectra because the most recent bursts of star formation dominate the stellar spectrum. For each object, the final burst alone makes up about 35–55 per cent of the total flux of the galaxy in the NIRCam bandwidth. Although this approach does not capture the individual galaxy progenitors, as in a method based on merger trees, the numerous mergers of low-mass galaxies, which have a large scatter in their star formation histories (e.g. Xu et al. 2016; Gutcke et al. 2022), average out, producing a more massive galaxy similar to the median values found in this work.

As we showed in Section 4, eight out of the 12 objects had log₁₀(f_BH) near the minimum value of –4. When the BH is this small, the stellar flux dominates the continuum, and therefore the offset in the AGN slope is not an important parameter. This is why we see wide distributions in f_BH with large uncertainty. We stress that in these cases, f_BH values have no significant meaning, regardless of their value. This is not the case for the four other objects where f_BH is larger, and the AGN contributes a noticeable amount to the SED, especially for 6115 where the AGN component is dominant.

5.1 AGN Spectra

In order to more accurately model the AGN, we attempted to include realistic AGN spectra from the SKIRTOR data base of a dusty torus at different inclination angles. Following the parameters given in fig. 4 in Stalevski et al. (2012), we ran emcee with the inclination angle as a free parameter. Since the dust-absorbed BPASS spectra match the shape of the photometry well when scaled by stellar mass, the best-fitting parameters still trended towards stellar-dominated spectra by predicting best-fitting inclinations greater than 40 deg, at which point a majority of the UV flux is attenuated into the IR, making AGN contributions negligible within our wavelength range of interest. We therefore decided to use the point-source emission model of AGNs and determine the AGN spectral slope in the UV with the δ_AGN parameter, assuming maximum possible contribution from the AGN.

As mentioned earlier to be consistent with previous studies on these objects, we did not include nebular emission lines from star-forming regions and AGN. Photometric detections of z ≥ 6 galaxies with Hubble Space Telescope (HST) showed that emission lines redwards of the Balmer break can contribute a significant amount of flux to broad-band filters (e.g. Smit et al. 2014). Emission-line spectra of H ii regions can indicate the likelihood of an AGN based on [O iii]/H β and [N ii]/H α ratios as shown in BPT diagrams (Baldwin, Phillips & Terlevich 1981) and will provide stronger evidence than colour–colour diagnostics of the existence of an AGN. Katz et al. (2022b) investigated the emission lines of galaxies identified in the JWST Early Release Observations and have used these ratios to determine whether these galaxies were likely to host AGNs. Given a theoretical model that is able to incorporate these elements into the spectra, we may be able to more easily fit a realistic AGN SED along with our realistic stellar SED and determine whether our high-redshift candidates host AGNs. However, this is highly dependent on the redshift of these galaxies. In practice, if the redshifts from A22 and H22 used in this work are correct, i.e. these galaxies are actually all at z > 9, the emission lines mentioned above would be redshifted to wavelengths greater than what is detectable by NIRCam. Therefore, these emission lines would have no impact on the photometry and including them in our spectra would not affect our results.

Given that we only have a few best-fitting models where an AGN possibly contributes to the continuum, we have too small of a sample to make any general predictions about the nature of the BHs and their seeding (Volonteri et al. 2022). If the AGN contribution is subdominant to the stellar component, colour–colour diagnostics alone are not sufficient to ascertain the existence of a central BH, and nebular emission lines detected with NIRSpec (Curtis-Lake et al. 2022) will be needed to determine their mass and accretion rate estimates and compare them with their host galaxy properties (e.g. Volonteri et al. 2022). With this additional information, constraints can be placed on the galaxy–BH relation during the Epoch of Reionization and possibly BH seeding mechanisms.

5.2 Stellar spectra and IMFs

We explored the possibility of an overabundance of massive stars, i.e. a top-heavy primordial IMF, by using Yggdrasil (Zackrisson et al. 2011). As discussed in Section 3.2, metal-free stars are thought to be generally massive because of the weak cooling rates associated with molecular hydrogen and from a lack of abundant transitions found in metals. The galaxies that had reasonable fits with Yggdrasil (6115 and GL-z12-1) may be indicative of a top-heavy IMF at high redshifts; however, we conclude this is generally not the case for the following reason.

No Balmer break exists in the younger stellar populations in Yggdrasil, but it is still not significant in the older available populations. This may contribute to the poor matching of Yggdrasil spectra by emcee, as a bump in emission at approximately 3600 Å in the rest frame is evident in the observed photometry. While the Pop II stellar populations available in Yggdrasil may have modelled such features more accurately, they do not have a top heavy IMF and would have been a redundant model, given what we have available in BPASS. There may not be an accurate enough top heavy IMF database publicly available to us to draw a conclusion whether this discrepancy in spectrum matching is due to the choice of IMF or a systematic error in the chosen spectra.

Compared to Yggdrasil, BPASS models more metal-enriched binary and SSPs and contains a strong Balmer break in its spectra. These features in conjunction with the applied absorption and attenuation functions allowed BPASS SEDs to match well with observed photometry. Almost every Yggdrasil model produced an AGN-dominated SED as shown in the middle panel of Fig. 13. Alternately, almost every BPASS model produced a stellar-dominated SED, except for 6115 which is completely AGN dominated.

As mentioned before, we did not include nebular emission (lines or continuum) in our SEDs. We recognize that this is a potentially important component of the spectra and we plan to include it in a future version of this work. To do this, we will process the stellar spectra with the radiative transfer code cloudy, similarly to how the Flexible Stellar Population Synthesis (fsps) code does. Byler et al. (2017) present useful information on the effect of the nebular emission and how they incorporate it with fsps using cloudy. They found that in the optical and NIR, nebular emission can contribute 30–50 per cent of the total flux at 10 Myr. Given that most of the optical band and NIR (4000 Å–5.3 |$\mu$|m) is redshifted beyond the NIRCam sensitivity, this would not have a significant effect on our results. At ages greater than 10 Myr, the nebular component for SSPs is negligible (see Fig. 12). However, for a constant star formation rate, the nebular emission contributes approximately 20 per cent or less to the total flux at wavelengths of 900–4000 Å (see Fig. 13 therein). By not including nebular emission with the BPASS spectra, we are slightly underestimating the flux. Correcting for this would increase the flux for a given stellar mass, and therefore reduce the stellar mass needed to fit the photometry. This may be the primary reason we see slightly higher stellar masses than A22 and H22.

6 CONCLUSIONS

We have developed a galaxy SED model tailored to matching high redshift photometry from JWST. We use this model to fit spectra to the JWST/NIRCam photometry of 12 high-redshift candidate galaxies from the CEERS and GLASS surveys. In comparison with the early JWST results we compared against (e.g. A22, H22, N22), our model add a few possibly important components to the spectral modeling.

First, we utilize high-redshift star formation history models extending to z = 15, which agree well with cosmological simulations. Secondly, we include AGN spectra in our fitting while many other works assume a stellar-only SED. Thirdly, we test multiple IMFs including a top-heavy IMF that uses the stellar spectra of Pop III.2 stars. The other two IMFs we tested include stellar spectra both with and without stellar binaries. These additions represent a more realistic environment for these very young galaxies assembling during the Epoch of Reionization and could differ from more evolved and massive galaxies at lower redshifts.

Our goal was to determine if the interplay of these features, particularly the inclusion of an AGN, could explain the unexpectedly large stellar masses and bright UV luminosity of galaxies from early JWST photometry. We found the following common trends in our sample:

• Our pipeline was able to fit the JWST photometry moderately well with values of reduced χ² comparable to other works. Although almost all our SEDs were dominated by the stellar component while the AGN was of little importance.

• Our pipeline was able to fit the JWST photometry moderately well with values of reduced χ² comparable to other works. Although almost all our SEDs were dominated by the stellar component while the AGN was of little importance.

• Only the top-heavy IMF given by the Yggdrasil models produced significant AGN contributions, but generally fit the photometry worse than the standard IMF for regular stars. The only exceptions were GL-z12-1 and 6115, where the top-heavy IMF was the best model, but only marginally. In general, models with a top-heavy IMF consistently produced lower stellar masses and larger BH fractions.

• There was no significant variation in halo mass, stellar mass, BH mass, or χ² when we included stellar binaries. This could indicate that the BPASS spectra do not differ significantly in our wavelength range of interest whether it is a single or binary stellar population. Additionally, because newly formed stars dominate the UV spectrum, the exponential halo growth rate was of little importance.

• Ultimately, we predict larger (smaller) best-fitting stellar masses than A22 and H22, yet still within the uncertainties, for eight (four) of the 12 galaxies. In these four cases, this could be caused by the presence of an AGN component, but emission-line diagnostics will be needed for their confirmation.

Although we could not confirm our initial hypothesis of an AGN component that moderates the stellar mass estimates, our approach can place constraints on the host halo growth history, star formation history, stellar populations, and central BH mass of these |$z \lower0.3em\hbox{$\, \buildrel\gt \over \sim \, $}10$| galaxies observed with JWST in its first data release. While photometry of these distant galaxies places large uncertainties on the intrinsic galactic properties, the first JWST/NIRSpec detections of |$z \lower0.3em\hbox{$\, \buildrel\gt \over \sim \, $}10$| galaxies are just being published, placing stronger constraints on their redshifts, stellar masses, metallicities, and dust attenuation (e.g. Curtis-Lake et al. 2022). The first six months of JWST high-redshift observations have been truly groundbreaking, and we, as a community, will continue to push the high-redshift boundary further uncovering the nature and origin of the first galaxies and their central BHs.

ACKNOWLEDGEMENTS

The authors thank Yuichi Harikane and Hakim Atek for providing additional photometric data that were not present in their publications. The authors thank Sandrine Ferrans and Rohan Srivastava for their contributions to an earlier version of this work. JHW thanks Raffaella Schneider, Rosa Valiante, and Marta Volonteri for useful discussions. DS is supported by the Future Investigators in NASA Earth and Space Science and Technology fellowship (FINESST) award 80NSSC20K1540. SS is supported by the NASA FINESST fellowship award 80NSSC22K1589. This work is supported by The National Science Foundation (NSF) grants OAC-1835213 and AST-2108020 and NASA grants 80NSSC20K0520 and 80NSSC21K1053. The Markov chain Monte Carlo calculations were performed with NSF’s LRAC allocation AST-20007 on the Frontera resource in the Texas Advanced Computing Center. The figures in this paper were constructed with the plotting library matplotlib (Hunter 2007).

DATA AVAILABILITY

The data relevant to each model and the best-fitting corner plots are available in the article’s online supplementary material in the form of tables and figures. The source code and data underlying this article and the corner plots of the models that are not the best-fitting models will be shared on reasonable request to the corresponding author.

Footnotes

1

A_j = (0.0036, 0.0017, 0.0011846, 0.0009410, 0.0007960, 0.0006967, 0.0006236, 0.0005665, 0.0005200, 0.0004817, 0.0004487, 0.0004200, 0.0003947, 0.000372, 0.000352, 0.0003334, 0.00031644)

2

Note Madau’s definition of x_c in footnote 3 contains a typo. We give the correct definition here.

3

See supplementary data for all best-fitting corner plots.

APPENDIX A: APPARENT MAGNITUDES TABLE

Table A1 lists the apparent magnitudes for the best-fitting model parameters and their differences from the observed photometry, along with their uncertainties.

Table A1.

Magnitudes for the best-fitting model parameters and their differences with observed magnitudes.

Galaxy	Model	F150W	F200W	F277W	F356W	F410M	F444W	MUV
		ΔF150W (σobs)	ΔF200W (σobs)	ΔF277W (σobs)	ΔF356W (σobs)	ΔF410M (σobs)	ΔF444W (σobs)	ΔMUV (σobs)
10234	B_chab100_bin	28.22	27.49	27.75	27.98	–	28.04	–
		–0.019 (0.07)	0.364 (0.05)	–0.047 (0.05)	–0.241 (0.05)		–0.064 (0.05)
1514	A_chab100_sin	28.01	27.68	27.41	27.14	–	26.35	–
		–0.004 (0.05)	–0.037 (0.05)	0.101 (0.05)	–0.064 (0.05)		0.009 (0.05)
1696	A_chab100_sin	26.71	26.76	26.85	26.84		26.23
		0.044 (0.05)	0.021 (0.05)	–0.099 (0.05)	–0.012 (0.05)		0.033 (0.05)
2462	A_chab100_sin	26.43	26.24	26.09	25.84	–	25.06	–
		0.070 (0.05)	–0.025 (0.05)	–0.106 (0.05)	–0.027 (0.05)		0.080 (0.05)
2779	B_chab100_bin	27.78	27.87	28.03	28.07	–	27.62	–
		–0.017 (0.05)	0.084 (0.05)	–0.095 (0.05)	0.003 (0.05)		0.035 (0.05)
6115	A_yggdrasil_sin	28.17	27.75	27.89	27.97	–	28.08	–
		0.025 (0.06)	0.004 (0.05)	–0.064 (0.05)	0.034 (0.05)		0.006 (0.05)
6878	A_chab100_sin	27.02	26.99	26.99	26.92	–	26.32	–
		–0.011 (0.05)	0.027 (0.05)	–0.027 (0.05)	0.010 (0.05)		0.002 (0.05)
CR2-z12 − 1	B_chab100_bin	28.76	27.76	27.94	28.09	28.16	28.08	–20.00
		–0.345 (0.4)	0.162 (0.1)	–0.057 (0.2)	–0.107 (0.1)	0.062 (0.1)	0.081 (0.2)	–0.096 (0.1)
CR2-z16 − 1	A_chab100_sin	–	27.90	26.35	26.41	26.45	26.48	–21.88
			–0.001 (0.1)	–0.048 (0.1)	0.008 (0.1)	0.151 (0.1)	–0.124 (0.1)	0.024 (0.1)
GL-z12 − 1	B_yggdrasil_sin	29.01	26.81	27.00	27.12	–	27.16	–21.02
		–0.091 (0.1)	0.111 (0.1)	0.103 (0.1)	–0.083 (0.1)		–0.043 (0.1)	–0.025 (0.1)
GL-z9 − 1	B_chab100_bin	26.93	26.62	26.61	26.59	–	26.30	–20.95
		0.025 (0.1)	0.023 (0.1)	0.009 (0.1)	–0.011 (0.1)		0.004 (0.1)	–0.054 (0.1)
GL-z9 − 2	B_chab100_bin	28.11	27.97	28.14	28.25	–	27.95	–19.68
		–0.086 (0.1)	0.175 (0.1)	–0.062 (0.1)	–0.046 (0.1)		0.050 (0.2)	0.018 (0.1)

Note. First row: Magnitudes for best fit. Second row: Difference between model magnitudes and observed magnitudes and the observed standard error. Observed values are from H22 and A22. Italicized numbers indicate where our filter magnitudes lie outside the observed standard error.

Author notes

These authors contributed equally to this work.