The Diversity and Variability of Star Formation Histories in Models of Galaxy Evolution2020The Diversity and Variability of Star Formation Histories in Models of Galaxy Evolution–2

We consider the star formation histories from a wide range of galaxy evolution models, ranging from hydrodynamical simulations (Illustris, IllustrisTNG, Mufasa, Simba, EAGLE), a semi-analytic model (Santa-Cruz SAM), an empirical model tuned to match observations across a range of observations (UniverseMachine), and three suites of zoom simulations (FIRE-2, g14 and Marvel/Justice League ) with a higher resolution and more explicit prescriptions for the interstellar medium (ISM) and stellar feedback (see reviews by Somerville & Davé 2015; Vogelsberger et al. 2020, for a summary of the individual components of these various models).

For simplicity, in the current analysis we limit ourselves to (i) considering only a fiducial run from each model, since some models have multiple runs varying the parameters of various sub-grid recipes, (ii) considering the SFHs of only central galaxies above a stellar mass threshold of $10^{9}$M_⊙ to partially mitigate resolution effects, and (iii) studying galaxies at $z\sim 0$, with model variants and redshift evolution to be considered in further work. Figure 2 shows the normalised distributions of stellar mass for galaxies from each model at $z\sim 0$, used in the current analysis. The FIRE-2, g14 and Marvel/JL zoom simulations have a much smaller sample of 14, 8 and 5 galaxies, respectively, spanning a range of stellar masses from $\sim 10^{9}$M_⊙ to $\sim 10^{11.5}$M_⊙. While these zoom simulations allow us to probe SFH fluctuations to shorter timescales compared to the large-volume models due to their much finer spatiotemporal resolution, which allows them to resolve GMC-scale structures and treat feedback more explicitly, these galaxies are not representative of a cosmological sample. Caution should therefore be employed in generalizing trends in their variability.

Each model of galaxy evolution is described briefly below, with references to relevant papers containing more detailed descriptions. Since the current analysis deals with galaxy SFHs, the descriptions focus on how each simulation implements star formation and feedback, and Table 1 contains a summary of the resolution, box size and number of galaxies from each simulation.

• •

  Illustris (Vogelsberger et al. 2014b; Genel et al. 2014)

  The Illustris project²² 2 https://www.illustris-project.org/ is a large-scale hydrodynamical simulation of galaxy formation using the moving mesh code AREPO (Springel 2010). The model includes recipes for primordial and metal-line cooling, stellar evolution and feedback, gas recycling, chemical enrichment, supermassive black hole (BH) growth and AGN feedback (Springel & Hernquist 2003; Vogelsberger et al. 2013). Given the spatial resolution of $\simeq 1$ kpc, giant molecular clouds are not resolved. A sub-resolution model for an effective equation-of-state (Springel & Hernquist 2003) is implemented where a star particle is stochastically produced above the critical hydrogen number density of $n_{\mathrm{SF}}=0.13\mathrm{cm}^{-3}$ on a density-dependent timescale that reproduces the observed Kennicutt-Schmidt relation (Schmidt 1959; Kennicutt 1989). Star formation results in supernovae, which release kinetic winds that expel gas from their surroundings and chemically enrich the ISM. These winds are implemented by launching hydrodynamically decoupled ‘wind particles’ that recouple to the gas when they leave the dense local ISM and reach a cell with a density $<0.05n_{\mathrm{SF}}$ (Springel & Hernquist 2003; Pillepich et al. 2018b). This results in a non-local coupling of the stellar wind feedback to the gas, in contrast to the local feedback from AGN. Feedback from AGN can be either thermal or kinetic, following the model of Springel et al. 2005; Sijacki et al. 2007. Galaxies in the simulation are quenched primarily due to radio mode feedback from AGN, with an expanding jet induced bubble transferring energy from the BH to the halo and heating the gas. Parameters of the Illustris model have been chosen to roughly reproduce the cosmic star formation rate density (SFRD), and the galaxy stellar mass function (SMF), the stellar mass-halo mass relation (SMHM), and the stellar mass-black hole mass relation (SMBH) at $z=0$.

• •

  IllustrisTNG (Pillepich et al. 2018b; Weinberger et al. 2017)

  A significantly updated version of the original Illustris project, IllustrisTNG³³ 3 https://www.tng-project.org/ carries over recipes for star formation and evolution, chemical enrichment, cooling, feedback with outflows, growth and multi-mode feedback from Illustris with substantial updates (Pillepich et al. 2018b; Weinberger et al. 2017; Nelson et al. 2018a). In addition to this, it incorporates new black hole driven kinetic feedback at low accretion rates, magnetohydrodynamics and improvements to the numerical scheme. Unlike Illustris, TNG injects winds isotropically with a modified wind speed that depends on the local 1D dark matter (DM) velocity dispersion, with a redshift dependence that matches the growth of the virial halo mass. AGN feedback is modeled using two modes: a pure thermal mode at high accretion rate (quasar mode) and a pure kinetic mode at low accretion rate (radio mode), with a kinetic wind feedback model (Weinberger et al. 2017) responsible for quenching galaxies (Weinberger et al. 2018). In addition to the observations used with Illustris, the TNG simulation parameters are also chosen to reproduce galaxy sizes and halo gas fractions at $z=0$.

• •

  Mufasa (Davé et al. 2016)

  The Mufasa meshless hydrodynamic simulations use the GIZMO code (Hopkins 2015), prescriptions for cooling and heating with Grackle (Smith et al. 2017), and star formation and feedback from massive stars using scalings from FIRE (Hopkins et al. 2014; Muratov et al. 2015). Star formation is implemented stochastically from gas particles using the Krumholz et al. 2009 formalism to estimate the H₂ formation at coarse resolution accounting for sub-grid clumping. Then, for densities $\geq 0.13\mathrm{cm}^{-3}$, stars are formed stochastically over local dynamical timescales ($t_{\mathrm{dyn}}=1/\sqrt{G\rho}$) with $\sim 2\%$ efficiency, following Kennicutt 1989. Sub-grid recipes for feedback from massive stars launch two-phase winds that drive material out of galaxies through a combination of type-II supernovae, radiation pressure and stellar winds. These winds are parametrized using a mass loading factor and wind speed, and scaling relations for these parameters based on galaxy properties are adopted from the FIRE simulations (Muratov et al. 2015) instead of being tuned to reproduce observations. Since Mufasa does not explicitly model AGN, quenching is accomplished by keeping all the gas in massive haloes heated (except gas that is self-shielded) to reproduce the effects of ‘maintenance mode’ feedback from long lived and AGB stars (Gabor & Davé 2015). Parameters in Mufasa have been chosen to reproduce the galaxy SMF at $z=0$.

• •

  Simba (Davé et al. 2019)

  The Simba cosmological galaxy formation simulations are built on the Mufasa simulations including black hole growth and feedback, using the GIZMO cosmological gravity+hydrodynamics code with its Meshless Finite Mass (MFM) solver (Hopkins 2014; Hopkins 2017). Similar to Mufasa, Simba uses a stochastic H₂ based star formation model, with the SFR given by the H₂ density divided by the local dynamical timescale. Simba also uses two-phase winds with updated mass loading factor scalings from FIRE (Anglés-Alcázar et al. 2017a), which is similar to those adopted by IllustrisTNG but with slightly lower wind velocities. Simba implements a torque limited BH accretion model along with a kinetic subgrid model for BH feedback similar to Anglés-Alcázar et al. 2017b, but with a variable outflow velocity. Wind particles are decoupled for a short amount of time ($10^{-4}\tau_{\mathrm{H}}$, where $\tau_{\mathrm{H}}$ is the Hubble time) from hydrodynamics and radiative cooling. The BH feedback is overall similar to the two-mode model in IllustrisTNG, with some differences detailed in Davé et al. 2019. The majority of galaxy quenching occurs due to the AGN jet mode feedback, with a bimodal distribution of quenching timescales found in Rodríguez Montero et al. 2019. Parameters in the Simba model have been chosen to reproduce the $M_{\mathrm{BH}}-\sigma$ relation and the galaxy SMF at $z=0$.

• •

  EAGLE (Schaye et al. 2015; Crain et al. 2015; Schaller et al. 2015; McAlpine et al. 2016)

  The Evolution and Assembly of GaLaxies and their Environments (EAGLE)⁴⁴ 4 http://icc.dur.ac.uk/Eagle/ is a set of cosmological hydrodynamic simulations of galaxy formation using a modified version of the Tree-PM smoothed particle hydrodynamics (SPH) code GADGET-3 (Springel 2005). EAGLE does not resolve molecular clouds for accurate modeling of the warm gas within galaxies, and implements sub-grid recipes for stellar evolution, cooling and heating of gas due to stars and other emission, metal enrichment of ISM gas and energy injection due supernovae, and the formation, accretion and feedback of AGN. Star formation occurs via gas particles that are stochastically converted into star particles at a pressure-dependent rate that reproduces the observed Kennicutt-Schmidt law (Schaye & Dalla Vecchia 2008). A metallicity-dependent density threshold (Crain et al. 2015) is adopted to ensure that star formation happens in cold, dense gas. The local ISM is heated stochastically due to feedback from massive stars and supernovae with a fixed temperature increment (Dalla Vecchia & Schaye 2012). At high SFR, this feedback can lead to large-scale galactic outflows (Crain et al. 2015). Similar to feedback from star formation, AGN feedback is implemented using a single-mode thermal feedback model. The fraction of radiated energy that couples to the ISM is calibrated to reproduce the stellar mass-black hole mass relation at $z=0$, and mimics the ‘radio’- and ‘quasar’-like modes depending on the BH accretion rate (Crain et al. 2015). Quenching is thought to happen on long timescales ($\sim 3-4$ Gyr) for low mass central galaxies due to stellar feedback, and high-mass centrals on shorter timescales due to AGN feedback and environmental quenching (Trayford et al. 2016; Wright et al. 2019). Parameters in the EAGLE suite are chosen to reproduce the galaxy SMF at $z=0.1$ and the disc galaxy size-mass relation.

• •

  Santa-Cruz SAM (Somerville et al. 2008; Somerville et al. 2015; Porter et al. 2014; Brennan et al. 2016)

  The Santa-Cruz Semi-Analytic Model contains a number of well motivated semi-analytic prescriptions for the hierarchical growth of structure, gas heating and cooling, star formation and stellar evolution, supernova feedback and its effect on the ISM and ICM, AGN feedback, starbursts and morphological transformations due to mergers and disc instabilities that are used in conjunction with the Bolshoi-Planck (Klypin et al. 2011; Rodríguez-Puebla et al. 2016; Klypin et al. 2016) dark matter simulation merger trees to construct populations of galaxies that are tuned to match observations at $z=0$. The model implements two modes of star formation: a ‘normal’ disc mode following the Schmidt-Kennicutt relation, along with exploding supernovae which drive outflows with recycling that occurs in isolated discs, and a ‘starburst’ mode that occurs as a result of a merger or internal disc instability. The SAM implements a multi-phase gas model for the ISM. Cold gas can be ejected from galaxies by winds driven by SN feedback. Heated gas is either trapped within the DM halo potential well, or ejected from the halo into the diffuse IGM. Brennan et al. 2016 and Somerville & Davé 2015 find that virial shock heating due to massive haloes alone is not enough to quench massive galaxies, with a significant role played by feedback from AGN activity, driven by galaxy mergers or in-situ processes like disc instabilities. Model parameters such as the strengths of stellar and AGN feedback are calibrated using the observed SMF at $z=0$, with further details in Porter et al. 2014.

• •

  UniverseMachine (Behroozi et al. 2013; Behroozi et al. 2019)

  The UniverseMachine is an empirical model that determines the SFRs of galaxies as a function of their host haloes’ potential well depths, assembly histories, and redshifts. The model uses halo properties and assembly histories from the Bolshoi-Planck dark matter simulation (Klypin et al. 2011; Rodríguez-Puebla et al. 2017; Klypin et al. 2016) in conjunction with a variety of observational constraints including the cosmic SFRD, observed SMFs, specific SFR functions, quenched fractions, UV luminosity functions, UV-stellar mass relations, IRX-UV relations, autocorrelation and cross-correlation functions, and the dependence of quenching on environment across $0<z<10$ to constrain its free parameters (see Table 1 in Behroozi et al. 2019). Star formation rates are parametrized in terms of redshift and halo properties, with the list of parameters in Table 2 of Behroozi et al. 2019, which include the scatter in the SFRs of star forming galaxies, a model for the SFR-$v_{M,\mathrm{peak}}$ relation⁵⁵ 5 Where $v_{M,\mathrm{peak}}$ is the maximum circular velocity of the halo at peak halo mass., quenched fraction properties and random errors in measuring stellar masses and star formation rates. The parameters are tuned using Markov Chain Monte Carlo optimization to match the observational constraints. In the current analysis we use SFHs from the public Data Release 1 of UniverseMachine.

• •

  FIRE-2 (Hopkins et al. 2018)

  The Feedback In Realistic Environments (FIRE)⁶⁶ 6 http://fire.northwestern.edu simulations considers a fully explicit treatment of the multi-phase ISM, and stellar feedback. The simulations in this work are specificially part of the “FIRE-2” version of the code; all details of the methods are described in Hopkins et al. 2018, Section 2. The simulations use the code GIZMO (Hopkins 2015),⁷⁷ 7 http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html, with hydrodynamics solved using the mesh-free Lagrangian Godunov “MFM” method. Gas dynamics and radiative cooling from a meta-galactic background and local sources are incorporated using tabulated cooling rates from CLOUDY (Ferland et al. 2017). Stars form by stochastically turning gas particles into stellar particles in dense, self-shielding molecular, self-gravitating regions above a density threshold. The stellar feedback prescription includes radiation pressure from massive stars, local photoionization and multi-wavelength photoelectric heating, core-collapse and type Ia supernovae with appropriate momentum and thermal energy injection, and stellar winds. The FIRE physics, source code, and all numerical parameters are identical to those in Hopkins et al. 2018. The higher resolution of the FIRE simulations resolves the ISM to a larger extent than the large-volume simulations. Hopkins et al. 2014 find that supernova feedback alone is not enough, radiative feedback (photo-heating and radiation pressure) is needed to destroy GMCs and enable efficient coupling of later supernovae to gas. Multiple feedback mechanisms are also responsible for regulating the ISM: supernovae regulate stellar masses/winds; stellar mass-loss fuels late star formation; radiative feedback suppresses accretion on to dwarfs and instantaneous star formation in discs. Feedback from supermassive black holes is not included in the simulations (Hopkins et al. 2018). While there are approximations for the momentum and energy deposition from SNe when the cooling radius is not resolved, the simulations are not explicitly tuned.

• •

  g14 (Governato et al. 2012; Munshi et al. 2013; Brooks & Zolotov 2014; Christensen et al. 2016; Brooks & Christensen 2016; Brooks et al. 2017; Christensen et al. 2018)

  The g14 suite of cosmological zoom simulations are run using the N-body+SPH code Gasoline (Wadsley et al. 2004) within a WMAP3 cosmology. The galaxies are chosen to have a range of merger histories and spin values. The g14 simulations follow the non-equilibrium formation and destruction of molecular hydrogen, and allow stars to form in the presence of H₂, with resolution high enough to resolve the disks of galaxies and the GMCs in which stars form (Christensen et al. 2012). Stars are born with a Kroupa et al. 1993 IMF, mass and metals are returned in stellar winds as star particles evolve and SN Ia and II return thermal energy to the surrounding gas (see Stinson et al. 2006 for details). For SN II, 10⁵¹ erg of energy are injected per SN. Metal diffusion occurs in the ISM (Shen et al. 2010), and a cosmic UV background is included following Haardt & Madau 1996. The g14 suite was calibrated to match the SMHM relation of Moster et al. 2013.

• •

  Marvel/Justice League (Bellovary et al. 2019, Munshi et al., in prep.)

  The Marvel/Justice League simulations are run using ChaNGa (Menon et al. 2015), the successor to Gasoline. The Marvel-ous dwarfs (henceforth Marvel) are a sample of field dwarfs (4-11 Mpc from a Milky Way-mass galaxy) at 65pc force resolution, while the DC Justice League (henceforth JL) are zooms of MW-mass disk galaxies and their surrounding environments at 170pc resolution. Many of the physics modules in ChaNGa remain the same as in Gasoline, such as the star formation and stellar feedback schemes, with the exception that 1.5$\times$10⁵¹ erg of thermal energy is injected per SN II. This increase is motivated by the fact that ChaNGa contains an improved implementation of Kelvin-Helmholtz instabilities compared to Gasoline (Wadsley et al. 2017), which leads to more efficient accretion onto the disk. An updated UV background is adopted, based on Haardt & Madau 2012. In addition, supermassive black hole growth and feedback is implemented using the models described in Tremmel et al. 2017. Parameters in the simulations were calibrated to reproduce the SMHM, SMBH, and SFRs of galaxies at $z=0$.

We compute SFHs for each galaxy in the hydrodynamical simulations under consideration (Illustris, IllustrisTNG, Mufasa, Simba, EAGLE, FIRE-2, g14 and Marvel/Justice League ) by performing a mass-weighted binning of the star particles with $\Delta t=100$ Myr. The choice of time bin is further explored in Section 2.3.1. For models where we only have access to the masses of the star particles at the time of observation, we account for mass-loss using the FSPS (Conroy et al. 2009; Conroy & Gunn 2010) stellar population synthesis code, adopting the initial mass function (IMF) of the stellar particles in the simulation. In this procedure, we consider all stellar particles belonging to a galaxy at $z\sim 0$ instead of tracing the gas-phase SFR as a function of time. The reason for this is twofold: (i) since the hydrodynamical simulations trace the times when star particles were formed, this gives us finer time-resolution than the snapshots saved for the different simulations, (ii) since the SFHs we observationally reconstruct are the sum over all the progenitors, this archaeological approach therefore allows us to compare directly with observations. Both UniverseMachine and the Santa Cruz SAM track the SFR for each galaxy, so we simply interpolate these to match the same time grid with $100$ Myr steps as the hydrodynamical simulations. In both cases, the resolution is fine enough that the interpolation does not need to up-sample the SFR. Since the UniverseMachine SFHs are stored in terms of scale factor instead of absolute time, an additional periodogram is computed using the uneven spacing to check that the PSDs are not significantly affected by the interpolation. Additional fine-resolution SFHs are also computed for the galaxies from the zoom simulations, with a timestep $\Delta t=1$ Myr.

Since the SFHs span a large dynamic range, we work in log SFR space in order to be able to better quantify the relative strengths of SFR fluctuations. Analyzing the SFHs in linear SFR space effectively amounts to a different weighting scheme. This choice is motivated by physical considerations, since the SFRs of star forming galaxies are often found to be distributed normally in log SFR space, with a tail towards low SFRs from passive galaxies that do not have ongoing star formation (Feldmann 2017; Hahn et al. 2019b; Caplar & Tacchella 2019).

The variability or ‘burstiness’ of galaxy SFRs is a topic of much interest, and has been studied in a variety of ways - using burstiness indicators based on the timescales of different SFR tracers (Guo et al. 2016; Emami et al. 2019; Broussard et al. 2019), fitting an exponential to the Pearson correlation coefficient of SFRs as a function of time-separation to quantify an ‘SFR evolution timescale’ (Torrey et al. 2018), quantifying the scatter in SFRs smoothing on different timescales (Hopkins et al. 2014; Matthee & Schaye 2019), using power spectral densities (PSDs) to quantify the variability in Fourier space (Caplar & Tacchella 2019; Wang & Lilly 2020b; Wang & Lilly 2020a), or performing a PCA decomposition of SFHs to get estimate the fraction of variance accounted for by different timescales (Matthee & Schaye 2019). In other studies involving timeseries data, the structure function (Hughes et al. 1992; MacLeod et al. 2010; Kozłowski 2016; Caplar et al. 2017) has also been used as a metric to quantify variability on different timescales in quasar and AGN studies.

In the current analysis, we choose to quantify the variability of galaxy SFHs using the PSD formalism, since

• •

  the PSD formalism is well studied and easily interpretable, and Fourier space provides an excellent domain to quantify and compare the variability of SFHs across different timescales;

• •

  the decomposition of variability into different frequencies, and therefore different timescales, allows us to understand the relative contribution to the overall burstiness from each timescale. This takes us one step closer toward relating this variability to the underlying physical processes responsible; and

• •

  evolving analysis techniques coupled with upcoming observational datasets will make it possible to obtain observational constraints in PSD space.

For a continuous time series $\psi(t)$, the PSD is defined in terms of the Fourier transform $f(k)=\int dt~e^{-ikt}\psi(t)$ as $\mathrm{PSD}(k)=|f(k)|^{2}$. In practice, we compute the PSD for each SFH using Welch’s method (Welch 1967), implemented in the scipy.signal.welch module.

The PSD corresponding to the SFH for an individual galaxy reports a phase-averaged estimate of the strength of SFR fluctuations at a given frequency⁸⁸ 8 Or, inverting it, at a given timescale..

For a sinusoidal signal with a frequency $\nu$, the corresponding PSD is given by a delta function at the frequency $\nu$, shown in the top panel of Figure 3. Generalized to more complicated timeseries, the PSD therefore provides a way to disentangle and interpret the strength of the fluctuations on different timescales, as previously done in studies of AGN variability and theoretically with SFHs (MacLeod et al. 2010; MacLeod et al. 2012; Caplar et al. 2017; Sartori et al. 2018; Caplar & Tacchella 2019). A sharp peak in the PSD would indicate strong SFR fluctuations at a given timescale, possibly driven by a physical process. However, physical processes acting over a range of timescales spread out the peaks and make it more difficult to isolate the effects of individual processes. An example of this is shown in the middle column of Figure 3, where the sum of three sinusoidal curves produces three peaks in the PSD, along with additional artifacts due to the finite length of the time series. Processes like hierarchical structure formation and correlated stochastic star formation additionally link short timescales to longer ones, creating an overall spectral slope to the PSD, shown in the bottom panel of Figure 3. Physical processes can additionally drive features at certain characteristic timescales, for e.g., the regulator model (Lilly et al. 2013, see also Bouché et al. 2010; Davé et al. 2012; Forbes et al. 2014b) predicts SFHs correlated below an equilibrium timescale of a galaxy’s gas reservoir, with the slope at timescales below the break steeper by 2 than the slope above it (Wang & Lilly 2020a; Tacchella et al. 2020). Such features can be seen as breaks in the PSD (Caplar & Tacchella 2019). The PSDs of galaxy SFHs therefore contain a wealth of information about the different physical processes responsible for its shape.

Examples of this procedure are shown in Figure 4, which shows SFHs for galaxies from the IllustrisTNG simulation (left column) as well as their corresponding PSDs (right column). The contribution to the PSDs at three different timescales due to the strength of SFH fluctuations are highlighted in different colours in the middle panels. Unlike the case for the sine wave, the power in these PSDs is spread over a large dynamic range, indicative of the stochastic nature of star formation and the wide range of timescales over which physical processes in galaxies induce variability in the star formation rates. With a thorough understanding of a galaxy’s evolution and merger history, it might be possible to interpret its individual PSD. However, in the current work we focus on studying the broader trends in a sample of galaxy SFHs and their corresponding PSDs as a way to compare different models of galaxy evolution on the same footing. In doing so, we examine the variability of galaxy SFHs on intermediate ($\sim 200$ Myr) to long ($\sim 10$ Gyr) timescales, and study the evolution in the PSDs as a function of stellar mass and star forming state (star forming vs quiescent). We choose stellar mass since it is a good tracer for the overall state of a galaxy, correlating well with a wide range of other physical properties including halo mass, SFR, metallicity and BH mass, and can be calculated self-consistently for all the models directly from the SFHs after accounting for mass-loss.

In this section, we study how the different models deviate from the overall sample behaviour (and from each other) by quantifying the overall SFH diversity as a function of time and stellar mass. We compute the median SFH in a given mass bin for individual models and compare it to the median SFH in a given mass bin across all the models. To account for the differing number of galaxies in a given mass bin across the different models, we randomly sample 1000 SFHs with replacement from the available SFHs at each step of the calculation. We repeat this sampling and calculation 100 times to adequately capture small ($\sim 0.02-0.1$ dex) fluctuations due to random seeds.

The result is shown in Figure 10, which shows the difference between the median SFHs in different bins of stellar mass. It should be noted that the median SFH across all models is not the ‘correct’ SFH, but merely a guide to the eye. Therefore, instead of comparing the deviation from the median for any given model, it is more instructive to (i) look at the differences between the models themselves, and (ii) use the median to get an idea of the overall variance among models at a given mass and epoch (shown as shaded grey regions). Although there is a considerable diversity across the different models, the largest differences occur when the SFR is low – at early epochs when galaxies are beginning to assemble their mass and when they are quenching. A locus of agreement across the various models exists in each mass bin, moving to higher redshifts with increasing mass. This is correlated with the epoch when the median SFHs peak in their SFR, as seen in Section 3.1.2. This means that despite these differences, the overall picture of galaxy mass assembly described by the models is similar.

At late times (low redshifts), there is an increase in the overall variance between the different models with increasing mass, ranging from $\sim 0.3$ dex at M${}_{*}\sim 10^{10.0}$M_⊙ to $\gtrsim 1$ dex for massive galaxies (M${}_{*}>10^{10.5}$M_⊙). The median SFHs across all models agree well at $z<1$ in the lowest mass bin. These trends are not easy to interpret since, as we discuss in Section 2.2, these SFHs are tracing the SFR of the main progenitor as well as of all the accreted systems. This means that this late time divergence is probably a combination of how the various models implement quenching as well as the SFH of the accreted systems. Although a full treatment studying the cause of these differences is outside the scope of this analysis, quantifying the differences between the PSDs for these SFHs begins to illustrate how differing strengths of SFR fluctuations across a range of timescales could shape the overall SFHs over the next few sections.

Additional plots showing the distributions of SFH parameters such as stellar mass, sSFR, SFH peak and width for the various models can be found in Appendix B.

Section 3.1 describes some of the overall trends in the PSDs - the distribution of power across a broad range of timescales, with an increase in power towards longer timescales / shorter frequencies. However, each model shows unique trends in how the PSDs evolve with stellar mass, as well as the actual strength of the PSD at different timescales.

Since there are a range of modeling assumptions and numerical recipes used across the various models we consider, a comparison in PSD space serves to highlight the differences in the resulting variability of their SFHs on different timescales. Figure 11 shows where the median PSDs of galaxies from the various models (as shown in Figures 5, 6, and 7) lie in PSD slope vs PSD power space at three representative timescales ($300$ Myr, $1$ Gyr and $10$ Gyr), and an interactive version of this plot spanning timescales ranging from $\sim 200$ Myr to $13$ Gyr can be found online⁹⁹ 9 https://kartheikiyer.github.io/psd_explorer.html. An equivalent plot showing individual galaxies from the zoom simulations is shown in Figure 12. The PSD power is the strength of SFR fluctuations or ‘burstiness’ at a given timescale. The local slope of the PSD at a given timescale is computed using the PSD within a log timescale of $\pm 0.1$ dex, and is a measure of how tightly coupled the PSD is to adjacent timescales. Changing this interval while computing the slope does not affect the overall trends across the models. A slope of 2 can be found in models of stochastic star formation described by random walks (Caplar & Tacchella 2019; Kelson et al. 2020). Tacchella et al. 2020 find this to emerge naturally within the framework of the gas regulator model (Lilly et al. 2013) and in modeling stochasticity due to GMC formation and destruction. Most high-mass and low-sSFR galaxies across the different models show a slope $\sim 2$, while UniverseMachine shows this at all stellar masses. Individual points for each model show the median slope and power of the PSDs in the same $0.5$ dex bins of stellar mass that are used in Figures 5, 6, and 7, highlighting evolution in PSD space as galaxies grow more massive.

A key point to note is that the various models are extremely diverse in (i) the region of PSD space they occupy at a given timescale, and (ii) their evolution with stellar mass. At $\sim 300$ Myr timescales, there seems to be an overall attractor toward increasing slope and decreasing power as the stellar mass increases, although this does not hold for all the models. The short timescale power generally increases as a function of decreasing stellar mass, indicating that lower mass galaxies are generally more bursty across a variety of models. This trend is not as prominent for the SAM and empirical model. Meanwhile, at $\sim 1$ Gyr timescales the models seem to follow a range of different behaviours, although most models seem to converge on a PSD slope of $\beta\sim 2$ at high stellar masses. On the longest timescales, both slope and power tend to increase with increasing stellar mass, in part due to the increased contribution from quenched galaxies to the long-timescale PSD power. Depending on how the individual models implement quenching, however, the rate and extent of this effect can vary greatly (see Section 3.4).

It should be noted that although these trends are shown using the median values for the PSD slope and power in a given mass bin, there is a large amount of variance in the range of slope and power values possible for individual galaxy PSDs due to features that may be present in individual galaxy SFHs based on stochastic events like halo accretion fueled star formation and major mergers. The variance in slopes is from $\sigma($PSD slope$)\sim 0.7-1.0$ (dex)² and in power is $\sigma($log PSD$)\sim 0.2-0.7$(dex)²Myr, corresponding to the shaded regions in the individual PSD plots and generally increasing with increasing stellar mass. While the large variance indicates that individual galaxies in a given mass range exhibit a large diversity in behaviour, trends across stellar mass are generally robust since they trace the behaviour of the entire population.

Given that the models span such a wide range in PSD slope and power at any given mass and timescale, observational constraints in this space (Caplar & Tacchella 2019; Wang & Lilly 2020b) would provide strong constraints on modelling galaxy physics.

Quenching becomes an increasingly important phenomenon as we consider galaxies with higher stellar masses. This phenomenon can be driven by a range of different physical processes acting on different timescales. Since the quenching of galaxies is an observably measurable phenomenon, it is therefore possible to get observational constraints on quenching timescales and connect them to the relevant physical processes. Here we explore the differences in the PSDs of actively star-forming and quiescent galaxies at $z=0$ to determine what, if any, differences they show at different stellar masses.

To perform this analysis, we first need to select galaxies that are quiescent at the time of observation. There exist multiple ways of performing this selection, depending on the definition of quenching (e.g., through a cut in UVJ space, in specific SFR, or a threshold distance from the SFR-M_∗ correlation, among others, see for example Donnari et al. 2019; Hahn et al. 2019b). In the current analysis, we separate galaxies into star forming vs quiescent using a commonly used threshold in sSFR ($\mathrm{sSFR}<0.2/\tau_{\mathrm{H}}\sim 10^{-10.83}yr^{-1}$ for quiescent galaxies at $z\sim 0$, see Pacifici et al. 2016; Carnall et al. 2019). This approach is motivated by two reasons: (i) since we already have access to the SFHs, this allows us to avoid the systematic assumptions of forward modeling rest-frame UVJ colours and the degeneracies of separating quiescent galaxies in that space, and (ii) we avoid the systematics of accounting for different SFR-M_∗ correlations across the different models (Hahn et al. 2019b) and use a uniform threshold for comparison across the models.

Having identified quiescent galaxies across the various models, we then compare the PSDs of quiescent galaxies to those of star-forming galaxies at different stellar masses. Since quenching distinctively alters the shape of a galaxy’s SFH, we expect the PSDs of quiescent galaxies to show more power on long timescales. Figure 13 shows the difference in the median log PSDs of quiescent and star forming galaxies in the same $0.5$ dex stellar mass bins used in the rest of this work. We exclude the highest mass bin (M${}_{*}\sim 10^{11.5}$M_⊙), since there are not enough star forming galaxies in all the models to perform this analysis.

We see that at low stellar masses, the quiescent galaxy PSDs generally show greater power on long timescales, with the exact timescale varying across models, ranging from $\sim 900$ Myr to $\gtrsim 3$ Gyr. As we go to higher stellar masses, we find that there is an excess of power across a range of shorter timescales in the IllustrisTNG, Simba, EAGLE, and SC-SAM models. This could be explained by processes driving quenching also driving variability in SFRs on other timescales. For example, multiple short episodes of feedback due to (i) AGN-driven outflows leading to the eventual quenching of galaxies, as seen in the implementation of jet mode AGN feedback (Rodríguez Montero et al. 2019) or (ii) X-ray feedback rapidly evacuating the star forming gas in the central regions (Appleby et al. 2020) could lead to increased short-timescale variability in Simba. In contrast, since Mufasa implements quenching primarily through a ’maintenance mode’ feedback, it does not show a strong evolution with stellar mass. Another explanation for this increase in power on short timescale for quiescent galaxies could be that the quiescent galaxies assemble their mass earlier, when SFHs in general were more bursty (Muratov et al. 2015; Hayward & Hopkins 2017). The phenomenon of quiescent galaxies assembling their mass earlier can be seen the correlation between t₅₀ and sSFR for quenched galaxies among the various models shown in Appendix B. However, the nature of this correlation is uniform across all the models and can not fully account for the variations in the difference between star-forming and quenched galaxy PSDs observed between the models.

In more detail, the Illustris and IllustrisTNG models both show sharp breaks above which the power in quiescent galaxies rises, with the break occurring on longer timescales in IllustrisTNG. The increase in power on short timescales with increasing mass is also more prominent in TNG compared to Illustris. The updated winds and AGN feedback in IllustrisTNG also show a noticeable increase in power on short timescales above masses $10^{10.5}$M_⊙, where AGN feedback becomes most effective. EAGLE shows a much broader range of timescales in comparison, similar to Simba albeit with high power at a given mass. The SAM shows a significant increase in power on timescales below $\sim 2$ Gyr, with this trend increasing with stellar mass. This seems to be primarily associated with stochastic starbursts on short timescales triggered by mergers, with more massive galaxies experiencing these events to a larger extent. UniverseMachine shows a moderate increase in power with quenching over all timescales.

In summary, PSDs across the different models show a range of behaviours when galaxies quench, with strong mass dependence in some models (IllustrisTNG, Simba, EAGLE, SC-SAM) and a range of timescale-specific breaks in the PSD ($\sim 900$ Myr in Illustris, $\sim 2-3$ Gyr in IllustrisTNG, $\sim 3-4$ Gyr in Simba, and $\sim 2$ Gyr in the SC-SAM). Observational constraints in PSD space for star forming and quiescent galaxy populations will provide sensitive probes of discriminating between the range of quenching mechanisms implemented across these models.

Upon examining the PSDs of galaxy SFHs across different models, we find that most of the power resides in the long timescales on which SFHs rise and fall. At early cosmic times, several models find the SFHs of galaxies to be correlated with the dark matter accretion histories (DMAHs) of their parent haloes (Wechsler & Tinker 2018). Diemer et al. 2017 model galaxy SFHs as log-normal curves, and find that the peak and width of SFHs in Illustris correlate strongly with the properties of their DMAHs with an offset between the formation times of haloes and galaxies that increases with stellar mass, along with a tight relation between the BH mass and peak time. Similarly, Qu et al. 2017 show that the SFHs of galaxies in EAGLE increasingly decorrelate from the halo accretion histories at increasing masses, by plotting the formation time vs accretion time for haloes and galaxies across stellar mass bins. They find this to be due to AGN feedback, which suppresses in-situ star formation and causes the stars in massive galaxies to form early and the galaxies to grow subsequently by mergers (i.e., the majority of star formation finished early), while haloes continue accreting mass until late times (i.e., massive haloes form late) (Neistein et al. 2006).

From an analytical standpoint, Kelson 2014 models star formation as a stochastic timeseries, with the ‘long-timescale memory’ encapsulated by a Hurst parameter of $\sim 0.98\pm 0.06$. In Kelson et al. 2016, this model is extended to derive stellar mass functions at early times, explicitly relating the variance of the SFRs for an ensemble of galaxies to the DM haloes and their ambient matter densities at the epoch when star formation begins. Kelson et al. 2020 analytically estimate the slope of the DMAH PSD to be $\sim 1$.

With this in mind, it would therefore be instructive to (i) compute the PSDs of DMAHs and study their behaviour, (ii) study the extent to which variability in galaxy SFHs is tied to the variability in the DMAHs of their parent haloes, and (iii) examine if SFHs and DMAHs are coherent, to understand if dark matter accretion drives star formation.

There exist a range of estimates in the literature for timescales associated with different physical processes, some of which are shown in Figure 1 and listed in Appendix C. In this Section, we briefly summarize the current state of our understanding regarding which physical processes can contribute to SFR fluctuations at various timescales. By doing this, we can begin to connect the different features seen in the median PSDs of SFHs in Section 3.1.2 to the underlying physical processes responsible. It also serves as a useful starting point for future analyses looking at these features in greater depth within specific models. Starting with processes that act on the shortest timescales, we gradually work our way to the longer timescales that are the focus of the bulk of this paper.

GMC formation and destruction: Star formation on small spatiotemporal scales occurs in GMCs, whose lifetimes are sensitive to a variety of factors including cloud collisions and mergers, feedback from supernovae, cosmic rays and photoionisation, turbulence in the ISM, and the growth of magnetic fields (Dobbs et al. 2012; Dobbs et al. 2015; Kim & Ostriker 2017; Semenov et al. 2017; Pakmor et al. 2017; Benincasa et al. 2019). Current upper bounds on theoretical predictions for GMC lifetimes range between $\sim 7-20$ Myr (Tasker 2011; Benincasa et al. 2019), with estimates for the timescales of individual processes that influence GMC lifetimes reported in Appendix C. Analytical models can also provide an understanding of when star formation in this regime can be bursty (Faucher-Giguère 2018).

Considering the rate of GMC formation and destruction to be a stochastic process, we would therefore expect a power-law PSD with slope $\beta\sim 2$ at these timescales (Kelson 2014; Tacchella et al. 2020). Although the large-volume models do not probe these timescales, the three suites of zoom simulations allow us to test this hypothesis. In fact, we do find the PSD in this timescale to be well-described by power-laws, and the g14 and Marvel/Justice League galaxies show slopes of $\beta\sim 1.6^{+0.4}_{-0.1}$ uncorrelated with stellar mass, while the FIRE-2 galaxies show slopes of $\beta\sim 1.8^{+0.5}_{-0.4}$, with a mild trend of increasing slope with stellar mass over timescales of $\sim 10-20$ Myr.

Dynamical processes within galaxies: A range of physical processes act to influence the state of the ISM on galaxy dynamical timescales ($\sim 10^{8}$ yr). These processes include turbulence in the ISM, molecular gas encountering spiral arms and bars, galactic winds, and the rapid cycling of ISM gas between star forming and non-star forming regions, in addition to the exponential growth of magnetic fields, and stochastic inflows of CGM gas¹¹¹¹ 11 The last two extend to longer timescales as well.. Analytical models account for these processes through a range of timescales, including timescales for gas accretion and cooling, as well as star formation, turbulent crossing and effective viscous timescales that describe how long it takes for accreted gas to reach the center of the galaxy (Dekel et al. 2009; Krumholz & Burkert 2010; Forbes et al. 2014a). For modeling these processes, resolution plays an extremely important role since resolving the ISM allows simulations to capture the effects of turbulence driven by feedback, as well as model the feedback self-consistently while relaxing the need for sub-grid recipes. Most of our knowledge in this regime comes from small-volume simulations (e.g., a slice of a galactic disk, Kim & Ostriker 2017, or an idealized disk Semenov et al. 2017) or zoom simulations focusing on individual galaxies (Hopkins et al. 2014; Ceverino et al. 2014; Christensen et al. 2016). Bursty star formation has been noted on timescales of $\sim 45$ Myr in the TIGRESS framework (Kim & Ostriker 2017), and on $\leq 100$ Myr timescales in FIRE-2 (Sparre et al. 2017; Hung et al. 2019).

Since there are many competing factors at play, we expect the PSDs in this regime (and beyond) to be complicated, and this is what we generally see in all the zoom simulation suites. Overall, while the PSDs can still be approximated with a power-law, several PSDs show minor peaks¹²¹² 12 $\sim 35-60$ Myr for Sandra in Justice League, Rogue in Marvel and in m12m, m11e, m11d, m11i, m11v, and m12i in FIRE-2 or breaks¹³¹³ 13 $\sim 60-90$ Myr for h986 from g14 and for m12b, m11f, m11i, m11g, m11c in FIRE-2 with average slopes in the $\sim 30-100$ Myr range of $\beta\sim 2.0^{+0.8}_{-0.7}$ for g14 and Marvel/Justice League and $\beta\sim 1.3^{+0.4}_{-0.5}$ for FIRE-2. All three suites of simulations show increased scatter in the power-law slopes, along with a trend of increasing slope with stellar mass in this timescale range, perhaps correlated with decreasing dynamical timescales as galaxies grow more massive.

Mergers: Mergers between galaxies bring in a combination of stars that have already formed and gas that can fuel a burst of subsequent star formation, with timescales ranging from $\sim 100-500$ Myr (Hernquist 1989; Barnes & Hernquist 1991; Barnes & Hernquist 1996; Mihos & Hernquist 1996; Robertson et al. 2006b; Hani et al. 2020). The effect on SFHs comes from mergers as a a primary mechanism for driving starbursts in galaxies (in addition to disk instabilities) and as a controversial trigger for quenching, depending on a variety of factors including the mass ratio, relative alignment, how gas-rich the merger is, and even if the merger triggers a central AGN (Hopkins et al. 2006; Governato et al. 2009). Zoom simulations also predict that mergers or counter-rotating streams can lower the angular momentum of the gas disk within galaxies, leading to a compaction of the gas phase, which results in an enhancement of the SFR (Zolotov et al. 2015). These phases can last for one to a few hundred Myr and move galaxies to the upper envelope of the star-forming sequence (Tacchella et al. 2016). Rodríguez Montero et al. 2019 find that major mergers cause enhanced SFR at all masses below a threshold of $\sim 10^{11}$M_⊙ in Simba. Tacchella et al. 2019 find trends consistent with centrally enhanced star formation due to ex-situ star formation for intermediate mass ($10^{10-11}$M_⊙) galaxies, with mergers responsible for over two thirds of the ex-situ component toward the high-mass portion of that range. A notable consequence of this is the increasing loss of coherence between the in-situ and full SFHs of galaxies in Figure 16 with increasing stellar mass.

In addition to the timescale of SFR enhancement following a merger, we also need to consider the fact that mergers themselves are stochastic events, and therefore carry an additional implicit timescale. Estimates of merger timescales are generally $\mathcal{O}(1)$ Gyr (Boylan-Kolchin et al. 2008; Lotz et al. 2011; Snyder et al. 2017), and can vary significantly depending on assumed definitions and factors like pair separation and angular momentum of the system. Due to these factors, it can be difficult to isolate the effects of mergers on galaxy PSDs.

Baryon Cycling: The global efficiency of how galaxies are able to convert their gas into stars is almost an order of magnitude different from local efficiencies in star forming regions. Semenov et al. 2017 tie this to the cycling of ISM gas between regions that are star forming and those that are not. In addition to this, gas that leaves the galaxy due to ejective feedback and returns also contributes to prolonging the period over which a galaxy continues to form stars (Christensen et al. 2016; Hopkins et al. 2018, also see review by Tumlinson et al. 2017 and references therein). The lifetimes and dynamics of cold clouds in the halo are also subject to a variety of timescales (Forbes & Lin 2019).

Estimated timescales for the cycling of baryons span a wide range, from $\sim 100$ Myr to about 3 Gyr (Oppenheimer et al. 2010; Christensen et al. 2016; Mitra et al. 2016; Anglés-Alcázar et al. 2017a; Grand et al. 2019). Some studies find the timescales to scale with halo or stellar mass (Oppenheimer et al. 2010; Mitra et al. 2016), while other studies find it to be largely independent of mass (Christensen et al. 2016). While we find evidence for peaks and breaks in the PSDs of individual galaxies on these timescales, especially in the zoom simulations (for example, in Figure 9), the broad range of timescales and the dependence on galaxy properties other than stellar mass results in these peaks being washed out in the median behaviour for an ensemble of galaxies. However, it is possible that breaks in the PSD could be correlated with baryon cycling processes, and bears further investigation in future work. In particular, the evolution of the break timescales with stellar mass in different models could help us understand why some studies show a significant mass-dependent trend while others do not. However, since mergers and other factors also play a role at these timescales, their effects also need to be accounted for in such an analysis.

The $\sim 1$ dex excess in the power of SFHs compared to DMAHs after accounting for quenching could correspond to contributions from baryonic processes like mergers and baryon cycling occurring on halo dynamical timescales as the galaxy grows, leading to imprints in the PSD on timescales $\propto 2\pi\tau_{\mathrm{dyn}}\sim 2\pi(0.1\tau_{\mathrm{H}})\sim 2.1-8.6$ Gyr over the past $\sim 10$ Gyr¹⁴¹⁴ 14 Although it is outside the scope of the current work, it would be an interesting exercise to model the excess in the SFH PSDs as an aggregate effect of baryonic processes across a range of redshifts using a broken power-law model with $\tau_{\mathrm{break}}\sim 2\pi\tau_{\mathrm{dyn}}(z)$, based on the formalism described in Caplar & Tacchella 2019 and Tacchella et al. 2020.. Since the dynamical time grows with decreasing redshift, the resulting contribution to the PSD would end up being smoothed out over a broad range of timescales. The plateau in simulations like Illustris and IllustrisTNG and individual galaxies in the zoom simulations at $\sim 1-3$ Gyr are also indicative of a decorrelation timescale that naturally arises in damped random walk models of star formation (Caplar & Tacchella 2019; Tacchella et al. 2020).

Quenching: Quenching in central galaxies can happen due to a lot of different factors - the shock heating of virial halo gas preventing cold-mode accretion (Dekel & Birnboim 2006), energy from AGN jets that heat gas and prevent it from forming stars (Somerville et al. 2008) and outflows that could remove cold gas from the galaxy (Di Matteo et al. 2005). Observational scaling relations like the $M_{\mathrm{BH}}-\sigma$ correlation tie behavior on large (galaxy-wide) scales to sub-kpc scales on which SMBHs grow, leading to a unique scenario where sub-grid recipes for implementing BH growth and feedback affect when and how galaxies quench. In addition to this, recipes for how simulations implement cooling and star formation, and the strength of winds that blow gas out of galaxies all contribute to the overall trends seen in galaxy quiescence. Finally, the haloes of galaxies set the inflow rate of gas into the central galaxy, as seen through the correlation on long timescales between the DMAHs and in-situ SFHs. This dependence could tie the fueling of the central AGN to that parent halo, ultimately determining when the onset of quenching occurs (Chen et al. 2019).

Figure 6 in Wright et al. 2019 shows a broad, unimodal distribution of quenching timescales in EAGLE galaxies extending out to $\tau_{\mathrm{H}}$ with a median of $\sim 2.5-3.3$ Gyr for low-mass centrals and at shorter timescales (median of $\sim 1.7-2.1$ Gyr) for high-mass centrals depending on the definition of quenching timescale. Longer quenching timescales at low masses are associated with stellar feedback prolonging star formation activity, while shorter timescales at high masses are associated with AGN activity. Simba, on the other hand, shows a bimodal distribution of quenching timescales (Rodríguez Montero et al. 2019), with a slow mode acting approximately over a dynamical time ($t_{Q}\sim 0.1\tau_{\mathrm{H}}$) that is more numerous overall for central galaxies, and a fast mode ($t_{Q}\sim 0.01\tau_{\mathrm{H}}$) that dominates at stellar masses of $M_{*}\sim 10^{10}-10^{10.5}$M_⊙. The fast quenching mode is associated with AGN jet quenching causing a rapid cessation of accretion, since it becomes active at this mass range, and merger rates are not preferentially elevated at these masses. Additionally, X-ray feedback can rapidly evacuate the central regions of galaxies (Appleby et al. 2020) and contribute to short-timescale variability. Sales et al. 2015 find a quenching timescale of $\sim 2-5$ Gyr for galaxies in Illustris. Nelson et al. 2018a find the colour-transition timescale, a tracer of the quenching timescale, to be $\sim 0.7-3.8$ Gyr for IllustrisTNG galaxies. Additionally, Joshi et al. 2020 find that morphological transformations in IllustrisTNG clusters occur on timescales of $\sim 0.5-4$ Gyr after accretion, with a control group showing a broader distribution. They also find that morphological transformation lags $\sim 1.5$ Gyr behind quenching for gas-poor disks, while it precedes quenching by $\sim 0.5$ Gyr for gas rich cluster galaxies, and by $\sim 2.5$ Gyr for gas-rich control galaxies.

In studying the excess PSD power on different timescales and stellar masses due to quenching, we find that the excess variability in IllustrisTNG on short timescales rises strongly at M${}_{*}\geq 10^{10.5}$M_⊙, correlated with the onset of strong kinetic-mode AGN feedback at M${}_{\rm BH}\sim 10^{8.5}$M_⊙ (Weinberger et al. 2018). While we are not in a position to speculate about timescales of $0.01\tau_{\rm H}$, we do find a tail of excess variability extending to the lowest timescales in Simba that could be related to the jet-mode AGN feedback. In EAGLE, Wright et al. 2019 find that galaxies at low masses quench primarily due to stellar feedback on long timescales, consistent with the excess power we see on timescales $\geq 2$ Gyr. As galaxies grow more massive (M${}_{*}\geq 10^{10.3}$M_⊙) mergers and black hole activity increase sharply, leading to overall shorter quenching times, and additional variability on all timescales, as seen in Figure 13. The excess power in the SAMs at high masses seems to be primarily due to increased stochastic starbursts triggered by mergers and AGN activity (Somerville et al. 2008). The increasingly featureless (scale-free) nature of the PSDs toward high stellar masses, where the fraction of quenched galaxies is the largest, could be due to the contribution to the PSD from quenching dominating all other contributions.

Since a combination of multiple processes is responsible for quenching at different stellar masses, it is difficult to constrain their relative strengths with observational measurements of quenching timescales. However, since these different processes also induce varying amounts of short-timescale variability, constraints in PSD space might be able to distinguish between processes and allow for better constraints on their relative strengths.

For the different galaxy evolution models we consider in Section 3, we see a large diversity in the amount of power in SFR fluctuations on a given timescale, the coupling between adjacent timescales, and the existence and location of breaks in the PSD. This makes the PSD a sensitive probe of both the strengths of physical processes and their numerical implementation in these models. Observational constraints in this space are therefore extremely important, and will allow us to better constrain the relative strengths of different processes for a population of galaxies at a given stellar mass and epoch. These observational constraints can come in three forms: (i) constraints in PSD space obtained by measuring the SFR variability of ensembles of galaxies, which can be compared to the models we study, (ii) constraints on the timescales for observed phenomena like quenching or rejuvenation, which can be tied to breaks or peaks in the PSD, and (iii) constraints on timescales of physical processes, which can be used to isolate the effects of different processes contributing to the PSD on a given timescale. In this section, we will briefly discuss each of these.

Ensemble constraints on SFR variability: The spectral energy distributions (SEDs) of galaxies are composed of spectrally distinct contributions from stellar populations formed at different ages relative to the time of observation. Interpreting these contributions gives us access to star formation rates averaged over different timescales. Nebular emission from the regions near short-lived O- and B-type stars provides constraints on SFR over the most recent $\sim 4-10$ Myr (Madau & Dickinson 2014). The rest-UV portion of the SED contains contributions from young stars that probe the SFR out to $\sim 30-100$ Myr, with a similar timescale probed by the rest-FIR portion of the SED, which contains the re-emitted light from the young stellar light absorbed by dust. In addition to this, features like the strength of H$\delta$ absorption and the 4000Å  break are sensitive to SFR within the last $\sim 1$ Gyr and to the light-weighted age within $\sim 2$ Gyr, respectively (Kauffmann et al. 2006; Wang & Lilly 2020b).

To get constraints in PSD space, it is useful to consider the SFR distributions for populations of galaxies and compare these distributions on different timescales to quantify a relative change in burstiness¹⁵¹⁵ 15 This makes an inherent assumption of ergodicity, that the PSDs obtained from a population of galaxies can be connected to the PSDs of individual galaxy SFHs over time. This assumption is explored in detail in Wang & Lilly 2020a.. While this straightforward to forward-model for the large-volume models, the small number of zoom galaxies make this a more involved procedure while considering those models. In these cases, a workaround is possible by realising samples from the PSDs of zoom galaxies, similar to the procedure followed in (Tacchella et al. 2020). In terms of the PSD formalism, this is equivalent to an observational constraint on the slope of the PSD between two timescales. This has been done for select timescales and populations of galaxies in Guo et al. 2016; Broussard et al. 2019; Emami et al. 2019. More recently, Caplar & Tacchella 2019 and Wang & Lilly 2020b; Wang & Lilly 2020a have performed analyses motivated by the PSD formalism to constrain the slope of the PSD and other features in its shape. Most relevant to the current work, Caplar & Tacchella 2019 fit broken power-law models to galaxy fluctuations around the star forming sequence at $z\sim 0$ with M${}_{*}=10^{10}-10^{10.5}$M_⊙. With degeneracies due to current observational uncertainties, they find that they cannot constrain both a slope and break timescale, but find a break timescale of $\sim 200$ Myr assuming a slope of $\beta=2$. Wang & Lilly 2020a extend this analysis in a spatially resolved direction and find PSD slopes of $\beta\sim 1-2$ in the timescale range $\sim 5$ Myr to $\sim 800$ Myr, assuming no break in the PSD (which implies that SFHs are correlated over the the age of the universe). They also find that the slopes generally decrease with stellar mass for M${}_{*}>10^{9}$M_⊙, and are correlated with estimated gas depletion timescales in galaxies. Going forward, these novel measurements can be used to constrain free parameters in the different models, and existing models can be used to make predictions for future observations with upcoming facilities like JWST and WFIRST.

2. Constraints on the timescales for observed phenomena: Combining the spectral features from distant galaxies across a range of wavelengths in a full SED fitting code allows us to estimate the star formation histories of individual galaxies with uncertainties (Heavens et al. 2000; Tojeiro et al. 2007; 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). While these observationally derived SFHs are not sensitive to variability on short timescales, they can be useful for measuring the timescales for morphological transformations, mergers, quenching and rejuvenation, and even recent starbursts. These timescales can then be linked to features in the PSDs of galaxy SFHs, such as peaks or breaks.

Pacifici et al. 2016 analyse a sample of quiescent galaxies from CANDELS at $0.2<z<2.1$ and find quenching timescales to be $\sim 2-4$ Gyr, with a strong mass dependence. Carnall et al. 2018 study a sample of quiescent galaxies from UltraVISTA at $0.25<z<3.75$ and find that the majority of galaxies quench on timescales of $\sim 0.4\tau_{\mathrm{H}}$, with a rising set of galaxies towards the lower redshift portion of their observations with quenching timescales $\sim 0.6\tau_{\mathrm{H}}$. Iyer et al. 2019 analysed a sample of CANDELS galaxies at $0.5<z<3.0$ and found that $\sim 15-20\%$ of galaxies showed evidence for multiple strong episodes of star formation, with the median timescale separating multiple peaks to be $\sim 0.4\tau_{\mathrm{H}}$, which matches the predictions using cosmological simulations by Tacchella et al. 2016. The study also found that the SFHs of galaxies were correlated with their morphological classification, with an elevation in SFR on timescales over the last $\sim 0.5$ Gyr in galaxies classified as mergers and interactions, and with a longer period of SFR decline for spheroids compared to disks.

A number of studies (Lotz et al. 2011; Snyder et al. 2017; Duncan et al. 2019) also use statistical estimates of the physical properties of galaxies to constrain merger rates and observability timescales. Pandya et al. 2017 uses a similar statistical approach to quantify the timescales on which galaxies experience quenching and rejuvenation by studying the relative number of galaxies that are star forming, quiescent, and transitioning between the two states at a given epoch.

Current observational techniques require a certain set of modeling assumptions, such as a choice of IMF, stellar population synthesis (SPS) model, and dust attenuation law. Combined with state-of-the-art observations, this leads to uncertainties of $\sim 0.2$ dex in estimating stellar masses, and $\sim 0.3$ dex in estimating star formation rates from SED fitting, with fractional uncertainties in SFR growing large as we go to lower values of SFR and older stellar populations. Caution should be exercised in analysing the variability across different timescales using these derived physical properties, with care taken in propagating measurement uncertainties and instrumental effects in observations to uncertainties on their estimated physical properties. One example of this procedure is in accounting for the difference between the observed and intrinsic scatter in the SFR-M_∗ correlation due to measurement uncertainties (Kurczynski et al. 2016; Boogaard et al. 2018), which would potentially affect the PSD slope described earlier in this section.

That being said, techniques to model and extract SFH information from galaxy SEDs are growing increasingly sophisticated (Leja et al. 2019b; Iyer et al. 2019), and are (i) better at estimating the older star formation in galaxies, (ii) using fully Bayesian techniques accounting for possible covariances between parameters, and (iii) implementing well motivated priors being used to break degeneracies between parameters like dust, metallicity and SFH. With this in mind, it is hoped that in addition to timescales, SFHs from upcoming surveys will also be able to provide direct constraints on PSD slope and power on the longest timescales. Functionally, this provides a way to infer the same information as the first class of constraints on these timescales, although in this case the timescales are estimated from the histories of individual objects as opposed to recent burstiness of ensembles of galaxies. This would, in principle, allow us to independently verify estimated timescales, and test the assumption of ergodicity inherent to constraints on the PSD obtained using ensembles of galaxies.

3. Constraints on the timescales of physical processes: In addition to the two approaches described above, observations can also directly measure timescales for gas depletion (Kennicutt Jr 1998; Wong & Blitz 2002; Bigiel et al. 2008), stellar winds (Sharp & Bland-Hawthorn 2010; Ho et al. 2016), disk formation (Kobayashi et al. 2007), bulge growth (Lang et al. 2014; Tacchella et al. 2015), black hole growth (Hopkins et al. 2005) and more, albeit for limited samples of galaxies. On short timescales, a large body of work also exists studying GMC lifetimes ($\sim 10-30$ Myr) (Zanella et al. 2019; Kruijssen et al. 2019; Chevance et al. 2020), measuring the extent to which this depends on environment, and the extent to which it is decoupled from galactic dynamics. Krumholz et al. 2017 also find episodic starbursts lasting $\sim 5-10$ Myr with intervals of $\sim 20-40$ Myr in a ring around the Milky-Way’s central molecular zone. Equivalent behaviour in the zoom galaxies would therefore manifest as a local peak in the PSD on those timescales. In the local universe, resolved observations of stellar populations allow us to constrain the SFHs of nearby galaxies using colour-magnitude diagrams (Weisz et al. 2011a). For galaxies where stellar populations can be resolved, additional timescale information can be obtained from chemical abundances, since the production of heavy elements by different types of supernovae trace a range of intermediate timescales (Kobayashi et al. 2007; Kobayashi & Nomoto 2009). However, the masses of these galaxies are often too low to compare against the large-volume models considered in the current work. Another interesting study along these lines uses the fact that supernovae are produced at a certain rate after an episode of star formation, to compute the delay time distributions of SN Ia using SN Ia yields in conjunction to observationally measured galaxy SFHs (Strolger et al. 2020). These constraints on the timescales of physical processes allow for better modeling of the individual components that contribute to the full PSD of a galaxy’s SFH, and sometimes provide an independent check of behaviour predicted using the PSDs.

Using the PSD formalism as our basis, it is therefore possible to constrain the PSD power on certain timescales or the PSD slope on certain timescale ranges using observations, with currently available data already starting to provide initial estimates of PSD slopes on $\sim 4-800$ Myr timescales. Taken together, the three types of constraints outlined above, i.e., (i) estimates of SFR variability using ensembles of galaxies, (ii) observationally measured timescales for phenomena like quenching, starbursts, and rejuvenation, and (iii) timescales for physical processes like gas depletion and GMC formation and destruction will allow us to compare features in the PSD such as slopes and peaks across different galaxy populations.

In order to investigate the effects of resolution of the numerical simulation, we consider three realizations of the TNG100 simulation (TNG100-1,2 and 3), which are identical except for resolution. These are described in further detail in Pillepich et al. 2018b, and contain star particles with initial masses of $9.44\times 10^{5}$M${}_{\odot}/h$, $7.55\times 10^{6}$M${}_{\odot}/h$, and $6.04\times 10^{7}$M${}_{\odot}/h$, respectively. Mufasa and Simba therefore fall somewhere between TNG100-2 and TNG100-3 in terms of resolution, while EAGLE is comparable to TNG100-1. We show the SFHs and corresponding PSDs for these runs in Figure 17.

In SFH space, we see that the different resolutions have a large impact on SFHs across all masses, especially in the portions with low SFRs. For the three simulations, with our adopted $100$ Myr time bins, the lowest SFRs we can probe are $\approx 10^{-2.02}$, $10^{-1.12}$ and $10^{0.68}~\mathrm{M}_{\odot}\mathrm{yr}^{-1}$ neglecting mass loss. We see this in effect as the median SFHs for low-mass galaxies grow increasingly dominated by shot-noise and in the case of TNG100-3, completely drop off the plot. This resolution effect also affects high-mass quiescent galaxies, which leads to the apparent more rapid quenching of the median SFHs in the highest-mass bins. This is simply because the SFRs can only drop to their minimum value from the quantum of SFR given the resolution, leading to a steeper apparent drop in the SFHs.

In PSD space, we see that the effect of lower resolution is to increase the amount of white-noise in the PSDs, which manifests as a flattening to spectral slopes of 0 towards shorter timescales. In addition to affecting the PSDs to higher masses, the white noise also increases in magnitude proportional to the mass of the star particles, leading to contamination at longer timescales in a given stellar mass bin. This effect is quantified in the analysis of Appendix A.

While resolution can be a limiting factor in any analysis of small-scale features in hydrodynamical models, convergence tests on individual simulations (Genel et al. 2018; Keller et al. 2019). and forward modeling the effects of discrete star particles as in Appendix A allow us to understand and account for these limitations. In this case, resolution effects mostly prevent us from studying the behaviour of the PSDs on small timescales, which can be circumvented using the zoom simulations, which have much higher resolution. It should also be noted that the location of the breaks in the PSD listed in Section 3 are robust to resolution, although their strengths can be affected by the amount of white noise. Therefore, breaks and peaks in the PSD of simulations should be carefully compared to observations (see also Section 4.2).

The considerable differences between the PSD slopes and power across the different models seen in Figure 11 are caused in part due to the different modeling assumptions for physical processes, for example, AGN seeding, growth and feedback, star formation and stellar feedback and processes governing galactic winds. These differences are also in part due to the resolution of the simulations, as seen in Section 4.3 and numerical techniques used to implement gravity and magnetohydrodynamics (MHD), ranging from no explicit treatment of MHD in empirical and semi-analytic models, to differences between smoothed particle hydrodynamics, adaptive mesh refinement schemes and other hybrid techniques in hydrodynamical simulations (e.g., Vogelsberger et al. 2012; Sijacki et al. 2012; Kereš et al. 2012) using codes like AREPO (Springel 2010) in Illustris and IllustrisTNG, GIZMO (Hopkins 2014) in FIRE-2, Mufasa and Simba, a forked version of GADGET-3 (Springel et al. 2005) in EAGLE, Gasoline (Wadsley et al. 2004) for the g14 suite, and ChaNGa (Menon et al. 2015) for the Marvel/Justice League suite of zoom simulations.

While the current work serves to illustrate the cumulative differences between models due to choices of numerical techniques and physical models, it is outside the scope of the current work to break down the individual contributions. Building on the current work, there are three directions in which we can begin to better connect individual physical processes to their relevant timescales in a model-independent way:

(i) Tacchella et al. 2020 propose an analytical model in PSD space using the gas regulator model of galaxy evolution (Lilly et al. 2013), extending the model to account for the creation and destruction of GMCs on short timescales. Using this model, they derive the PSD as a broken power-law with multiple breaks that characterise the equilibrium timescale of gas inflow and the average lifetime of GMCs. Applied to PSDs from the different models we consider, this can explain the effective timescales for these processes across the various models.

(ii) In a slightly different direction, many of the models we consider have run additional simulations varying the input physics. For example, there is a set of $25$ Mpc³ boxes run for Illustris-TNG varying a single parameter per run, including stellar and black hole feedback mechanisms, galactic wind scalings, and aspects of star formation (Pillepich et al. 2018a; Nelson et al. 2018b). The Simba model contains additional simulations varying the AGN feedback model (Davé et al. 2019). Crain et al. 2015 describe model variations within the EAGLE suite varying stellar and AGN feedback. Choi et al. 2017 contains a suite of zoom simulations that are run with and without AGN feedback. The Santa Cruz SAM and other semi-analytic models are also capable of being run multiple times varying model parameters.

Using all of this data, it should be possible to characterise the effects of varying physical modeling assumptions with individual models, and use this across several models to understand the general trends and timescales for physical processes like stellar and AGN feedback and baryon cycling. However, as Pillepich et al. 2018a note, ‘the optimal choices for wind as well as black hole feedback strongly depend on the whole ensemble of galaxy formation mechanisms incorporated into the model.’ What holds for a given model need not generalise across all models, and extreme caution should be exercised while extrapolating trends from individual models, using the full available range of observational constraints described in Section 4.2 as benchmarks.

(iii) Further studies will also be needed to investigate the link between the well-studied effects of spatial turbulence on star formation (Larson 1981; Nakamura & Li 2005; Krumholz & McKee 2005; Padoan & Nordlund 2011) and the natural emergence of power-law spatial correlation functions (Guszejnov et al. 2018) to its temporal manifestations studied in this work. Studies like di Leoni et al. 2015, which look at the joint spatio-temporal power spectra for numerical simulations of turbulent flows to identify the signatures of physical mechanisms, provide a useful starting point in this regard.