# Modelling Stochastic Star Formation History of Dwarf Galaxies in GRUMPY

Yue Pan,1,2 Andrey Kravtsov1,3,4

1Department of Astronomy & Astrophysics, The University of Chicago, Chicago,
IL 60637, USA
2Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA
3Enrico Fermi Institute, The University of Chicago, Chicago, IL 60637 USA
4Kavli Institute for Cosmological Physics, The University of Chicago, Chicago,
IL 60637 USA
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We investigate the impact of bursty star formation on several galaxy scaling
relations of dwarf galaxies using the GRUMPY galaxy formation model. While
this model reproduces the star formation rate (SFR)–stellar mass, stellar
mass–gas mass, and stellar mass–metallicity relations, the scatter of these
relations in the original model is smaller than observed. We explore the
effects of additional stochasticity of SFR on the scaling relations using a
model that reproduces the level of SFR burstiness in high-resolution zoom-in
simulations. The additional SFR stochasticity increases the scatter in the
SFR–stellar mass relation to a level similar to that exhibited by most nearby
dwarf galaxies. The most extreme observed starbursting dwarfs, however,
require higher levels of SFR stochasticity. We find that bursty star formation
increases the scatter in the colour–magnitude distribution (CMD) for brighter
dwarf galaxies $(M_{V}<-12)$ to the observed level, but not for fainter ones
for which scatter remains significantly smaller than observed. This is due to
the predominant old stellar populations in these faint model galaxies and
their generally declining SFR over the past 10 Gyrs, rather than quenching
caused by reionization. We examine the possibility that the colour scatter is
due to scatter in metallicity, but show that the level of scatter required
leads to an overestimation of scatter in the metallicity–mass relation. This
illustrates that the scatter of observed scaling relations in the dwarf galaxy
regime represents a powerful constraint on the properties of their star
formation.

###### keywords:

Galaxies: formation – Galaxies: evolution – Galaxies: star formation –
Galaxies: dwarf

††pubyear: 2022††pagerange: Modelling Stochastic Star Formation History of
Dwarf Galaxies in GRUMPY–References

##  1 Introduction

Scaling relations between galaxy properties are a valuable window into
processes driving galaxy formation (Kennicutt, 1989; Kauffmann et al., 2003;
Kormendy et al., 2009; Lilly et al., 2013; Tacchella et al., 2016; Mason et
al., 2023). Over the past two decades, the relation between galaxy stellar
mass, $M_{\star}$, and its star formation rate, $\rm SFR$, was studied very
actively (see, e.g., Popesso et al., 2023, and references therein). The nearly
linear relation, often referred to as the star-forming main sequence (SFMS),
is widely considered to be instrumental to understanding how stellar mass of
galaxies is built up and the processes that regulate star formation in
galaxies (e.g., Leitner & Kravtsov, 2011; Leitner, 2012; Lee et al., 2015).

Observations revealed that the SFMS exists both at $z=0$ (e.g., Brinchmann et
al., 2004; Salim et al., 2007; Renzini & Peng, 2015; Cano-Díaz et al., 2016)
and at higher redshifts (see e.g., Noeske et al., 2007; Daddi et al., 2007;
Pannella et al., 2009; Rodighiero et al., 2010; Wuyts et al., 2011; Whitaker
et al., 2012; Mancuso et al., 2016; D’Silva et al., 2023).

In general, both the form of the SFMS (its amplitude and slope) and the
scatter of the $\log_{10}\rm SFR$ around the mean relation,
$\Delta_{\rm{MS}}$, are important. The scatter is directly related to
seemingly random variations (stochasticity) of star formation rate and can
provide valuable insights into the physical processes that affect star
formation histories (SFHs) on different mass scales (e.g., Matthee & Schaye,
2019).

For bright galaxies, $\Delta_{\rm{MS}}\approx 0.20-0.35$ dex is estimated in
both observations and simulations (e.g., Daddi et al., 2007; Whitaker et al.,
2012; Speagle et al., 2014; Torrey et al., 2014; Sparre et al., 2015; Matthee
& Schaye, 2019). Moreover, it has been found that the scatter
$\Delta_{\rm{MS}}$ of the SFMS can be accurately described by the log-normal
distribution (e.g., Guo et al., 2013; Chang et al., 2015; Schreiber et al.,
2015; Davies et al., 2019). There are also indications that scatter increases
with decreasing stellar mass of galaxies (e.g., Santini et al., 2017) and that
it is particularly prevalent in galaxies at $z\gtrsim 6$ (Dressler et al.,
2023; Looser et al., 2023; Atek et al., 2023; Endsley et al., 2023; Ciesla et
al., 2023).

Some studies have used the fact that $\Delta_{\rm MS}$ is moderate to suggest
that galaxies’ SFRs on the MS persist over long periods, regulated by various
physical processes such as gas accretion and feedback-driven outflows (Bouché
et al., 2010; Daddi et al., 2010; Genzel et al., 2010; Davé et al., 2011;
Lilly et al., 2013; Dekel & Mandelker, 2014; Rodríguez-Puebla et al., 2016;
Tacchella et al., 2016). However, there is no consensus on the redshift
evolution of $\Delta_{\rm{MS}}$. While some studies concluded that it is
independent of redshift (e.g., Daddi et al., 2007; Noeske et al., 2007;
Whitaker et al., 2012; Speagle et al., 2014; Ciesla et al., 2014; Pessa et
al., 2021), others find that it evolves with redshift (e.g., Kurczynski et
al., 2016; Santini et al., 2017; Katsianis et al., 2019; Tacchella et al.,
2020; Davies et al., 2022; Shin et al., 2023).

Building on these observations and models, Caplar & Tacchella (2019) utilized
the power spectral density (PSD) distribution to characterize the
stochasticity of SFHs for all star-forming galaxies with $M_{\star}\approx
10^{10}M_{\odot}$. They found that the SFHs of galaxies tend to be correlated
on a timescale of $\approx 200$ million years. This timescale is shorter than
the dynamical time of the dark matter halo, implying that the baryonic
effects, such as feedback and reincorporation of galactic winds, which operate
on the dynamical timescale of galaxies, play a crucial role in shaping the
SFHs of galaxies. In several follow-up studies, the PSD was used to
investigate the effect of gas cycling through ISM and giant molecular clouds
(GMCs) on the SFR of galaxies (Tacchella et al., 2020), analysis of the level
of SFR stochasticity in different galaxy formation models and simulations
(Iyer et al., 2020), and the impact of different observational SFR tracers on
the inferred SFR stochasticity (Iyer et al., 2022). These studies elucidated
the contribution of different physical processes to the variation of SFR on
different time scales.

On short time scales ($\lesssim 100$ Myr), the amplitude of the SFR
fluctuations contains valuable information about the creation and destruction
of individual GMCs (Scalo & Struck-Marcell, 1984, 1986; Kruijssen & Longmore,
2014; Krumholz & Kruijssen, 2015; Semenov et al., 2017; Faucher-Giguère, 2018;
Orr et al., 2019) caused by supernova explosions, cosmic rays, and
photoionization (Gnedin et al., 2008; Parrish et al., 2009; Hopkins et al.,
2014; Faucher-Giguère, 2018; Tacchella et al., 2020; Semenov et al., 2021b;
Semenov et al., 2021a). In the regime of extremely low SFR in dwarf galaxies,
the stochasticity can even be caused by the formation of individual massive
stars (see e.g. Fumagalli et al., 2011; da Silva et al., 2012, 2014).

On intermediate time scales of $\sim 0.1-1$ Gyr, various dynamic processes
such as galactic winds, disc formation, bulge growth, environmental effects,
bar-induced inflows, and galaxy mergers are thought to influence the
fluctuations in star formation (Oppenheimer & Davé, 2008; Tacchella et al.,
2016; Sparre et al., 2017; Torrey et al., 2018; Wang & Lilly, 2020a, b).

On longer time scales ($>1$Gyr), galaxy quenching and the overall interplay
between gas accretion and outflows are known to play significant roles
(Behroozi et al., 2013; Moster et al., 2018; Birrer et al., 2014; Rodríguez-
Puebla et al., 2016; Weinberger et al., 2017; Anglés-Alcázar et al., 2017;
Tacchella et al., 2018; Behroozi et al., 2019).

Although the stochasticity of SFR in massive galaxies was explored
extensively, our understanding of the scatter in the SFMS in the dwarf galaxy
regime remains rather limited, even though studies of star formation in nearby
dwarf galaxies indicate that their SFRs are bursty (e.g., Emami et al., 2019).
While the approximately linear relation between SFR and stellar mass extends
to galaxies of stellar mass as low as $\approx 10^{5}M_{\odot}$, the scatter
growth dramatically with decreasing stellar mass for galaxies with
$M_{\star}\lesssim 10^{8}\,M_{\odot}$ (see, e.g., Kravtsov & Manwadkar, 2022).
In particular, the latter authors showed that scatter stemming from the
variations of mass assembly histories and structural properties of dwarf
galaxies is not sufficient to explain the observed scatter of specific SFRs in
the dwarf galaxy regime.

This is not entirely surprising, given that star formation in dwarf galaxies
is modulated by a few individual star forming regions, which leads to burstier
star formation. At the same time, outflows are expected to be more efficient
in dwarf galaxies which may result in periods of low SFR following bursts
thereby leading to stronger variations of SFR. These processes are confined
within localized regions in galaxies and require detailed modeling of
interstellar medium and star-forming regions in galaxies, as well as their
destruction and accompanying feedback effects.

However, many galaxy formation models and simulations do not model or resolve
the relevant scales and model relevant processes. For example, Kravtsov &
Manwadkar (2022) show that a simple regulator-based GRUMPY model reproduces
the mean observed correlations between properties of dwarf galaxies, as well
as luminosity function and radial distribution of Milky Way’s dwarf satellites
(Manwadkar & Kravtsov, 2022). However, the model does not model small-scale
star-forming regions and thus predicts a much smaller scatter than observed in
the colour–magnitude and stellar mass–SFR relations. Similarly, Pan et al.
(2023) found that dwarf galaxy colours in the Auriga cosmological
hydrodynamical simulations (Grand et al., 2017, 2021) have a significantly
smaller scatter than observed. This is likely due to the fact that ISM in
these simulations is pressurized and ISM and formation and destruction of
individual star-forming processes are not modeled, but overall star formation
is modeled statistically with a subgrid model.

The goal of this study is to evaluate the degree of SFR stochasticity in dwarf
galaxies required to reproduce observed scatter in the colour-magnitude and
stellar mass–SFR relations. To this end, we use the regulator-based GRUMPY
galaxy formation model mentioned above, but add scatter to the SFR predicted
by the model in a controlled manner.

This paper is organized as follows. We describe the GRUMPY model in Section 2.
We summarize the PSD formalism used in the model of SFR stochasticity in
Section 2.2 and the revised model with stochasticity in Section 2.3. We
investigate the effects of stochasticity on the SFR–mass relation in Section
3.1, the colour–magnitude diagram in Section 3.2, and the metallicity–stellar
mass relation in Section 3.3. In Section 4, we compare our results with
results from other theoretical models and simulations in Section 4.1and 4.2,
explore the sensitivity of different SFR indicators to the stochasticity of
SFR in observations in Section 4.3, and discuss several caveats in Section
4.4. Finally, we present our conclusions in Section 5.

The cosmological parameters adopted in this study are those of the Caterpillar
simulation suite: $h=H_{0}/100=0.6711$, $\Omega_{\rm m0}=0.32$,
$\Omega_{\Lambda}=0.68$.

##  2 Galaxy formation model with stochastic SFR

We begin with description of how stochasticity is modeled in this study in
conjunction with the SFR produced by the fiducial GRUMPY model (Kravtsov &
Manwadkar, 2022). We first review the GRUMPY model with a focus on the
modeling of star formation, molecular hydrogen, and metallicity – quantities
most relevant for the present study. We then describe the stochastic SFR
framework based on the power spectral density approach of Tacchella et al.
(2018). Finally, we combine the two components and describe how stochastic SFR
is modelled in this study.

###  2.1 The GRUMPY galaxy formation model

####  2.1.1 Halo mass tracks and GRUMPY galaxy formation model

We model the population of dwarf galaxies around the Milky Way using tracks of
haloes from the Caterpillar suite of $N$-body simulations (Griffen et al.,
2016), which simulated 32 MW-sized haloes.111The Caterpillar project website
can be found at https://www.caterpillarproject.org. We use the highest
resolution set of simulations, LX14, which allows to maximize the size of our
model galaxy sample.

The halo tracks were produced by a modified version of the Rockstar halo
finder and the Consistent Trees Code (Behroozi et al., 2013), as described in
Section 2.5 of Griffen et al. (2016). As shown by Manwadkar & Kravtsov (2022,
see their Fig. 1), the subhalo peak mass function in the LX14 simulations is
complete at $\mu=M_{\rm peak}/M_{\rm host}\gtrsim 4\times 10^{-6}$ or $M_{\rm
peak}\approx 4\times 10^{6}M_{\odot}$ for the host halo mass $M_{\rm
host}\approx 10^{12}\,M_{\odot}$, even in the innermost regions of the host
($r<50$ kpc). This degree of completeness is sufficient to model the full
range of luminosities observed in Milky Way satellites since the faintest
ultra-faint dwarfs (UFDs) are hosted in haloes of $M_{\rm peak}\gtrsim
10^{7}\,M_{\odot}$ in our model (Kravtsov & Manwadkar, 2022; Manwadkar &
Kravtsov, 2022).

The halo mass evolution tracks from the Caterpillar suite are used as the
basis for the GRUMPY galaxy formation model Kravtsov.Manwadkar.2022 based on a
regulator-type galaxy formation framework (e.g., Krumholz & Dekel, 2012; Lilly
et al., 2013; Feldmann, 2013). Namely, the model solves a system of coupled
differential equations to follow the evolution of key galaxy properties,
including the effects of UV heating after reionization, gas accretion
suppression onto small mass haloes, galactic outflows, a model for the gaseous
disk and its size, molecular hydrogen mass, star formation, and more. The
model also accounts for the evolution of the half-mass radius of the stellar
distribution. The galaxy model parameters used in this study are identical to
those in the fiducial model of Manwadkar & Kravtsov (2022).

As the primary of this study is on the stochasticity of star formation rate,
we briefly review how star formation is modeled in the GRUMPY model. The model
uses the evolved exponential radial gas profile, $\Sigma_{g}(R)$, to estimate
the mass of molecular hydrogen using the model of Gnedin & Draine (2014):

$M_{\rm
H_{2}}=2\pi\int_{0}^{\rm{R_{HI}}}f_{\rm{H_{2}}}(\Sigma_{g})\Sigma_{g}(R)RdR$
(1)

where $R_{\rm{HI}}$ is the radius corresponding to the assumed self-shielding
threshold and $f_{\rm H_{2}}$ is molecular fraction computed using current gas
metallicity and mean star formation rate.

The star formation rate is computed using $M_{\rm{H_{2}}}$ assuming a constant
depletion time of molecular gas, $\tau_{\rm{dep,H_{2}}}$, and instantaneous
recycling of the gas:

$\dot{\mathcal{M}}_{\star}=(1-\mathcal{R})\dot{M}_{\star}$ (2)

where

$\dot{M}_{\star}=\frac{M_{\rm{H_{2}}}}{\tau_{\rm{dep,H_{2}}}}=\rm{SFR}$ (3)

$\mathcal{R}=0.44$ represents the proportion of gas that is converted into
stars and subsequently returned to the ISM assuming the Chabrier (2003)
initial mass function of stars (e.g., Leitner & Kravtsov, 2011; Vincenzo et
al., 2016). Meanwhile, $\tau_{\rm{dep,H_{2}}}$ is assumed to be
$\tau_{\rm{dep,H_{2}}}=2.5$ Gyr in the fiducial model consistent with
measurements in nearby galaxies (e.g., Bigiel et al., 2008; Bigiel & Blitz,
2012; Bolatto et al., 2011; Rahman et al., 2012; Leroy et al., 2013).

Observations indicate that $\tau_{\rm H_{2}}$ may vary between galaxies and
with redshift. Thus, the assumption of constant $\tau_{\rm H_{2}}$ is one of
the reasons the fiducial model may underestimate SFR stochasticity.

Kravtsov & Manwadkar (2022) and Manwadkar & Kravtsov (2022) showed that the
model successfully reproduces the general trends observed in the SFR–mass
relation, mass–metallicity relation, and colour–magnitude relation of Milky
Way satellite galaxies, as well as luminosity function and radial distribution
of the Milky Way satellites. However, as noted in the Introduction, the model
predicts a significantly smaller scatter in star formation rates and colours
of dwarf galaxies than what is actually observed in the Universe.

This is not surprising because the GRUMPY framework does not include modeling
of the processes that can result in significant stochasticity of star
formation, such as formation and destruction of individual star-forming
regions (e.g., Iyer et al., 2020). Here we introduce additional SFR
stochasticity in a controlled manner using PSD formalism, as described below.
However, we first describe the sources of stochasticity that are included in
the model.

####  2.1.2 Sources of stochasticity in the original model

The fiducial GRUMPY model includes two sources of stochasticity. The first
arises from the variations of the mass assembly histories (MAHs) of haloes of
a given final mass. The diversity in halo MAHs can result in a range of
scatter in galaxy properties and can be particularly significant in the ultra-
faint dwarf (UFD) regime due to the suppression of gas accretion and star
formation after the epoch of reionization (Rey et al., 2019; Katz et al.,
2020).

The second source of stochasticity arises from the assumed stochasticity in
the size of gas distribution in galaxies. The latter is modeled using an
exponential profile:

$\Sigma_{g}(R)=\Sigma_{0}\exp\left(-\frac{R}{R_{d}}\right)$ (4)

where $\Sigma_{0}$ is the central surface density, $R_{d}$ is the disc size.
The model assumes that at each redshift $z$ during evolution, the disc size is
proportional to the halo radius enclosing density contrast of 200 times the
critical density, $R_{\rm{200c}}(z)$:

$R_{d}=\chi_{d}R_{\rm{200c}}(z)$ (5)

The value of $\chi_{d}$ has a mean of 0.06 (Mo et al., 1998) and a random log-
normal Gaussian scatter of 0.25 dex to account for the fact that there is
significant scatter in disc sizes at a given stellar mass in observations. The
same random $\chi_{d}$ is used throughout its evolution for a given galaxy,
but different galaxies have different $\chi_{d}$ values.

Although these sources of stochasticity between galaxies produce some scatter
in galaxy scaling relations and in particular reproduce scatter in observed
size–magnitude relation (Kravtsov & Manwadkar, 2022; Manwadkar & Kravtsov,
2022), the scatter in the $M_{\star}-\rm SFR$ and colour–magnitude relations
is considerably smaller than in the corresponding observed relations, hence
motivating additional sources of SFR scatter.

###  2.2 Power Spectral Density of stochastic SFR

We follow Caplar & Tacchella (2019) and model the specific star formation rate
(SFR) of a galaxy with stellar mass $M_{\star}$ relative to the main sequence
(MS) relation as:

$\Delta_{\rm{MS}}=\log_{10}\left(\frac{\rm{SFR}}{\langle\rm{{SFR}}\rangle(M_{\star}),}\right)$
(6)

where $\langle\rm SFR\rangle(M_{\star})$ represents the SFR of the mean MS at
the mass $M_{\star}$ and $\Delta_{\rm{MS}}$ is a Gaussian random number.
Consequently, the distribution of $\rm{SFR}/\langle\rm SFR\rangle(M_{\star})$
is log-normal, which is consistent with observations (Guo et al., 2013; Davies
et al., 2019; Caplar & Tacchella, 2019).

Figure 1: An illustration of how $\Delta_{\rm{MS}}$ and its PSD depend on
$\alpha$, $\tau_{\rm{break}}$, and $\sigma_{\rm{int}}$. We show $\alpha=1,2,3$
and $\tau_{\rm{break}}=0.1,1$ Gyr. In each panel we show
$\sigma_{\rm{int}}=0.2$ (grey) and 0.6 (purple). As $\alpha$ increases and
$\tau_{\rm{break}}$ increases, adjacent points become more correlated. A
larger $\sigma_{\rm{int}}$ value elevates the PSD on all timescales.

We thus model stochasticity of star formation relative to the mean
$M_{\star}-{\rm SFR}$ relation by producing correlated Gaussian random number
for each time step during galaxy evolution integration. The numbers are
generated assuming a given power spectral density (PSD) of the time series in
the Fourier domain.

The time series, represented by a discrete set of values $s(t_{n})=s_{n}$ at
evenly spaced intervals $t_{n}=t_{\rm min}+n\Delta t$ for $n=[0,N-1]$, covers
a time span of $T=N\Delta t$. The Fourier transform of the series, $S_{k}$,
where $k$ represents the wave number and ranges from $-N/2$ to $N/2$:

$\displaystyle s_{n}$
$\displaystyle=\frac{1}{N}\sum\limits_{k=0}^{N-1}S_{k}\cdot\exp\left(\frac{i2\pi
n}{N}k\right)$ (7)
$\displaystyle=\frac{1}{N}\sum\limits_{k=0}^{N-1}S_{k}\cdot\exp\left[\frac{i2\pi(t_{n}-t_{\rm
min})}{\Delta tN}k\right]$
$\displaystyle=\frac{1}{N}\sum\limits_{k=0}^{N-1}S_{k}\cdot\exp\left[i2\pi(t_{n}-t_{\rm
min})f_{k}\right]$

The DC signal has a frequency $f_{0}=0$ and a corresponding $S_{0}$ value of
1. The frequency $f_{k}$ is defined as $k/T$.

According to the Parseval’s theorem, for a real series $s_{n}$:

$\frac{1}{N}\sum\limits_{n=0}^{N-1}s^{2}_{n}=\frac{1}{N^{2}}\sum\limits_{k=0}^{N-1}|S_{k}|^{2}$
(8)

For a series $s_{n}$ with a mean of zero, the left-hand side represents the
variance, and the right-hand side indicates that the sum of the squares of
Fourier components of $s_{n}$ divided by $N^{2}$ should equal the variance of
$s_{n}$.

The power spectrum $P(k)$ is defined as:

$P(k)=\frac{1}{N^{2}}|S_{k}|^{2}$ (9)

Note that $P(k)$ is dimensionless, and according to the Parseval’s theorem
above:

$\sigma^{2}(s_{n})=\sum\limits_{k=0}^{N-1}P(k)$ (10)

The power spectral density (PSD) is defined as

${\rm PSD}(k)=\frac{(\Delta
t)^{2}}{T}|S_{k}|^{2}=\frac{T}{N^{2}}N^{2}P(k)=TP(k)=N\Delta tP(k).$ (11)

and has unit of time.

To generate a sequence of Gaussian random numbers with a given power spectrum,
we use a grid of Fourier components at frequencies $f_{\rm grid}$ with
amplitudes set to random numbers drawn from the Gaussian pdf with zero mean
and unit variance and multiplied by $\sqrt{P(k)}=\sqrt{{\rm PSD}(f)/T}$.

We follow (Caplar & Tacchella, 2019) and use PSD of the form

$\rm{PSD}(\emph{f})=\frac{\sigma^{2}}{1+(\tau_{\rm{break}}\
\emph{f}_{\rm{grid}})^{-\alpha}},$ (12)

where $\sigma=\sigma_{\rm{int}}\tau_{\rm{decor}}$ characterizes the
variability over a long time scale and
$\tau_{\rm{break}}=2\pi\tau_{\rm{decor}}$ characterizes the timescale over
which the random numbers are no longer correlated.

We generate a time series of $\Delta_{\rm MS}$ perturbations this way with
$\Delta t=10\,\rm Myr$ – the value chosen to correspond to the typical
lifetime of star forming regions in galaxies (Fall et al., 2010; Kruijssen et
al., 2019; Semenov et al., 2017; Tacchella et al., 2020; Semenov et al.,
2021b).

Figure 1 shows an illustration of how $\Delta_{\rm{MS}}$ and its PSD changes
with $\alpha$, $\tau_{\rm{break}}$, and $\sigma_{\rm{int}}$. In general, a
larger $\alpha$ and $\tau_{\rm{break}}$ increases the correlation between
adjacent points, and a larger $\sigma_{\rm{int}}$ increases the power on all
timescales.

Figure 2: An illustration of non-stochastic SFH, stochastic SFH, and
stochastic+perturbed SFH (left) along with their corresponding PSD (right).
The solid green line represents the SFH of a typical halo track that exhibits
high star formation during the first 2 Gyr, followed by a steady SFR decline.
Due to the rather smooth SFH, the power on short-term fluctuations (thick
green line on the right) is low. To create a more bursty SFH, we add
stochasticity to the SFH of this sample galaxy, following the process
described in Section 2, which yields a red curve on the left. Since we
manually introduce short-term fluctuations, the power on short-term
fluctuations is much greater than in the non-stochastic case (red curve on the
right compared to green curve). Moreover, as shown in Section 3.1, we test
another model where we introduce an additional PSD corresponding to a delta
function on the timescale of $\sim 13$ Myr that increases SFR by a random
number. This is a test to show how much additional perturbation is needed to
match the high SFR of starburst galaxies like those in Lin et al. (2022). As
demonstrated, adding such a PSD results in significantly greater power on
short fluctuation timescales (blue curve on the right). Figure 3:
$M_{\star}-\rm SFR$ (left panel) and $M_{\rm g}-\rm SFR$ (right panel)
relations for model galaxies in the non-stochastic run (green) and stochastic
run (red) as a comparison to observations such as the ALFALFA survey (Durbala
et al., 2020, grey dots), STARBIRDS survey (McQuinn et al., 2019, purple
circles), McGaugh et al. (2017, pentagons), Jimmy et al. (2015, squares), the
SHIELD I & II surveys (McQuinn et al., 2015a, McQuinn et al. 2022, submitted;
pink circles), and Lin et al. (2022, diamonds). For Lin et al. (2022), we only
include points that have SFR fractional error $<20$%. The median SFR values in
each mass bin for the non-stochastic run (green), stochastic run (red), and
stochastic+post-process perturbation (blue) are shown by solid lines. The
transparent colour band illustrates the 2-sigma variation, and any points
outside the 2-sigma width of the median are plotted as individual points. We
find that the linear relation between SFR and $M_{\star}$ extends down to
$M_{\star}\sim 10^{5}M_{\odot}$ dwarf galaxies. The 2-sigma width of the
SFR-$M_{\star}$ and SFR-$M_{g}$ relation in the non-stochastic run is
relatively smaller than that for the stochastic run. When we add additional
perturbations to the model SFR using another form of PSD as described in
Section 3.1, we obtain a much higher SFR (blue) that is within the same range
as observed galaxies in SHIELD I & II, although we are unable to reach the
same high SFRs as those in Lin et al. (2022). As discussed in Section 4,
different SFR indicators can cause significant differences in the SFR value,
with indicators such as H$\alpha$ used in Lin et al. (2022) tending to produce
higher values. Overall, our model results agree well with observations, with
future studies warranted to explore the physics that cause star bursts. Figure
4: $M_{g}$-$M_{\star}$ relation for model galaxies in the non-stochastic and
stochastic runs compared with observations. Similar to Figure 3, we plot the
median SFR for each mass bin and the 1-sigma and 2-sigma width for the
stochastic run. Switching from a non-stochastic to a stochastic SFR model does
not significantly alter the median relationship between gas and stellar mass.

###  2.3 Modelling stochastic SFR in GRUMPY

In what follows, we use the PSD given by eq. 12 with the values $\alpha=2$,
$\tau_{\rm{break}}=2$ Gyr, and $\sigma_{\rm{int}}=0.6$. We find that using the
PSD given by equation 12 with these parameter values produces results similar
to those of the EAGLE, Illustris-TNG, and FIRE-II simulations (Iyer et al.,
2020). An example of a galaxy SFH and corresponding PSD in such a model is
presented in Figure 2. This validates our approach to generating stochasticity
and demonstrates that a simple analytical model can produce SFR similar to the
SFR in state-of-the-art galaxy formation simulations.

For each individual galaxy, we generate a sequence of $\Delta_{\rm MS}$
Gaussian random numbers spaced by $\Delta t=10\,\rm Myr$, as described in the
previous section. This sequence is interpolated to a given $t$ during the
integration of galaxy evolution. To account for the fact that the mean of a
log-normal distribution is biased high compared to the unperturbed value, we
adjust the mean of the $\Delta_{\rm{MS}}$ by subtracting the variance of the
sequence divided by 2.

Stochasticity is then introduced by perturbing the SFR outputs with
$\Delta_{\rm{MS}}$:

$\rm{SFR_{stoch}}=10^{\Delta_{\rm{MS}}}\cdot\rm{SFR}$ (13)

We assume that the PSD does not change with redshift, since observed values of
$\Delta_{\text{MS}}$ depend only mildly on redshift. We tested that varying
$\Delta_{\text{MS}}$ by an amount significantly larger than the observed
redshift dependence does not change our results.

Given that galactic outflows in the model are assumed to be proportional to
SFR, $\dot{M}_{\rm out}=\eta\dot{M}_{\star}$ stochasticity should also induce
stochasticity of the outflow rate. Physically, however, outflows may not
necessarily respond to each local burst of star formation in a given star-
forming region, but should arise from a collective effect of such local
bursts. We thus assume that the outflow rate is proportional to the stochastic
SFR rate averaged over 100 Myr time scale:

$\dot{M}_{\rm{g,out}}=\eta_{w}\dot{M}_{\star,\rm{ave}},$ (14)

where $\eta_{w}$ is the mass-loading factor, and $\dot{M}_{\star,\rm{ave}}$ is
the running average of the instantaneous SFR over the past 100 Myr.

Figure 5: _Upper_ : Colour–magnitude diagram of g-r colour vs $M_{r}$ relation
for model galaxies in the non-stochastic (upper left) and stochastic runs
(upper right) compared to observations such as the MegaCam survey Muñoz et al.
(2018), the SAGA survey (Geha et al., 2017; Mao et al., 2021), and the ELVES
survey (Carlsten et al., 2020, 2022). Similar to Figure 3, we plot the median
SFR for each mass bin (thick lines) and the 1-sigma and 2-sigma width with the
2-sigma outliers (green dots for the left two columns and red dots for the
right two columns). We also show the distribution of dwarf galaxies in the
GAMA survey as a 2D histogram in grey. The non-stochastic model exhibits a
tightly clustered and almost linear colour-magnitude relation across
magnitudes, which contrasts with observations that display more scatter in the
$-20<M_{r}<-8$ magnitude range. The addition of stochastic star formation rate
(SFR) increases the scatter in CMD, particularly for brighter magnitudes
$M_{r}<-12$, but improves the comparability with observations for galaxies
brighter than $M_{r}=-12$, while still falling short of reproducing the
comparable scatter in faint dwarf galaxies. _Lower:_ the $r$-band effective
surface brightness as a function of $M_{r}$ for both non-stochastic and
stochastic runs, with negligible differences in scatter between the two cases,
and both able to reproduce the mean trend and scatter of the relation similar
to observations. Further discussion on the findings is presented in Section 4.
Figure 6: Colour–magnitude diagram of g-r colour vs $M_{r}$ relation for model
galaxies in the non-stochastic (left), stochastic SFR (middle), and stochastic
SFR + stochastic metallicity runs (right) compared to the same observations as
in Figure 5. Although perturbing the SFR does not increase the scatter in
colour, especially for faint dwarf galaxies with $M_{r}>-12$, perturbing the
metallicity significantly increases the scatter in colour across all
magnitudes. This aligns with observations that show some quenched galaxies
having blue colours and UV emission (Greene et al., 2023), where stochasticity
in metallicity likely plays a role in their blue colour instead of star
formation. Figure 7: The gas metallicity–mass (right) and stellar
metallicity–mass (left) relations for the non-stochastic and stochastic runs
with observations from various sources such as Lee et al. (2006); Berg et al.
(2016); Kojima et al. (2020); McQuinn et al. (2015b, 2020); Kirby et al.
(2013); Simon (2019). Similar to Figure 3, we plot the median SFR for each
mass bin and the 1-sigma width. In the non-stochastic case, the model
accurately predicts the mean trend of the metallicity-mass relation, but the
scatter is tighter than observed. The stochastic run shows slightly increased
scatter, but it maintains the metallicity–mass relation without any breaks.
The introduction of stochasticity in the SF model is seen to be effective in
preserving the metallicity-mass relation.

##  3 Effect of stochasticity on galaxy scaling relations

In this section, we explore effects of SFR stochasticity implemented described
in the previous section on the $M_{\star}-{\rm SFR}$, $M_{\rm g}-\rm SFR$
relations (Section 3.1), and on the colour-magnitude relation (Section 3.2).

###  3.1 Effect of stochasticity on galaxy SFR

Figure 3 shows a comparison of the $M_{\star}-\rm SFR$ and $M_{\rm g}-\rm SFR$
relation for model galaxies and observed galaxies from the ALFALFA survey
Durbala et al. (2020) and a number of other samples (Jimmy et al., 2015; Teich
et al., 2016; McGaugh et al., 2017; James et al., 2017; McQuinn et al., 2015a;
McQuinn et al., 2019, 2021; Lin et al., 2022). Compared to results shown in
Figure 11 of Kravtsov & Manwadkar (2022) the widths of $M_{\star}-\rm SFR$ and
$M_{\rm g}-\rm SFR$ relations are considerably larger for the model with added
SFR stochasticity. However, the increase is not sufficient to account for
extreme starbursting dwarfs reported by Lin et al. (2022) or even for galaxies
with the highest SFR in other samples (e.g., SHIELD and STARBIRDS galaxies).
This indicates that additional physical processes might be at play to boost
the SFR to such a high level. One possibility is that SFR in the smallest
dwarf galaxies in these samples is dominated by the formation and destruction
of one or a couple of massive star clusters, that temporarily boost SFR by a
large factor (see, e.g., Zick et al., 2018).

We explore this possibility by adding an additional stochasticity to our model
on the time scale of 10 Myr corresponding to a typical lifetime of a star-
forming region that gives rise to individual clusters. This corresponds to
adding the delta function to the PSD at this time scale. The amplitude of the
delta function corresponds to the perturbation of the actual SFR of

$\rm{SFR_{\rm{pert.}}}=\rm{SFR_{\rm{stoch.}}}\cdot r\cdot 20$ (15)

where $r$ is a random variable drawn from the exponential distribution
$3e^{-3x}$.

The result is shown by the blue line and blue shaded region in Figure 3. Such
additional stochasticity allows the model to match most observed galaxies,
although not the most extreme starbursts. Such starbursts are very rare in
observations, however, and only a dozen is found in the entire Dark Energy
Survey.

Figure 4 shows a comparison between non-stochastic and stochastic runs for the
$M_{g}-M_{\star}$ relation. Introducing stochasticity does not change the mean
trend and scatter of this relation.

In Section 4.3, we delve into the sensitivity of the stochasticity of SFR to
the choice of SFR indicators, as demonstrated by previous studies (Lee et al.,
2009; da Silva et al., 2014; Emami et al., 2019).

###  3.2 Effect of stochasticity on galaxy colours

Figure 5 illustrates that the scatter in the $g-r$ colours for bright galaxies
($M_{r}<-12$) modeled with additional SFR stochasticity agrees well with
observation, whereas the scatter for fainter galaxies ($M_{r}>-12$) has a
scatter much smaller than observed. This discrepancy may be due to the fact
that dwarf galaxies in the model remain quenched with after reionization and
thus stochasticity does not affect their magnitudes. Observed dwarf galaxies
with masses $M_{\star}\gtrsim 10^{5}\,M_{\odot}$ often do show indications of
recent star formation or star formation at intermediate epochs (e.g., McQuinn
et al., 2023; Jones et al., 2023), and such rejuvenated dwarf galaxies are
also present at least in some models of galaxy formation (e.g., Rey et al.,
2020).

To test the effect of reionization on the quenching of dwarf galaxies in the
model, we lowered the critical mass. $M_{c}(z)$, at which baryonic fraction is
suppressed by a factor of 2 relative to the universal value due to post-
reionization UV heating of the intergalactic medium (IGM). Thus, lowering
$M_{c}(z)$ results in more low-mass haloes that are able to accrete gas and
prolong star formation after the epoch of reionization. However, lowering
$M_{c}(z)$ by a factor of four does not produce significantly more blue low-
mass dwarf galaxies. Changing the form of $M_{c}(z)$ to a step function with a
filtering halo mass as low as $10^{8}M_{\odot}$ also does not change the
colour distribution for faint dwarf galaxies. Thus, a declining SFH and lack
of scatter in colour for faint dwarf galaxies is not due to accretion
suppression resulting from the post-reionization UV heating of the IGM.
Rather, it must be due to a reason internal to the galaxies themselves.

One possibility is that the optical light of observed dwarf galaxies includes
emission lines arising from ionized HII regions, while these are not included
in calculations of colours of the model galaxies. Indeed, Martín-Manjón et al.
(2008) demonstrated that the presence of emission lines from the interstellar
medium (ISM) can significantly impact the colours of starbursting dwarf
galaxies (see their Figures 7 and 10). To assess this effect, we employed the
Flexible Stellar Population Synthesis (FSPS) code (Conroy et al., 2009; Conroy
& Gunn, 2010) in conjunction with its Python bindings, PyFSPS
222https://github.com/dfm/python-fsps, to generate single stellar population
instances. We compared the colours of these instances with and without
incorporating the contribution of nebular emission lines and found only minor
differences. Therefore, the difference in the scatter of colours is unlikely
to be due to the omission of nebular emission lines in the colour calculation.

Recently, Greene et al. (2023) showed that many early-type dwarf galaxies
($M_{\star}\lesssim 10^{7}M_{\odot}$) exhibit excess UV emissions and blue
colours. These blue colours are not due to ongoing star formation, but due to
significant variations of metallicity between galaxies of a given $M_{\star}$.
Lower metallicity galaxies have bluer colours and scatter in metallicity thus
should result in the scatter of colours. To test this, we perturbed each model
galaxy’s metallicity $Z$ by a log-normal Gaussian random number:

$Z_{\rm{stoch}}=Z\cdot 10^{\Delta_{\rm{Z}}},$ (16)

where $\Delta_{\rm{Z}}$ is a Gaussian random number centered on 0 with
standard deviation of $\sigma=0.5$.

The effect of such additional metallicity scatter on colours is shown in the
right panel of Figure 6. Comparison of the colour–magnitude relations in the
three panels of this figure shows that additional scatter in metallicity
indeed leads to a significant increase in scatter of colour, especially for
dwarf galaxies with $M_{r}>-12$. Thus, the tail of blue colours of quenched
low-mass dwarf galaxies in observations could be due to metallicity scatter.
At the same time, as we show below, inclusion of such large metallicity
scatter overestimates the scatter in the metallicity–stellar mass relation.

###  3.3 Effect of SFR stochasticity on the stellar mass-metallicity relation

Given the potentially significant effects of metallicity scatter on dwarf
galaxy colours we investigate the impact of stochastic SFH on the stellar
mass–metallicity relation of dwarf galaxies, which is widely utilized as a
diagnostic tool for modeling galaxy feedback, chemical evolution, inflows, and
outflows (e.g. Garnett, 2002; Peeples & Shankar, 2011; De Lucia et al., 2020;
Finlator & Davé, 2008). In the case of dwarf galaxies, observations have been
conducted to measure stellar metallicity using absorption lines in stellar
emission (see, for instance, the review by Simon, 2019) and to determine gas
metallicity using emission lines from HII regions in galaxies with gas and
ongoing star formation (Lee et al., 2006; Berg et al., 2012; Berg et al.,
2016).

Figure 7 illustrates the relation between galaxy stellar mass and gas and
stellar phase metallicities, $Z_{\rm g}$ and $Z_{\star}$, of model and
observed dwarf galaxies. The scatter observed in the $Z_{\star}-M_{\star}$ and
$Z_{g}-M_{g}$ relations is comparable or slightly higher in the stochastic run
when compared to the non-stochastic run.

The width of the 1-sigma band indicates that the overall scatter in the gas
metallicity–gas mass relationship is higher than the scatter in the stellar
metallicity–stellar mass relationship. This can be attributed to the fact that
$Z_{\star}$ is the average metallicity in the cumulative stellar mass and
increases monotonically, whereas $Z_{\rm g}$ is influenced by star formation
and metal enrichment outflows and inflows, which govern the instantaneous gas
mass (Finlator & Davé, 2008; Lilly et al., 2013; Torrey et al., 2019; van Loon
et al., 2021). It is worth noting that the model $Z_{\rm g}$ and $Z_{\star}$
values do not consider the observational uncertainties that contribute to the
scatter of the relations of observed galaxies.

Figure 8, on the other hand, shows that additional stochasticity in the
metallicity explored in the previous subsection (eq. 8) results in a
metallicity scatter that is significantly larger than observed.

Figure 8: The stellar metallicity–mass relation for the non-stochastic and
stochastic SFR + stochastic metallicity runs compared with the same
observations as in Figure 7. Similar to Figure 7, we plot the median
$Z_{\star}$ for each mass bin and the 1-sigma and 2-sigma width. The
stochastic SFR + stochastic metallicity run shows an overestimated scatter
compared to observations.

##  4 Discussion

Results presented above indicate that star formation rates in dwarf galaxies
are very stochastic. Furthermore, we find that adding SFR stochasticity to the
SFHs of model galaxies results in increased $g-r$ colour scatter in bright
dwarf galaxies $(M_{r}<-12)$ bringing their colour distribution in a better
agreement with observations. In contrast, colours of faint dwarf galaxies
$(M_{r}>-12)$ are largely unaffected by the addition of the SFR stochasticity.

In this section, we discuss interpretation of these results, compare our
results with those of other theoretical models and observational studies,
examine how observed SFR stochasticity depends on the use of different SFR
indicators, and investigate the link between short-term and long-term
fluctuations and galaxy evolution. Lastly, we highlight the limitations of
interpretation of our results and outline potential avenues for future
research.

###  4.1 Comparison with other models

It is interesting to compare the level of SFR stochasticity required to match
observed scatter in the $M_{\star}-\rm SFR$ relation in the dwarf galaxy
regime with stochasticity in various galaxy formation models and simulations
presented by Iyer et al. (2020). These authors computed the PSD of galaxies
with $10^{9}M_{\odot}<M_{\star}<10^{11.5}M_{\odot}$ in cosmological
simulations (Illustris (Genel et al., 2014; Vogelsberger et al., 2014),
Illustris TNG (Weinberger et al., 2017; Pillepich et al., 2018), Mufasa (Davé
et al., 2016), Simba (Davé et al., 2019), EAGLE (Crain et al., 2015; Schaller
et al., 2015; Schaye et al., 2015)), zoom-in simulations (FIRE-2 (Hopkins et
al., 2018), g14 (Brooks & Zolotov, 2014; Brooks & Christensen, 2016;
Christensen et al., 2016; Christensen et al., 2018; Brooks et al., 2017),
Marvel/Justice League) (Bellovary et al., 2019), semi-analytical models
(Somerville et al., 2008; Porter et al., 2014; Somerville et al., 2015;
Brennan et al., 2017) and semi-empirical UniverseMachine model (Behroozi et
al., 2013, 2019).

Iyer et al. (2020) found that different models exhibit a wide range of SFR
stochasticity. The largest stochasticity levels are reached in the simulations
with efficient, explosive feedback that do not rely on subgrid ISM models,
such as EAGLE and FIRE-2 simulations. Simulations that do use a subgrid ISM
model exhibit considerably less stochasticity on time scales $\lesssim 100$
Myr. Semi-analytic models that do not model local ISM structure and do not add
stochasticity explicitly, such as the Santa Cruz SAM and Universe Machine,
exhibit lower SFR stochasticity still. This is because the SFR stochasticity
arising due to variations in the mass assembly histories and internal global
parameters of galaxies are sub-dominant to the stochasticity on $\lesssim 100$
Myr time scales induced by the rise and fall of individual star-forming
regions in galaxies (Tacchella et al., 2020; Shin et al., 2023).

The inherent level of stochasticity in the GRUMPY model is comparable to that
of simulations with subrid physics and semi-analytic models. On the other
hand, the addition of explicit SFR stochasticity at the level required to
match observed scatter in the $M_{\star}$-SFR relation brings the SFR
stochasticity to the level similar to that exhibited by low-mass galaxies in
the FIRE-2 and EAGLE simulations.

###  4.2 Scatter of the $g-r$ colours compared to observations and simulations

The main challenge that we identified is reproducing a large observed scatter
of $g-r$ colours of faint dwarf galaxies ($M_{r}\gtrsim-15$). Chaves-Montero &
Hearin (2021) examined effects of star formation rate burstiness on the
colours of massive galaxies ($\log_{10}M_{\star}>9.25$) in the Illustris TNG
simulations and model galaxies in the UniverseMachine. They found that
burstiness can affect optical colours of individual galaxies at the $\approx
0.2-0.3$ level, but does not significantly change the colour distribution of
galaxy population as a whole. This is different from our finding in the dwarf
galaxy regime, likely because the level of burstiness implied by the observed
scatter in the SFR-$M_{\star}$ regime is significantly larger than the level
considered in that study.

Comparison of observations with the Auriga simulation suite of MW-sized
galaxies (Grand et al., 2017) revealed that the colour-magnitude distribution
for $M_{\star}\lesssim 10^{6}\,M_{\odot}$ in simulations is considerably
narrower than in observations, and there are notably fewer faint blue dwarf
galaxies in simulations than in observations (Pan et al., 2023). These authors
attributed this discrepancy to the absence of low-temperature/molecular gas
cooling in the stellar feedback and ISM model used in the Auriga simulations.
In the case of faint dwarf galaxies, the molecular gas self-shielding effect
is likely underestimated when the UV radiation perturbs the gas due to the
limitations of the model. Therefore, most smallest dwarf galaxies are quenched
and become red, with almost no star-forming blue galaxies forming following
reionization.

This is similar to what we find in the GRUMPY model for smaller galaxies. The
reasons may also be similar as star formation in this model is based on the
model of molecular hydrogen of Gnedin & Draine (2014), which was not tested at
the metallicities of the ultra-faint galaxies. At the same time, star
formation in low-metallicity galaxies may not be well traced by molecular
hydrogen to begin with (Nobels et al., 2023, Polzin et al. in prep.). Thus,
the star formation rate may be underestimated in the faintest dwarf galaxies
due to the adopted star formation model in GRUMPY.

Field galaxies of $M_{\star}\gtrsim 10^{6}\,M_{\odot}$ in both Auriga and
APOSTLE simulations continue forming stars to low redshifts, although
satellite galaxies do quench after their accretion onto MW-sized system (Digby
et al., 2019). Likewise, Garrison-Kimmel et al. (2019) also find that galaxies
with $M_{\star}\lesssim 10^{6}\,M_{\star}$ are quenched by UV heating after
reionization in the FIRE-2 zoom-in simulations of dwarf galaxies. However,
they also find that galaxies of larger mass generally continue forming stars
until $z=0$.

In the GRUMPY model galaxies of larger mass also continue to form stars and
their colours are affected by SFR stochasticity as was demonstrated above.

###  4.3 Contribution of uncertainties in observational SFR indicators to the
apparent SFR stochasticity

Although there are physical drivers of bursty star formation in dwarf
galaxies, part of the apparent scatter of the SFR in these galaxies could be
due to uncertainty in the observational SFR indicators. Different observed
bands contain distinct information about star formation, with shorter
wavelengths probing most recent star formation, while longer wavelengths
probing star formation over longer time scales (e.g., Boquien et al., 2014).
Additionally, spectral features such as $H\alpha$ and UV-based SFR indicators,
$H\alpha$ and Ca-H,K absorption line equivalent widths, the Dn(4000) spectral
break, and broadband galaxy colours are sensitive to the model parameters in
the PSD formalism (Iyer et al., 2022). Therefore, comprehending the
variability of SFRs across different indicators can provide information about
stochasticity on different time scales and inform comparisons with model
results.

The most commonly used SFR indicators are $H\alpha$ and UV+IR fluxes
(Kennicutt, 1998; Kennicutt & Evans, 2012; Boquien et al., 2014; Davies et
al., 2016). $H\alpha$ photons originate from gas ionized by young stars ($<$20
Myr), providing a measure of the current SFR $(\sim 4-10)$ Myr in galaxies. In
dwarf galaxies, however, Hα can significantly underpredict the star formation
rate due to paucity of massive stars (e.g., Lee et al., 2009). UV emission
directly comes from star-forming regions, while some is absorbed and
reprocessed by dust, re-emerging in the IR. Thus, the sum of UV and total IR
luminosities provides the total SFR over the timescale of $\approx 30-100$
Myr.

The distribution of SFRs measured using two indicators can be used to
investigate the change in SFRs over time and provide a statistical view of the
behavior of the galaxy population. However, such analysis is challenging due
to several uncertainties and assumptions, including the stochastic IMF and
dust properties, the monotonicity of SFRs over different timescales, and the
inherent difficulty in measuring SFRs accurately (Johnson et al., 2013;
Shivaei et al., 2018).

At low SFRs, inferring an SFR from a star formation indicator (SFI) in an
individual galaxy is problematic. When calculating the sum of the mass and
ages of all the stars that comprise a population using stellar population
synthesis (SPS) models to map observed luminosity to underlying physical
properties, it is common to assume that the IMF and SFH are well populated.
However, if the IMF is not fully sampled, the luminosity’s highly non-linear
dependence on stellar mass can have significant consequences for the
population’s luminosity. Thus, there is no longer a deterministic relation
between the total mass and age of the population to the total luminosity and
colour of its integrated light. This makes the inverse problem of determining
the mass or age of a simple stellar population from its photometric properties
ill-posed. When the IMF and SFH are sparsely sampled in the low SFR regime,
there is no unique mapping between SFRs and SFIs, and thus no unique way to
infer an SFR from an SFI in an individual (da Silva et al., 2014).

For SFRs $\lesssim 10^{-4}M_{\odot}$/yr, the scatter in the $H\alpha$
indicator, which is used for most dwarfs at these masses, is expected to be
very high due to insufficient sampling of IMF and SFH. In addition, Lee et al.
(2009) demonstrated that for lower luminosity dwarf galaxies (roughly less
active than the Small Magellanic Cloud), H$\alpha$ tends to increasingly
underpredict the total SFR relative to the FUV. By SFR $\sim 0.003M_{\odot}$
yr-1, the average H$\alpha$-to-FUV flux ratio is lower than expected by a
factor of two, and at the lowest SFRs probed, the ratio exhibits an order of
magnitude discrepancy for the handful of galaxies that remain in the sample
due to the IMF deficiency in dwarf and low surface brightness galaxies. The
upper panels in Figure 4 of da Silva et al. (2014) illustrate these results.
According to Cignoni et al. (2019), the level of stochasticity on a 10
million-year time scale is quite significant, whereas it becomes much less
significant on a 100 million-year time scale.

###  4.4 Caveats

The interpretation of our results should be approached with caution due to the
constraints of our simple regulator-type model.

First, we make the assumption of a smooth, radial-averaged exponential profile
for both the gas surface density $\Sigma_{g}$ and the molecular gas surface
density $\Sigma_{\rm{HI}}$. However, in reality, the interstellar medium (ISM)
can be highly structured (Kravtsov, 2003; Semenov et al., 2017), and a radial
average may not accurately capture the overdensities and underdensities within
the ISM. This is particularly relevant for low-mass galaxies, where star
formation occurs in a small number of bright star-forming regions, making it
difficult to effectively average out fluctuations. In blue compact dwarf
galaxies, some giant molecular clouds are massive enough to produce
significant torques on surrounding disk material, leading to kinematic
decoupling and stirring up the inhomogeneity of the ISM. Thus, a closer
examination of the detailed structure of the ISM is necessary in future
studies.

Second, we applied the same level of stochasticity to all dwarf galaxies with
stellar masses between $10^{5}M_{\odot}$ and $10^{9}M_{\odot}$. However,
previous studies have demonstrated that the intrinsic scatter of the main
sequence (MS) for star-forming galaxies with stellar masses between
$10^{8}M_{\odot}$ and $10^{11}M_{\odot}$ does not remain constant, but
exhibits a "U-shape" trend (Davies et al., 2022; Huang et al., 2023).
Specifically, for low-mass galaxies with
$10^{8}M_{\odot}<M_{\star}<10^{10}M_{\odot}$, stochastic starbursts and
stellar feedback events are the main drivers of intrinsic MS scatter. For
intermediate-mass galaxies around $10^{10}M_{\odot}$, the effect of star
formation and feedback becomes less pronounced, resulting in a decrease in
intrinsic MS scatter. In contrast, high-mass galaxies
($M_{\star}>10^{10.3}M_{\odot}$) exhibit an increased intrinsic MS scatter due
to AGN feedback. For even lower-mass dwarf galaxies with
$M_{\star}<10^{8}M_{\odot}$, the star formation rate is subject to small-
number statistics and becomes more bursty. Therefore, to accurately model
these variations, the level of stochasticity incorporated in the model should
be a function of stellar mass, with higher stochasticity for low-mass galaxies
and decreasing stochasticity as we move towards intermediate masses around
$10^{10}M_{\odot}$. However, as this project is a proof of concept, we defer a
more detailed study of the dependence of stochasticity on galaxy stellar mass
to future work, focusing instead on demonstrating the non-negligible impact of
incorporating stochasticity on several galaxy scaling relations.

##  5 Conclusions

We present the impact of bursty star formation on several galaxy scaling
relations in dwarf galaxies using the GRUMPY model. This simple regulator-type
model reproduces the average trend of the SFR–mass, stellar mass–gas mass,
stellar mass–metallicity relation, and CMD. However, the scatter of these
relations in the original model is smaller than observed. To address this, we
utilize a formula for the power spectral density (PSD) to generate a series of
random numbers with a log-normal distribution to perturb the original SFR. We
adjust the parameters in the PSD formalism to control the amplitude and
variability of burstiness in the model. We experiment with the effects of
long-term and short-term SFR variations in the model and compare them to
observations. The main results can be summarized as follows:

  1. 1.

We find that the PSD of equation 12 with parameters $\alpha=2$,
$\tau_{\rm{break}}=2$Gyr, and $\sigma=0.6$ produces both short-term and long-
term variations, and the amplitude of burstiness comparable to zoom-in galaxy
formation simulations, such as FIRE-2.

  2. 2.

Bursty star formation increases the scatter in the SFR-$M_{\star}$ relation
compared to the non-bursty case. We find, however, that rare extreme starburst
galaxies, such as those in Lin et al. (2022) imply even larger levels of
stochasticity.

  3. 3.

Adding the SFR stochasticity increases the scatter in the CMD (Figure 5) in
brighter dwarf galaxies $(M_{V}<-12)$. For fainter dwarf galaxies
$(M_{V}>-12)$, the scatter is comparable to that of the non-bursty model. We
confirmed that this result is not due to low-mass dwarf galaxies being
quenched by reionization.

  4. 4.

We show that metallicity scatter can increase the scatter of colours of faint
dwarfs $(M_{V}>-12)$, however this leads to the overestimate of scatter in the
mass-metallicity relation compared to observed Local Group dwarf satellite
galaxies.

## Acknowledgements

This research made use of data from the SAGA Survey (saga- survey.org), which
was supported by NSF collaborative grants AST-1517148 and AST-1517422 and by
Heising-Simons Foundation grant 2019-1402.

The software packages NUMPY (van der Walt et al., 2011), SCIPY (Jones et al.,
2001), Astropy (Astropy Collaboration et al., 2018), MATPLOTLIB (Hunter,
2007), and GITHUB were invaluable in conducting the analyses presented in this
paper. Additionally, this research relied heavily on the Astrophysics Data
Service (ADS) and the preprint repository arXiv.

We thank Anirudh Chiti and Nick Gnedin for constructive feedback on an earlier
version of this paper.

## Data Availability

The data underlying this article will be shared on reasonable request to the
authors.

## References

  * Anglés-Alcázar et al. (2017) Anglés-Alcázar D., Faucher-Giguère C.-A., Quataert E., Hopkins P. F., Feldmann R., Torrey P., Wetzel A., Kereš D., 2017, MNRAS, 472, L109
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Atek et al. (2023) Atek H., et al., 2023, arXiv e-prints, p. arXiv:2308.08540
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A., Primack J. R., 2013, ApJ, 763, 18
  * Behroozi et al. (2019) Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143
  * Bellovary et al. (2019) Bellovary J. M., Cleary C. E., Munshi F., Tremmel M., Christensen C. R., Brooks A., Quinn T. R., 2019, MNRAS, 482, 2913
  * Berg et al. (2012) Berg D. A., et al., 2012, ApJ, 754, 98
  * Berg et al. (2016) Berg D. A., Skillman E. D., Henry R. B. C., Erb D. K., Carigi L., 2016, ApJ, 827, 126
  * Bigiel & Blitz (2012) Bigiel F., Blitz L., 2012, ApJ, 756, 183
  * Bigiel et al. (2008) Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2846
  * Birrer et al. (2014) Birrer S., Lilly S., Amara A., Paranjape A., Refregier A., 2014, ApJ, 793, 12
  * Bolatto et al. (2011) Bolatto A. D., et al., 2011, ApJ, 741, 12
  * Boquien et al. (2014) Boquien M., Buat V., Perret V., 2014, A&A, 571, A72
  * Bouché et al. (2010) Bouché N., et al., 2010, ApJ, 718, 1001
  * Brennan et al. (2017) Brennan R., et al., 2017, MNRAS, 465, 619
  * Brinchmann et al. (2004) Brinchmann J., Charlot S., White S. D. M., Tremonti C., Kauffmann G., Heckman T., Brinkmann J., 2004, MNRAS, 351, 1151
  * Brooks & Christensen (2016) Brooks A., Christensen C., 2016, in Laurikainen E., Peletier R., Gadotti D., eds, Astrophysics and Space Science Library Vol. 418, Galactic Bulges. p. 317 (arXiv:1511.04095), doi:10.1007/978-3-319-19378-6_12
  * Brooks & Zolotov (2014) Brooks A. M., Zolotov A., 2014, ApJ, 786, 87
  * Brooks et al. (2017) Brooks A. M., Papastergis E., Christensen C. R., Governato F., Stilp A., Quinn T. R., Wadsley J., 2017, ApJ, 850, 97
  * Cano-Díaz et al. (2016) Cano-Díaz M., et al., 2016, ApJ, 821, L26
  * Caplar & Tacchella (2019) Caplar N., Tacchella S., 2019, MNRAS, 487, 3845
  * Carlsten et al. (2020) Carlsten S. G., Greco J. P., Beaton R. L., Greene J. E., 2020, ApJ, 891, 144
  * Carlsten et al. (2022) Carlsten S. G., Greene J. E., Beaton R. L., Danieli S., Greco J. P., 2022, ApJ, 933, 47
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chang et al. (2015) Chang Y.-Y., van der Wel A., da Cunha E., Rix H.-W., 2015, ApJS, 219, 8
  * Chaves-Montero & Hearin (2021) Chaves-Montero J., Hearin A., 2021, MNRAS, 506, 2373
  * Christensen et al. (2016) Christensen C. R., Davé R., Governato F., Pontzen A., Brooks A., Munshi F., Quinn T., Wadsley J., 2016, ApJ, 824, 57
  * Christensen et al. (2018) Christensen C. R., Davé R., Brooks A., Quinn T., Shen S., 2018, ApJ, 867, 142
  * Ciesla et al. (2014) Ciesla L., et al., 2014, A&A, 565, A128
  * Ciesla et al. (2023) Ciesla L., et al., 2023, arXiv e-prints, p. arXiv:2309.15720
  * Cignoni et al. (2019) Cignoni M., et al., 2019, ApJ, 887, 112
  * Conroy & Gunn (2010) Conroy C., Gunn J. E., 2010, ApJ, 712, 833
  * Conroy et al. (2009) Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * D’Silva et al. (2023) D’Silva J. C. J., Lagos C. D. P., Davies L. J. M., Lovell C. C., Vijayan A. P., 2023, MNRAS, 518, 456
  * Daddi et al. (2007) Daddi E., et al., 2007, ApJ, 670, 156
  * Daddi et al. (2010) Daddi E., et al., 2010, ApJ, 713, 686
  * Davé et al. (2011) Davé R., Oppenheimer B. D., Finlator K., 2011, MNRAS, 415, 11
  * Davé et al. (2016) Davé R., Thompson R., Hopkins P. F., 2016, MNRAS, 462, 3265
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * Davies et al. (2016) Davies L. J. M., et al., 2016, MNRAS, 461, 458
  * Davies et al. (2019) Davies L. J. M., et al., 2019, MNRAS, 483, 1881
  * Davies et al. (2022) Davies L. J. M., et al., 2022, MNRAS, 509, 4392
  * De Lucia et al. (2020) De Lucia G., Xie L., Fontanot F., Hirschmann M., 2020, MNRAS, 498, 3215
  * Dekel & Mandelker (2014) Dekel A., Mandelker N., 2014, MNRAS, 444, 2071
  * Digby et al. (2019) Digby R., et al., 2019, MNRAS, 485, 5423
  * Dressler et al. (2023) Dressler A., et al., 2023, arXiv e-prints, p. arXiv:2306.02469
  * Durbala et al. (2020) Durbala A., Finn R. A., Crone Odekon M., Haynes M. P., Koopmann R. A., O’Donoghue A. A., 2020, AJ, 160, 271
  * Emami et al. (2019) Emami N., Siana B., Weisz D. R., Johnson B. D., Ma X., El-Badry K., 2019, ApJ, 881, 71
  * Endsley et al. (2023) Endsley R., et al., 2023, arXiv e-prints, p. arXiv:2306.05295
  * Fall et al. (2010) Fall S. M., Krumholz M. R., Matzner C. D., 2010, ApJ, 710, L142
  * Faucher-Giguère (2018) Faucher-Giguère C.-A., 2018, MNRAS, 473, 3717
  * Feldmann (2013) Feldmann R., 2013, MNRAS, 433, 1910
  * Finlator & Davé (2008) Finlator K., Davé R., 2008, MNRAS, 385, 2181
  * Fumagalli et al. (2011) Fumagalli M., da Silva R. L., Krumholz M. R., 2011, ApJ, 741, L26
  * Garnett (2002) Garnett D. R., 2002, ApJ, 581, 1019
  * Garrison-Kimmel et al. (2019) Garrison-Kimmel S., et al., 2019, MNRAS, 489, 4574
  * Geha et al. (2017) Geha M., et al., 2017, ApJ, 847, 4
  * Genel et al. (2014) Genel S., et al., 2014, MNRAS, 445, 175
  * Genzel et al. (2010) Genzel R., et al., 2010, MNRAS, 407, 2091
  * Gnedin & Draine (2014) Gnedin N. Y., Draine B. T., 2014, ApJ, 795, 37
  * Gnedin et al. (2008) Gnedin N. Y., Kravtsov A. V., Chen H.-W., 2008, ApJ, 672, 765
  * Grand et al. (2017) Grand R. J. J., et al., 2017, MNRAS, 467, 179
  * Grand et al. (2021) Grand R. J. J., et al., 2021, MNRAS, 507, 4953
  * Greene et al. (2023) Greene J. E., Danieli S., Carlsten S., Beaton R., Jiang F., Li J., 2023, ApJ, 949, 94
  * Griffen et al. (2016) Griffen B. F., Ji A. P., Dooley G. A., Gómez F. A., Vogelsberger M., O’Shea B. W., Frebel A., 2016, ApJ, 818, 10
  * Guo et al. (2013) Guo K., Zheng X. Z., Fu H., 2013, ApJ, 778, 23
  * Hopkins et al. (2014) Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Huang et al. (2023) Huang R., Battisti A. J., Grasha K., da Cunha E., Lagos C. d. P., Leslie S. K., Wisnioski E., 2023, MNRAS, 520, 446
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Iyer et al. (2020) Iyer K. G., et al., 2020, MNRAS, 498, 430
  * Iyer et al. (2022) Iyer K. G., Speagle J. S., Caplar N., Forbes J. C., Gawiser E., Leja J., Tacchella S., 2022, arXiv e-prints, p. arXiv:2208.05938
  * James et al. (2017) James B. L., Koposov S. E., Stark D. P., Belokurov V., Pettini M., Olszewski E. W., McQuinn K. B. W., 2017, MNRAS, 465, 3977
  * Jimmy et al. (2015) Jimmy Tran K.-V., Saintonge A., Accurso G., Brough S., Oliva-Altamirano P., 2015, ApJ, 812, 98
  * Johnson et al. (2013) Johnson B. D., et al., 2013, ApJ, 772, 8
  * Jones et al. (2001) Jones E., Oliphant T., Peterson P., 2001, SciPy: Open Source Scientific Tools for Python, http://www.scipy.org
  * Jones et al. (2023) Jones M. G., et al., 2023, arXiv e-prints, p. arXiv:2310.01478
  * Katsianis et al. (2019) Katsianis A., et al., 2019, ApJ, 879, 11
  * Katz et al. (2020) Katz H., et al., 2020, MNRAS, 494, 2200
  * Kauffmann et al. (2003) Kauffmann G., et al., 2003, MNRAS, 341, 54
  * Kennicutt (1989) Kennicutt Robert C. J., 1989, ApJ, 344, 685
  * Kennicutt (1998) Kennicutt Robert C. J., 1998, ARA&A, 36, 189
  * Kennicutt & Evans (2012) Kennicutt R. C., Evans N. J., 2012, Annual Review of Astronomy and Astrophysics, 50, 531
  * Kirby et al. (2013) Kirby E. N., Cohen J. G., Guhathakurta P., Cheng L., Bullock J. S., Gallazzi A., 2013, ApJ, 779, 102
  * Kojima et al. (2020) Kojima T., et al., 2020, ApJ, 898, 142
  * Kormendy et al. (2009) Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009, ApJS, 182, 216
  * Kravtsov (2003) Kravtsov A. V., 2003, ApJ, 590, L1
  * Kravtsov & Manwadkar (2022) Kravtsov A., Manwadkar V., 2022, MNRAS, 514, 2667
  * Kruijssen & Longmore (2014) Kruijssen J. M. D., Longmore S. N., 2014, MNRAS, 439, 3239
  * Kruijssen et al. (2019) Kruijssen J. M. D., et al., 2019, Nature, 569, 519
  * Krumholz & Dekel (2012) Krumholz M. R., Dekel A., 2012, ApJ, 753, 16
  * Krumholz & Kruijssen (2015) Krumholz M. R., Kruijssen J. M. D., 2015, MNRAS, 453, 739
  * Kurczynski et al. (2016) Kurczynski P., et al., 2016, ApJ, 820, L1
  * Lee et al. (2006) Lee H., Skillman E. D., Cannon J. M., Jackson D. C., Gehrz R. D., Polomski E. F., Woodward C. E., 2006, ApJ, 647, 970
  * Lee et al. (2009) Lee J. C., et al., 2009, ApJ, 706, 599
  * Lee et al. (2015) Lee N., et al., 2015, ApJ, 801, 80
  * Leitner (2012) Leitner S. N., 2012, ApJ, 745, 149
  * Leitner & Kravtsov (2011) Leitner S. N., Kravtsov A. V., 2011, ApJ, 734, 48
  * Leroy et al. (2013) Leroy A. K., et al., 2013, AJ, 146, 19
  * Lilly et al. (2013) Lilly S. J., Carollo C. M., Pipino A., Renzini A., Peng Y., 2013, ApJ, 772, 119
  * Lin et al. (2022) Lin Y.-H., et al., 2022, arXiv e-prints, p. arXiv:2211.02094
  * Looser et al. (2023) Looser T. J., et al., 2023, arXiv e-prints, p. arXiv:2306.02470
  * Mancuso et al. (2016) Mancuso C., Lapi A., Shi J., Cai Z. Y., Gonzalez-Nuevo J., Béthermin M., Danese L., 2016, ApJ, 833, 152
  * Manwadkar & Kravtsov (2022) Manwadkar V., Kravtsov A. V., 2022, MNRAS, 516, 3944
  * Mao et al. (2021) Mao Y.-Y., Geha M., Wechsler R. H., Weiner B., Tollerud E. J., Nadler E. O., Kallivayalil N., 2021, ApJ, 907, 85
  * Martín-Manjón et al. (2008) Martín-Manjón M. L., Mollá M., Díaz A. I., Terlevich R., 2008, MNRAS, 385, 854
  * Mason et al. (2023) Mason C. A., Trenti M., Treu T., 2023, MNRAS, 521, 497
  * Matthee & Schaye (2019) Matthee J., Schaye J., 2019, MNRAS, 484, 915
  * McGaugh et al. (2017) McGaugh S. S., Schombert J. M., Lelli F., 2017, ApJ, 851, 22
  * McQuinn et al. (2015a) McQuinn K. B. W., Mitchell N. P., Skillman E. D., 2015a, ApJS, 218, 29
  * McQuinn et al. (2015b) McQuinn K. B. W., et al., 2015b, ApJ, 815, L17
  * McQuinn et al. (2019) McQuinn K. B. W., van Zee L., Skillman E. D., 2019, ApJ, 886, 74
  * McQuinn et al. (2020) McQuinn K. B. W., et al., 2020, ApJ, 891, 181
  * McQuinn et al. (2021) McQuinn K. B. W., et al., 2021, ApJ, 918, 23
  * McQuinn et al. (2023) McQuinn K. B. W., Mao Y.-Y., Buckley M. R., Shih D., Cohen R. E., Dolphin A. E., 2023, ApJ, 944, 14
  * Mo et al. (1998) Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319
  * Moster et al. (2018) Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822
  * Muñoz et al. (2018) Muñoz R. R., Côté P., Santana F. A., Geha M., Simon J. D., Oyarzún G. A., Stetson P. B., Djorgovski S. G., 2018, ApJ, 860, 66
  * Nobels et al. (2023) Nobels F. S. J., Schaye J., Schaller M., Ploeckinger S., Chaikin E., Richings A. J., 2023, arXiv e-prints, p. arXiv:2309.13750
  * Noeske et al. (2007) Noeske K. G., et al., 2007, ApJ, 660, L43
  * Oppenheimer & Davé (2008) Oppenheimer B. D., Davé R., 2008, MNRAS, 387, 577
  * Orr et al. (2019) Orr M. E., Hayward C. C., Hopkins P. F., 2019, MNRAS, 486, 4724
  * Pan et al. (2023) Pan Y., et al., 2023, MNRAS, 519, 4499
  * Pannella et al. (2009) Pannella M., et al., 2009, ApJ, 698, L116
  * Parrish et al. (2009) Parrish I. J., Quataert E., Sharma P., 2009, ApJ, 703, 96
  * Peeples & Shankar (2011) Peeples M. S., Shankar F., 2011, MNRAS, 417, 2962
  * Pessa et al. (2021) Pessa I., et al., 2021, A&A, 650, A134
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Popesso et al. (2023) Popesso P., et al., 2023, MNRAS, 519, 1526
  * Porter et al. (2014) Porter L. A., Somerville R. S., Primack J. R., Johansson P. H., 2014, MNRAS, 444, 942
  * Rahman et al. (2012) Rahman N., et al., 2012, ApJ, 745, 183
  * Renzini & Peng (2015) Renzini A., Peng Y.-j., 2015, ApJ, 801, L29
  * Rey et al. (2019) Rey M. P., Pontzen A., Agertz O., Orkney M. D. A., Read J. I., Saintonge A., Pedersen C., 2019, ApJ, 886, L3
  * Rey et al. (2020) Rey M. P., Pontzen A., Agertz O., Orkney M. D. A., Read J. I., Rosdahl J., 2020, MNRAS, 497, 1508
  * Rodighiero et al. (2010) Rodighiero G., et al., 2010, A&A, 518, L25
  * Rodríguez-Puebla et al. (2016) Rodríguez-Puebla A., Primack J. R., Behroozi P., Faber S. M., 2016, MNRAS, 455, 2592
  * Salim et al. (2007) Salim S., et al., 2007, ApJS, 173, 267
  * Santini et al. (2017) Santini P., et al., 2017, ApJ, 847, 76
  * Scalo & Struck-Marcell (1984) Scalo J. M., Struck-Marcell C., 1984, ApJ, 276, 60
  * Scalo & Struck-Marcell (1986) Scalo J. M., Struck-Marcell C., 1986, ApJ, 301, 77
  * Schaller et al. (2015) Schaller M., Dalla Vecchia C., Schaye J., Bower R. G., Theuns T., Crain R. A., Furlong M., McCarthy I. G., 2015, MNRAS, 454, 2277
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schreiber et al. (2015) Schreiber C., et al., 2015, A&A, 575, A74
  * Semenov et al. (2017) Semenov V. A., Kravtsov A. V., Gnedin N. Y., 2017, ApJ, 845, 133
  * Semenov et al. (2021a) Semenov V. A., Kravtsov A. V., Caprioli D., 2021a, ApJ, 910, 126
  * Semenov et al. (2021b) Semenov V. A., Kravtsov A. V., Gnedin N. Y., 2021b, ApJ, 918, 13
  * Shin et al. (2023) Shin E.-j., Tacchella S., Kim J.-h., Iyer K. G., Semenov V. A., 2023, ApJ, 947, 61
  * Shivaei et al. (2018) Shivaei I., et al., 2018, ApJ, 855, 42
  * Simon (2019) Simon J. D., 2019, ARA&A, 57, 375
  * Somerville et al. (2008) Somerville R. S., Hopkins P. F., Cox T. J., Robertson B. E., Hernquist L., 2008, MNRAS, 391, 481
  * Somerville et al. (2015) Somerville R. S., Popping G., Trager S. C., 2015, MNRAS, 453, 4337
  * Sparre et al. (2015) Sparre M., et al., 2015, MNRAS, 447, 3548
  * Sparre et al. (2017) Sparre M., Hayward C. C., Feldmann R., Faucher-Giguère C.-A., Muratov A. L., Kereš D., Hopkins P. F., 2017, MNRAS, 466, 88
  * Speagle et al. (2014) Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15
  * Tacchella et al. (2016) Tacchella S., Dekel A., Carollo C. M., Ceverino D., DeGraf C., Lapiner S., Mandelker N., Primack Joel R., 2016, MNRAS, 457, 2790
  * Tacchella et al. (2018) Tacchella S., Bose S., Conroy C., Eisenstein D. J., Johnson B. D., 2018, ApJ, 868, 92
  * Tacchella et al. (2020) Tacchella S., Forbes J. C., Caplar N., 2020, MNRAS, 497, 698
  * Teich et al. (2016) Teich Y. G., et al., 2016, ApJ, 832, 85
  * Torrey et al. (2014) Torrey P., Vogelsberger M., Genel S., Sijacki D., Springel V., Hernquist L., 2014, MNRAS, 438, 1985
  * Torrey et al. (2018) Torrey P., et al., 2018, MNRAS, 477, L16
  * Torrey et al. (2019) Torrey P., et al., 2019, MNRAS, 484, 5587
  * Vincenzo et al. (2016) Vincenzo F., Matteucci F., Belfiore F., Maiolino R., 2016, MNRAS, 455, 4183
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * Wang & Lilly (2020a) Wang E., Lilly S. J., 2020a, ApJ, 892, 87
  * Wang & Lilly (2020b) Wang E., Lilly S. J., 2020b, ApJ, 895, 25
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Whitaker et al. (2012) Whitaker K. E., van Dokkum P. G., Brammer G., Franx M., 2012, ApJ, 754, L29
  * Wuyts et al. (2011) Wuyts S., et al., 2011, ApJ, 742, 96
  * Zick et al. (2018) Zick T. O., Weisz D. R., Boylan-Kolchin M., 2018, MNRAS, 477, 480
  * da Silva et al. (2012) da Silva R. L., Fumagalli M., Krumholz M., 2012, ApJ, 745, 145
  * da Silva et al. (2014) da Silva R. L., Fumagalli M., Krumholz M. R., 2014, MNRAS, 444, 3275
  * van Loon et al. (2021) van Loon M. L., Mitchell P. D., Schaye J., 2021, MNRAS, 504, 4817
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science and Engineering, 13, 22