Stochastic Modeling of Star Formation Histories. III. Constraints from Physically Motivated Gaussian Processes

Galaxy formation and evolution involve a variety of effectively stochastic processes that operate over different timescales. The extended regulator model provides an analytic framework for the resulting variability (or “burstiness”) in galaxy-wide star formation due to these processes. It does this by relating the variability in Fourier space to the effective timescales of stochastic gas inflow, equilibrium, and dynamical processes influencing giant molecular clouds' creation and destruction using the power spectral density (PSD) formalism. We use the connection between the PSD and autocovariance function for general stochastic processes to reformulate this model as an autocovariance function, which we use to model variability in galaxy star formation histories (SFHs) using physically motivated Gaussian processes in log star formation rate (SFR) space. Using stellar population synthesis models, we then explore how changes in model stochasticity can affect spectral signatures across galaxy populations with properties similar to the Milky Way and present-day dwarfs, as well as at higher redshifts. We find that, even at fixed scatter, perturbations to the stochasticity model (changing timescales vs. overall variability) leave unique spectral signatures across both idealized and more realistic galaxy populations. Distributions of spectral features including Hα and UV-based SFR indicators, Hδ and Ca H and K absorption-line strengths, D n (4000), and broadband colors provide testable predictions for galaxy populations from present and upcoming surveys with the Hubble Space Telescope, James Webb Space Telescope, and Nancy Grace Roman Space Telescope. The Gaussian process SFH framework provides a fast, flexible implementation of physical covariance models for the next generation of spectral energy distribution modeling tools. Code to reproduce our results can be found at https://github.com/kartheikiyer/GP-SFH.


INTRODUCTION
Large galaxy surveys like SDSS (York et al. 2000;Strauss et al. 2002;Abazajian et al. 2009), GAMA (Driver et al. 2009) and COSMOS (Scoville et al. 2007;Lilly et al. 2007;Weaver et al. 2022) reveal an enormous diversity in galaxy demographics across different epochs, and many studies in modern galaxy evolution have devoted considerable energy in attempting to explain this diversity from a physical standpoint, usually through analytical models and cosmological simulations (see review by Somerville & Davé 2015).
A part of this picture involves the stochasticity in star formation, which is regulated by a variety of physical processes acting over many orders of magnitude in spatial and temporal scales.This stochasticity or 'burstiness' can be physically observed through short-timescale star formation rate (SFR) indicators such as Hα or UV-based SFRs, which probe recent star formation averaged over the most recent ∼ 4 − 10 and ∼ 10 − 100 Myr respectively (Madau & Dickinson 2014;Flores Velázquez et al. 2021;Tacchella et al. 2022a).It can also be observed by studying resolved star formation in galaxies, where most star formation appears to to occur in discrete clumps traced by the rest-UV that are spatially offset from the bulges of galaxies where most older stellar populations live (Guo et al. 2014;Huertas-Company et al. 2020), with the creation and destruction of these clumps due to ISM physics and feedback leading to stochasticity on the timescales of the clump lifetimes (Semenov et al. 2018(Semenov et al. , 2021)).On a population level, this could leave signatures in the scatter of scaling relations like the SFR-M * correlation1 (Noeske et al. 2007;Daddi et al. 2007;Elbaz et al. 2007) and the mass-metallicity relation (Tremonti et al. 2004), which show overall coherent behaviour for the 'average' galaxy that is tied to the growth of their parent dark-matter halos, but significant variation (at the ∼ 0.3 − 0.5 dex level) in the star formation rates of individual galaxies that seem to fluctuate around these average relations (Kauffmann et al. 2006;Tacchella et al. 2016;Matthee & Schaye 2019).
These fluctuations in SFR are regulated by a variety of physical processes ranging from the local creation and destruction of stars in giant molecular clouds (GMCs), to dynamical processes like disk formation and bulge growth, to galaxy-wide processes that include mergers, galactic winds from stellar and AGN feedback, and baryon cycling that couples a galaxy to its surrounding circumgalactic medium (CGM) (Tacchella et al. 2020;Iyer et al. 2020).On the largest scales, however, a galaxy's growth is tied to the reservoir of fuel available to it to form stars, which are in turn tied to the accretion rates of their parent halos, galactic depletion-times and outflows, and the large-scale structure of the environment they live in (see review by Wechsler & Tinker 2018).
Analyzing the stochasticity of star formation across a range of timescales 2 therefore provides us with a way to constrain the relative strengths of these physical processes.A particularly effective way to quantify and assess this is by quantifying the 'burstiness' or stochasticity of a galaxy's past star formation as a function of timescale, and comparing these to theoretical predictions.Fourier space has proved instructive in this regard, with the power spectral density (PSD) of galaxy star formation over time being used to construct analytical models of galaxy stochasticity due to different physical processes (Forbes et al. 2014;Tacchella et al. 2020), study the regulation of star formation across different cosmological simulations (Iyer et al. 2020), and constrain simple models of stochasticity using observational data from large galaxy surveys (Caplar & Tacchella 2019;Wang & Lilly 2020a).
These studies hit on a fundamental aspect of galaxy evolution, that the growth of pure dark matter halos is essentially a scale-free process that leads to a power-law PSD (Guszejnov et al. 2018;Kelson et al. 2020).This is the reason why (to first order) the halo mass is such a good predictor of galaxy properties and methods like sub-halo abundance matching have met with such remarkable success.Other aspects including stellar feedback, baryon cycling and multiphase ISM astrophysics decouple star formation from this hierarchical build-up and add additional stochasticity to this on a range of (often interwoven) timescales, leading to an overall complex power spectrum that can be studied and understood with careful analysis.Caplar & Tacchella (2019, Paper I in this series, hereafter CT19) defines the PSD of a galaxy's star formation history (SFH) and assuming the shape of a simple broken-power law, linked it to SFR distributions averaged over different timescales.Tacchella et al. (2020, Paper II in this series, hereafter TFC20) builds on this, using the widely successful 3 gas regulator model (Lilly et al. 2013) coupled with the stochastic inflow of gas (Kelson et al. 2020, TFC20) to derive a more general form for the PSD of galaxy SFHs.
The PSD of the Extended Regulator model depends only on an overall level of stochasticity for gas inflow and GMCs, and characteristic timescales for effective gas inflow, equilibrium, and GMC lifetimes.TFC20 formulates this PSD and proposes the timescales for a few galaxy populations (Milky-way analogues, dwarf galaxies, galaxies at high redshift and massive galaxies at cosmic noon) given our current knowledge of burstiness in these galaxy populations.However, to verify these models and observationally probe these timescales, we need a framework in which the PSDs can be tested observationally.
To observationally measure the PSD (i.e. the variability in SFR over different timescales), we can leverage the fact that a range of spectral features measure changes in the star formation rate averaged over different timescales.Some of the commonly measured spectral features include the nebular Hα line, rest-UV+IR SED that traces emission from young massive stars & re-emitted radiation from dust, Hδ absorption from the photospheres of smaller (mostly Atype) stars, and the 4000 Å break from the accumulation of ionized metal absorption lines for older stellar populations (Gonçalves et al. 2012;Flores Velázquez et al. 2021).Taken together, these indicators probe fluctuations in SFRs on timescales ranging from ∼ 5 Myr to 10 Gyr, and their ratios have previously been used to obtain estimates of the burstiness of galaxy populations (Guo et al. 2016;Emami et al. 2019a;Broussard et al. 2019a;Faisst et al. 2019;Wang & Lilly 2020b).
Spectral energy distribution (SED) fitting methods go a step further and estimate the full star formation histories (SFHs) of individual galaxies using the full range of available multiwavelength spectral information, whether it is photometry, spectroscopy, or a combination of the two (Heavens et al. 2000;Tojeiro et al. 2007;Dye 2008;Pacifici et al. 2012;Pacifici et al. 2016;Smith & Hayward 2015;Leja et al. 2017;Iyer & Gawiser 2017;Carnall et al. 2018;Leja et al. 2019a;Iyer et al. 2019).A few of these methods that implement non-parametric SFHs (Pacifici et al. 2012;Leja et al. 2017;Iyer et al. 2019) also allow users to implement priors on SFR burstiness on one or more timescales.The non-parametric Dense Basis method (Iyer & Gawiser 2017;Iyer et al. 2019) for SFH reconstruction that allows us to incorporate physical priors on SFR stochastcity through a Gaussian process covariance function (also called the kernel).These kernels are related to the PSDs through the Weiner-Khinchin theorem4 (Wiener 1930;Khinchin 1938), which allows us to relate the frequency-domain PSDs to the time-domain auto-covariance functions (ACFs).
For the first time, this opens up a way to explicitly incorporate an analytical framework for correlated SFRs over a range of timescales into an SED modeling framework.This is a crucial development for multiple reasons: 1.It provides a physical explanation for how SFRs vary over time through the three timescales -the timescale over which stochastic gas inflow is correlated, the mass-loss or equilibrium timescale on which gas is consumed/removed from the reservoir, and and average GMC lifetime timescales.While a range of different physical processes are responsible for regulating SF in galaxies, these three timescales have been shown to capture a significant portion of the effective dynamics of galaxy populations5 (TFC20, Tacconi et al. 2020).
2. It offers a framework to forward-model galaxy observations based on a stochasticity model to compare against existing models and to determine sensitive spectral features for future observations (e.g.Whitaker et al. 2014;Parul & Bailin 2021).
3. It allows us to explicitly compare against existing priors e.g., uncorrelated SFRs in adjacent timesteps, or the Dirichlet & continuity type priors (Leja et al. 2017;Iyer et al. 2019) that assume an arbitrary amount of stochasticity and/or correlation between adjacent SFR bins in a model to test their efficacy and their relation to the effective timescales in the regulator model.
4. It provides an intuitive framework for modeling SFH priors, or incorporating SFH priors from cosmological simulations by estimating the extended regulator model parameters directly from their SFHs.
In this paper, we incorporate the physically motivated SFR stochasticity model proposed in TFC20 within the framework of a Gaussian process (GP; Rasmussen & Williams 2006;Iyer et al. 2019).Using this, we then use the flexibility of the Extended Regulator model of TFC20 to define GPs corresponding to various regimes of stochasticity that we might find in galaxy populations -ranging from the bursty behaviour of galaxies at high redshifts to the long-timescale correlated behaviour of Milky Way analogues at lower redshifts.We then use these GPs to generate mock SFHs in a computationally inexpensive manner, which is crucial if these are to be used in SED fitting.By running these SFHs through a stellar population synthesis framework (FSPS; Conroy et al. 2009;Conroy & Gunn 2010;Johnson et al. 2021), we then model spectra corresponding to these SFHs and use them to identify observables that can be used to tell the models apart, laying the foundations of future work where this can be incorporated into SED fitting packages.A schematic representation of our main approach is shown in Figure 1, and the code used to implement the GPs is publicly available at https://github.com/kartheikiyer/GP-SFH.
This paper is structured as follows.In §2, we go through the formalism and provide a description of the Extended Regulator model first presented in TFC20, derive the associated auto-covariance function (ACF), and highlight specific cases likely to correspond to various galaxy populations of interest.In §3, we describe the implementation of the derived ACF as a "physical kernel" in a Gaussian process (GP), and describe how it can be used as a prior in SED fitting codes that have flexible models for galaxy SFHs.In §4, we investigate how differences in the underlying ACFs can translate over to spectral signatures using stellar population synthesis (SPS) models, and investigate how these differences manifest in populations of simulated galaxies.We discuss our findings in §5 and conclude in §6.
Note that while it is most natural to have a process using the base-e logarithm ln SFR(t), we convert to the base-10 logarithm log SFR(t) in most plots to facilitate comparisons with other quantities and measurements in the literature.Throughout this paper magnitudes are in the AB system; we use a standard ΛCDM cosmology with Ω m = 0.3, Ω Λ = 0.7 and H 0 = 70 km s −1 Mpc −1 .

MODELLING STOCHASTICITY IN STAR FORMATION RATES
We start by building physical intuition of how different physical processes related to galaxies can affect stochasticity and correlations in the SFRs of individual galaxies across cosmic time, summarized through their power spectral densities (PSDs) and associated auto-covariance functions (ACFs).In §2.1, we provide a brief set of definitions for the PSD and ACF and their relationship with each other.In §2.2, we derive results for the Extended Regulator model presented in TFC20.Finally, in §2.3 we highlight the four special cases from TFC20 -Milky Way analogues, dwarfs, galaxies at cosmic noon, and galaxies at high-redshift -along with a description of their general expected behavior.A more detailed set of derivations are presented in §A.
Following TFC20, we will assume that the stochastic processes described by the Extended Regulator model correspond to variability in the natural log of the star formation rate, ln SFR around a time-dependent mean µ(t) and are captured entirely by the corresponding PSD/ACF.This will allow us to model SFHs as log-Gaussian Processes, which we will return to in Section 3.

Terms and Definitions
We start by informally defining a stochastic process as something that can generate infinite realizations of a time series {x 1 , x 2 , . . .x n } ≡ {x t } n 1 ≡ x n at any times t = 1, . . ., n (i.e. the x t values change every time we simulate from the process).The collection of x n values will then follow some joint probability distribution P (x n ) which is defined by the stochastic process.
The simplest way to explore the correlation structure in a given stochastic process is to compute the auto-covariance function (ACF)6 (1) between x t and x t at two times t and t , where µ(t) is the time-dependent mean and P (x t , x t ) is the joint distribution of x t and x t defined by the process.Assuming our process is stationary so that the auto-covariance function only depends on the separation (i.e.time lag) between any two given times τ ≡ t − t rather than the individual times t and t themselves, we can instead write the auto-covariance function as In addition to defining correlation structure as a function of time t, we can also do the same as a function of frequency f .We first define a "windowed" version of x(t) (3)  Both the PSD and ACF highlight that variability on short timescales is dominated by giant molecular clouds (GMCs) while variability on longer timescales is dominated by a combination of gas inflows and cycling between atomic and molecular gas.
Using the ACF, we define a Gaussian Process (GP; bottom left) to quickly and easily sample many realizations of the log SFR(t) over time relative to some mean log SFR base (t), which here we take to be 1M /yr.We again highlight both long-term and short-term (see inset) variability driven by gas inflows/cycling and dynamical processes, respectively.Finally, we generate SEDs at z ∼ 0 by running the GP-sampled SFHs through a stellar population synthesis model (bottom right).The variability in the SFHs leads to differences in the overall normalization (from differences in total stellar mass formed) as well as particular spectral features (from differences in relative contributions of various stellar populations).Properties apart from the SFH assumed while generating the spectra are detailed in Section 3.2, and are held constant to highlight differences from the SFHs.Understanding how the distributions of spectral features corresponding to various galaxy populations correspond to the properties of their PSDs will allow us to put constraints on the timescales on which gas inflows, baryon cycling, and GMC physics drive galaxy evolution.
for a window function w T (t) with some width (duration) T centered around t. Taking its Fourier transform the power spectral density (PSD) is then where the limit T → ∞ assumes the stochastic process is not localized in time.We can interpret the PSD as the relative amount of variance as a function of frequency, where larger values indicate stronger correlations at particular frequencies.
While the ACF and PSD can be computed directly from a given stochastic process, they can also be directly computed from each other.Based on the Wiener-Khinchin theorem, in the continuous-time limit the PSD S(f ) and ACF C(τ ) are Fourier pairs and we can convert between the two using: This property is extremely useful, as many stochastic processes (such as the Extended Regulator model in §2.2) can be much easier to describe in frequency rather than in time (and vice versa).

Extended Regulator Model
TFC20 introduced the Extended Regulator model as a way to characterize how stochastic processes that drive 1. gas inflow rates, 2. gas cycling (between atomic and molecular gas) in equilibrium, and 3. the formation and disruption of giant molecular clouds (GMCs) relate to ln SFR(t).Assuming that the behavior of each component follows a damped random walk with some decorrelation timescale τ dec and variability σ, each PSD can be shown to have a broken power-law PSD of the form where s 2 is the absolute normalization (scatter squared) for f = 0. Making the well-justified assumptions that 1. the process describing the behavior of GMCs is largely independent of those describing gas inflow and cycling in equilibrium and 2. the processes describing gas inflow and equilibrium gas cycling are coupled, the full PSD of the Extended Regulator model can be written as where s 2 gas = s 2 in s 2 eq is the total variability in gas inflows and equilibrium cycling, s 2 dyn is the variability in dynamical processes including the creation and destruction of GMCs, and τ in , τ eq , and τ dyn are the de-correlation timescales associated with gas inflows, cycling in equilibrium, and GMC formation/disruption respectively, and we use the values of β l = 0, β h = 2 for the power-law slopes of the gas inflow term as defined in Eqn.23 and Table 2 of TFC20.
Using the Wiener-Kninchin Theorem, the corresponding ACF is then 3).To be comparable across different timescales, the four ACFs have all been scaled such that C(τ = 0) ≈ 1.0 dex (left panel ).We then generate a single realisation of the SFH log SFR(t) using a GP with each of the four ACFs with a fixed random number seed to highlight differences that arise from varying the ACF (right panel ).ACFs with longer-timescale correlations tend to spend larger stretches of time above and below the mean SFH, while short-timescale kernels tend to show more "bursty" behaviour that generate larger (but short-lived) excursions from the mean.The Milky Way analogue model in particular includes correlations at both long and short timescales from gas cycling and inflows, respectively, giving it the smoothest, most correlated SFH.
for all τ in = τ eq and where we've replaced s → σ to emphasize that S(f = 0) = s 2 = σ 2 = C(τ = 0).See §A for further details and more general results.
Compared to the PSD, the ACF offers different insights into the correlation structure.In particular, it shows that a damped random walk leads to correlations that decay exponentially with time (∝ e −|τ |/τ dec )7 .Note also that the variance of the Extended Regulator model now becomes since there are now multiple independent stochastic processes involved.
We can quantify the extent and strength of the correlations using the auto-correlation time which is a measure of the relative correlation contributed by all possible time lags τ .For the Extended Regulator model, evaluating this expression gives where we've defined a gas as the fractional contribution to the variance from the gas component.In the limit where σ gas σ dyn so GMC formation/disruption is the dominant source of variability, this gives which is just twice the GMC formation/disruption time (since τ can range from −∞ to +∞).In the limit where σ gas σ dyn , we instead have τ A,gas = 2τ gas = 2(τ in + τ eq ) ( 14) which is related instead to the combined timescales involved in gas inflows τ in and cycling τ eq .This final case, where the GMC contribution to the variability is assumed to be negligible, is what TFC20 refer to as the Regulator model.
The PSD and ACF for the Extended Regulator model are shown in Figure 1 while comparisons between the Regulator model (two damped random walks) and a single damped random walk are shown in Figure 2.

Special Cases
Following TFC208 , we consider the following special illustrative cases to highlight the behavior of the Extended Regulator model in four different regimes: 1. Milky Way Analogue (MWA): (τ in , τ eq , τ dyn ) = (150, 2500, 25) Myr.Based on the long-term secular evolutionary trends seen in the MW, this includes an extremely large τ eq (2.5 Gyr) along with approximate order of magnitude (> 5) differences between various timescales, with τ eq τ in τ dyn .SFHs will be dominated by the long-running equilibrium timescale, with small perturbations from changes to the inflow rate along with small amounts of additional white noise from GMCs on much shorter timescales.
2. Dwarf : (τ in , τ eq , τ dyn ) = (150, 30, 10) Myr.Although it has the same τ in as MW, τ eq has substantially decreased to account for the much more rapid gas cycling (and τ dyn to a lesser extent for similar reasons) expected in these low-mass galaxies.This has the effect of making the SFHs substantially burstier on short timescales ( 100M yr).
It also includes smaller changes in scale (∼ 5), leading to somewhat larger impacts from τ in > τ eq > τ dyn .SFHs will be dominated by the variability timescales associated with gas inflows, but with larger perturbations from equilibrium and white noise from GMCs compared to our Milky Way analogue.
3. Cosmic Noon: (τ in , τ eq , τ dyn ) = (100, 200, 50) Myr.This case is designed to simulate a typical 10 9 M galaxy around z ∼ 2. The equilibrium time is larger by an order of magnitude relative to the dwarf case due to the larger overall mass, with smaller changes in scale (∼ 2) and longer-lived GMCs.As τ eq τ eq τ dyn , all timescales remain quite relevant, leading to larger and more correlated fluctuations.4. High-z: (τ in , τ eq , τ dyn ) = (16, 15, 6) Myr.Our last case is designed to simulate the conditions for a galaxy at z ∼ 4 − 6.Here, τ eq ≈ τ in with both only ∼ 2τ dyn , with extremely short timescales due to the lower masses involved along with the more disruptive environments that many galaxies find themselves in.Since the gasrelated timescales are almost identical, we expect this case to have behavior most similar to a Matern32 kernel (see §A.4.1) but with additional perturbations caused by GMCs on somewhat similar timescales.
We consider two classes of models when deciding on the values of the scatter σ gas and σ dyn : 1. Fixed (0.4 dex): We normalize our results such that σ gas = 0.39 dex and σ dyn = 0.07 dex, so that the relative contribution from gas inflows/cycling versus GMC formation/disruption are always fixed.This allows us to fix the scatter in log SFR at 0.4 dex, and isolate the impact that relative changes in timescales may have on SFHs and associated observables.We choose 0.4 dex since it is close to the commonly measured value for the scatter in the SFR-M * correlation (Kurczynski et al. 2016;Iyer et al. 2018;Leja et al. 2021).

Variable (TFC20):
We normalize our results to the values reported in TFC20 (see their Figure 9) of σ MWA = 0.17 dex, σ Dwarf = 0.53 dex, σ Noon = 0.24 dex and σ High−z = 0.27 dex.This involves relative changes in both the overall scatter and the relative contributions from gas inflows/cycling versus GMC formation/disruption.
The general behaviour of each case for fixed and variable scatter is highlighted in Figures 3 and 4, respectively.

GAUSSIAN PROCESS IMPLEMENTATION
We now describe how we can use the ACFs computed in §2 to quickly generate realizations of galaxy star formation histories (SFHs).In §3.1, we provide a brief overview of Gaussian Processes (GPs).In §3.2, we describe how we use them to generate realizations of synthetic galaxy spectra.Our implementation is publicly available at https: //github.com/kartheikiyer/GP-SFHand summarized in Figure 1.

Brief Overview of Gaussian Processes
A Gaussian Process (GP) is a generalization of the Gaussian distribution to the space of functions (Rasmussen & Williams 2006).Similar to our definition of a stochastic process, this means that our GP can generate an infinite set of values y n ≡ {y t } t=n t=1 at any time t = 1, . . ., n whose joint probability distribution will always follow a multidimensional Gaussian distribution where µ n is the n vector of mean values generated from some mean function µ(t) at times t = 1, . . ., n and C n,n is the n × n covariance matrix evaluated at each pair of times t = 1, . . ., n and t = 1, . . ., n.If we have some values y m that are known and want to predict a set of new possible values at our given times t = 1, . . ., n, this can be done by exploiting the fact that the joint distribution is which gives a conditional distribution of Together, these properties along with the functional nature of GPs are often summarized using the following notation: where the ∼ indicates "is a realization of" rather than the "is of the same order-of-magnitude as" definition usually used in the astrophysics literature.
Taken together, the above results imply that we can use GPs to quickly and easily generate realizations of our data y n either completely from scratch or conditioning on some known values y m based on some input mean µ(t) and covariance C(t, t ) functions that we can easily replace with any of the ACFs derived in §2.In particular, switching over to our variable of interest (ln SFR) and the model of interest (the Extended Regulator model discussed in §2.2), the model we explore in our paper takes the form ln SFR(t) ∼ GP(ln SFR base (t), C ExReg (τ )) where ln SFR base (t) is some baseline SFH we are interested in studying and again τ = t − t .In other words, any given realization of the SFH will depend on both the "baseline" (mean) SFH, ln SFR base (t), as well as the particular ACF C ExReg (τ ) defined by the Extended Regulator model.In practice, our GP is implemented following a similar procedure to Iyer et al. ( 2019) using a multidimensional Gaussian prior initialized at every point of an input time array.This is done through an instance of the GP_SFH() class, which is initialized with a user-determined kernel at a particular redshift, along with an astropy.cosmology()object and a fsps.stellarpopulation()object for generating spectra and other observables.Upon initialization, the instance computes the covariance matrix specified by the ACF at a range of time values ranging from 0 to the t univ at the specified redshift.Once this matrix is computed, realisations of SFHs can be generated simply by sampling a multivariate normal distribution at each time in the array with the covariance structure determined by the kernel.This can then be conditioned on observable constraints using eqn.17.

Implementation with Stellar Population Synthesis Models
To generate spectra corresponding to draws from the GP, we pass the star formation histories though the Flexible Stellar Population Synthesis (FSPS) code (Conroy et al. 2009;Conroy & Gunn 2010).To highlight the differences in galaxy spectra that arise solely from differences in the ACF, we choose a simple set of modeling assumptions (listed in Table 1) when generating spectra while keeping everything else fixed to their default FSPS values.We also demonstrate the effects of varying some of these parameters on our observables of interest in Section 5.1.
In practice, changes to the stellar population parameters can be made simply by reassigning the input fsps.stellarpopulationobject linked to the GP-SFH class instance.This modular implementation allows for a pre-computed covariance matrix to be rapidly associated with many different stellar population parameters while modeling and fitting SEDs, since that is the rate-limiting step to generating SFH realisations.For spectral features, we choose to consider features sensitive to star formation on a range of timescales (Kauffmann et al. 2003).Hα emission from O-and B-stars probes star formation on 4 − 10 Myr timescales (Flores Velázquez et al. 2021;Tacchella et al. 2022a).Hδ absorption traces star formation over the last 0.1 − 1 Gyr (Worthey & Ottaviani 1997;Wang & Lilly 2020b).Finally, the 4000 Å break strength (D n (4000)) provides a reliable tracer of the median age of the stellar populations that make up a galaxy, as can be seen in Figure 2 of Kauffmann et al. (2003), who point out that these indices are largely insensitive to dust attenuation effects, and the distribution of galaxies in this space is sensitive to stochasticity in star formation over the most recent ∼ 2 Gyr in a galaxy's past.This happens because galaxies that form stars at a steady rate tend to occupy a very narrow locus in Hδ EW -D n (4000) space.In addition to these, we also consider the equivalent width of the Ca H, K lines at λ 3933.6,3968.5 Å.The Ca II K line traces older stellar populations, while the Ca H absorption line is blended with H and [Ne III] and effectively traces intermediate age populations (Mayya et al. 2004;Zhu & Ménard 2013).The resulting equivalent widths probe a range of timescales as seen in Figure 5 and Appendix B.
For our spectra, we get the Hα line luminosity directly from the FSPS outputs and following a procedure similar to Kauffmann et al. (2003), we use the 3850−3950 and 4000−4100 Å continuum bands introduced by Balogh et al. (1999) to compute the strength of the D n (4000) break, and compute the Hδ EW using the bandpasses of 4083.50 − 4122.25 Å, 4041.60 − 4079.75Å, and 4128.50 − 4161.00Å for the index and blue/red continuum bands respectively, defined in Table 1 of Worthey & Ottaviani (1997).For the Ca-K EW , we use 3929.51− 3941.22Å, 3907.01 − 3929.51Å, and 3941.22 − 3961.02Å and for Ca-H EW , we use 3961.02− 3980.83Å, 3941.22 − 3961.02Åand 3980.82 − 3997.03Å for the index, blue and red bandpasses respectively.For better comparison across the different stellar masses that could be produced due to bursts and troughs in individual realisations of SFHs, we normalize these quantities by the stellar masses of each realisation, effectively reporting e.g., the distribution of Hα luminosity (in L ) per solar mass for the different cases discussed in Section 4. The Hδ EW , D n (4000) and broad-band galaxy colors remain unaffected by this normalisation.

SPECTROPHOTOMETRIC SIGNALS OF CHANGING STOCHASTIC BEHAVIOR
Following the procedure highlighted in Figure 1, in this section we identify the particular spectrophotometric signals that can help to distinguish different types of stochastic behavior (i.e.varying correlation timescales τ in , τ eq , and τ dyn and fluctuation strengths σ gas and σ dyn ).We generate 10000 realizations of various SFHs (log SFR(t)) for each of the cases outlined in §2.3 with the scatter fixed (Figure 3) and the scatter matched to TFC20 (Figure 4).We then feed these into FSPS to generate a set of UV-to-IR galaxy SEDs as described in Section 3.2.
To highlight the behavior of our model, we first investigate overall effects that the parameters in the Extended Regulator model may have on a few key observables.Our results are highlighted in Figure 5, where we vary single parameters in the Extended Regulator model while holding the rest fixed at fiducial values (σ gas = 1.0, τ eq = 500 Myr, τ eq = 150 Myr, σ dyn = 0.1, τ dyn = 10 Myr).We find that the exact mechanism for adjusting the "burstiness" of star formation -whether through the overall level of variability in the gas inflows/cycling (σ gas ) or GMC formation/disruption (σ dyn ), or through the duration of the (gas equilibrium) correlation time τ eq -leaves different imprints on various observables even at fixed SFR scatter.In particular, while varying σ gas affects both the stellar mass and (s)SFR distributions, varying τ eq only affects the stellar mass distribution.This gives at least one way to distinguish populations with differing amounts of variability about the same base set of SFHs.
Since changing correlation timescales and the level of scatter can lead to large differences in the final stellar mass formed, we choose to normalize all SEDs based on their final stellar mass before investigating possible (relative)   3 with fixed scatter.Using 10000 SFH realizations for each case, we highlight differences in the median stellar mass-normalized full UV-to-IR galaxy SEDs with respect to the reference SED of a galaxy with the base SFH (i.e., constant SFR of 1M /yr).As in Figure 5, increasing the correlation timescales tends to decrease the final stellar mass and increase the variability due to longer bursts/quiescent periods above/below the base SFH.On a broad scale, this results in a decreased relative sSFR for models with more variability, leading to signatures in the rest-UV.More subtle signatures of the stochasticity are also visible in the relative strengths of the 4000 Å break and absorption lines.
differences.This helps to highlights trends as a function of specific star formation rate (sSFR) rather than just SFR, and helps to account for the increasing (expected) variance in the total stellar mass formed for more bursty SFHs.

Fixed (0.4 dex) scatter
The results of this procedure assuming a fixed (0.4 dex) scatter are the SFR is shown in Figures 6 and 7 for all four of the models we consider.Figure 6 shows the median spectra for each model compared to a reference spectrum corresponding to the base SFH, and Figure 7 shows a corner plot comparing the full distributions of the spectral features described in Section 3.2.As expected based on our results in Figure 5, since the gas inflow/cycling physics dominates the main behavior of the model, decreasing the associated gas timescale (i.e.τ gas = τ in + τ eq ) to make SFHs more bursty leads to both larger stellar masses and a tighter overall distribution for a given scatter.While the massnormalized spectra have very similar medians across the four archetypes, they diverge in the optical and ultraviolet.This, combined with subtle differences in the Hδ EW (rising with the decreasing τ gas as SFHs form more of their stars in recent bursts) and the slope of the Hα-D n (4000) correlation (becoming shallower and more dispersed as SFHs become less correlated) indicate that, even in the case where the only thing varying are the timescales, constraining them should be possible with a large enough sample size.
Since obtaining spectroscopic data for large ensembles of galaxies is expensive and may not always be feasible, we also consider distributions of broad-band colors corresponding to the different archetypes.We consider three different colorcolor spaces -(i) the commonly used UVJ diagram (Wuyts et al. 2007;Williams et al. 2009;Muzzin et al. 2013), where the different archetypes can be differentiated based on their sSFR distributions, (ii) the NUV-r-K diagram (Arnouts et al. 2013;Moutard et al. 2016), which traces the differences in stellar mass distributions since the K S band probes the rest-frame 1.6µ feature, and the N U V − r probes the effects of dust and SFR, and (iii) the recently introduced  3 with fixed scatter.Using 10000 SFH realizations for each case, we highlight differences in the distributions of spectral features sensitive to star formation across a range of timescales, including Hα, Hδ, and Dn(4000).While the differences in the median value of observables are subtle, the joint distributions do show some differences that arise from changes in the underlying sSFR distribution.
wide-baseline FUV-V-(Wise)W3 diagram following Leja et al. (2019b), which is more sensitive to lower sSFRs than the UVJ diagram, where galaxies with low sSFRs tend to populate the top-left portion of the space.As shown in Figure 8, all three color-color spaces provide a means to differentiate between the different regimes of stochasticity typified by the four archetypes with sufficient sample sizes and SNR.In addition to this, upcoming observations with JWST will help push rest-frame colors out to higher redshifts, and provide extremely high SNR probes of differences in the stochasticity across the four archetypes considered here.

Variable (TCF20) scatter
In the more realistic case from TFC20 where the scatter among each archetype also varies, we see larger differences in the corresponding observables (Figures 9 and 10).To first order, this is associated with the large variations in σ gas , leading to shifts and/or broadening of the stellar mass and SFR distributions that propagate down to Hα, Hδ EW , and D n (4000).These can be seen most clear by looking at the distribution of values for the dwarf model (which has the largest σ gas ) and the MW analogue (which has the longest τ gas ).Perhaps the biggest change in terms of the spectral features is the Hδ EW distribution which, while still varying in width, is much more centered around a small mean amount of emission from recent SFR.This mainly happens due to the reduced scatter in all the models, which provides less opportunity for the bursty episodes of SF to form large amounts of stellar mass and thereby reduce the sSFRs on average.In terms of the broad-band colors, changing the scatter has the effect of scaling the corresponding distribution of colors by a proportional amount, while maintaining the differences in shape due to varying stochasticity.

DISCUSSION
Having established the GP-SFH formalism and demonstrated sensitivity to the extended regulator model parameters in Section 2, we consider the implications and predictions we can make with JWST observations using this formalism in 5.1.We consider the effects of other stellar population variables such as dust and metallicity in Section 5.2.Sec 5.3 demonstrates how the GP-SFH can be used as a stochasticity prior for binned SFHs, and 5.4 considers how populations of galaxies with specific spectral features can provide novel constraints on stochasticity.Sec 5.5 considers relaxing the assumption of a fixed base SFH, and Sec 5.6 considers the assumption of stationarity and shows an example of relaxing it with a time-varying kernel and 5.7 considers some caveats and challenges of the GP-SFH formalism.

Implications for JWST
The release of JWST data from the early release observations (ERO, aka "Webb's First Deep Field"; Pontoppidan et al. 2022), GLASS (Merlin et al. 2022) and CEERS (Finkelstein et al. 2022) have demonstrated the incredible potential for probing star formation in galaxies across an incredible range of redshifts, environments and stellar masses.This also enables studies of star formation stochasticity using the colors and spectral features in this work to constrain the ExReg model parameters -i.e.stochasticity amplitudes and timescales.
At redshifts around cosmic noon, we will see a major improvement in being able to directly measure rest-frame colors.In contrast to HST ACS+WFC3, which can measure rest-frame UVJ colors at 0.21 z 0.29, JWST's NIRCam filters will extend this to 1.51 z 2.55, and a combination of HST+JWST will span the entire 0.21 z 2.55 range.This is also similar to the redshift range that slitless spectroscopy using NIRISS will be able to measure the spectral features discussed in this work for large populations of galaxies (Willott et al. 2022).
At higher redshifts, as we attempt the challenging task of measuring the star formation rates and histories of these galaxies, care must be taken to account for the dependence of the results on the assumed priors (Tacchella et al. 2022b;Whitler et al. 2022), and use a combination of spectroscopic and photometric data where available (Tacchella et al.  6, but now for the four cases highlighted in Figure 4 with variable (TFC20) scatter.The variation in scatter dramatically expands the set of sSFR distributions due to the changing distribution of total stellar mass formed after a fixed time given a constant SFH.These differences are most prominent when comparing the model with the smallest scatter and longest timescales (MW analogue; red) with the one with the largest scatter (dwarf; orange).2022c).The GP-SFH formalism described here can help motivate SFH priors for the next generation of SED fitting based studies based on estimates of the stochasticity from lower-redshift analogs or simulations.
Additionally, although we did not consider its effects in this work, the evolution of gas-phase and stellar metallicity in galaxies is also tied to their star formation, and further studies of their correlated properties (as in Camps-Fariña et al. 2022;Zhou et al. 2022) could provide further observable tests of the ExReg model timescales, and priors for SED fitting codes that explicitly allow for evolution in chemical enrichment over time (Pacifici et al. 2013;Thorne et al. 2021).

Effects of varying other SED parameters
SED modeling depends on a host of assumptions about the stellar populations that make up a galaxy, in addition to dust attenuation and emission from dust heating, nebular regions and AGN.In our current analysis we have held most of these constant in order to isolate and study the effects of perturbing the SFH stochasticity model, but it is informative to consider the extent to which varying these additional parameters will broaden the distributions we expect to observe.
Figure 11 shows the effect of varying the stellar metallicity and the dust attenuation (assuming a Calzetti dust law; Calzetti et al. 2000) on the spectral indices we consider.This analysis is done for a galaxy with a fixed SFH of 1M /yr and other parameters corresponding to Table 1.The effects of a distribution of values in either of these stellar population parameters would correspond to a broadening in the distribution of spectral indices by an amount proportional to the mean and width of the dust/metallicity distribution.For example, a distribution of ∼ 1 dex in metallicity centered around solar metallicity would correspond to a spread of ∼ 0.04 in D n (4000) and ∼ 0.3 dex in log Hα luminosity.While convolving the distributions in Figure 9 does make the different models harder to discriminate between, it is still distinct enough to be possible with a large enough sample size.This is additionally helped by the fact that the broadening of distributions in the spectral indices is not homogeneous, and in fact displays quite different  7, but now for the four cases highlighted in Figure 4 with variable (TFC20) scatter.The variation in scatter dramatically expands the set of sSFR distributions due to the changing distribution of total stellar mass formed after a fixed time given a constant SFH.These differences lead to greater distinguishing power using the distributions of the spectral features we consider.
signatures across the three indices for dust and metallicity -notably that dust attenuation does not affect the Hδ EW .We have not shown the effects of varying SPS models or the IMF, since that would correspond to an overall shift in the indices rather than a broadening of the distribution.In the rest-optical part of the SED that we study in this work, we are also not significantly affected by AGN, dust re-emission and other factors that manifest in the mid-to-far IR portions of the SED.These effects have also been studied in relation to SFR stochasticity in the literature (Broussard et al. 2019b;Emami et al. 2019b While we expect the distribution of these spectral features to be broadened due to variations in these quantities, the differences are still sufficient to tell the different toy models in Figure 10 apart.
One additional factor to note is that the broadening predicted by Figure 11 assumes that variations in dust and metallicity are independent of SFR stochasticity and star formation history.However, given that stochasticity determines the frequency of sharp bursts of star formation, it is likely that it will be correlated with the chemical enrichment of the galaxy.Although this is outside the scope of this work, cosmological simulations of galaxy evolution could shed light on the link between these parameters and help develop correlated priors for use with future observations.

A prior for binned SFHs
The Gaussian process implementation described in this paper can also be adapted as a prior for binned SFHs that are used by spectrophotometric fitting codes like Prospector (Leja et al. 2017;Johnson et al. 2021) or CMD based methods like Match (Dolphin 2002;Weisz et al. 2014).Prospector in particular uses either a Dirichlet prior or a continuity prior which parametrizes the log SFR ratios between bins using a Student's-t distribution.To incorporate the covariance models described in this work, it suffices to replace these priors with the covariance values C ExReg (t i , t j ) where t i , t j correspond to the centers of the individual bins.In Appendix D, we verify that this procedure yields SFHs identical to sampling from the high-resolution GP-SFH and degrading the resolution to match the logarithmically spaced time bins used in these codes.While we do not reproduce the distribution of observational metrics in Figures 6-9 using this formalism, we expect them to be very similar due to the lack of information encoded in galaxy SEDs about short-term variability at large lookback times.Figure 12 shows samples of binned SFHs corresponding to the four stochasticity regimes used in this work.5.4.How do we choose a 'population of galaxies' to study observationally?
An underlying assumption in comparing distributions of spectral features or colors is that galaxies in the sample belong to the same underlying population and thus can be described using the same model for star formation stochasticity.Caution should therefore be exercised when creating a galaxy sample to minimize contamination by other populations.Methods like clustering in colors, physical properties, or SFHs can be used to select galaxies that are likely to belong to the same population.Additional methods like selecting specific galaxy sub-populations using self-organising maps (Teimoorinia et al. 2022;Holwerda et al. 2022;Davidzon et al. 2022) or in the latent space of variational auto-encoders (Portillo et al. 2020) can also be used for this purpose.
It is currently not well understood whether properties of SF stochasticity (i.e.properties like σ gas or τ eq in the Extended Regulator model) correlate strongly with other physical properties like galaxy size, environment, morphology or dynamics.Since the variability in these quantities can span a wide range of timescales, that may or may not relate to the timescales on which SFR fluctuations are correlated.It would be interesting to study this further, since the stochasticity model can also influence the presence of certain galaxy populations.
For example, in Figure 8, we notice that the Dwarf and Cosmic noon populations show a slight excess of galaxies with (r − K < 0.75) and (N U V − r ∼ 3).This region is highlighted in the left panel of Figure 13.Since spectral sensitivity falls off as a function of age at different rates depending on the wavelengths under consideration, an assumed model for SFH stochasticity can produce unique spectral signatures depending on a combination of broadband filters.In the NUV-r-K color-color space for example, the r-K color has a mild linear dependence with age, except for a short period between ∼ 5 − 60Myr during which the color sharply decreases.In complement to this, the NUV-r color is relatively flat until ∼ 20Myr, after which it shows a linear dependence with age.Because of this combination of sensitivities, a portion of NUV-r-K color space (≈ (N U V − r > 2) & (r − K < 0.5)) is uniquely sensitive to galaxies that recently experienced a sudden recent rise and fall in their star formation histories9 .
Since the probability of such an event is directly proportional to the amount of burstiness and effective timescales over which SFR is correlated, the four stochasticity models considered above make differing predictions for the probability that a galaxy can have such an event (and therefore on the number of galaxies in a given sample).We examine this Figure 13.Post starburst galaxy SFHs are naturally produced by the different models, although the fraction of these galaxies with a recent burst on a ∼ 100 Myr timescale in any given sample is higher for the dwarf and cosmic noon stochasticity regimes.This is picked up by a color-color selection shown in the bottom panels, similar to Figure 8 where we select galaxies with R − Ks < 0.75 and N U V − R > 2. Median SFHs along with the 68% variance for these galaxies is shown in the top panel, where the dwarf and noon samples show a clear post-starburst feature in the last ∼ 1 Gyr.This is unlikely in the MWA case due to the much longer correlation timescales of SFHs, and in the High-z case because of the much shorter effective correlation timescales.Statistical samples of galaxies can therefore allow us to use the relative abundance of post-starburst galaxies to constrain the overall sample stochasticity.
portion of the NUV-r-K color space better in Figure 13, finding that the SFHs of galaxies with these colors tend to show a strong post-starburst feature in their SFHs, with the Dwarf and Cosmic noon ACFs resulting in a higher number of these galaxies compared to the MW or high-z populations.Additional quantities like the fraction of PSB galaxies and the timescales of the recent burst could therefore be useful tracers of SF stochasticity in future studies.

Varying the base SFHs
The analysis in Sections 4.1 and 4.2 assumes that the individual SFHs can be described as perturbations around a base SFH, which we assume to a constant SFH with SF R = 1M /yr.While it is possible to relax this assumption in our current framework, caution needs to be exercised when the variability of the base SFHs is comparable to that of the timescales in the extended regulator kernel, since this will modify the ACF by adding power on the longer timescales when the autocorrelation time of the base SFH is of the same order as that for the extended regulator model τ A,baseSFH τ A,ExReg .Appendix C presents a detailed discussion of the effects on spectral features with an implementation where the base SFHs themselves vary across the sample and are drawn from a distribution.While this leads to a broadening of the distributions of individual spectral features, we find that it is still possible to differentiate between the models by comparing distributions of spectral features.

The assumption of stationarity and ergodicity
For simplicity, the derivation in Section 2 assumes that the star formation histories of a population of galaxies are stationary and ergodic.The assumption of stationarity requires that the PSD or ACF of a galaxy SFH does not have an explicit time dependence.However, it is not necessary that SFHs in the real universe follow this, with either the stochasticity or timescales of the PSD model evolving with time.However, (i) for most science cases that discriminate between different models of stochasticity, the evolution is slow enough that this assumption is expected to hold (see the discussion in §3.2 of Wang & Lilly 2020a), and (ii) if/when we decide to relax this assumption of stationarity, the kernel in our Gaussian process formalism can be updated to account for that.Indeed, non-stationary kernels are an open topic of research in Gaussian processes (Rhode 2020) and models for the time-evolution of the ACF are an important part of the future work enabled by this formalism.
As observational data with future telescopes unlocks new timescales and large populations of galaxies across different cosmic epochs, our models can be updated to include variations in the extended regulator model parameters as a function of time.In this case, it is possible to relax the assumption of stationarity (i.e., Eqn. 2) and implement a more general Gaussian process for galaxy (log) SFRs, as shown in Figure 14.A simple extension in this direction would be to make one or more of the parameters in the extended regulator model time-dependent.As an example of this case, we consider a simple time-dependence to both σ gas and τ eq given by This form is chosen for the σ to allow for increased variability at higher redshifts, while the increasing τ eq is linked to the increasing dynamical times of galaxies with decreasing redshifts.Ergodicity is a more subtle issue, and concerns the fact that we perform our analysis using populations of galaxies.For a dynamical system, ergodicity implies that the variability of a single galaxy's SFH over time is equivalent to the variability of an ensemble of SFHs observed at a given epoch, after accounting for selection functions and completeness.This is considered in Sec 3.2 of Wang & Lilly (2020a), which also considers the need for ergodicity while working with galaxy PSDs.While it is unlikely that galaxy SFHs are fully ergodic, due to the changing conditions in which galaxies form stars at different epochs, the extent to which this is violated is expected to be minor and can be tested in the future using cosmological simulations.

Challenges and caveats
Techniques studying the stochasticity of star formation across timescales rely on a host of modeling assumptions, and these should be kept in mind while using inference to determine the SFH kernel parameters.This paper aims to provide a framework for forward modeling observables that can be combined with realistic distributions in other properties like the mean SFH model, stellar and gas-phase metallicity and dust attenuation, along with assumptions of stellar tracks, isochrones, IMF, dust law, binary fraction, TP-AGB contributions, and other factors that could systematically bias results in any SED-based analysis.This is mitigated somewhat by the fact that the code developed here is modular and runs off FSPS.As a result, the SFH stochasticity model can be folded into existing SED fitting codes like dense basis (Iyer et al. 2019) and Prospector (Leja et al. 2017;Johnson et al. 2021), or combined with realistic models for physical properties like dust (Hahn et al. 2022;Nagaraj et al. 2022) in a population-level forward modeling setup.
Another concern is that the assumed form for the GP covariance does not take into account all the possible physical processes that induce variability in galaxy SFRs.While this is undoubtedly true, the processes modeled here are 'effective' timescales as a result of a wide range of physical processes -τ in can include changes in the gas reservoir due to pristine inflows, outflows, mergers, and the re-integration of gas blown out in previous cycles; τ eq accounts for the fact that multiphase gas undergoes complex transformations across a galaxy as a result of stellar and AGN feedback, and τ dyn accounts for dynamical processes that create and destroy GMCs, which can be due to stellar feedback but also gas compression by spiral arms or compaction-induced starbursts.The next generation of large-volume, high-resolution simulations (Pillepich et al. 2019;Nelson et al. 2019;Feldmann et al. 2022) and studies that systematically vary feedback prescriptions (Villaescusa-Navarro et al. 2021, 2022) are needed to robustly understand these effective timescales, and further refine the GP kernel as needed.
A minor concern is the scalability of the GP-formalism in the limit of high-resolution SFHs and/or computationally intensive forward-modeling comparisons to large datasets of observations.However, fast GP-based methods that include automatic differentiation (e.g.tinyGP; Foreman-Mackey et al. 2022) can provide workarounds to these issues as they arise.

CONCLUSION
Modeling the stochasticity of star formation over a range of timescales provides a way to connect observations of galaxy star formation histories to the underlying physical processes driving galaxy growth and quenching.
In this paper, we propose a fast, modular, Gaussian process-based (GP) formalism implementing the extended regulator model based on TFC20 as an autocovariance function (ACF) or kernel.This model provides a parametric form for SFR stochasticity as a combination of different physical processes, and is completely characterised by three effective timescales corresponding to stochastic gas inflows, equilibrium and dynamical processes influencing GMC creation and destruction.
Implementing a GP with this kernel allows us to make extremely fast draws of galaxy SFHs with a particular SFH autocovariance structure and to use them to forward-model galaxy spectra and dependent observables.Studying these observables as a function of the kernel parameters allows us to quantify differences as a function of the extended regulator model's timescales and thus differentiate between different regimes of stochasticity.This is illustrated by considering four toy-models for galaxy populations -Milky Way analogs, dwarf galaxies, massive galaxies at cosmic noon, and galaxies at high redshifts.We model the spectra of these galaxies using FSPS and study distributions of spectral features including Hα and UV-based SFR indicators, Hδ and Ca-H,K absorption line equivalent widths, the D n (4000) spectral break and broadband galaxy colors, finding that these distributions are sensitive to the Extended Regulator model parameters, and their distributions and covariances can be used to discriminate between the models.
Since increasing the amount of stochasticity leads to greater stellar masses formed in intense bursts of star formation, the sSFR distribution, and therefore the flux in SFR tracers like Hα and rest-U V , are sensitive to the overall level of stochasticity.Complementary to this, the Hδ EW vs D n (4000) space traces star formation over longer timescales, and is extremely sensitive to both the overall level of variability as well as the timescales on which SFR is correlated.We also find that the rest-frame broadband colors reveal populations of galaxies such as post-starbursts that are preferentially found in models that allow SFR correlations on the timescales that the colors are sensitive to, and can thus be used as additional constraints on the ExReg model parameters.
The GP-SFH formalism can also be easily incorporated into existing SED fitting codes to provide realistic priors for SFH covariance or infer them from future high S/N spectrophotometric observations using JWST, used to study the effective timescales in cosmological simulations, and further expanded to include factors like non-stationarity.Code to reproduce our results can be found at kartheikiyer/GP-SFH.   (Caswell et al. 2019), scipy (Virtanen et al. 2020), numpy (Walt et al. 2011), corner (Foreman-Mackey 2016), FSPS (Johnson et al. 2021), dense basis (Iyer et al. 2019) and hickle (Price et al. 2018).
Figure 16.Implementing the ExReg kernel with the realistic SFHs described in Iyer et al. (2019).The top panel shows draws from a Dirichlet distribution, with perturbations from the various toy models added on to the base SFHs.The bottom panels follow the same format as Figure 6, and show that despite the increased scatter due to SFH variations, the differences between the models are still distinguishable.
the underlying population itself possesses some intrinsic variability in their base SFHs around some population mean, we still observe some differences in Figure 17 in the joint distribution of observables, noticeably Hα and Hδ.The former difference comes from the varying sSFR distributions due to the perturbations of the different stochasticity archetypes, and the latter due to the varying amount of SFH burstiness in the last ∼ 1 − 2 Gyr, with the median of the distribution rising with decreasing τ gas .

D. VERIFYING THE COVARIANCE FOR BINNED SFHS
In practice, many non-parametric SFH codes use SFHs binned in roughly logarithmic bins in lookback time with varying priors on stochasticity.In this Section, we verify that directly sampling from the covariance function using bin-centers is equivalent to sampling high-resolution SFHs and degrading them to the same coarse time-bins.
We start with 1000 samples of high-resolution SFHs from a GP corresponding to each of the four galaxy regimes described in §2.We then degrade them to binned SFHs with 10 equally spaced bins in log lookback time such that the SFR in a given bin is the average of the SFR in that interval, as shown in the middle panels of Figure 18.As expected, the covariance matrix computed from the binned SFHs matches the analytical estimate, with small differences due to spectral leakage from the finite length of the time-series.Repeating this analysis now drawing from the SFHs with the coarse time array corresponding to the bin centers with the same kernels confirms that the coarse SFHs are consistent with the distribution obtained by coarsening the high-resolution SFHs.This also results in a significantly faster GP  2019).The top panel shows draws from a Dirichlet distribution, with perturbations from the various toy models added on to the base SFHs.The bottom panels follow the same format as Figure 7, and show that despite the increased scatter due to SFH variations, the differences between the models are still distinguishable.
that may be more suitable for forward-modeling large ensembles of observations that require repeated computation of the covariance matrix.

Figure 1 .
Figure1.An illustration of the overall strategy used in this paper to test the impact that stochasticity in star formation histories (SFHs) can have on galaxy spectral energy distributions (SEDs).We start with power spectral density (PSD; top left) of a dwarf-like galaxy population defined within the framework of the Extended Regulator model of TFC20 (see §2.3), with individual components highlighted.We then convert this PSD to the corresponding auto-covariance function (ACF; top right).Both the PSD and ACF highlight that variability on short timescales is dominated by giant molecular clouds (GMCs) while variability on longer timescales is dominated by a combination of gas inflows and cycling between atomic and molecular gas.Using the ACF, we define a Gaussian Process (GP; bottom left) to quickly and easily sample many realizations of the log SFR(t) over time relative to some mean log SFR base (t), which here we take to be 1M /yr.We again highlight both long-term and short-term (see inset) variability driven by gas inflows/cycling and dynamical processes, respectively.Finally, we generate SEDs at z ∼ 0 by running the GP-sampled SFHs through a stellar population synthesis model (bottom right).The variability in the SFHs leads to differences in the overall normalization (from differences in total stellar mass formed) as well as particular spectral features (from differences in relative contributions of various stellar populations).Properties apart from the SFH assumed while generating the spectra are detailed in Section 3.2, and are held constant to highlight differences from the SFHs.Understanding how the distributions of spectral features corresponding to various galaxy populations correspond to the properties of their PSDs will allow us to put constraints on the timescales on which gas inflows, baryon cycling, and GMC physics drive galaxy evolution.

Figure 2 .
Figure2.The effects of varying the gas inflow timescale τin and equilibrium gas cycling time τeq on galaxy star formation histories.We highlight four ACFs: two damped random walks with long and short correlation timescales τeq, respectively, and two Regulator model realizations (see §A.5) corresponding to a Milky Way analogue and a dwarf galaxy (see §2.3).To be comparable across different timescales, the four ACFs have all been scaled such that C(τ = 0) ≈ 1.0 dex (left panel ).We then generate a single realisation of the SFH log SFR(t) using a GP with each of the four ACFs with a fixed random number seed to highlight differences that arise from varying the ACF (right panel ).ACFs with longer-timescale correlations tend to spend larger stretches of time above and below the mean SFH, while short-timescale kernels tend to show more "bursty" behaviour that generate larger (but short-lived) excursions from the mean.The Milky Way analogue model in particular includes correlations at both long and short timescales from gas cycling and inflows, respectively, giving it the smoothest, most correlated SFH.

Figure 3 .
Figure3.An illustration of the ACF (left panel ) and several corresponding realizations of SFHs (right panels) for the four cases discussed in §2.3: Milky Way analogues at z = 0.1 (MWA; red), dwarf galaxies at z = 0.1 (orange), massive galaxies at cosmic noon at z = 2 ("noon"; green), and galaxies at high redshifts (high-z; blue).SFHs are generated assuming a base SFR of 1M /yr for the age of the universe (i.e.ignoring formation times) and with a fixed total scatter of 1 dex and a fixed relative contribution of GMC formation/disruption to the variance of f dyn = 0.03.These highlight the relative changes in behavior that arise only from varying correlation timescales τin, τeq, and τ dyn .Within each panel, SFHs are plotted with various intensities to help the eye distinguish them.

Figure 4 .
Figure 4. Similar to Figure3, but the overall scatter for each case is set to the values used in TFC20.These highlight overall changes in behavior that arise both from varying correlation timescales as well as changes in the overall burstiness of SFHs.

Figure 5 .
Figure 5. Effects of varying the Extended Regulator model parameters on the distributions of stellar mass and spectral features spanning a range of timescales: (in increasing order) Hα & NUV fluxes, equivalent widths of Hδ & Ca-H,K, u-band flux, and finally the strength of the Dn(4000) Å break.Histograms at the bottom of each panel show the fiducial distribution, corresponding to (σgas, σ dyn , τin, τeq, τ dyn ) ≡ (1.0, 0.1, 500, 150, 10), while the individual bars show the change in the median (thick solid line) and the width of the distribution, shown using the 16-84 th percentiles (thin error bars) upon changing each of the ExReg model parameters.The non-degenerate changes in the observables upon perturbing the ExReg model allow us to use a combination of spectral features sensitive to a range of timescales to test and constrain the model parameters observationally.

Figure 6 .
Figure6.Observational signatures caused by only varying correlation timescales τin, τeq, and τ dyn based on the four cases highlighted in Figure3with fixed scatter.Using 10000 SFH realizations for each case, we highlight differences in the median stellar mass-normalized full UV-to-IR galaxy SEDs with respect to the reference SED of a galaxy with the base SFH (i.e., constant SFR of 1M /yr).As in Figure5, increasing the correlation timescales tends to decrease the final stellar mass and increase the variability due to longer bursts/quiescent periods above/below the base SFH.On a broad scale, this results in a decreased relative sSFR for models with more variability, leading to signatures in the rest-UV.More subtle signatures of the stochasticity are also visible in the relative strengths of the 4000 Å break and absorption lines.

Figure 7 .
Figure 7. Observational signatures caused by only varying correlation timescales τin, τeq, and τ dyn based on the four cases highlighted in Figure3with fixed scatter.Using 10000 SFH realizations for each case, we highlight differences in the distributions of spectral features sensitive to star formation across a range of timescales, including Hα, Hδ, and Dn(4000).While the differences in the median value of observables are subtle, the joint distributions do show some differences that arise from changes in the underlying sSFR distribution.

Figure 8 .
Figure8.The distribution of galaxy colors corresponding to the populations described in Figure6in U-V-J (left), NUV-r-K (middle), and NUV-V-W3 (right).While the broad distributions are similar among all populations due to their similar SFHs, we see meaningful shifts in both the wings of the distributions and the overall general centers.

Figure 9 .
Figure9.As in Figure6, but now for the four cases highlighted in Figure4with variable (TFC20) scatter.The variation in scatter dramatically expands the set of sSFR distributions due to the changing distribution of total stellar mass formed after a fixed time given a constant SFH.These differences are most prominent when comparing the model with the smallest scatter and longest timescales (MW analogue; red) with the one with the largest scatter (dwarf; orange).

Figure 10 .
Figure10.As in Figure7, but now for the four cases highlighted in Figure4with variable (TFC20) scatter.The variation in scatter dramatically expands the set of sSFR distributions due to the changing distribution of total stellar mass formed after a fixed time given a constant SFH.These differences lead to greater distinguishing power using the distributions of the spectral features we consider.

Figure 11 .
Figure11.The impact that varying dust and metallicity has on the spectral features shown in Figure7and 10.To highlight the size of the effect, we keep the axis ranges in the corner identical while highlighting the variability in the top inset panels.While we expect the distribution of these spectral features to be broadened due to variations in these quantities, the differences are still sufficient to tell the different toy models in Figure10apart.

Figure 12 .
Figure 12.Log-binned SFHs sampled from Gaussian processes using the Extended Regulator model kernel, corresponding to the four stochasticity regimes described in Section 2.

Figure 14 .
Figure 14.Implementing a non-stationary kernel in the GP-SFH formalism.While following most of the behaviour of the MW analogs, this kernel has a time-varying component to (i) σgas, which rises to its maximum value at z ∼ 1, after which it falls again, and (ii) τeq, which keeps rising with time, leading to smoother, less 'bursty' SFHs at lower redshifts.(Left) The time-dependent 'kernel' for the GP with the dotted lines showing a deviation of 1 Gyr to the future and past, and (right) the median behaviour of the galaxy sample (black line and shaded region) and individual realisations (red lines) of galaxy SFHs.
JSS & KGI jointly led the work and drafted the majority of the text.NC to ST contributed ideas and helped to draft and edit the text.The authors are grateful to Harry Ferguson, Camilla Pacifici, and Vince Estrada-Carpenter for their comments and suggestions.JSS is grateful to Rebecca Bleich for continuing to tolerate his excessive patronage of McDonald's.KGI is grateful to the nonzero intersection in the Venn diagram of Noopur Gosavi, 5e, sci-fi and potatoes.The Dunlap Institute is funded through an endowment established by the David Dunlap family and the University of Toronto.EG acknowledges support from NASA ADAP grant 80NSSC22K0487.

Figure 17 .
Figure17.Implementing the ExReg kernel with the realistic SFHs described inIyer et al. (2019).The top panel shows draws from a Dirichlet distribution, with perturbations from the various toy models added on to the base SFHs.The bottom panels follow the same format as Figure7, and show that despite the increased scatter due to SFH variations, the differences between the models are still distinguishable.

Figure 18 .
Figure 18.Top: Covariance function calculated for logarithmically spaced time bins for the four cases considered in this paper.Middle: Individual high-resolution SFHs degraded to the 10 log-spaced bins in lookback time.Bottom: Comparison of the computed correlation function to the analytic ACF for the four cases.Differences are due to shot-noise from a limited sample size and finite length for the time-series.

Table 1 .
Modeling assumptions for generating spectra corresponding to different ACF cases.