Bursty Star Formation Naturally Explains the Abundance of Bright Galaxies at Cosmic Dawn

In this Letter, we analyze the same set of simulations as recently studied by Sun et al. (2023a), which is a subset of the High-Redshift suite (Ma et al., 2018a, b, 2019) of the FIRE-2 cosmological zoom-in simulations (Hopkins et al., 2018). The FIRE-2 simulations use the GIZMO code with its meshless-finite mass (MFM) hydro solver (Hopkins, 2015), and include multiple channels of stellar feedback to regulate star formation. Star formation occurs in dense molecular gas ($n_{\mathrm{H}}>1000\,\mathrm{cm^{3}}$) that is self-gravitating and self-shielding. The stellar feedback mechanisms implemented include: (1) energy, momentum, mass, and metal injection from core collapse and Type Ia supernovae and winds from OB and AGB stars, (2) photoionization and photoelectric heating, and (3) radiation pressure. A redshift-dependent but homogeneous ionizing background is also included following Faucher-Giguère et al. (2009).³³3The version of the ionizing background used in these simulations reionizes the universe at $z_{\rm reion}\approx 10$, which is earlier than the mid-point of reionization of $z_{\rm reion}\approx 8$ favored by more recent observational constraints (e.g., Planck Collaboration et al., 2020; Faucher-Giguère, 2020). However, our main results focus on the bright end of the UVLF, which arises from relatively massive halos, whereas the suppression of galaxy formation due to heating by the ionizing background primarily affects low-mass halos ($M_{\rm h}\lesssim 10^{9}$ $M_{\odot}$; e.g. Gnedin, 2000; Noh & McQuinn, 2014). Moreover, an earlier reionization redshift implies that in the present simulations, galaxy formation is suppressed starting earlier in the small halos, so adopting a more up-to-date reionization model would (if anything) enhance the predicted UV luminosity density. Similar arguments apply to other IGM heating processes. The baryonic (dark matter) mass resolution of the set of simulations considered in this work is $m_{\mathrm{b}}=7\times 10^{3}\,M_{\odot}$ ($m_{\rm DM}=4\times 10^{4}\,M_{\odot}$), except for the simulations z5m11a and z5m11b, which have $m_{\mathrm{b}}\approx 1\times 10^{3}\,M_{\odot}$ ($m_{\rm DM}=5\times 10^{3}\,M_{\odot}$). The gravitational softenings are fixed in physical units to $\epsilon_{\mathrm{DM}}=42$ pc for the dark matter and $\epsilon_{\mathrm{star}}=2.1$ pc for stars. The gravitational softenings are adaptive for gas, with a minimum of $\epsilon_{\mathrm{b}}=0.42$ pc. This is, again, with the exception of z5m11a and z5m11b (see Figure 4 for a list of simulation IDs considered in this work), which have $\epsilon_{\mathrm{DM}}=21$ pc, $\epsilon_{\mathrm{star}}=1.4$ pc, and $\epsilon_{\mathrm{b}}=0.28$ pc.

Part of the High-Redshift suite of simulations was presented and analyzed in detail by Ma et al. (2018a, b) for the predicted properties of the simulated galaxy population at $5\leq z\leq 12$, including sizes, morphologies, scaling relations, and number statistics measured by the stellar mass and luminosity functions. In this follow-up analysis of Ma et al. (2018b) motivated by recent JWST observations of the abundance of galaxies at $z\gtrsim 10$, we follow closely the methodology adopted in Ma et al. (2018b) for fair comparisons, but the sample size of high-$z$, massive galaxies has been substantially increased to better determine the bright-end behavior of galaxy UVLFs at cosmic dawn. Below, we will only briefly summarize the key information about the sample of simulated galaxies pertinent to the analysis presented here. We refer interested readers to the aforementioned papers for further details about the FIRE-2 simulations and the High-Redshift suite.

For a robust analysis of UVLFs at their bright end, we build a maximum possible sample size of massive galaxies by making use of all the zoom-in simulations available at each redshift above the ending redshift $z_{\mathrm{end}}$. In each zoom-in region, we consider all the well-resolved halos⁴⁴4Halos containing at least $10^{4}$ particles in total and uncontaminated by low-resolution particles are considered “well-resolved”. that host a central galaxy, rather than the one hosting just the most massive, primary galaxy (typically near the center of the zoom-in region). Following Ma et al. (2018b), we define galaxies based on catalogs of halos identified with the Amiga Halo Finder (AHF; Knollmann & Knebe 2009). The radius $R_{\mathrm{max}}$ at which the halo rotation curve reaches maximum is used to define a galaxy by incorporating star particles within $R_{\mathrm{max}}/3$ and excluding the contamination from subhalos outside $R_{\mathrm{max}}/5$. We restrict the scope of our UVLF analysis to halos with mass $M_{\mathrm{h}}>10^{7.5}\,M_{\odot}$ in snapshots at $8\leq z\leq 12$ because most of the recent UVLF measurements at $z<8$ with JWST can be well explained by previous theoretical predictions and a sufficiently constraining sample of spectroscopically-confirmed galaxies is not available at $z>12$ (Harikane et al., 2023a). In Appendix A, we illustrate how the halo/galaxy sample is constructed with (snapshots of) the 26 individual zoom-in simulations, which build up a total sample of $\approx 25,000$ galaxy snapshots over $8\leq z\leq 12$. For all simulations, snapshots are saved at a cadence of every 10–20 Myr.

We process the simulated galaxy sample in order to arrive at their 1600 Å UV magnitudes $M_{\mathrm{UV}}$, following Sun et al. (2023a). Templates of binary, single-stellar-population (SSP) spectra from BPASS v2.1 (Eldridge et al., 2017) are interpolated and applied to star particles according to their stellar age and metallicity, assuming a Kroupa IMF (Kroupa, 2001). Including nebular (continuum) emission can in principle augment both the UV emissivity and variability (Byler et al., 2017), although we opt to ignore it here as nebular emission is not expected to strongly affect the measurement of $M_{\mathrm{UV}}$, especially when compared with effects of SFR variations. Two notable differences from Sun et al. (2023a) exist, though, for the treatment of (1) the connection between $M_{\mathrm{UV}}$ and the SFH and (2) dust attenuation, on which we elaborate below.

Using UV magnitudes derived for the sample of simulated galaxies binned into redshift bins of width $\Delta z=\pm 0.5$, we calculate the UVLF through a convolution with the halo mass function (HMF) following the “HMF-weighting” method introduced by Ma et al. (2018b). This method has been verified to provide robust estimates of the UVLF from galaxy samples drawn from zoom-in simulations, so we only summarize briefly here. First, in narrow halo mass and redshift bins, we count the number of simulated halos $N_{\mathrm{S}}$ from the sample and compute the expected number of halos $N_{\mathrm{E}}$, which scales with the HMF, $\mathrm{d}n/\mathrm{d}\log M_{\mathrm{h}}$, calculated using the hmf code (Murray et al., 2013) for the fitting function from Behroozi et al. (2013). A common weight $w=N_{\mathrm{E}}/N_{\mathrm{S}}$ is assigned to all the halos in the same bin, such that a summation of halo weights in a given mass bin yields the expected number of halos in the universe. These weights are then applied to sample galaxies binned in $M_{\mathrm{UV}}$ to obtain the UVLF, which is essentially a convolution between the HMF and $M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation including the full, $M_{\mathrm{h}}$-dependent distribution (see Section 2.4 of Ma et al., 2018b). Finally, we stress that, compared with Ma et al. (2018b) where only a subset of the High-Redshift suite was analyzed, we substantially increase the number of samples of massive halos/bright galaxies in this work (a factor 8 increase of halos with $M_{\mathrm{h}}>10^{10}\,M_{\odot}$ at $z=10$) by considering the full High-Redshift suite as in Ma et al. (2019), thereby extending the magnitude down to which the UVLF at $z>10$ can be reliably determined to $M_{\mathrm{UV}}<-20$, overlapping with the bright-end UVLF probed by JWST.

Following the methods outlined in Sections 2.2 and 2.3, we first use our samples of simulated galaxies to quantify the $M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation in different redshift regimes, assuming either “bursty” or “smoothed” SFH. A comparison of the $M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relations at $z=8$, 10, and 12 from our simulations is shown in the top row of Figure 1. Overall, galaxies become more UV-bright at higher $M_{\mathrm{h}}$ and, at a given $M_{\mathrm{h}}$, $M_{\mathrm{UV}}$ decreases modestly with increasing redshift as a result of more rapid halo growth at higher redshift. A significant scatter in $M_{\mathrm{UV}}$ around the median relation that gradually increases towards lower masses exists, which is a sign of increasing star formation burstiness at low masses, given the proportionality between $L_{\mathrm{UV}}$ and the SFR. At a fixed $M_{\mathrm{h}}$, we find a modest trend for the scatter in $M_{\mathrm{UV}}$ to decrease with decreasing redshift that continues to $z<8$ (not shown). This tentative evidence for the redshift evolution of the UV variability might be testable using comparisons of SFR indicators sensitive to different star formation timescales or high-precision measurements of the halo–galaxy connection with galaxy clustering (see Section 4). From the comparison between the “bursty” and “smoothed” cases shown by the 5–95th percentiles (especially in the top middle panel where three “smoothed” cases with varying $\tau_{\mathrm{SF}}$ are shown), it can be seen that evaluating $M_{\mathrm{UV}}$ from a smoothed SFH leads to a shallower $M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation with a reduced scatter in $M_{\mathrm{UV}}$ at higher masses, which effectively suppresses the population of UV-bright galaxies at a given $M_{\mathrm{h}}$.

In the bottom row of Figure 1, we show the UVLF at $z=8$–12 implied by the $M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation. From the comparisons against recent observational constraints and between the two SFH cases, several key results are immediately apparent. First, in the fiducial, “bursty” SFH scenario, the predicted UVLFs agree remarkably well with the observational constraints available. In particular, our $z\gtrsim 10$ predictions lie safely above the firm lower bounds set by the dust-uncorrected, spectroscopically-confirmed samples recently compiled by Harikane et al. (2023a), and they are also broadly consistent with the variety of measurements based on photometrically selected candidates (see the caption for details)⁶⁶6We have verified by bootstrapping 1000 times the simulated galaxy samples that the statistical uncertainty on the UVLF, especially at the bright end, is small enough that it does not affect the bright-end comparisons of interest to this study. In the brightest bin, the 1$\sigma$ statistical uncertainties in $\log\phi$ estimated from bootstrapping are approximately 0.15 dex, 0.15 dex, and 0.3 dex at $z=8$, 10, and 12, respectively.. Unlike some other theoretical predictions (e.g., Mason et al., 2023; Yung et al., 2023), for which a clear tension with the spec-$z$ lower bounds exists without modifications, our bursty-case predictions do not require any additional tuning of UV variability or production efficiency to match observations. Despite uncertainties associated with the treatment of dust, this good agreement implies that the UVLFs observed by JWST at $z\gtrsim 10$ are consistent with generally “normal” SFE and UV production efficiency as predicted by the FIRE-2 simulations. As demonstrated in Ma et al. (2018b), the relation between $M_{*}$ and $M_{\mathrm{h}}$ in these simulations is broadly consistent with extrapolations from lower $z$ where empirical analyses show that the SFE is strongly suppressed by stellar feedback in low-mass halos (Behroozi et al., 2013; Tacchella et al., 2018).

Second, in the contrasting, “smoothed” SFH scenario, a clear deficit of UV-bright galaxies is seen as a result of suppressed up-scattering in $M_{\mathrm{UV}}$ of low-mass halos when the SFR is averaged over a long timescale $\tau_{\mathrm{SF}}$. The underestimated abundance of UV-bright galaxies reveals the important role played by the burstiness of star formation in determining the number statistics of galaxies at cosmic dawn. As also shown by the comparison of different $\tau_{\mathrm{SF}}$ values at $z=10$, smoothed SFHs with $\tau_{\mathrm{SF}}\gtrsim 100\,$Myr result in bright-end UVLFs that are too steep compared with observations, especially the photometrically selected samples, for which the bright-end UVLF can be underpredicted at $>2\sigma$ level in some cases (Castellano et al., 2023; Donnan et al., 2023). At $z=12$, predictions of the smoothed SFH are in tension with even the most conservative lower limits derived from only the spectroscopically confirmed samples (Harikane et al., 2023a). It is therefore clear that the UVLF serves as a useful probe of the burstiness in the SFH, as been noted in e.g., Furlanetto & Mirocha (2022) and Shen et al. (2023), although in practice it can be challenging to extract the burstiness information from only the UVLF measurements (see Section 4). The overall shallower $M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation when the SFH is smoothed also leads to slightly steeper slope at the faint end, although the effect is much smaller than the suppression at the bright end.

The binned UVLFs without dust attenuation extracted from our simulations at $z=8$, 10, and 12 are summarized in Table 1. As has been demonstrated in Figure 1, dust attenuation only modestly affects the UVLF at the very bright end, reducing $\phi$ (in the brightest bin) by approximately 0.4, 0.25, and 0.01 dex at $z=8$, 10, and 12, respectively. The binning scheme is chosen such that the brightest $M_{\mathrm{UV}}$ bin contains more than ten simulated galaxies for robust statistics. Meanwhile, we fit the dust-free UVLF at $8\leq z\leq 12$ assuming a universal double-power law (DPL) in $M_{\mathrm{UV}}$,

We specify the redshift-dependent DPL parameters $\phi_{*}$, $M^{*}_{\mathrm{UV}}$, $\alpha$, and $\beta$ in the form of a single power law as $\phi_{*}(z)=\phi_{*,0}[(1+z)/10]^{\phi_{*,1}}$, $M^{*}_{\mathrm{UV}}(z)=M^{*,0}_{\mathrm{UV}}[(1+z)/10]^{M^{*,1}_{\mathrm{UV}}}$, $\alpha_{*}(z)=\alpha_{*,0}[(1+z)/10]^{\alpha_{*,1}}$, and $\beta_{*}(z)=\beta_{*,0}[(1+z)/10]^{\beta_{*,1}}$, where the best-fit parameters are found to be $\phi_{*,0}=-2.01$, $\phi_{*,1}=0.68$, $M^{*,0}_{\mathrm{UV}}=-17.26$, $M^{*,1}_{\mathrm{UV}}=-0.08$, $\alpha_{*,0}=-0.31$, $\alpha_{*,1}=-0.93$, $\beta_{*,0}=0.68$, and $\alpha_{*,1}=0.93$.

Figure 2 shows a comparison between the binned and best-fit UVLFs predicted by our simulations and other theoretical predictions in the literature based on cosmological hydrodynamical simulations (Ocvirk et al., 2020; Vijayan et al., 2021; Dawoodbhoy et al., 2023; Kannan et al., 2023; Wilkins et al., 2023). Overall, our predicted UVLFs show a weaker redshift evolution beyond $z=8$ compared with the predictions from the MillenniumTNG (Kannan et al., 2023) and CoDa II (Ocvirk et al., 2020; Dawoodbhoy et al., 2023) simulations, which results in a higher abundance of bright ($M_{\mathrm{UV}}\lesssim-20$) galaxies at $z\gtrsim 10$. Our bright-end predictions are generally comparable to those from the FLARES simulations (Vijayan et al., 2021; Wilkins et al., 2023) in both normalization and slope, despite the vastly different nature of the simulations and methods to evaluate the UVLF. It is noteworthy, though, that the FIRE-2 simulations analyzed in this work have significantly higher resolution ($m_{\mathrm{b}}\approx 7\times 10^{3}\,M_{\odot}$ in FIRE-2 vs. $m_{\mathrm{b}}\approx 2\times 10^{6}\,M_{\odot}$ in FLARES), which allows us to predict the UVLFs at $8\leq z\leq 12$ down to $M_{\mathrm{UV}}\sim-10$ vs. the FLARES predictions down to $M_{\mathrm{UV}}\sim-18$. We have also verified that UVLFs in this work and from Ma et al. (2018b, 2019) are in good agreement in the overlapping regime.

By integrating the predicted UVLFs, we can derive the UV luminosity density, $\rho_{\mathrm{UV}}$, as a function of time, which traces the cosmic star formation rate density (SFRD). Since at $z\gtrsim 10$ only the brightest end ($M_{\mathrm{UV}}\ll M_{\mathrm{UV,*}}$) of the UVLF has been probed, we follow Harikane et al. (2023a) to compare the UV luminosity density contributed by galaxies brighter than $M_{\mathrm{UV}}=-18$, namely $\rho_{\mathrm{UV,bright}}=\rho_{\mathrm{UV}}(M_{\mathrm{UV}}<-18)$, which corresponds to the contribution from halos with $M_{\mathrm{h}}\gtrsim 10^{10}\,M_{\odot}$ at $z=10$. The unconstrained contribution by fainter, lower-mass galaxies is highly sensitive to the faint-end slope of the UVLF and might even outweigh $\rho_{\mathrm{UV,bright}}$ (Sun & Furlanetto, 2016), but the comparison restricted to $M_{\mathrm{UV}}<-18$ galaxies still serves as a useful test of the overall abundance of bright, massive galaxies and their SFE at cosmic dawn⁷⁷7Results from this work, Harikane et al. (2018, 2023b), and Bouwens et al. (2023) are integrated down to $M_{\mathrm{UV,lim}}=-18$, whereas the rest are down to $M_{\mathrm{UV,lim}}=-17$. Figure 3 thus shows conservatively that our simulations without smoothing predict enough total UV emission compared with observations, regardless of the modest difference in $M_{\mathrm{UV,lim}}$..

Figure 3 shows a comparison of the cumulative UV luminosity density between the dust-attenuated predictions from our simulations and a compilation of constraints from observations and theoretical forecasts in the literature. Throughout, dust-attenuated predictions from models/simulations (curves) are compared with observations (data points), which are dust-uncorrected. Over $8\leq z\leq 12$, dust-attenuated luminosity densities predicted by our simulations without smoothing the SFH are fully consistent with observations of both photometric galaxy candidates and spectroscopically-confirmed galaxies that provide firm lower limits. Due to the integrated nature of $\rho_{\mathrm{UV}}$, the “smoothed” case appears more consistent with the spec-$z$-only lower limits here than at the bright end of the UVLF as shown in Figure 2. In both cases with dust attenuation, a power-law evolution of $\rho_{\mathrm{UV}}\propto(1+z)^{-0.3}$ over $8\leq z\leq 12$ is implied, which appears more gradual compared with the predictions by some previously proposed semi-analytic/semi-empirical models, such as in Mason et al. (2015) and Harikane et al. (2018).