How bursty star formation shapes spectral diversity at high redshifts: breaking degeneracies in early galaxy formation

Introduction

Among the most striking discoveries from the James Webb Space Telescope (JWST) are the numerous luminous galaxies at $z > 10$ revealed by early NIRCam campaigns (e.g. [20, 102, 1, 9, 38, 94, 117, 143]. In the past two years, the discovery and subsequent spectroscopic confirmation of galaxies at $z \gtrsim14$ (e.g. [18, 103]) has pushed the observational frontier to only $\lesssim300$ Myr after the Big Bang. The high abundance of $z > 10$ galaxies poses one of the most significant challenges to current models of galaxy formation. Standard models developed before JWST, tracing the evolution of the halo mass function, predicted a much steeper decline in the number density of luminous galaxies toward these early epochs

Corresponding author: Viola Gelli

Email: [email protected] (e.g., C. A. Mason et al. [91]; S. Tacchella et al. [134]; L. Y. A. Yung et al. [146]). Even new models specifically tuned to match the observed galaxy population at $z \sim10$ still struggle to reproduce the latest UVLF constraints at $z \sim14$ (e.g. A. Weibel et al. [141]; L. Whitler et al. [143]; I. Chemerynska et al. [23]).

The majority of explanations proposed to explain the UVLF fall into two classes: either galaxies are intrinsically brighter at fixed halo mass relative to pre-JWST expectations; or the scatter in UV luminosity increases, allowing galaxies hosted by lower-mass halos to temporarily populate the bright end of the UV luminosity function. The first case includes scenarios with increased star formation efficiency (SFE) due to inefficient feedback (A. Dekel et al. [35]; Z. Li et al. 2024; H.-Z. Chen et al. 2026) or efficient cloud scale star formation (e.g. S. H. Menon et al. [96]; R. S. Somerville et al. [125]), a top-heavy initial mass function (e.g. K. Inayoshi et al. [73]; E. R. Cueto et al. [31]; A. Trinca et al. [138]; A. Hutter et al. [72]), or reduced dust attenuation due to radiation pressure (e.g. A. Ferrara et al. [52]; F. Fiore et al. [54]) or changes in the grain size distribution (e.g. J. McKinney et al. [93]; D. Narayanan et al. [105]; L. Sommovigo et al. [126]). The second class invokes bursty star formation histories, where feedback driven SFR fluctuations can temporarily boost UV luminosities in the numerous low-mass halos at high redshift into the bright end of the UV luminosity function (e.g. C. A. Mason et al. [92]; A. Pallottini & A. Ferrara [109]; X. Shen et al. [120]; G. Sun et al. [132]; A. Kravtsov & V. Belokurov [81]; I. Nikolić et al. [108]). Recent models including halo mass dependent burstiness (e.g. V. Gelli et al. [57]; J. B. Muñoz et al. [101]), motivated by predictions of feedback becoming increasingly efficient in lower-mass systems with shallow potential wells (e.g. C.-A. Faucher-Giguere 2018; X. Ma et al. 2018; V. Gelli et al. [61]; S. R. Furlanetto & J. Mirocha [56]), showed that this can further boost the UV luminosity function, bringing models predictions into good agreement with JWST measurements up to $z \sim12$. However, reproducing the observed UVLF at $z \sim14$ is still a major challenge. Crucially, models with increased SFE and increased burstiness both enhance the number densities of galaxies at fixed $M_{\mathrm{UV}}$, making it impossible to distinguish between these scenarios with currently observed luminosity functions alone.

Spectroscopy offers a promising way to break this degeneracy. Different spectral features are sensitive to star formation on different timescales, providing direct constraints on the underlying star formation histories (SFH) and its burstiness. As the first JWST NIRSpec prism spectra of luminous $z > 9$ galaxies were taken, they revealed a remarkable diversity in their physical properties. Some UV-luminous systems exhibit extremely strong UV emission lines, potentially linked to their compact morphologies and extreme surface brightness densities (e.g., GN-z11, GHZ2, MoM-z14 A. J. Bunker et al. [13]; M. Castellano et al. [21]; R. P. Naidu et al. [103]). These galaxies are often interpreted as compact starbursts undergoing rapid recent increases in star formation, or possibly hosting active galactic nuclei (e.g. Y. Harikane et al. [65]; V. Kokorev et al. [79]). In contrast, other luminous galaxies show weak emission lines and more extended morphologies (e.g., Gz9p4, GS-z14 K. Boyett et al. [10]; S. Carniani et al. [18]; Z. Chen et al. [25]). While the origin of weak emission line systems has been challenging to interpret in low resolution prism spectra, recent studies with additional constraints from ultra-deep NIRSpec grating spectroscopy (e.g. Z. Chen et al. [25]) and MIRI follow-up to access rest-frame optical nebular emission lines and Balmer breaks at $z > 10$ (e.g. J. M. Helton et al. [68]; Y. Harikane et al. [66]) suggest that, at least some, UV-bright weak-emission line galaxies may be observed after a recent downturn phase in their star formation histories. However, while detailed studies of individual objects are critical for understanding the physical conditions in early galaxies, they provide only single snapshots of star formation histories. Progress in understanding the physical origin behind the $z \sim14$ UVLF excess now requires large statistical samples to characterize the full galaxy population properties.

In the low redshift universe, spectral diagnostics such as combination of UV and Balmer-lines have long been used to constrain burstiness and SFH in large spectroscopic samples (e.g. D. R. Weisz et al. [142]; Y. Guo et al. [62]; N. Emami et al. [42]). JWST has now made the same approach possible at $z \sim6$, where observations reveal large scatter in the star-forming main sequence (e.g. J. W. Cole et al. [29]; L. Ciesla et al. [28]; C. Simmonds et al. 2024; C. Di Cesare et al. [36]; A. Stoffers et al. [128]), broad distributions of emission-line equivalent widths at fixed UV luminosity (e.g. R. Endsley et al. [44]; R. Navarro-Carrera et al. [106]), and substantial diversity in inferred ionizing photon production efficiencies (e.g. G. Prieto-Lyon et al. [113]; C. Simmonds et al. [124]). These broad distributions are naturally expected if galaxies are observed at different phases of highly time-variable SFHs, including intense starbursts, downturn phases, and temporary quenching episodes following feedback-driven suppression of star formation (e.g. T. J. Looser et al. [86]; V. Strait et al. [129]; T. J. Looser et al. [87]; J. Trussler et al. [139]; W. M. Baker et al. [2]). We are now at a stage where these analyses can be extended to $z > 9$. Thanks to on-going spectroscopic efforts, particularly with NIRSpec prism (e.g. JADES, CEERS, RUBIES, CAPERS; D. J. Eisenstein et al. [41]; A. J. Bunker et al. [14]; S. L. Finkelstein et al. [53]; A. de Graaff et al. 2025; C. T. Donnan et al. [40]), JWST is now delivering increasingly large spectroscopic samples of bright galaxies at $z \gtrsim9$. These samples are enabling statistical studies and revealing a huge variety in galaxies’ spectral properties, including broad distributions of UV slopes and Balmer-line equivalent widths (EW, e.g. G. Roberts-Borsani et al. [116]; A. Saxena et al. [119]; M. W. Topping et al. [137]; M. Tang et al. [136]), suggesting highly stochastic SFHs in the early universe. These spectra reveal not only significant diversity, but also evidence of evolution of spectral properties with redshift at fixed UV luminosity. For example, M. Tang et al. [136] find that the fraction of galaxies with strong emission lines increases from $z \sim6$ to $z \sim10$. Understanding and interpreting these trends requires a framework that connects star formation histories to the observed statistical distributions of spectral diagnostics and predicts their evolution toward $z \sim14$.

In this work, we develop such a framework by combining bursty star formation histories modeling with spectral synthesis. We create a dark matter halo-based star formation history model, calibrated at $z \sim6$, and explore its evolution at $z \sim10 - 14$. We test whether different scenarios for high-$z$ star formation histories that reproduce the $z \sim14$ UV luminosity function, such as enhanced star formation efficiency and increased burstiness, predict distinct distributions of spectral diagnostics, providing a way to identify the physical origin of the unexpectedly abundant luminous galaxies at the highest redshifts.

The paper is structured as follows: in § 2 we describe the model for populations of galaxies with bursty star formation histories; in § 3 we calibrate and anchor the model based on observations at $z = 6$; in § 4 we explore the predicted evolution from $z = 6$ to 10; in § 5 we explore different physical scenarios at $z = 14$ and their impact on a number of observables distributions; in § 6 and § 7 we present a discussion of the results and summary of the main conclusions.

Throughout the paper we assume a flat $\Lambda$CDM cosmology with $\Omega_m = 0.3$, $\Omega_\Lambda= 0.7$ and $h = 0.7$.

Methods

We aim to build a framework to model bursty star formation histories and emission for large populations of high-redshift galaxies. In the following, we describe how we model the time evolution of star formation in individual halos, then how we compute their spectral energy distributions (SED) to predict observable properties, and finally how we construct galaxy populations across different redshifts.

Halo growth and baseline star formation rate

We first derive the mass-assembly history of each dark-matter halo. Starting from its mass at a given redshift, we reconstruct its growth history using an analytic halo-growth model. The halo growth follows the Einstein–de Sitter high-redshift approximation¹, as derived in e.g. A. Dekel et al. [35]:

We use this halo-growth history to define a baseline (smooth) star-formation history, onto which stochastic fluctuations will be added. We assume that the star formation rate traces the baryonic accretion rate onto the halo (e.g. C. A. Mason et al. [91]; S. Tacchella et al. [134]; X. Shen et al. [120]; L. Y. A. Yung et al. [147]):

where $f_b = \Omega_b/\Omega_m$ is the baryon fraction and $\epsilon_{\star}(M_h)$ the mass-dependent but redshift-independent instantaneous star formation efficiency. This efficiency controls how accreted baryons are converted into stars and is therefore a key parameter of the model. Following previous models (e.g. B. P. Moster et al. [99]; X. Shen et al. [120]; J. B. Muñoz et al. [101]), we parameterize the star formation efficiency as a double power law:

where we fix $\alpha= -0.6$, $\beta= 0.5$, $M_c = 10^{11.8}M_{\odot}$. These values correspond to the best-fitting $z \sim6$ star formation efficiency inferred by J. B. Muñoz et al. [101], also consistent with observational constraints by M. Shuntov et al. [122] and Y. Harikane et al. [64]. The normalization of the star formation efficiency $\epsilon_{\star,0}$ is a free parameter of our model.

Stochasticity

Both simulations and observations point to a picture where the star formation histories of galaxies undergo significant fluctuations over time (e.g. C.-A. Faucher-Giguère 2018; X. Ma et al. 2018; V. Gelli et al. [61]; K. G. Iyer et al. [75]; S. Tacchella et al. [135]; T. J. Looser et al. [86]; A. Bhagwat et al. [5]; L. Ciesla et al. [28]; R. Endsley et al. [46, 43]). These fluctuations arise from the complex interplay between gas accretion, gas depletion, and feedback from radiation, massive stars, and supernovae, which can heat and eject gas from the galaxy. This interplay leads to alternating phases of enhanced and suppressed star formation, giving rise to bursty star formation histories. Modeling this variability is essential for interpreting high-z observations, as such fluctuations can lead to substantial variations in the observed spectral properties of galaxies over time (e.g. G. Sun et al. [131]).

Simulations show that SFR are typically found to be distributed normally in log-SFR space (e.g. R. Feldmann [49]; N. Caplar & S. Tacchella [17]; K. G. Iyer et al. [75]; A. Pallottini & A. Ferrara [109]). To analytically describe such realistic bursty star formation histories, we introduce a stochastic component to the baseline smooth SFR of Equation (2). For each galaxy, we assume that the instantaneous SFR fluctuates around the median baseline value:

where $\Delta$ describes the SFR fluctuations in log space. We note that, by adopting this parametrization in log-space, the burstiness of SFH is responsible for driving higher mean SFR with respect to the baseline median value, which will lead to higher average luminosities whenever stochasticity is present. In other words, the $\epsilon_\star$ parameter determines the median SFR trend, while $\Delta$ determines the mean SFR and fluctuations around the median.

We describe the time-variability $\Delta$ in the frequency domain through a power spectral density (PSD). Following previous works (e.g., N. Caplar & S. Tacchella [17]; S. Tacchella et al. [135]; K. G. Iyer et al. [74]; A. Kravtsov & V. Belokurov [81]; G. Sun et al. [133]; E. Burnham et al. [15]), we adopt a broken power-law PSD form:

where $f$ is frequency, $\sigma$ is the amplitude, $\tau$ is a characteristic timescale over which the fluctuations become uncorrelated. The parameter $\alpha$ sets the high-frequency slope of the PSD, such that larger values of $\alpha$ suppress short-timescale fluctuations more strongly, resulting in smoother star formation histories. Different theoretical models and simulations consistently find $\alpha= 2$ (N. Caplar & S. Tacchella [17]; S. Tacchella et al. [135]), which we fix throughout this work.

The characteristic timescale $\tau$ is a free parameter in our model and it is expected to be tightly related to the physical processes regulating star formation in galaxies. For instance, different feedback mechanisms typically suppress star formation on different characteristic timescales, e.g. Active Galactic Nuclei and radiation pressure on $\lesssim10$ Myr as they act efficiently after their onset, supernovae on $\sim10 - 50$ Myr as they reflect the delay due to lifetimes of different stars (e.g. C. Leitherer et al. [84]; P. F. Hopkins et al. [69]; S. R. Furlanetto & J. Mirocha [56]), and mergers on $> 100$ Myr (e.g. M. Boylan-Kolchin et al. [11]; M. H. Hani et al. [63]).

The amplitude of the fluctuations $\sigma_{\log\mathrm{SFR}}$ is also expected to strongly depend on the galaxies’ physical properties. Several models and simulations suggest a significant halo mass dependence: low-mass halos, with their shallow gravitational potential wells, are very sensitive to feedback processes and easily experience strong variations in the SFR, whereas more massive systems are typically more stable against both internal and external perturbations (e.g. C.-A. Faucher-Giguère 2018; X. Ma et al. 2018; V. Gelli et al. [61]).

For this reason, we model the SFR fluctuations’ amplitude as a mass-dependent quantity in the form:

where the normalization $\sigma_{10}$ corresponds to the amplitude at $M_{\mathrm{h}} = 10^{10}M_\odot$ and is treated as a free parameter in our modeling. We fix the mass dependent slope to a negative value of $\gamma= -0.15$ to produce higher burstiness in lower mass halos, reflecting the increasing ability of more massive halos to retain gas against feedback processes. The value was set following constraints obtained by J. B. Muñoz et al. [101] when fitting $\gamma$ as a free parameter, and it is also consistent with the results of I. Nikolić et al. [108] and A. Basu et al. [3]. We set a minimum amplitude of $\sigma_{\log\mathrm{SFR}} = 0.3$ to account for the level of burstiness expected from stochastic dark matter accretion rates alone (K. Ren et al. [115]).

We note that, because the amplitude of SFR fluctuations depends on halo mass, the level of burstiness naturally evolves with time as the halo grows: while the baseline SFR increases in time with halo mass, the level of burstiness gradually decreases, such that an individual halo experiences larger SFR fluctuations at earlier times than at later times. While this effect is often neglected in similar models (e.g. K. G. Iyer et al. [75]; E. Burnham et al. [15]; J. B. Muñoz et al. [101]), which assume a time-independent fluctuation amplitude for a given halo, our model captures the evolution of burstiness within individual halos. As a result, galaxies experience systematically burstier star-formation histories at earlier times, potentially affecting observables sensitive to star formation on longer timescales, such as the UV luminosity and Balmer-break strength (see the next § 2.3).

Spectral modeling

As we aim to infer the impact of bursty star formation histories into high-$z$ observables, a key ingredient of our framework is the forward modeling of the spectral energy distributions (SED) of galaxies.

We generate SEDs using STARBURST99 (C. Leitherer et al. [84]) with PADOVA stellar evolutionary tracks (G. Bertelli et al. [4]), computing both stellar and nebular continuum emission for each galaxy based on its SFH. The nebular continuum and emission-line equivalent widths are computed assuming that ionizing photons are fully absorbed by the surrounding gas and adopting standard case-B recombination. The SFHs are sampled in logarithmic timesteps, with the final $\sim50$ Myr resolved at a very fine resolution ($\Delta t \le1$ Myr) to accurately capture the effects of short-term variability on emergent emission. We adopt a standard Kroupa IMF, a double power-law with slopes of 1.3 and 2.3 in the ranges [0.1, 0.5] M$_\odot$ and [0.5, 100] M$_\odot$ respectively (P. Kroupa [83]). We discuss the potential implications of varying this assumption in § 6. We fix the stellar metallicity to $Z = 0.1$ $Z_\odot$, appropriate for the high-redshift systems considered (e.g. [33, 104, 118]). While exploring a full metal-enrichment history would be interesting and will be addressed in future work, metallicity variations are expected to be sub-dominant to SFH variability on the distributions of spectral properties we explore in this paper (e.g. [121, 60, 37, 46]).

To estimate the UV magnitude $M_{\mathrm{UV}}$, defined as the rest-frame absolute AB magnitude at 1500 Å, we add dust attenuation following C. A. Mason et al. [91] and adopting G. R. Meurer et al. [97] attenuation law $A_{\mathrm{UV}} = 4.43 + 1.99\beta$, including the UV slope $\beta(M_{\mathrm{UV},z})$ empirically from observations by R. J. Bouwens et al. [8]. We do not consider dust attenuation at $z > 10$, as recent results suggest its effect to be negligible in very high-$z$ galaxies ([32, 98, 137]).

From the modeled SED, we use a set of key observables, commonly used to characterize star formation in high-$z$ galaxies, that are sensitive to different timescales: UV magnitude, Balmer line equivalent

widths, UV slope, and the strength of the Balmer break. To build intuition, in Figure 1 we show the evolution of these observables, following a single burst of star formation (of $M_\star= 10^8 M_\odot$) for 300 Myr.

The top panel shows the evolution of the intrinsic UV magnitude $M_{\mathrm{UV}}$, i.e. prior to dust attenuation. The UV emission peaks at $\sim2-5$ Myr due to the post-ZAMS evolution of massive ($\gtrsim15M_\odot$) stars, including the onset of Wolf–Rayet phases, and subsequently gradually declines as these stars end their lives. At later times ($>5$ Myr), the UV continuum still receives contribution from longer lived lower mass stars, causing a gradual decline of the emission. The UV magnitude requires $\sim50-100$ Myr to fade by $\sim3$ mag, thus remaining sensitive to star formation over relatively long timescales.

Hydrogen emission line EW (H$\alpha$ and H$\beta$, second panel), defined as the ratio between line luminosity and the underlying continuum emission, are powerful tracers of recent star formation activity. Hydrogen recombination lines are produced in HII regions ionized by short-lived massive OB stars, and trace the SFR over the past $\lesssim5$ Myr, as the ionizing photon output drops rapidly once these stars evolve off the main sequence. In contrast, the optical continuum is primarily sensitive to the accumulated stellar mass and longer-lived stellar populations. As a result, hydrogen line EWs broadly trace the specific SFR. Following a starburst, the EW remains high for the first $\sim5$ Myr, before declining rapidly as the massive stars die. For example, the H$\alpha$ EW decreases by a factor $\sim10$ between 3 and 10 Myr².

The UV continuum slope (third panel) is defined as $\beta_{\mathrm{UV}}$, where $f_\lambda\propto\lambda^{\beta_{\mathrm{UV}}}$. We measure $\beta_{\mathrm{UV}}$ by fitting the slope of the rest-frame continuum between 1600–2500 Å, following the standard definition (e.g. D. Calzetti et al. [16]). The UV slope is sensitive to the presence of young stellar populations, nebular continuum emission, and dust attenuation, which can redden the intrinsic spectrum. Here we show the evolution of the intrinsic $\beta_{\mathrm{UV}}$, i.e. without any dust attenuation. Following a starburst, galaxies exhibit very blue slopes ($\beta_{\mathrm{UV}} \lesssim-2.5$) during the first $\sim10$ Myr, when the emission is dominated by young massive stars and nebular continuum emission. Around $\sim5-10$ Myr, the nebular contribution, which makes the UV slope less steep than the stellar continuum, rapidly declines. As a result the overall UV slope becomes temporarily even bluer. At later times, as the stellar population ages and the short-lived massive O/B stars die, the UV continuum becomes dominated by longer-lived lower mass stars and the spectrum progressively reddens, reaching $\beta_{\mathrm{UV}} \sim-2$ after $\sim100$ Myr.

The fourth panel shows the Balmer break, defined as the flux ratio $F_\nu(4200\text{\AA})/F_\nu(3500\text{\AA})$ (e.g. C. Binggeli et al. [6]; A. Vikaeus et al. [140]). Following a starburst, the Balmer break first exhibits a bump around $\sim5$ Myr due to strong nebular continuum emission powered by massive stars. After $\sim10$ Myr, it increases steadily as the stellar population ages, exceeding unity after $\sim20-30$ Myr when A-type stars begin to dominate the spectrum. The Balmer break is therefore primarily driven by the relative contribution of older stellar populations compared to young, massive stars: a weak break indicates a spectrum dominated by recent star formation, while a stronger break reflects an increasing contribution from evolved stars.

While these trends provide useful intuition, bursty star formation histories complicate the interpretation of galaxy observables, since different tracers retain memory of star formation over distinct and overlapping timescales. To illustrate this, we present two representative examples of bursty SFH and their emission in Figure 2 for two galaxies with same $M_{\mathrm{UV}}=-19$. The SFHs of two $z=6$ galaxies are constructed using the method described above, with parameters that will be calibrated in § 3. The left panel shows their star formation histories over the last 60 Myr, highlighting that systems with identical UV luminosities can exhibit very different SFHs. The galaxy in red is in a declining phase, with its SFR decreasing from $\sim5$ M$_\odot$ yr$^{-1}$ to nearly zero over the past $\sim10$ Myr, while the blue galaxy is experiencing a recent burst of star formation. A commonly used way to characterize this burstiness is to quantify the ratio between the star formation rate averaged over the last 3 Myr and that averaged over the last 50 Myr, SFR$_{3}$/SFR$_{50}$ (see e.g. J. W. Cole et al. [29]; T. J. Looser et al. [86]; R. Endsley et al. [43]; V. Kokorev et al. [79]; M. Tang et al. [136]; K. E. Heintz et al. [67]; A. Stoffers et al. [128]). Values above unity indicate a recent increase in star formation relative to the longer-term average, while values below unity trace a downturn or quenching phase. In the two examples shown here, this ratio has values of 0.05 for the declining system and 20.5 for the bursting one.

The diversity of SFHs at fixed $M_{\mathrm{UV}}$ propagates directly into differences in spectral properties. The right panel of Figure 2 shows the corresponding spectra, highlighting variations in Balmer break, Balmer line EW and UV slope driven by the different recent star formation activity. The bursting galaxy shows stronger emission-line equivalent widths, a Balmer jump, and slightly bluer UV slope, all reflecting the dominance of young stellar populations. Conversely, the galaxy in a declining phase shows weaker line emission, a well developed Balmer break and slightly redder UV slope. Studying the distributions of these observables as a function of $M_{\mathrm{UV}}$ across large samples can therefore be key to constraining burstiness in high-redshift galaxies.

Generating galaxy populations

Having described how we model the evolution and emission of individual halos, the next step is to simulate a population of galaxies to predict global properties, scaling relations, and the distributions of observables across different redshifts.

We assume a standard $\Lambda$CDM cosmology with dark-matter halos distributed according to the halo mass function (HMF), for which we adopt the form of D. S. Reed et al. [114]³. At each redshift we draw $\gtrsim10^7$ halos from the HMF. For each halo, we generate a random SFH with the PSD method described above, and forward model the corresponding SED.

We impose a lower limit on the halo masses capable of forming stars, set by the atomic cooling threshold corresponding to virial temperatures of $T_{\mathrm{vir}} \sim10^4$ K (V. Bromm & N. Yoshida [12]). We note that it is crucial to simulate halos down to this limit to be mass complete at fixed UV luminosity⁴. In fact, due to stochasticity in their SFH, low-mass halos may experience upward fluctuations in star formation and contribute to the galaxy population at high luminosities.

Bursty star formation at $z = 6$

To investigate how bursty star formation shapes galaxy properties across cosmic time, we here first establish a fiducial model calibrated using the well-constrained galaxy population at $z = 6$ (§ 3.1). We then test its predictions against multiple observational diagnostics at the same redshift (§ 3.2).

Model calibration at $z = 6$

We aim to identify a set of parameters that best reproduce the observed UV luminosity function and the inferred stochasticity of star formation histories at $z = 6$. For this purpose, we use the UV luminosity function from R. J. Bouwens et al. [7] together with estimates of the $\mathrm{SFR}_{3}/\mathrm{SFR}_{50}$ ratio distributions from R. Endsley et al. [43], which provide direct constraints on star formation variability in $z = 6$ galaxies.

R. Endsley et al. [43] analyzed a sample of 368 galaxies at $z \sim6$ across a wide magnitude range ($-22 \lesssim M_{\mathrm{UV}} \lesssim-16$), using JWST/NIRCam photometry to infer recent SFH from H$\alpha$-to-UV luminosity ratios combined with SED fitting. They find that fainter galaxies have broader distributions of H$\alpha$-to-UV, hinting at increased burstiness in lower mass halos. We note that the nature of the population of weak emission line sources is not fully established: low nebular emission line EW can arise not only from phases of suppressed star formation but also from other factors, such as high ionizing photon escape fractions or strong differential dust extinction (e.g. see R. Endsley et al. [43]). However, SFH are likely to play a dominant role in suppressing EW (e.g., I. Shivhaei et al. 2018; R. Endsley et al. [45]; Z. Chen et al. 2026a; Y. Harikane et al. [66]. Throughout this paper, we thus assume that the observed emission-line EWs, and the inferred SFR ratios, primarily trace the SFH. We discuss this assumption further in § 6.

We explore a grid of models varying the burstiness and star formation efficiency parameters introduced in § 2, in the ranges $\sigma_{10} \in[0.3, 1.5]$, $\tau\in[0.1, 100]$ Myr, and $\epsilon_{\star,0} \in[0.01, 1]$. We find a fiducial model with $\sigma_{10} = 0.6$, $\tau= 50$ Myr, and $\epsilon_{\star,0} = 0.03$ provides the best match to the $z \sim6$ observational constraints. Notably, the selected fluctuation timescale is consistent with values found in simulations, where SFHs are modulated by stellar feedback processes such as supernovae (e.g. X. Ma et al. 2018; V. Gelli et al. [61, 60]; A. Pallottini & A. Ferrara [109]). The values are also consistent with recent results by J. B. Muñoz et al. [101], who used a similar parametrization of the time-variable SFH.

We show the fiducial mass-dependent model resulting from the calibration in Figure 3. The left panel shows the UVLF compared with data from R. J. Bouwens et al. [7]. To assess the impact of the assumed mass-dependent scatter, we compare the results from the fiducial mass-dependent model with a model with no mass dependence, with $\sigma_{\log\mathrm{SFR}} = 0.5$, shown with dashed curve. Both models reproduce the observed UVLFs reasonably well. Their shapes differ slightly due to the different halo mass dependence of burstiness: in the mass-dependent model, the larger scatter in lower-mass halos enhances the faint end of the UVLF more than the bright end, resulting in a slightly steeper luminosity function.

The inset shows the resulting UV scatter $\sigma_{\mathrm{UV}}$ as a function of halo mass, defined as the standard deviation of the $M_{\mathrm{UV}}$ distribution at fixed $M_{\mathrm{h}}$ and computed in 0.5 dex mass bins. The UV scatter reflects, but is not identical to, the assumed mass-dependent scatter in log SFR (Equation (6)), since the UV luminosity is sensitive to star formation variations over longer timescales. We find $\sigma_{\mathrm{UV}}(10^{10}M_{\odot}) \sim0.9$ mag for our mass-dependent model, consistent with results from J. B. Muñoz et al. [101].⁵ The mass-independent model instead yields an approximately flat relation with $\sigma_{\mathrm{UV}} \sim0.8$ mag for $M_{\mathrm{h}} \lesssim10^{11}M_{\odot}$. At high halo masses, we see a decline of $\sigma_{\mathrm{UV}}$ driven by dust attenuation, which preferentially suppresses UV-bright galaxies during starburst phases, reducing the amplitude of upward luminosity fluctuations and hence the UV scatter. In the absence of dust attenuation, for the mass independent model the $\sigma_{\mathrm{UV}}$} -- $M_{\mathrm{h}}$ relation remains flat over the full halo mass range.

The two panels on the right show the fractions of galaxies in recent upturns (left) and downturns (right) of star formation, defined as those with ratios of $\mathrm{SFR}_{3}/\mathrm{SFR}_{50} > 5$ and $< 0.2$ respectively, following the definition by R. Endsley et al. [43] to compare with their data, as a function of $M_{\rm UV}$. In the mass-dependent model, we find an upturn fraction of $\sim10 - 15\%$, with a mild increase toward fainter UV magnitudes, while the downturn fraction increases more strongly toward the faint end. These trends with $M_{\rm UV}$ are driven by mass-dependent burstiness: lower-mass halos, which on average host fainter galaxies, exhibit stronger stochasticity and are therefore more likely to experience both enhanced and suppressed periods of star formation.

We find our mass-dependent model reproduces the general trends inferred by R. Endsley et al. [43], with very good agreement for the fraction of upturning systems across the full magnitude range and for downturning systems down to $M_{\rm UV} > -18$. In the faintest magnitude bin ($M_{\rm UV} \sim-17$, shown in lighter shade) our model underproduces downturn systems relative to the fraction inferred by R. Endsley et al. [43]. With our current prescriptions, we find no parameter combination can simultaneously reproduce this faint bin together with the upturn and downturn fractions at brighter luminosities and the UVLF shape. This may indicate that, particularly in faint galaxies, which are more likely to reside in the lowest mass halos, star formation may undergo longer suppression phases or more rapid post-burst quenching than captured by our single-timescale stochastic model (e.g. V. Gelli et al. [60]). We also note, however, that a high ionizing photon escape fraction can also reduce Hα EW (used to infer short timescales SFR), and may therefore contribute to an enhanced inferred fraction of galaxies in downturn phases (R. Endsley et al. [46]). We discuss this further in § 6.

By contrast, the mass-independent model (shown in gray) does not describe the observations well. This model yields upturn and downturn fractions that are independent of UV magnitude, and fails to reproduce the increase in downturn fraction toward fainter galaxies. We also show, for comparison, a model with a lower level of stochasticity, $\sigma_{\log{\rm SFR}} = 0.3$ (dotted gray curve), comparable to the level of scatter expected from stochastic halo accretion alone (e.g. K. Ren et al. [115]). This model predicts upturn and downturn fractions below $\lesssim5\%$ across all magnitudes, substantially underestimating the observed values. This demonstrates that the stochasticity associated with halo accretion alone is insufficient to reproduce the observed data. A substantial level of burstiness in star formation histories and its mass dependence are required to reproduce the observed fractions of galaxies in recent upturns and downturns.

Testing $z=6$ spectral observables predictions

To further validate the $z = 6$ model, we compare its predictions against two observables which were not used in the calibration. We compare our predictions for the distribution of Hα equivalent widths (EW, § 3.2.1) and intrinsic UV slope, $\beta_{\rm UV}$ (§ 3.2.2), both as a function of UV magnitude, to recent JWST observations. In each case, we expect the predictions of the models as a function of UV magnitude to be most important for understanding the impact of mass-dependent stochasticity.

H equivalent width

Figure 4 shows the distributions of H$\alpha$ equivalent widths in two UV magnitude bins ($M_{\mathrm{UV}} = -20$ to $-17.4$). The top panel shows the observed H$\alpha$ EW distributions as derived by R. Endsley et al. [46] from deep NIRCam imaging of more than 250 galaxies from JADES. The middle and bottom panels show the distributions as predicted by our fiducial mass-dependent model and by the mass-independent model respectively. The observations (top) show two trends in the H$\alpha$ EW distributions from the UV-bright to faint bin: a slight shift toward lower EW and a broadening of the distribution.

Our mass-dependent model (middle) successfully reproduces both trends. The shift of the median to lower EW reflects the lower SFR of less massive halos (due to lower accretion rates), which on average host fainter galaxies. The broadening of EW distributions for fainter galaxies arises naturally in the presence of mass-dependent burstiness: lower-mass halos, hosting on average fainter galaxies, exhibit stronger stochasticity and star formation rates variations, hence giving rise to broader distributions of H$\alpha$ EW. In contrast, the mass-independent model (bottom panel) reproduces a shift to lower EW at fainter UV-magnitudes, but does not reproduce the broadening of EW distribution. Because the level of burstiness is the same for all halo masses in this model, the shape of the EW PDF remains approximately unchanged between the two magnitude bins. Any small differences in shape between the two $M_{\mathrm{UV}}$ bins are consistent with statistical fluctuations due to finite sample size.

UV slopes

In Figure 5 we explore the predicted intrinsic UV slopes, i.e. prior to dust attenuation, as a function of UV magnitude. By excluding dust, we isolate the contribution of star formation variability to the UV slope. We also compare to observations from JWST data from M. W. Topping et al. [137] and M. C. Jecmen et al. [76].

Both the mass-dependent (blue solid line) and mass-independent models (dashed gray line) predict an almost flat relation for $M_{\mathrm{UV}} \lesssim-15$, with brighter galaxies having slightly bluer intrinsic UV slopes due to their higher average accretion and star formation rates. Compared to recent JWST observations at $z \sim6$ M. W. Topping et al. [137]; M. C. Jecmen et al. [76], our models both produce bluer median UV slopes than observed for galaxies brighter than $M_{\mathrm{UV}} < -18$. At the bright-end ($M_{\mathrm{UV}} < -18$) the observed behavior is driven by stronger dust attenuation in more luminous galaxies.

In the absence of dust attenuation, our model predicts an intrinsically nearly flat $\beta_{\mathrm{UV}}$-$M_{\mathrm{UV}}$ relation over this magnitude range.

At the faint end ($M_{\mathrm{UV} } \gtrsim-15$), the mass-dependent model predicts a progressive reddening, with $\beta_{\mathrm{UV}}$ increasing more rapidly toward lower luminosities. This behavior arises naturally from mass-dependent burstiness: faint galaxies reside in more bursty low-mass halos, where an increasing contribution from downturn phases shifts the median toward redder intrinsic UV slopes. Such faint galaxies have long been out of reach observationally. Excitingly however, recent JWST observations of strongly lensed galaxies from the GLIMPSE program show evidence for a similar faint-end reddening, with a turnover at $M_{\mathrm{UV}} \gtrsim-17$ (M. C. Jecmen et al. [76]). A flattening between $M_{\mathrm{UV}} \sim-18$ and $-17$ is also reported from deep JADES observations by M. W. Topping et al. [137]; and G. Felicioni et al. [51] find a consistent trend across multiple lensed fields extending to $M_{\mathrm{UV}} \sim-19$. Although the observed $\beta_{\mathrm{UV}}$ values are systematically redder than our intrinsic predictions, as expected due to dust attenuation, the detection of an enhanced reddening toward faint magnitudes is qualitatively consistent with mass-dependent stochastic star formation histories.

In the mass-dependent model, the intrinsic $\beta_{\mathrm{UV}}$ distribution shows a significant dispersion around the median (shaded region shows 16–84 percentile range). The scatter increases toward fainter magnitudes, as a consequence of mass-dependent burstiness: lower-mass halos hosting fainter galaxies are more strongly stochastic and therefore exhibit a wider range of UV slopes. By contrast, the mass-independent model shows approximately constant scatter and no systematic broadening toward faint magnitudes, since the burst amplitude does not depend on halo mass. A large scatter in the $\beta_{\mathrm{UV}}$–$M_{\mathrm{UV}}$ relation is measured in $z \sim6$ observations (e.g. [116, 119, 137]), with $\beta_{\mathrm{UV}}$ values spanning $-2.8 \lesssim\beta_{\mathrm{UV}} \lesssim-1.5$. This range is much broader than our model predictions: likely due to additional sources of variations such as dust attenuation and high escape fractions. However, our model suggests that burstiness alone contributes for part of the observed scatter, producing variations of $\beta_{\mathrm{UV}}$ of about $\sim0.2$ at fixed $M_{\mathrm{UV}}$.

Having established that the mass-dependent model successfully reproduces many statistical properties of galaxies at $z \sim6$, we now test whether the same model naturally reproduces the observed evolution of burstiness toward $z = 10$.

Redshift z=6 to z=10 evolution

We now apply the $z = 6$-calibrated model parameters at $z = 10$ to explore the predicted redshift evolution and test whether the model reproduces JWST burstiness measurements at these epochs. We once again focus on the fraction of upturns ($\mathrm{SFR}_{3}/\mathrm{SFR}_{50} > 5$) and downturns ($\mathrm{SFR}_{3}/\mathrm{SFR}_{50} < 0.2$) as a function of $M_{\mathrm{UV}}$, here focusing on their evolution from 6 to 10.

Figure 6 shows the prediction for the fiducial model (top row) and for the mass-independent model (bottom row). As increasingly large JWST NIRSpec prism spectroscopy samples are now making it possible to infer star formation histories and estimate these fractions at high redshifts, we can compare with recent observational results. In particular, we compare our predictions with M. Tang et al. [136], who analyzed a sample of 61 $-22 < M_{\mathrm{UV}} < -17$ galaxies with prism spectroscopy across $z = 9$–14 and compared their properties against a sample of 401 galaxies at $6 < z < 9$.

The fiducial model at $z = 6$ shows a fraction of upturns (top left) between $\sim8$–14% across the magnitude range $-21 < M_{\mathrm{UV}} < -15$, increasing to $\sim20\%$ at $z = 10$. This redshift evolution is primarily driven by mass dependent burstiness: at fixed $M_{\mathrm{UV}}$, galaxies are hosted by lower-mass halos toward higher redshift, which are intrinsically more bursty and therefore more likely to exhibit high $\mathrm{SFR}_{3}/\mathrm{SFR}_{50}$. The increase in accretion rates at earlier times further contributes to this trend, driving more rapidly rising star formation histories at higher-$z$ and therefore also higher average $\mathrm{SFR}_{3}/\mathrm{SFR}_{50}$ ratios. This latter effect is stronger for brighter galaxies, which are hosted by more massive halos with higher accretion rates at fixed redshift, resulting in a flatter dependence of the upturn fraction with $M_{\mathrm{UV}}$.

Comparing with observations, M. Tang et al. [136] inferred star formation histories by fitting the observed spectra and used these to estimate the $\mathrm{SFR}_{3}/\mathrm{SFR}_{50}$ ratios, adopting the same definition as in our model. They find that the fraction of upturns increases from $\sim10\%$ to $\sim25\%$ from $z \sim6$ to $z \sim10$, in agreement with our model’s prediction. This agreement supports our picture in which galaxies at higher redshift are progressively burstier, as they are hosted by lower-mass, more stochastic halos, resulting in higher upturn fractions, an effect further amplified by the rising accretion rates at high-$z$.

In contrast, the mass independent model (bottom) predicts a much weaker increase in the upturn fraction from $z = 6$ to $z = 10$. In this model, the increase is driven purely by the rise in accretion rates, pushing the average $\mathrm{SFR}_{3}/\mathrm{SFR}_{50}$ ratios to higher values. However, this increase is not sufficient to match the observed upturn fraction at $z \sim10$.

The panels on the right show the fraction of galaxies in downturn phases, those that have experienced a recent suppression in star formation and therefore undergo periods of reduced UV luminosity. In the fiducial mass-dependent model (top right), the downturn fraction increases toward fainter luminosities, reflecting the mass dependence of burstiness. However, in contrast to upturn fractions, the downturn fraction shows little evolution from $z = 6$ to $z = 10$. This is because, at fixed $M_{\mathrm{UV}}$, galaxies in downturn phases are preferentially systems that have been scattered down in UV luminosity from more massive halos. At higher redshift, a fixed UV luminosity corresponds to lower-mass and more bursty halos, increasing the probability of suppression episodes. At the same time, higher accretion rates at higher redshift lead to higher average $\mathrm{SFR}_{3}/\mathrm{SFR}_{50}$ ratios, reducing the probability of having low $\mathrm{SFR}_{3}/\mathrm{SFR}_{50}$ systems. These two effects approximately compensate, leading to a nearly redshift-independent downturn fraction evolution. The observed downturn fractions inferred by M. Tang et al. [136] are broadly consistent with our predicted trends, though they may hint at a mild increase in downturns at earlier epochs compared to our model. Conversely, the mass-independent model (bottom right) does not reproduce the data. The downturn fraction is constant with $M_{\mathrm{UV}}$, and it decreases from $z = 6$ to 10 by a factor $\sim1.5$. This is because, without mass-dependence, the $\mathrm{SFR}_{3}/\mathrm{SFR}_{50}$ distribution shifts to higher values from $z = 6$ to 10, while retaining the same scatter, hence reducing the fraction of galaxies below the $\mathrm{SFR}_{3}/\mathrm{SFR}_{50} < 0.2$ threshold. Importantly, this shows the downturn fractions are mainly determined by the amplitude and mass-dependence of stochasticity, rather than accretion rates. As we will discuss in § 6, this makes the downturn fraction a particularly sensitive tracer of bursty star formation.

Overall, the results presented in this and the previous section demonstrate that: (i) a high level of burstiness is required to reproduce the observed fractions of upturns and downturns and the distributions of H$\alpha$ EW, and (ii) a mass-dependent burstiness is necessary to reproduce the observed trends with $M_{\mathrm{UV}}$ and redshift. Having validated the fiducial model with mass-dependent burstiness at $z = 6 - 10$, we now explore the model’s predictions for key observables at $z = 14$.

Star formation models at $z > 10$

The JWST discovery of luminous galaxies at $z = 14$ poses a significant challenge to models of early galaxy formation that struggle to reproduce the high observed abundances. Motivated by the success of our model at $z \sim6-10$, here, we test whether the burstiness inferred from our $z = 6$ calibration can account for the observed UVLF at $z = 14$, and investigate possible modifications to the bursty star formation histories required to reconcile the predictions with the JWST observations.

Applying the $z = 6$ calibrated fiducial model at higher redshifts, we find that it reproduces the observed UVLF remarkably well up to $z = 11$. At $z = 12$–$13$, the model continues to match the faint end but begins to underpredict the abundance of bright galaxies at $M_{\mathrm{UV}} \lesssim-19$ (Appendix § 9). By $z = 14$, this discrepancy becomes substantial (solid gray curve in the left panel of Figure 7). As also found in previous studies (V. Gelli et al. [57]; J. B. Muñoz et al. [101]), mass-dependent scatter enhances the abundance of UV-bright galaxies relative to models without stochasticity, such as the model of C. A. Mason et al. [91] (dashed line). This enhancement arises because lower-mass halos are more abundant relative to higher-mass halos at high redshift and also exhibit stronger burstiness, allowing them to contribute significantly to the bright end of the UVLF. However, the scatter calibrated at $z = 6$ remains insufficient to reproduce the observed abundance of UV-bright galaxies at $z = 14$ (e.g. S. L. Finkelstein et al. [53]; Y. Harikane et al. [65]; A. Weibel et al. [141]; L. Whitler et al. [143]; K. C. Kreilgaard et al. [82]; D. J. McLeod et al. [95]).

We therefore explore different scenarios designed to reproduce the UVLF at $z = 14$, inspired by solutions proposed in recent theoretical work. As discussed in the introduction, these can be broadly divided into two classes according to how they impact galaxies’ star formation histories: (i) increasing the median SFR at fixed halo mass (e.g. via higher star formation efficiency as shown by e.g. A. Dekel et al. [35]; R. S. Somerville et al. [125]), and (ii) increasing the scatter in luminosity at fixed halo mass through more bursty star formation (e.g. C. A. Mason et al. [92]; J. B. Muñoz et al. [100])⁶.

Within our framework, we model these two scenarios by varying the parameters for the normalization of the star formation efficiency $\epsilon_{\star,0}$ and for the bursts amplitude $\sigma_0$. Starting from the fiducial $z = 6$ calibrated model, we set each of the two parameters to reproduce the observed $z \sim14$ UVLF from JWST, exploring the two scenarios:

• Increased SFE⁷: the normalization of the efficiency $\epsilon_{\star,0}$ is increased by a factor $\times5$, while the stochasticity is the same as in the fiducial model $\sigma_{\mathrm{UV}}(M_{\mathrm{h}} = 10^{10}M_{\odot}) \approx1$ mag;

• Increased burstiness: we increase the amplitude of stochastic fluctuations to $\sigma_{10} = 0.95$, with $\sigma_{\mathrm{UV}}(M_{\mathrm{h}} = 10^{10}M_{\odot}) \approx1.4$ mag, while leaving the efficiency unchanged.

These models are illustrated in Figure 7. On the left we show the predicted UVLF at $z = 14$ for the increased-SFE and increased-burstiness cases in orange and purple, that are tuned to match JWST data from S. L. Finkelstein et al. [53]; Y. Harikane et al. [65]; L. Whitler et al. [143]; K. C. Kreilgaard et al. [82]. More recent constraints from A. Weibel et al. [141] and D. J. McLeod et al. [95], based on wide areas and pure parallel surveys, are also shown and, although they present lower number densities than the other studies, they remain above the prediction of the fiducial model calibrated at $z = 6$, suggesting that some modification to the star formation histories is still required to match the observed abundance of bright galaxies at $z \sim14$.

The middle panel shows the UV scatter as a function of the halo mass for the fiducial model at $z = 6$, $10$ and $14$, and for the two $z \sim14$ models. The fiducial model, shown in grey, predicts $\sigma_{\mathrm{UV}}$ increases for lower masses, going from $\sim0.5$ mag to $\sim1.2$ mag in $M_{\mathrm{h}} = 10^{12}M_{\odot}$ to $10^{8}M_{\odot}$ halos at $z = 6$–$10$. We note the model also predicts the amplitude of $\sigma_{\mathrm{UV}}$ increases slightly at fixed halo mass at higher redshifts. This weak redshift dependent evolution of the UV scatter arises because UV luminosity traces star formation over $\sim50$ Myr, a timescale over which halos grow in mass. A halo at a given mass at $z \sim14$ grows more over the past 50 Myr than a halo of the same mass at $z \sim6$. As a result, its UV light includes contributions from earlier phases when the halo was less massive than the same final-mass halo at $z = 6$, and therefore more bursty. This naturally leads to a higher UV scatter at higher-$z$, despite our prescription of redshift independent SFR burstiness at fixed halo mass.

The $\sigma_{\mathrm{UV}}$ for the $z = 14$ increased-SFE model is identical to that of the fiducial $z = 14$ model, as the UV scatter is not affected by the increase of the star formation efficiency. In contrast, the increased-burstiness model has similar mass dependence of the UV scatter, but a systematic shift to higher values, ranging from $\sim1.5$ at low masses to $\sim1$ at the high-mass end. This shows

pendix C for how our mass-dependent burstiness model affects the mean SFR. that, despite the fact that the $z \sim14$ UV luminosity functions predicted by the increased-SFE and increased-burstiness models are essentially indistinguishable, the scatter in UV luminosity at fixed $M_{\rm h}$ differs between the two models.

The difference between the two models is driven by the underlying star formation histories. In the right panel of Figure 7, we show examples of two typical SFHs for galaxies with $M_{\rm UV} \sim-21$, i.e. selected to be representative of systems hosted by halos near the median halo mass at that UV luminosity in each model.

In the increased-burstiness model, a typical $M_{\rm UV} \sim-21$ galaxy has a star formation history consisting of a very strong young burst in the last $\sim5$ Myr, with little to no star formation previous to that event. It is hosted in a relatively low mass halo with temporarily boosted UV luminosity driven by the high SFR at the time of observation. Conversely, in increased-SFE model the galaxy has instead lower SFR at the time of observation, but it has been forming stars from earlier on, the emission from slightly older stellar populations ($\sim20$ Myr) still contributing to the observed UV light.

Given the similarity of the predicted UVLFs, our goal is to establish which observables can best distinguish between the models. As shown in § 3, the distributions of key observables across large samples can differ systematically between models, providing a powerful way to constrain the burstiness of star formation histories and the physical origin of the UVLF excess. We therefore investigate the predicted evolution of distributions of key observables from $z = 6$ to $z = 14$, and as a function of $M_{\rm UV}$ in these models. In particular, we focus on: SFR$_3$/SFR$_{50}$ ratios (§ 5.1), nebular emission line equivalent widths (§ 5.2), Balmer breaks (§ 5.3), UV slopes (§ 5.4) and $z = 17$ UVLF (§ 5.5).

SFR ratios

In § 3 we introduced the SFR$_3$/SFR$_{50}$ ratio as a powerful diagnostic for studying the statistical properties of recent star formation histories of galaxies. We here analyze the predictions of the SFR$_3$/SFR$_{50}$ ratio at $z = 14$, comparing the increased-burstiness and the increased-SFE scenarios introduced in the previous section.

We classify galaxies according to their SFR$_3$/SFR$_{50}$ ratio as experiencing strong upturns ($> 10$), upturns ($> 5$) and downturns ($< 0.2$). Our aim here is to quantify how likely galaxies at different redshifts and luminosities are to be caught in phases of enhanced or suppressed star formation. To this end, each panel of Figure 8 shows the fraction of galaxies in strong upturns, upturns and downturns respectively as a function of UV magnitude. When interpreting these trends, it is important to keep in mind that, at fixed halo mass, galaxies undergoing an upturn have boosted UV luminosities, while galaxies in downturns appear fainter than average.

As shown in § 4, the fiducial model (here shown with grey curves) predicts a strong increase in the upturn fractions from $z = 6$ to $10$, consistent with recent JWST measurements [136], while the downturn fraction remains approximately constant over the same redshift range. We now examine how these trends continue toward $z = 14$ in the increased-burstiness and increased-SFE scenarios, shown with the purple and orange curves respectively. We note that the increased-SFE scenario has the same level of stochasticity as the fiducial model (see middle panel of Figure 7). Therefore, the evolution of the SFR ratios at $z = 14$ is identical in the two models and, for this reason, we do not show the fiducial $z = 14$ predictions in the following plots.

In both scenarios, the fraction of galaxies in strong-upturn (top) and upturn (middle) phases increases relative to $z = 10$, with a stronger rise toward brighter magnitudes. This evolution reflects the combined effect of mass-dependent stochasticity and higher accretion rates at higher redshift. At fixed $M_{\mathrm{UV}}$, galaxies at higher-$z$ reside in lower-mass, more bursty halos, making upturns more likely. In addition, halo accretion rates increase more rapidly at earlier epochs, producing more rapidly rising star formation histories and therefore higher average $\mathrm{SFR}_{3}/\mathrm{SFR}_{50}$ ratios, further increasing the fraction of galaxies with recent upturns. This implies that, at fixed luminosity, galaxies are more likely to be observed in upturn phases at higher-$z$ than at lower-$z$. The effect becomes stronger at higher luminosities, where bursty phases become increasingly important. Thus, the bright end of the UVLF is increasingly dominated by galaxies undergoing recent upturns in star formation at higher redshift.

Comparing the two models, the increased-burstiness scenario predicts a significantly larger fraction of galaxies in upturn phases than the increased-SFE model. In particular, $\sim30$-$40\%$ of galaxies brighter than $M_{\mathrm{UV}} \lesssim-19$ are in strong upturn phases at $z = 14$ in the increased-burstiness case, compared to only $\sim20\%$ in the increased-SFE scenario. The differences between the two models are less pronounced when considering the broader upturn definition ($\mathrm{SFR}_{3}/\mathrm{SFR}_{50} > 5$). Meanwhile the downturn fractions (bottom panel) show a much stronger separation between the two models. The increased-SFE model shows no evolution in the downturn fraction relative to the $z = 6$ and $10$ fiducial model, as the higher accretion rates at high redshift and the increased stochasticity of lower-mass halos have competing effects that approximately compensate, as discussed in § 4. In contrast, the increased-burstiness model predicts downturn fractions higher by a factor of $\sim2$ relative to $z = 6$ and $10$ at all magnitudes, reaching $\sim20$-$30\%$ for galaxies with $M_{\mathrm{UV}} \gtrsim-18$. This highlights the importance of detecting galaxies in suppressed phases to constrain burstiness, as their fraction as a function of $M_{\mathrm{UV}}$ directly traces the mass dependence of SFR fluctuations $\sigma_{\log\mathrm{SFR}}(M_h)$. At $z = 6$ there is already growing evidence for such downturn systems (e.g. J. Trussler et al. [139]; R. Endsley et al. [43]) and even quiescent phases (e.g. “lulling” or “mini-quenched” galaxies; T. J. Looser et al. [86]; V. Strait et al. [129]; T. J. Looser et al. [87]). These sources become increasingly challenging to detect toward higher redshift, as they are expected to be more common at faint UV magnitudes (see also V. Gelli et al. [58]). However, our models show that, in the case of strongly bursty scenarios, a significant fraction of downturn galaxies is expected even among bright systems at $z = 14$, accounting for $\sim20\%$ of galaxies at $M_{\mathrm{UV}} \sim-21$. Interestingly, bright galaxies with evidence for recent suppression of star formation have now been observed at $z > 9$ (S. Carniani et al. [19]; Z. Chen et al. [25]; Y. Harikane et al. [66]), as we will discuss in § 6.

Our prediction for the $z = 14$ evolution of downturn fractions provides a promising observational test to discriminate between models: an increase of downturn galaxies with redshift, at any fixed $M_{\mathrm{UV}}$, would indicate enhanced burstiness, while no evolution would disfavor an increase in the burstiness amplitude. In the following subsections we turn to more direct observables that are sensitive to different timescales for star formation, to establish empirical tracers that may discriminate between these models.

Nebular line equivalent widths

Nebular emission line equivalent widths provide a sensitive probe of short-timescale SFH variations: hydrogen recombination lines trace star formation over $\lesssim5$ Myr, while the underlying continuum reflects the accumulated stellar population. Their ratio therefore approximately traces the galaxy’s sSFR. H$\beta$ can be observed with NIRSpec up to $z \sim10$, while H$\gamma$ extends to $z \sim12$. At higher redshifts, H$\alpha$ becomes accessible with MIRI (e.g. T. Y.-Y. Hsiao et al. [70]; J. M. Helton et al. [68]; J. A. Zavala et al. [149]). All three lines equivalent widths show similar trends, as they trace the same star-formation timescale. Here, we focus on H$\beta$ EW to facilitate comparison with the available data from M. Tang et al. [136] at $z \sim6$-$10$.

The two panels on the left in Figure 9 show the probability distribution functions of H$\beta$ EW at fixed magnitude of $M_{\mathrm{UV}} \sim-21$. We begin by examining the gray distributions (identical in both panels), which illustrate the evolution of the fiducial model from $z = 6$ to $z = 10$. The median EW shifts to higher values toward earlier times, increasing from $\sim100$ Å at $z = 6$ to $\sim150$ Å at $z = 10$. This increase reflects the higher accretion rates at earlier epochs, which lead to higher typical SFR. At the same time, galaxies at fixed UV luminosity have lower stellar masses, resulting in higher sSFR and therefore larger EW. The same increase in halo accretion rates toward higher redshift also drives the evolution of the SFR$_3$/SFR$_{50}$ ratios shown in Figure 8, leading to larger recent-to-past SFR ratios.

Simultaneously, the widths of the distributions are set by the burstiness of star formation. The width of the EW distribution increases at higher redshift. This EW widening with increasing redshift arises naturally from the mass dependence of burstiness: at fixed $M_{\mathrm{UV}}$, galaxies reside in lower-mass halos toward higher-$z$, and these halos exhibit stronger SFR fluctuations producing broader EW distributions. As a result, even in the fiducial model at $z = 6$ and $10$, at fixed $M_{\mathrm{UV}}$ we predict both an enhanced high-EW tail (systems observed during recent starbursts) and an enhanced low-EW tail (systems caught in downturn phases) toward higher redshift, even among luminous $M_{\mathrm{UV}} \sim-21$ galaxies.

The top-left and bottom-left panels also show the predictions of H$\beta$ EW distributions for the two $z = 14$ star formation models, increased-SFE and increased-burstiness respectively. In both scenarios, the trends already visible between $z = 6$ and $z = 10$ continue toward earlier epochs, with the distributions shifting to higher median EWs (due to an increase in accretion rates) and becoming broader (due to mass-dependent burstiness).

However, the models differ in the extent of this broadening. In particular, the increased-burstiness scenario produces an extremely wide EW distribution, with both very high-EW ($> 500$ Å) and very low-EW ($< 20$ Å) galaxies. This is a direct consequence of the large amplitude and time-variability of the SFR. During starburst phases, the SFR can reach typical values higher than in the fiducial and increased-SFE models, leading to very high EW, up to 600 Å for $M_{\mathrm{UV}} \sim-21$ galaxies. Conversely, during downturn phases the SFR can drop to very low values. Given that we adopt variability on timescale $\sim50$ Myr, typical of SN feedback, the downturns can last longer than the 10 Myr timescale over which the EW of the preceding starburst declines (see Figure 1). As a result galaxies can spend extended periods in very low-SFR state after the starburst signature has faded, leading to extremely low EW.

To quantify the broadening of the EW distributions, in the right panels of Figure 9 we show the fraction of strong and weak-line galaxies, defined as H$\beta$ EW $> 240$ Å and $< 50$ Å respectively, as a function of $M_{\mathrm{UV}}$ and their redshift evolution. The top-right panel shows the fraction of strong-line galaxies. At fixed UV magnitude, the fraction of strong line emitters increases with redshift because galaxies have higher typical star formation rates at earlier times, shifting the EW distribution toward larger values. This evolution closely mirrors the increase in the fraction of galaxies undergoing recent star formation upturns (Figure 8, left). The predicted increase of high-EW galaxies with $z$ is consistent with observations by M. Tang et al. [136], who report an increasing fraction of strong emitters from $z=6$ to $z=10$. In particular, for galaxies with $-20<M_{\mathrm{UV}}<-19$, they find a fraction of $\sim10\%$ of sources with EW $>240$ Å at $z\sim6$, rising to $25\%$ at $z\sim10$, in broad agreement with our predicted evolution. The measurement based on the $6<z<9$ sample lies slightly above our $z=6$ prediction, as expected since the observed sample also includes galaxies at higher redshifts, where our model predicts somewhat higher fractions.

Both the increased-SFE and increased-burstiness models predict a further increase in the fraction of strong-line galaxies from $z=10$ to $z=14$. However, the increase is substantially larger in the increased-burstiness model, reaching $\sim55\%$ for galaxies with $M_{\mathrm{UV}}\sim-21$, compared to $\sim45\%$ in the increased-SFE model.

The bottom-right panel of Figure 9 shows the fractions of galaxies with low Hβ EW. At all redshifts, we predict a higher fraction of weak-line systems among fainter galaxies, as the faint end is more strongly populated by lower-mass galaxies undergoing downturns. From $z=6$ to $z=10$, the fraction of weak-line galaxies decreases by a factor of $\sim2$ at all magnitudes, because the median $H\beta$ EW is higher at earlier epochs. However, the behavior at $z = 14$ depends strongly on the assumed star formation model. In particular, the increased-burstiness scenario yields higher fraction than at $z \sim10$, whereas the increased-SFE model predicts a further decline. In other words, if burstiness increases toward earlier times, the fraction of low-EW galaxies is expected to increase from $z = 10$ to $14$.

Measuring the fractions of strong and weak line emitters and their redshift evolution would provide a powerful way to break degeneracies between models. In particular, an increasing population of weak line emitters at fixed $M_{\mathrm{UV}}$ with redshift is a distinctive signature of increasingly bursty star formation histories.

Balmer break

Having examined how bursty star formation histories shape the nebular emission line EWs, which trace star formation over the past $\lesssim5$ Myr, we now turn to the Balmer break. This provides a complementary test of our models by tracing the emergence of spectra dominated by older ($\gtrsim20$-30 Myr) stellar populations.

Figure 10 illustrates the distributions of Balmer break strengths at fixed magnitude $M_{\mathrm{UV}} \sim-21$. From redshift $z = 6$ to $z = 10$ the typical median Balmer break values become lower. This trend is expected, as galaxies at higher redshift are on average found at earlier stages of their assembly and therefore host younger stellar populations with fewer old stars.

Focusing on the different $z = 14$ star formation models in the two panels, the increased-SFE model shows no significant evolution in the shape of the Balmer break strength distribution with respect to $z = 10$, with only a shift to lower values due to the average younger stellar ages. The scenario with increased-burstiness produces a broader distribution with respect to $z = 10$, spanning both very low and very high Balmer break values, from 0.3 to 1.6. In particular, galaxies with strong Balmer jumps (Balmer break $< 0.5$) represent 12% of the increased-burstiness population, compared to only 4% in the increased-SFE case. Very strong ($> 1.3$) Balmer break galaxies correspond instead to galaxies dominated by older stellar populations ($\gtrsim30$ Myr, see Figure 1) from previous starburst phases. These are very rare in both scenarios at $z = 14$, but twice more common in the increased-burstiness model (2%) compared to the increased-SFE (1%). Detections of Balmer breaks at such high redshift would therefore provide a signature of strong burstiness in star formation histories at these epochs. Such observations are feasible thanks to the JWST/MIRI instrument, which can access the continuum redward of the Balmer break for $z > 10$

galaxies (e.g. [e.g. 68]. Recent examples of Balmer break detections are now reported (e.g. [e.g. 66]), providing direct evidence of evolved stellar populations that may have formed during earlier starburst episodes, highlighting the power of this diagnostic to test early star formation histories.

UV slope

We now focus on the UV continuum slope $\beta_{\mathrm{UV}}$, which responds to variations in the SFH on up to $\sim100$ Myr timescales. Following a burst, $\beta_{\mathrm{UV}}$ becomes redder as the burst fades, with variations of $\Delta\beta_{\mathrm{UV}} \sim0.5$ over $\sim100$ Myr (see Figure 1). As in § 3 we show the intrinsic UV slopes, before dust attenuation, as $z \gtrsim9$ galaxies appear to be relatively unattenuated (e.g. F. [32, 98, 137]. We focus in particular on the redshift evolution of the $\beta_{\mathrm{UV}}$–$M_{\mathrm{UV}}$ relation. Figure 11 shows the trend for $z = 6$, $10$ and $14$, using the same color coding of Figure 8 for the two $z = 14$ star formation models. The redshift evolution reveals that galaxies at all magnitudes become progressively bluer towards earlier epochs. The general trend of $\beta_{\mathrm{UV}}$ decreasing with $z$ is commonly observed (e.g. [119, 137]) and usually attributed to the presence of dust which, building up at later times and in more evolved galaxies, plays a great role in reddening the spectra at lower $z$. However, here we show that, even when not including dust, spectra are expected to be bluer at high-$z$. Galaxies with $M_{\mathrm{UV}} = -20$ go from median intrinsic UV slopes of $\sim-2.5$ at $z \sim6$ to $\sim-2.65$ at $z \sim14$. This evolution is primarily driven by the increasing halo accretion rates at early times, which lead to higher average SFR and systematically younger stellar populations.

At $z = 14$, the two star-formation models both predict very blue UV slopes on average, $\beta_{\mathrm{UV}} \sim-2.65$ for $M_{\mathrm{UV}} \sim-20$, reflecting the extremely young ages of the stellar populations. However, the increased-burstiness model tends to have broader distributions of possible $\beta_{\mathrm{UV}}$ values at all magnitudes, including a substantial population of faint galaxies (30% in the range $-15 < M_{\mathrm{UV}} < -10$) with intrinsic UV slope $\beta_{\mathrm{UV}} \gtrsim-2.4$, compared to 20% in the increased-SFE model. These

systems correspond to galaxies undergoing downturn phases in their star formation histories, with the last strong burst having occurred more than $\sim40$ Myr earlier. In practice, detecting such faint sources at $z > 10$ and robustly measuring their UV slopes is extremely challenging even for JWST, although recent observations are pushing the current limits, with reported examples of red sources at $z > 10$ (e.g. [40]). Differences between the two $z = 14$ models also arise at the bright end, where the increased-burstiness model predicts a larger fraction of very blue galaxies with $\beta_{\mathrm{UV}} < -2.7$ at $M_{\mathrm{UV}} < -20$: 9% compared to only 2% in the increased-SFE model. We recall that these predictions do not include dust attenuation, which would likely redden the galaxies, particularly the more evolved and typically brighter ones. Nevertheless, an excess of very blue galaxies at the bright end together with very red galaxies at the faint end could hint at enhanced burstiness at high redshift.

z>15 luminosity functions

At the time of writing, the highest-redshift spectroscopically confirmed galaxy lies at $z = 14.44$ ([103]), with $M_{\mathrm{UV}} = -20.2$. While several galaxy candidates have been reported at even higher redshifts ($15 < z < 20$; e.g. [22, 80, 112]), none have yet been spectroscopically confirmed.

The confirmation of galaxies at these extreme epochs, or non-detections providing more stringent upper limits on the UVLF, may provide additional constraints on star formation models. In particular, models that are calibrated to reproduce the UVLF at $z \sim14$ may begin to diverge significantly at higher redshifts, offering a way to break degeneracies.

To explore this, in Figure 12 we show the $z \sim17$ UVLFs predicted by the fiducial model calibrated at $z = 6$ (gray curve) and the two $z = 14$ models, i.e. increased-burstiness and increased-SFE. Both $z = 14$ models predict higher number densities at $z = 17$ than the fiducial model by a factor $> 10$. The increased-burstiness model predicts higher number densities compared to the increased-SFE scenario. This difference arises from the growing contribution of low-mass halos, which are both more abundant and more strongly affected by burstiness due to the mass-dependent stochasticity. As discussed in Section 5 and in [57], this effect becomes increasingly important at higher redshift: the combination of the steep halo mass function and $\sigma_{\log\mathrm{SFR}}(M_{\mathrm{h}})$ decreasing with $M_{\mathrm{h}}$ boosts the relative contribution of the numerous low-mass systems to the UVLF.

In the figure, we also show current photometric constraints on the $z \sim17$ UVLF. Results from P. G. Pérez-González et al. [112] and V. Kokorev et al. [80] show an excess of galaxy candidates relative to pre-JWST predictions at $z = 17$, suggesting that the over-abundance observed at $z = 14$ may persist at higher redshifts. While V. Kokorev et al. [80] data seem to tentatively favor the increased-burstiness model, the current uncertainties cannot robustly rule out either scenario. D. J. McLeod et al. [95] instead finds number densities that lie below constraints from other data, indicating a more rapid dearth of galaxies above $z > 15$, potentially disfavoring both the increased-burstiness and increased-SFE scenarios as currently calibrated. However, the inferred upper limits on the galaxy number density remain approximately one order of magnitude above the $z = 17$ fiducial prediction, so we cannot rule out that a modification to the star formation histories may still be required at $z = 17$.

At present, the observational picture at $z \gtrsim15$ remains highly uncertain, with different analyses reaching substantially different inferred abundances. This regime is based on photometric candidates and small-number statistics, with a limited number of independent deep fields, making the inferred abundances particularly sensitive to contamination and selection effects and cosmic variance. Larger samples and confirming the nature of the $z > 14.5$ candidates will therefore be crucial to more

reliably establish the $z > 14.5$ population and constrain star formation models.

Discussion

Our results show that halo mass-dependent burstiness in star formation rate is essential to reproduce the observed properties of the $z = 6$ galaxy population and their evolution to $z = 10$, in particular the fractions of galaxies undergoing star-formation upturns and downturns. At $z = 14$ we showed that scenarios with increased-burstiness and increased-SFE, both built to reproduce the $z = 14$ UVLF observed by JWST, predict distinct distributions of spectral observables, offering a path to interpret the spectral diversity observed in JWST data and break the UVLF degeneracy. We now explore the most promising observational avenues to constrain stochasticity using galaxies in star-formation downturn phases (§ 6.1) and discuss future prospects (§ 6.2).

Constraining stochasticity with star formation downturns

Our results show that the fraction of galaxies caught in star formation downturns provides one of the most promising probe of the stochasticity of star formation (§ 5.1,§ 5.2). In particular, the dependence of the downturn fraction on $M_{\mathrm{UV}}$ reflects the mass dependence of the stochasticity, while its overall amplitude traces the strength of the fluctuations. Measuring the downturn fraction over a range of $M_{\mathrm{UV}}$ and redshifts can therefore constrain how stochasticity and its mass dependence evolve with cosmic time. The fraction of galaxies in recent upturns is similarly sensitive to stochasticity, but its evolution with redshift is degenerate with the increasing average gas accretion rate. Downturns are less subject to this degeneracy and therefore provide a cleaner probe of stochasticity.

To illustrate how downturn phases, here defined as $\mathrm{SFR}_3/\mathrm{SFR}_{50} < 0.2$, directly correlate with burstiness, we realize a set of models at $z = 14$ with increasingly high amplitude of bursts (scaling star formation efficiency to match the observed UVLF). In Figure 13 we show the fraction of galaxies in downturns as a function of the UV magnitude for models spanning a range of stochasticity, quantified by the UV scatter at $M_h = 10^{10}M_\odot$. Because the scatter is mass-dependent, downturns are more common at fainter $M_{\mathrm{UV}}$ in every model. For comparison, a model with mass independent scatter and $\sigma_{\mathrm{UV}}(M_h = 10^{10}M_\odot) = 1$ is also shown, highlighting that in this case the fraction of downturns does not vary with $M_{\mathrm{UV}}$. Crucially, at fixed $M_{\mathrm{UV}} = -21$ the downturn fraction increases by a factor of $\sim4$ as

$\sigma_{\mathrm{UV}}(M_{\mathrm{h}}=10^{10}M_{\odot})$ increases from 1 to 1.6 mag. In the most extreme model ($\sigma_{\mathrm{UV}=1.6}$ mag), more than 35% of $M_{\mathrm{UV}}=-20$ galaxies are caught in a downturn at any given time. The downturn fraction at fixed $M_{\mathrm{UV}}$ is thus a direct, quantitative estimator of burstiness.

Measuring the downturn fraction observationally requires robustly recovering the SFH of high redshift galaxies. As shown in § 5.2 and Figure 9, weak recombination line EWs are a key tracer of downturns: the fraction of weak-EW systems increases by $\sim1.5\times$ from $z=10$ to $z=14$ in the increased-burstiness model, while it declines over the same interval in the increased-SFE model. Although the predicted change is relatively small, the opposite evolution provides a clear qualitative distinction between the two scenarios. However, low EWs are not a unique signature of suppressed star formation as they could also arise when high fractions of ionizing photons escape the galaxy (e.g. R. Endsley et al. [46]). Ongoing and upcoming JWST programs provide promising opportunities to break this degeneracy. Deep NIRSpec grating spectroscopy enables more accurate reconstruction of star formation histories through resolved low-ionization absorption lines that can provide priors on the ionizing photon escape fraction $f_{\mathrm{esc}}$, breaking degeneracies present in prism spectra (e.g. Z. Chen et al. [25]; C. Keerthi Vasan G. et al. [78]).

Recovering $\mathrm{SFR}_{3}/\mathrm{SFR}_{50}$ directly would be an even stronger diagnostic for constraining stochasticity and distinguishing $z>10$ models, but requires combining emission-line EWs with information on the UV continuum and Balmer break strengths. In this respect, MIRI observations can probe Balmer breaks and H$\alpha$ emission at $z>7$, providing complementary diagnostics for both recent and older star formation activity (e.g. T. Y.-Y. Hsiao et al. [71]; J. M. Helton et al. [68]; Y. Harikane et al. [66]).

Intriguingly, recent observations are beginning to provide examples of bright high-redshift systems consistent with galaxies in post-burst downturn phases. Z. Chen et al. [25] report a $z=9.3$ $M_{\mathrm{UV}}=-21.66$ source with a Balmer break and low H$\beta$ EW of $\sim25$ Å, and weak UV emission and stellar wind features, consistent with a galaxy in a downturn phase in the past 10–20 Myr. Similarly, a recent study by Y. Harikane et al. [66] reported the detection of a Balmer break ($>1.1$) in a $M_{\mathrm{UV}}\sim-20$ galaxy at $z=11$. The detected Balmer break, together with the lack of strong H$\alpha$ and [OIII] emission lines, implies a star formation history in which this galaxy experienced an intense burst of star formation approximately $\sim10$ Myr prior to the time of observation, followed by a phase of suppressed activity. Current samples with deep grating resolution spectroscopy or MIRI coverage that enable robust estimates of past star-formation are still small, but these detections show that downturn phases are already identifiable at $z>9$. Expanding these samples to determine the fraction of $z>10$ galaxies in downturn phases, as a function of $M_{\mathrm{UV}}$ and redshift, will be important for constraining the role of stochastic star formation at high redshifts. In particular, if increased burstiness is responsible for the high abundance of luminous galaxies at $z\sim14$, our model predicts that downturns become more common with increasing redshift at fixed $M_{\mathrm{UV}}$.

Future prospects

Our framework provides a flexible and predictive approach to exploring a wide range of bursty star formation models and the signatures they leave on SEDs. Here, we discuss several physical ingredients that can impact bursty star formation and consider which assumptions could be relaxed to further explore their role and identify additional observables.

Despite the success of our fiducial mass-dependent model at $z\sim6-10$, departures from observations point toward important additional physics shaping galaxies’ spectra. In particular, in § 3 we showed that the model underproduces the downturn fraction in the faintest magnitude bin reported by R. Endsley et al. [43] by a factor $\sim1.5\times$. The higher observed downturn fraction could point to suppressed phases lasting longer in low-mass galaxies than our model currently allows. In our model, we characterize SFH fluctuations with a single timescale of $\tau\sim50$ Myr, consistent with values found in high-redshift galaxy simulations. However, multiple physical processes operating over different timescales could in principle give rise to more complex SFHs characterized by multiple timescales (e.g. S. Tacchella et al. [135]). In addition, the duty cycle of starbursts and downturns may not be symmetric, particularly in low-mass galaxies where strong feedback episodes can easily evacuate gas and delay the recovery of star formation (e.g. S. R. Furlanetto & J. Mirocha [56]; T. Dome et al. [37]; V. Gelli et al. [59, 58]). Incorporating these refinements could resolve the faint-end discrepancy with R. Endsley et al. [43]. However, because this effect is again driven by a halo-mass dependence, it is not expected to alter the redshift evolution of the downturn fraction, which retains a similar slope across epochs (§ 5.1). This highlights the importance of constraining the luminosity dependence of downturn fractions, which can also provide key information on the characteristic timescales and duty cycles governing stochastic star formation.

The faint-end downturn discrepancy could also reflect a higher ionizing photon escape fraction in the faintest galaxies, which would suppress the inferred EW independently of the underlying SFH, however we expect this may be a sub-dominant effect to SFH. We do not explicitly model ionizing photon escape fractions or their temporal variability. Variations in $f_{\mathrm{esc}}$ could affect observables such as $\beta_{\mathrm{UV}}$ and nebular emission line equivalent widths, especially immediately following starburst episodes. In particular, high escape fractions can temporarily produce very blue UV slopes over short timescales ($< 5$ Myr; M. W. Topping et al. [137]), introducing additional degeneracies with SFH-driven signatures. However, this effect is expected to impact only a small fraction of galaxies: as shown by M. W. Topping et al. [137], very blue galaxies ($\beta_{\mathrm{UV}} < -2.8$), which require very high escape fractions to explain, account only for $\lesssim4\%$ of their JADES sample of galaxies at $z \gtrsim6$ with $M_{\mathrm{UV}} > -20$. Deep JWST spectroscopy providing detections of multiple low-ionization absorption features may independently place constraints on $f_{\mathrm{esc}}$ in high-redshift galaxies (e.g. Z. Chen et al. [25]; C. Keerthi Vasan G. et al. [78]), helping to disentangle these effects.

Our results show that increased burstiness could provide an explanation for the diverse spectral properties of galaxies at $z > 10$, raising the question of which physical processes can drive such variation in SFR and suggesting strong feedback processes at play. One interesting possibility could be variation of the stellar IMF, as a larger population of massive stars could, in addition to increasing the UV luminosity produced per unit stellar mass (e.g. E. R. Cueto et al. [31]; A. Trinca et al. [138]; A. Hutter et al. [72]), also enhance stellar feedback, and give rise to more bursty star formation histories (e.g. E. Cueto & A. Hutter [30]). Recent ultra-deep rest-UV spectroscopy of $z \gtrsim9$ galaxies, including GN-z11, has revealed stellar-wind features consistent with a contribution from very massive stars with masses $> 100\,M_{\odot}$ (Z. Chen et al. [26]; R. Marques-Chaves et al. [90]). From a theoretical perspective, the IMF may depend on the physical conditions of star formation, with studies suggesting that low-metallicity environments could favor a more top-heavy IMF (e.g. S. Chon et al. [27]). Our model currently assumes a Kroupa IMF with an upper mass limit of $100\,M_{\odot}$, and therefore does not directly capture these effects. Future work exploring IMF variations and the contribution of very massive stars and their feedback will be important to understand the interplay with bursty star formation and the resulting spectral signatures at $z > 10$.

Beyond the stellar emission and spectral features explored here, bursty star formation may also affect other observable properties. In particular, bursty star formation may impact gas-phase metallicities, as episodes of intense star formation may drive both rapid metal production and feedback-driven outflows. This can introduce scatter in the mass-metallicity relation (e.g. M. Curti et al. [33]; A. Pallottini et al. [110]). Although metallicity variations are expected to have a smaller impact on $M_{\mathrm{UV}}$ than variations in SFR, they could strongly affect metal emission-line strengths and ratios, providing additional observables with which to test star formation stochasticity. Our model assumes a fixed metallicity of $Z_{\star} = 0.1\,Z_{\odot}$, but in future work, exploring the interplay between metallicity and bursty star formation could provide additional constraints on the physical processes shaping early galaxies.

Overall, our results show that spectral diagnostics provide a powerful avenue to distinguish between different physical scenarios for early galaxy formation. Combining this approach with additional probes, such as metallicity-sensitive observables and galaxy clustering (see e.g., J. B. Muñoz et al. [100]; V. Gelli et al. [57]; I. Nikolić & C. A. Mason [107]; G. Sun et al. [130]), will provide complementary constraints on the physical processes shaping the first galaxies.

Conclusions

We developed a framework to model bursty star formation histories across cosmic time, with the goal of assessing the impact of stochasticity on multiple galaxy observables. We model stochastic fluctuations of the star formation histories of galaxies, characterized by the halo mass-dependent instantaneous star formation efficiency $\epsilon_\star(M_{\rm h})$ and the burst amplitude $\sigma_{\log{\rm SFR}}(M_{\rm h})$. We generate galaxy populations by sampling halos from the halo mass function at each redshift and assigning bursty star formation histories, which we then forward-model into spectral energy distributions to predict key observables ($\beta_{\rm UV}$, hydrogen emission-line EW, and Balmer break strengths). We calibrate the model at $z = 6$ to reproduce both the UV luminosity function and the observed ratios of ${\rm SFR}_{3}/{\rm SFR}_{50}$. Our main findings are:

• A significant level of burstiness is required at $z = 6$ to reproduce the distributions of ${\rm SFR}_{3}/{\rm SFR}_{50}$ inferred from UV-to-optical SED fitting by R. Endsley et al. [43], corresponding to a UV scatter of $\sigma_{\rm UV}(10^{10}\,M_\odot) = 0.9$ mag.

• The burstiness must be mass-dependent, with low-mass halos more bursty than high-mass ones, to reproduce the increased fractions of galaxies in downturn phases (${\rm SFR}_{3}/{\rm SFR}_{50} < 0.2$) toward fainter $M_{\rm UV}$ at $z = 6$ observed by R. Endsley et al. [43]. The same model also reproduces the broader distributions of H$\alpha$ EW toward fainter $M_{\rm UV}$ observed by R. Endsley et al. [46] and the evolution in upturn and downturn fractions from $z = 6$ to $z = 10$ inferred from NIRSpec spectra by M. Tang et al. [136], which a mass-independent model fails to reproduce. This mass dependence is consistent with that found by J. B. Muñoz et al. [101], and reflects the greater susceptibility of lower-mass halos to star formation suppression.

• When extrapolated from $z = 6$ to higher redshifts, at fixed $M_{\rm UV}$, the model predicts that the median observables evolve toward: higher H$\beta$ EW, bluer UV slopes, and weaker Balmer breaks, mainly driven by increasing halo accretion rates. Simultaneously, the distributions of ${\rm SFR}_{3}/{\rm SFR}_{50}$ and H$\beta$ EW broaden toward high redshift, as galaxies at fixed $M_{\rm UV}$ reside in lower-mass, more bursty halos.

• Our mass-dependent model predicts that the fraction of galaxies in upturns (${\rm SFR}_{3}/{\rm SFR}_{50} > 5$) at fixed $M_{\rm UV}$ naturally rises with redshift, with a larger increase at brighter magnitudes. This is in agreement with the increase from 10% to 25% between $z \sim6$ and 10 estimated by M. Tang et al. [136] from JWST spectroscopy. Thus, we predict the bright end of the UVLF becomes increasingly dominated by galaxies undergoing recent star formation upturns at higher redshift.

• The mass-dependent model predicts that the downturn fraction (${\rm SFR}_{3}/{\rm SFR}_{50} < 0.2$) rises toward fainter $M_{\rm UV}$, directly reflecting the mass dependence of burstiness, but remains nearly constant with redshift unless stochasticity increases.

When extrapolated beyond the $z = 6$ calibration, the mass-dependent model remains in very good agreement with the observed UVLF up to $z = 10$--11. However, it starts to fall short at $z \gtrsim12$ at the bright end, and by $z = 14$ it underpredicts the observed abundance of UV-luminous galaxies by $\sim1.5$ dex, demonstrating that the level of mass-dependent burstiness inferred from the $z = 6$ calibration is insufficient to explain the abundance of UV-luminous galaxies at the earliest epochs.

Motivated by this discrepancy between the observed JWST $z \sim14$ UVLF and the predictions from the $z = 6$-calibrated model, we explore two modifications to the star formation histories that reproduce the observed $z = 14$ UVLF: increased burstiness ($\sigma_{\rm UV}(M_{\rm h} = 10^{10}M_\odot) = 1.4$ mag) and increased star formation efficiency (a factor $\times5$ in the $\epsilon_{\star,0}$ parameter). While these two scenarios are fully degenerate in their predicted number densities, they produce distinct and testable signatures in the distributions of spectral observables:

• The abundance of galaxies undergoing a recent downturn of star formation strongly depends on the level of burstiness, with downturns becoming more common toward fainter $M_{\rm UV}$ as a consequence of the mass dependence of burstiness. In particular, in the fiducial and increased-SFE model, the fraction of systems with ${\rm SFR}_{3}/{\rm SFR}_{50} < 0.2$ remains approximately constant from $z = 6$ to $z = 14$ at all $M_{\rm UV}$, while it increases by a factor of $\sim2$ in the increased-burstiness scenario, as a consequence of the more common suppression phases in the presence of strong stochasticity;

• Stronger SFR variations during galaxy evolution cause the increased-burstiness model to produce significantly broader distributions of hydrogen emission line equivalent widths, UV slopes and

Balmer breaks compared to the increased-SFE model;

• The increased-burstiness model predicts an increasing population of low-equivalent-width galaxies due to more frequent strong downturn phases: the fraction of galaxies with H$\beta$ EW < 50 Å is expected to increase by a factor of $\sim1.5$ from $z \sim10$ to $z \sim14$ at $-22 < M_{\rm UV} < -18$, while it decreases in the increased-SFE scenario.

Overall, our results show how the striking spectral diversity observed among UV-luminous galaxies at $z > 10$ is a natural consequence of bursty star formation histories, and that its statistical characterization provides a powerful way to constrain models of early galaxy formation. In particular, we predict that the abundance of galaxies in downturn phases represents a key discriminant between competing scenarios. While individual post-burst systems have already been identified at $z > 9$ [136, 25, 66], building statistically significant samples and tracking their abundance as a function of redshift and UV magnitude will be key for distinguishing enhanced burstiness from increased star formation efficiency models at $z \sim14$. Deep NIRSpec spectroscopy, combined with MIRI follow-up to access Balmer breaks and H$\alpha$, makes this an achievable goal in the near future.

Acknowledgments

We thank Ryan Endsley for providing useful comments on a draft of the manuscript, and Robert Feldmann and Andrei Mesinger for useful discussions. VG acknowledges support by the Carlsberg Foundation (CF22-1322). CAM acknowledges support by the European Union ERC grant RISES (101163035), Carlsberg Foundation (CF22-1322), and VILLUM FONDEN (37459). JBM was supported by the National Science Foundation under Grants AST-2307354 and AST-2408637, and NASA through grant JWST-GO-03224. Views and opinions expressed are those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

UV magnitude distributions

The relation between the halo mass and UV luminosity, commonly described through the conditional distribution $P(M_{\rm UV}|M_h)$, is a key ingredient in galaxy formation models for connecting UV luminosity functions to the underlying dark matter halo population.

In § 3, we find a scatter of the relation of $\sigma_{\rm UV}(M_h = 10^{10}M_\odot) = 0.9$ mag in our fiducial model, slightly lower than the value inferred by V. Gelli et al. [57]. However, in V. Gelli et al. [57], the conditional distribution $P(M_{\rm UV}|M_h)$ is assumed to be Gaussian, as done in many empirical and analytical models, whereas in our current framework it is not imposed a priori and instead emerges self-consistently from the stochastic star formation histories, and can therefore deviate from a Gaussian. Since $P(M_{\rm UV}|M_h)$ is key in determining how galaxies are redistributed across the UV luminosity function, here we examine the shape of the distribution as predicted by our model.

Figure A1 shows the resulting distribution for halos with $M_h = 10^{10}M_\odot$, compared to the Gaussian approximation adopted in V. Gelli et al. [57]. The model distribution is clearly non-Gaussian, with a pronounced skewness toward higher luminosities. Thus, while $\sigma_{\rm UV} = 0.9$ mag provides a useful measure of its overall width, it does not fully characterize the effect of stochasticity on the UV luminosity function. This asymmetry arises naturally from the non-linear mapping between bursty star formation histories and UV emission. Following a burst, the UV luminosity declines gradually over tens of Myr (Figure 1), so galaxies spend an extended period in moderately UV-bright states. Conversely, the suppression of UV luminosity during downturn phases is limited because stars formed in previous episodes continue to contribute to the UV continuum. The observed UV emission therefore reflects the superposition of stellar populations with different masses and ages, producing an intrinsically asymmetric distribution. Such behavior is also found in high-resolution simulations that explicitly resolve bursty star formation (e.g. [77, 132] (X. Ma et al. 2018), indicating that it is a generic consequence of stochastic SFHs rather than a feature specific to our model.

The skewness of $P(M_{\rm UV}|M_h)$ therefore modifies the impact of scatter on the UV luminosity function relative to the standard Gaussian scatter typically assumed in analytical and empirical models. The relative contribution of up-scattered and down-scattered halos differs from the Gaussian case, affecting both the normalization and the shape of the UV luminosity function.

Redshift $z=10$-$13$ UVLF

In Section 3 we calibrated the fiducial model at $z = 6$ to reproduce the observed UVLF (R. J. Bouwens et al. [7]) and the fractions of galaxies in star formation upturns and downturns (R. Endsley et al. [43]). Here, we test how the model performs in reproducing UVLFs at higher redshift when we apply the same parameters, without further tuning.

Figure A2 shows the UVLFs predicted by the fiducial model, with the $z = 6$ calibrated parameters compared with photometric and spectroscopic JWST constraints (R. J. Bouwens et al. [9]; P. G. Pérez-González et al. [111]; C. T. Donnan et al. [39]; S. L. Finkelstein et al. [53]; B. Robertson et al. [117]; C. J. Willott et al. [144]; M. Franco et al. [55]; Y. Harikane et al. [65]; A. Weibel et al. [141]; L. Whitler et al. [143]; I. Chemerynska et al. [23]; K. C. Kreilgaard et al. [82]; D. J. McLeod et al. [95]).

At $z = 10$ the model agrees well with the data across the full probed magnitude range. At $z = 11$–$12$ the agreement remains good at the faint end, but the model starts underpredicting the data at the bright end, $M_{\mathrm{UV}} \lesssim-18$. The discrepancy between the fiducial model and the data grows with redshift. By $z = 13$ the faint end of the model still agrees with the constraints from I. Chemerynska et al. [23] but lies $\sim0.5$ dex below those of L. Whitler et al. [143], while the bright end falls short of the data by $\sim0.5$–$1$ dex, anticipating the larger discrepancy at $z = 14$ (Figure 7).

Star formation efficiency

Star formation efficiency is a central quantity linking galaxy growth to dark matter halo assembly, and is widely used to interpret both stellar mass buildup and UV luminosity functions. Our model describes star formation as stochastic fluctuations around a median efficiency. It is therefore useful to clarify how these fluctuations affect the effective star formation efficiency, particularly as their amplitude varies with halo mass. A common definition of star formation efficiency is SFE = SFR$/\left(f_b\dot{M}_h\right)$. However, in the presence of bursty star formation histories, this definition is affected by stochastic variability in the star formation rate and its interpretation becomes non-trivial. In this section, we therefore quantify how burstiness modifies both the effective SFE and its scaling with halo mass.

In our framework, the star formation efficiency parameter $\epsilon_\star$ sets the normalization of the median SFE at fixed halo mass, as galaxies undergo log-normal stochastic fluctuations around this median. Physically, this implies that star formation proceeds through intermittent bursts and downturns around a well-defined baseline efficiency. Because of the lognormal nature of the fluctuations, the arithmetic mean SFR is shifted above the median by a factor $e^{\sigma_{\ln\mathrm{SFR}}^2/2}$. This reflects the fact that high-SFR episodes contribute more to the total stellar mass formed than low-SFR phases suppress it. As a result, bursty star formation increases the cumulative stellar mass growth at fixed halo mass, leading to a higher effective star formation efficiency and UV luminosity.

The integrated efficiency $M_\star/\left(f_bM_h\right)$ traces this cumulative effect of star formation over cosmic time. In Figure C1, we compare this quantity for our fiducial mass-dependent model and for a mass-independent burstiness model at $z = 6$.

We find that mass-dependent burstiness leads to a significantly shallower SFE–$M_h$ relation, similar to what found in FIREbox simulations by R. Feldmann et al. [50]. Here, this arises because low-mass halos experience stronger stochastic fluctuations in log-SFR, which enhance their cumulative stellar mass growth through strong bursts more efficiently than in massive systems. In contrast, in the mass-independent model the burst-induced boost is uniform across halo masses, preserving a steeper SFE–$M_h$ relation.