Are the ultra-high-redshift galaxies at $z>10$ surprising in the context of standard galaxy formation models?

The semi-analytic model (SAM) developed by the Santa Cruz group is a versatile galaxy formation model that implements a suite of carefully curated physical processes to track the formation and evolution of galaxies via dark matter halo merger trees (Somerville & Primack, 1999; Somerville et al., 2008, 2012; Somerville et al., 2015, 2021; Popping et al., 2014). We refer the reader to these papers for a full description of the model components and Yung et al. (2022) for a schematic flowchart of the internal workflow of the model. The well-established model has been shown to reproduce the observed evolution of distribution functions of galaxy properties at $0<z<6$ (Somerville et al., 2015) and $4<z<10$ (Yung et al., 2019a, b), with free parameters in the model only calibrated to a subset of observational constraints at $z\sim 0$.

The Santa Cruz SAM includes a suite of standard physical processes that are also adopted by other semi-analytic models or incorporated in cosmological hydrodynamic simulations. These processes include cosmological accretion and atomic cooling, suppression of gas accretion after reionization by the meta-galactic photoionizing background, star formation and stellar-driven winds, chemical evolution, black hole feedback, and mergers. The cold gas disc is partitioned into atomic, molecular, and ionized components depending on the gas-phase metallicity of the disc using fitting functions based on numerical hydrodynamic simulations by Gnedin & Kravtsov (2011, denoted as GK). In addition, the model adopts a H₂-based star formation recipe. In a similar spirit as the Kennecutt-Schmidt star formation (SF) relation where SFR is proportional to the surface density of cold gas (e.g. Kennicutt, Jr., 1998), we adopt a SF relation that scales with the surface density of molecular hydrogen (e.g. $\Sigma_{\text{SFR}}\propto\Sigma_{\text{H${}_{2}$}}^{\alpha}$), where the slope $\alpha$ steepens from 1 to 2 above a critical value of $\Sigma_{\text{H${}_{2}$}}$. This steepening is motivated by observations (Sharon et al., 2013; Rawle et al., 2014; Hodge et al., 2015; Tacconi et al., 2018) as well as theory (Ostriker et al., 2010), and has been shown to be essential in order for the SAM to reproduce the observed evolution in the $z>4$ UV luminosity functions, stellar mass functions and star formation rate distribution functions (Somerville et al., 2015; Yung et al., 2019a).

The model parameters are calibrated as shown in Gabrielpillai et al. (2022), such that the SAM outputs match the observed $z\sim 0$ stellar mass function (Baldry et al., 2012; Bernardi et al., 2013; Moustakas et al., 2013), stellar-to-halo mass ratio (Rodríguez-Puebla et al., 2017), cold ISM gas fraction vs. stellar mass (Calette et al., 2018; Catinella et al., 2018), stellar metallicity (Gallazzi et al., 2005; Kirby et al., 2011), and $M_{\text{BH}}$–$M_{\text{bulge}}$ relation (McConnell & Ma, 2013; Kormendy & Ho, 2013). These physical parameters are not re-‘tuned’ at higher redshift, and we use identical values of the parameters in this study. Note that these parameters are slightly different from the ones used in the Yung et al. JWST forecast paper series, due to the switch from Extended Press-Schecter based merger trees to N-body based trees.

The full star formation and chemical enrichment history is convolved with publicly available results from the stellar population synthesis (SPS) model Binary Population and Spectral Synthesis (BPASS¹¹1https://bpass.auckland.ac.nz/, v2.3; Eldridge et al., 2017; Byrne et al., 2022) to produce a high-resolution SED for each galaxy. The rest-frame UV magnitude is calculated from the SED with a tophat filter centred around 1600Å. Yung et al. (2020a) and Yung et al. (2020b) explored various SPS models and have shown that the binary star populations in BPASS yield a higher ionizing photon budget which may help reproduce the reionization history implied by various constraints derived from intergalactic medium and cosmic microwave background observations (see Robertson et al. 2013 and Robertson et al. 2015 for a thorough review on the various types of constraints and how they are derived). However, at rest 1600Å the BPASS models yield very similar predictions to single star models such as those of Bruzual & Charlot (2003), and we do not expect our predictions to be very sensitive to the assumed stellar population models. The version of the BPASS model adopted in this work assumes a fiducial broken power law stellar Initial Mass Function, with an upper slope $a_{1}=-1.30$ between $0.1$ – $0.5$  M_⊙ and a lower slope $a_{2}=-2.35$ between $0.5$ – $300$  M_⊙. In this work, we utilize SEDs with stellar metallicity spanning $Z=10^{-5}$ to 0.040. We note that dust attenuation is omitted in this work, as we wish to obtain the most optimistic predictions for the ultra-high-$z$ galaxy populations. In addition, we include only stellar continuum emission, and neglect nebular emission, although we will include the latter in a forthcoming work.

Gadget at Ultrahigh Redshift with Extra-Fine Timesteps (gureft, pronounced graft) is a suite of dark matter-only cosmological simulations designed to have mass and temporal resolution sufficient to capture the merger history of dark matter halos at extreme redshifts (Yung et al., 2023). The merger trees extracted from this simulation suite will ultimately be grafted together to paint a holistic picture of the earliest emergence of halos that are massive enough to host galaxies, which are not available from most publicly accessible cosmological simulations due to insufficient mass resolution and sparsely spaced snapshots (see Table 1). The suite consists of four simulated volumes, with 5, 15, 35, and 90 Mpc $h^{-1}$ on a side, containing $1024^{3}$ dark matter particles, with mass resolutions specified in Table 1. These boxes are labelled gureft-05, gureft-15, gureft-35, gureft-90 accordingly. The initial conditions for the suite are generated with the MUlti-Scale Initial Conditions (Music, Hahn & Abel, 2011) code, with second-order Lagrangian perturbation theory (2LPT) at $z=200$. The evolution of dark matter is carried out with the publicly available version of Gadget-2 code²²2https://wwwmpa.mpa-garching.mpg.de/gadget/ (Springel, 2005) down to $z=6$.

For each simulated volume, 170 snapshots are saved between $z=40$ to 6 with spacing roughly one-tenth of the typical halo dynamical time. For comparison, the Bolshoi-Planck simulation only has 30 snapshots and the IllustrisTNG simulations (Nelson et al., 2019) has 14 snapshots in the same redshift range. Halo and merger tree catalogues are generated with the rockstar (Behroozi et al., 2013a) and consistent-tree codes (Behroozi et al., 2013b), based on the virial mass definition of Bryan & Norman (1998).

Fig. 1 shows the halo mass functions (HMFs) between $z\sim 6$ to 20 from the four gureft boxes (colour-coded). We also show predictions from the Bolshoi-Planck simulations at $z=6$, 10, and 15 for comparison, as well as the $z=0$ halo mass function to guide the eye. We note that no mass cutoff is applied in the making of these HMFs. Thus, the turnover observed on the low-mass end reflects where the halo populations become incomplete due to the mass resolution of the simulation, which is typically $\sim 100$ times larger than the mass of DM particles. Tabulated HMFs constructed by combining well-resolved halos from the four gureft volumes are available in Appendix A. For more details on the gureft simulations, as well as a comparison of the halo mass functions and other properties with past work and commonly used analytic models and fitting functions, please see Yung et al. (2023).

One-point distribution functions are a frequently used summary statistic for characterizing and comparing galaxy populations from observations and simulations. In Fig. 2, we present the stellar mass functions (SMFs) for galaxies across the full mass range that we can resolve between $z=8$ to 17 by combining the results from the four gureft volumes. One can see from the excellent agreement in the regions where the different boxes overlap that, as shown previously in Gabrielpillai et al. (2022), the convergence is excellent (i.e. the SAM results are insensitive to the resolution of the merger trees as long as the merger history is well resolved). These results are compared to SMFs derived from past Hubble and Spitzer observations (Song et al., 2016). The flattening or turnover seen at the low-mass end of gureft-35 and gureft-90 is caused by incompleteness due to the mass resolution limit of the underlying cosmological simulations. The gureft-05 and gureft-15 boxes are designed to resolve low-mass halos and their merger histories, and thus here the turnover reflects the halo mass limit where atomic cooling becomes inefficient (we do not include processes that can cool gas below $10^{4}$ K). This has been previously explored using the extended Press-Schechter (EPS)-based merger trees (e.g. Somerville & Kolatt, 1999) and a grid of halo masses as presented in Yung et al. (2019b), which is also shown in this comparison. As previously shown, our model predictions are in very good agreement with the observationally derived stellar mass functions at $z=8$. We refrain from showing stellar mass functions for higher redshifts or for any of the recent JWST based studies, as the uncertainties on these stellar mass estimates are still extremely large.

Fig. 3 presents the rest-frame UV luminosity functions (UV LFs) between $z=8$ to 17 for galaxies simulated by the SC SAM. These predictions are compared to the recent JWST observations by Finkelstein et al. (2022b, 2023b, 2023a) from CEERS (Bagley et al., 2023b), Leung et al. (2023) from NGDEEP (Bagley et al., 2023a), Casey et al. (2023) from COSMOS-Web (Casey et al., 2022), Donnan et al. (2022), Harikane et al. (2023), and Bouwens et al. (2023a). In addition, we include pre-JWST measurements from Bouwens et al. (2014) and Finkelstein et al. (2015), and a compilation of results from Finkelstein (2016). These include observations from lensed fields from Livermore et al. (2017) and Bouwens et al. (2022), as well as the latest results derived from past Hubble and Spitzer Space Telescope observations presented by Oesch et al. (2018), Stefanon et al. (2019), Bowler et al. (2020), and Finkelstein et al. (2022a). We show the $z\sim 17$ luminosity function estimates of Bouwens et al. (2022) and Harikane et al. (2023) as open points to illustrate how these published results compare with our model predictions. However, the single galaxy from the CEERS field (first discovered by Donnan et al., 2022) that went into the Bouwens et al. (2022) estimate has now been shown to be at $z\sim 5$ (Arrabal Haro et al., 2023). The Harikane et al. (2023) estimate is based on two galaxies, the same CEERS galaxy now known to be at lower redshift, and a candidate from the Stephan’s Quintet field (S5-z16-1). Similar to the CEERS object, S5-z16-1 has a bimodal photometric redshift solution with one peak at $z=16.4$ and another at $z=4.9$. The UV luminosity implied by the $z=16.4$ solution would be extremely high ($M_{\rm UV}=-21.6$). One can see that even one object of this luminosity at this redshift would be in very strong tension with our model predictions.

In Appendix A, we provide the tabulated data for the SMFs and UV LFs by combining the predicted galaxies across the four simulated volumes within the mass (and luminosity) range where the galaxy populations are complete. This approach is similar to the one adopted by Vogelsberger et al. (2020) to produce JWST predictions from the suite of IllustrisTNG simulations (e.g. Nelson et al., 2018; Pillepich et al., 2018; Nelson et al., 2019).

As also shown in previous works Yung et al. (2019b), the rest-UV luminosity functions predicted by our models are in good agreement with observations at $z=8$–10. Note that the predictions shown here do not include any dust attenuation, so the small excesses in our predicted UV LFs at the bright end at $z\sim 8$ and 9 could be cured by invoking a moderate amount of extinction (see Yung et al. (2019b), which did include modelling of dust attenuation). We also note that the free parameters used in this work are based on the re-calibration presented in Gabrielpillai et al. (2022). In particular, the slope for the stellar feedback-mass outflow relation (see Yung et al. (2019a) for full description) has been slightly increased from $\alpha_{\text{rh}}=2.8$ to 3.0, resulting in higher mass outflow rates in low mass halos, and hence lower masses and luminosities for the galaxies hosted by those halos. As a result, the faint-end (low mass-end) slopes for galaxies presented in this work are slightly shallower than those presented in previous works. This is expected and is in agreement with the controlled experiment shown in Yung et al. (2019a, b). By $z=11$, we see a moderate discrepancy between our predicted UV LF and the observational measurements, with the number densities predicted by our models an order of magnitude or more lower than the observational estimates. The discrepancy is about a factor of 30 at $z=13$.

The number density constraints estimated for observed galaxies presented in previous figures suffer from various uncertainties, including photometric errors, redshift errors, and field-to-field variance due to Poisson sampling statistics as well as galaxy clustering (field-to-field variance). Photometric errors can cause a so-called “Eddington bias”, where the much more numerous intrinsically low-luminosity galaxies are “scattered” into the higher luminosity bins, causing a flattening of the high-luminosity slope of the luminosity function). The full uncertainty budget has not been fully accounted for in the published error bars. Similarly, the simulated outputs also carry uncertainties due to cosmic variance associated with our relatively small simulated volumes.

In order to better assess the level of tension between our theoretical predictions and the observational measurements, in Fig. 4 we show a “zoomed in” version of the $z=11$ and $z=13$ UV LF comparisons, where we focus on the luminosity range of the detected galaxy populations. Here we have selected the luminosity range where the simulated galaxies are best captured for each gureft box and combine them. For instance, at $z=12$ and 13, galaxies with $M_{\text{UV}}\leq-18$ are sourced from gureft-90, and galaxies between $-18\leq M_{\text{UV}}\leq-13$ are from gureft-35.

Note that we do not attempt to correct for redshift uncertainties (which may shift the galaxies in redshift within these bins) or catastrophic redshift errors (which may cause galaxies that are at much lower redshifts to be counted in this ultra high-z population), as these corrections remain very uncertain and are not straightforward to model. We experimented with including photometric noise in our predicted luminosity functions, and find that it introduces a minor source of uncertainty in this comparison, so we do not show it here.

We explore the potential impact on the number density of simulated galaxies due to cosmic variance. We show cosmic variance estimates from the cosmological simulation bluetides (Feng et al., 2016; Bhowmick et al., 2020), which is 400 h^-1 on a side, $\sim 88$ and $\sim 1492$ times larger than the box sizes of gureft-90 and gureft-35, respectively. We also show estimates from the cosmic variance calculator of Trenti & Stiavelli (2008). A detailed discussion of how we obtained these estimates for each panel, and of uncertainties associated with them, can be found in Section 4. We have chosen the CEERS survey as a representative field for these estimates. We show a representative “total error” estimate as the sum in quadrature of the quoted error from FL23 and the larger of our cosmic variance error estimates. In their $z\sim 11$ bin, the uncertainty due to Poisson sampling and photometric errors combined (a fractional error of $\sim$0.44 at $M_{\rm UV}=-20$) is comparable to that from cosmic variance (fractional error of $\sim 0.2-0.3$). At $z\sim 14$ and $M_{\rm UV}=-20$, the quoted fractional error is 100%, and the uncertainty due to cosmic variance is $\sim 0.3-0.77$.

In Fig. 4, we also show the impact of theoretical uncertainties in the stellar initial mass function. For our fiducial calculations, we have adopted an IMF that is similar to the standard Chabrier (2003) IMF which is typical for galaxies in the nearby Universe. However, there are multiple reasons to think that the IMF might be top heavy in the early Universe (as discussed in more detail in Section 4). It has been suggested that stellar populations could be 2-10 times brighter in the UV if the IMF is top heavy (e.g. Raiter et al., 2010). The blue shaded regions in Fig. 4 show what happens to our predictions if we simply boost the luminosity of each galaxy in our simulation by a factor of 2–10. One can see that, at $z=11$–13, a fairly modest “IMF boost” of a factor of four would bring our model predictions into reasonable agreement with the observational estimates.

Although our model includes merger-triggered bursts, there could be short timescale star formation stochasticity due to ISM processes that we do not resolve or attempt to represent in our models. Given there are fewer bright objects, adding stochasticity causes more galaxies to scatter into the bright end of the UVLF, in an effect similar to Eddington bias. In order to explore the effect of UV stochasticity in a simple way, we add a stochastic component to the predicted UV magnitudes in post-processing. The stochastic component is assumed to follow a Gaussian distribution centred at zero for specified values of the standard deviation $\sigma_{\rm UV}$. We repeat the process of adding the stochastic component and computing UV LFs for 2,500 independent realizations and report the median of the outcome. In Fig. 5, we show the effect of $\sigma_{\rm UV}=0.5$, 1.0, 1.5, and 2.0; and compare them to the same observational constraints shown in Fig. 3. We find that a relatively large value of $\sigma\simeq 2$ is required to match the observations at $z\sim 11$–13.

In Fig. 6 we show how the composite stellar mass function and rest-UV luminosity function evolve from $z=6$ to 17. These are compiled by combining all of the gureft volumes. This figure highlights the relatively rapid predicted evolution of both the SMF and UVLF. This appears to be in tension with the current observational results, which suggest weaker evolution over cosmic time at the highest redshifts. Fig. 7 shows the same quantity but sliced in a different way – this shows the number density in a bin of $M_{\text{UV}}$ as a function of redshift. The number density of luminous galaxies is predicted to evolve more rapidly than that of lower luminosity galaxies, as has also been predicted (and observed) at lower redshifts (sometimes called ‘downsizing’).

Scaling relations are conditional relationships between two physical or observable properties of galaxies, and they are important diagnostics of the physics that shapes galaxy properties. In addition to rest-frame $M_{\text{UV}}$ and $M_{*}$, physical properties such as star formation rate (SFR) and UV spectral slope ($\beta_{\text{UV}}$) can be inferred from the measured multi-band photometry with SED fitting techniques (e.g. Johnson et al., 2021; Carnall et al., 2018; Iyer et al., 2019; Boquien et al., 2019). In this subsection, we compare outputs from the SC SAM to scaling relation constraints inferred by observations. In these studies, physical properties such as stellar masses are obtained through SED fitting. It is important to keep in mind that these techniques rely on a complex modelling procedure, and the estimates can have large uncertainties due to corresponding uncertainties in the underlying simple stellar population models, and assumed priors regarding star formation history, chemical evolution, dust, etc. (e.g. Papovich et al., 2011; Tacchella et al., 2022b). Furthermore, these inferred properties are subject to uncertainties in distance.

Fig. 8 shows the $M_{*}$–$M_{\text{UV}}$ relation from our predicted SAM galaxies, where individual galaxies are shown by the data points and the median and 16th and 84th percentiles are shown with the cyan line and shaded area. These results are compared to results presented by Whitler et al. (2022) derived from NIRCam imaging in the CEERS field and the flexible Bayesian SED fitting code Prospector with a built-in ‘continuity’ prior. In addition, we include constraints derived from past Hubble and Spitzer observations by Tacchella et al. (2022b). These quantities are also based on the Prospector tool, assuming a ‘bursty’ version of a continuity prior. We note that while these galaxies span a wide range of redshift ($8<z<11$ for Whitler et al. (2022) and $9<z<11$ for Tacchella et al. (2022a)) with relatively large uncertainties in photometric redshift, we show all observed galaxies across all panels in Fig. 8 for discrete redshift snapshots at $z=8$, 11, and 13 from our simulations. The median of the $M_{*}$–$M_{\text{UV}}$ relation are overlaid in the last panel. It is intriguing that the model predictions for the relationship between stellar mass and rest-UV luminosity are in good agreement with the observational estimates at $z=8$ and $z=11$. In the $z=13$ bin, there are very few model galaxies in the luminosity range of the observations, but the observed galaxies seem to lie along a plausible extrapolation of the simulated relation for lower luminosity galaxies. Also notable is the rather weak evolution in the $M_{*}$–$M_{\text{UV}}$ over this redshift interval from $8<z<13$, with galaxies predicted to be slightly more luminous for a given stellar mass at higher redshift.

Fig. 9 shows a comparison between the modelled and observed $M_{\text{UV}}$ and rest-UV slope $\beta_{\text{UV}}$. For the simulated galaxies, we compute $\beta_{\text{UV}}=-0.4(M_{\text{FUV}}-M_{\text{NUV}})/\log(\lambda_{\text{FUV}}/\lambda_{\text{NUV}})-2$ with magnitude in far UV (FUV, centred at 1600 Å) and near UV (NUV, centred at 2300 Å) bands. These results are compared to results from JWST observations reported by Topping et al. (2022) and Cullen et al. (2023). In addition, we also show results from Hubble and Spitzer observations reported by Tacchella et al. (2022b). The median of the $M_{\text{UV}}$–$\beta_{\text{UV}}$ relations from all three redshift snapshots are overlaid in the last panel. Some of the observed galaxies have significantly redder values of $\beta_{\text{UV}}$ than our models predict, which could be due to the presence of dust in some of these galaxies. Additionally, our models assume a metallicity floor of $10^{-3}$ Z_sun, and the BPASS spectral synthesis models do not include Pop III type tracks or SEDs. The extremely blue slopes reported by Topping et al. (2022) are an exciting hint that we may be witnessing the formation of Pop III stars in some of these objects. The broad range of UV slopes seen in the observations suggests that the galaxy population could be a mixture of Pop II and Pop III stars at these epochs.

Fig. 10 shows the relation between $M_{*}$ and $\cal{K}_{\text{UV}}$. The factor $\cal{K}_{\text{UV}}$ is an empirically measured conversion factor that relates the rest-frame UV luminosity to the star formation rate $\text{SFR}={\cal{K}_{\text{UV}}}\times L_{\nu}(\text{UV)}$, where we consider the FUV luminosity and SFRs averaged over 100 Myr for this comparison. Similar to previous figures, we show the individual galaxies with data points and the median and 84th and 16th percentiles with cyan lines and shaded regions. We compare our simulations with values reported by Madau & Dickinson (2014), ${\cal{K}}_{\text{UV}}=0.72\times 10^{-28}$ M_⊙ yr^-1 erg^-1 s Hz and ${\cal{K}}_{\text{UV}}=0.88\times 10^{-28}$ M_⊙ yr^-1 erg^-1 s Hz as adopted by Kennicutt, Jr. (1998). We note that these values are adjusted from a Salpeter (1955) to Chabrier (2003) IMF by multiplying by a factor of 0.63 (Madau & Dickinson, 2014). Since $\cal{K}_{\text{UV}}$ is effectively the efficiency of producing observed FUV light for a given SFR, it is expected that its value will be sensitive to star formation and chemical enrichment histories. The medians of the $M_{*}$–$\cal{K}_{\text{UV}}$ relations at $z=8$ to 17 are overlaid in the last panel. As expected, galaxies in the early universe contain younger stellar populations and are more efficient at producing UV light. Similarly, Fig. 11 shows the relation between metallicity of cold ISM gas ($Z_{\text{cold}}$) and $M_{*}$.

Additional scaling relations with dark matter halo mass are presented in Appendix B.

Halo assembly and merger histories form the backbone of galaxy formation simulations, especially for the semi-analytic modelling approach. Cosmological simulations are subject to tension between simulated volume and mass resolution, subject to available computational resources. While extremely large volume is required to capture overdense environments where the very rare, massive halos are found, it is also critical to have sufficient mass resolution to properly resolve halos’ assembly histories. While this tension is always present in modelling galaxy formation (since we can observe a much larger dynamic range of galaxy environments than we can currently simulate), there has been less attention directed at creating a suitable suite of N-body simulations for studying structure formation in the “ultra-high redshift” Universe ($z>6$) — up until now. Establishing connections between halos identified across multiple simulation snapshots and constructing the full merger histories of halos also requires high temporal resolution. However, most state-of-the-art cosmological N-body simulations store a limited number of snapshots at ultrahigh redshift (likely a thoughtful decision made to optimize storage usage), and these are insufficient to reconstruct merger histories of these halos. For these reasons, many publicly available cosmological simulations do not deliver halo merger trees that can be used by SAMs to explore predictions for galaxies in the exciting redshift frontier that has recently been opened up by JWST ($10<z<17$). We designed the gureft simulation suite specifically to overcome these barriers. In a companion paper (Yung et al., 2023), we will present updated fitting functions describing the halo mass function and other useful halo property distributions in the ultra-high redshift universe from the gureft suite along with a comparison with previous halo mass functions that have been used in the literature.

We add that the mass range probed by the higher resolution boxes, gureft-05 and gureft-15, may be probing a regime that can be affected by uncertainties in the physical nature of dark matter. In this work, we assumed ‘vanilla’ $\Lambda$CDM. If the dark matter has a wavelike or “fuzzy” nature, this could substantially impact the low-mass end slope of the HMF (Menci et al., 2017; Hui et al., 2017; Kulkarni & Ostriker, 2022, see Hui 2021 for a review). Additionally, if the dark matter is warm or self-interacting, this could also impact the abundances of low-mass halos at early times and their internal structures and formation histories (e.g. Menci et al., 2016; Lovell et al., 2018; Adhikari et al., 2022). The nature of dark energy could also impact structure formation at early times. Models with Early or Dynamical Dark Energy can predict higher abundances of low-mass halos at early times (Menci et al., 2020; Klypin et al., 2020; Menci et al., 2022).

A comprehensive comparison between our predictions and those of the large number of other published models and simulations is beyond the scope of this work, but we comment briefly here on the general theory landscape and on a few key comparison points. First, we comment on the predictions of the “Semi-analytic Forecasts for JWST” paper series (Yung et al., 2019a, b, 2020a), which adopted the same semi-analytic model presented here with slightly different parameters, but implemented within merger trees constructed using the Extended Press-Schechter formalism. We show in Fig. 2 and Fig. 3 that there is very good agreement between the Yung et al. predictions and those of the present work from $z\sim 8$–13, with the small differences on the faint end of the UVLF being largely due to the re-calibration of the stellar feedback parameters as discussed in Section 3. At $z>13$, our new and more robust calculations based on our new gureft suite predict number densities that are higher by an order of magnitude or more than the EPS based calculations. We note that this has little to no impact on any of our published results or conclusions, as we deliberately focussed on redshifts $z\lesssim 12$ specifically because we did not trust the EPS merger trees at higher redshifts.

Comparisons with other model predictions are shown in Harikane et al. (2023) for $z=12$ and $z=14$–16 (see their Figure 16) and FB23 (see their Figure 14) for $z\sim 8$–14. The general impression from these comparisons is that our predicted number densities at $z\sim 8$–10 are generally consistent (within a factor of a few) with those of large cosmological hydrodynamic simulations such as SIMBA (Wu et al., 2020), Millennium-TNG (Kannan et al., 2023), and THESAN (Kannan et al. 2022; see also Yung et al. (2019b) for a detailed comparison of our model predictions with simulation predictions available at that time out to $z=10$). The hydrodynamic simulation FLARES and the semi-analytic model DELPHI predict number densities at $m<28.5$ and $z\sim 11$ that are higher than ours by about an order of magnitude and a factor of 3, respectively (F23). Our predicted number densities are also very similar to those of the semi-empirical model UniverseMachine (see Behroozi et al., 2020, for a detailed comparison). In summary, the significant under-prediction of $z\gtrsim 10$ galaxies observed by JWST in physics-based models is a fairly generic result and certainly not specific to our models.

With robust halo catalogs and merger trees spanning halos over a mass range of $\sim 10^{5.5}$–$10^{12.5}$ from $z\sim 20$–6, we are able to compute predictions for the properties of ultra-high redshift galaxies over an unprecedented dynamic range. Our predictions are based on the well-established Santa Cruz semi-analytic models, which adopt a very standard set of physical processes, including atomic cooling, partitioning of galaxies’ ISM into multiple phases, star formation and stellar feedback, chemical evolution, and black hole seeding, growth, and feedback. Although these models contain simplified representations of these physical processes, we have shown in the past that the predictions for many key galaxy properties are remarkably similar to those from more detailed (and much more computationally expensive) numerical hydrodynamic simulations (Yung et al., 2019a; Gabrielpillai et al., 2022, 2023).

One of the puzzles of galaxy formation manifested in observations of nearby galaxies is why star formation is so inefficient and why galaxies and halos capture such a small fraction of available baryons (e.g. Moster et al., 2010; Behroozi et al., 2013c). The standard solution to this puzzle is that heating and winds from massive stars and supernovae eject baryons and make star formation inefficient in low mass halos, and similar processes associated with accreting black holes eject gas and prevent it from cooling in massive halos (e.g. Somerville & Davé, 2015). Our models, like all large volume cosmological simulations, contain parameterized phenomenological treatments of these key “feedback” processes, as it is infeasible to simulate the full range of scales and physical processes. These parameters have been tuned by hand to reproduce key global observable (or semi-observable) properties of galaxies in the nearby Universe, such as their stellar-to-halo mass relation, cold gas fractions, stellar metallicities, etc. There is absolutely no guarantee that the parameterizations of key physical processes such as stellar feedback scale correctly in terms of the physical properties that are actually relevant. Therefore it is actually a bit surprising that the models perform as well as they do out to such high redshifts ($z\sim 10$). In this section, we ask what we learn from the apparently growing tension between the model predictions and the observational estimates at $z\gtrsim 10$.

We note that we have verified in the past that under our current implementation, black hole feedback has a negligible impact on the galaxy populations at these very high redshifts, so the main processes shaping galaxy properties in our models are the rate of gas accretion and cooling, the efficiency with which cold gas can form stars, and the efficiency of stellar driven outflows, which eject gas from the ISM and suppress future star formation. Taking the observational estimates at face value for purpose of this discussion, we consider four hypotheses for why our models might be predicting lower number densities of galaxies than the observations suggest:

1. 1.

   The standard $\Lambda$CDM model does not predict early enough formation of halos massive enough to host the observed galaxies.

2. 2.

   Cooling is inefficient

3. 3.

   Star formation is inefficient

4. 4.

   Stellar feedback is too strong

In the above sub-section, we attempted to quantify the roles of different physical processes within the framework of our existing models. However, there are several other theoretical uncertainties that we have not touched on, which we summarize here.

We note that the metallicity of the circumgalactic medium impacts how rapidly it cools, so chemical yields and the efficiency with which metals are ejected from, and circulate back into, the CGM will also impact the predicted cooling rates. Similarly, the chemical enrichment of the ISM and the formation of dust (not modelled here self-consistently) will have a strong impact on the modelling of star formation. We emphasize that the star formation prescription we have adopted was based on radiation-hydrodynamic simulations that are more than a decade old at this point, and on local universe conditions. In reality, H₂ may simply trace the dense gas where star formation is efficient rather than being a pre-condition for star formation (e.g. Krumholz et al., 2011). At gas metallicities less than about 0.1 Z_⊙, additional drivers of the thermal state of the ISM become important, including reduced photoelectric heating, cosmic ray ionization heating, and H₂ cooling (Bialy & Sternberg, 2019).

A related issue is the stochasticity, or burstyness, of star formation. If the process of star formation at early times is highly stochastic, then galaxies could experience short periods of time when their UV luminosities are boosted due to starbursts. In fact, there is mounting observational evidence that ultra-high redshift galaxies did have extremely bursty star formation histories (e.g. Dressler et al., 2023; Endsley et al., 2023a; Looser et al., 2023). Several modelling studies have shown that very bursty star formation could help to increase the number of luminous objects that would be observable at very high redshift and reduce the tension between model predictions and ultra-high redshift observations (Mason et al., 2023; Shen et al., 2023; Sun et al., 2023a, b). We also find that adding a stochastic burst component in post-processing appears to be able to reconcile our fiducial model predictions with observations out to $z\sim 13$. However, in agreement with other studies (e.g. Shen et al., 2023), we find that a rather large stochastic component ($\sigma_{\rm UV}\sim 2$, where $\sigma_{\rm UV}$ is the root variance of a Gaussian random deviate in UV magnitude) is needed. This is much larger than the stochasticity produced in the high-resolution radiation-hydrodynamic cosmological simulations of Pallottini & Ferrara (2023), which yield typical $\sigma_{\rm UV}\sim 0.6$. We are therefore inclined to think that stochasticity alone is not responsible for resolving the discrepancy. However, weaker feedback, higher star formation efficiency, and bursty star formation are all expected to be consequences of the same physical changes in the ISM in ultra-high redshift galaxies. In future work, we plan to attempt to model these effects self-consistently in a physically motivated manner.

Another important process is the suppression of accretion and cooling by the meta-galactic photoionizing background after the Universe is reionized. We adopt a highly simplified treatment of photoionization feedback in our models, in which we assume that the Universe was reionized everywhere instantaneously at $z=8$. After reionization, the baryon fraction that accretes into low mass halos is reduced, according to the the results of Okamoto et al. (2008). However, in the real Universe, early forming, luminous objects in dense environments will create bubbles which will start to overlap and form larger reionized regions. We do not properly model the complex interplay between photo-ionization feedback from both local and meta-galactic radiation fields and stellar feedback. This will be increasingly important to model accurately as JWST begins to push to lower luminosity galaxies. We also note here that Harikane et al. (2023) propose that the high observed number density of ultra-high redshift galaxies could be explained by the lack of suppression of star formation via photo-ionization before reionization. Our models do not include any suppression of star formation by the UV background at $z>8$, yet they still underpredict the observed counts — therefore we suggest that this cannot be a complete explanation.

Another open question is how dust forms and is destroyed in the ISM at these early epochs. As noted above, dust is important for the chemistry and cooling of the ISM, and also impacts the observed SED of galaxies through attenuation and re-emission. Models and simulations have been presented that more self-consistently track the formation and destruction of dust through multiple channels (e.g. Popping et al., 2017; Vogelsberger et al., 2020; Dayal et al., 2022), but this physics is not incorporated into our current models. This will be an important ingredient to include as we update the realism of our modelling of star formation and the ISM in future work.

Another possibility is that some of the UV light in some of the ultra-high-z galaxies could be contributed by an accreting supermassive black hole (AGN), which would provide a luminosity boost (Pacucci et al., 2022). Our current models include black hole seeding and accretion, and do not predict a significant contribution from AGN at these redshifts, but this is subject to the details of our seeding and accretion modelling. Volonteri et al. (2023) investigated this possibility, finding that a massive black hole seed (such as might form through the “direct collapse” mechanism) growing at close to the Eddington rate could be detectable with JWST, and would be difficult to disentangle from star formation powered galaxies via the available JWST colours. Candidate AGN at ultra-high redshift have been reported from early JWST observations (e.g. Ono et al., 2023). However, in current samples, the majority of the detected objects have extended (rather than point source) morphologies, suggesting that an AGN is not the dominant source of luminosity (e.g. Harikane et al., 2023). Future JWST observations will help to clarify the contribution of AGN to the detected ultra-high-z population, and place constraints on models of black hole seeding and accretion.

Perhaps one of the largest uncertainties at ultra high redshifts is the properties of the stellar populations at these epochs. There is a large literature on the first stars (Pop III) which formed out of primordial gas and how their properties may have differed from present day stellar populations (e.g. Bromm et al., 1999; Abel et al., 2000; Zackrisson et al., 2011; Visbal et al., 2018). Although there is general agreement that the IMF was likely to have been more top-heavy in the early Universe, due both to the higher CMB temperature and lower metallicity, there is no detailed consensus on the expected Initial Mass Function (IMF) of Pop III stars, and even less on how the IMF may have evolved as the Universe is gradually polluted by metals (e.g. Pop III.2 stars). Clearly this will have a potentially profound impact on the galaxy SED, the chemical enrichment, and the stellar feedback — none of the existing cosmological simulations that have been compared with the JWST observations self-consistently incorporate these effects. Harikane et al. (2023) suggest that adopting a top-heavy IMF in early galaxies could help reconcile the theoretical predictions with observed number densities. They show that Pop III.1 and Pop III.2 populations with top-heavy IMFs modelled with the stellar population code YGGDRASIL (Zackrisson et al., 2011) can produce 3-4 times more UV photons for a given SFR than standard assumptions, in part due to the strong emission lines excited by massive stars. Raiter et al. (2010) similarly show the expected change in the ratio of UV luminosity to SFR for several IMF variants that might arise in the conditions typical of the high redshift Universe, also finding an expected increase of a factor of 2-10. This possibility has also been discussed by FB23, who showed using abundance matching that a factor of $\sim 2$ boost in the UV luminosities could reconcile the differences with models. We showed a simple representation of how a top-heavy IMF would affect our model predictions in Fig. 4, by simply boosting the luminosity of each galaxy in our models by a factor of 2–10. A boost of a factor of four brought our predicted number densities into good agreement with the observations. Clearly, modelling the formation of Pop III stellar populations and the transition through Pop III.2 to Pop II, and self-consistently incorporating appropriate stellar population models, is an important next step in modelling the ultra-high redshift Universe. Some early studies have found extremely blue UV slopes in ultra-high redshift galaxies, which could be an indication of the presence of extremely metal poor or even Pop III stars (Topping et al., 2022), although we note that these can also be indicative of high escape fraction (e.g. Endsley et al., 2023b; Jung et al., 2023; Mascia et al., 2023a, b). Emission line diagnostics from future JWST spectroscopic observations will be able to shed light on the physical nature of the UV radiation sources (Cleri et al., 2023).