It's More Complicated Than You Think: A Forward Model to Infer the Recent Star Formation History, Bursty or Not, of Galaxy Populations

TL;DR

This work tackles the challenge of constraining galaxy star-formation histories (SFHs) at the population level by inferring the power and timescales of SFR fluctuations from JWST-like spectroscopic observables. It introduces a simulation-based inference framework that forward-models galaxy populations through two SFH models: a simple single-frequency oscillator and a flexible power spectral density (PSD) model spanning 1 Myr–10 Gyr, augmented by population-level recent-SFH slope parameters. Using forward-modeled rest-UV to rest-optical features measured with FSPS/Prospector and realistic noise, the authors train neural density estimators (normalizing flows) to recover posterior distributions of PSD components and slope parameters, achieving precise recovery for both bursty and smooth histories and demonstrating the framework's ability to distinguish FIRE-2-like versus Illustris-like feedback prescriptions. The study also reveals limitations: long-timescale power is harder to recover in bursty populations due to outshining, and systematic uncertainties in SPS and dust require careful treatment. Overall, this method provides a robust, scalable pathway to constrain feedback physics governing star formation from large, uniformly selected spectroscopic samples with JWST.

Abstract

Observations of the early Universe (z > 4) with the James Webb Space Telescope reveal galaxy populations with a wide range of intrinsic luminosities and colors. Bursty star formation histories (SFHs), characterized by short-term fluctuations in the star formation rate (SFR), may explain this diversity, but constraining burst timescales and amplitudes in individual galaxies is challenging due to degeneracies and sensitivity limits. We introduce a population-level simulation-based inference framework that recovers the power and timescales of SFR fluctuations by forward-modeling galaxy populations and distributions of rest-UV to rest-optical spectral features sensitive to star formation timescales. We adopt a stochastic SFH model based on a power spectral density formalism spanning 1 Myr-10 Gyr. Using simulated samples of N=500 galaxies at z~4 with typical JWST/NIRSpec uncertainties, we demonstrate that: (i) the power of SFR fluctuations can be measured with sufficient precision to distinguish between simulations (e.g., FIRE-2-like vs. Illustris-like populations at >99% confidence for timescales < 100 Myr); (ii) simultaneously modeling stochastic fluctuations and the recent (t_L < 500 Myr) average SFH slope is essential, as secular trends otherwise mimic burstiness in common diagnostics; (iii) frequent, intense bursts impose an outshining limit, and bias inference toward underestimating burstiness due to the obscuration of long-timescale power; and (iv) the power of SFR fluctuations can be inferred to 95% confidence across all timescales in both smooth and bursty populations. This framework establishes a novel and robust method for placing quantitative constraints on the feedback physics regulating star formation using large, uniformly selected spectroscopic samples.

Figures (16)

• Figure 1: Effect of the single-frequency illustrative model on SFR indicator distributions. Left panel: Three SFH model realizations, with $\alpha=0.1$, $\sigma=-0.1$, and $\delta t = 20$ Myr (blue), $\alpha=-0.3$, $\sigma=1.0$, and $\delta t = 40$ Myr (red), and $\alpha=-0.3$, $\sigma=1.0$, and $\delta t = 20$ Myr (teal dashed) in a burst phase. The gray dashed line indicates the SFMS at $z=4$ to be sSFR$_\text{SFMS} = 10^{-8.8} \ \mathrm{yr^{-1}}$. Right panel: The marginal and joint distribution of stellar mass normalized H$\alpha$ flux and the Balmer break strength for 500 samples of each SFH model realization. We marginalize over stellar metallicity and the SFH phase for each sample. The distribution of H$\alpha$/UV flux from each population model is shown in the bottom right panel.
• Figure 2: Effect of the flex-PSD model on observable distributions. Upper left panel: Two realizations of the flex-PSD model, a bursty model (blue, power-law slope of $\beta = 1$), and a smooth model (red, power-law slope of $\beta= 2$). For the smooth model, we display the log-normal Gaussian components as dashed lines with each normalization annotated with the flex-PSD model parameter name. We additionally show the approximate timescales that are probed by the dominant stellar spectral type in the integrated light of an SED: O stars on 1-10 Myr timescales, B stars on 10-100 Myr timescales, A stars on 100-500 Myr timescales, F stars on 0.5-2 Gyr, and G stars $> 2$ Gyr. Upper Right Panel: Two example SFHs sampled from the two models. Lower panels: Marginal distributions of observed spectral features from mock populations simulated from each PSD. We include stellar-mass normalized H$\alpha$ flux, stellar-mass normalized FUV flux density, Balmer break strength, and rest-frame U-V color.
• Figure 3: Workflow of the proposed method. A training set is generated by sampling parameters $\theta$ from a chosen population-level SFH model, including a distribution of stellar metallicities and redshifts in which to marginalize over. For a given sample of parameters, a population of $N$ SFHs is sampled from the SFH model, and the nuisance parameters are marginalized over. These parameters are passed to an SPS framework to model the observed spectrum of each object. Spectral features of interest are measured, with observational and counting noise applied. The SBI model is then trained to learn the mapping from these observable distributions to the posterior distribution of $\theta$. To apply this method to observed galaxy populations, $N$ observed spectra are corrected for dust attenuation, the same spectral features are measured, and the resulting distributions are passed to the trained model to infer posteriors on the selected population parameters.
• Figure 4: Effect of flex-PSD and SFH slope parameters on observable distributions. Each row displays the change in the PSD/SFHs and the observable distributions as we decrease (blue) or increase (red) the parameter by some amount. The black curve is the same population model in all panels.
• Figure 5: SBI model residuals for the single-frequency SFH model. Upper Row: Posterior median $\pm 1\sigma$ versus the true value of the hyperparameter. From left to right, $\sigma$, $\delta t$, and $\alpha$. A one-to-one line is displayed as a black dashed line for reference. $\sigma_\text{NMAD}$, the normalized median absolute deviation, is shown in the top left corner. Lower Row: Distribution of the standardized residuals for each simple SFH hyperparameter. A unit-Gaussian is displayed as a black dashed line for reference.