Stochastic modeling of star-formation histories I: the scatter of the   star-forming main sequence

Neven Caplar; Sandro Tacchella

arXiv:1901.07556·astro-ph.GA·June 5, 2019

Stochastic modeling of star-formation histories I: the scatter of the star-forming main sequence

Neven Caplar, Sandro Tacchella

PDF

TL;DR

This paper introduces a stochastic framework to model galaxy star-formation histories, linking the scatter in the star-forming main sequence to underlying temporal correlations and decorrelation timescales.

Contribution

It develops a novel stochastic modeling approach using a broken power-law power spectrum to interpret the scatter in galaxy star-formation rates.

Findings

01

Estimated decorrelation timescale of ~200 Myr for star-formation histories.

02

Constraints on the power spectrum parameters from multi-band observational data.

03

Highlighting the importance of baryonic effects over dynamical timescales.

Abstract

We present a framework for modelling the star-formation histories of galaxies as a stochastic process. We define this stochastic process through a power spectrum density with a functional form of a broken power-law. Star-formation histories are correlated on short timescales, the strength of this correlation described by a power-law slope, $α$ , and they decorrelate to resemble white noise over a timescale that is proportional to the timescale of the break in the power spectrum density, $τ_{break}$ . We use this framework to explore the properties of the stochastic process that, we assume, gives rise to the log-normal scatter about the relationship between star-formation rate and stellar mass, the so-called galaxy star-forming main sequence. Specifically, we show how the measurements of the normalisation and width ( $σ_{MS}$ ) of the main sequence, measured in several…

Equations72

A C F (δ t) = \frac{lim _{T - > \infty} \frac{1}{T} \int _{(T)} Δ ( t ) Δ ( t + δ t ) d t}{σ _{Δ}^{2}}

A C F (δ t) = \frac{lim _{T - > \infty} \frac{1}{T} \int _{(T)} Δ ( t ) Δ ( t + δ t ) d t}{σ _{Δ}^{2}}

\frac{d Δ ( t )}{d t} + \frac{1}{τ _{decor}} Δ (t) = ϵ (t),

\frac{d Δ ( t )}{d t} + \frac{1}{τ _{decor}} Δ (t) = ϵ (t),

P S D (f) = \frac{σ _{int}^{2}}{1/ τ _{decor}^{2} + ( 2 π f ) ^{2}} = \frac{σ ^{2}}{1 + ( τ _{break} f ) ^{2}} .

P S D (f) = \frac{σ _{int}^{2}}{1/ τ _{decor}^{2} + ( 2 π f ) ^{2}} = \frac{σ ^{2}}{1 + ( τ _{break} f ) ^{2}} .

A C F (t) = exp (- t / τ_{decor}) .

A C F (t) = exp (- t / τ_{decor}) .

P S D (f) = \frac{σ ^{2}}{1 + ( τ _{break} f ) ^{α}} .

P S D (f) = \frac{σ ^{2}}{1 + ( τ _{break} f ) ^{α}} .

B = \frac{σ / μ - 1}{σ / μ + 1},

B = \frac{σ / μ - 1}{σ / μ + 1},

SFR_{ind} (t_{i}) = \frac{1}{N _{j}} t_{j} = 1 \sum N_{j} SFR (t_{i} - t_{j}) .

SFR_{ind} (t_{i}) = \frac{1}{N _{j}} t_{j} = 1 \sum N_{j} SFR (t_{i} - t_{j}) .

Δ (t_{i}) = lo g SFR (t_{i}) - lo g SFR_{MS} (M_{⋆}, t_{i}),

Δ (t_{i}) = lo g SFR (t_{i}) - lo g SFR_{MS} (M_{⋆}, t_{i}),

Δ_{ind} (t_{i}) = lo g SFR_{ind} (t_{i}) - lo g SFR_{MS} (M_{⋆}, t_{i}) .

Δ_{ind} (t_{i}) = lo g SFR_{ind} (t_{i}) - lo g SFR_{MS} (M_{⋆}, t_{i}) .

Δ_{ind} (t_{i}) = lo g \frac{1}{N _{j}} t_{j} = 1 \sum N_{j} SFR (t_{i} - t_{j}) - lo g SFR_{MS} = lo g \frac{1}{N _{j}} t_{j} = 1 \sum N_{j} 1 0^{Δ (t_{i} - t_{j}) + l o g SFR_{MS}} - lo g SFR_{MS} = lo g \frac{1}{N _{j}} t_{j} = 1 \sum N_{j} 1 0^{Δ (t_{i} - t_{j})} .

Δ_{ind} (t_{i}) = lo g \frac{1}{N _{j}} t_{j} = 1 \sum N_{j} SFR (t_{i} - t_{j}) - lo g SFR_{MS} = lo g \frac{1}{N _{j}} t_{j} = 1 \sum N_{j} 1 0^{Δ (t_{i} - t_{j}) + l o g SFR_{MS}} - lo g SFR_{MS} = lo g \frac{1}{N _{j}} t_{j} = 1 \sum N_{j} 1 0^{Δ (t_{i} - t_{j})} .

P [Δ (t_{i} - t_{j})] = N (Δ (t_{j}) exp (- \frac{t _{i} - t _{j}}{4 τ _{decor}}), σ_{MS} 1 - exp (- \frac{t _{i} - t _{j}}{2 τ _{decor}})),

P [Δ (t_{i} - t_{j})] = N (Δ (t_{j}) exp (- \frac{t _{i} - t _{j}}{4 τ _{decor}}), σ_{MS} 1 - exp (- \frac{t _{i} - t _{j}}{2 τ _{decor}})),

E [1 0^{N (μ, s)}] = 1 0^{μ + \frac{s ^{2} l n ( 10 )}{2}} .

E [1 0^{N (μ, s)}] = 1 0^{μ + \frac{s ^{2} l n ( 10 )}{2}} .

E [1 0^{Δ (t_{i} - t_{j})}] = 1 0^{Δ (t_{j}) e x p (- \frac{t _{i} - t _{j}}{4 τ _{decor}}) + σ_{MS}^{2} (1 - e x p (- \frac{t _{i} - t _{j}}{2 τ _{decor}})) \frac{l n ( 10 )}{2}} .

E [1 0^{Δ (t_{i} - t_{j})}] = 1 0^{Δ (t_{j}) e x p (- \frac{t _{i} - t _{j}}{4 τ _{decor}}) + σ_{MS}^{2} (1 - e x p (- \frac{t _{i} - t _{j}}{2 τ _{decor}})) \frac{l n ( 10 )}{2}} .

⟨ Δ_{ind} (t_{i})⟩ = lo g \frac{1}{N _{j}} t_{j} = 1 \sum N_{j} 1 0^{Δ (t_{j}) e x p (- \frac{t _{i} - t _{j}}{4 τ _{decor}}) + σ_{MS}^{2} (1 - e x p (- \frac{t _{i} - t _{j}}{2 τ _{decor}})) \frac{l n ( 10 )}{2}} .

⟨ Δ_{ind} (t_{i})⟩ = lo g \frac{1}{N _{j}} t_{j} = 1 \sum N_{j} 1 0^{Δ (t_{j}) e x p (- \frac{t _{i} - t _{j}}{4 τ _{decor}}) + σ_{MS}^{2} (1 - e x p (- \frac{t _{i} - t _{j}}{2 τ _{decor}})) \frac{l n ( 10 )}{2}} .

σ_{ind, MS}^{2} (t_{i}) = \frac{σ _{MS}^{2}}{N _{j}^{2}} t_{j} = 1 \sum N_{j} t_{j^{^{'}}} = 1 \sum N_{j} A C F (t_{j} - t_{j^{^{'}}}),

σ_{ind, MS}^{2} (t_{i}) = \frac{σ _{MS}^{2}}{N _{j}^{2}} t_{j} = 1 \sum N_{j} t_{j^{^{'}}} = 1 \sum N_{j} A C F (t_{j} - t_{j^{^{'}}}),

σ_{ind, MS} (t_{i}) = \frac{σ _{MS}}{N _{j}} 1 + 2 t_{j} = 1 \sum N_{j} (1 - \frac{t _{j}}{N _{j}}) A C F (t_{j})^{1/2} .

σ_{ind, MS} (t_{i}) = \frac{σ _{MS}}{N _{j}} 1 + 2 t_{j} = 1 \sum N_{j} (1 - \frac{t _{j}}{N _{j}}) A C F (t_{j})^{1/2} .

A C F (t) = \frac{1}{2} [- (t)^{2 H} g (\frac{1}{t}) + (t + 1)^{2 H} + (t - 1)^{2 H} g (\frac{1}{t - 1}) - t^{2 H}]

A C F (t) = \frac{1}{2} [- (t)^{2 H} g (\frac{1}{t}) + (t + 1)^{2 H} + (t - 1)^{2 H} g (\frac{1}{t - 1}) - t^{2 H}]

g (t^{^{'}}) = 1 + 2 H \int_{0}^{t^{^{'}}} ((1 + s)^{H - 1/2} - s^{H - 1/2})^{2} d s,

g (t^{^{'}}) = 1 + 2 H \int_{0}^{t^{^{'}}} ((1 + s)^{H - 1/2} - s^{H - 1/2})^{2} d s,

A C F (t) = \frac{1}{2} ((t + 1)^{2 H} - 2 (t)^{2 H} + (t - 1)^{2 H}) .

A C F (t) = \frac{1}{2} ((t + 1)^{2 H} - 2 (t)^{2 H} + (t - 1)^{2 H}) .

SFR_{ind} (t) = SFR (t^{'}) f_{ind} (t - t^{'}),

SFR_{ind} (t) = SFR (t^{'}) f_{ind} (t - t^{'}),

SFR_{ind} (t) = \int_{- \infty}^{t} SFR (t^{'}) f_{ind} (t - t^{'}) d t^{'},

SFR_{ind} (t) = \int_{- \infty}^{t} SFR (t^{'}) f_{ind} (t - t^{'}) d t^{'},

SFR_{ind} (t_{i}) = t_{j} = 0 \sum \infty SFR (t_{i} - t_{j}) w_{ind} (t_{j}),

SFR_{ind} (t_{i}) = t_{j} = 0 \sum \infty SFR (t_{i} - t_{j}) w_{ind} (t_{j}),

w_{ind} (t_{j}) = \int_{t_{j} - δ /2}^{t_{j} + δ /2} f_{ind} (t^{'}) d t^{'},

w_{ind} (t_{j}) = \int_{t_{j} - δ /2}^{t_{j} + δ /2} f_{ind} (t^{'}) d t^{'},

Δ = lo g SFR - lo g SFR_{MS} .

Δ = lo g SFR - lo g SFR_{MS} .

Δ_{ind} (t_{i}) = lo g t_{j} = 0 \sum \infty 1 0^{Δ (t_{i} - t_{j})} w_{ind} (t_{j}) .

Δ_{ind} (t_{i}) = lo g t_{j} = 0 \sum \infty 1 0^{Δ (t_{i} - t_{j})} w_{ind} (t_{j}) .

1 0^{x} = 1 + ln (10) x + \frac{1}{2} ln (10)^{2} x^{2} + \frac{1}{6} ln (10)^{3} x^{3} + ...

1 0^{x} = 1 + ln (10) x + \frac{1}{2} ln (10)^{2} x^{2} + \frac{1}{6} ln (10)^{3} x^{3} + ...

Δ_{ind} (t_{i}) = lo g j = 0 \sum \infty (1 + ln (10) Δ (t_{i} - t_{j}) + \frac{1}{2} ln (10)^{2} Δ^{2} (t_{i} - t_{j}) + ....) w_{ind} (t_{j}) = lo g (j = 0 \sum \infty w_{ind} (t_{j})) + ln (10) t_{j} = 0 \sum \infty Δ (t_{i} - t_{j}) w_{ind} (t_{j}) + \frac{1}{2} ln (10)^{2} t_{j} = 0 \sum \infty Δ^{2} (t_{i} - t_{j}) w_{ind} (t_{j}) + .... = lo g 1 + ln (10) t_{j} = 0 \sum \infty Δ (t_{i} - t_{j}) w_{ind} (t_{j}) + \frac{1}{2} ln (10)^{2} t_{j} = 0 \sum \infty Δ^{2} (t_{i} - t_{j}) w_{ind} (t_{j}) + .... .

Δ_{ind} (t_{i}) = lo g j = 0 \sum \infty (1 + ln (10) Δ (t_{i} - t_{j}) + \frac{1}{2} ln (10)^{2} Δ^{2} (t_{i} - t_{j}) + ....) w_{ind} (t_{j}) = lo g (j = 0 \sum \infty w_{ind} (t_{j})) + ln (10) t_{j} = 0 \sum \infty Δ (t_{i} - t_{j}) w_{ind} (t_{j}) + \frac{1}{2} ln (10)^{2} t_{j} = 0 \sum \infty Δ^{2} (t_{i} - t_{j}) w_{ind} (t_{j}) + .... = lo g 1 + ln (10) t_{j} = 0 \sum \infty Δ (t_{i} - t_{j}) w_{ind} (t_{j}) + \frac{1}{2} ln (10)^{2} t_{j} = 0 \sum \infty Δ^{2} (t_{i} - t_{j}) w_{ind} (t_{j}) + .... .

lo g (1 + x) = \frac{x}{ln ( 10 )} - \frac{x ^{2}}{2 ln ( 10 )} + \frac{x ^{3}}{3 ln ( 10 )} + ...

lo g (1 + x) = \frac{x}{ln ( 10 )} - \frac{x ^{2}}{2 ln ( 10 )} + \frac{x ^{3}}{3 ln ( 10 )} + ...

Δ_{ind} (t_{i}) = t_{j} = 0 \sum \infty Δ (t_{i} - t_{j}) w_{ind} (t_{j}) + \frac{ln ( 10 )}{2} t_{j} = 0 \sum \infty Δ^{2} (t_{i} - t_{j}) w_{ind} (t_{j}) - t_{j} = 0 \sum \infty Δ (t_{i} - t_{j}) w_{ind} (t_{j})^{2} + ...

Δ_{ind} (t_{i}) = t_{j} = 0 \sum \infty Δ (t_{i} - t_{j}) w_{ind} (t_{j}) + \frac{ln ( 10 )}{2} t_{j} = 0 \sum \infty Δ^{2} (t_{i} - t_{j}) w_{ind} (t_{j}) - t_{j} = 0 \sum \infty Δ (t_{i} - t_{j}) w_{ind} (t_{j})^{2} + ...

σ_{ind, MS}^{2} (t_{i}) = \mbox V A R t_{j} = 0 \sum \infty Δ (t_{i} - t_{j}) w_{ind} (t_{j}) = t_{j} = 0 \sum \infty t_{j^{^{'}}} = 0 \sum \infty E (Δ (t_{i} - t_{j}) w_{ind} (t_{j}) Δ (t_{i} - t_{j^{^{'}}}) w_{ind} (t_{j^{^{'}}})) - E (Δ (t_{i} - t_{j}) w_{ind} (t_{j})) E (Δ (t_{i} - t_{j^{^{'}}}) w_{ind} (t_{j^{^{'}}})) = t_{j} = 0 \sum \infty t_{j^{^{'}}} = 0 \sum \infty w_{ind} (t_{j}) w_{ind} (t_{j^{^{'}}}) E (Δ (t_{i} - t_{j}) Δ (t_{i} - t_{j^{^{'}}})) = t_{j} = 0 \sum \infty t_{j^{^{'}}} = 0 \sum \infty w_{ind} (t_{j}) w_{ind} (t_{j^{^{'}}}) \cdot C O V (Δ (t_{i} - t_{j}), Δ (t_{i} - t_{j^{^{'}}})) = σ_{MS}^{2} t_{j} = 0 \sum \infty t_{j^{^{'}}} = 0 \sum \infty w_{ind} (t_{j}) w_{ind} (t_{j^{^{'}}}) A C F (t_{j} - t_{j^{^{'}}}),

σ_{ind, MS}^{2} (t_{i}) = \mbox V A R t_{j} = 0 \sum \infty Δ (t_{i} - t_{j}) w_{ind} (t_{j}) = t_{j} = 0 \sum \infty t_{j^{^{'}}} = 0 \sum \infty E (Δ (t_{i} - t_{j}) w_{ind} (t_{j}) Δ (t_{i} - t_{j^{^{'}}}) w_{ind} (t_{j^{^{'}}})) - E (Δ (t_{i} - t_{j}) w_{ind} (t_{j})) E (Δ (t_{i} - t_{j^{^{'}}}) w_{ind} (t_{j^{^{'}}})) = t_{j} = 0 \sum \infty t_{j^{^{'}}} = 0 \sum \infty w_{ind} (t_{j}) w_{ind} (t_{j^{^{'}}}) E (Δ (t_{i} - t_{j}) Δ (t_{i} - t_{j^{^{'}}})) = t_{j} = 0 \sum \infty t_{j^{^{'}}} = 0 \sum \infty w_{ind} (t_{j}) w_{ind} (t_{j^{^{'}}}) \cdot C O V (Δ (t_{i} - t_{j}), Δ (t_{i} - t_{j^{^{'}}})) = σ_{MS}^{2} t_{j} = 0 \sum \infty t_{j^{^{'}}} = 0 \sum \infty w_{ind} (t_{j}) w_{ind} (t_{j^{^{'}}}) A C F (t_{j} - t_{j^{^{'}}}),

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Stochastic modeling of star-formation histories I:

the scatter of the star-forming main sequence

Neven Caplar,1 Sandro Tacchella,2

1Department of Astrophysical Sciences, Princeton University, 4 Ivy Ln., Princeton, NJ 08544, USA

2Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA E-mail: [email protected]: [email protected]

(Accepted XXX. Received YYY; in original form ZZZ)

Abstract

We present a framework for modelling the star-formation histories of galaxies as a stochastic process. We define this stochastic process through a power spectrum density with a functional form of a broken power-law. Star-formation histories are correlated on short timescales, the strength of this correlation described by a power-law slope, $\alpha$ , and they decorrelate to resemble white noise over a timescale that is proportional to the timescale of the break in the power spectrum density, $\tau_{\rm break}$ . We use this framework to explore the properties of the stochastic process that, we assume, gives rise to the log-normal scatter about the relationship between star-formation rate and stellar mass, the so-called galaxy star-forming main sequence. Specifically, we show how the measurements of the normalisation and width ( $\sigma_{\rm MS}$ ) of the main sequence, measured in several passbands that probe different timescales, give a constraint on the parameters of the underlying power spectrum density. We first derive these results analytically for a simplified case where we model observations by averaging over the recent star-formation history. We then run numerical simulations to find results for more realistic observational cases. As a proof of concept, we use observational estimates of the main sequence scatter at $z\sim 0$ and $M_{\star}\approx 10^{10}~{}M_{\odot}$ measured in H $\alpha$ , UV+IR and the u-band, and show that combination of these point to $\tau_{\rm break}=178^{+104}_{-66}$ Myr, when assuming $\alpha=2$ . This implies that star-formation histories of galaxies lose “memory” of their previous activity on a timescale of $\sim 200$ Myr, highlighting the importance of baryonic effects that act over the dynamical timescales of galaxies.

keywords:

galaxies: evolution – galaxies: star formation – galaxies: statistics

††pubyear: 2019††pagerange: Stochastic modeling of star-formation histories I: the scatter of the star-forming main sequence–D

1 Introduction

Recent observational studies have produced a wealth of data on galaxy populations at different epochs. This has enabled us to study galaxies in a statistical manner and explore scaling relations between a wide range of galaxy properties, such as luminosity, color, stellar mass, size, velocity dispersion, and central stellar mass density, with cosmic time (e.g., Faber & Jackson, 1976; Kormendy, 1977; Djorgovski & Davis, 1987; Baldry et al., 2004; van der Wel et al., 2014; Mosleh et al., 2017; Barro et al., 2017). However, how and why individual galaxies change their physical properties and evolve to create such scaling relations is still unclear, since we observe each individual galaxy only once during its cosmic lifetime.

A possible way forward is the use of analytical and numerical models. A large body of theoretical work, including cosmological simulations, has been developed aiming to explain the evolving galaxy properties and constrain physical mechanisms that give rise to the observations (e.g., Vogelsberger et al., 2014; Hopkins et al., 2014; Ceverino et al., 2014; Porter et al., 2014; Henriques et al., 2015; Schaye et al., 2015; Pillepich et al., 2018). Currently, such models still have relatively limited predictive power because of the inadequate spatial/mass/temporal resolution, producing small galaxy samples and/or the uncertainties associated with the “sub-grid” physics (e.g., Somerville & Davé, 2015).

Another way forward is based on observations: precise evaluation of the number densities of different types of galaxies, coupled with continuity equations, can be used to construct high-fidelity phenomenological models that explain a large number of observed features in the Universe (Peng et al., 2010; Leitner, 2012; Carollo et al., 2013; Damjanov et al., 2014; Kelson, 2014; Caplar et al., 2015; Kelson et al., 2016; Caplar et al., 2018). Additionally, galaxy stellar ages can be used to understand how galaxies relate to their precursors at earlier epochs to constrain, for example, the mechanism responsible for the change in morphology (Fagioli et al., 2016; Williams et al., 2017; Tacchella et al., 2017; Zahid & Geller, 2017; Damjanov et al., 2018).

More general than simple “ages” are star-formation histories (SFH), which are one of the critical ingredients needed to understand the evolution of galaxies. Having the full information on how galaxies change their star-formation rate (SFR) as a function of cosmic time would allow us to project the galaxy back in time and understand on which timescales the galaxy grew its mass as a consequence of star formation. Furthermore, SFHs encode information about the variability on short and long timescales that arise from different physical mechanisms. There are several physical processes – acting on different timescales – that can increase the SFR, such as gas accretion or mergers, and that can decrease the SFR, such as ram pressure stripping of gas or feedback from black holes and stars.

SFHs are related to the scaling relation between the SFRs and stellar masses ( $M_{\star}$ ) of galaxies: observational studies have produced a consensus that since at least redshift 3, the majority of star-forming galaxies have a SFR that is strongly correlated with their mass ( $\mathrm{SFR}\propto M_{\star}^{\beta}$ , with $\beta\sim 1$ ), creating the so-called star-forming main sequence (MS). It has been established that the MS normalization increases with the look-back time as $(1+z)^{\sim 2.2}$ (Brinchmann et al., 2004; Daddi et al., 2007; Elbaz et al., 2007; Noeske et al., 2007; Whitaker et al., 2012; Speagle et al., 2014; Schreiber et al., 2015; Boogaard et al., 2018). This increase of the overall normalization of the MS with lookback time can be understood by the higher dark matter accretion rate onto halos, which leads to higher gas fractions in high-redshift galaxies (Lilly et al., 2013; Tacconi et al., 2013; Tacchella et al., 2013, 2018; Genzel et al., 2015; Rodríguez-Puebla et al., 2016).

The most noticeable feature is that the MS relation at any given redshift shows a rather small scatter of $\sigma_{\rm MS}\sim 0.2-0.4~{}\mathrm{dex}$ , depending on the exact definition of star-forming galaxies and the method used to infer SFRs and stellar masses. This rather small scatter has been used to argue that the SFRs of galaxies on the MS are sustained for extended periods of time in a quasi-steady state of gas inflow, gas outflow, and gas consumption (Bouché et al., 2010; Daddi et al., 2010; Genzel et al., 2010; Davé et al., 2012; Lilly et al., 2013; Dekel & Mandelker, 2014; Tacchella et al., 2016; Rodríguez-Puebla et al., 2016). Clearly, the MS scatter is directly related to the variability in the SFH and hence encodes physical mechanisms acting on different timescales. In general terms, Abramson et al. (2015, see also ) highlights that if the MS scatter arises due to short-term fluctuations, star-forming galaxies with similar mass mostly grew-up together, while if the MS scatter arises due to long-term fluctuations (Peng et al., 2010; Behroozi et al., 2013) then star-forming galaxies with similar mass did not grow-up together and key physics lies in the process that diversifies star formation histories (Gladders et al., 2013; Kelson, 2014).

Based on cosmological zoom-in simulations (VELA simulations; Ceverino et al. 2014; Zolotov et al. 2015), Tacchella et al. (2016) suggested that star-forming galaxies with $M_{\star}\approx 10^{9}-10^{11}~{}M_{\odot}$ oscillate about the MS ridgeline on timescales of $\leq 2$ Gyr at $z\sim 1-4$ . The propagation of galaxies upwards toward the upper envelope of the MS is due to gas compaction, triggered by e.g., mergers, counter-rotating streams, and/or violent disk instabilities. The downturn at the upper envelope is due to central gas depletion by peak star formation and outflows while inflow from the shrunken gas disk is suppressed. Rodríguez-Puebla et al. (2016) argue, based on advanced abundance matching, that the MS scatter could be set by the halo mass accretion rate, which itself has a scatter of $\sim 0.3$ dex, matching the observed dispersion of the MS. Matthee & Schaye (2018) show, based on the EAGLE simulations, that the scatter of the MS at $z\sim 0$ originates from a combination of short timescales ( $\la 1~{}\mathrm{Gyr}$ ) that are presumably associated with self-regulation from cooling, star formation and outflows, and long time-scale ( $\sim 10~{}\mathrm{Gyr}$ ) variations related to differences in halo formation times. Torrey et al. (2018) show, based on the IllustrisTNG simulations, that galaxy offsets from the MS and mass-metallicity relation oscillate over similar time-scales (roughly the dynamical time of the dark matter halo) and are often anticorrelated. They speculate that the SFR and metallicity evolution tracks may become decoupled in galaxy formation models dominated by feedback-driven globally bursty SFHs, which could weaken the FMR residual correlation strength.

It is well-known that measurements of star-formation rate in different observed bands carry different information, with shorter wavelengths carrying more information about recent star-formation, and longer wavelengths being more sensitive to star-formation on longer timescales. However, measuring SFHs from the observed spectral energy distribution (SED) of galaxies is still very challenging (Conroy, 2013). State-of-the-art SED fitting codes are nowadays able to fit non-parametric SFHs to the galaxy SEDs and provide constraints of the SFHs on timescales of 100 Myr to a few Gyrs, i.e., on rather long timescales (Thomas et al., 2005; Graves & Schiavon, 2008; Davies et al., 2015; Pacifici et al., 2016; Iyer & Gawiser, 2017; Leja et al., 2018a; Carnall et al., 2018). Short timescales can in principle be inferred from comparing nebular lines to the ultraviolet (UV) continuum. For instance, Weisz et al. (2012) and Guo et al. (2016) inferred a systematic decline and increased scatter in nebular line to FUV ratios toward low-mass systems. They interpreted their data with simple SFH models that included bursts lasting tens of Myr, and with periods of $\sim 250$ Myr, finding that these low mass galaxies are rather “bursty”. However, large uncertainties remain, including variations of the initial mass function (IMF), differential dust attenuation distribution in front of young and old stars, and metallicity effects (e.g., Shivaei et al., 2018).

All this research implies that galaxies exhibits wide range of SFHs, exhibiting bursts, drops and periods non-changing star-formation. In this work, we propose to model the time-dependence of star formation and the movement of star-forming galaxies around the MS ridgeline as a purely stochastic process. We aim to describe the stochastic behaviour in very general terms and therefore define its properties in the frequency domain through a power spectrum density (PSD) that is modelled as a broken power-law. The high-frequency slope of the PSD determines how quickly the SFR changes on short timescales and is going to be connected with physical drivers of star-formation. The broken power-law form allows for a break in the correlation, i.e., sets up a timescale on which the SFR in a single galaxy loses “memory” of previous star-formation.

We note that this timescale can, in principle, be longer than the timescale of the Universe. In that case, the PSD is effectively a pure power-law for all practical purposes and the star-formation in a given galaxy has a “long-term dependence”, i.e., it is correlated over the cosmic timescales. This process, with a single power-law slope, is used by Kelson (2014) to model stellar mass growth in the Universe. He described the SFHs using a so-called “Hurst parameter”, which is very closely related to the slope of the PSD used in this work. We discuss the connection between describing the stochastic process via the PSD with a broken power law (and its high-frequency slope) and the description using “Hurst parameter” more extensively in Appendix A.

In this paper we focus on the Myr to a Gyr timescale and exclude the longer timescales in the analysis, which we believe set the overall evolution of the MS. Characterisation of the PSD on these “short” scales is critical for a couple of reasons. Firstly, it provides a test for simulations, which aim to match various distributions (e.g., galaxy mass function) over a span of redshifts, but do not necessarily aim to match stochasticity of star-formation with cosmic time. Secondly, the recovered parameters of the underlying stochastic process are connected to the physical drivers of star-formation and can provide valuable insight into the properties of star-formation. In particular, the timescale of the break in the PSD should be connected with the timescale on which star-formation due to internal galaxy properties becomes less important compared to external feeding of the galaxy. For an example in a related field, stochastic modelling of the short-term variability ( $\sim$ 1 year) of active galactic nuclei (AGN) is one the main methods to explore physical processes powering AGN (Kelly et al., 2009; MacLeod et al., 2010; Dexter & Agol, 2011; MacLeod et al., 2012; Simm et al., 2016; Caplar et al., 2017), and is providing serious challenges to standard viscous accretion disk models commonly used to model AGN accretion (see review by Lawrence 2018). Finally, understanding the variability of the SFR on short timescale is crucial for proper SED fitting, as the characteristic timescale on which SFR varies is directly related to the size of parameter confidence intervals (Leja et al., 2018b).

In this work we lay the basic groundwork of our approach and conduct a simple analysis of the observations in order to constrain the PSD. We focus our investigation on rather massive star-forming galaxies with $M_{\star}\approx 10^{10}-10^{10.5}~{}M_{\odot}$ . We work in this rather narrow mass range because measurement errors are the small, fraction of galaxies that are quenching is negligible, and result are being less affected by possible variations of dust attenuation with galaxies mass (as elaborated in Section 5). In a companion paper (Iyer et al. in prep.), we analyse and measure the PSDs for a wide range of numerical (cosmological and zoom-in simulations) and semi-analytical models. We find that PSDs are well described by a broken power-law, where the short timescales (high frequencies) follow $\propto f^{\alpha}$ and above a certain decorrelation timescale, $\tau_{\rm break}$ , the PSD is constant. Nearly all of the different theoretical models produce $\alpha\approx-2$ , while having different values for $\tau_{\rm break}$ , ranging from $\tau_{\rm break}\sim 100~{}\mathrm{Myr}$ in the FIRE simulations to $\tau_{\rm break}\sim 1000~{}\mathrm{Myr}$ in IllustrisTNG. We show that the PSD is set by physical ingredients in the models and, hence, the PSD contains a wealth of information that has not yet been exploited. In a separate, upcoming paper we will address many issues that complicate the precise determination of the PSD in observations, including possible variations in the IMF and dust attenuation. There, we also wish to use observational constraints to analyse properties of galaxies at higher redshifts, where the proper consideration of these complications is critically important.

This paper is structured as follows. In Section 2 we introduce our model and give a short overview of the properties of stochastic processes. In Section 3 we show how we can infer parameters describing movement of galaxies about the MS. We show both analytical and numerical results in the simplified case in which response of each star-formation indicator is given by a step function. In Section 4 we show predictions for the observational results given the more realistic response functions. In Section 5 we use results developed in previous sections and determine, as a proof of concept, parameters describing stochastic processes by using observations of the measured scatter of the MS in the local universe. We conclude in Section 6. In the appendix we expand on several more technical aspects of the paper.

Throughout the paper we use the terms “scatter of MS”, “standard deviation of MS” and “width of MS” interchangeably, and denote it by $\sigma_{\rm MS}$ . We use the term “intrinsic” MS, or “momentary” MS to denote width of the MS that would be detected if we could measure the star-formation over the previous 1 Myr. We expand on our choice for 1 Myr in Section 2.3. We make our codes available at github.com/stacchella/ and github.com/nevencaplar/.

2 Description of star-formation histories with stochastic processes

In this section, we introduce our approach to modelling SFHs with stochastic processes. In particular, we start by highlighting and justifying the assumptions we make in our model. We then explain what a stochastic process is and how it can be described mathematically. Finally, we show how we generate SFHs with a stochastic process and discuss their properties such as burstiness.

2.1 Assumptions

We propose to model time-dependence of the SFR relative to the MS ridgeline of star-forming galaxies as a purely stochastic process. We denote the SFR relative to the MS ridgeline for a galaxy with stellar mass $M_{\star}$ by $\Delta_{\rm MS}=\log(\mathrm{SFR}/\mathrm{SFR}_{\rm MS}(M_{\star}))$ , where $\mathrm{SFR}_{\rm MS}(M_{\star})$ is the SFR of the MS at mass of $M_{\star}$ . Throughout this work, we assume that:

•

there exists a MS, i.e., there is a power-law relation between SFR and stellar mass for star-forming galaxies that are separate from the quiescent galaxy population;

•

the distribution of SFRs of galaxies about the MS ridgeline is well described by a log-normal function;

•

time-variability of $\Delta_{\rm MS}$ of all star-forming galaxies of a given mass can be described with the same stochastic process given by a broken power-law PSD.

The first two assumptions are well justified from observations. Specifically, the existence of the MS have been put forward by several different observing groups (e.g., Noeske et al., 2007; Elbaz et al., 2007; Daddi et al., 2007; Peng et al., 2015). Furthermore, several groups claim that the scatter can be well described by a log-normal (Guo et al., 2013; Schreiber et al., 2015; Chang et al., 2015; Davies et al., 2019). Feldmann (2017) shows that the whole galaxy population, i.e., star-forming and quiescent galaxies, can be well fitted by a zero-inflated negative binomial distribution, which is only a small modification from the log-normal distribution describing star-forming galaxies, but has the advantage that no division between star-forming and quiescent galaxies is needed. There is still debate in the literature about the exact shape of the MS, including its slope and turnover at high mass at low redshifts, and also about the existence of a MS at high redshifts ( $z>3$ ; Speagle et al. 2014; Whitaker et al. 2014; Iyer et al. 2018). Since we focus in this paper on moderately massive ( $\approx 10^{10}~{}M_{\odot}$ ) galaxies at $z\sim 0$ and model the SFRs relative to the MS, this uncertainty on the exact shape is not critical.

The third assumption that all massive star-forming galaxies are described with the same stochastic processes is the simplest reasonable assumption and it enables us to maintain analytical tractability of the model. In this work we do not explore further variables that could influence the PSD, such as environment. It is also possible that lower-mass galaxies (including dwarfs) follow a stochastic process with different properties; this will be investigated in the future. Finally, as also mentioned in the introduction, we assume that the stochastic behaviour can be described through a broken power-law PSD, i.e., we assume that the 3-point correlation function (and any higher order terms) is zero. This is again a reasonable assumption given the current observational constraints, as further discussed in Section 5.

2.2 General description of stochastic processes

In this section, we will discuss several different stochastic processes and develop an intuition for parameters that we will use in the rest of paper. We characterise different stochastic processes by their PSD, a quantity denoting the amount of power that a process has as a function of frequency, per unit frequency. The PSD gives a quantitative estimate of the importance of each mode, i.e., the importance of each frequency interval for a given stochastic process. We will also extensively use the Fourier pair of the PSD, the autocorrelation function (ACF). The autocorrelation function provides the same information as the PSD, but in the time domain instead of the frequency domain. The definition of ACF, for a process that leads to a time-series $\Delta(t)$ is:

[TABLE]

where $\sigma^{2}_{\Delta}$ is the long-term variance of the process. We can see that highly correlated timeseries, where $\Delta(t)\sim\Delta(t+\Delta t)$ , give ACF $\sim 1$ , while uncorrelated timeseries will lead to ACF $\sim 0$ .

The simplest stochastic process is white noise. This process has equal power at all frequencies. In this work we denote frequency with $f$ . The statement about equal power on all frequencies can be mathematically written as $PSD(f)\propto f^{\alpha}$ , with $\alpha=0$ . This process is perfectly uncorrelated, i.e., for a white noise process it holds that $ACF(t)=0$ , for all times when $t>0$ .

Other processes can be defined by modifying the slope $\alpha$ of the PSD. By decreasing the slope $\alpha$ the nature of process changes and values separated by a given $\Delta t$ become more similar. Some well known examples are given by $\alpha=-1$ (pink noise) and by $\alpha=-2$ (red noise); the latter corresponding to a random walk or Brownian motion.

All natural process have some characteristic time $\tau$ , after which a process becomes uncorrelated, i.e., loses its “memory” of the previous values and behaves like a white noise process. We can express that statement in terms of the ACF by saying that the ACF at times longer than the decorrelation time, $\tau_{\rm decor}$ has to be very small, $ACF(t>\tau_{\rm decor})\ll 1$ , or, in the closely related terms of the PSD, by requiring that the slope at small frequencies below “break” frequency, $f_{\rm break}$ , has to be close to zero, $\alpha(f<f_{\rm break})\approx 0$ . In that case, the PSD has a broken power-law shape, with two distinct slopes at large and small frequencies.

To illustrate this point further, let’s consider the simplest and most instructive case where the PSD is described by a broken power-law, which corresponds to a damped random walk. The damped random walk is defined by the solution to the stochastic differential equation (Kelly et al., 2015):

[TABLE]

where $\epsilon(t)$ is a continuous white-noise process with zero mean and variance $\sigma_{\rm int}^{2}$ . In this form, one can recognise $\sigma_{\rm int}^{2}$ as an intrinsic amount of “randomness” entering the system, while $\tau_{\rm decor}$ describes the timescale by which the process gravitates towards a stable-state solution. This process is also known as the Orstein-Uhlenbeck process in physics, or the Vasicek model in financial literature.

The damped random walk process leads to the following PSD:

[TABLE]

where $\sigma=\sigma_{\rm int}\tau_{\rm decor}$ is the “long-term” variability and $\tau_{\rm break}=2\pi\tau_{\rm decor}$ . The connection of the description of the process in frequency and time domain can be achived via the Wiener-Khinchin theorem (Emmanoulopoulos et al., 2010). Application of this procedure for this particular choice of PSD leads to

[TABLE]

We see that this process has all the properties that we mentioned above. At short timescales (high frequencies), it is very well correlated (ACF $\sim$ 1) and described with a PSD with a slope that equals to 2. Then, the process becomes rapidly decorrelated at a timescale $\tau_{\rm decor}$ , and at long timescales (low frequencies), the ACF drops quickly (ACF $\sim$ 0) and the PSD is well described with a flat slope. Observationally, high values of ACF imply similar measurements, while decorrelation implies little to no similarity between two measurements. As such, while at short timescales star-formation measurements will be very similar, at timescales larger than $\tau_{\rm decor}$ the observations of the same galaxy will have star-formation measurements which are largely independent of each other.

In this work, we consider a larger group of processes, described with

[TABLE]

This choice offers additional flexibility to encompass a wide range of physical processes that are not necessarily described by random walks, providing a more general description of the possible underlying stochastic process. We show examples of generated SFHs relative to the MS, $\Delta_{\rm MS}$ , in Figure 1. The generation of these SFHs is described in the next section. Figure shows simulated SFHs as a function of high-frequency slopes of the PSD ( $\alpha$ in Equation 5)111In the rest of the paper $\alpha$ will appear as a positively defined quantity, i.e., $\alpha$ is defined with Equation (5), and not with $PSD(f)\propto f^{\alpha}$ . and break timescales ( $\tau_{\rm break}$ in Equation 5). The figure highlights large qualitative differences between SFHs generated with different parameters. When the slope of the PSD is shallow and/or $\tau_{\rm break}$ is small, the SFH rapidly oscillates around the mean and is very quickly decorrelated. On the other hand, steeper slopes and longer $\tau_{\rm break}$ lead to slower oscillations and more correlated behaviour, as described above. In Section 2.4, we will return to this point in order to quantify the “burstiness” of these processes.

All of the examples shown in Figure 1 have been created with $\sigma=0.4~{}\mathrm{dex}$ . We use this value motivated by the measurement of the MS width in Davies et al. (2019) in H $\alpha$ , and we show in Section 4 that the measurement in H $\alpha$ is largely representative of the intrinsic width of MS, $\sigma_{\rm MS}$ . A number of results that we will derive depend somewhat on the intrinsic width of MS. In the rest of the paper, we work with $\sigma_{\rm MS}=0.4~{}\mathrm{dex}$ , but note at each point where our results depend on the exact choice of this value.

Finally, we note that we are implicitly assuming that the stochastic process that we are using is stationary. A stationary process has the property that its mean, variance and autocorrelation function are not time dependent. That the variance, i.e., the width of MS is not time dependent over last several Gyr is consistent with observational data and simulations (Noeske et al. 2007; Rodighiero et al. 2010; Katsianis et al. in prep.). In other words, we assume that the properties of stochastic process do not change with time and that is valid to use the same parameters to describe observations of galaxies over the whole time range probed. It is important to note that this does not mean that we assume that the MS is not allowed to change with time. The main stochastic process modelled in this work is movement of galaxies about the MS ridgeline ( $\Delta_{\rm MS}$ ). The MS itself evolves over cosmological timescales, but because we consider galaxies at $z=0$ , this evolution is negligible. We further discuss this effect in Section 5, when comparing the results from our modelling with the observational data.

2.3 Generating star-formation histories

Here we describe how we generate numerically a sample of SFHs that we use in the rest of the paper when we analyze properties of stochastic processes and observational consequences. For each choice of break timescale $\tau_{\rm break}$ and high-frequency slope $\alpha$ we simulate 1000 SFHs using the algorithm by Timmer & Koenig (1995) and implemented in Python by Connolly (2015). For each galaxy, we only consider the SFR as an integrated quantity over the whole galaxy, i.e., we measure only total SFRs over the whole galaxy without deeper considerations of possible difference in different galactic (sub-)components.

When simulating SFHs we need to assume some minimal timescale at which to separate our measurements of SFR (time bins). In this work we set this timescale to 1 Myr. This choice is motivated by the free-fall time and the lifetime of giant molecular clouds which are at least of that order (Murray et al., 2010; Hollyhead et al., 2015; Freeman et al., 2017; Grudić et al., 2018). Operationally this means that we are assuming that below the scale of 1 Myr, the SFR is perfectly correlated, i.e., remains constant. This is a reasonable assumption since we are considering only total SFRs integrated over the whole galaxy, which itself is made up of several star-forming events (giant molecular clouds). In terms of PSD that means that the domain of the PSD is defined up to scales of $f<1/\mbox{Myr}$ . In terms of ACF this leads to a discontinuity in our assumed ACF at 1 Myr, where we set ACF to be equal to unity. The effect of this assumptions is less and less pronounced when the stochastic process has larger $\alpha$ and $\tau_{\rm break}$ , since in this case the process is already highly correlated on 1 Myr scales, and our assumption is satisfied by construction. In Section 5, when we compare with observational data at low redshifts, the deduced parameters are outside the range that is critically sensitive to this assumption, but this might not be the case for studies at higher redshifts where deeper investigation of this assumption should be conducted.

For each SFH we simulate 4282 points, corresponding to 4.282 Gyr, i.e., to a galaxy evolving from $z=0.4$ to $z=0$ . We note that choice of the length is somewhat aribtary and carries no great physical significance. We often use much shorter SFHs and our results are not dependent on the actual duration of the simulated histories, provided they are significantly longer than the $\tau_{\rm break}$ of the considered process. In Figure 2 we illustrate the results of this procedure. For each $\tau_{\rm break}$ and $\alpha$ , we show a selection of three SFHs relative to the MS (left panel), together with their PSDs and ACFs (middle and right panels). We also show median relations in black, which have been deduced from all 1000 simulated galaxies. We also explicitly show the effect of a cut at 1 Myr, seen as the grey region in the plots of the PSD, and seen in the fact that ACFs extended “only” to $\log(\mathrm{time}/\mathrm{Myr})=0$ .

2.4 Quantifying the burstiness of star-formation histories

Here we want to connect the parameters that we have introduced ( $\alpha$ and $\tau_{\rm break}$ ) with a more intuitive quantity which describes the tendency of the SFH to have “bursts”, i.e., periods of star formation which are much higher (or lower) than the mean of the population. From Figures 1 and 2 it is apparent that the power-law slope $\alpha$ and the break timescale $\tau_{\rm break}$ are connected with the variability, or “burstiness”, of the SFH. Following Goh & Barabási (2008, see also ), we define the burstiness parameter, for a single galaxy, as

[TABLE]

where $\sigma$ is the standard deviation of $\Delta_{\rm MS}$ , and $\mu$ is the mean of $\Delta_{\rm MS}$ of an individual SFH, measured over a certain time-interval. It is important to note that the exact value of $B$ depends on the interval over which it is measured. We choose this time interval to be 30 Myr in order to be able to describe the large diversity of variability in the SFHs seen in Figure 1. This choice of time interval is also consistent with the one used by Applebaum et al. (2018). In general, $B$ has a value in the bounded range $(-1,1)$ , and its magnitude correlates with the signal’s burstiness. The average deviation from the MS for one galaxy, $\mu$ , will approach zero as the SFH becomes more bursty. Therefore, $B=1$ is a maximally bursty signal, $B=0$ is neutral, and $B=-1$ corresponds to a completely regular (periodic or constant) signal.

To show trends as a function of $\alpha$ and $\tau_{\rm break}$ , we compute the mean of $B$ over our sample of 1000 numerically generated model galaxies for each point of the $\alpha$ - $\tau_{\rm break}$ parameter space. Figure 3 shows how the burstiness $B$ varies as a function of the PSD parameters $\alpha$ and $\tau_{\rm break}$ . Large values of $\alpha$ and $\tau_{\rm break}$ correspond to slowly varying, nearly constant SFHs (see Figure 1), corresponding to $B=-1$ , i.e., highly non-bursty. Most individual galaxies rarely venture to the edges of the MS and therefore experience no episodes of very high (or very low) SFR in comparison to the average over 30 Myr. On the other hand, PSDs with small $\alpha$ and short $\tau_{\rm break}$ produce SFHs with $B=1$ (highly bursty). These SFH rapidly change their SFR and have very varied SFHs over the considered time range of 30 Myr. Furthermore, Figure 3 highlights the degeneracy between $\alpha$ and $\tau_{\rm break}$ : contours of constant $B$ extend along well-defined curves in the $\alpha$ - $\tau_{\rm break}$ plane. This is something that we will encounter again in sections 3 and 4, where we investigate how the scatter of the MS as seen from different indicators depends on the PSD parameters $\alpha$ and $\tau_{\rm break}$ .

3 Finding parameters of stochastic star-formation history in a toy model

In the previous section we have shown how different choices for parameters describing the PSD (namely the slope $\alpha$ and the break timescale $\tau_{\rm break}$ ) lead to different SFHs of star-forming galaxies. Estimating these stochastic parameters from observations is difficult, since it is still a big challenge to infer the SFH from a galaxy’s SED (see Introduction). By measuring the SFR using different indicators we are able to estimate how the SFR varied in the galaxy’s recent past (e.g., Kennicutt, 1998). Specifically, measurements of the SFR with indicators that are based on nebular lines will probe very recent star-formation ( $\sim 10^{6}~{}\mathrm{yr}$ ), while measurement in the UV and optical bands will be sensitive to much longer timescales ( $\sim 10^{7}~{}\mathrm{yr}$ and $\sim 10^{8}~{}\mathrm{yr}$ , respectively; see also Section 4). As we will show below, by studying a large sample of galaxies for which the SFRs have been estimated with different indicators, the parameters of the underlying stochastic process can be deduced.

To develop understanding, we present the simplest possible case in which measurements with different indicators correspond exactly to averaging the recent SFR over a given timescale. We will develop analytical expressions for the offset as well as scatter of the MS as a function of the averaging timescale. This approach will help us to see connections and correlations between different parameters more clearly. Afterwards, in Section 4, we discuss realistic response functions for different SFR indicators and redo much of the calculations with more complex and more realistic models. In contrast with mainly analytical approach in this section, when dealing with more realistic cases we will derive results with greater emphasis on numerical calculations.

3.1 Impact of averaging of the star-formation histories on the shape of the measured MS

In this section we show how the measured normalization and width of the MS change as a function of the timescale over which the SFH is averaged. These changes depend on the the properties of the stochastic process (i.e., PSD) and, hence, can also be used to infer these properties. We first show analytically derived relations between the measured offset from the MS and then the expression for the width of the MS. We then confirm these results by our numerically generated model galaxies (Section 2.3) and show how parameters of stochastic process can be deduced from measurements.

3.1.1 Mean measured offset from the main sequence

We first consider mean offsets from the MS, as a function of the parameters of the stochastic process, when measured by indicators probing different timescales. We denote the SFR that has been measured with a given indicator by $\mathrm{SFR_{\rm ind}}$ . As already mentioned, we assume now the simplest possible case in which the response of each indicator is given by a step function each with different lengths, i.e., different response times. In that case, the response step function can be divided into $N_{j}$ equal time steps. The measurement of the SFR with a given indicator at time $t_{i}$ is then equal to the averaged “true” SFR of the previous $N_{j}$ steps, where $t_{j}$ are labels for these equidistant time steps:

[TABLE]

Let’s now define the logarithmic distance of the SFR from the ridgeline of the intrinsic MS222We henceforth use just $\Delta$ as shorthand for $\Delta_{\rm MS}$ in equations, in order to reduce the size and complexity of our expressions.

[TABLE]

and its equivalent quantity, logarithmic distance of measured SFR from the ridgeline of intrinsic MS

[TABLE]

In this work we are also assuming that the time and mass dependence of $\log\mathrm{SFR}_{\rm MS}(M_{\star},t)$ can be neglected, i.e., $\log\mathrm{SFR}_{\rm MS}(M_{\star},t)\equiv C$ , where C is a constant in units of $M_{\odot}~{}\mathrm{yr}^{-1}$ . This is done to simplify the calculation and is indeed true at low redshifts (see further discussion in Section 5.2.2).

By combining Equations (7), (8) and (9) we find that

[TABLE]

The above expression is true for a single galaxy, but in observations we are looking for a statistical quantity $\langle\Delta_{\rm ind}(t_{i})\rangle$ , where we use $\langle\rangle$ to denote the average value of the measured offset from MS, i.e., to denote that the measurement has been averaged over a sample of observed galaxies.

$\langle\Delta_{\rm ind}(t_{i})\rangle$ can be explicitly derived for the case of damped random walk (equations (2) and (3)). When the galaxy is at a MS distance of $\Delta(t_{j})$ at some time $t_{j}$ , its expected distribution at later time $t_{i}-t_{j}$ is given by (Kelly et al., 2009; MacLeod et al., 2010):

[TABLE]

where with $P$ we denote expected distribution and $\mathcal{N}(\mu,s)$ denotes a normal distribution with mean $\mu$ and standard deviation of $s$ . We are actually interested in the quantity that appears in Equation (10), i.e., the expectation value of the quantity $10^{\Delta(t_{i}-t_{j})}$ . The expectation value, $E$ , for a lognormal distribution is given by

[TABLE]

Therefore, the expectation value of the process at time $t_{i}-t_{j}$ is given by

[TABLE]

Combining together the expressions from equations (10) and (13) we obtain:

[TABLE]

We see that the measured offset depends on the duration of the response of an indicator ( $N_{j}$ ), the properties of underlying stochastic process (only $\tau_{\rm decor}$ in case of the damped random walk) and the intrinsic width of the main sequence ( $\sigma_{\rm MS}$ ). We can also see that larger offsets are realized for larger values of $N_{j}$ , $\tau_{\rm decor}$ and $\sigma_{\rm MS}$ . We will investigate the interplay of these parameters more closely in Section 3.2.

3.1.2 Width of the MS

We now turn to the width of the MS and explore this property as a function of parameters of the stochastic process. The derivation is somewhat involved; see Appendix B for details. There we show that the lowest order expansion, as a function of $\sigma_{\rm MS}$ , of the measured variance of a sample of galaxies, $\sigma^{2}_{\rm ind,MS}$ , is given by

[TABLE]

where $\sigma_{\rm MS}$ is the intrinsic width of the MS. This double summation can be simplified, by using symmetry properties of the auto-correlation function, to yield the measured $\sigma_{\rm ind,MS}$ (Priestley, 1981; Ivezić, Ž., Connolly, A. J., Vanderplas, J. T. and Gray, A., 2014):

[TABLE]

We note that this equation has the same form as the equation for the uncertainty of the mean for the dataset consisting of $N_{j}$ measurements with a homoscedastic measurement error of $\sigma_{\rm MS}$ . This is another intuitive way to think about this quantity that we are measuring. In essence, we are exploiting the ergodic property of our sample. Measuring the width of the distribution created by galaxies which are moving stochastically around the MS with an indicator over $N_{j}$ times is equivalent to measuring the mean of a sample of galaxies on MS at $N_{j}$ different times!

We see that in case of $ACF(t_{j})=0$ , the expression in the brackets of Equation (16) is unity, and therefore the measured scatter of the MS drops very quickly as $\sigma_{\rm MS}/\sqrt{N_{j}}$ , which is the familiar expression for the uncertainty of the mean from any analysis that uses independent measurements. On the other hand, in case of highly correlated data ( $ACF(t_{j})\approx 1$ ), the measured scatter of the MS is essentially unchanged and corresponds to the initial $\sigma_{\rm MS}$ .

3.2 Comparison to numerically generated model galaxies

We now consider the consequences of the theory developed above, i.e., we study the interplay between the stochastic process of star-formation and the timescale of the indicators used to probe the SFR. We show in Figure 4, as a function of PSD parameters, how the shape of the MS (i.e., distribution of $\Delta_{\rm MS}$ at fixed mass) depends on the different averaging timescale. The grey line, which is the same in all panels, shows the intrinsic distribution of the 1000 model galaxies about the MS ridgeline ( $\Delta_{\rm MS}$ ). By construction, this intrinsic distribution is defined in 1 Myr time bins and is a normal distribution with $\sigma_{\rm MS}=0.4~{}\mathrm{dex}$ .

Let us first consider the effect of averaging the SFH, i.e., mimicking an observation by averaging the SFR over a certain time interval $t$ (which we denote with $\langle t\rangle$ ). Note that this is fundamentally different than considering an average of SFH relative to the MS, i.e., considering an average of $\Delta_{\rm MS}(t)$ (see Equation (7)). The blue, orange, green, and red histograms show the distribution about the MS ridgeline after averaging over 5 Myr, 10 Myr, 100 Myr, and 1000 Myr, respectively. As expected, we find that measurements over the shorter averaging timescales give measured distributions that are closer to the intrinsic distribution. On the other hand, long averaging timescales produce a narrower distribution with a median value that is biased to higher values than in the original distribution. This is a consequence of arithmetic averaging of logarithmic quantity (see Equation (14)) and serves as a useful reminder that the measured median of the MS is not necessarily at the same value as the intrinsic median of the MS.

We now consider different processes measured with the same indicator, i.e., consider the differences between histograms in different panels with the same colour. We find that the observed distribution matches the intrinsic distribution better for larger values of $\alpha$ and $\tau_{\rm break}$ . This is because the SFR of individual galaxies changes less during the time covered by an indicator and is therefore more similar to its momentary value.

Therefore, we see that there is an interplay between the parameters of a stochastic process ( $\alpha$ and $\tau_{\rm break}$ ) and the averaging timescale (timescale of the observational indicator). In particular, for a given slope $\alpha$ , there is a strong connection between between the timescale of the process and the averaging timescale. This can be seen in Figure 4 as a similarity of the MS distributions that have the same $\tau_{\rm break}$ and averaging timescale in a given row of the figure (compare orange histograms in the panels on the left hand side, green in the middle and red on the right hand side).

We note that our formalism is able to predict further moments of the distribution, such as skewness and kurtosis. However, given that those increasingly depend on the details on the assumed intrinsic distribution and are less readily observed quantities, we choose not to extend that part of our analysis.

We present the information about offset ( $\Delta_{\rm MS,\langle t\rangle}$ ) and width ( $\sigma_{\rm MS,\langle t\rangle})$ of the MS contained in Figure 4 in a more condensed form in Figure 5. In the panels on the left we show changes of the measured width of the MS as a function of (a) a fixed slope $\alpha$ and changing $\tau_{\rm break}$ , and (b) a fixed $\tau_{\rm break}$ and changing slope $\alpha$ . On the right hands side we show the result of the same exercise on the offset of the measured MS ridgeline from the intrinsic MS ridgeline.

These plots confirm our intuitive findings from above. For a given indicator corresponding to a single averaging timescale $\langle t\rangle$ , the larger $\alpha$ and shorter $\tau_{\rm break}$ leads to measuring a tighter MS and a larger offset from the median of the intrinsic MS.

Until now we have discussed how measurements (averaging in this simple toy case) and parameters of stochastic processes influence the statistical parameters of the final distribution such as median and standard deviation. In Figure 6 we show a full example of the effect of averaging the SFR (i.e., measuring the SFR with a given indicator) on the observed MS. For all three panels we assume an averaging timescale of 100 Myr and the slope $\alpha=2$ , while varying the $\tau_{\rm break}$ of the process. The blue points show the individual model galaxies. The dashed red line indicates the predicted connection between the measured and actual (intrinsic) offset from the MS according to Equation (14). This prediction is in excellent agreement with the relation measured from simulated model galaxies, which is shown in orange. In the bottom right corner of each panel we also show the measured scatter of the MS and the predicted value from Equation (16).

In the left panel of Figure 6, we can now explicitly see the central limit theorem in action: because the SFRs vary rapidly, the measured and intrinsic SFRs are almost completely unrelated. This is a consequence of the fact that the observational timescale is much longer than the decorrelation time of the process. This means that during the timescale of the observation, each galaxy has moved up and down the MS ridgeline and hence the measured SFR gives us only an estimate of the mean of the whole population. In the middle and right panels of Figure 6, we increase $\tau_{\rm break}$ : the decorrelation times become longer and the agreement between the intrinsic and measured SFRs improve as galaxies change their SFRs less significantly during the time that measurement is sensitive to.

3.3 Measuring parameters of stochastic process from the mock data

In the previous section we have extensively discussed how averaging, i.e., measuring, SFRs changes the inferred properties of the MS, such as position of ridgeline and width. Furthermore, we have shown how this effect depends on the parameters of the stochastic process itself. Here we want to show how we can use these insights to deduce the parameters of the stochastic process from mock observations. We will then apply the same procedure, with more realistic response functions of different indicators, when discussing real data in Section 5.

Here we consider an example in which we have simulated 10000 galaxies with $\alpha=2$ and $\tau_{\rm break}=400$ Myr. We have then convolved these galaxies with three step functions of 1, 10 and 100 Myr respectively, mimicking 3 different indicators. Given that we 1 Myr is the shortest time-scale that we consider, “measurment” with the short indicator is the equivalent to the intrinsic SFH. We then measure the ratio of the widths of MS measured in each of these indicators, with associated error bars, and then find parts of the parameter space that are consistent with these measurements, both from considering Equation (16) and from comparing with results from a simulation which spans the whole parameters space. We consider ratios of the measurements in different indicators to minimise the dependence of the result on the intrinsic width of the MS, which is of course known in the simulation, but is not necessarily available when considering real data. We see from Equation (16), which is first order approximation to the actual measured width, that considering ratios of measured widths eliminates dependence on intrinsic $\sigma_{\rm MS}$ . Some dependence comes from higher order terms to Equation (16), which we discuss in Appendix B.

We show the result of this procedure in Figure 7. The dashed lines show the result deduced from analytic expression in Equation (16) while the differently coloured areas show 1- $\sigma$ constraints from the numerical simulations. We see that, although a single ratio of MS width measurements of two indicators is not sufficient to break the degeneracy between $\alpha$ and $\tau_{\rm break}$ , combination of several indicators can break this degeneracy and provide a unique solution for the parameters of the process. This is not surprising, and is a restatement of the fact that we need two different “measurements” in order to uniquely determine two free parameters. It is also advantageous if the indicators cover substantially different range of timescales, in order to be sensitive to different parts of the PSD. Given that we aim to be as agnostic as possible to the intrinsic width of the MS and therefore consider ratios of the MS width, this effectively requires a measurement of the width of the MS with three indicators. Addition of any more indicators or combinations of indicators will not bring more information in this plane, but can be used as a consistency check.

Furthermore, we see that there is a small offset from the input values and the values deduced from using Equation (16). This is a consequence of the fact that Equation (16) is only the first order approximation to the full calculation of the width, as already discussed in Section 3.1.2 and Appendix B. We use full numerical simulations to avoid this effect here and when we compare with the real data in Section 5.

4 Finding parameters of stochastic star-formation history with realistic response functions

We highlighted in the previous section that the MS scatter for the SFR average over different timescales can give a tight constraint on the PSD of the underlying stochastic process. In this section, we present a detailed investigation quantifying how does the normalisation and width of the MS change when probed with different SFR indicators, which themselves trace different timescales. In Section 5, we will use this information to show example of constraining the PSD from observations.

4.1 SFR tracers

There are multiple methods for determining SFRs from observational tracers such as from the UV, IR and emission lines. In the end, all the methods – irrespective of whether the observable is directly starlight, dust-reradiated starlight or nebular line emission – quantitatively link the observable to the intrinsic emission from hot young stars. Therefore, a crucial part of each method deals with how to correct for dust attenuation. Similarly, the IMF, in particular the high-mass end slope, and the star-formation history play a role in the exact luminosity-SFR calibration.

In the scope of this work, we are primarily interested in the distribution of SFRs probed by different tracers and hence timescales. We assume that measurements of dust-corrected SFRs as provided in the literature are accurately dust corrected, i.e., we do not discuss any dust effects. Similarly, we assume a Chabrier (2003) IMF and solar metallicity throughout. Variations in metallicity, dust attenuation (in particular differential dust effect between young and old stars) and IMF could be, in principle, important. We will address this in an upcoming publication, where we will study the effect of variability of SFHs, IMF and dust attenuation in the observational plane of luminosities. In general, these effects would increase the scatter, implying that our results are basically an upper limit on the variability. In Section 5, we focus on a galaxy sample in a small mass range ( $M_{\star}=10^{10}-10^{10.5}~{}M_{\odot}$ ) in order to minimize variations in IMF, metallicity and dust content.

We are considering the following standard SFR tracers (see Kennicutt, 1998; Kennicutt & Evans, 2012; Davies et al., 2016):

•

H $\alpha$ : H $\alpha$ photons arise from the gas ionized by young ( $<20$ Myr) stars. Thus, H $\alpha$ provides a measure of the current SFR in galaxies.

•

FUV and NUV: UV continuum emission arises from hot, massive O and B stars with $M_{\star}>3~{}\mathrm{M}_{\odot}$ , and hence is a good tracer of more recent star formation, though on longer timescales than H $\alpha$ .

•

$u$ -band: The rest-frame $u$ -band ( $\sim 3500$ Å) emission arises from the photospheres of young, massive stars, tracing the star formation on $\la 100$ Myr timescales. This emission is less strongly affected by dust obscuration than the UV, but the possible large contribution from older stellar populations makes it more difficult to interpret.

•

WISE W3: Probing infrared (IR) flux from star-forming galaxies in the $3-100~{}\mu\mathrm{m}$ range gives a reliable estimate for the ongoing star-formation: the amount of flux emitted in the IR is directly related to the UV emission from newly formed stars (Calzetti et al., 2007, e.g.). We consider the Wide-field Infrared Survey Explorer (WISE) passband W3 at rest-frame $12~{}\mu\mathrm{m}$ to estimate the SFR (Cluver et al., 2014).

•

UV+IR: As discussed above, UV emission arises directly from star-forming regions, while some fraction of this emission is absorbed and reprocessed by dust, being re-emitted in the IR. We therefore sum both UV and total IR luminosities to obtain the total SFR, based on the bolometric luminosity of OB stars.

In order to now predict what the SFH, as seen by a certain indicator looks like, we need to convolve the intrinsic SFH with the response function of that indicator. Response functions for the aforementioned indicators are shown and discussed in Appendix C. Briefly, we predict the luminosity evolution for a Simple Stellar Population (SSP) using the Flexible Stellar Population Synthesis code (FSPS333https://github.com/cconroy20/fsps; Conroy et al., 2009). FSPS has been extensively calibrated against a suite of observational data (for details see Conroy & Gunn, 2010). Throughout this work, we adopt the MILES stellar library and the MIST isochrones. Since the W3 and UV+IR SFR tracers are based on emission from a large wavelength range, these indicators depends on the dust attenuation prescription. We assume a uniform screen with Cardelli et al. (1989) attenuation curve and an opacity of $\tau=0.5$ at $5500~{}$ Å, motivated by findings by Salim et al. (2018) in the local universe. Changing this prescription does not affect our findings here significantly. As expected, and shown in Appendix C, the H $\alpha$ SFR indicator shows a rapid response function, followed by the FUV SFR indicator and the other indicators.

Figure 8 shows the effect on the SFH by convolving the SFH with the response functions for the different indicators. In particular, the upper panel shows the distance from the MS, while the lower panel shows the absolute SFR under the assumption that the MS has a SFR of $1~{}M_{\odot}~{}\mathrm{yr}^{-1}$ at all times and masses. This SFH has PSD described with a high frequency power-law slope of $\alpha=2$ and a decorrelation timescale of $\tau_{\rm break}=100~{}\mathrm{Myr}$ .

The SFR inferred from H $\alpha$ traces the intrinsic SFH well, only missing fluctuations on the shortest timescale ( $\sim 4$ Myr). All other tracers, which probe longer timescales, show a much larger differences compared to the intrinsic SFH. We find that the SFHs are not only smoother, but because of the long-timescale tail in the response function, the SFRs are biased high. There is an overall mean offset of the MS. Additionally, we find that all the variability on short timescales ( $<20$ Myr) are washed out, giving rise to a narrower MS width. This is overall consistent with the simple analysis presented in the previous section.

4.2 Shape of the MS as seen by different SFR indicators: normalization and scatter

We go now one step further and look how the effect on the MS depends on the PSD itself. Figure 9 is a remake of Figure 4, showing the distribution of the distances from the MS ridgeline for different SFR indicators. We see the same qualitative behavior that we discussed extensively in Figure 4, when we discussed the simplified case in which we simply averaged SFRs over a certain timescale. Measurements in H $\alpha$ , which probe star-formation on the shortest timescale, follows the intrinsic distribution closely for all but the most “bursty” set of parameters. On the other hand, all other indicators show some level of deviation, and produce a narrower distribution that is offset from the intrinsic value of the MS. Indicators that act on longer time-scales tend do produce narrower distributions, which is in very good aggrement with observations (e.g., Figure 7 from Davies et al. 2019). The differences in the shapes, even if not so pronounced as in the simple case shown in Figure 4, can be again used to infer the intrinsic properties of the underlying stochastic process.

Although both the normalisation as well as the width of the MS as seen by different indicators can in principle be used to derive the PSD, we are focusing now on the width of the MS. We do this because the overall normalisation depends on the exact conversion used to derive SFRs from fluxes, which are uncertain and depend on several assumption like the IMF, metallicity, as well as the assumed SFH. On the other hand, the width of the MS is a relative measure that depends only on the rank-ordering in the galaxy sample and as such it is largely independent on the exact flux-SFR conversion.

Figure 10 shows how the scatter of the MS changes when traced by different SFR indicators. This figure assumes that the intrinsic width of the MS is 0.4 dex (probed on a 1 Myr timescale). Specifically, each panel shows the MS scatter as a function of the PSD parameters $\alpha$ and $\tau_{\rm break}$ . We see a similar pattern that we have also seen when we have discussed “burstiness" (Figure 3). Consistent with the previous figure, measurements in H $\alpha$ follow the intrinsic width closely for large parts of the parameter space. Measurements in other indicators that depend more heavily on the previous SFH provide a measurement which is less readily interpretable. In general, as discussed in Section 2, processes with low $\alpha$ and $\tau_{\rm break}$ that oscillate wildly, lead to tighter observed MS, courtesy of the central limit theorem. Although measured width in different indicators are relatively similar, the differences between them enable us to constraint parameters of the underlying process in Section 5.

4.3 Ratio of SFRs with different indicators

We now move forward to a MS-independent discussion and directly compare different SFR tracers with each other in Figure 11. Specifically, we plot the ratio of FUV and H $\alpha$ SFRs versus the ratio of $u$ -band and H $\alpha$ SFRs. The color-coding corresponds to the ratio of the SFR measured over 1 Myr and the SFR measured over 100 Myr. There are three trends visible from this figure. Firstly, the tightness of the relation increases with increasing burstiness. This arises because all SFR indicators probe the same SFR if the SFR varies on shorter timescales than the timescales of the SFR indicators. Secondly, the slope decreases with decreasing burstiness. This can be explained by the fact that the H $\alpha$ and FUV SFR indicators are able to probe short-term variability better than the $u$ -band SFR indicator, which, in the case of slowly varying SFHs, leads to a FUV-to-H $\alpha$ ratio close to 1 while the $u$ -band-to-H $\alpha$ ratio shows significantly more scatter. Thirdly, the color-coding shows that if the SFH quickly increases, H $\alpha$ is able to probe this, leading to ratios smaller than 1, while if the SFH quickly decreases, the ratios are all larger than 1. Overall, this figure highlights that – in the framework presented in this work – there are quantifiable consequences which arise when simply investigating the relation between the SFR measured with different indicators.

5 Observational constraints

In this section, we wish to apply our insights about different stochastic processes and their influence on observed MS properties. The main goal is to use the MS width measurements to determine parameters of the stochastic process, i.e., slope $\alpha$ and decorrelation timescale $\tau_{\rm break}$ of the PSD, that govern the evolution of star-forming galaxies about the MS ridgeline. The analysis in this section should be regarded only as a proof of concept as we work with the simplest possible assumptions, in order to highlight the power of our approach. We show that the result for the parameters that we derive are consistent when using different indicators, which is highly encouraging.

5.1 Main result

The ideal dataset should contain statistically significant sample of galaxies and be covered in a several observational bands, in order to check consistency between different SFR indicators. Additionally, the analysis is simplified in the case of a galaxy sample at low redshift. This is because of (a) the slow evolution of the MS as a function of time, which simplifies interpretation of measured SFRs as an offset from the virtually non-changing MS; and (b) the slow growth of mass, meaning that measurement of total stellar mass and star-formation are independent quantities. We discuss these issues more in Section 5.2. Additionally, low redshift sample also avoids potential problems and systematics that might arise due to K-correction.

We find the measurement of the MS width presented in Davies et al. (2019) to be well suited for this analysis. Davies et al. (2019) used data from Galaxy and Mass Assembly Survey (GAMA; Driver et al., 2011; Liske et al., 2015; Driver et al., 2016; Baldry et al., 2018) and conducted an analysis which was explicitly measuring width of the MS with different SFR indicators. The sample consists of 9005 galaxies and is restricted to $z<0.1$ , satisfying our requirement to work with a large sample at low redshift. It includes only galaxies which are not classified as being a member of a group or a pair in the GAMA group catalogue, as described in Robotham & Driver (2011). This means that our final result will be more suited to describing star-formation fluctuations due to “secular” evolution instead of mergers.

We consider the following measurements: (a) extinction-corrected H $\alpha$ SFRs, measured with the process outlined in Gunawardhana et al. (2011); Hopkins et al. (2013); Gunawardhana et al. (2015) and using line measurements from Gordon et al. (2017), (b) UV+IR SFRs, derived from the galaxy spectra as described in Brown et al. (2014), and (c) extinction-corrected $u$ -band SFRs from rest-frame u-band luminosity Davies et al. (2016). SFRs have been derived using Chabrier IMF. We refer the reader to Davies et al. (2019) for more details.

We consider MS width measurements for SFR selected galaxies between $10^{10.0}$ and $10^{10.5}~{}M_{\odot}$ as we expect that the measurement errors are the smallest for this subset of galaxies. Furthermore, more massive galaxies are transitioning from the star-forming to the quiescent population (quenching, e.g., Peng et al., 2010), which could lead to an artificially large scatter of the MS, which itself would depends heavily on the selection of star-forming galaxies. We estimate, using Equation (27) from Peng et al. (2010) with reasonable quenching time of 1 Gyr, that only around 2 $\%$ of star-forming galaxies at this mass are being “quenched” at redshift 0. Additionally, considering the rather narrow mass range makes our result less affected by possible variations of dust and IMF within the sample. In order to then derive the observational constraints, we assume a relative error on the ratios of $7\%$ to include possible observational uncertainties in the error budget (Davies et al., 2019).

The final result is shown in Figure 12. As elaborated in Section 3, we consider ratios of the MS widths measured with different indicators to reduce the influence of the unknown intrinsic width of the MS on our measurement. We see the similar structure as in Figure 7, with the general degeneracy between $\alpha$ and $\tau_{\rm break}$ . We wish to highlight the consistency between the parameter space covered by all three combination of the indicators used, albeit the ratio $\sigma_{\rm u,MS}/\sigma_{\rm UV+IR,MS}$ is quite uninformative as it covers large part of the considered parameters space. Nevertheless, we find the consistency between all of the indicators very encouraging, especially given the simplicity of our assumptions and minimal consideration of the dust and IMF effects.

Unfortunately, with this data and simplicity of our analysis is not possible to resolve the degeneracy between $\alpha$ and $\tau_{\rm break}$ , as it has been possible in the simplified case considered in Section 3.2. This is due to combination of two factors. Firstly, this is due to larger uncertainties in the observational measurement of the $\sigma_{\rm MS}$ , driven by the smaller number of galaxies in the sample and by the observational uncertainties. Secondly, fact that realistic indicators are not a simple step functions but have non-trivial time dependence that mixes many timescales reduces our ability to differentiate between the influence of $\alpha$ and $\tau_{\rm break}$ in observations (see Figure 10 and Appendix C).

In order to quantify our final result in one number, we report $\tau_{\rm break}$ after setting the high-frequency slope at $\alpha=2$ , i.e., for the case of damped random walk. We emphasize however that the data itself does not prefer $\alpha=2$ over any other particular value for $\alpha$ and our choice for reporting $\tau_{\rm break}$ for damped random walk is driven by theoretical and empirical reasons outside of the scope of this paper. The first reason for considering the damped random walk is its theoretical and physical simplicity. Secondly, the slope of $\alpha=2$ is seen when studying time-variability properties of galaxies in simulations (Iyer et al. in prep.). In the case of $\alpha=2$ , we find that the region of the $\tau_{\rm break}$ parameter space that is within 1- $\sigma$ contours for all combinations of indicators is given with $\tau_{\rm break}=178^{+104}_{-66}$ Myr. We intentionally do not conduct further statistical consideration of the allowed region, given that we believe that systematic uncertainties, due to variation of dust, IMF and observational uncertainties dominate our error budget and would make any such analysis highly uncertain. We do note that measurement of the width of MS by other indicators available in Davies et al. (2019) are in broad agreement with the estimated value that is reported above.

5.2 Discussion

5.2.1 Mass dependence

The result of $\tau_{\rm break}=178^{+104}_{-66}~{}\mathrm{Myr}$ that we have derived is only valid for the sample that we considered, consisting of isolated galaxies with $\log(M_{\star}/M_{\odot})=10.0-10.5$ . As such, the derived timescale is not necessarily the same for the whole range of masses and environments in the Universe. We do note that there is a mass dependence in the measurements of the width of the MS in Davies et al. (2019). We do not conduct the mass dependence analysis in this work given the possible mass dependence in the dust and IMF properties. We do however note that the measured widths of the MS suggest similar or somewhat shorter timescale for $\tau_{\rm break}$ for the smaller masses within the dataset ( $8<\log(M/M_{\odot})<10$ ), with some dependence on the exact set of indicators used.

5.2.2 Caveats

Of course, our analysis comes with several caveats. One complication that we consider is the intrinsic correlation between the stellar mass and measured SFR. As we have elaborated, the measured SFR can be considered as an imperfect averaging over the intrinsic SFR, where the length of this averaging is dependent on the indicator used for the measurement. On the other hand, the stellar mass of a galaxy is also obviously related with the previous SFR, given by the integral over the complete SFH of the galaxy. Hence, the measurement of the MS width is only meaningful if the total stellar mass added over the period that a given indicator is measuring is small compared to the total stellar mass of a galaxy. Given that the mass doubling time for star-forming galaxies today is very long, i.e., comparable to Hubble time, the mass added over the time effectively probed by the different indicators is negligible (see also Appendix D). This effect, however, needs to be accounted for at higher redshifts (e.g., $z>1$ ) where the timescale of the measurement of star-formation in slow-response indicators starts being a non-negligible fraction of the lifetime of the Universe (e.g., Tacchella et al., 2018).

Similarly, when conducting our analysis of observational results we have assumed that the overall change of the mean SFR of the MS is negligible. Consider the timescale of 400 Myr, which is more than twice as long as the $\tau_{\rm break}$ derived above. We show in Appendix C that by this timescale all of the considered indicators have outputted more than 60% of their complete luminosity, i.e., are very weakly sensitive to the star-formation on longer timescales. If we assume relatively steep evolution of the sSFR of the MS, such as $(1+z)^{3}$ (e.g., Lilly et al., 2013, and references therein), we find that the ratio between $\mathrm{sSFR}_{\rm MS}(400\mbox{ Myr})/\mathrm{sSFR}_{\rm MS}(\mbox{today})\approx 1.1$ . As such this change will have minimal effect on the inferred result in this study, but it needs to be considered at higher redshifts.

Finally, observational uncertainties introduce scatter and, hence, increase the MS width obtained by Davies et al. (2019). These uncertainties arise, beside other things, from the aperture correction, the dust attenuation correction, and the assumption that all galaxies have the same IMF and metallicity. It is beyond the scope of this paper to assess the impact of those uncertainties in detail, but these all contribute to relatively large error on measured $\sigma_{\rm ind,MS}$ that we have used.

6 Summary and conclusions

The knowledge of how the SFR changes as a function of time for individual galaxies is a critical ingredient needed for understanding the evolution of galaxies. SFHs encode information about the variability on short and long timescales that arise from different physical processes, such as gas accretion, mergers, and feedback from stars and black holes. The best method to model, and mathematical describe, SFHs is is still being debated.

In this paper, we present a framework for modelling SFHs of galaxies as a stochastic process. We focus on star-forming galaxies and assume that a stochastic process describes how these galaxies move about the MS ridgeline. We also assume that at any moment the intrinsic distribution of SFRs around the MS ridgeline is log-normal. We define this stochastic process through a power spectrum density. The functional form of the PSD is given by a broken power-law, with high-frequency slope $\alpha$ that flattens out at a frequency corresponding to $\tau_{\rm break}$ . The physical meaning of the PSD is that SFHs are correlated on short timescales, where the strength of this correlation is described by the slope $\alpha$ , and they decorrelate to resemble white noise at a timescale that is proportional to $\tau_{\rm break}$ . Specifically, when the slope of the PSD is shallow and/or $\tau_{\rm break}$ is small, the SFH violently oscillates around the mean and is very quickly decorrelated. On the other hand, steeper slopes and longer $\tau_{\rm break}$ lead to slower oscillations and a more correlated behavior.

We demonstrate that, even though we can not observationaly follow the time-evolution of individual galaxies, we can deduce the properties of the stochastic process. This can be done by measuring properties of the MS, such as the MS normalisation and the MS scatter, with different observational indicators that are sensitive to different timescales. There exists a degree of degeneracy between the influence of the parameters $\alpha$ and $\tau_{\rm break}$ , that can be roughly described by the “burstiness”, i.e., by the fact that galaxies described with smaller $\alpha$ and shorter $\tau_{\rm break}$ tend to more often have SFRs which are outside of the 1- $\sigma$ width of the MS. We show that this degeneracy can be broken, in principle, by measuring the properties of the MS in several observational indicators.

We then consider realistic observational SFR indicators, and show results for the measured properties of the MS in H $\alpha$ , FUV, NUV, $u$ -, W3 and UV+IR bands. We use those results to, as a proof of concept, deduce parameters of the underlying stochastic process from the data from the GAMA survey, presented in Davies et al. (2019). For a sample of isolated galaxies with $\log(M_{\star}/M_{\odot})=10.0-10.5$ and using MS widths measured from H $\alpha$ , $u$ - and UV+IR bands, we find the area of the parameter space that is consistent with this data. With an assumption that the high-frequency slope is given with $\alpha=2$ , we obtain $\tau_{\rm break}=178^{+104}_{-66}$ Myr. The motivation for $\alpha=2$ stems from numerical simulations as well as the theoretical attractiveness of the damped random walk process. The quoted errobars are deduced from the observational 1- $\sigma$ errors of all the available indicators.

Our result of $\tau_{\rm break}\approx 200$ Myr indicates that the SFHs of galaxies decorrelate, i.e., lose memory of their previous SFH, on roughly this timescale. This is shorter than the dynamical time of the dark matter halo, which is roughly 10% of the Hubble time (i.e., a bit more than a Gyr at $z\sim 0$ ). Therefore, we conclude that baryonic effects within the halos, which act on dynamical timescales of galaxies, play an important role in shaping the SFHs of galaxies. Examples of such effects could include feedback and reincorporation of galactic winds.

In the future, we will apply this framework to a wide range of numerical and semi-analytical models, where we can measure the PSD directly of the models’ SFHs. Preliminary results show that the obtained PSDs are well described by a broken power-law with nearly all theoretical models producing $\alpha\approx 2$ , while having different values for $\tau_{\rm break}$ . We show that the PSD is set by physical ingredients in the models and hence that the PSD contains a wealth of information that has not yet been researched. On the observational side, we will address many issues that complicate precise determination of the PSD in observations, including possible variations in the IMF and dust attenuation. There we also wish to use these insights to analyse properties of galaxies at higher redshifts, where it is especially important to consider these complications.

Acknowledgements

We thank Luke Davis for sharing his measurements from the GAMA survey with us. Furthermore, we are thankful to Charlie Conroy, Luke Davis, Daniel Eisenstein, Kartheik Iyer, Ben Johnson, Daniel Kelson, Joel Leja, and Michael Strauss for very useful comments. We thank Sophie Reed and Hassan Siddiqui for proofreading the manuscript. This research made use of NASA’s Astrophysics Data System (ADS), the arXiv.org preprint server, the Python plotting library matplotlib (Hunter, 2007), astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013), and the python binding of FSPS (Foreman-Mackey et al., 2014). S.T. is supported by the Smithsonian Astrophysical Observatory through the CfA Fellowship.

Appendix A Fractional Brownian motion and Hurst parameter

In this appendix, we shortly discuss the connection between the description of the stochastic processes that we have used in this work and the fractional Brownian motion described with the Hurst parameter, as used in Kelson (2014)444Here we discuss “type 2” fractional Brownian motion, which is more suited for descriptions of physical process as in that description the future events are only dependent on the previous events. “Type 1” processes are dependent on both future and past events, violating causality.. Fractional Brownian motion is process described with the ACF (Marinucci & Robinson, 1999):

[TABLE]

where

[TABLE]

and $H$ is called the “Hurst exponent” or “Hurst parameter” and it is bound between 0 and 1. In the limit when $t\rightarrow\infty$ the relation simplifies to (Marinucci & Robinson, 1999):

[TABLE]

Tarnopolski (2016) shows that when $H>1/2$ , which describes a physically interesting case where observations are more correlated at shorter timescales, the autocorrelation function behaves as $ACF(t)\propto t^{2-2H}$ . The fact that the ACF decreases slower than $t^{-1}$ is often expressed as a statement that the process has a “long-range” dependence. One can compare this to ACFs that we present in right panel of Figure 2 (see also Equation (4)), which all reach $ACF(t)\sim 0$ at timescales correlated with $\tau_{\rm break}$ . In other words, by modelling the PSD as broken power law, where the small-frequency slope is zero (mimicking white noise), we have explicitly allowed for the system to decorrelate.

Given that the processes that we use in this work have an explicit timescale built into them through the $\tau_{\rm break}$ parameter, we see that any connection that we want to derive between these two different approaches (fractional Brownian motion and broken-law PSD) will fundamentally depend on the timescale over which the connection is measured. In any case, it is crucial to note that any such a connection is essentially only an approximation, given the different definitions of both processes. In Figure 13 we show the Hurst parameter which best approximates the SFHs generated with different $\alpha$ and $\tau_{\rm break}$ . The SFHs used here have a length of 4.2 Gyr as described in Section 2.3. We see similar dependence as in when we have discussed “burstiness” and widths of the MS measured in different indicators. As expected, a larger Hurst parameter is found for SFHs that are well correlated (larger $\tau_{\rm break}$ and $\alpha$ ) while less correlated, “bursty” SFHs (smaller $\tau_{\rm break}$ and $\alpha$ ) lead to lower Hurst parameters.

Appendix B Connection between star-formation measurements in different indicators

In this appendix we extend our analytical calculation of the measured width of MS, as seen after measurement in different indicators. We start by noticing that the star formation rate measured by a given indicator is connected to the previous star formation and is given by the intensity of that star-formation episode and a response of that particular indicator to the star formation, i.e.,

[TABLE]

where $\mathrm{SFR}_{\rm ind}$ is the measured star-formation, $\mathrm{SFR}$ is the actual star-formation, and $f_{\rm ind}$ are response functions (the ones that we have used in this work are presented in Appendix C). In the case of continuous star-formation the measured star-formation is given by

[TABLE]

where we have convolved the total star-formation history with an appropriate response function up to the time of the measurement. When discretized the above integral is written as a sum

[TABLE]

where $SFR_{\rm ind}(t_{i})$ is the measured star-formation in bin $t_{i}$ and $w_{\rm ind}(t_{j})$ is the discretized version of the response function:

[TABLE]

where $\delta$ is the duration of the discretization bin. In this work we use $\delta$ = 1 Myr (see Section 2 for justification of this time-scale).

Now, consider the ensemble of galaxies moving around Main-Sequence in stochastic fashion with a log-normal distribution function. Lets denote with $\Delta$ 555As mentioned in the main body of the text, we use just $\Delta$ as shorthand for $\Delta_{\rm MS}$ in order to reduce the size and complexity of our expressions. $\log$ -distance from the median value of MS i.e., (see Equations (8) and (9)):

[TABLE]

We rewrite the Equation (22) in terms of $\Delta$ :

[TABLE]

We note that that factors in the exponent $\Delta(i-j)$ are in general expected to be small as they are bound by the width of the MS. In that case we can use the MacLauren expansion:

[TABLE]

Lets continue by taking first two orders of expansion above in Equation (26):

[TABLE]

where we have used the property that the response function has to sum to unity, i.e., $\sum^{\infty}_{t_{j}=0}w_{\rm ind}(t_{j})=1$ . We now expand the $\log$ function using

[TABLE]

and we find the expression

[TABLE]

Let us now find variance of the measured MS after considering only the first term of the expansion above :

[TABLE]

where we have used the fact that the co-variance of the process is given by $COV(\Delta(t_{i}),\Delta(t_{j})$ = $\sigma^{2}_{\rm MS}ACF(t_{i}-t_{j})$ and where we have used symbol $E$ to denote expectation value. We have also used the definition of the covariance together with the fact that the expectation value of a stationary process equals zero:

[TABLE]

Taking into the account next order terms in Equation (29) will lead to extra terms in the summation (30), with the smallest order being proportional to $\sigma^{4}_{\rm MS}$ . This is because all $\sigma^{3}_{\rm MS}$ terms will inevitably be connected with terms $E(\Delta^{2}(t_{i}))\cdot E(\Delta(t_{j}))$ or $E(\Delta(t_{i}))\cdot E(\Delta^{2}(t_{j}))$ which are equal to zero, given the stationary of our problem. For the realistic widths of MS ( $\sigma_{\rm MS}\sim$ 0.4 $\pm 0.1$ ) these $\sigma^{4}_{\rm MS}$ terms are sub-dominant and therefore lead only to relatively small corrections (e.g., see the relatively small difference between the input and recovered parameters in Figures 6 and 7).

In the simplest case, discussed in Section 3, in which we assume that a given SFR indicator is given by a step function lasting $N_{j}$ time units, i.e.,

[TABLE]

the last expression in Equation (30) can be written as

[TABLE]

and we can now recognize Equation (15) used in the main body of the text.

Appendix C Luminosity evolution of a Simple Stellar Population (SSP)

As we have extensively discussed in the main text of the paper, the ability to differentiate different stochastic processes that govern star-formation is dependent on being able to measure the SFR with different indicators that probe different timescales. To be sensitive to a wider range of different PSDs, these response functions of the SFR indicators would ideally be well differentiated and cover a large range of timescales. We show in Figure 14 the luminosity evolution, seen in different passbands, of a Simple Stellar Population (SSP). These are effective response functions, i.e., the observed SFH in any of these bands can be considered as a convolution of the intrinsic SFH with these luminosity evolutions. In the left panel, we see that all of the traces exhibit peak at early times ( $\sim$ 5 Myr), followed by an exponential-like decline. The right panel shows the cumulative luminosity evolution, indicating more clearly the long term behaviour of the different indicators. H $\alpha$ responds very quickly to star-formation – therefore, the measurement of the MS width in H $\alpha$ is very similar to the intrinsic MS width for almost the whole range of parameter space that we have considered, as seen in Figure 10. As expected, the indicators corresponding to longer wavelengths tend to have slower responses, showing long “memory” of the previous SFH. The fact that all of these indicators are a mixture of the recent and long-term star-formation complicates matters and makes the interpretation of the measurements, such as width of MS, quite unintuitive, but they can be modelled through simulations (Figure 10).

Appendix D Timescales of different SFR tracers

We have shown in the main body of the manuscript how the key relations are simplified and well understood in the case where the SFR response function is a simple step function. On the other hand, in Appendix C we show that the real response functions have non-trivial time dependence and that they are clearly not well described by a step function. Crucially, the timescales that best approximate different SFR tracers actually depend on the property of the SFH itself. In practice, similar to the derivation of the luminosity-SFR calibration, one usually assumes a constant SFR and then determines the luminosity-weighted mean age of the stellar population (e.g., Kennicutt, 1998; Kennicutt & Evans, 2012; Madau & Dickinson, 2014). A more sophisticated approach is to consider more realistic SFHs from numerical and analytical simulations, with the caveat that the variability of the SFHs is simulation-dependent and it might bias the measurement (e.g., Tacchella et al., 2018). Already Johnson et al. (2013) showed – based on empirical SFHs – that fluctuations in the SFHs alone can cause factor of $\sim 2$ variations in the UV luminosities relative to the assumption of a constant SFH over the past 100 Myr. Thanks to our framework of stochastic SFHs introduced in Section 2, we are able to generate and understand SFHs that have different variability properties according to a power spectral density.

Here we are looking to derive a singular timescale over which to average the realistic star-formation response to best approximate complex and realistic functions with a step function of a given duration. We will first briefly express our problem analytically and afterwards derive these timescales numerically for the different indicators that we used in Section 4 (H $\alpha$ , FUV, NUV, u, W3 and UV+IR).

The statement about approximating realistic star-formation responses as a step function can be expressed as a search for the best $N_{\rm approx}$ to approximate Equation (25):

[TABLE]

where

[TABLE]

Comparing the expressions within the logarithms we see that the condition simplifies to the following equation for $N_{\rm approx}$ :

[TABLE]

From here we can see that in order to derive $N_{\rm approx}$ we have to consider both the properties of stochastic process, that govern time dependence of $\Delta(t_{i}-t_{j})$ term, and various indicators with their response function described with $w_{\rm ind}(t_{j})$ . In other words, as stated before, the timescale associated with each SFR traces depend on the previous SFH (in addition to the properties of each indicator).

In order to derive the timescales of the different SFR tracers numerically, we generate SFHs with different $\alpha$ and $\tau_{\rm break}$ , as described in Section 2. We estimate the timescale of the SFR tracer by solving Equation (33) in a statistical sense, by minimizing the RMS difference of the SFR of the tracer and the intrinsic SFR averaged over different timescales. In other words, we search for $N_{\rm approx}$ that minimizes the variance of the quantity $\Delta_{\rm ind}(t_{i})-<\Delta(t_{i})>_{N_{\rm approx}}$ .

In Figure 15, we present the result of this analysis for H $\alpha$ , FUV, NUV, u, W3 and UV+IR. As expected, H $\alpha$ corresponds to shortest timescale while UV+IR and u-band correspond to longer timescales (e.g., see Figure 14). We however note large variations across the parameter space of $\tau$ and $\alpha$ , confirming our intuition about the importance of considering the underlying stochastic process on making a statement about timescales probed by different indicators. The differences are obviously larger for the longer indicators. Intuitively, when the indicator is very quick, as is the case for H $\alpha$ or FUV, effect of “mixing in” of the underlying process is reduced and changes of $\tau_{\rm break}$ and $\alpha$ that affect longer timescale are not probed in this analysis.

For the set of parameters that we have deduced in Section 5, $\alpha=2$ and $\tau_{\rm break}=178^{+104}_{-66}$ Myr, we find that the effective times of the indicators are $t_{\rm eff}(\mbox{H}\alpha)=2.0_{-0.1}^{+0.4}~{}\mathrm{Myr}$ , $t_{\rm eff}(\mbox{FUV})=63_{-1}^{+2}~{}\mathrm{Myr}$ , $t_{\rm eff}(\mbox{NUV})=119_{-2}^{+8}~{}\mathrm{Myr}$ , $t_{\rm eff}(\mbox{u-band})=164_{-16}^{+15}~{}\mathrm{Myr}$ , $t_{\rm eff}(\mbox{W3})=96_{-14}^{+7}~{}\mathrm{Myr}$ , and $t_{\rm eff}(\mbox{UV+IR})=84_{-5}^{+6}~{}\mathrm{Myr}$ . Of course, these estimates suffer from the same caveats that we have discussed in Sections 5 and just therefore only be taken as a general guidance, rather than strict quantitative measurement.

Bibliography121

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Abramson et al. (2015) Abramson L. E., Gladders M. D., Dressler A., Oemler Jr. A., Poggianti B., Vulcani B., 2015, Ap J , 801, L 12 · doi ↗
2Applebaum et al. (2018) Applebaum E., Brooks A. M., Quinn T. R., Christensen C. R., 2018, preprint, ( ar Xiv:1811.00022 )
3Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A , 558, A 33 · doi ↗
4Baldry et al. (2004) Baldry I. K., Glazebrook K., Brinkmann J., Ivezić Ž., Lupton R. H., Nichol R. C., Szalay A. S., 2004, Ap J , 600, 681 · doi ↗
5Baldry et al. (2018) Baldry I. K., et al., 2018, MNRAS , 474, 3875 · doi ↗
6Barro et al. (2017) Barro G., et al., 2017, Ap J , 840, 47 · doi ↗
7Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Conroy C., 2013, Ap J , 770, 57 · doi ↗
8Boogaard et al. (2018) Boogaard L. A., et al., 2018, A&A , 619, A 27 · doi ↗