Semi-analytic forecasts for JWST - I. UV luminosity functions at z = 4 - 10
L. Y. Aaron Yung Rachel S. Somerville Steven L. Finkelstein Gergö Popping Romeel Davé


firstpage–lastpage
Development of a SEMPA system for magnetic imaging with ns time resolution and phase-sensitive detection
    Mathias Kläui
    December 30, 2023
========================================================================================================




In anticipation of the upcoming deployment of the James Webb Space Telescope (JWST), we present high-redshift predictions by the well-established Santa Cruz semi-analytic model. We update the models by re-calibrating them after adopting cosmological parameters consistent with recent constraints from Planck. We provide predictions for rest-frame UV luminosity functions for galaxy populations over a wide range of M_UV from ∼-6 to ∼-24 between z = 4 – 10. In addition, we present the corresponding predictions for observed-frame galaxy number counts in different redshift bins in the full set of NIRCam filters. We provide predictions of the quantitative effect on these observables of varying the physical recipes implemented in the models, such as the molecular gas depletion time (star formation efficiency) scalings or the scalings of outflow rates driven by stars and supernovae with galaxy circular velocity. Based on these results, we discuss what may be learned about the physical processes that shape galaxy formation from JWST observations of galaxy number densities at different intrinsic luminosities. All data tables for the results presented in this work are available at <https://www.simonsfoundation.org/semi-analytic-forecasts-for-jwst/>.




galaxies: evolution – galaxies: formation – galaxies: high-redshifts – galaxies: star formation – galaxies: statistics





§ INTRODUCTION


The soon-to-be-launched James Webb Space Telescope (JWST) will possess the unprecedented infrared sensitivity and spatial resolution required for detecting faint, distant galaxies that are extremely difficult or impossible to detect with any current facilities. These observations will provide significant insights into the statistical properties of the galaxy population near cosmic dawn. This is of great interest for determining whether the physical processes that shape galaxy formation are very different in the early universe from locally, as well as for constraining which objects are responsible for reionizing the Universe. 

Hubble Space Telescope (HST) has detected nearly 2000 galaxy candidates at high redshifts (z ∼ 4 – 10) from both blank and gravitationally lensed fields <cit.>. Additionally, brighter objects have also been discovered with ground-based facilities, such as the United Kingdom Infra-Red Telescope (UKIRT) and the Visible and Infrared Survey Telescope for Astronomy (VISTA) <cit.>. These observations have provided constraints on the space density of relatively bright galaxy populations up to z ∼ 10 <cit.>. However, the constraints on the faintest populations currently rely solely on fields that are lensed by massive foreground galaxy clusters, and the corrections for magnification are quite uncertain <cit.>.
JWST will probe much deeper down the luminosity function in unlensed fields, as well as providing more secure redshift measurements for high-z candidates that are currently selected by the Lyman-break technique. Many members of the community are currently engaged in planning the optimal strategies to take advantage of the limited lifetime of JWST to achieve these breakthroughs in our understanding of the very early Universe. Theoretical predictions may be able to aid in designing the observing strategy and trade-offs in imaging area, depth, and wavelength coverage. Perhaps more importantly, it is critical to assess the uncertainties on current theoretical predictions and the plausible range in galaxy properties at these extreme redshifts.

With the cosmological parameters that govern early structure formation fairly well constrained by experiments such as the Wilkinson Microwave Anisotropy Probe <cit.> and Planck satellite <cit.>, combined with other constraints from Baryon Acoustic Oscillations, Supernovae, and weak lensing <cit.>, the major uncertainties in forecasting the abundance of galaxies in the early Universe arise from our lack of a rigorous theory for how dense, cold molecular clouds form, how stars form within these clouds, and how thermal energy, momentum, and radiation from stars and supernovae affect the subsequent efficiency of star formation (“stellar feedback”; see <cit.>, <cit.>, and <cit.>). Although feedback from radiation and jets produced by accretion onto supermassive black holes is likely to be critical for regulating galaxy formation in the lower redshift Universe ( and references therein), at very early times z ≳ 6 most galaxies probably have not had time to form massive black holes in their nuclei, so it is likely (although not certain), that this form of feedback is sub-dominant at these epochs. 

Some theorists have suggested that star formation remains very efficient out to very high redshifts <cit.>, while others have suggested that the low metallicity environments will make it difficult to form molecular hydrogen and lead to inefficient star formation <cit.>. The efficiency of star formation in low-mass halos at early times is also critical to determining which objects reionized the Universe, and how reionization proceeded in space and time. Numerous studies have shown that the already detected galaxy population is likely insufficient to reionize the Universe by z∼ 6–10 as required by observations <cit.>, and it is generally assumed that the shortfall in the ionizing photon budget is made up by faint galaxies that form in low-mass halos <cit.>. Many published theoretical predictions for reionization assume a fixed efficiency for converting baryons into stars in halos <cit.>. However, it is known that in the very nearby Universe, star formation is extremely inefficient in low-mass halos, leading to questions about whether the local and high redshift observations are in tension. 

Computational methods have long been used in modeling the formation and evolution of galaxies. However, given that these processes operate over a vast range of scales, both spatially and temporally, modeling galaxy formation in a cosmological context remains one of the greatest challenges in astrophysics today . One possible approach to predicting galaxy properties at high redshift is to constrain the relationship between galaxy and dark matter halo properties at lower redshift, and then assume that these relationships hold out to higher redshift <cit.>. Another is to implement physical recipes for the key processes that are thought to shape galaxy formation (e.g. gas accretion and cooling, star formation and stellar feedback, chemical enrichment, black hole formation and feedback) within a cosmological framework (see ,  for a recent review). These physical recipes may be implemented within numerical simulations, which attempt to explicitly follow the equations of gravity, hydrodynamics, and thermodynamics for particles or grid cells <cit.>. However, current cosmological simulations are unable to directly resolve the multiphase structure of the interstellar medium, the formation and evolution of individual stars and supermassive black hole (SMBH), or their interaction with their surroundings. Thus, “sub-grid” recipes must be used to represent these processes on scales below those that can be explicitly resolved in the simulations. In most current numerical simulations, these sub-grid recipes are generally phenomenological and contain free parameters that must be calibrated to match global galaxy observations (see discussion in ). Due to the high computational expense, rather limited dynamic range can be achieved with these techniques, and the ability to explore different sub-grid recipes and parameter values is also limited.

An alternative is semi-analytic simulations, which apply simplified recipes for these same physical processes within dark-matter halo “merger trees” either extracted from dissipationless N-body simulations or created using semi-analytic Monte Carlo techniques. Despite the absence of a spatially defined “grid” in semi-analytic models (SAMs), the essence of using phenomenological or physically driven recipes to model galaxies is very similar to the motivation of sub-grid recipes in conventional numerical simulations. One way to view SAMs is as, essentially, a book-keeping scheme that tracks the movement of mass between different reservoirs, as gas accretes from the diffuse intergalactic medium (IGM) into galactic halos, flows from halos into the interstellar medium (ISM) of galaxies, is converted from ISM into stars, and potentially is ejected from the ISM back into the hot diffuse galactic halo or back out into the IGM. In addition to gas and stars, heavy elements can also be tracked as they are produced by stars and circulate through the cosmic baryon cycle. Like the numerical simulations, the physical recipes in SAMs are phenomenological, and contain free parameters that must be calibrated to global galaxy observations. 

This method has been widely used to explore a very broad range of galaxy properties <cit.> and it has shown that the predictions for the statistical properties of galaxies obtained using these methods are in very good agreement with those from numerical simulations (). Due to the much greater computational efficiency and flexibility of these methods, however, they provide several advantages over the numerical approach, namely by allowing researchers to explore a broader range of parameter space and different physical recipes. Moreover, SAMs are also able to span a larger dynamical range in dark matter halo mass, from the smallest halos that are believed to be able to form stars to the most massive and rare objects. 

Recently, many other groups have made theoretical predictions for high-z JWST observations using numerical (e.g. ), phenomenological(e.g. ), and semi-analytic (e.g. ) methods. These works collectively probe a broad range of aspects, including the physical and photometric properties of high-z galaxies <cit.>, population synthesis and synthetic spectra <cit.>, and the impact of broadband filter choice on estimating photometric redshifts and on recovering their physical properties from observations <cit.>. Despite the great interest in using JWST to uncover the physical nature and the assembly histories of high-z galaxies, however, only a handful of studies have systematically varied the underlying physics or model parameters in a controlled way within the same modeling framework. We find that these kinds of predictions will be quite informative for the interpretation of JWST observations.

In this work, we use the well-established Santa Cruz SAM, which has been shown to successfully reproduce key galaxy observations at lower redshifts (z≲6), to make predictions for galaxy populations at intermediate to high redshifts (z = 4–10). Importantly, we make use of the Santa Cruz SAM version that tracks multiple phases of gas in the ISM (atomic, molecular, and ionized) using a suite of recipes based on empirical considerations or predictions from detailed numerical simulations containing treatments of molecular chemistry and radiative transfer (, hereafter ; , hereafter ). Our fiducial model is based on the empirically grounded assumptions that stars form in environments that are dominated by molecular gas, and that the formation of molecular gas depends on the gas surface density, metallicity, and the local UV radiation field. However, we also explore how sensitive our predictions at ultra-high redshift are to our recipes for star formation and stellar feedback. In this paper we focus on first-order, directly observable quantities, namely one-point distribution functions for rest-frame UV luminosity and observed-frame magnitudes in JWST NIRCam filters. In a series of planned papers, we will explore physical properties of galaxies such as stellar masses, molecular gas content, metallicities, and dark matter halo masses, and will also make direct predictions for the implications for reionization. 

The key components of this work are summarized as follows: the basic elements of the Santa Cruz SAM used in this work are summarized in <ref>. We then present the rest-frame UV LFs in <ref> and counts in observed JWST filters in <ref>. We discuss our results in <ref>, and summarize and conclude in <ref>.




§ THE SEMI-ANALYTIC FRAMEWORK



The SAM used in this study is the same one outlined in <cit.>. Hence, we will not describe the model in full here and we refer the reader to the following works for full detail for the  modeling framework developed by the Santa Cruz group: <cit.>; <cit.>; <cit.>;  and . Throughout this work, we adopt cosmological parameters that are consistent with the ones reported by  in 2015: Ω_m = 0.308, Ω_Λ = 0.692, H_0 = 67.8Mpc^-1, σ_8 = 0.831, and n_s = 0.9665.

Dark matter halo merger histories, more commonly known as `merger trees', are the backbone of the semi-analytic modeling framework. These merger trees can either be extracted from dissipationless N-body simulations or constructed using semi-analytic methods based on the Extended Press-Schechter (EPS) formalism <cit.>. In this work, in order to maximize the dynamic range and computational efficiency, we adopt the EPS-based method of <cit.> with updates as described in <cit.>. We have also run our models using halo merger trees from the Bolshoi Planck N-body simulations, and find  very similar results over the dynamical range spanned by that simulation. However, in addition to having inadequate mass resolution to resolve the low-mass halos that host the faint galaxy population that will be detectable by JWST, Bolshoi Planck has only 29 snapshots stored above z=6, which is inadequate to construct accurate merger histories for z>6 galaxies. In addition, robust identification of halos at high redshift becomes tricky. Friends-of-friends methods, such as used in e.g. the BlueTides simulation <cit.> and the DRAGONS simulation suite <cit.>, have been shown to artificially link together distinct halos at high redshift <cit.>. Although we are in the process of generating a new suite of simulations, halo catalogs, and merger trees that are carefully designed to model faint galaxies in the high redshift Universe (Yung et al. in prep), we are confident that the EPS-based approach has sufficient accuracy for the somewhat qualitative exploration that we present here.

At each output redshift, we set up a grid of one hundred root halos with masses spanning the range in virial velocity V_vir≈ 20 - 500, which covers halos ranging from close to the atomic cooling limit to the rarest objects expected to be detected in high-redshift surveys. We then weight each of these root halos by the expected abundance of dark matter halos of the given mass at the respective redshift, using the fitting functions  provided in <cit.> based on results from the MultiDark suite <cit.> of N-body simulations. Although this involves extrapolating the fitting functions to lower halo masses and higher redshifts than those that are directly probed by the MultiDark Suite, we have have validated these results using an unpublished suite of very high resolution, small box simulations kindly made available to us by Eli Visbal <cit.>. The assembly history is then traced down to a minimum progenitor mass of M_res, which we refer to as the mass resolution of our simulations; here we set M_res = 10^10 or 1/100th of the root halo mass, whichever is smaller. For each root halo in the grid, one hundred Monte Carlo realizations of the merger histories are generated.

After constructing these semi-analytic merger trees, our SAM implements a suite of fairly standard recipes such as cosmological accretion and cooling, star formation and stellar-driven winds, chemical evolution, black hole feedback, and  mergers. We again refer readers to  and  for full details.




 §.§ Gas Partitioning and Star Formation


In the most recent iteration of the Santa Cruz SAM (, ), disks are subdivided into annuli and the cold gas in each annulus is partitioned into an atomic (HI), ionized (HII), and molecular () component. In  and , several different recipes for gas partitioning are investigated, including an empirical recipe in which the molecular fraction is determined by the disk mid-plane pressure, and several variants of recipes in which the molecular fraction is determined by the gas surface density, metallicity, and the intensity of the local UV radiation field. Somewhat surprisingly, these studies found that most results were not very sensitive to which gas partitioning recipe was used. However, they did find that the recipe based on the prescription of <cit.>, in which the dependence on the UV background was not taken into account, failed to reproduce sufficient numbers of low-mass galaxies and was disfavored. Overall, the metallicity and UV-background dependent recipe based on simulations by <cit.> was found to perform the best, and was adopted as the “fiducial” model. As we are concerned that some assumptions contained in the empirical pressure-based model may no longer hold at high redshift, we do not explore it here.

The usual picture is that the first stars form out of “primordial” molecular hydrogen, form Pop III stars, and pollute early halos with metals <cit.>. In the current work we do not model the formation of metal-free Population III  stars explicitly. Instead (as in previous works), we set a metallicity floor of  Z_pre-enrich to represent pre-enrichment from Pop III stars. We adopt Z_pre-enrich = 10^-3Z_⊙ <cit.>. In SPT15, we show that our results do not depend sensitively on the choice of this value.

Many earlier generations of SAMs adopt what we refer to as the “classic” Kennicutt-Schmidt () relation to model the rate at which cold gas is converted into stars. In this approach, the star formation rate (SFR) is assumed to scale as a power of the total (cold; HI+) gas density <cit.>, where in some implementations (as in ), only gas above a fixed surface density is assumed to participate in star formation. This may approximate to first order the transition from predominantly molecular to predominantly atomic gas that the above multiphase modeling attempts to capture, but there is observational evidence that this critical surface density depends on gas metallicity <cit.>.

A more recent generation of SAMs adopts a -based star formation recipe <cit.>. Observations of nearby spirals have shown that the SFR surface density is nearly linearly proportional to the molecular hydrogen surface density <cit.>. However, these observations only probe  surface densities up to about 50-80 2 . Mounting evidence from both observation and theory suggests that the slope of the SF relation may steepen to ∼ 2 at higher surface densities <cit.>. We have explored both a “single slope” SF relation (which we refer to as Big1), in which the molecular gas depletion time is effectively invariant with both galaxy properties and redshift, and a “two slope” relation (referred to as Big2), in which the molecular gas depletion time decreases (and star formation efficiency increases) with increasing  surface density. The surface density of SFR is given by the expression:

    Σ_SFR = A_SF/τ_*,0(Σ_/10M_⊙pc^-2) (1+Σ_/Σ_,crit)^N_SF,

where the critical  surface density Σ_,crit = 70 2 , A_SF is the SF relation normalization, and τ_*,0 is a tunable normalization parameter. We adopt a value of A_SF from observational determinations of the relevant Kennicutt-Schmidt relation (as given in Table <ref>), and allow τ_*,0 to vary by about a factor of 50% up or down, reflecting the observational uncertainty in the true normalization of the Kennicutt-Schmidt relation (although here we find that τ_*,0=1 produces good results).  In , we found that although at z∼ 0, the predictions for galaxy properties were not very sensitive to which SF relation was adopted, an increasingly large discrepancy was seen at higher redshifts, up to the highest redshifts z∼ 6 explored in that paper. Accordingly, in this work we explore the implications of adopting these different SF recipes at even higher redshifts. The model variants are summarized in Table <ref>. 






 §.§ Photoionization Feedback

Gas accretion and, ultimately, star formation activity in a galaxy can be reduced in the presence of a strong photoionizing background <cit.>; this effect is sometimes known as photoionization `squelching' <cit.>. Studies have shown that this effect is especially effective in low-mass halos, and it is thought to be significant in suppressing the collapse of gas into small mass halos and preventing the overproduction of dwarf galaxies in the local universe <cit.>. 

As in our previously published models, we adopt the approach proposed by <cit.> to model photoionization squelching, in which the fraction of baryons that can collapse into halos of a given mass M_halo at redshift z in the presence of a photoionizing background is computed using the following function:

    f_b(M_halo, z) = ⟨ f_b⟩{1+(2^α/3 - 1) [M_halo/M_char(z)]^-α}^-3/α,

where ⟨ f_b⟩ is the cosmic mean baryon fraction, and M_char(z) is the mass at which halos retain half of the universal baryon fraction. In the present work, we adopt the redshift dependent characteristic mass obtained by <cit.> from hydrodynamic simulations including a uniform metagalactic background (based on tabulated results kindly provided to us by T. Okamoto), and adopt α = 2 as favored by those authors. We switch squelching on at a fixed redshift  z_squelch = 8, which is consistent with the redshift at which instantaneous reionization occurred as reported by the <cit.>.

The characteristic mass M_char that we adopt ranges from a value of 1.4 × 10^7  just after reionization to a value of ∼ 9.3 × 10^9  (corresponding to a virial velocity of ∼ 25 km/s) at z=0. Note that these currently favored values of the characteristic mass <cit.> are significantly smaller than the values obtained by <cit.>, and used in previous versions of the Santa Cruz SAM.

Our treatment of squelching is rather crude and of course is not self-consistent. Some previous semi-analytic studies have attempted to model photoionization squelching in a more self-consistent manner, based on the predicted UV emissivity of the galaxy population that emerges in the SAM <cit.>, and some numerical simulations have also included photoionization squelching self-consistently <cit.>. However, these studies have found that the self-consistent modeling of squelching compared with the simpler approach adopted here does not have a significant effect on predictions for observable galaxy properties or reionization <cit.>. Indeed, as we show in fig. <ref>, the  filtering mass is so low that squelching has an almost undetectable effect on any observable galaxy properties at the redshifts that we investigate here. This is in agreement with the result found by <cit.>.



 §.§ Black Hole Growth and Feedback

Top-level halos are seeded with black holes with an initial mass of 10^4. Black hole accretion and growth may occur in two modes. Mergers and disk instabilities trigger relatively rapid, radiatively efficient accretion (see S08 and <cit.> for details). In addition, a less rapid, radiatively inefficient mode of fueling is associated with Bondi accretion from the hot halo (again, see S08 for details). Two modes of AGN feedback are also implemented, associated with each accretion mode. The radiatively efficient mode drives winds that eject material from the cold gas reservoir and drive it out of the halo. The radiatively inefficient mode is associated with relativistic jets that are assumed to heat the hot halo gas, reducing or even quenching cooling. The two models of AGN feedback have associated parameters, ϵ_ wind (see Eqn. 17 of ) and κ_AGN (equivalent to κ_radio in Eqn. 20 of ). The parameter ϵ_ wind=0.5 is set based on results from hydrodynamic simulations of binary galaxy mergers (see ). The value of the parameter controlling the “jet mode” feedback (κ_AGN) is adjusted to fit constraints on the relationship between stellar mass and halo mass from abundance matching, and the stellar mass function derived from observations (see Appendix <ref>).

In practice, we find that the implementation of AGN feedback currently included in our model has no noticeable impact on the predictions presented in this paper – our results remain unchanged when we switch off both modes of AGN feedback. However, AGN feedback does impact the calibration of the model shown in the Appendix, and so we describe it here for completeness. 




 §.§ Stellar Feedback

Star formation and the baryon fraction in galaxies is thought to be regulated by large-scale outflows driven by massive stars and supernova explosions. The details of how these winds are driven and what determines their efficiency are poorly understood, and cosmological simulations are in general unable to drive these outflows directly without adopting various “tricks” (see  for a discussion). A frequently adopted assumption, loosely motivated by the expectations of momentum-conserving or energy-conserving winds (see ) is that the mass outflow rate scales with a power of the galaxy potential well depth, represented by the internal velocity dispersion or rotation velocity, times the star formation rate:


    ṁ_out = ϵ_SN( V_0/V_c)^α_rhṁ_*,

where ṁ_* is the star formation rate, V_c is the circular velocity of the galaxy, normalized by the arbitrary constant V_0 = 200, and  ϵ_SN and α_rh are tunable free parameters.

Variants of this recipe are almost universal in semi-analytic models, and several numerical cosmological hydrodynamic simulations effectively adopt it by administering “kicks” to particles according to these or similar scaling relations <cit.>. Intriguingly, several simulations in which winds are driven by attempting to implement stellar feedback “directly” also find that the wind mass loading factors η≡ṁ_out/ṁ_* follow this type of scaling. <cit.> find that their zoom-in simulations with the GASOLINE code and a “blastwave” approach to driving winds yield α_ rh∼ 2, and <cit.> find that at high redshift η depends strongly on V_c in the FIRE simulations, where α_ rh≳ 2 for low-V_ c halos and α_ rh∼ 1 at higher V_ c. 




 §.§ Stellar Populations and Dust Attenuation


For each galaxy, we store a two-dimensional grid recording the mass of stars formed with a given age and metallicity. Unattenuated synthetic spectral energy distributions (SED) are created for each galaxy by convolving this grid with simple stellar population models (SSP). Here we use the SSP models of <cit.> with the Padova1994 <cit.> isochrones and a <cit.> initial mass function (IMF). 

We include dust attenuation using a simple approach similar to the one presented in <cit.>. We assume that the face-on extinction optical depth in the V-band is given by

    τ_V,0 = τ_dust,0 Z_cold m_cold/(r_gas)^2 ,

where the optical depth normalization of dust, τ_dust,0, is a free parameter, Z_cold is the metallicity of the cold gas, m_cold is the mass of the cold gas in the disc, and r_gas is the radius of the cold gas disc. We then calculate  the attenuation in the V-band for a given inclination using a “slab” model, which assumes that the radiation sources (stars) are embedded in a slab of dust:

    A_V = -2.5log_10[1-exp[-τ_V,0/cos(i)]/τ_V,0/cos(i)] .

Here, i is the inclination, which is chosen randomly for each galaxy. We then translated this to the attenuation in the UV-band (1600), A_UV, assuming a fixed attenuation curve <cit.>. We experimented with using the two-component <cit.>-type model (as implemented in ), in which an age-dependent attenuation curve is effectively adopted by assuming that young stars are enshrouded in higher optical-depth “birth clouds”. However, we found that this type of model predicted very large amounts of extinction in high redshift galaxies, which was incompatible with observations. 

Several previous works have found that simple models for dust attenuation, in which the dust extinction optical depth scales with gas metallicity and column density with a fixed normalization, overpredict the amount of attenuation in high redshift galaxies <cit.>. One possible interpretation is that the dust-to-metal ratio is not constant, and depends on galaxy properties that effectively cause it to change with cosmic time. Indeed, recent observations have shown that the dust-to-metal ratio differs significantly across different galaxy mass scales <cit.>, and models that attempt to self-consistently track the production and destruction of dust by various channels find that the dust-to-metal ratio may change over cosmic time <cit.>. We plan to develop more detailed and physical models of dust attenuation in high redshift galaxies in future work (Popping et al. in prep), but for the present work we adopt a very simple and ad hoc approach in which we allow the dust normalization parameter τ_dust,0 to be a function of redshift <cit.>. We tune τ_dust,0 by hand and find the values summarized in Table <ref> provide a good match to available observational constraints on the bright-end of the UV LFs at z = 4 – 10. This is then implemented in our model with a functional form τ_dust,0 = exp(az+b) with parameter values a = -0.31 and b = -2.726. A comparison of the values adopted in this work and in  can be found in fig. <ref>. Note that the cosmological parameters and SAM ingredients and parameters adopted in  are a little bit different than the ones adopted in this work, which is why the values do not match up perfectly in the redshift range of overlap.





We obtained the latest version of the published NIRCam filter response functions[<https://jwst-docs.stsci.edu/display/JTI/NIRCam+Filters>] with optical telescope element (OTE). For each galaxy, we calculate its apparent magnitude by convolving the synthesized SED with the NIRCam filters. We also include absorption due to the intervening IGM. The effective optical depth of the IGM along the line-of-sight, to a source at some redshift, at wavelength λ is calculated using the expression given in <cit.>.





 §.§ Chemical Evolution
 
The production of metals is modeled using a simple approach that is commonly adopted in semi-analytic models (see e.g. ; ; ). In a given timestep where a parcel of new stars d m_* is created, a mass of metals dM_Z = y dm_* is also formed, where y is the `effective' chemical yield, or mean mass of metals produced per mass of stars. Here we assume that the chemical yield is constant. In principle, y could be obtained from stellar evolution models, but these model yields are uncertain by a factor of ∼2, and the single-element instantaneous recycling approach to chemical evolution that we are using here is somewhat crude, so we instead treat the chemical yield as a free parameter while restricting it to an expected range. Once created, metals are assumed to be mixed instantaneously with the cold gas in the disk. We track the mean metallicity of the cold gas Z_cold, and new star parcels created out of this gas are assumed to have the same stellar metallicity Z_* as the mean metallicity of the cold gas in that timestep. Supernova feedback ejects metals from the disk, along with cold gas. These metals are either mixed with the hot gas in the halo, or ejected from the halo into the “diffuse” Intergalactic Medium (IGM), in the same proportion as the reheated cold gas(see ). The ejected metals in the `diffuse gas' reservoir are also reaccreted into the halo in the same manner as the gas.

Throughout this paper, the yield y and all metallicities are given in solar units, which we take to be Z_⊙ = 0.02. Although this formally represents the total metallicity, we note that as we track only the enrichment associated with Type II supernovae, our metallicity estimates probably correspond more closely with α-type elements. Note that because enriched gas may be ejected from the halo, and primordial gas is constantly being accreted by the halo, this approach is not equivalent to a standard “closed box” model of chemical evolution.




 §.§ Calibrating the free parameters

We calibrate our models to a standard set of z∼ 0 observables, and then leave all free parameters (except the dust normalization, as noted above) fixed. Relative to the WMAP5 cosmology used in  and , the Planck cosmology adopted here results in significantly different predictions for the abundance of dark matter halos as a function of cosmic time (see ), and as a result the free parameters of the SAM need to be re-calibrated. We show the calibration quantities, along with some other diagnostic quantities that are not used in calibration, in Appendix <ref>. The details of the calibration procedure are also presented in the Appendix. The model parameters used throughout this work are summarized in Table <ref>. 






§ REST-FRAME UV LUMINOSITY FUNCTIONS



Perhaps the most basic statistical characterization of the galaxy population is the one-point distribution function of an observable quantity, such as luminosity (commonly referred to as the “luminosity function”  (LF)). The change of the LF in different redshift bins constrains the evolution of the galaxy population over cosmic time. In this section, we present the UV LFs predicted by our SAM at z = 4 - 10. With the large dynamic range of dark matter halo masses probed by our models, our predictions cover a wide UV luminosity range between M_UV∼ -6 to -24. All binned luminosity functions presented in this work are available for download online[ <https://www.simonsfoundation.org/semi-analytic-forecasts-for-jwst/>]. Throughout this work we use a tophat filter of width of 400Å centered at 1600Å to calculate the rest-frame UV luminosity.





In fig. <ref>, we show the distribution functions for the intrinsic rest-frame UV luminosities, without accounting for the effect of dust attenuation. Our results show that the choice of star formation recipe can significantly alter the number density of bright galaxies. Recall that in the GK-Big1 model, the star formation efficiency (molecular gas depletion timescale) is effectively fixed at a constant value. Keeping in mind that the  depletion timescale in nearby spirals is 1–2 Gyr <cit.>, it is perhaps unsurprising that this becomes the limiting factor in forming stars at times when the age of the Universe is significantly less than this. In the GK-Big2 model, the “steepening” of N_SF→ 2 effectively leads to a higher star formation efficiency and shorter t_ dep, mol in higher density gas. As galaxies are much more compact and gas rich at high redshift, this effectively leads to higher SF efficiencies at high redshift. The super-linear dependence of the “classic” Kennicutt-Schmidt SF relation () goes in the same direction, but to a lesser extent, as it assumes a slightly shallower slope N_ SF=1.5. Interestingly, as already suggested in , we find that the formation of molecular gas is not a significant limiting factor for star formation even at these very high redshifts, in contrast to the suggestions of <cit.>. On the other hand, the faint end of the LF is insensitive to the choice of SF model (within the limited range of models that we have tested here). We will discuss the reasons for this in <ref>.

We show our results alongside with a compilation of UV LF constraints on the bright end presented in <cit.>, which consist of both ground- and space-based observations from <cit.>. Note that the observational UVLFs are potentially impacted by attenuation by dust in these high-z galaxies, but dust attenuation is not included in the intrinsic model UV LF predictions shown in fig. <ref>. We notice that the resultant LFs from the GK-Big1 model seem to fit observations fairly well while others overpredict the population of galaxies at z ≲ 6 in the absence of dust.  At higher redshifts z ≳ 9 the results from the KS and GK-Big2 models seem to in better agreement with observations.



In fig. <ref>, we illustrate the effects of dust attenuation and photoionization squelching by comparing the outputs from our fiducial models with these effects switched on and off. The free parameter τ_dust,0 in our dust recipe is calibrated to match available constraints from HST observations by <cit.> and <cit.>. These observations along with constraints on the faint end from the Frontier Field by <cit.> are presented here for comparison. The dust model predicts that the effect of dust attenuation should be much stronger in bright, massive galaxies given their higher fraction of cold gas and the higher metallicity therein, whilst the effect on faint galaxies is minuscule. We fitted the UV LFs predicted by our fiducial model, with attenuation by dust, with the Schechter function using the least-squares method, where the Schechter function is given by

    ϕ(M)= 0.4ln(10)ϕ^* 10^-0.4(M-M^*)(α+1) e^-10^-0.4(M-M^*)

<cit.>. The best-fit Schechter parameters are presented in Table <ref>.  




In both cases, where dust attenuation is or is not included, our results show that the demographics of luminous galaxies seem to have evolved more rapidly than their faint counterparts. For instance, the “knee” feature in the LFs is rather modest at high redshifts but this feature quickly develops and becomes distinct at z ∼ 7. Moreover, we show that the position of the knee evolves as a function of redshift due to the rapidly evolving dust content. The extinction in the UV-band due to the presence of dust is estimated based on the physical properties of individual galaxies. The UV extinction as a function of M_UV is illustrated in fig. <ref>. The two-dimensional histograms are color-coded according to the conditional number density (Mpc^-3) of galaxies in each bin, which is normalized to the sum of the number density in its corresponding (vertical) rest-frame UV magnitude bin. The 50th, 16th, and 84th percentiles are marked in each panel to illustrate the statistical distribution. We have verified that the scatter in A_UV is dominated by the scatter in physical properties that are used to calculate τ_V,0 (see Eqn. <ref>), whilst the randomly assigned inclination i is sub-dominant. Our results show that even with our simple approach to modeling dust, the scatter can be quite large.

Fig. <ref> shows the β_UV-M_UV relation, where the rest-frame UV luminosity presented here includes dust attenuation, and the photometric rest-frame UV spectral slope β_UV is calculated using the following expression 

    β_UV   = log(f_λ,FUV / f_λ,NUV)/log(λ_FUV / λ_NUV)
       = -0.4 (m_FUV-m_NUV)/log(λ_FUV / λ_NUV) - 2

from <cit.>, where λ_FUV and λ_NUV are the central wavelengths of the far- and near-UV bands from the Galaxy Evolution Explorer (GALEX) survey, where we adopted λ_FUV≃ 1530  and λ_NUV≃ 2300  for our calculations. Our results are compared to the range of β_UV spanned by typical local starburst galaxies and the median values of β_UV from a compilation of observations, including GOODS-S DEEP, GOODS-S WIDE, HUDF09, and WFC3 Early Release Science <cit.>. We also compare our results to measurements reported by <cit.> at z = 7 – 9, <cit.> at z = 4 – 8, <cit.> at z = 5, and <cit.> at z = 10. The predicted β_UV-M_UV relation is is surprisingly good agreement with the observations.  

In fig. <ref>, we zoom into the bright end of the UV LFs and explore the sensitivity of the bright end behavior to the SF timescale τ_*,0, where larger τ_*,0 means less efficient star formation, and vice versa. This parameter also effectively multiplies the gas consumption timescale, and therefore larger τ_*,0 would result in higher gas fractions (again, see ). Here we increase or decrease τ_*,0 from our fiducial value by a factor of 2 and find that this results in galaxy populations with M_UV≳ -18 mildly deviating from our fiducial model, with the strength of the deviation seeming to scale with luminosity. Joint constraints on the gas and dust content and SFR in high redshift galaxies from ALMA and JWST will be extremely valuable for breaking the degeneracy between dust and star formation efficiency. With regard to the bright end, we also note that although AGN feedback is responsible for shaping the bright end of the galaxy luminosity function in our models at z≲ 2, we have checked that switching the AGN feedback on and off has no noticeable effect on any of the results presented here.





 §.§ The faint galaxy populations







Fig. <ref> shows a zoomed-in figure that focuses on the faint end of the UV LFs. Note that the effect of dust is predicted to be negligible in this regime due to the generally low metallicity in these galaxies. The turnover in the LF near M_UV∼-9 is not due to resolution, but corresponds to the sharp cutoff in the atomic cooling function at 10^4 K (which corresponds to V_ vir≃ 17 ).  The slope of the UV LFs remains fairly steep down to this limit. In addition to the estimates of the observed LF from the lensed Frontier Fields, we show constraints on the LF at z = 4 - 7 from multiple studies. <cit.> estimate the evolution of the faint end of UV LFs by measuring the star formation histories of 37 Local Group galaxies out to z∼5 using “fossil evidence” from resolved stellar populations. Interestingly, their estimates are a factor of ∼10 higher than our predictions. The "fossil" method provides a very interesting complementary approach for estimating the faint end of the high redshift luminosity function, but currently it involves some rather uncertain corrections (D. Weisz, private communication). The largest uncertainties likely arise from two factors. The first is the volumetric correction, which attempts to correct for the fact that star formation histories were not available for all known Local Group galaxies, and for Local Group galaxies that are not currently detected. The second issue is that deep resolved color magnitude diagrams (that reach the oldest main sequence turnoff) were not at the time available for many Local Group galaxies. This can lead to a bias that yields increased SF estimates at early times (D. Weisz, priv. comm.). Establishing better links between high-redshift observations and those derived from Local Universe galaxies with resolved stellar population studies is an exciting ongoing area of research, which JWST will also help to advance.

At z∼6 – 8, we show additional studies based on the Frontier Fields. Here we show both the published data from <cit.> and the unpublished, Eddington bias-corrected number densities for comparison (R. Livermore, priv. comm.). <cit.> provide constraints on z∼6 galaxies by reanalyzing the Frontier Field observations with a more comprehensive treatment for the magnification and systematic uncertainties. <cit.> estimated constraints for LFs using a combined analysis of three lensed fields with their associated parallel fields.

In addition to the fiducial model parameters, we explored the sensitivity of the faint-end slope to the efficiency of stellar driven winds by comparing four additional cases of SN feedback slope α_rh = 2.0, 2.4, 3.2, and 3.6, where larger values of α_rh imply a steeper dependence of the mass loading factor η on galaxy circular velocity, so more gas is ejected from low-mass galaxies. Although these values of α_rh are not consistent with observations at z ∼ 0, we present these results as an attempt to quantify the effects of strong stellar feedback on the low-mass galaxy populations and to illustrate how our predictions would change if the effective scalings vary with cosmic time. The results show that the faint end is indeed sensitive to the velocity scaling of the mass loading factor. Higher α_rh leads to stronger suppression of star formation in low-mass halos, and hence to a flattening of the luminosity function. The adopted value of α_rh also shifts the luminosity where the LF “turns over” from ∼ -8 (for the “strongest”) feedback to ∼ -11 (for the weakest), because α_rh also affects the slope of the relationship between luminosity (or stellar mass) and halo mass (so the luminosity of galaxies forming in halos with mass corresponding to the atomic cooling limit is shifted up or down). Note that varying the parameter ϵ_SN simply shifts the normalization of the whole luminosity function below the knee up or down, as shown in <cit.>.

Fig. <ref> shows the evolution of the cumulative number density of galaxies above some rest-frame luminosity threshold M_UV,lim. First we focus on the evolution of the bright end of the UV LFs, and present the results calculated with M_UV,lim = -19. We show the bracketing cases from our fiducial and GK-Big1 models, which provide the most and least optimistic scenario regarding forming luminous galaxies at high redshift. We compare our results to a compilation of studies from observations, empirical extrapolations, and numerical hydrodynamic simulations <cit.>. This comparison shows that results from other studies are more or less bracketed by our two star formation scenarios.

Additionally, we probe the faint end of the LFs by showing the cumulative number density for galaxies with M_UV < -15.0. Aside from our fiducial model, we include two bracketing cases where we let α_rh = 2.0 and 3.6, which correspond to weaker and stronger stellar feedback efficiency, respectively. This kind of test will place constraints on the efficiency of stellar winds and whether there are additional dependencies in the wind mass loading scalings beyond the simple assumptions adopted here. The comparison shows that the prediction from our fiducial model is similar to other studies. However, some studies predict fewer galaxies than we predicted at lower redshifts due to a shallower faint-end slope of the UV LFs, which in our model framework would require stronger stellar feedback. 

In fig. <ref>, we show the cosmic star formation rate (CSFR) for our fiducial and GK-Big1 models, with star formation rate integrated down to M_UV = -17 ± 1.0. Our result is compared to <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>. The CSFR predicted by the KS model is very similar to the ones from our fiducial model, as the UV LF results hinted, and we therefore omitted the KS model from this figure.

Overall, it is intriguing that our models, which were previously calibrated and tested only at much lower redshifts, agree so well with the existing observations all the way out to z∼ 10.  Our models predict that the number density of galaxies evolves quite rapidly between z = 4 - 10, and that the bright end evolves more rapidly than the faint end. Our models predict that the slope of the UV LFs will remain fairly steep below the current detection limit until M_UV∼-9. 
This slope is not sensitive to the choice of star formation model but is very sensitive to the scaling of the mass loading factor of stellar-driven winds with galaxy circular velocity, as is the luminosity where the LF turns over at the faint end. The bright end is sensitive to the star formation efficiency or gas depletion timescale and its dependence on gas surface density. However, the effects of dust and the SF efficiency on the bright end of the LF are degenerate, so independent probes of these quantities are needed. 





 §.§ Comparison with other models


In addition to comparing our results with observational constraints, we also compare our predictions with a collection of theoretical studies, including empirical models, semi-analytic models, and cosmological hydrodynamic simulations. Note that studies that are fully or partially numerical are subject to the inevitable tension between simulated volume and spatial and time resolution. For example, simulations with high resolution are only feasible to run over small volumes and rare, massive objects are not well sampled. On the other hand, simulations with coarser resolution and larger volumes do well at capturing formation of large-scale structure. However, reliable predictions for small objects that fall near or below the resolution limit are not possible. For this reason, the galaxy mass and luminosity range covered varies among these studies. And note that among all models compared here, our model has the widest dynamic range and is capable of carrying the simulation to z = 0 for comparison with local observations.

Overall, our results agree extremely well with other studies in regimes where observational constraints are available. On the other hand, our predictions for faint galaxy populations, for which observational constraints are very limited due to the lack of direct detections, are in close proximity with a number of other studies, so long as the simulations have adequate resolution to resolve these objects. The faint-end slopes predicted by our fiducial model tend toward the shallower end, especially at higher redshifts, among the range collectively predicted by these models.

Both empirical models included in this comparison employ some empirical relations between halo mass and star formation efficiency that are calibrated using halo abundance matching at redshifts where observational constraints are available. These relations are then extrapolated to make predictions at higher redshifts. The  model, which is based on models introduced by <cit.>, uses a redshift independent star formation efficiency that depends on the halo mass and assembly time, calibrated at z ∼ 5. The  study explores the outcome of various star formation efficiency models, including one that assumes a power-law extrapolation below a cutoff halo mass limit M_h = 2 × 10^10  and a best-fit redshift independent model. Note that their models are calibrated to bright galaxies at z = 6 – 8,  and results at z = 9 and 10 are predictions that are yet to be published. 

On the other hand, the details for the numerical hydrodynamic simulations included in this comparison are summarized in the following. BlueTides is a large-volume cosmological hydrodynamic simulation that focuses on the high-redshift (z ≳ 8) universe, with a box that is 400 Mpc h^-1 on a side, resolving galaxies with M_* ≳ 10^8  toward the end of their simulation <cit.>. Vulcan is a high-resolution (dark matter particle ∼ 10^5 ), relatively small volume (25 comoving Mpc on a side) simulation that aims to quantify the contribution of faint galaxies to cosmic reionization <cit.>. The Cosmic Dawn (CoDa) simulation is a large-scale hydrodynamic simulation coupled with radiative transfer modeling, resolving galaxies down to ∼ 10^8  in box that is ∼ 100 Mpc on a side <cit.>. Cosmic Reionization On Computers (CROC) is a cosmological simulation based on the adaptive refinement tree (ART) method <cit.>. The CROC results compared here are from their fiducial simulation box with 20 Mpc h^-1 on a side and 512^3 initial grid cells, with spacial resolution of 100 pc <cit.>. The DRAGONS project <cit.> consists of a SAM meraxes that is built on top of the Tiamat suite of n-body simulations. We also include results from a high-resolution hydrodynamic simulation with embedded self-consistent radiative transfer that is similar to the one described in <cit.>. The three runs included here have box sizes of 6.0, 7.5, and 10.0 Mpc h^-1 on a side, with 2×256^3, 2×320^3, and 2×512^3 particles; the UV background is discretized into 32, 40, and 64 voxels. 

In fig. <ref>, we compare the bright end of our UV LFs to the predictions of , , DRAGONS, CROC, CoDa, Vulcan, and BlueTides. In fig. <ref>, we compare our predictions for the faint-end behavior of the UV LFs  to , , DRAGONS, CROC, and . The  model assumed a M_UV = -12 cutoff for the atomic cooling limit in their work and the  model considered a limiting magnitude of M_UV = -13 in their CSFR calculation. For this comparison, we plot the results from these analytic studies extrapolated down to fainter magnitudes to illustrate the faint-end slopes predicted by these models and how they stack up with other numerical results. On the other hand, numerical simulations provide predictions with physical processes traced self-consistently as far down as resolution permits, although they are still limited by the assumptions inherent in their sub-resolution recipes for star formation and stellar feedback. Moreover, we note that the very faint populations predicted by these numerical studies might be somewhat resolution dependent, rather than physical. For instance, some subtle differences can be spotted between the UV LFs based on the coarser Tiamat (at z = 6 – 10) and the finer Tiny Tiamat (at z = 6, 8, and 10) N-body simulations from the DRAGONS simulation suite in our comparison <cit.>. Similar behavior can be seen among the three different runs in the  simulations. However, heating from a photoionizing background is tracked self-consistently with on-the-fly radiative transfer in this model. The Renaissance Simulations <cit.>, SPHYNX <cit.>, and FIRE <cit.> are theoretical studies that address similar issues but are not included in our comparison here.








§ PREDICTED APPARENT MAGNITUDE FUNCTIONS WITH NIRCAM FILTERS







In this section, we provide the distribution functions of apparent magnitude (AMF) for high-redshift galaxy populations, where the latest published JWST NIRCam filter response functions are used to compute the apparent magnitude for the galaxies predicted by our model. We also provide estimates of the number of detected objects for several example survey configurations.


https://jwst-docs.stsci.edu/display/JTI/NIRCam+Filters

In fig. <ref>, we showcase the evolution of the AMFs over redshift as seen in the eight NIRCam broadband filters: F070W, F090W, F115W, F150W, F200W, F277W, F356W, and F444W. The central wavelength and the detection limit of these filters are summarized in Table <ref>. Detection limits are 10-σ point source limiting AB magnitudes from the JWST Exposure Time Calculator[<https://jwst.etc.stsci.edu/>, v1.2.2] for an exposure time of 10^4 second, measured in 0.04” diameter apertures for the short-wavelength camera, and 0.08” diameter apertures for the long-wavelength camera, assuming a low background level. An approximate detection limit m_lim∼ 29 is marked in the figure with a black dashed line. Objects to the left of the line are too faint to be detected at this depth. We also marked where one object is expected per NIRCam field of view (2×2.2×2.2 arcmin^2) with a horizontal dashed line in each panel. Objects with counts below this line are too rare to be found on average in a single pointing of NIRCam. The seven panels together demonstrate how the high-redshift galaxy populations quickly drop out from the filters with shorter wavelengths due to absorption by the intervening IGM. Our results show that the F200W filter, our choice of nominal filter, is best suited for detecting objects across all redshifts of interest, and might even be able to pick up bright objects beyond z ∼ 10 with an extended exposure time. The Schechter function fitting parameters for the F200W filter are presented in Table <ref>. Redundant detections from multiple filters also improve the accuracy of estimates of dust attenuation, which is essential to properly uncover the underlying intrinsic UV luminosity. Galaxies at high redshifts rapidly drop out of the shorter wavelength filters F070W, F090W, and F115W for z < 10, and F150W at z ∼ 12, whilst F200W, F277W, F356W, and F444W remain advantageous for possible detections up to z ∼ 15 given sufficient exposure time. Although MIRI is designed to detect longer wavelength radiation, it will not be able to detect such high redshift galaxies due to its low sensitivity. 



In fig. <ref>, we further investigate the population of galaxies that are expected to be found by JWST in different surveys. We show the same AMFs from the F200W filter and provide the detection limits from a number of simulated surveys. The detection limits of past HST surveys and upcoming JWST surveys are summarized in Table <ref>. First, we show the detection limit of the HST CANDELS-Wide survey for reference. Then, from bright to faint, we show the detection limits estimated for a JWST investment comparable to CANDELS-W and HUDF, and the very optimistic case where we gain a factor of 10× magnification from gravitational lensing from a massive galaxy cluster. We indicated where ∼10 objects are expected to be found in survey areas similar to that of the HUDF field (2.4×2.4arcmin^2) and CANDELS-W (∼100 arcmin^2), which are both estimated assuming dz = 1 slices centered at the given redshift. We also show the cases where we deviate from the fiducial value and let α_rh = 2.8 ± 0.4 and ± 0.8, where stronger feedback gives shallower faint-end slopes and vice versa. This comparison shows that JWST will be able to reliably constraint the faint-end slope up to z∼8, or even z∼10 in lensed fields. This will provide important constraints on stellar feedback physics.





In fig.<ref>, we show the number of objects per Mpc^3 above some detection limits, m_lim, expected to be found using the NIRCam F200W filter. We consider two cases where we let m_lim,wide = 29 and m_lim,deep = 31.5, representing the expected detection limit for JWST wide-field surveys (e.g. CANDELS-Wide) and deep-field surveys (e.g. HUDF). In the upper panel, we illustrate the difference among the three SF models considered in this work, and in the lower panel we show the cases where we let α_rh = 2.0 and α_rh = 3.6. Here, a more direct comparison with observations can be made, without the need for applying K-corrections to the observations to get from observed to rest-frame magnitudes. One can see that, similarly, observations of the bright population will constrain the star formation efficiency, while observations of the faint end will constrain stellar driven winds.  








§ DISCUSSION



The formation and evolution of galaxies are governed by a complex network of intertwined physical processes that operate over many physical scales. Although many details remain to be worked out, the community seems to have reached a broad consensus regarding the main physical processes that shape galaxy properties at z ≳ 6; these include cosmological accretion, strong stellar-driven winds that are more efficient at low masses, and (more controversially), black hole feedback that preferentially suppresses star formation at high masses. However, these processes can sometimes have degenerate effects on galaxy properties. Although differentiating the effects of these processes is extremely challenging, being able to do so is crucial for interpreting galaxy observations at all redshifts. Observational constraints at z ≳ 6 are currently quite limited, leaving significant uncertainties in the current theories of galaxy formation and their predictions at very high redshifts. One of the main highlights of this work is that we explored how variations in several uncertain physical processes in theoretical models will impact the global physical properties of galaxies and what may be seen by future observations. We also discuss how these processes may be disentangled with the significant insights JWST observations have to offer.

High-z star-forming galaxies are most luminous in the observed frame infrared, making them observable with NIRCam and MIRI onboard JWST, and it has been shown that the choice of broadband filters can have an impact on the photometric redshift estimation and the follow-up interpretation <cit.>. In fig. <ref>, we demonstrated how objects drop out from shorter wavelength filters starting at z ∼ 5, making the F200W filter the nominal choice for blind surveys for these deep-field objects across z = 4 – 10 with investments comparable to CANDELS-Wide and HUDF. For galaxies beyond z ∼ 10, observations in the longer wavelength filters, namely F356W and F444W, will be required with exposure times ≳ 10^5 seconds. 

With its unprecedented sensitivity, JWST will be able to probe large populations of galaxies in deep space and provide constraints for the faint-end slope of the UV LFs at z ≲ 8 and for the evolution of the bright-end up to z ∼ 10. These observations will put the empirical relations derived from low-redshift observations to the test in extreme conditions. Moreover, the redundancy in multiple NIRCam broadband filters at lower redshifts will enable multiwavelength measurements, which is crucial to recovering the physical properties of high-z galaxies and breaking degeneracies in the underlying physical processes.

In this work, we utilized a computationally efficient SAM to make quantitative estimates of observable properties for galaxy populations expected to be detectable with JWST. The physically motivated empirical recipes adopted in the model are well-tested at lower redshifts (z ≲ 6, see ). We present model predictions at redshifts where these models have never been tested before, and where observational constraints are relatively limited. By exploiting the efficiency of our model, we also provide forecasts for various model variations, including exploring multiple SF recipes and parameterizations of outflow rates for stellar-driven winds.

Physical processes that shape the formation of galaxies are degenerate, and yet each of them evolves slightly differently over time and can simultaneously affect multiple physical properties. Traditionally, there has been tension for galaxy formation models to simultaneously match the observed gas fraction, stellar metallicity, and stellar fraction. In order to calibrate our model, we carefully balance the model parameters for multiple physical processes to match observations at z ∼ 0. Fig. <ref> summarizes how the outputs of our calibrated model compare with observational constraints. See appendix for details of the calibration process and tests.

Massive star forming galaxies in the lower redshift Universe are obscured by dust in the line of sight, which makes it extremely difficult to determine the true underlying stellar content. In our model, the attenuation effect due to dust in the galactic disc is estimated based on galaxy physical properties, incorporated with a dust-to-metal ratio guided by observations. In agreement with previous studies, we find that the optical depth of dust, or effectively the dust-to-metal ratio, is required to evolve with redshift. However, we also note that our very simple dust recipe does not accurately represent the complex geometries of dust relative to stars that may be common in high redshift galaxies <cit.>.

We find that star formation remains fairly efficient in low-mass halos, in conflict with the <cit.>  model, which argues that star formation is heavily suppressed at z > 2 in dark matter halos with masses < 10^11  due to the low metallicity of gas in these halos. However, recent lensed and very deep fields observations provide strong support for LF remaining fairly steep down to low halo.  masses. We are aware of the fact that the underlying assumptions of these metallicity-based multiphase-gas partitioning recipes tend to break down in extremely low-metallicity environments (e.g. Z_cold < 0.05 Z_⊙ <cit.>, which might cause our model to over-predict the stellar content in our least massive halos and at higher redshifts.

Cosmic reionization in our models is treated with a rather ad-hoc prescription, in which the whole universe is reionized uniformly and instantaneously. Our models currently do not attempt to model the formation of metal-free Population III stars explicitly. Instead, we set a metallicity floor of  Z_pre-enrich to represent pre-enrichment in the initial hot gas in halos and the gas accreted onto halos due to cosmological infall.




 §.§ Physical processes that shape the bright end of UV LFs


The strong evolution in the bright end of the UV LFs indicates a rapid growth of massive galaxies in the early universe. However, the physical processes that drive or regulate this evolution are still unclear. Due to the limited direct observations at z ≳ 6, it is extremely difficult to disentangle the effects of multiple processes.

Each observed property of galaxies is resulting from a unique combination of underlying physical processes. Although the underlying fundamental physics should remain the same at all times, since the physical processes that drive the evolution of these objects depend on a number of effectively redshift dependent conditions, the global properties of galaxies may also be effectively redshift dependent. The evolution of these resulting global physical properties change subtly depending on the strength of the impact of each of the physical processes. Therefore, being able to disentangle the impact of these driving physical processes, as well as the evolution of the resultant observable properties, is most important for understanding galaxy formation and for making accurate predictions for the galaxies that are yet to be directly observed. 

In this work, we have shown that the choice of star formation recipe and the gas depletion time both have a very strong impact on the abundance of bright galaxies. We have also provided model outputs where we systematically vary some of the physical recipes and associated parameters. For instance, changing the gas depletion timescale by a factor of two in either direction seems to mainly impact galaxies of M_UV≲  -19.  Under this variation, the bright end of the UV LFs deviates mildly from the fiducial results but still remains mostly within the uncertainties in the observations. On the other hand, altering the SF recipes would significantly change the predicted number density for luminous galaxy populations, where the distinction occurs at M_UV≲ -14 at z = 10, and this limit gradually evolves to M_UV≲ -20 at z = 4. The abundance of faint galaxies is largely insensitive to either of these processes due to self-regulation of SF in low-mass halos, which will be discussed in detail in the next subsection.

To focus on the impact of these processes on the underlying stellar populations, in fig. <ref> we compare the intrinsic UV LFs in the absence of dust. The free parameters in the empirical SF recipes that are calibrated to z ∼ 0 observations and the SF relation do not evolve over time, with the fiducial GK-Big2 and GK-Big1 representing the more and less optimistic SF scenarios, respectively. Both our fiducial model and the classic  model are able to match high-z observations quite well without further parameter retuning. Given that the general consensus is that high-z galaxies are much less dusty than their low-z counterparts, these empirical SF recipes might hold to a certain extent even in extreme environments. Interestingly, GK-Big1 is able to match observations quite well at low redshifts (z ≲ 6) when dust attenuation is omitted, but it underpredicts the number of bright galaxies at higher redshifts. This implies that if the dust attenuation at these redshifts is significantly lower than what we have assumed, a mildly evolving Σ_SFR–Σ_ relation may be able to produce sufficient numbers of bright galaxies at high redshifts.

It is well known that at z∼ 2 – 4, the most rapidly star-forming galaxies are heavily obscured by dust. The dust content of galaxies at higher redshift is uncertain. Since dust is built up over generations of star formation, the general expectation is that high redshift galaxies should be less dusty than lower redshift ones <cit.>. Moreover, there is mounting evidence that the dust geometry may be different in high redshift galaxies, giving rise to differences in attenuation for a given dust mass <cit.>. Most theoretical models that include the effects of dust assume a relationship between dust optical depth and the metallicity and gas surface density, similar to the one we have adopted here. Our simple empirical model for dust is calibrated guided by available observations at z = 4 – 10, which requires a redshift dependent dust-to-metal ratio in order to be able to simultaneously match constraints at all redshifts within the framework of our simple model. Our results show that dust attenuation seems to be less effective at M_UV≳ -19, however, with a much stronger redshift evolution comparing to bracketing cases for gas depletion times we showed.  Future multiwavelength observations  with multiple NIRCam filters may help constrain the dust content (e.g. using the UV continuum slope β as in , also see ).

Using similar SAMs and a similar dust model,  found that their models with a Calzetti attenuation curve underproduced UV-luminous galaxies at high redshift (z∼ 3–5). Our updated models perform better in this regard, as a result of the updated Planck cosmology and associated recalibration. As shown in fig. <ref>, the amount of dust required to match observations between = 4 – 6 has been reduced by roughly a factor of two. However, the UV LFs remain quite sensitive to τ_dust,0 and, as pointed out in previous works, using a fixed dust parameter that is normalized to observations at low redshifts would systematically underproduce UV-luminous galaxies at high redshift. In future works we plan to explore models with a more self-consistent treatment of dust formation and destruction and more complex treatment of dust geometry, and the implications for Far-IR and mm and sub-mm observations (Popping et al. in prep).

Although we find that black hole feedback plays a negligible role in shaping galaxy properties at these redshifts, this may be due to the specific manner in which black hole seeding, growth, and feedback are implemented in these models.  AGN feedback is implemented in our models via two different modes: “jet mode” (also called “radio mode”) and “radiative mode” (also called “quasar mode” or “bright mode”). “Radio mode”  feedback (, ) is implemented as a heating term which (in the Santa Cruz SAMs) scales as a power-law function of the black hole mass, based on observations of radio jets in nearby galaxy groups and clusters (see ). The jets are assumed to be able to couple efficiently with the hot gas only when the cooling time is longer than the dynamical time, which tends to be the case in massive halos at late times. As a result, jet mode feedback becomes effective only in massive halos at redshifts below about z∼ 1.

The Santa Cruz SAMs also include “radiative mode” feedback, in which cold gas can be removed from galaxies via radiation pressure driven winds associated with radiatively efficient accretion onto a black hole. Although this mode can act at high redshift, it tends to have little effect on galaxy properties, because it is assumed that the winds can only act “ejectively” on the cold interstellar gas, and new gas tends to cool rapidly and replenish the cold gas reservoir, especially at high redshift. However, the treatment of radiative mode feedback in the existing models is based on an older suite of hydrodynamic simulations of binary mergers between idealized galaxies with no initial hot gas halo. More recent cosmological zoom-in simulations including thermal and kinetic feedback from radiatively efficient black hole accretion find that there is also a strong preventative feedback effect, as these winds can significantly reduce the density of gas near the centers of halos and thereby suppress cooling for much longer timescales <cit.>. In this picture, we might expect quenching via radiative mode feedback to be more effective at higher redshifts. 




 §.§ Physical processes that shape low-luminosity galaxies


In this work, we updated the Santa Cruz SAM treatment of photoionization squelching by implementing the fitting function presented in <cit.> to model heating from an photoionizing background. Our results in fig. <ref> show that the effect is negligible at all redshifts. Some studies suggest that photoionization squelching plays a  significant role in shaping the faint end of the galaxy LF, suppressing the faint end of the UV LFs in the range M_UV∼ -12 – -10 <cit.>. Some others have shown that dense clumps are extremely hard to penetrate and become ionized due to a high recombination rate. Simulations have shown that cooling in gas clumps should be fairly efficient when densities are comparable to the virial density of a halo <cit.>. <cit.> used radiative hydrodynamic simulations with radiative transfer to show that a photoionizing background is devastating only for low-mass systems of V_c ≲ 20 , which are halos close to the atomic cooling limit. Similar results are also presented by <cit.>, where only extremely low density gas clumps are completely evaporated and the overall effect of a photoionizing background on galaxy formation is much weaker than previously thought.

In our models, stellar-driven winds play a dominant role in shaping the slope of the faint-end of the UVLF as well as the location of the turnover. This type of simple empirical model for stellar-wind feedback has been widely adopted in both SAMs and some cosmological numerical hydrodynamic simulations. The recipe has been shown to be quite successful in reproducing observations at lower redshifts for galaxies across a wide range of masses. As discussed in , SF in low-mass halos is self-regulated by stellar feedback. For instance, when SF becomes more efficient, more gas is ejected by energetic stellar winds, and hence reducing the supply of cold gas and yielding less efficient SF, and vice versa. This effect has also been examined and demonstrated in a number of studies <cit.>. Here we show that this remains the case in low-mass halos up to very high redshift, even in models with metallicity dependent, -based star formation. 

The free parameters in our physically motivated stellar feedback recipe have only been tuned to match observed SMF and stellar-to-halo mass ratio at z ∼ 0. In addition to the calibrated fiducial values, we also show example cases where we systematically vary the SN feedback slope α_rh slightly by ± 0.4 and ± 0.8 while keeping the rest of the model unchanged. Our results clearly demonstrate that this parameter indeed plays a dominant role even at extreme redshifts. However, we also show that the more massive galaxies are insensitive to variations in this parameter, because stellar-driven winds cannot efficiently escape the deep potential wells of the halos that host these objects.

Our model is one of the very few simulations that is capable of resolving objects as tiny as V_c ∼ 20 in a cosmological context, close to the atomic cooling limit. We find that the turnover in the UVLF due to the atomic cooling limit occurs at around M_UV∼ -8 at z∼ 4, moving slightly brighter to about M_UV∼ -9 at z∼ 10. We find that the magnitude where the turnover occurs shifts brighter or fainter by up to ∼ 1 magnitude under variations in the stellar feedback parameter α_rh. This is because α_rh changes the slope of the median relationship between halo mass and stellar mass or SFR, shifting the rest-UV magnitude corresponding to a halo circular velocity of V_c ∼ 20 .




 §.§ Probing ultra high-redshift galaxies beyond JWST
 

The galaxy populations predicted by our model span a wide range of luminosities and redshifts. Even though JWST will be able to detect many objects during the epoch of reionization, our models predict that there will still be a significant population of objects too faint to be detected even by JWST.  Although these galaxies are unlikely to play a significant role in cosmic reionization, they are thought to have hosted Pop III stars and polluted the ISM with the first heavy elements.  Therefore, constraints on these objects are a very important missing piece in the current formation theory for stars and galaxies. Probing these objects directly requires instruments with sensitivity many times higher than what JWST has achieved. However, the bulk effects of these objects can be studied via metal absorption (e.g. ) or intensity mapping <cit.>. Many ongoing and planned intensity mapping pathfinders (e.g. CHIME, HERA, HIRAX, Tianlai, BINGO, LOFAR, MeerKat, CONCERTO, STARFIRE) are paving the way to future high-z large-scale multiline intensity mapping surveys. 

On the other hand, planned wide-field surveys, such as those that will be carried out with Euclid and the Wide-Field Infrared Survey Telescope (WFIRST), will probe unprecedented areas, providing  constraints on the massive, bright galaxy populations <cit.>.

We plan to exploit our model framework to forecast and provide an interpretive framework for these and other observations in future projects. We plan to further explore the progression of cosmic reionization arising from our predicted galaxy populations, and compare that to the latest observational constraints from the Lyα forest and the Thomson scattering optical depth for the cosmic microwave background (Yung et al., in prep).




§ SUMMARY AND CONCLUSIONS



In this work, we presented predictions for galaxy populations that are expected to be detected in upcoming JWST observations, and showed how they can  potentially constrain the physical processes that govern the formation and evolution of these objects. Our galaxies are modeled using the well-established Santa Cruz semi-analytic model with the recently updated multiphase gas partitioning and -based SF recipes. We also used semi-analytic dark matter halo merger trees that are constructed based on the EPS formalism to achieve the very wide dynamic range and computational efficiency required for forecasting observable properties of galaxy populations over a wide range of masses and redshifts. Moreover, we adopted the updated Planck cosmology and the model parameters were recalibrated to match the latest observational constraints near z∼ 0.

By exploiting the high efficiency of our model, we were also able to systematically vary the SF recipes, as well as the sub-grid physical parameters for gas depletion timescale and stellar feedback relation slope, in a controlled manner. We use these results to explore the physical processes that have degenerate effects on galaxies, which create tension in matching galaxy properties and other cosmological observables. We also  discuss whether these processes can be disentangled and what we expect to learn from the upcoming deep-field observations. 

Predictions for rest-frame UV luminosity functions at z = 4 – 10 are presented and are compared to existing observations and other models. We include the effects of dust attenuation using an empirical dust recipe with a redshift dependent dust-to-metal ratio. Although the free parameters in our model are only calibrated to match observations at z ∼ 0, our results matches surprisingly well with UV LF constraints at z > 4. In addition, our results agree extremely well with previous theoretical studies, particularly for more luminous galaxies. We predicted that the faint end of the UV luminosity functions will remain steep below the current detection limit until M_UV∼ -9. We showed that the gas depletion time and the choice of star formation recipe have strong influences on star formation in luminous galaxies. Conversely, the effect of feedback from AGN is found to be negligible, although this may be due to shortcomings in our modeling of black hole seeding, accretion, and/or feedback.  However, starlight from the most intrinsically UV-luminous galaxies is also heavily obscured by dust, making it extremely difficult to disentangle these degenerate effects from multiple physical processes. On the other hand, star formation in low-mass halos seems to be most strongly affected by stellar feedback, with photoionization feedback having a negligible effect on the populations that we studied. 

We estimated the apparent magnitudes of our predicted galaxy population utilizing the published JWST NIRCam broadband filters and presented them in the form of one-point distribution functions, which may be compared directly with observations. We also estimated the effects of dust attenuation and illustrated the sensitivity to stellar feedback efficiency. We show that, with a simple dust model, the effect of dust is only significant for rapidly star forming, metal rich galaxies and is more important at lower redshifts. 

We summarize our main conclusions below. 

(i) A relatively simple and computationally efficient semi-analytic model, which has been calibrated only to z ∼ 0 observations, produces predictions that agree remarkably well with observed UV luminosity functions from z∼ 4–10.

(ii) Star formation physics and gas depletion time are dominant in determining the abundance of bright, massive galaxy populations, and these physical processes are degenerate with the attenuation due to dust. Therefore, it is critical to obtain independent probes on the dust content of high redshift galaxies in order to be able to disentangle the underlying physics.

(iii) The faint-end slope of UV LFs is mainly sensitive to the scaling of the mass-loading factor for stellar driven winds with halo or galaxy properties. We find that the effect of photoionization squelching on galaxies that will be detectable with JWST is negligible. JWST observations will be able to place important constraints on stellar feedback. However, JWST will not be able to probe down to the atomic cooling limit.

(iv) In our models, the absolute magnitude at which the “turnover” in the UV LF occurs due to the atomic cooling limit is also sensitive to the adopted recipe for stellar feedback. This is because stellar feedback changes the halo mass that hosts galaxies of a given UV luminosity.



§ ACKNOWLEDGEMENTS

The authors of this paper would like to thank Aldo Rodríguez-Puebla, Viraj Pandya, Peter Behroozi, Takashi Okamoto, Daniel Weisz, Eli Visbal, and David Spergel for useful discussions. We also thank Kristian Finlator, Rachael Livermore, and Guochao Sun for providing unpublished results for comparison. We also thank Rebecca Sesny and Dylan Simon for the creation of the project website for online data release. AY and RSS  thank the Downsbrough family for their generous support, and gratefully acknowledge funding from the Simons Foundation.



mnras








§ TABULATED VALUES FOR SELECTED UV LFS AND AMFS


 

Tabulated UV LFs from our fiducial model, with and without dust attenuation, and AMFs calculated with our nominal NIRCam F200W filter are provided in the tables <ref> – <ref>. Other UV LFs shown in this work are available online.










§ RE-CALIBRATION FOR PLANCK COSMOLOGY


 

In this appendix, we briefly summarize the calibration of our model after updating the cosmological parameters to values consistent with recent constraints from Planck. The free parameters in our model are calibrated by hand such that the outputs from our fiducial model match observations at z ∼ 0. The calibration quantities are the stellar-to-halo mass ratio and stellar mass function, stellar mass-metallicity relation, cold gas fraction versus stellar mass relation for disk-dominated galaxies, and the black hole mass vs. bulge mass relation. We adopt observational constraints for these quantities from <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>. We show our model outputs compared with these observational constraints in fig. <ref>. In addition, we show two quantities used as a cross-check but not used directly in the calibration: cold gas phase mass-metallicity relation (lower-left panel in fig. <ref>) and the  mass function (fig. <ref>), along with observational constraints from <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>. For scaling relations, the blue solid line marks the median and the dashed lines mark the 16th and 84th percentile for central galaxies. Both central and satellite galaxies are included in the distribution functions (such as the SMF and  mass function).  has shown that the different SF recipes produce results that converge at z ∼ 0. The re-calibration for updated cosmology have uniform effects across the models and therefore should not qualitatively change the results among models. Hence, we refer reader to such work for the relative differences among SF models at low redshifts.  

Traditionally, there has been tension in the predicted physical properties in galaxy formation models. In one case, matching the stellar fraction closely tends to lead to an excess in gas fraction. Since the underlying star formation efficiency that has immediate effects on stellar fraction is directly related to the production of metals, which our multiphase-gas partitioning, star formation, and dust attenuation recipes are extremely sensitive to, calibrating our model to simultaneously match the observed gas fraction, stellar metallicity, and stellar fraction is very challenging. The uncertainty of about a factor of two in the observed normalization of the Kennicutt SF law and the observed gas fractions provides us with some leeway for calibrating our model. After carrying out multiple tests with different configurations, we find that leaving the SF timescale τ_*,0 close to unity yields the best results, with gas fractions slightly higher than values reported by <cit.>. If we strictly enforce the gas fraction to match observations by decreasing τ_*,0, we will also need to compensate for that by increasing AGN feedback to keep the model from overproducing massive galaxies. 

While we configure the SF timescale and AGN feedback to fit the massive end of the stellar-to-halo mass ratio and stellar mass function, the faint populations seem to be more sensitive to the SN feedback slope α_rh. Structure forms earlier in the Planck cosmology relative to the WMAP-5 cosmology used in our previous work, leading to a higher number density of low-mass dark matter halos at early times. This requires us to  increase α_rh to suppress the formation of low-mass galaxies in order to match both the stellar-to-halo mass ratio and stellar mass function, which comes with a side effect of further steepening the mass-metallicity relation. As discussed in previous works , this tension, along with other discrepancies at intermediate redshifts <cit.> hint that our rather simple recipe for stellar feedback needs to be revised. Moreover, since the abundance of low-mass galaxies is very sensitive to the choice of α_rh, deviating from the calibrated value would certainly would certainly raise tension with observational constraints at z ∼ 0. The impact of alternative α_rh values at z ∼ 0 has been examined in <cit.>, which they found that making alternate assumptions for ejected mass would lead to a ∼ 0.25 dex changes relative to the fiducial model. Some of the values, including the ones disfavored by local observations are explored in this work in effort to quantify its effect on the low-mass galaxy populations. The differences among the three star formation recipes presented in this work have been explored in . 








§ TESTING THE HALO MASS FUNCTION

 
In this work, we use the fitting functions for halo mass functions provided by <cit.> that is fitted to the Bolshoi-Planck simulations <cit.>. However, the mass resolution from Bolshoi-Planck simulation is well above V_vir∼ 20 which our grid of root halos reaches. To check the validity of these fitting functions for HMF, we used are compared to very high resolution, small box simulations that are similar to the ones presented in <cit.>. A comparison between the halo mass functions and the numerical simulations is shown in fig. <ref>.






