Population Models for Star Formation Timescales in Early Galaxies: The First Step Towards Solving Outshining in Star Formation History Inference

The growth of galaxies is regulated by physical processes operating across different timescales, ranging from the creation and destruction of giant molecular clouds at below $\sim$ 10 Myr to quenching at over $\sim 10^{3}$ Myr (see C. F. McKee & E. C. Ostriker 2007 for a theoretical overview). Understanding the observed diversity of the galaxy populations—a cumulative effect of these intertwined individual processes—has long been a focus in the field of galaxy formation and evolution.

Theoretical studies based on cosmological simulations offer a promising approach to dissect the formation and evolution of galaxies, but face the challenge of the wide dynamic range required to describe the physics involved. While models agree on a core set of critical processes shaping galaxy properties, their phenomenological implementations (sub-grid physics) differ, leading to a dispersion in model predictions (R. S. Somerville & R. Davé, 2015; T. Naab & J. P. Ostriker, 2017). Observationally, the imprints of these individual processes are anticipated to shape the star formation histories (SFHs) of galaxies, with each process exerting its influence over distinct characteristic timescales (e.g., K. G. Iyer et al. 2020). This means that by constraining the variability of SFHs, it would be possible to disentangle the various mechanisms driving or suppressing star formation.

Recent breakthroughs facilitated by the James Webb Space Telescope (JWST) suggest a staggering diversity of inferred SFHs within galaxy populations at early times, reigniting interest in the study of its variability (e.g., T. J. Looser et al. 2025; C. Witten et al. 2025). The discovery of mini-quenched galaxies at early epochs (T. J. Looser et al., 2024; V. Strait et al., 2023) and the over-abundance of luminous high-redshift galaxies have challenged theoretical models (e.g., Y. Harikane et al. 2023; C. A. Mason et al. 2023). Large amplitude, short timescale fluctuations in star formation rates (SFRs) have been suggested in $z>1$ studies (A. van der Wel et al., 2011; H. Atek et al., 2011; M. V. Maseda et al., 2018). This phenomenon may be prominent in the gas-rich $z>4$ universe, and has been proposed to explain these exotic objects (A. Pallottini & A. Ferrara, 2023; X. Shen et al., 2023; G. Sun et al., 2023b; A. L. Faisst & T. Morishita, 2024), in addition to several other theoretical possibilities, including an evolving initial mass function (IMF; E. R. Cueto et al. 2024; P. van Dokkum & C. Conroy 2024; L. Y. A. Yung et al. 2024), a decrease in dust attenuation with increasing redshift (V. Mauerhofer & P. Dayal, 2023; A. Ferrara, 2024), a high efficiency of gas conversion into stars (A. Dekel et al., 2023; A. Renzini, 2023; J. Sipple & A. Lidz, 2024), and black hole contribution to the UV luminosity (Y. Ono et al., 2018; F. Pacucci & A. Loeb, 2022). This stochasticity is predicted by theoretical models of galaxy formation at high redshift as the feedback timescale approaches or becomes smaller than the dynamical time of the system, disrupting the otherwise simple correlation between smooth gas accretion and SFRs (G. S. Stinson et al., 2007; S. Tacchella et al., 2016; R. Feldmann et al., 2017; C.-A. Faucher-Giguère, 2018; L. Bassini et al., 2023; P. F. Hopkins et al., 2023a; T. Dome et al., 2024). Short-lived and frequent deviations from an “average” SFR would naturally explain both unexpected low-mass quiescent systems as post-starbursts and overly luminous galaxies as starbursts.

The standard methodology for constraining SFH variability is via H$\alpha$-to-UV flux ratios, as the two are expected to trace star formation on different timescales (D. R. Weisz et al., 2012; B. D. Johnson et al., 2013; M. Sparre et al., 2017; N. Caplar & S. Tacchella, 2019; A. L. Faisst et al., 2019; J. A. Flores Velázquez et al., 2021). H$\alpha$ emission arises from recombination in nebular clouds, primarily driven by O and B-type stars on timescales of $\lesssim$ 5 or 10 Myr. UV emission is generated by O and B-type stars with lifetimes of $\lesssim 300$ Myr (J. Kennicutt, 1998), although due to outshining of the older stellar light by recent star formation (e.g., C. Papovich et al. 2001), the UV emission in star-forming galaxies typically traces $30-70$ Myr timescales. Pre-JWST studies at $z<3$ have demonstrated the feasibility of this method (e.g., J. C. Lee et al. 2009; N. Emami et al. 2019), but it faces two main challenges: (i) a purely rising SFH would lead to maximal H$\alpha$ and UV fluxes similar to a bursty SFH (V. Mehta et al., 2023); (ii) translating the H$\alpha$-to-UV flux ratios to timescales is difficult as the timescale probed by the UV is a strong function of the intrinsic SFH, meaning that every object would need to be uniquely and explicitly modeled (e.g., B. D. Johnson et al. 2013; J. A. Flores Velázquez et al. 2021).

In principle, one can solve both these problems by performing individual spectral energy distribution (SED) fits instead. Indeed, with advancements in panchromatic SED fitting techniques, numerous studies attempt to estimate SFH variability by fitting photometric observations from deep JWST surveys (e.g., L. Ciesla et al. 2024), in addition to using the H$\alpha$-to-UV flux ratios (Y. Asada et al., 2024; R. Endsley et al., 2024). However, the challenge here is that the inference of SFHs from individual SED fitting of photometric data, at least in broad-band photometry, cannot distinguish between an “extremely bursty” and a “smooth” model. B. Wang et al. (2024a) has demonstrated that both models can produce statistically acceptable solutions, as the Bayes factor shows no preference for either model.

Fundamentally, the challenge is in solving outshining. To infer SFH from integrated light is to infer the past from the present, but with asymmetric information: a recent burst of star formation dramatically outshines the older, fainter stellar populations, even in the near-infrared where old stars are relatively brightest. This prevents an accurate inference of the full SFH, introducing biases in the stellar mass and SFR estimates. C. Papovich et al. (2001) analyzed the HST photometric data of a sample of Lyman-break galaxies at $2\lesssim z\lesssim 3.5$, and found that a hypothetical old stellar population could in principle contain up to $\approx 3-8$ times the stellar mass of the young stars that dominate the observed SED; B. Wang et al. (2024a) showed that highly stochastic SFHs constrained only by broadband photometry produce $\sim 0.8$ dex systematics in SFR averaged over a short timescale and $\sim 0.3$ dex systematics in average stellar age. Recent studies performed on cosmological simulations also hint at systematic biases in inferred stellar masses (e.g., D. Narayanan et al. 2024; R. K. Cochrane et al. 2025). Owing to the expected increasing prevalence of burstiness at $z\gtrsim 4$ as outlined previously, and paired with the fact that the SFH is degenerate with most key parameters measured from galaxy SEDs, outshining is perhaps the biggest challenge in interpreting the light from galaxies in the early universe.

In addition, galaxy brightness (and therefore detectability in a flux-limited survey) and burstiness are inherently interconnected, as a recent burst makes galaxies bright. It follows that the duty cycle of the SFH variation would be biased in a flux-limited sample; that is, galaxies without a burst at the low-mass end (over a range of $10-100$ in stellar mass according to a study done on the FIRE simulation; G. Sun et al. 2023a) are preferentially missed. This selection effect has broad implications for statistical studies of the galaxy population such as the mass assembly history and global SFH.

All of the above factors motivate the need for a thorough examination of the extent to which burstiness can be observationally constrained. In this paper, we take a simple model of burstiness—essentially, repeated up-and-down fluctuations in the SFH—and evaluate the accuracy of recovered SFHs for bursty systems inferred with various methods, including state-of-the-art SED fitting of simulated high signal-to-noise (S/N) spectra, and population-level analyses of spectral features. The latter is motivated by recent developments in making the population distributions of the parameters of interest, such as photometric redshift or stellar population parameters, as the inference objective, bypassing individual fits (J. Alsing et al., 2024; J. Li et al., 2024). All tests have a quantitative target for recovery, enabled by our parameterization of the SFH. This work provides a comprehensive assessment of burstiness constraints, and critically, sets out a potential path forward for solving outshining.

The structure of this paper is as follows. Section 2 presents the SFH model and mock observations. Section 3 focuses on inferring SFH variations from individual galaxies, using common setups in SED fitting. Section 4 proceeds to show the improvement on these constraints given a correct SFH prior. Section 5 assesses the performance of the widely used population-level burstiness indicator, the H$\alpha$-to-UV flux ratio. Section 6 proposes our population-level approach to complement the H$\alpha$-to-UV method, and demonstrates its additional constraining power. Section 7 ties all the pieces together to address outshining, and also provides additional discussions on the individual SED fits including the prior influences, and the observability and selection effects due to SFH variations. We conclude in Section 8.

Where applicable, we adopt the best-fit cosmological parameters from the WMAP 9 yr results: $H_{0}=69.32$ ${\rm km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{M}=0.2865$, and $\Omega_{\Lambda}=0.7135$ (G. Hinshaw et al., 2013), and a G. Chabrier (2003) IMF over the mass range of $0.08-120~{\rm M}_{\odot}$. Unless otherwise mentioned, we report the median of the posterior, and 1$\sigma$ error bars are the 16th and 84th percentiles.

From Figure 3, it appears that different SFHs can have noticeable impacts on the resulting spectra. The most straightforward approach to estimating SFH variations then is via SED fitting of individual galaxies. We thus begin our analysis by joint photometric and spectral fitting of simulated galaxies, using a setup that follows the most common practices. For ease of comparison, this section contains the methods and results of this individual approach, with further discussions supplied in Section 7.3.

We note that the goal is to infer the population-level SFH parameters, as opposed to the usual integrated stellar population properties. For this reason, we conduct this experiment in three stages, described as follows.

Thus far, we have adopted the common setup (the continuity and bursty continuity priors) used in SED fits, but these state-of-the-art priors do not specifically weight towards the oscillating behavior of the simple SFHs; i.e., the model allows for but does not prefer the oscillating SFHs used as inputs. Having demonstrated that even high-quality spectra with a neutral prior cannot recover the detailed SFH, it is natural to ask whether individual SFHs can be accurately inferred with a prior tuned to the underlying population SFH. Such priors may be needed to solve the outshining problem, and can be empirically measured: B. Wang et al. (2023b) demonstrated the utility of physically motivated, informative priors in breaking degeneracies among stellar population parameters. More recently, Z. Gao et al. (2025) recovered cosmic SFR density and stellar mass density by constructing a prior based on a correlation between photometric colors and SFHs from the IllustrisTNG simulation (R. Weinberger et al., 2017; A. Pillepich et al., 2018), whereas J. T. Wan et al. (2024) showed the potential for a hierarchical modeling approach with a flexible stochastic SFH prior that aims to describe a wide range of SFH behaviors.

In this section, we test the hypothesis that, if the short-term SFH variability is known a priori, the SFH of individual galaxies can indeed be recovered via SED fitting. To this end, we fit mock observations using an oscillating SFH model as parameterized in Section 2.1. While this choice is intentionally avoided in Section 3.1—given that real galaxies are not expected to strictly follow a simple SFH form—it serves as a clear test here to address the question: if the model is informed that the SFH follows an oscillating pattern, can we recover the SFH from fitting individual spectra? Here, we use uniform priors on all four population-level SFH parameters: $\phi\in[0,2\pi]$, $\sigma\in[0,1.0]$ dex, $\alpha\in[-1,0]$, and $\delta t\in[10,70]$ Myr.

The summary statistics quantifying the stellar population parameter recovery are included in Figure 6. With the correct model for population-level burstiness, the distribution of the residuals in the inferred mass deviates from 0 by $<0.05$ dex, the scatter is $0.1$ dex, and the distribution of uncertainty-normalized residuals becomes close to a unit Gaussian with an outlier fraction of $0.2$.

Furthermore, the intricate variations in the SFH can also be recovered. An example is shown in Figure 10, where the inferred SFH is consistent with the true SFH. We note that the good recovery is contingent upon an accurate measurement of dust attenuation. This is because dust attenuation is needed to decode the Balmer emission lines, tracing SFR in the most recent $\sim 10$ Myr, and UV fluxes tracing SFR in the most recent $\sim 30-70$ Myr scale. The degeneracy between dust attenuation and intrinsic emission-line luminosity, given an observed flux, can cause a factor of $>2$ (up to orders of magnitudes) deviation from the true recent SFRs. Uncertainty in the deblending of H$\alpha$ and [N ii] introduces a second order effect where the recent SFH can be incorrectly inferred at a level of $\sim 20$%, since [N ii] typically is $\sim 10$% of H$\alpha$ flux.

Figure 10-b additionally includes two examples of inferring the SFHs with photometry alone. Generally, for cases where clear age indicators (e.g., a Balmer break) are captured by the medium bands, the input SFH can be recovered without spectroscopy. However, as alluded to earlier, in some cases accurate and precise measurements of the Balmer emission lines are necessary for recovering the SFH.

The above findings suggest that given the correct SFH model, it is indeed possible to infer the SFH by performing SED fitting. Put in another way, outshining can be solved if the right SFH model is known. A key assumption here is that the current burstiness observable in the population is consistent with the past SFH behavior. For cases where star formation is primarily self-regulated and stationary—rather than driven by external events like mergers, or evolving sharply with time due to e.g., gas availability—this assumption is likely reasonable. We revisit this point in Section 7.

To summarize, the key to address outshining is to empirically measure the population distribution of bursty SFHs. Possible ways forward to achieve this is the focus of the second half of the paper.

The most well-studied measurement of burstiness is likely the scatter in H$\alpha$-to-UV flux ratios. There is rich literature on this variability indicator (K. Glazebrook et al., 1999; J. C. Lee et al., 2009; D. R. Weisz et al., 2012; Y. Guo et al., 2016; N. Emami et al., 2019; T. Nanayakkara et al., 2020), and it can be considered as a prototype of a population-level fit. We therefore start the presentation of population-level constraints of burstiness with H$\alpha$-to-UV flux ratios, and test its constraining power using our mock observations.

In Figure 11, we first show the H$\alpha$-to-UV flux ratios as functions of the observed phases for two SFH families that differ in fluctuation amplitudes. For these cases, the H$\alpha$-to-UV flux ratios differ markedly. The larger amplitude leads to an increased scatter in H$\alpha$/UV, as expected. However, this difference becomes significantly less noticeable for models with the same fluctuation amplitude, but different durations and/or slopes.

The significant changes in the sSFRs for the more steeply rising SFH models drive the sudden change in H$\alpha$/UV near $\phi=\pi$. However, some scatter is present in this transition, making it appear smoother than one might expect from the abrupt changes in the simulated SFHs at $\phi=\pi$. The varying metallicities partly contribute to the smoothness. More importantly, the abrupt changes in the simulated SFHs are on the instantaneous sSFRs. As shown in J. Choi et al. (2017), the number of ionization photons—and therefore—the H$\alpha$-to-UV flux ratios, is proportional to SFR averaged over 10 Myr, as opposed to instantaneous SFRs.

The simulated H$\alpha$/UV ratios are broadly consistent with the observationally inferred values reported in Y. Asada et al. (2024), although the maximum value of $\sim 0.07$ lies at the extreme end of their distribution. This upper limit also exceeds the maximum value found in R. Endsley et al. (2024), likely due to differences in modeling assumptions, as R. Endsley et al. (2024) derived their H$\alpha$/UV ratios from BEAGLE (J. Chevallard & S. Charlot, 2016) fits to the photometry. We stress that our primary goal is not to precisely match the observations but rather to encompass a reasonable range of burstiness, since these simulations are designed to simply serve as a controlled baseline for evaluating the constraining power of different methods.

To quantify the constraining power of H$\alpha$-to-UV flux ratios, we compute the Wasserstein distance between the observed, noisy distributions to the model distributions. Also known as the earth mover’s distance, the Wasserstein distance is a measure of the similarity between two probability distributions. It is similar to the perhaps more widely known Kullback-Leibler divergence, but is more robust to outliers and is sensitive to the spatial arrangement of the probability mass (V. M. Panaretos & Y. Zemel, 2019), and has been applied to astronomical data (e.g., J. Li et al. 2024). It is possible to explore alternative distance metrics in the future, but here the best-fit model is found by simply identifying the shortest Wasserstein distance.

We note that while the observed, noisy distributions are from the same model grid used to test the individual SED fits, the model distributions being compared to are generated from a denser grid ($\sigma\in$ [0.2, 0.9] dex with 0.1 dex increments, $\delta t\in$ [15, 45] Myr with 5 Myr increments, $\alpha\in$ [-0.3, -0.7] dex with 0.05 increments). We repeat the mock observations with different realized noises and the resulting inference procedure 1,000 times for each SFH family sampled on the denser model grid as a robustness test. We then calculate the fraction of population-level SFH parameters correctly identified as the best-fit, and consider the modal best-fit parameters as the overall best-fit. For example, suppose we have a noisy distribution of H$\alpha$/UV corresponding to the SFH family of $\sigma=0.3,\delta t=40,\alpha=-0.4$. In the 1000 trials, 80% of the time the shortest Wasserstein distance gives $\sigma=0.4$, 10% of the time gives $\sigma=0.3$, and 10% of the time gives $\sigma=0.5$. We take $\sigma=0.4$ as the best-fit, and plot it as a hexagon in Figure 8. The shade indicates the possible range (i.e., from 0.3 to 0.5 in this example).

As seen from Figure 8, the H$\alpha$-to-UV flux ratio almost always gets the large fluctuating amplitude right, and gets the small fluctuating amplitude quite close ($\sigma\sim 0.5$ dex when the truth is 0.4 dex). Notably, it never mistakenly infer one model with true $\sigma=0.4$ dex as $\sigma\sim 0.8$ dex, and vice versa. This performance is markedly better than the individual SED fits, where the spread in the posterior medians can be larger than the difference between these two amplitudes. However, such discerning power vanishes when considering the other two parameters, namely duration and slope, as the best-fits only differ marginally for the different cases. Put another way, the H$\alpha$-to-UV flux ratio only has constraining power on the fluctuating amplitude: it is not very constraining for either the timescale of bursts or the underlying SFH slope.

These results demonstrate that the H$\alpha$-to-UV flux ratio is sensitive to recent bursts, which in turn can be linked to many physical processes driving the SFH, and that it is difficult to disentangle them based on the distribution of H$\alpha$-to-UV flux ratio alone. This challenge has also been implied in V. Mehta et al. (2023), where they found that a smoothly rising or declining SFR over the long term can produce similar scatter in H$\alpha$-to-UV flux ratios as short-term bursts.

In summary, the H$\alpha$-to-UV flux ratio is a useful indicator for burstiness, but going beyond inferring a general characteristic of SFH variability (e.g., the fluctuation amplitude in our case) to constrain all the population-level SFH parameters requires additional information. This motivates us to explore a simultaneous population-level approach, as detailed below.

The physical processes of various characteristic timescales influence certain spectral features in different ways. Figure 3 illustrates this with four examples of model SFHs and their corresponding spectra, and the spectral features of Balmer breaks, Balmer emission lines, near- and far-ultraviolet flux densities (NUV, FUV) are highlighted. While isolated analysis of individual spectra is insufficient to robustly infer the population-level SFH parameters, simultaneous analysis of the distribution of the relevant spectral features may give greater discriminating power. This is because while determining the behavior of a single galaxy is challenging, but many galaxies living within the same space can effectively trace out phase-space to provide a comprehensive view of the population-level characteristics.

The above serves as a motivation for our proposed population-level approach; at the same time, it is also a key assumption that the formation histories of most galaxies vary on similar timescales. Here we show that galaxies drawn from the same SFH family can be used to identify that family. In real observations, one must somehow also first group galaxies into different SFH families.

Simulations suggest that this may be challenging, as the probability of a burst (at fixed stellar mass) depends on a stability parameter, which involves gas mass and its specific angular momentum (E. Cenci et al., 2024). In other words, even at fixed mass the burstiness can depend on other galaxy properties such as their structure (see also P. F. Hopkins et al. 2023b). However, as a first-order analysis, we may assume that galaxies in similar redshift and mass ranges likely experience comparable environmental conditions, leading to similar SFH variations.

Having discussed the individual and population-level fits separately in Sections 3.4, 4, 5, and 6.3, we now tie all the pieces together. We begin this section by a brief summary, presenting the implications from this work for solving the longstanding challenge of outshining. Key lingering questions and future improvements are discussed in Section 7.2.

The remaining two subsections supply additional discussion on the effects of priors placed on the SFH in individual SED fitting, and the implications of burstiness for studies assessing completeness.

This paper establishes the first step toward solving outshining—a longstanding challenge in interpreting the light from galaxies—by empirically measuring and then applying a model for short-term galaxy burstiness. This issue has become especially pertinent in light of recent observations from JWST of the early mini-quenched galaxies (e.g., T. J. Looser et al. 2024; V. Strait et al. 2023; R. Endsley et al. 2024; A. Weibel et al. 2025; W. M. Baker et al. 2025; J. A. A. Trussler et al. 2025) and the over-abundance of luminous high-redshift galaxies (e.g., R. P. Naidu et al. 2022; H. Atek et al. 2023; C. M. Casey et al. 2023; S. L. Finkelstein et al. 2023; B. E. Robertson et al. 2023).

We present a thorough examination of the extent of constraints attainable regarding bursty star formation, defined here as recurrent up-and-down fluctuations on tens of Myr timescales over the past 500 Myr of formation history. We begin our study by assessing the state-of-the-art SED-fitting methods for inferring burstiness, and find that standard techniques with neutral priors fail to recover detailed recent SFHs and often under-estimate the total stellar mass; however, a successful recovery is achievable with the correct model for burstiness. We then introduce a population-level technique leveraging the distributions of timescale-sensitive spectral features, and demonstrate its potential for accurately inferring burstiness. While motivated by high-redshift findings, our results are also applicable to studies of lower-redshift galaxies. The main findings are summarized as follows.

First, in the presence of bursty star formation, the population-level SFH parameters, particularly the fluctuation amplitude, lead to strong selection effects in flux-limited surveys. The simulated brightness from our SFH models, for a galaxy population of the same redshift and mass, span over one magnitude in rest near-infrared, where the fluxes are expected to be the most stable. The variability in rest-optical is greater—more than two magnitudes. This may explain the prevalence of bursty star formation at high redshift in objects typically selected in the rest-optical, or in spectroscopic samples with more complex selection functions.

Second, even with exquisite data with no modeling systematics, individual SED fits are unable to place meaningful constraints on the population-level SFH parameters. The spread in posterior medians in the inferred fluctuation amplitude, duration, and slope, from fitting an ensemble SEDs simulated from one SFH family, range from 0.5 to 1.0 dex, 10 to 120 Myr, and 0.5 to 1.2, respectively. These spreads are greater than the spacing of the model grid, meaning that the sampled SFH parameters cannot be distinguished. Critically, if these were real data, we would interpret the wide spread in posterior medians as evidence of a diversity in galaxy SFH timescales. This is especially concerning when the true SFH parameters fall outside the posterior coverage (Figures 4–5).

Third, the above findings are consistent for the tests done using the continuity and the bursty continuity priors. Surprisingly, the continuity prior performs marginally better than the bursty continuity prior. Assuming the continuity prior, the mean biases in the inferred masses, and SFRs averaged over the most 100 Myr are approximately $-0.06$ dex, and $0.1$ dex, respectively; assuming the bursty continuity prior, the biases increase to $-0.3$ dex and $0.2$ dex. This is because by preferring smooth transitions between time bins in the SFH, the continuity prior is less prone to outshining.

Fourth, the population distributions of spectral features—Balmer break strength, Balmer emission lines, NUV/FUV flux densities—are sensitive to at least some of the population-level SFH parameters. Considering these four observables simultaneously, we can always identify the correct population-level SFH parameters via the shortest Wasserstein distance. This work serves as a proof of concept, and the full methodology and validation of our proposed population-level method will be presented in a future work.

Finally, under the assumption that current population-level burstiness predicts past SFH, the complete formation history, even when the light from recent bursts dominates the SED, can be accurately inferred given the correct burstiness model. In other words, the proposed population-level method for constraining burstiness is the key to address outshining.

To conclude, bursty star formation has emerged as a potential unifying explanation for the high-redshift temporarily quenched galaxies and overly luminous galaxies discovered by JWST; yet observational constraints on burstiness remain uncertain. Outshining represents perhaps the greatest challenge in interpreting the light from galaxies experiencing bursty star formation. The simultaneous population-level method proposed in this paper has the strong potential to constrain burstiness, and, thereby, addressing the issue of outshining. This work paves the way for a complete understanding of the observed diverse populations of galaxies.

We thank Elia Cenci for providing the FIREbox SFHs. We also thank Hiranya Peiris, Chris Hayward, and Kartheik Iyer for valuable discussions. We are grateful to the anonymous referee for the helpful comments. B.W. and J.L. acknowledge support from JWST-GO-02561.022-A. Computations for this research were performed on the Pennsylvania State University’s Institute for Computational and Data Sciences’ Roar supercomputer. This publication made use of the NASA Astrophysical Data System for bibliographic information.

In the main text, we adopt the same time bins as used in the input SFH when fitting SEDs. However, as it is difficult to know the timescales that the SFH varies a priori, certain form of rebinning is likely required.

J. Leja et al. (2019) has tested different number of time bins extensively on smoothly varying SFHs, and concluded that rebinning has no impact on the SED-fitting results. We revisit this issue here given the more complex bursty SFHs tested in this work. We rebin the input SFH into 5 Myr fine bins in the first 100 Myr, and the rest into a uniform logarithmic grid with the last bin fixed to 90% of the age of the universe. This choice is motivated by the need to be sufficiently flexible to describe SFHs potentially varying on a range of timescales, while keeping the total number of bins (and hence the total number of free parameters) to be within a manageable number. As shown in Figure A.1, the input and the rebinned SFHs predict different model spectra. While alternative non-parametric SFH implementations exist (e.g., K. A. Suess et al. 2022; K. G. Iyer et al. 2024), we demonstrate, in this abbreviated manner, that the binning of SFH would need to be carefully selected when inferring SFH variability in SED fitting.

We adopt the same general Prospector settings as outlined in Section 3.1, but using flexible time bins (K. A. Suess et al., 2022), which allows the bin widths to change, so that the period where the SFH is inferred to have the most variations can be more described in greater resolution.

Four observing strategies are considered: photometry only (including broad and medium bands), photometry + low-resolution Prism spectroscopy, photometry + medium-resolution G235M spectroscopy, and photometry + medium-resolution G395M spectroscopy. At our redshift range of interest ($z\sim 4-7$), the G235M grating covers a critical age indicator—the Balmer break—and also H$\beta$, whereas the G395M grating covers H$\alpha$ and H$\beta$.

The results, where we assume the ideal scenario of dust attention being known perfectly, are presented in Figure B.1. This rejuvenating case—where short periods of star formation and quiescence occur in the recent past—produces a model spectrum in which the most recent burst dominate the observed SED. Specially, this galaxy experiences a period of star formation over the most recent 340–540 Myr, followed by a 200 Myr long quiescence, and then forms stars again during the most recent 40 Myr.

Without the complication of the degeneracy between dust attention and emission-line strength, the most recent SFHs are inferred accurately in all cases. However, outshining induces significant uncertainty in inferring the past SFHs.

We also perform the same fits but letting dust attention to be fit for. The inferred becomes all consistent with a single burst, and the most recent SFRs are overestimated by a factor of $\sim 2$.

The challenge being put forth by the complex likelihood surface is manifested in the underestimated uncertainties in Figures 6 and D.1. To verify the hypothesis that the sampler fails to locate the global maximum on the likelihood surface, we fit the same galaxy multiple times. The only randomness in this process is the noise—while the S/N is always set to 20, we let the realized noise to be generated randomly for each trial. As seen from Figure C.1, the inferred masses as well as their uncertainties can change significantly. In addition, we repeat the same exercise using nautilus, which is also a nested sampler but with enhanced efficiency aided by neural networks (J. U. Lange, 2023). Interestingly, nautilus tends to consistently identify the underestimated mass as the final solution.

A rigorous test of various samplers is out of the scope of this paper; however, our findings suggest that the sampling problem here is indeed a challenging one. The bias in mass recovery is perhaps primarily driven by the challenge in finding the global maximum on the likelihood surface. High-dimensional SED modeling is known to have a complex likelihood space, and the complex behavior of bursty and rising SFHs exacerbates this challenge. As discussed in B. Wang et al. (2024b), the nested sampling method (J. Skilling, 2004), which is also used in this paper, is already better suited to sample multi-modal posteriors than other traditional techniques such as Markov chain Monte Carlo (J. Goodman & J. Weare, 2010). In principle, the problem of locating the global maximal likelihood modes can be mitigated by substantially increasing in the accuracy settings, but the resulting increase in the CPU-hours would quickly become prohibitively expensive.

Possible improvements on the current state-of-art sampling process used in astronomy (F. Feroz et al., 2009; D. Foreman-Mackey et al., 2013; W. J. Handley et al., 2015; J. S. Speagle, 2020; J. Buchner, 2021) include gradient-based samplers (M. D. Hoffman & A. Gelman, 2011), combining with machine-learning techniques (M. Karamanis et al., 2022; J. U. Lange, 2023), or new inference workflows (C. Hahn & P. Melchior, 2022; G. Khullar et al., 2022; B. Wang et al., 2023a; K. Zhang et al., 2023; M. Ho et al., 2024). It is worth noting that this missing-mode issue does not undermine the main conclusion in the paper that the individual spectrum does not contain enough information for complex recent SFH recovery. We have shown that all the spectra are well fit, but the inferred SFH can be very different from the input SFH. A sampler that is able to better map likelihood space can improve the uncertainty calibration, but unlikely to solve the fundamental problem regarding the information content.

The long-wavelength coverage enabled by MIRI improves the recovery of mass and averaged SFR. The bias typically is improved by $0.1-0.2$ dex for the inferred masses, and $<0.06$ dex for the SFRs. This is likely because the mass-to-light ratio shows minimal variation in this wavelength range, particularly at rest 2 $\mu$m (e.g., T. Kettlety et al. 2018). However, the MIRI bands are still insufficient to accurately trace the SFH variations from individual SED fits.