# Stochastic star formation and the abundance of $z>10$ UV-bright galaxies

Andrey Kravtsov1,2,3*  — and Vasily Belokurov4*  —
1Department of Astronomy and Astrophysics, The University of Chicago, Chicago,
IL 60637 USA
2Kavli Institute for Cosmological Physics, The University of Chicago, Chicago,
IL 60637 USA
3Enrico Fermi Institute, The University of Chicago, Chicago, IL 60637
4Institute of Astronomy, Madingley Rd, Cambridge, CB3 0HA, UK
<EMAIL_ADDRESS>
###### Abstract

We use a well-motivated galaxy formation framework to predict stellar masses,
star formation rates (SFR), and ultraviolet (UV) luminosities of galaxy
populations at redshifts $z\in 5-16$, taking into account stochasticity of SFR
in a controlled manner. We demonstrate that the model can match observational
estimates of UV luminosity functions (LFs) at $5<z<10$ with a modest level of
SFR stochasticity, resulting in the scatter of absolute UV luminosity at a
given halo mass of $\sigma_{M_{\rm UV}}\approx 0.75$. To match the observed UV
LFs at $z\approx 11-13$ and $z\approx 16$ the SFR stochasticity should
increase so that $\sigma_{M_{\rm UV}}\approx 1-1.3$ and $\approx 2$,
respectively. Model galaxies at $z\approx 11-13$ have stellar masses and SFRs
in good agreement with existing measurements. The median fraction of the
baryon budget that was converted into stars, $f_{\star}$, is only
$f_{\star}\approx 0.005-0.05$, but a small fraction of galaxies at $z=16$ have
$f_{\star}>1$ indicating that SFR stochasticity cannot be higher. We discuss
several testable consequences of the increased SFR stochasticity at $z>10$.
The increase of SFR stochasticity with increasing $z$, for example, prevents
steepening of UV LF and even results in some flattening of UV LF at $z\gtrsim
13$. The median stellar ages of model galaxies at $z\approx 11-16$ are
predicted to decrease from $\approx 20-30$ Myr for $M_{\rm UV}\gtrsim-21$
galaxies to $\approx 5-10$ Myr for brighter ones. Likewise, the scatter in
median stellar age is predicted to decrease with increasing luminosity. The
scatter in the ratio of star formation rates averaged over 10 and 100 Myr
should increase with redshift. Fluctuations of ionizing flux should increase
at $z>10$ resulting in the increasing scatter in the line fluxes and their
ratios for the lines sensitive to ionization parameter.

###### keywords:

high-redshift – galaxies: luminosity function, mass function – galaxies:
stellar content – galaxies: formation – galaxies

##  1 Introduction

The era of the James Webb Space Telescope (JWST) is revolutionizing our
understanding of the earliest stages of galaxy evolution due to the unique
capabilities of its instruments to probe stellar populations at $z>10$ (e.g.,
see Robertson, 2022, for a review). Indeed, deep observations with JWST during
the first two years of its operation detected significant numbers of $z>9$
galaxies bright in the rest-frame ultraviolet (UV) band (e.g., Naidu et al.,
2022; Castellano et al., 2022; Finkelstein et al., 2023; Harikane et al.,
2023, 2024; Donnan et al., 2023; Robertson et al., 2023b; McLeod et al., 2024;
Casey et al., 2024; Hainline et al., 2024), which opened a window into the
earliest stages of galaxy evolution. Deep surveys are also finding increasing
numbers of fainter galaxies (e.g., Hainline et al., 2024; Austin et al., 2023;
Simmonds et al., 2024; Atek et al., 2023, 2024) which will lead to reliable
measurements of the UV luminosity functions form (UV LFs) at these high
redshifts.

Abundance of bright galaxies was estimated in several recent studies
(Finkelstein et al., 2023; Wilkins et al., 2023; Robertson et al., 2023a; Yung
et al., 2024), although robust measurement of redshift is often difficult and
thus debate about the abundance of the brightest galaxies at $z>7$ is ongoing
(e.g., Desprez et al., 2023). Comparisons of these estimates with model
predictions led to the conclusions that theoretical models significantly
underestimate the abundance of UV bright galaxies at $z>10$.

Theoretical analyses indicated that to match observed estimates of the
high-$z$ UVLF requires high star formation efficiency and suppressed feedback
in the gas-rich, starbursting environments of these galaxies (Boylan-Kolchin,
2023; Dekel et al., 2023; Qin et al., 2023) or a top-heavy IMF (Inayoshi et
al., 2022; Trinca et al., 2024; Yung et al., 2024; Wang et al., 2023a; Qin et
al., 2023; Ventura et al., 2024). At the brightest galaxy tail, Ferrara et al.
(2023) suggest that a rapid decrease of dust obscuration may play a role in
reconciling models with observations. Furthermore, the star formation rate
(SFR) and UV luminosties of $z>9$ galaxies may be overestimated due to the
contribution of active galactic nuclei (e.g., D’Silva et al., 2023).

At the same time, Keller et al. (2023), McCaffrey et al. (2023), Sun et al.
(2023) and Kannan et al. (2023) showed that results of high-resolution galaxy
formation simulations are generally consistent with observed abundance
galaxies in the JWST surveys at $z\sim 10-12$ for the range of luminosities
and star formation rates (SFRs) where simulations and observations overlap.
This seemingly contradictory conclusion provides a hint of why some of the
simpler models may have failed to match the abundance of the high-$z$
galaxies. Star formation in high-resolution simulations of galaxies at these
redshifts tends to be very bursty, while simple models only account for the
average star formation of galaxies of a given mass.

Indeed, several studies used simple analytic models for star formation in dark
matter halos to point out that strong UV variability that can be induced by
the variations in the SFR and/or initial mass function (IMF) of stars can help
explain the abundance of the UV-bright galaxies by up-scattering low-mass
galaxies that happen to be UV-bright at a given time (Ren et al., 2019;
Mirocha & Furlanetto, 2023; Mason et al., 2023; Shen et al., 2023; Muñoz et
al., 2023). At the same time, Sun et al. (2023) showed that the SFR
stochasticity in the FIRE-2 galaxy formation simulations is sufficient to
produce UV LF consistent with measurements of the UV LF at $z\approx 10-12$.

On the observational side, recent JWST observations of spectroscopically
confirmed galaxies also indicate that star formation in galaxies of $\lesssim
10^{9}\,M_{\odot}$ at $z>5$ is very bursty (e.g., Looser et al., 2023; Atek et
al., 2023), which is consistent with the inferences from the broad-band JWST
NIRCam photometry modelling (e.g., Dressler et al., 2024; Endsley et al.,
2023; Simmonds et al., 2024; Asada et al., 2024; Ciesla et al., 2023; Cole et
al., 2023). These observational studies show that burstiness of star formation
increases with increasing redshift and with decreasing galaxy luminosity.
Thus, both theoretical models and observations that SFR is quite stochastic
during the early stages of galaxy evolution.

At the same time, it is not yet clear whether model predictions agree with
observations in detail or what the physical drivers of the SFR stochasticity
are. In galaxy formation simulations compelling mechanisms have been
identified to explain the burstiness of SFR in dwarf galaxies. For example, Su
et al. (2018) showed that feedback due to individual supernovae and hypernovae
leads to significant bursts of SFR in the FIRE-2 simulations. Sugimura et al.
(2024), on the other hand, show that suppression of cooling due to H2
dissociation by the FUV radiation delays star formation until gas accumulates
and becomes gravitationally unstable, which leads to a strong burst.
Subsequent removal of gas by the feedback-driven outflow leads to a lull in
star formation until the cycle repeats. In both studies, however, the effects
are specific to dwarf galaxies and were shown to be significant in halos of
mass $M_{200}\lesssim 10^{9}\,M_{\odot}$.

In higher mass halos bursty SFR can be driven by the imbalance between gas
accretion time scale and time scales for star formation and feedback that
prevents equilibrium star formation regime in the interstellar medium.
Faucher-Giguère (2018) argued that all high-redshift galaxies should be in
such a regime due to high accretion and merger rates expected during early
stages of galaxy evolution.

In this paper, we use a regulator-type galaxy formation model (Kravtsov &
Manwadkar, 2022; Manwadkar & Kravtsov, 2022) to explore the role of bursty
star formation in shaping galaxy UV LFs at $z\sim 5-16$. The goal is to test
and elucidate the conclusions of the previous studies discussed above. We show
that our results are consistent with the results of Sun et al. (2023) who
showed that the level of SFR stochasticity in FIRE-2 simulations is sufficient
to match observed UV LFs at $z\approx 8-12$. Here we extend the analyses to a
wider range of redshifts of $z\in 5-16$. We show that the amounts of
stochasticity required to match observed UV LF at $z>10$ are smaller than
deduced in some of the previous studies. At the same time, the stochasticity
is required to increase with increasing redshifts from $\sigma_{M_{\rm
UV}}\approx 0.75$ at $z<9$ to $\sigma_{M_{\rm UV}}\approx 1.-1.3$ at $z\approx
10-13$ and to $\sigma_{M_{\rm UV}}\approx 2$ at $z=16$. We further show that
SFR stochasticity plays the key role in shaping galaxy UV LF.

The paper is organized as follows. In Section 2 we describe the methodology
used to construct a representative sample of halos at $z\in 5-16$ and galaxy
formation model used to populate the halos with galaxies. We also present an
approach we used to add stochasticity to the baseline model SFR in a
controlled manner. We present our results in Section 3 and discuss them in
Section 4. Our results and conclusions are summarized in Section 5. In the
appendix A we present a demonstration of how UV LF changes for different
assumptions of star formation depletion time, while in the appendix B we
present comparisons of the halo mass accretion rate approximation we use to
model the evolution of halo mass in our framework to other calibrations. We
show that our approximation is accurate across a wide redshift range of
$0<z<20$ and that small differences in calibrations of the accretion rate at
$z>10$ have a negligible effect on the predicted galaxy UV LF at $z\approx
12.$

In calculations throughout this paper we assume flat $\Lambda$+Cold Dark
Matter ($\Lambda$CDM) cosmology with the mean density of matter in units of
the critical density of $\Omega_{\rm m}=0.32$, mean density of baryons of
$\Omega_{\rm b}=0.045$, Hubble constant of $H_{0}=67.11\,\rm
km\,s^{-1}\,Mpc^{-1}$, the amplitude of fluctuations within tophat spheres of
$R=8h^{-1}$ Mpc of $\sigma_{8}=0.82$, and the primordial slope of the power
spectrum of $n_{s}=0.95$. Halo virial masses throughout this study are defined
within the radius enclosing density contrast of 200 relative to the critical
density at the corresponding redshift and are denoted as $M_{\rm 200c}$.

##  2 Methods

The galaxy formation framework we use in this study is applied to predict
galaxy population properties for representative samples of model galaxies at
all relevant luminosities and redshifts. This is done using samples of halos
that follow the expected halo mass function at each considered redshift and
halo mass evolution tracks constructed using an accurate approximation for the
halo mass accretion rate, as described in Section 2.1.

The galaxy formation model used here was shown to reproduce observed
properties of $\lesssim L_{\star}$ galaxies at $z=0$ down to the faintest
ultra-faint dwarf galaxies (Kravtsov & Manwadkar, 2022; Manwadkar & Kravtsov,
2022; Kravtsov & Wu, 2023) and is briefly outlined in Section 2.2 below. The
basic model is supplemented with a simple prescription for a controlled
stochasticity in the star formation rate described in Pan & Kravtsov (2023,
see Section 2.3 below) to explore its effects on the UV LF.

###  2.1 Halo evolution model

To model the luminosity function of galaxies accurately at a given redshift
$z$ we first construct large samples of model halos using the following
procedure. We use an accurate cubic spline approximation of the cumulative
halo mass function computed using Tinker et al. (2008) approximation of halo
mass function and the inverse transform sampling method to generate a random
sample of halo masses in a given volume. The halo sample thus accurately
follows the Tinker et al. (2008) halo mass function by construction.
Experiments showed that halo samples in the mass range $\log_{10}M_{\rm
200c}\in[9,13.5]$ and cubic volumes of size $L\in[500-2000]h^{-1}\,\rm Mpc$
provide an accurate representation of the abundance of halos hosting galaxies
at $z\in[5-16]$.

Large volumes are required to model the bright end of the galaxy luminosity
function accurately. Given that the halo mass function is steep, halo samples
over the entire mass range in such volumes are large. For example, the number
of halos in the above mass range in a $L=500h^{-1}$ Mpc cube at $z=5$ is
$\approx 1.74\times 10^{8}$ and this number is dominated by halos near the
lower mass limit of the sample. Modeling galaxy formation in all of the halos,
however, is not necessary. Instead, we select a random fraction of halos for
modeling which is given as a function of halo mass $M_{\rm 200c}$ using:

$f(M_{\rm 200c})=n\left(\frac{M_{\rm 200c}}{10^{9}\,M_{\odot}}\right)^{\eta},$
(1)

where $n\approx 10^{-5}-5\times 10^{-6}$ and $\eta\approx 1.5-2.35$ provide
sufficiently large samples of model galaxies to reliably measure luminosity
functions at different $z\in[5,16]$. When the luminosity function is computed
at a given $z$, each model galaxy is weighted by $f^{-1}(M_{\rm 200c})$. This
procedure ensures good sampling of all relevant decades of the luminosity
function, while keeping the number of model galaxies reasonably small.

Once a halo sample is drawn at a given redshift, $z_{\rm f}$, we use each halo
mass as a starting point and construct a halo mass evolution track. The track
is computed by integrating the equation of halo mass evolution $\dot{M}_{\rm
200c}=\mu(M_{\rm 200c})$ back in time using an accurate approximation for the
mean halo mass accretion rate $\dot{M}_{\rm 200c}$:

$\dot{M}_{\rm 200c}=0.3606\,M_{\odot}\,{\rm Gyr}^{-1}\,M_{\rm
200c}(t)^{1.091}\,t^{-1.8},$ (2)

where halo mass $M_{\rm h}$ is in $M_{\odot}$ and time $t$ is in Gyrs. We
present comparisons of this simple equation to other approximations proposed
in the literature, including the high-$z$ calibration of Yung et al. (2023),
in the Appendix B. This approximation for mass accretion rate was derived
using analyses of the mass evolution histories of halos formed in cosmological
$\Lambda$CDM simulations. The comparisons show that the above approximation
accurately describes the mean accretion rate of halos of a given mass and has
a convenient simple form for integration in time. We also show in the Appendix
B (see Figure 17 and related discussion) that small differences between
different approximations have a negligible effect on the model UV luminosity
function.

Note that the equation 2 approximates the average rate of accretion, while
mass assembly histories of halos of a given mass have a scatter which
corresponds to the scatter in the halo mass accretion rates. Although it would
be straightforward to model this scatter, we choose not to do so in this study
because the scatter in the star formation in our fiducial model is completely
dominated by the scatter introduced by the stochastic SFR model described
below. The scatter in the halo mass accretion rate produces a sub-dominant
scatter of only $\lesssim 0.1$ dex in stellar mass at a given halo mass.

In our analysis we thus consider two model variants: the base model that has
no sources of SFR scatter and which describes only the average evolution of
galaxy properties in halos of a given mass and the model in which different
amounts of SFR stochasticity are introduced.

Figure 1: Left panel: Example star formation history of one of the model
galaxies produced by the fiducial model (blue long dashed line) and star
formation histories with different levels of stochasticity added
($\sigma_{\Delta}=0.08,0.15,0.30$ shown as violet, purple, and blue lines,
respectively). Right panel: the median power spectrum densities of model
galaxies in the fiducial model (thick blue long dashed line) and of the model
star formation histories with different levels of stochasticity added (solid
violet, purple and blue lines). The dotted, short-dash, and dot-dashed color
lines show power spectra of $z>5$ galaxies in high-resolution zoom-in
simulations by the FIRE-2 project (Ma et al., 2018, 2019, 2020) and by Garcia
et al. (2023, see text for details). Although the scatter is significant,
there is a general trend of the PSD amplitude increasing with decreasing
stellar mass of galaxies.

###  2.2 The model of galaxy formation

The halo mass evolution tracks computed as described above used as the basis
for the GRUMPY galaxy formation model (Kravtsov & Manwadkar, 2022) built upon
a regulator-type galaxy formation framework (e.g., Krumholz & Dekel, 2012;
Lilly et al., 2013; Feldmann, 2013). The model solves a system of differential
equations that describe the evolution of gas mass, stellar mass, size, and
stellar and gas-phase metallicities. It also includes galactic outflows, a
model for the gaseous disk and its size, molecular hydrogen mass, star
formation, and more. The model parameters used in this study are identical to
those in the fiducial model of Manwadkar & Kravtsov (2022), which was shown to
reproduce properties of observed $L\lesssim L_{\star}$ galaxies at $z=0$ down
to the ultra-faint galaxies (cf. also Kravtsov & Wu, 2023).

The GRUMPY model assumes that at all times the ISM follows the exponential
radial gas profile, $\Sigma_{g}(R)$ with a half-mass radius proportional to
the parent halo virial lradius. The mass of molecular hydrogen treated as a
proxy for star-forming gas is estimated as

$M_{\rm H_{2}}=2\pi\int_{0}^{\rm{R_{HI}}}f_{\rm
H_{2}}(\Sigma_{g})\Sigma_{g}(R)RdR$ (3)

where $R_{\rm{HI}}$ is the radius corresponding to the assumed self-shielding
threshold and $f_{\rm H_{2}}$ is the molecular fraction computed using current
gas metallicity and UV flux proportional to the mean star formation rate using
the model of Gnedin & Draine (2014). We note in passing that modeling of
molecular gas at low metallicities is uncertain and star formation may not be
traced by molecular gas (Polzin et al., 2023). We use the above prescription
for molecular gas as a proxy for star-forming gas. As we noted, the model with
this star formation prescription reproduces a broad range of properties of
$L\lesssim L_{\star}$ galaxies and thus likely captures realistic scaling of
the star-forming gas fraction with metallicity.

The star formation rate is computed using $M_{\rm H_{2}}$ assuming a constant
depletion time of molecular gas, $\tau_{\rm sf}$, and instantaneous recycling
of the gas:

$\dot{\mathcal{M}}_{\star}=(1-\mathcal{R})\frac{M_{\rm H_{2}}}{\tau_{\rm sf}}$
(4)

where $\mathcal{R}=0.44$ is the gas fraction that is returned to the ISM via
mass loss in winds and supernovae explosions assuming the Chabrier (2003)
initial mass function of stars (e.g., Leitner & Kravtsov, 2011; Vincenzo et
al., 2016). We use the fiducial value of $\tau_{\rm sf}=2$ Gyr consistent with
measurements in nearby galaxies (e.g., Bigiel et al., 2008), but we explore
models in which $\tau_{\rm sf}$ decreases with increasing redshift, as
indicated by observations of high $z$ galaxies. In addition, $\tau_{\rm sf}$
varies significantly from galaxy to galaxy. We do not model this scatter
explicitly, because we will include different sources of scatter in SFR in a
controlled manner using an explicit model of SFR stochasticity described
below.

The model includes outflows of the ISM gas that are assumed to be proportional
to the mean SFR, $\dot{M}_{\rm out}=\eta_{w}M_{\rm H_{2}}/\tau_{\rm sf}$ with
the mass-loading factor $\eta_{w}$ dependent on the current stellar mass of
the galaxy in a way expected in the energy-driven wind models. With the
outflow inclusion, the model successfully reproduces the general trends
observed in the SFR–stellar mass relation, mass–metallicity relation, and
several other properties of $<L_{\star}$ galaxies (Manwadkar & Kravtsov,
2022). However, it does not include modeling of the processes that can result
in significant SFR stochasticity, such as formation and destruction of
individual star-forming regions (e.g., Tacchella et al., 2020; Iyer et al.,
2020; Sugimura et al., 2024). As described above, we also do not introduce
stochasticity due to scatter in halo assembly histories and depletion times.
Instead, we introduce the collective SFR stochasticity due to these various
processes in a controlled manner using the PSD formalism described below.

###  2.3 The SFR stochasticity model

In addition to the evolution of galaxy properties computed in the base model,
we model stochasticity of star formation relative to the mean SFR estimated by
the base model by producing correlated Gaussian random number for each time
step during galaxy evolution integration, as described in Pan & Kravtsov
(2023). The numbers are generated assuming a given power spectral density
(PSD) of the time series in the Fourier domain. Namely, the mean SFR in the
model at a time moment $t_{n}$ is perturbed as

${\rm SFR}_{\rm stochastic}={\rm SFR}_{\rm mean}\times 10^{\Delta},$ (5)

where $\Delta$ is a correlated random number drawn from the Gaussian pdf with
zero mean and unit variance and multiplied by $\sqrt{P(k)}=\sqrt{{\rm
PSD}(f_{k})/T}$, where wavenumber $k$ corresponding to frequency $f_{k}$ is
defined as $k=f_{k}T$ and where $T$ is the duration of galaxy evolution track.

We follow Caplar & Tacchella (2019) and use the PSD of the form

${\rm PSD}(f)=\frac{\sigma^{2}_{\Delta}}{1+(\tau_{\rm break}f)^{\alpha}},$ (6)

where $\sigma_{\Delta}$ characterizes the amplitude of the SFR variability
over long time scales and $\tau_{\rm break}$ characterizes the timescale over
which the random numbers are effectively uncorrelated. Parameter $\alpha$
controls the slope of the PSD at high frequencies (short time scales).

In our models, we fix the slope $\alpha$ and $\tau_{\rm break}$ to the values
$\alpha=2$ and $\tau_{\rm break}=100$ Myr that are physically motivated by the
time scales of gas evolution and star formation in giant molecular clouds in a
typical ISM (see Tacchella et al., 2020, for a detailed discussion). In the
analyses presented below we only vary the amplitude of the PSD,
$\sigma_{\Delta}$. Note that the SFR stochasticity is introduced after the
model evolution is computed to predict the distribution of star formation
rates and corresponding UV luminosities. It does not affect the outflows in
the model that are assumed to be proportional to the mean star formation rate
computed by the model.

Physically, the stochasticity can be driven by fluctuations in gas mass
accretion rate, feedback-driven outflows (Tacchella et al., 2020), and other
modulation processes, such as far UV suppression of molecular gas formation in
high-density regions (Sugimura et al., 2024) and imbalance between gas
accretion rate and the time scale of star formation and feedback in gas-rich
galaxies (Faucher-Giguère, 2018).

The right panel of Figure 1 shows examples of the PSDs of the SFRs modeled at
$z>5$ using our galaxy formation model with stochasticity corresponding to
different $\sigma_{\Delta}$ values with the thick solid (for
$\sigma_{\Delta}>0$) and thick blue dashed (for $\sigma_{\Delta}=0$) line. The
corresponding star formation histories are shown in the left panel of the
figure. Note that the PSD for $\sigma_{\Delta}=0$ simply reflects the average
shape of the star formation history ${\rm SFR}(t)$. The thin dotted and dot-
dashed lines show PSDs of star formation histories of the individual galaxies
in the FIRE-2 simulations of galaxies forming in halos
$M_{200}=10^{9}\,M_{\odot}$ at $z=5$ (Ma et al., 2018, 2019, 2020) and in
galaxies at $z>7$ in the simulations of Garcia et al. (2023). The thin short-
dashed lines show median PSDs of star formation histories of galaxies forming
in halos of mass $M=10^{10}\,M_{\odot}$, $10^{11}\,M_{\odot}$,
$10^{12}\,M_{\odot}$ at $z=0$ (Wetzel et al., 2023). The left panel of the
figure illustrates the variation of SFR for one of the model galaxies with
these different levels of SFR stochasticity.

Figure 1 shows that the stochasticity of the SFR in these high-resolution
zoom-in simulations can be well described by the model PSD described above
with $\sigma_{\Delta}\approx 0.08-0.4$. In the analyses below will use the
values in this range and will show that the stochasticity of this level is
sufficient to match the observed UV LFs at $z\in[5,16]$.

Figure 2: UV luminosity functions at redshifts $z=5,6,7,8,9$ (from the top
down) measured from the HST and JWST observations (points) and predicted by
the model (lines). LF at each redshift is offset by the factor of $10^{6-z}$
for clarity. The dashed lines show model UV LF with the fiducial parameters
and depletion time $\tau_{\rm sf}=2$ Gyr when stochasticity is not included,
while solid lines show results for the same model but with SFR stochasticity
added with the PSD parameters $\sigma_{\Delta}=0.08$, $\alpha=2$, $\tau_{\rm
break}=100$ Myr.

###  2.4 Computing the UV luminosities of model galaxies

To compute the monochromatic UV luminosity of model galaxies at
$\lambda=1500\,\rm\mathring{A}$, we tabulated the grid of luminosities,
$L_{1500}$, for stellar populations of a given age and metallicity using the
Flexible Stellar Population Synthesis model v3.0 (FSPS, Conroy et al., 2009;
Conroy & Gunn, 2010b) and its Python bindings, Python-FSPS. The table was used
to construct an accurate bivariate cubic spline approximation to compute
$L_{1500}$ for stellar populations of a given age and metallicity produced by
the model.

Specifically, we use finely spaced time outputs of the models to compute the
integral $L_{1500}$ due to all stars formed by the current time using the
evolution of stellar mass and stellar metallicity computed by the model. We do
not include any nebular emission in the calculation of $L_{1500}$. We also do
not include dust effects, which are expected to affect the brightest galaxies
in the UV at $z\lesssim 7$. At $z\gtrsim 7$ these effects are expected to be
minor as indicated by the values of the UV spectral slope $\beta$ in observed
galaxies and by theoretical models (e.g., Dayal et al., 2022; Lewis et al.,
2023) and are considerably smaller than the uncertainty of the UV luminosity
function predictions due to uncertainty in the amount of SFR stochasticity.

##  3 Results

Table 1: SFR stochasticity normalization and resulting scatter of $M_{\rm UV}$ in different model parameter choices at | $z=10.7$ | $z=12$ | $z=13$ | $z=16$
---|---|---|---|---
$\tau_{\rm sf}=2$ Gyr, $\sigma_{\Delta}=0.08$, UV LF is not matched
$\sigma_{M_{\rm UV}}$ | 0.69 | 0.83 | 0.84 | 0.85
$\tau_{\rm sf}=2$ Gyr, variable $\sigma_{\Delta}$, UV LF is matched
$\sigma_{\Delta}$ | 0.15 | 0.15 | 0.16 | 0.25
$\sigma_{M_{\rm UV}}$ | 1.1 | 1.4 | 1.5 | 2.2
variable $\tau_{\rm sf}$, variable $\sigma_{\Delta}$, UV LF is matched
$\tau_{\rm sf}$ (Gyr) | 1.0 | 0.5 | 0.5 | 0.25
$\sigma_{\Delta}$ | 0.13 | 0.13 | 0.15 | 0.22
$\sigma_{M_{\rm UV}}$ | 0.97 | 1.2 | 1.4 | 1.9
Figure 3: UV luminosity functions at redshifts $z=10.7,12,13,16$ (from the top
left panel clockwise) measured using JWST observations (points) and predicted
by the model (lines). The dashed lines show model of UV LF with the fiducial
parameters and decreasing depletion time $\tau_{\rm sf}=2$ Gyrs with
stochasticity amplitude $\sigma_{\Delta}=0.08$ (PSD normalization), while dot-
dashed lines show results for the same model but with $\sigma_{\Delta}$
increasing with increasing redshift from $\sigma_{\Delta}\approx 0.15$ at
$z=10.7-13$ to $\sigma_{\Delta}\approx 0.25$ at $z=16$ (see Table 1). The
solid lines show a similar model but with star formation depletion time
decreasing from $\tau_{\rm sf}=1$ Gyr at $z=10.7$ to $\tau_{\rm sf}=0.25$ Gyr
at $z=16$; in this model, the amount of SFR stochasticity ($\sigma_{\Delta}$
values) required to match observed LFs are lower (see Table 1). The dotted
line in the $z=12$ panel (upper right corner) shows the UV LF estimate using
the FIRE-2 simulations from Sun et al. 2023. The observational estimates using
spectroscopic samples are shown by the brighter red, while the estimates from
photometric samples are shown with fainter symbols. Specific measurements used
at each redshift are: $z=10.7$ Harikane et al. (2023, 2024, diamonds), Donnan
et al. (2023, squares), Leung et al. (2023, circles), Chemerynska et al.
(2023, downward triangles); $z=12$ Robertson et al. (2023a, stars), Bouwens et
al. (2023b, crosses), Finkelstein et al. (2022, triangle); $z=13$ Bouwens et
al. (2023a, crosses), Donnan et al. (2023, squares), Robertson et al. (2023a,
stars); $z=16$ Harikane et al. (2023, faint circle and upper limit), Harikane
et al. (2024, upper limit at $M_{\rm UV}=-21.9$), Bouwens et al. (2023a,
cross).

###  3.1 UV luminosity functions at $5<z<16$

We first present the results of the galaxy formation model described above for
the UV luminosity functions of galaxies at $z<10$ in Figure 2. The dashed
lines show predictions of the model with star formation depletion time of
$\tau_{\rm sf}=2$ Gyrs and without any sources of stochasticity of SFR, while
the solid lines show the model with a moderate amount of stochasticity added
to star formation rate with the power spectrum amplitude of
$\sigma_{\Delta}=0.08$ (see eq. 6 in Section in 2.3). This level of SFR
stochasticity is comparable to the model shown by the magenta line (second
lowest thick line) in the right panel of Figure 1 and is on the low side of
the stochasticity observed in high-resolution galaxy formation simulations
shown in that figure. It corresponds to the scatter of $M_{\rm UV}$ at a given
halo mass of $\sigma_{M_{\rm UV}}\approx 0.75$, which is comparable to
$\sigma_{M_{\rm UV}}$ in high-resolution simulations at these redshifts
(Pallottini & Ferrara, 2023; Sun et al., 2023). It is also consistent with the
conclusions of Shen et al. (2023) about the amount of $M_{\rm UV}$ scatter
required to reproduce UV LFs at $z\in[5,9]$.

Figure 2 shows that the model with moderate SFR stochasticity provides a
remarkably good match to observational measurements across the entire
luminosity range probed by observations at these redshifts. Some minor
discrepancies at luminosities $M_{1500}\lesssim-21$ and $z<7$ are likely due
to dust attenuation not accounted in the model LF. As expected, the effect of
the SFR stochasticity on the LF increases at brighter luminosities where LF is
steeper. Likewise, the effect of stochasticity increases with increasing
redshift for the same reason.

Figure 3 compares the same model with SFR stochasticity of
$\sigma_{\Delta}=0.08$ by the long-dashed lines and observational estimates of
UV LF at $z>10$. At these redshifts, the model underpredicts the abundance of
observed galaxies and the difference increases with increasing redshift. Such
discrepancy of the model predictions and observations at these high redshifts
is generic to most pre-JWST galaxy formation models (e.g., Finkelstein et al.,
2023; Harikane et al., 2023, 2024; Wilkins et al., 2023; Yung et al., 2024;
Chemerynska et al., 2023). However, as we discussed above in Section 2.3, the
levels of stochasticity in very high-resolution simulations of galaxy
evolution at $z>5$ is large and most previous models of UV LF do not take this
into account. Indeed, Figure 3 shows that if the SFR stochasticity amplitude
is increased to $\sigma_{\Delta}=0.15-0.25$, comparable to those measured in
the simulations of Ma et al. (2018, 2019, 2020) and Garcia et al. (2023), the
model (shown by the dot-dashed lines) can match the observational estimates of
UV LFs.

Figure 3 shows that the increase of SFR stochasticity with increasing $z$
required to match observations leads to the gradual flattening of the UV LF at
higher $z$, which is a distinct prediction of this model. This can be seen in
Figure 4 which shows model UV luminosity functions at redshifts
$z=10.7,12,13,16$ side-by-side. It shows that the increase of SFR
stochasticity with increasing $z$ results in the LF that does not steepen at
higher $z$ but becomes somewhat shallower at $z>12$.

Figure 4: Model UV luminosity functions at redshift $z=10.7,12,13,16$ (solid,
dot-dashed, short-dashed, and long-dashed lines, respectively) shown side-by-
side to illustrate the fact that increase of SFR stochasticity with increasing
$z$ results in the LF that does not steepen at higher $z$ but becomes somewhat
shallower at $z>12$.

The current observational measurements are consistent with such flattening
(see, e.g., Fig. 13 in Harikane et al., 2023), but measurements at $z\gtrsim
12$ are still too uncertain to probe the slope of the UV LF reliably. This
prediction can be tested as the number of reliably detected galaxies at these
redshifts increases. The model also predicts that the luminosity function
extends to higher luminosities than the range probed by current observations.
Future observations should thus find galaxies with higher UV luminosities than
the maximum luminosities of the current samples.

Note that the effect of SFR stochasticity on the model UV LF is twofold. The
main effect is the upscattering of numerous lower-mass galaxies to higher
luminosities due to SFR and UV luminosity bursts. Additionally, the lognormal
scatter of SFR assumed in the model increases the mean SFR in galaxies of a
given mass. The latter effect can be compensated by other parameters of the
model affecting the average level of star formation. In particular, the
depletion time of star-forming gas, $\tau_{\rm sf}$, directly affects the star
formation efficiency, and thus SFR, for a given amount of gas. The solid lines
in Figure 3 show UV LFs of the model in which $\tau_{\rm sf}$ is allowed to be
different at different $z$: specifically, $\tau_{\rm sf}=1$ Gyrs at $z=11$ and
$z=12$, $\tau_{\rm sf}=0.5$ Gyrs at $z=11$ and $\tau_{\rm sf}=0.25$ Gyr at
$z=16$. This model can produce UV LF similar to that of the model with
$\tau_{\rm sf}=2$ Gyr, but with smaller SFR stochasticity. The specific
normalization of the power spectrum of stochastic SFR required to match
observed UV LF at each redshift is presented in Table 1.

Although these depletion times values are smaller than the value of $\tau_{\rm
sf}=2$ Gyr typical for $z=0$ galaxies (e.g., Bigiel et al., 2008), they are
comparable to those observed in starbursting galaxies (e.g., Kennicutt &
Evans, 2012; Huang & Kauffmann, 2014; Hunt et al., 2015; Díaz-García & Knapen,
2020) and generally in galaxies at higher $z$ (e.g., Tacconi et al., 2020). In
particular, direct estimates of molecular depletion time in galaxies at
$z\gtrsim 6$ indicate values of measured $\approx 0.05-1$ Gyr (Pavesi et al.,
2019; Vallini et al., 2024).

Nevertheless, the effect of decreasing $\tau_{\rm sf}$ alone on the UV LF is
insufficient. As we show in Appendix A, a factor of ten decrease of depletion
time does increase the amplitude of LF by an order of magnitude but this is
still two orders of magnitude short of matching the observational estimates of
LF amplitude at this redshift. Thus, in our model LF at $z\gtrsim 13$ cannot
be matched by a realistic increase of star formation efficiency. This is
because the increase of star formation rate is partly offset by the increase
in the gas outflow rate, which is proportional to the SFR.

###  3.2 UV luminosity–stellar mass relation of $z>10$ model galaxies

Figure 5 shows the UV absolute magnitude–stellar mass correlation of model
galaxies at $z=10.7,13,16$ in models with and without SFR stochasticity. The
median correlation can be described approximately by a relation
$\log_{10}M_{\star}=\log_{10}M_{-20\star}+s(M_{1500}+20)$, where
$\log_{10}M_{-20\star}=8.2,8.04,7.8$ and $s=0.25,0.33,0.4$ at $z=10.7,13,16$,
respectively. The scatter around the correlation in the model with SFR
stochasticity is moderate because scatter in stellar mass and absolute
magnitude are correlated. Nevertheless, the scatter is both non-negligible and
asymmetric which can make a conversion of the UV luminosity to $M_{\star}$
subject to biases unless a correct form of the pdf $p(M_{\star}|M_{1500})$ is
adopted and scatter is taken into account properly.

Figure 5: UV absolute magnitude–stellar mass correlation of model galaxies at
$z=10.7,13,16$. The dashed lines show the relation in the model without SFR
stochasticity. The solid lines show the median $M_{\star}$ in bins of
$M_{1500}$ in the model with SFR stochasticity that matches observed UV LF at
$z>10$ (shown by solid lines in Fig. 3). The shaded regions show the $68.2$
and $95.4$ percentiles of the stellar mass distribution.

Figure 6: Star formation rate as a function of stellar mass in model and
observed galaxies at $z=10.7,12,13,16$. The solid lavender lines show the
median relations for the fiducial model but with SFR stochasticity added with
the PSD parameters $\alpha=2$, $\tau_{\rm break}=100$ Myr and
$\sigma_{\Delta}=0.13,0.13,0.15,0.22$ at $z=10.7,12,13,16$. The UV LF for this
model is shown by the solid lines in Figure 3. The $1\sigma$ and $2\sigma$
scatter in SFR at a given $M_{\star}$ is shown by the dark and light blue
bands. The three dotted lines show linear relation ${\rm SFR}={\rm sSFR}\times
M_{\star}$ for ${\rm sSFR}=10^{-7},10^{-8},10^{-9}\,\rm yr^{-1}$. The lower
subpanels in each of the four redshift panels show the $1\sigma$ and $2\sigma$
scatter in $\log_{10}\rm SFR$ around the median. The red points in the first
three panels show available observational measurements for galaxies at the
corresponding redshifts (squares, Atek et al. 2023; triangle in the $z=10.7$
panel, Roberts-Borsani et al. 2023; circles, Franco et al. 2023; stars,
Robertson et al. 2023a; upward triangle in the $z=12$ panel, Castellano et al.
2024; downward triangles in the $z=12$ and $z=13$ panels, Wang et al. 2023b);
diamonds, Morishita et al. 2024 and Morishita & Stiavelli 2023).

.

###  3.3 Star formation rates of $z>10$ model galaxies

Figure 6 shows the evolution of the $M_{\star}$-SFR relation from $z=16$ to
$z=11$ in the model with variable $\sigma_{\Delta}$ and $\tau_{\rm sf}$. At
$z=11-13$ we also plot star formation rates and stellar masses inferred for
individual observed galaxies at the corresponding redshifts. The figure shows
that the model with stochastic SFR is broadly consistent with observations,
which means that model galaxies of a given stellar mass have correct star
formation rates. It also shows that matching UV LFs at $z>10$ requires that
the scatter of SFR at higher $z$ increases dramatically with increasing
redshift. At $z=10.7$ the scatter in $\log_{10}\rm SFR$ at a given $M_{\star}$
is $\approx 0.4-0.5$ dex, which is consistent with a recent observational
estimate at that redshift by Cole et al. (2023). It remains to be seen if the
required scatter in the SFR is consistent with observations at higher $z$. As
noted above for the UV LF, the model predicts the existence of galaxies with
larger stellar masses and SFRs than in the population detected by the current
surveys. These galaxies are rare but should be found as the volume probed by
surveys increases.

Figure 7: The ratio of star formation rates averaged over the past 10 Myr and
100 Myr as a function of the absolute magnitude at $1500\,{\rm\AA}$ for model
galaxies at $z=10.7$, $z=13$, $z=16$. The dashed lines show the median $\rm
SFR_{10}/SFR_{100}$ in the model without SFR stochasticity, while solid lines
show the median for the fiducial model with SFR stochasticity at the level
that matches observational estimates of UV LF at these redshifts. The shaded
bands show the $68\%$ and $95\%$ of the distribution for the latter model.

.

Figure 8: The ratio of the Hα luminosity to the luminosity per unit frequency
at $\lambda=1500\,\rm\AA$ as a function of the absolute magnitude at
$1500\,{\rm\AA}$ for model galaxies at $z=10.7$, $z=13$, $z=16$. Hα luminosity
is proportional to the flux of hydrogen ionizing photons produced by the
youngest stellar populations of model galaxies (age $\lesssim 3-5$ Myr), while
UV luminosity reflects star formation history on the time scale of tens of
Myr. The dashed lines show the median ratio in the model without SFR
stochasticity, while solid lines show the median ratio for the fiducial model
with SFR stochasticity at the level that matches observational estimates of UV
LF at these redshifts. The shaded bands show the $68\%$ and $95\%$ of the
distribution for the latter model. Figure 9: The median age of stars of model
galaxies, $\tau_{\rm sf,50}$, in millions of years as a function of their
$\lambda=1500\,{\rm\AA}$ absolute magnitude, $M_{1500}$ at redshifts
$z=10.7,13,16$. The dashed lines show the median $\tau_{\rm sf,50}$ in the
model without SFR stochasticity. The stellar ages in this model are quite
small due to rapid growth of galaxies at these redshifts. The solid lines show
the median $\tau_{\rm sf,50}$ in the model with SFR stochasticity that matches
observed UV LF at $z>10$ (shown by solid lines in Fig. 3). The shaded regions
show the $68.2$ and $95.4$ percentiles of the $\tau_{\rm sf,50}$ distribution.

An indicator of the star formation histories and stochasticity of the current
star formation is the ratio of star formation rates averaged over the past 10
Myr and 100 Myr. Figure 7 shows this ratio as a function of the absolute
magnitude at $\lambda=1500\,\rm\mathring{A}$ for model galaxies at $z=10.7$,
$z=13$, $z=16$. The dashed lines in the figure panels show the median ratio in
the model without SFR stochasticity which have values of $\rm
SFR_{10}/SFR_{100}\approx 2$ at $z=10.7$ and $5$ at $z=16$. These values are
due to the rapid growth and rapidly increasing SFR of model galaxies at these
redshifts, which make most recent star formation larger than the star
formation rate averaged over a longer period.

Figure 7 shows that SFR stochasticity results in a significant scatter around
the median $\rm SFR_{10}/SFR_{100}$ ratio which is predicted to increase with
increasing $z$. Cole et al. (2023) estimated the scatter of this ratio of
$\approx 0.4-0.5$ dex at $9<z<12$, consistent with the scatter exhibited by
galaxies in our stochastic SFR model. We note that the scatter is significant
both for the values of $\rm SFR_{10}/SFR_{100}$ above and below the median
value. SFR stochasticity thus produces galaxies with enhanced and depressed
recent star formation at the same $M_{1500}$. This means there should be a
wide range of ionizing luminosities produced by galaxies of a given UV
luminosity. This prediction is potentially testable using observational
estimates of the Hα and $L_{1500}$ luminosities. Figure 8 shows this ratio
predicted by our model at the same three redshifts. The Hα luminosity was
computed assuming Case B ionization as (e.g., eq. 2 in Kennicutt, 1998)

$L_{\rm H_{\alpha}}=1.37\times 10^{-12}\dot{N}_{\rm ion}\ \rm erg\,s^{-1},$
(7)

where $\dot{N}_{\rm ion}$ is the emission rate of photons with
$\lambda<912\,\rm\mathring{A}$ per second. This flux was computed self-
consistently for model galaxies using their star formation histories and the
tables of ionizing flux for single-age stellar populations with binaries from
the BPASS stellar population synthesis package version 2.3 (Byrne et al.,
2022).

In addition to modulating the $L_{\alpha}/L_{1500}$ ratio, the scatter in the
ionizing flux should result in a wide range of line ratios for the ions
sensitive to ionizing radiation, such as OII (Kewley et al., 2019). The highly
bursty star formation will result in a wide range of the [OIII]/[OII] ratios
in $z>10$ galaxies.

Finally, Figure 9 shows that the median age of stars in model galaxies,
$\tau_{\rm sf,50}$, is $\approx 10-40$ Myr at $z\in[11,16]$ broadly consistent
with estimates for observed galaxies at these redshifts (e.g., Franco et al.,
2023; Morishita & Stiavelli, 2023; Robertson et al., 2023b, a; Curtis-Lake et
al., 2023; Carnall et al., 2023; Zavala et al., 2024). Note that the median
ages are $\approx 10-30$ Myr even without SFR stochasticity due to extremely
rapid growth of galaxies at these redshifts.

The SFR stochasticity introduces a distinct trend of decreasing median
$\tau_{\rm sf,50}$ and decreasing scatter with increasing galaxy luminosity.
This is because with the levels of SFR stochasticity required to match
observed UV LF most of the bright galaxies become bright due to a strong burst
of star formation, which reduces the median age of their stars. Indications of
such a trend are found for galaxies at $z<9$ (Endsley et al., 2023).

These trends represent a potential test of the scenario in which UV LF
function at $z>10$ is shaped by the significant SFR stochasticity. Figure 9
also shows that the scatter of stellar median ages is quite large at
$M_{1500}\lesssim-20$ due to the variation of periods of specific bursts and
lulls in star formation.

###  3.4 Galaxy-halo connection

Given that the model with SFR stochasticity reproduces the observed estimates
of the UV LF across $z\in 11-16$, it is interesting to consider the galaxy-
halo connection it implies in the form of the UV luminosity-halo mass
relation.

Figure 10 shows the UV absolute magnitude at $\lambda=1500\,\rm\mathring{A}$
as a function of halo mass $M_{\rm 200c}$ at redshifts $z=10.7,14,16$. For
reference, the dotted line in each panel shows the relation

$M_{1500}=-17-3.3\left(\log_{10}M_{200}-10\right),$ (8)

which approximates the median relation of the stochastic SFR model at
$z=10.7$.

The median luminosities at a given halo mass increase with increasing redshift
and their scatter increases due to the assumed increase in the SFR
stochasticity at higher $z$. The corresponding increase in the luminosity of
galaxies in the model without SFR stochasticity is much milder. In the
presence of significant scatter the distribution pdf $p(M_{1500}|M_{\rm
200c})$ is generally different from the pdf $p(M_{\rm 200c}|M_{1500})$ and
Figure 11 shows that this is the case here. The relation between the median
$M_{\rm 200c}$ in bins of $M_{1500}$ and absolute magnitude is no longer
linear and it becomes progressively flatter with increasing $z$. This
flattening is the consequence of the lack of halos of mass $M_{\rm
200c}\gtrsim 3\times 10^{11}\,M_{\odot}$ that are rare at these redshifts.

Figure 11 shows that most model galaxies in the observed luminosity range
($-21\lesssim M_{1500}\lesssim-17$) are hosted by halos in a narrow mass range
of less than two orders of magnitude. The stellar masses of model galaxies in
such halos are in the range $M_{\star}\approx 10^{6}-5\times
10^{9}\,M_{\odot}$. The median stellar mass – halo mass relation for this
range of $M_{\star}$ is well described by a power law

$M_{\star}=f_{\star}\times 0.15\times
10^{10}\,M_{\odot}\left(\frac{M_{200}}{10^{10}\,M_{\odot}}\right)^{1.4}$ (9)

where $f_{\star}$ is the fraction of the baryon budget that was converted into
stars in halos of mass $M_{200}$.

Figure 12 shows that the median $f_{\star}$ varies from $f_{\star}\approx
0.005-0.02$ at $z=10.7,12,13$ to $f_{\star}\approx 0.01-0.05$ at $z=16$. Thus,
halos are expected to convert only a few percent of their available baryons
into stars on average to match observed UV LFs at $z>10$ in agreement with
conclusions of Mason et al. (2023).

However, the scatter around the median increases dramatically with increasing
redshift, as required by the condition to match the observational estimate
using UV LFs. The tail of the distribution at $z=13$ reaches $f_{\star}\approx
0.1-0.3$. These values are still very reasonable and are comparable to the
stellar fractions of $L_{\star}$ galaxies at $z=0$, such as the Milky Way. At
$z=16$, however, the $95.4$ percentile of the distribution extends to
$f_{\star}\approx 1$, which indicates that the SFR stochasticity of the model
at this redshift is as high as it can be.

Figure 10: Absolute magnitude at $1500\,\rm\mathring{A}$, $M_{1500}$, as a
function of halo mass, $M_{\rm 200c}$, of model galaxies at $z=10.7,13,16$.
The dashed lines show the results of the model with the fiducial parameters
and decreasing molecular depletion time with increasing $z$ from $\tau_{\rm
sf}=1$ Gyr at $z=10.7$ to $0.25$ Gyr at $z=16$ and when stochasticity is not
included. The scatter in this model shown by the narrow dark blue bands around
the dashed lines is negligible (a small amount of scatter results from the
different gas disk sizes assumed in the model for each galaxy). The solid
lavender lines show the median relations for the same model but with SFR
stochasticity added with the PSD parameters $\alpha=2$, $\tau_{\rm break}=100$
Myr and $\sigma_{\Delta}=0.13,0.15,0.22$ at $z=10.7,13,16$. The UV LF for this
model is shown by the solid lines in Figure 3.

.

Figure 11: Distribution of halo mass, $M_{\rm 200c}$, of model galaxies as a
function of absolute magnitude at $1500\,\rm\mathring{A}$, $M_{1500}$, at
$z=10.7,13,16$. The plots show the results of the model with the fiducial
parameters and decreasing molecular depletion time with increasing $z$ from
$\tau_{\rm sf}=1$ Gyr at $z=10.7$ to $0.25$ Gyr at $z=16$ and SFR
stochasticity added with the PSD parameters $\alpha=2$, $\tau_{\rm break}=100$
Myr and $\sigma_{\Delta}=0.13,0.15,0.22$ at $z=10.7,13,16$ with which this
model matches observed UV LF. The solid lavender lines show the median $M_{\rm
200c}$ in bins of $M_{1500}$, while the $68.2$ and $95.4$ percentiles of the
distribution are shown by the darker and lighter shaded bands. The figure
shows that most model galaxies in the observed luminosity range
($M_{1500}\lesssim-17$) are hosted by halos of a fairly narrow mass range. The
average mass corresponding to observed luminosities decreases with increasing
$z$. The increasing SFR stochasticity at higher $z$ leads to significant
flattening of the relation.

.

Figure 12: Fraction of the baryon mass that should be available to halos,
$\Omega_{\rm b}/\Omega_{\rm m}\,M_{200}$, converted into stars in model
galaxies at redshifts $z=10.7$, $z=13$, and $z=16$. The solid lines show the
median of the distribution, while dark and light-shaded regions show the
$68.2\%$ and $95.6\%$ regions of the distribution of $f_{\star}$ at a given
$M_{200}$. The horizontal dashed line shows $f_{\star}=1$ when galaxies
convert their entire available baryon budget into stars.

##  4 Discussion

###  4.1 Comparisons with previous studies and alternative scenarios to
explain $z>10$ UV LF

During the last two years, several scenarios have been put forth to explain
the high abundance of UV-bright galaxies at $z>10$ (see, e.g., Section 6.1 of
Harikane et al., 2023).

For example, Dekel et al. (2023) posited an idea that galaxies at these
redshifts undergo starbursts in which usual feedback processes do not have the
time to operate and gas is converted into stars in bursts with high
efficiency. Several studies have also proposed to explain the high abundance
of UV bright galaxies by invoking an increase of star formation efficiency
with increasing redshift (Boylan-Kolchin, 2023; Qin et al., 2023; Mason et
al., 2023).

Our results indicate that increased star formation efficiency is not necessary
to explain UV LF at $z\leq 13$ when SFR is stochastic, as the highest
fractions of baryons converted into stars are reasonable $f_{\star}\lesssim
0.3$. Moreover, our model results show that if outflows are as efficient as
they need to be to explain properties of $z=0$ dwarf galaxies higher star
formation efficiency results in only a modest increase of UV LF amplitude.
This is because higher star formation efficiency is partly compensated by
higher gas outflow rate, which results in smaller gas mass. In simulations
with weak feedback such regulation does not occur and higher UV luminosities
can be achieved by increasing star formation efficiency. However, galaxies in
such simulations may have incorrect properties at lower redshifts.

Another discussed idea is the top-heavy IMF that increases the UV luminosity
per unit stellar mass formed (Inayoshi et al., 2022; Trinca et al., 2024; Yung
et al., 2024; Ventura et al., 2024). There are good reasons to expect top-
heavy IMF at these high redshifts, where the minimum temperature of the ISM
gas is higher due to higher cosmic microwave background temperature (e.g.,
Jermyn et al., 2018; Chon et al., 2022) or other star-formation related
processes (e.g., Riaz et al., 2020; Steinhardt et al., 2022) and where a
significant fraction of stars may be forming in massive bound clusters (e.g.,
Belokurov & Kravtsov, 2023; Adamo et al., 2024; Mowla et al., 2024). The
latter are environments conducive to the formation of massive and supermassive
stars via stellar collisions (e.g., Katz et al., 2015; Gieles et al., 2018;
Katz, 2019). Rasmussen Cueto et al. (2023), however, explored evolving IMF
within their galaxy formation model and found that top-heavy IMF by itself is
unlikely to explain the bright end of UV LF at $z>10$ because the increase of
the UV luminosity per unit mass of stars born is counteracted by the
suppression of star formation due to increased feedback.

The stochasticity of SFR was proposed as a natural explanation of the high
amplitude of UV LF at $z>10$ in several recent studies (Mason et al., 2023;
Shen et al., 2023; Sun et al., 2023; Muñoz et al., 2023). Qualitatively, our
results are consistent with their conclusions. However, we find that the
amount of stochasticity required to match the observed UVLFs at $z=11-13$ is
considerably smaller than was concluded in most previous studies. Shen et al.
(2023) presented results of a simple model assuming a constant efficiency of
star formation and SFR proportional to the halo mass accretion rate. They
concluded that to explain UV LF at $z\lesssim 9$ a modest level of scatter in
$M_{\rm UV}$ of $\sigma_{M_{\rm UV}}\approx 0.75$ consistent with our finding.
At $z>10$ Shen et al. (2023) deduce a much larger required scatter of
$\sigma_{M_{\rm UV}}\approx 2.0$ at $z=12$ and $\sigma_{M_{\rm UV}}\approx
2.5$ at $z=16$. Mason et al. (2023), on the other hand, concluded that
$\sigma_{M_{\rm UV}}\approx 1.5$ is sufficient to match UV LFs at $z\approx
11-13$, the number closer to our estimate.

We show that at $z=11-13$ the observed UV LF can be reproduced with a
considerably smaller scatter of $\sigma_{M_{\rm UV}}\approx 1-1.4$. As shown
in Figure 13, these values of $\sigma_{M_{\rm UV}}$ are consistent with the
scatter measured in the FIRE-2 simulations. The figure shows the standard
deviation of the UV absolute magnitude at $\lambda=1500\,\rm\mathring{A}$,
$\sigma_{M_{\rm UV}}$ as a function of the logarithm of halo mass in the
FIRE-2 simulations (Sun et al., 2023) and SPHINX galaxy formation simulations
(Katz et al., 2023a). The solid horizontal lines show $\sigma_{M_{\rm UV}}$
for our fiducial model that matches the observational estimates of UV LF at
different redshifts with the extent of the lines showing the effective halo
mass range contributing to galaxies in the observed luminosity range. Note
that our model matching observational estimates of UV LF corresponds to
$\sigma_{M_{\rm UV}}$ similar to those measured in the FIRE-2 simulations and
we also reproduce the UV LF estimated using FIRE-2 simulations (see the $z=12$
panel of Figure 3).

The scatter in the SPHINX simulation is consistent with the FIRE-2 simulations
at $M_{\rm 200c}\gtrsim 3\times 10^{10}\,M_{\odot}$ but is smaller at smaller
masses. At $z=10.7$ the level of $\sigma_{M_{\rm UV}}$ in the SPHINX is close
to the level required to match observed UV LF, but is smaller at higher $z$.
The scatter in the SPHINX does reach values of $\sigma_{M_{\rm UV}}\approx
1.5-2$ but at $\log_{10}M_{\rm 200c}\lesssim 9$. There are thus differences in
the predicted degree of SFR stochasticity in current simulations, likely due
to different implementations of star formation and feedback.

Figure 13: Standard deviation of the UV absolute magnitude at
$\lambda=1500\,{\rm\AA}$, $\sigma_{M_{1500}}$ as a function of the logarithm
of halo mass in the FIRE-2 simulations (Sun et al., 2023) and in the SPHINX
galaxy formation simulations (Katz et al., 2023a). The solid horizontal lines
show $\sigma_{M_{1500}}$ for our fiducial model that matches the observational
estimates of UV LF at different redshifts. The horizontal extent of each line
shows the $1\sigma$ spread of halo mass that host model galaxies contributing
to UV LF probed by current observations at each $z$.

At $z=16$ we find that $\sigma_{M_{\rm UV}}\approx 1.9-2$ is required by
existing estimates of galaxy abundance at $-22\lesssim M_{\rm UV}\lesssim-20$.
This is also smaller than the previous estimate of Shen et al. (2023). We
believe the difference arises because Shen et al. (2023) first compute the UV
LF without scatter and then convolve it with a Gaussian of a given width.
However, in the presence of a large scatter this procedure is not accurate and
UV LF needs to be computed by a convolution of the halo mass function with the
Gaussian pdf $p(M_{\rm UV}|M_{\rm 200c})$. If we convolve our $z=16$ UV LF in
the model without SFR stochasticity with a Gaussian as was done by Shen et al.
(2023) we do indeed find that $\sigma_{M_{\rm UV}}\approx 2.5$ is needed to
match observational UV LF estimate. Yet, the average scatter in our fiducial
model that matches observations at a given halo mass is only $\sigma_{M_{\rm
UV}}\approx 1.9$.

In fact, we show in Section 3.4 that $\sigma_{M_{\rm UV}}\approx 2$ at $z=16$
we deduce should be considered an upper limit because it requires that a small
fraction of the model galaxies should form more stars than the mass of
available baryons. This shows that higher levels of SFR stochasticity than
those used in our model are unrealistic.

###  4.2 Other sources of scatter of UV luminosity

SFR stochasticity is not the only process that can contribute to the scatter
of UV luminosity. Fluctuations of $M_{\rm UV}$ can be enhanced for a given SFR
stochasticity if the initial mass function of stars becomes top-heavy during
SFR bursts. The surface density of star formation, $\Sigma_{\rm SFR}$ should
increase dramatically during such bursts and it is known that in low-$z$
galaxies the fraction of stars formed in massive bound clusters increases with
increasing $\Sigma_{\rm SFR}$ (e.g., see Adamo et al., 2020, for a review). At
$\Sigma_{\rm SFR}\gtrsim 10-30\,M_{\odot}\,\rm yr^{-1}\,kpc^{-2}$,
characteristic of extreme starbursts at $z=0$ and observed in a significant
fraction of galaxies at $z>10$ (e.g., Morishita & Stiavelli, 2023; Robertson
et al., 2023b; Morishita et al., 2024), $\gtrsim 50\%$ of star formation is
expected to form in massive bound clusters according to the observed trend for
local galaxies. As $\Sigma_{\rm SFR}$ can increase to even higher values
during strong SFR bursts (in fact $\Sigma_{\rm SFR}$ for many $z>10$ are lower
limits due to limited angular resolution of observations), close to $100\%$ of
stars may form in massive bound clusters.

This is now directly observed in galaxies at $z\gtrsim 7$ (Adamo et al., 2024;
Mowla et al., 2024; Castellano et al., 2024) and in the low-metallicity
stellar population of the Milky Way that should have formed at similar
redshifts (Belokurov & Kravtsov, 2023). As noted above, massive bound clusters
are environments conducive to the formation of massive and supermassive stars
(e.g., Gieles et al., 2018). When close to $100\%$ of stars form in massive
clusters the IMF can become top-heavy due to the formation of such stars. The
impact of the SFR bursts on the UV luminosity can thus be considerably
enhanced compared to the standard IMF assumed in our analysis. On the
observational side, potential signatures of top-heavy IMF have been indeed
indicated by recent observational analyses (Cameron et al., 2023).

###  4.3 Observational indications of bursty star formation at $z>5$

A growing number of observational analyses indicate a large degree of
burstiness in $z>7$ galaxies. Endsley et al. (2023) show that the population
of galaxies at $z\approx 7-9$ exhibits large fluctuations of SFR and that
brightest galaxies tend to be undergoing a starburst. Simmonds et al. (2024)
show that bursty dwarf galaxies at $4\lesssim z\lesssim 9$ contribute more
ionizing photons than non-bursty galaxies of the same stellar mass and thus
likely play a key role in reionizing the Universe. Moreover, they find strong
evidence that ionizing luminosity increases with decreasing galaxy luminosity
and increasing redshift indicating that burstiness increases with increasing
$z$.

Likewise, Ciesla et al. (2023) used JWST JADES survey data to estimate the
amount of SFR stochasticity in galaxies at $6<z<12$ and found that
stochasticity is not only ubiquitous but is increasing towards higher
redshifts (see also Cole et al., 2023). They find that at $z>9$ up to $\approx
90\%$ of galaxies of mass $M_{\star}\gtrsim 10^{9}\,M_{\odot}$ show signs of
bursty star formation rate in the past 100 Myr, while this fraction drops to
$\approx 15\%$ at $z<7$. This is consistent with our finding that the amount
of SFR stochasticity required to match the UVLFs at $z>5$ decreases rapidly
around $z\approx 10$.

Ciesla et al. (2023) also estimated that the dispersion of UV luminosity
$\sigma_{M_{\rm UV}}$ stays approximately constant as a function of $z$ at
$\sigma_{M_{\rm UV}}\approx 1.2$. Thus, their results imply that it is the
fraction of starbursting galaxies that rapidly declines with decreasing
redshift, not the dispersion $\sigma_{M_{\rm UV}}$. We note, however, that at
$z>9$ where almost all galaxies are bursty according to their analysis, our
results indicate that $\sigma_{M_{\rm UV}}\approx 1.2$ is sufficient to match
the observed UVLF at $z\lesssim 13$, especially if lower molecular gas
depletion time is considered (see Figure 3).

###  4.4 Observational signatures of bursty star formation

Our results show that SFR stochasticity at the level required to match UV LF
estimates at $z>10$ has several potentially testable predictions. First, this
scenario predicts that UV LF does not steepen with increasing $z$ at $z>10$,
but gets even somewhat shallower at $z>12$ (see Section 3.1 and Figures 3 and
4). This prediction should be testable with future JWST observations.

The model also predicts a decrease in the median stellar age of the stellar
population with increasing UV luminosity (see also Ren et al., 2019; Mason et
al., 2023). A significant scatter of median stellar ages for galaxies of
$M_{1500}\gtrsim-21$ is expected with the significant decrease of scatter for
brighter galaxies (see Fig. 9). The increasing SFR stochasticity at $z>10$
should result in a larger fraction of galaxies with extremely young ages
$\lesssim 5-10$ Myr.

Significant SFR stochasticity modifies the distribution of virial masses of
halos that host observed galaxies within a given $M_{1500}$ range (see Fig.
11). In this scenario, UV-bright galaxies occupy predominantly low-mass halos
and become UV bright during SF bursts. In contrast, in the scenario in which
the UV bright galaxies are bright due to high SF efficiency, the bright
galaxies should occupy massive halos. As pointed out by Muñoz et al. (2023),
this difference should in principle be manifested in the different clustering
amplitude of UV-bright galaxies predicted by these two scenarios.

SFR stochasticity results in a substantial scatter of the star formation rates
averaged over 10 and 100 Myr, SFR10/SFR100 (Fig. 7). The scatter of this ratio
in model galaxies is dominated by the fluctuations of SFR10, which is in
qualitative agreement with measurements of Cole et al. (2023) at $z\approx
10-11$. To explain the observational estimates of UV LF at higher $z$, SFR
stochasticity must increase at $z>11$ and this leads to the increase the
scatter around the median SFR10/SFR100 from $\approx 0.5$ dex at $z=10.7$ to
$\approx 0.7$ dex at $z=13$ and $\approx 1$ dex at $z=16$.

A related quantity is the ratio of Hα to UV luminosity (Fig. 8). This ratio is
essentially the ionization production efficiency, $\xi_{\rm ion}$:

$\xi_{\rm ion}=\frac{\dot{N}_{\rm ion}}{L_{1500}}\ {\rm
erg^{-1}\,Hz}\approx\frac{7.28\times 10^{11}\,L({\rm H}_{\alpha})}{L_{1500}}\
\rm erg^{-1}\,Hz$ (10)

where $\dot{N}_{\rm ion}$ is the number of the Lyman continuum photons emitted
by a galaxy per second, and $L_{1500}$ is the monochromatic luminosity at
$\lambda=1500\,{\rm\AA}$ in units of $\rm erg\,s^{-1}Hz^{-1}$. Note that the
$L({\rm H}_{\alpha})/L_{1500}$ ratio predicted by our model assumes the
standard IMF (the IMF adopted in BPASS v2.3). This ratio should be different
if top-heavy IMF is used to explain the high UV luminosities of $z>10$
galaxies.

The fluctuations of $\xi_{\rm ion}$ due to SFR burstiness will result in the
fluctuations of the ionization factor $U$, which is proportional to
$\dot{N}_{\rm ion}$. The stochastic SFR scenario predicts that the brightest
galaxies $M_{1500}\lesssim-21$ should have larger Hα equivalent widths due to
the smaller median age of their stellar population. The ratio of Hα to UV
luminosities should have a smaller scatter than in fainter galaxies (see Fig.
8). Indications of such trend have been found in $z<9$ galaxies (Endsley et
al., 2023).

The fluctuations of the ionization factor will also result in fluctuations of
the line flux rations for the lines sensitive to ionizing radiation, such as
[OIII]/[OII], [NIII]/[NII], etc (see, e.g., Fig. 7 in Kewley et al., 2019).
Given that SFR stochasticity is required to increase with redshift the effects
of the ionizing flux fluctuations should also increase. We thus can expect an
increase in the typical [OIII]/[OII] line flux ratios with increasing $z$.

Also, strong fluctuations of SFR and corresponding ionizing flux are expected
to lead to suppression of $158\,\mu\rm m$ [CII] emission in high-$z$ galaxies
(Katz et al., 2023b), which means that increasing SFR burstiness should result
in a lower number of [CII] detections at $z>10$ in ngVLA observations (Carilli
et al., 2018). Katz et al. (2023b) show that combined information from the
[CII], Hα, and $L_{1500}$ can be used to identify galaxies that experienced
extreme SFR fluctuations recently.

##  5 Summary and conclusions

We use the GRUMPY galaxy formation framework to predict stellar masses, star
formation rates, and UV luminosities of galaxy populations at $z\in 5-16$. The
galaxy formation model was previously shown to reproduce observed properties
of $\lesssim L_{\star}$ galaxies at $z=0$ down to the faintest ultra-faint
dwarf galaxies (Kravtsov & Manwadkar, 2022; Manwadkar & Kravtsov, 2022;
Kravtsov & Wu, 2023). Here we use the model with halo mass tracks obtained by
integration of a new simple approximation to the halo mass accretion rate
accurate across a wide redshift interval of $0<z<20$ (eq. 2). We add SFR
stochasticity to the star formation rate predicted by the model to explore its
effects on the observed UV LF in a controlled manner and its potential
observational implications. Our results and conclusions are as follows.

  * 1.

We show that the model can match observational estimates of observed UV LFs at
$5<z<10$ with a modest level of SFR stochasticity resulting in $\sigma_{M_{\rm
UV}}\approx 0.75$ (Figure 2). Such SFR stochasticity is consistent with the
expectation of high-resolution galaxy formation simulations.

  * 2.

To match the observed UV LFs at $z\approx 11-13$ the SFR stochasticity should
increase by a factor of $\approx 1.5-2$ producing $\sigma_{M_{\rm UV}}\approx
1-1.3$ (see Table 1). We show that such a stochasticity level is consistent
with the SFR stochasticity observed in high-resolution zoom-in simulations and
is significantly smaller than the level of stochasticity deduced by previous
studies.

  * 3.

We show that if the SFR stochasticity increases with increasing redshift at
the rate suggested by our results, the observed UV LF at $z>10$ does not
steepen and even becomes somewhat flatter at $z>12$ (Figures 3 and 4).

  * 4.

We show that our model matching $z>10$ UV LFs is in good agreement with
existing measurements of the SFR and stellar masses of galaxies at these
redshifts (Figure 6). The model also predicts the scatter of the star
formation rates averaged over 10 and 100 Myr, SFR10/SFR100 in agreement with
the current estimate at $z\approx 10-11$. This scatter is predicted to
increase at higher $z$ due to the increase of SFR stochasticity required to
match the observational estimates of UV LFs at $z>10$ (Figure 7).

  * 5.

The typical median age of stellar populations in the model with SFR
stochasticity is $\approx 20-30$ Myr with a large scatter and the model
predicts a systematic decrease of the median stellar age and scatter around it
at the absolute magnitudes of $M_{1500}\lesssim-21$ (see Figure 9).

  * 6.

The scatter of short-term SFR and stellar ages results in a substantial
scatter in the ionizing flux at a given UV luminosity. This has several
consequences that can be potentially probed by observations. For example, the
increase of SFR stochasticity with redshift results in a larger fraction of
galaxies experiencing a strong burst or a significant lull in their star
formation. The associated fluctuations of the ionization parameter will result
in the increasing scatter in the line fluxes and their ratios with increasing
$z$ for the lines sensitive to this parameter (see Discussion in Section 4.3).

  * 7.

According to our model, galaxies of $M_{\rm UV}\in[-17,-22]$ are hosted by
halos of mass $M_{\rm 200c}\approx 5\times 10^{9}-3\times 10^{11}\,M_{\odot}$
at $z\approx 10.7$ but are hosted in halos of progressively smaller mass at
higher $z$.

  * 8.

Stellar masses of model galaxies at observed luminosities at $z>10$ are in the
range $\approx 5\times 10^{6}-5\times 10^{9}\,M_{\odot}$ with the median
stellar mass-halo mass relation well described by $M_{\star}=f_{\star}\times
0.15\times 10^{10}\,M_{\odot}M_{200,10}^{1.4}$, where $M_{200,10}$ is halo
mass in units of $10^{10}M_{\odot}$ and $f_{\star}$ is the fraction of the
baryon budget that was converted into stars in halos of mass $M_{200}$, which
varies from $f_{\star}\approx 0.005-0.03$ at $z=10.7-13$ to $f_{\star}\approx
0.01-0.05$ at $z\approx 16$ (Figure 12). On average the model galaxies thus
convert only a few percent of the available baryon budget into stars and
extremely high star formation efficiency is not required to match observed UV
LF at $z\leq 13$.

  * 9.

The scatter of $f_{\star}$ is substantial but most galaxies have
$f_{\star}\lesssim 0.4$ at $z\approx 11-13$. At $z=16$, on the other hand, the
level of SFR stochasticity required to match the current estimate of the
abundance of UV bright galaxies predicts that a small fraction of galaxies
have $f_{\star}>1$ (Figure 12). Thus, the level of SFR stochasticity required
to match current estimates of UV LF at $z\approx 16$ are as high as they can
be.

Our results and results of previous related theoretical studies strongly
motivate exploration of the physical processes that can drive the burstiness
of star formation in both dwarf and massive galaxies at $z\gtrsim 5-10$ as
these processes may be prevalent during early stages of galaxy evolution and
can thus shape the properties of low-metallicity stellar populations of
galaxies. On the observational side, several different observational probes
discussed in Section 4.3 can be used to measure or constrain the level of SFR
stochasticity at $z>10$, thereby providing a test of this explanation for high
abundance of UV-bright galaxies at these redshifts.

## Acknowledgments

We are grateful to Brant Robertson, Harley Katz, and the UChicago structure
formation group for useful discussions during this project. We are also
grateful to Massimo Ricotti for sharing star formation histories of the
simulated galaxies from Garcia et al. (2023) with us. AK was supported by the
National Science Foundation grants AST-1714658 and AST-1911111 and NASA ATP
grant 80NSSC20K0512.

We use simulations from the FIRE-2 public data release (Wetzel et al., 2023).
The FIRE-2 cosmological zoom-in simulations of galaxy formation are part of
the Feedback In Realistic Environments (FIRE) project, generated using the
Gizmo code (Hopkins, 2015) and the FIRE-2 physics model (Hopkins et al.,
2018). Analyses presented in this paper were greatly aided by the following
free software packages: NumPy (Oliphant, 2015), SciPy (Jones et al., 01 ),
Matplotlib (Hunter, 2007), FSPS (Conroy & Gunn, 2010a) and its Python bindings
package Python-FSPS111https://github.com/dfm/python-fsps, BPASS stellar
population synthesis tables for ionizing luminosity Byrne et al. (2022), and
Colossus cosmology package (Diemer, 2018). We have also used the Astrophysics
Data Service (ADS) and arXiv preprint repository extensively during this
project and the writing of the paper.

## Data Availability

The codes and data used in this study are available upon request.

## References

  * Adamo et al. (2020) Adamo A., et al., 2020, Space Sci. Rev., 216, 69
  * Adamo et al. (2024) Adamo A., et al., 2024, arXiv e-prints, p. arXiv:2401.03224
  * Asada et al. (2024) Asada Y., et al., 2024, MNRAS, 527, 11372
  * Atek et al. (2023) Atek H., et al., 2023, MNRAS, 524, 5486
  * Atek et al. (2024) Atek H., et al., 2024, Nature, 626, 975
  * Austin et al. (2023) Austin D., et al., 2023, ApJ, 952, L7
  * Belokurov & Kravtsov (2023) Belokurov V., Kravtsov A., 2023, MNRAS, 525, 4456
  * Bigiel et al. (2008) Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2846
  * Bouwens et al. (2023a) Bouwens R., Illingworth G., Oesch P., Stefanon M., Naidu R., van Leeuwen I., Magee D., 2023a, MNRAS, 523, 1009
  * Bouwens et al. (2023b) Bouwens R. J., et al., 2023b, MNRAS, 523, 1036
  * Boylan-Kolchin (2023) Boylan-Kolchin M., 2023, Nature Astronomy, 7, 731
  * Byrne et al. (2022) Byrne C. M., Stanway E. R., Eldridge J. J., McSwiney L., Townsend O. T., 2022, MNRAS, 512, 5329
  * Cameron et al. (2023) Cameron A. J., Katz H., Witten C., Saxena A., Laporte N., Bunker A. J., 2023, arXiv e-prints, p. arXiv:2311.02051
  * Caplar & Tacchella (2019) Caplar N., Tacchella S., 2019, MNRAS, 487, 3845
  * Carilli et al. (2018) Carilli C. L., Murphy E. J., Ferrara A., Dayal P., 2018, arXiv e-prints, p. arXiv:1810.07536
  * Carnall et al. (2023) Carnall A. C., et al., 2023, MNRAS, 518, L45
  * Casey et al. (2024) Casey C. M., et al., 2024, ApJ, 965, 98
  * Castellano et al. (2022) Castellano M., et al., 2022, ApJ, 938, L15
  * Castellano et al. (2024) Castellano M., et al., 2024, arXiv e-prints, p. arXiv:2403.10238
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chemerynska et al. (2023) Chemerynska I., et al., 2023, arXiv e-prints, p. arXiv:2312.05030
  * Chon et al. (2022) Chon S., Ono H., Omukai K., Schneider R., 2022, MNRAS, 514, 4639
  * Ciesla et al. (2023) Ciesla L., et al., 2023, arXiv e-prints, p. arXiv:2309.15720
  * Cole et al. (2023) Cole J. W., et al., 2023, arXiv e-prints, p. arXiv:2312.10152
  * Conroy & Gunn (2010a) Conroy C., Gunn J. E., 2010a, FSPS: Flexible Stellar Population Synthesis, Astrophysics Source Code Library, record ascl:1010.043
  * Conroy & Gunn (2010b) Conroy C., Gunn J. E., 2010b, ApJ, 712, 833
  * Conroy et al. (2009) Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486
  * Correa et al. (2015) Correa C. A., Wyithe J. S. B., Schaye J., Duffy A. R., 2015, MNRAS, 450, 1521
  * Curtis-Lake et al. (2023) Curtis-Lake E., et al., 2023, Nature Astronomy, 7, 622
  * D’Silva et al. (2023) D’Silva J. C. J., Driver S. P., Lagos C. D. P., Robotham A. S. G., Summers J., Windhorst R. A., 2023, ApJ, 959, L18
  * Dayal et al. (2022) Dayal P., et al., 2022, MNRAS, 512, 989
  * Dekel et al. (2023) Dekel A., Sarkar K. C., Birnboim Y., Mandelker N., Li Z., 2023, MNRAS, 523, 3201
  * Desprez et al. (2023) Desprez G., et al., 2023, arXiv e-prints, p. arXiv:2310.03063
  * Díaz-García & Knapen (2020) Díaz-García S., Knapen J. H., 2020, A&A, 635, A197
  * Diemer (2018) Diemer B., 2018, ApJS, 239, 35
  * Donnan et al. (2023) Donnan C. T., McLeod D. J., McLure R. J., Dunlop J. S., Carnall A. C., Cullen F., Magee D., 2023, MNRAS, 520, 4554
  * Dressler et al. (2024) Dressler A., et al., 2024, ApJ, 964, 150
  * Endsley et al. (2023) Endsley R., et al., 2023, arXiv e-prints, p. arXiv:2306.05295
  * Faucher-Giguère (2018) Faucher-Giguère C.-A., 2018, MNRAS, 473, 3717
  * Feldmann (2013) Feldmann R., 2013, MNRAS, 433, 1910
  * Ferrara et al. (2023) Ferrara A., Pallottini A., Dayal P., 2023, MNRAS, 522, 3986
  * Finkelstein et al. (2022) Finkelstein S. L., et al., 2022, ApJ, 940, L55
  * Finkelstein et al. (2023) Finkelstein S. L., et al., 2023, ApJ, 946, L13
  * Franco et al. (2023) Franco M., et al., 2023, arXiv e-prints, p. arXiv:2308.00751
  * Garcia et al. (2023) Garcia F. A. B., Ricotti M., Sugimura K., Park J., 2023, MNRAS, 522, 2495
  * Gieles et al. (2018) Gieles M., et al., 2018, MNRAS, 478, 2461
  * Gnedin & Draine (2014) Gnedin N. Y., Draine B. T., 2014, ApJ, 795, 37
  * Griffen et al. (2016) Griffen B. F., Ji A. P., Dooley G. A., Gómez F. A., Vogelsberger M., O’Shea B. W., Frebel A., 2016, ApJ, 818, 10
  * Hainline et al. (2024) Hainline K. N., et al., 2024, ApJ, 964, 71
  * Harikane et al. (2023) Harikane Y., et al., 2023, ApJS, 265, 5
  * Harikane et al. (2024) Harikane Y., Nakajima K., Ouchi M., Umeda H., Isobe Y., Ono Y., Xu Y., Zhang Y., 2024, ApJ, 960, 56
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Huang & Kauffmann (2014) Huang M.-L., Kauffmann G., 2014, MNRAS, 443, 1329
  * Hunt et al. (2015) Hunt L. K., et al., 2015, A&A, 583, A114
  * Hunter (2007) Hunter J. D., 2007, Computing In Science & Engineering, 9, 90
  * Inayoshi et al. (2022) Inayoshi K., Harikane Y., Inoue A. K., Li W., Ho L. C., 2022, ApJ, 938, L10
  * Iyer et al. (2020) Iyer K. G., et al., 2020, MNRAS, 498, 430
  * Jermyn et al. (2018) Jermyn A. S., Steinhardt C. L., Tout C. A., 2018, MNRAS, 480, 4265
  * Jones et al. (01 ) Jones E., Oliphant T., Peterson P., et al., 2001--, SciPy: Open source scientific tools for Python, http://www.scipy.org/
  * Kannan et al. (2023) Kannan R., et al., 2023, MNRAS, 524, 2594
  * Katz (2019) Katz H., 2019, in Latif M., Schleicher D., eds, , Formation of the First Black Holes. pp 125--143, doi:10.1142/9789813227958_0007
  * Katz et al. (2015) Katz H., Sijacki D., Haehnelt M. G., 2015, MNRAS, 451, 2352
  * Katz et al. (2023a) Katz H., et al., 2023a, The Open Journal of Astrophysics, 6, 44
  * Katz et al. (2023b) Katz H., et al., 2023b, MNRAS, 518, 270
  * Keller et al. (2023) Keller B. W., Munshi F., Trebitsch M., Tremmel M., 2023, ApJ, 943, L28
  * Kennicutt (1998) Kennicutt Robert C. J., 1998, ARA&A, 36, 189
  * Kennicutt & Evans (2012) Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531
  * Kewley et al. (2019) Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, ARA&A, 57, 511
  * Kravtsov & Manwadkar (2022) Kravtsov A., Manwadkar V., 2022, MNRAS, 514, 2667
  * Kravtsov & Wu (2023) Kravtsov A., Wu Z., 2023, MNRAS, 525, 325
  * Krumholz & Dekel (2012) Krumholz M. R., Dekel A., 2012, ApJ, 753, 16
  * Leitner & Kravtsov (2011) Leitner S. N., Kravtsov A. V., 2011, ApJ, 734, 48
  * Leung et al. (2023) Leung G. C. K., et al., 2023, ApJ, 954, L46
  * Lewis et al. (2023) Lewis J. S. W., Ocvirk P., Dubois Y., Aubert D., Chardin J., Gillet N., Thélie É., 2023, MNRAS, 519, 5987
  * Lilly et al. (2013) Lilly S. J., Carollo C. M., Pipino A., Renzini A., Peng Y., 2013, ApJ, 772, 119
  * Looser et al. (2023) Looser T. J., et al., 2023, arXiv e-prints, p. arXiv:2306.02470
  * Ma et al. (2018) Ma X., et al., 2018, MNRAS, 478, 1694
  * Ma et al. (2019) Ma X., et al., 2019, MNRAS, 487, 1844
  * Ma et al. (2020) Ma X., et al., 2020, MNRAS, 493, 4315
  * Manwadkar & Kravtsov (2022) Manwadkar V., Kravtsov A. V., 2022, MNRAS, 516, 3944
  * Mason et al. (2023) Mason C. A., Trenti M., Treu T., 2023, MNRAS, 521, 497
  * McCaffrey et al. (2023) McCaffrey J., Hardin S., Wise J. H., Regan J. A., 2023, The Open Journal of Astrophysics, 6, 47
  * McLeod et al. (2024) McLeod D. J., et al., 2024, MNRAS, 527, 5004
  * Mirocha & Furlanetto (2023) Mirocha J., Furlanetto S. R., 2023, MNRAS, 519, 843
  * Morishita & Stiavelli (2023) Morishita T., Stiavelli M., 2023, ApJ, 946, L35
  * Morishita et al. (2024) Morishita T., et al., 2024, ApJ, 963, 9
  * Mowla et al. (2024) Mowla L., et al., 2024, arXiv e-prints, p. arXiv:2402.08696
  * Muñoz et al. (2023) Muñoz J. B., Mirocha J., Furlanetto S., Sabti N., 2023, MNRAS, 526, L47
  * Naidu et al. (2022) Naidu R. P., et al., 2022, ApJ, 940, L14
  * Oliphant (2015) Oliphant T. E., 2015, Guide to NumPy, 2nd edn. CreateSpace Independent Publishing Platform, USA
  * Pallottini & Ferrara (2023) Pallottini A., Ferrara A., 2023, A&A, 677, L4
  * Pan & Kravtsov (2023) Pan Y., Kravtsov A., 2023, arXiv e-prints, p. arXiv:2310.08636
  * Pavesi et al. (2019) Pavesi R., Riechers D. A., Faisst A. L., Stacey G. J., Capak P. L., 2019, ApJ, 882, 168
  * Polzin et al. (2023) Polzin A., Kravtsov A. V., Semenov V. A., Gnedin N. Y., 2023, ApJ in press, p. arXiv:2310.10712
  * Qin et al. (2023) Qin Y., Balu S., Wyithe J. S. B., 2023, MNRAS, 526, 1324
  * Rasmussen Cueto et al. (2023) Rasmussen Cueto E., Hutter A., Dayal P., Gottlöber S., Heintz K. E., Mason C., Trebitsch M., Yepes G., 2023, arXiv e-prints, p. arXiv:2312.12109
  * Ren et al. (2019) Ren K., Trenti M., Mason C. A., 2019, ApJ, 878, 114
  * Riaz et al. (2020) Riaz R., Schleicher D. R. G., Vanaverbeke S., Klessen R. S., 2020, MNRAS, 494, 1647
  * Roberts-Borsani et al. (2023) Roberts-Borsani G., et al., 2023, Nature, 618, 480
  * Robertson (2022) Robertson B. E., 2022, ARA&A, 60, 121
  * Robertson et al. (2023a) Robertson B., et al., 2023a, arXiv e-prints, p. arXiv:2312.10033
  * Robertson et al. (2023b) Robertson B. E., et al., 2023b, Nature Astronomy, 7, 611
  * Rodríguez-Puebla et al. (2016) Rodríguez-Puebla A., Behroozi P., Primack J., Klypin A., Lee C., Hellinger D., 2016, MNRAS, 462, 893
  * Shen et al. (2023) Shen X., Vogelsberger M., Boylan-Kolchin M., Tacchella S., Kannan R., 2023, MNRAS, 525, 3254
  * Simmonds et al. (2024) Simmonds C., et al., 2024, MNRAS, 527, 6139
  * Steinhardt et al. (2022) Steinhardt C. L., et al., 2022, ApJ, 934, 22
  * Su et al. (2018) Su K.-Y., et al., 2018, MNRAS, 480, 1666
  * Sugimura et al. (2024) Sugimura K., Ricotti M., Park J., Garcia F. A. B., Yajima H., 2024, arXiv e-prints, p. arXiv:2403.04824
  * Sun et al. (2023) Sun G., Faucher-Giguère C.-A., Hayward C. C., Shen X., Wetzel A., Cochrane R. K., 2023, ApJ, 955, L35
  * Tacchella et al. (2020) Tacchella S., Forbes J. C., Caplar N., 2020, MNRAS, 497, 698
  * Tacconi et al. (2020) Tacconi L. J., Genzel R., Sternberg A., 2020, ARA&A, 58, 157
  * Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
  * Trinca et al. (2024) Trinca A., Schneider R., Valiante R., Graziani L., Ferrotti A., Omukai K., Chon S., 2024, MNRAS,
  * Vallini et al. (2024) Vallini L., et al., 2024, MNRAS, 527, 10
  * Ventura et al. (2024) Ventura E. M., Qin Y., Balu S., Wyithe J. S. B., 2024, MNRAS, 529, 628
  * Vincenzo et al. (2016) Vincenzo F., Matteucci F., Belfiore F., Maiolino R., 2016, MNRAS, 455, 4183
  * Wang et al. (2023a) Wang Y.-Y., Lei L., Yuan G.-W., Fan Y.-Z., 2023a, ApJ, 954, L48
  * Wang et al. (2023b) Wang B., et al., 2023b, ApJ, 957, L34
  * Wechsler et al. (2002) Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52
  * Wetzel et al. (2023) Wetzel A., et al., 2023, ApJS, 265, 44
  * Wilkins et al. (2023) Wilkins S. M., et al., 2023, MNRAS, 519, 3118
  * Yung et al. (2023) Yung L. Y. A., Somerville R. S., Nguyen T., Behroozi P., Modi C., Gardner J. P., 2023, arXiv e-prints, p. arXiv:2309.14408
  * Yung et al. (2024) Yung L. Y. A., Somerville R. S., Finkelstein S. L., Wilkins S. M., Gardner J. P., 2024, MNRAS, 527, 5929
  * Zavala et al. (2024) Zavala J. A., et al., 2024, arXiv e-prints, p. arXiv:2403.10491

##  Appendix A High-$z$ UV LF and variations of depletion time

Star formation depletion time is one of the parameters of galaxy formation
model used in this study. Figure 14 illustrates the effect of changing
depletion time on the model UV LF. It compares the $z=16$ UV LF in the models
with $\tau_{\rm sf}=2$ Gyr and $\tau_{\rm sf}=0.25$ Gyr without SFR
stochasticity (thick and thin dashed lines, respectively). Although the
amplitude of the LF is increased by an order of magnitude in the model with
$\tau_{\rm sf}=0.25$ Gyr, its amplitude is still more than two orders of
magnitude smaller than the corresponding models with SFR stochasticity at the
level required to match the observational estimate of UV LF at this redshift
(solid and dot-dashed lines). Thus, in our model LF at $z\gtrsim 13$ cannot be
matched by a realistic increase of star formation efficiency. At the same
time, the fiducial model that matches observational estimates of UV LF is
consistent with the inferred star formation rates of $z>10$ galaxies.

Figure 14: Model UV luminosity functions at redshift $z=16$. The dashed lines
show model of UV LF with the fiducial parameters without any SFR stochasticity
in the $\tau_{\rm sf}=2$ Gyr (thick dashed line) and $\tau_{\rm sf}=0.25$ Gyr
(thin dashed line). The thick dot-dashed and thin solid lines show results for
the same models but with SFR stochasticity added with the PSD parameters
$\alpha=2$, $\tau_{\rm break}=100$ Myr and $\sigma_{\Delta}=0.25$ and 0.22,
respectively (corresponding to $\sigma_{\rm M_{\rm UV}}=2.2$ and $1.9$, see
Table 1). The figure shows that the amplitude of UV LF can be increased with
increasing efficiency of star formation (smaller $\tau_{\rm sf}$) or with SFR
stochasticity. For more efficient star formation smaller amplitude of SFR
stochasticity is required to match UV LF.

##  Appendix B Tests of the halo mass accretion rate approximation

In our model calculations, we use halo tracks constructed by integrating a
simple expression for the halo mass accretion rate of equation 2. Here we
present a comparison of this expression with other calibrations of halo mass
accretion rate and demonstrate its accuracy.

Figure 15: Comparison of the halo mass evolution obtained by integration of
the approximation used in this study (eq. 2) and individual mass accretion
histories of the MW-sized halos from the Caterpillar simulation suite Griffen
et al. (2016). The dashed line shows exponential mass accretion rate
approximation of Correa et al. (2015, dashed line).  Figure 16: Comparison of
the halo mass accretion rates using approximation used in this study (eq. 2)
and approximations of Correa et al. (2015, dashed line), Rodríguez-Puebla et
al. (2016, dotted line) and Yung et al. (2024, dot-dashed line) for the mean
halo evolution tracks computed using integration of eq. 2 for halos with $z=0$
$M_{\rm 200c}$ mass of $10^{11}$, $10^{12}$, $10^{13}$, $10^{14}$ in panels
from left to right.

Figure 15 shows halo tracks obtained by integrating different approximations
of the halo mass accretion rate for halos of $M_{\rm 200c}=10^{12}\,M_{\odot}$
at $z=0$. Specifically, we use the mass accretion rate of equation 2 and
approximations of Rodríguez-Puebla et al. (2016) and the mass accretion rate
of Correa et al. (2015) corresponding to the exponential approximation to the
halo mass as a function of $z$ (Wechsler et al., 2002). The $M_{\rm 200c}(z)$
obtained with these approximations for $\dot{M}_{\rm 200c}$ is compared to the
halo mass tracks for halos in the Caterpillar suite of high-resolution
simulations of $M_{\rm 200c}\approx 10^{12}\,M_{\odot}$ halos (Griffen et al.,
2016). All approximations produce similar tracks at $z>10$, but the
approximation of the equation 2 provides a better description of the mass
evolution at higher $z$.

Figure 16 shows a direct comparison of the $\dot{M}_{\rm 200c}$ for halos of a
broad range of $z=0$ mass to the halo mass accretion rate approximation of
Correa et al. (2015), Rodríguez-Puebla et al. (2016), and the recent
calibration of Yung et al. (2023) focusing specifically on halo evolution at
$z>5$. The figure shows that our approximation is close to the approximations
of Rodríguez-Puebla et al. (2016) and Yung et al. (2024).

The differences of $\lesssim 30-40\%$ at $z<17$ exist but they have a small
effect on the UV LF prediction. Figure 17 shows a comparison of the UV LFs at
$z=12$ for our fiducial model and model in which the mass accretion rate was
lowered by 30% at all redshifts. The figure shows that the differences induced
by such modifications of the halo mass accretion rate have a negligible effect
on the UV LF prediction. Although this may seem surprising given the common
assumption that $L_{1500}\propto{\rm SFR}\propto\dot{M}_{\rm 200c}$ in the
calculations of UV LF for a halo population, in our model the star formation
rate is set by the amount of star-forming gas, which constitutes only a
fraction of the total gas and is not directly linked to the instantaneous halo
mass accretion rate. This indicates that residual uncertainties in the halo
mass accretion rate are not a significant source of uncertainty in galaxy
formation model predictions at these redshifts.

Figure 17: Comparison of the fiducial model prediction for the UV LF at $z=12$
with and without SFR stochasticity (long-dashed and solid line) with the model
calculations in which halo mass accretion rate was reduced by 30% at all
redshifts (short-dashed and dotted lines), which corresponds approximately to
the typical difference of our approximation for $\dot{M}_{\rm 200c}$ (eq. 2)
and approximations of Yung et al. (2024) and Rodríguez-Puebla et al. (2016) at
these $z$. The figure shows that differences of this level have a negligible
effect on the UV LF at $z\sim 12$.