Log-normal star formation histories in simulated and observed galaxies

TL;DR

This study demonstrates that log-normal functions effectively model the overall shape of galaxy star formation histories in simulations and observations, outperforming traditional models and revealing insights into galaxy evolution.

Contribution

It provides a comprehensive analysis showing log-normal functions fit galaxy SFHs well, compares simulated and observed data, and explores the physical implications of the fitted parameters.

Findings

Log-normal fits reproduce 85% of cumulative SFHs within 5% error.

Log-normal performs better than delayed-tau models for most SFHs.

Simulated and observed SFHs occupy similar parameter spaces, with some differences.

Abstract

Gladders et al. have recently suggested that the star formation histories (SFHs) of individual galaxies are characterized by a log-normal function in time, implying a slow decline rather than rapid quenching. We test their conjecture on theoretical SFHs from the cosmological simulation Illustris and on observationally inferred SFHs. While the log-normal form necessarily ignores short-lived features such as starbursts, it fits the overall shape of the majority of SFHs very well. In particular, 85% of the cumulative SFHs are fitted to within a maximum error of 5% of the total stellar mass formed, and 99% to within 10%. The log-normal performs systematically better than the commonly used delayed-tau model, and is superseded only by functions with more than three free parameters. Poor fits are mostly found in galaxies that were rapidly quenched after becoming satellites. We explore the log-normal parameter space of normalization, peak time, and full width at half maximum, and find that the simulated and observed samples occupy similar regions, though Illustris predicts wider, later-forming SFHs on average. The ensemble of log-normal fits correctly reproduces complex metrics such as the evolution of Illustris galaxies across the star formation main sequence, but overpredicts their quenching timescales. SFHs in Illustris are a diverse population not determined by any one physical property of galaxies, but follow a tight relation, where $\mathrm{width}\propto\mathrm{(peak\ time)}^{3/2}$. We show that such a relation can be explained qualitatively (though not quantitatively) by a close connection between the growth of dark matter halos and their galaxies.

Equations

$$\mathrm{SFR}(\mathrm{t})=\frac{\mathrm{A}}{\sqrt{2\pi {\tau }^2}\times \mathrm{t}}\exp \left(-\frac{(\ln \mathrm{t}-{\mathrm{T}}_0{)}^2}{2{\tau }^2}\right)$$

$$M_{\mathrm{final}}=A\times 10^9\times f_{\mathrm{ret}}$$

$$\mathrm{sSFR}({\mathrm{z}}_{\mathrm{obs}})\le \frac{0.2}{10^9\times t_{\mathrm{obs}}}$$

$$\mathrm{cSFR}(\mathrm{t})\equiv {\int }_0^{\mathrm{t}}\mathrm{SFR}({\mathrm{t}}^{\prime })\mathrm{d}{\mathrm{t}}^{\prime }$$

$$D\equiv \frac{\max \left(\left|\mathrm{cSF}{\mathrm{R}}_{\mathrm{sim}}(\mathrm{t})-\mathrm{cSF}{\mathrm{R}}_{\mathrm{fit}}(\mathrm{t})\right|\right)}{\mathrm{cSFR}(t_0)}$$

$$M(z)=M_0{e}^{-\alpha z}$$

$${\sigma }_{\mathrm{SFR}}=0.83t_{\mathrm{peak}}^{3/2},$$

$$\frac{dM}{dt}\propto e^{-\alpha z}(1+z{)}^{5/2}.$$

$$t_{\mathrm{peak}}=e^{T_0-{\tau }^2},$$

$$\mathrm{SF}{\mathrm{R}}_{\mathrm{peak}}=\frac{\mathrm{A}}{\sqrt{2\pi {\tau }^2}}{\mathrm{e}}^{-{\mathrm{T}}_0+{\tau }^2/2}.$$

$$\mathrm{cSFR}(\mathrm{t})\equiv {\int }_0^{\mathrm{t}}\mathrm{SFR}({\mathrm{t}}^{\prime })\mathrm{d}{\mathrm{t}}^{\prime }=\frac{\mathrm{A}}{2}\left[1-\mathrm{erf}\left(-\frac{\ln \mathrm{t}-{\mathrm{T}}_0}{\tau \sqrt{2}}\right)\right]$$

$$\mathrm{erf}(\mathrm{x})=\frac{1}{\sqrt{\pi }}{\int }_{-x}^x{e}^{-t^2}dt.$$

$$\Delta (\ln t)=2\sqrt{2\ln (\beta )}\tau .$$

$$\Delta t_{\beta }$$

$$=2t_{\mathrm{peak}}\sinh \left(\sqrt{2\ln (\beta )}\tau \right).$$

$${\sigma }_{\mathrm{SFR}}=2t_{\mathrm{peak}}\sinh \left(\sqrt{2\ln (2)}\tau \right).$$

$$\mathrm{SFR}(t)=\frac{A}{{\tau }^2}(t-t_{\mathrm{i}})\exp \left(-\frac{t-t_{\mathrm{i}}}{\tau }\right),$$

$$\mathrm{cSFR}(t)=A\times 10^9\left[\left(-\frac{t-t_{\mathrm{i}}}{\tau }-1\right)\exp \left(-\frac{t-t_{\mathrm{i}}}{\tau }\right)+1\right].$$

$$\mathrm{SFR}(t)=A{\left[{\left(\frac{t}{\tau }\right)}^B+{\left(\frac{t}{\tau }\right)}^{-C}\right]}^{-1}.$$

$$t_{\mathrm{peak}}=\tau {\left(\frac{B}{C}\right)}^{-\frac{1}{B+C}}$$

$${\mathrm{SFR}}_{\mathrm{peak}}=A{\left[{\left(\frac{B}{C}\right)}^{-\frac{B}{B+C}}+{\left(\frac{B}{C}\right)}^{\frac{C}{B+C}}\right]}^{-1}.$$

$$\mathrm{cSFR}(t)=A\times 10^9\times \tau \times f(x,B,C)$$

$$f(x,B,C)\equiv \frac{x^{1-B}\left[(B-1)B{x}^{B+C}{f}_1+C(C+1)f_2-C(C+1)\right]}{(B-1)(C+1)(B+C)}$$

$$\displaystyle f_{1}\equiv\$$

$$\displaystyle f_{2}\equiv\$$

Full text

††journal: The Astrophysical Journal, 839:26 (20pp), 2017 April 10

\SIthousandsep

,

1. Introduction

One of the most fundamental aspects of galaxy formation is the star formation history (SFH), both in individual galaxies and in the universe overall. Unfortunately, the time-resolved SFHs of individual galaxies are difficult to measure observationally in all but the most local galaxies where we have access to resolved stellar populations (Weisz et al., 2011, 2014; Skillman et al., 2014; Williams et al., 2011; Lewis et al., 2015). In more remote systems, we need to rely on stellar archaeology, i.e. measurements of stellar ages based on photometric or spectroscopic observations and stellar population synthesis models (Tinsley, 1968; Gallagher et al., 1984; Sandage, 1986; Kauffmann et al., 2003; Thomas et al., 2005; Kennicutt & Evans, 2012; McDermid et al., 2015; Leja et al., 2017). The number of independent time bins in such SFH measurements is usually small, though better time resolution can be achieved with more flexible parameterizations of the SFH (Tojeiro et al., 2007, 2009; Pacifici et al., 2013, 2016a, 2016b). Instead of focusing on individual galaxies, one can try to measure the star formation rate (SFR) of galaxy populations at different redshifts and connect them to their progenitors statistically, but such inferences are complicated by scatter and merging in a $\Lambda$CDM universe (Behroozi et al., 2013a; Torrey et al., 2015, 2017; Wellons & Torrey, 2017).

In a global sense, however, many fundamental aspects of star formation in the universe are now well established. The star formation rate density (SFRD) peaks around $z\approx 2$ and declines thereafter (Lilly et al., 1996; Madau et al., 1998; Hopkins & Beacom, 2006; Madau & Dickinson, 2014). At any redshift along this global trajectory, star-forming galaxies exhibit a correlation between their stellar mass and SFR, leading to the concepts of a “star formation main sequence” and a quiescent population (Brinchmann et al., 2004; Noeske et al., 2007; Elbaz et al., 2007; Karim et al., 2011; Whitaker et al., 2012; Speagle et al., 2014). Massive galaxies tend to form their stars earlier, a trend known as downsizing (Cowie et al., 1996; Heavens et al., 2004; Treu et al., 2005; Bundy et al., 2006; Neistein et al., 2006; Kriek et al., 2007; Conroy & Wechsler, 2009). Finally, the SFR has been connected to a number of physical properties of galaxies such as morphology (Postman & Geller, 1984; Wuyts et al., 2011) and environment, which turns out to be closely related to whether a galaxy is a satellite or central (Oemler, 1974; Dressler, 1980; Peng et al., 2010, 2012).

One might hope that these global observations would strongly constrain the SFHs of individual galaxies, but this connection is not easily established. For example, average SFHs can be inferred by integrating the main-sequence SFR over time (e.g., Leitner, 2012), but this approach leads to inconsistencies (Leja et al., 2015). Instead, the most successful theoretical models link the growth of stellar mass to the growth of the dark matter halos that galaxies inhabit, for example, via subhalo abundance matching (e.g., Kravtsov et al., 2004; Conroy et al., 2006; Moster et al., 2013; Behroozi et al., 2013b), halo occupation distributions (e.g., Peacock & Smith, 2000; Seljak, 2000; Hearin et al., 2016), semi-analytic models (Kauffmann et al., 1993; Somerville et al., 2001; Guo et al., 2011), or other assumptions (Bouché et al., 2010; Davé et al., 2012; Lilly et al., 2013; Tacchella et al., 2013; Mitra et al., 2017). One important conclusion from these models is that there has to be significant scatter between halo and galaxy masses (and thus growth histories) in order to explain observations (More et al., 2009; Behroozi et al., 2013b; Reddick et al., 2013; Gu et al., 2016). As a result, even models that agree on global constraints can lead to orthogonal interpretations of the evolution of individual galaxies. A good example for such disagreement is the “rapid quenching111The term “quenching” is somewhat ambiguous. In this paper, we use it to mean the cessation of star formation, without any presumption as to whether the decrease happens quickly or slowly, and whether it happens due to a diminishing gas supply or other physical processes.” framework where galaxies follow the main sequence until they sharply fall below the main sequence (Peng et al., 2012; Wetzel et al., 2013; Tinker et al., 2016; Tacchella et al., 2016b) and the “stochastic” framework where correlated scatter and the central limit theorem lead to the main sequence (Kelson, 2014; Kelson et al., 2016).

Recently, another rather different picture has been proposed. Inspired by the fact that the global SFRD is well fit by a log-normal in time, Gladders et al. (2013, hereafter G13) suggested that this form might also describe the SFHs of individual galaxies (G13; Dressler et al., 2013, 2016; Oemler et al., 2013a; Abramson et al., 2015, 2016). The log-normal SFR is given by the expression

[TABLE]

where $A$, $T_{0}$, and $\tau$ are free parameters (throughout the paper, $t$ refers to the time since the big bang, not lookback time). G13 emphasize that the key assumption need not be the exact functional form of the log-normal, but rather its steep rise and slow decline in linear time, suggesting a physical picture different from main sequence star formation interrupted by sudden quenching. While there is some evidence that the majority of galaxies cease their star formation gradually (Noeske et al., 2007; Schawinski et al., 2014; Peng et al., 2015; Eales et al., 2017; Gutcke et al., 2017), the log-normal SFR is no more than an assumption. The appeal of this assumption is that it allows for wide-ranging predictions if the log-normal parameters for a sample of galaxies can be inferred. This procedure was implemented by G13 who found surprisingly good agreement of the predicted stellar mass functions and their evolution, the star formation main sequence, downsizing, and many more complicated metrics (G13; Abramson et al., 2016). These successes cannot trivially be reproduced with symmetric forms of the SFH such as a Gaussian in linear time (G13).

Given the scarcity of reliable, time-resolved SFH observations, cosmological simulations of galaxy formation can help to differentiate between the different physical pictures. These simulations have recently reached a reasonable agreement with a number of observables, providing some level of confidence in their predictions for individual galaxies (Vogelsberger et al., 2014a; Torrey et al., 2014; Schaye et al., 2015; Davé et al., 2016). More specifically, Sparre et al. (2015) showed that the galaxy population of the Illustris simulation broadly matches the observed main sequence, and investigated individual SFHs using principal component analysis (see also the related analyses of Simha et al. 2014 and Cohn & van de Voort 2015).

In this paper, we systematically investigate the fundamental assumption of G13, namely whether the log-normal functional form is a good fit for SFHs in the Illustris simulation and for the inferred SFHs of Pacifici et al. (2016b, hereafter P16). We find that log-normals fit the majority of Illustris galaxies very well, particularly in the mass range studied by G13. We investigate the log-normal parameter space of normalization, peak time, and width as a common language for simulated and observed SFHs. Our goal is not to test the fitting procedure of G13 in detail or to compare the galaxy population in Illustris to that of G13. Instead, we study which physical properties of Illustris galaxies translate into particular values of the log-normal parameters. We also compare the log-normal to other fitting functions and discuss the implications of the log-normal framework in terms of quenching and the global star formation properties of the universe.

The paper is organized as follows. In Section 2, we describe the simulated and observed data our analysis is based on. In Section 3, we investigate the quality of the log-normal fits and the resulting parameter space. In Section 4, we discuss the implications of our findings for the global star formation properties of the universe and the quenching of star formation. We summarize our results in Section 5. In Appendices A and B we give mathematical expressions of various fitting functions and discuss additional details regarding fits and fit results.

Throughout the paper, we extensively use the language of log-normal functions (Equation 1). This functional form imposes the constraint that the rise and decline of star formation are symmetric in logarithmic time. Hence, the parameters $T_{0}$ and $\tau$ are in units of logarithmic time, and can be interpreted as the time when half the stars have formed and the duration of the galaxy’s star formation. However, the log-normal does not peak at $t=\exp(T_{0})$, and we thus substitute the more intuitive parameter space of the SFR’s peak time, $t_{\rm peak}$, and its full width at half maximum, $\sigma_{\rm SFR}$. Figure 1 illustrates how these parameters determine the shape of the log-normal. Finally, the normalization $A$ corresponds to the total integrated star formation, a quantity we re-cast as the total stellar mass formed,

[TABLE]

where $f_{\rm ret}=0.6$ is the retention factor due to stellar evolution (see Section 2.3). The numerical factor accounts for the fact that SFRs are measured in $\>{M_{\odot}/{\rm yr}}$ whereas times are in gigayears. Detailed expressions for a number of useful properties of log-normal SFHs are given in Appendix A.1.

2. Simulation, Data, and Methods

In this section, we introduce the observed and simulated datasets our analysis is based on, as well as our method for extracting and fitting SFHs from simulations.

2.1. The Gladders et al. Galaxy Sample

Our primary observed dataset consists of the stellar masses, SFRs, and best-fit log-normal parameters of G13. The underlying galaxy observations were taken from a number of galaxy surveys. In particular, the stellar masses and SFRs for galaxies with $M_{*}>4\times 10^{10}\>{M_{\odot}}$ were drawn from the Sloan Digital Sky Survey (see Oemler et al., 2013b, and references therein), whereas the PG2MC survey (Calvi et al., 2011) was used for galaxies with $M_{*}<4\times 10^{10}\>{M_{\odot}}$. Due to the smaller volume of the latter survey, the samples were re-normalized to create equal weights. The resulting sample contains $2094$ galaxies with a mean redshift of $0.0678$, and is complete above $M_{*}=10^{10}\>{M_{\odot}}$. We refer the reader to G13 and the references therein for details on how the stellar masses and SFRs were computed.

The fundamental assertion of G13 is that the log-normal, a functional form that describes the global SFRD of the universe, might also be a good description of the SFHs of individual galaxies. However, for each observed galaxy, only two variables are known ($M_{*}$ and SFR) whereas the log-normal has three free parameters (Equation 1). Thus, G13 used additional global constraints, namely the SFRD (accounting for the contribution from galaxies with $M_{*}<10^{10}\>{M_{\odot}}$) and SFR distributions back to $z\approx 1$. Some of these observations were drawn from the IMACS Cluster Building Survey (ICBS, Oemler et al., 2013b; Dressler et al., 2013; Oemler et al., 2013a). In a combined fit over the global constraints and the stellar masses and SFRs of each galaxy, they found the best-fit $T_{0}$ and $\tau$ for the $2094$ galaxies (see their Figure 9). The normalization of the log-normals was set such that the integrated SFR matches the observed $M_{*}$, meaning that one has to assume a retention factor, $f_{\rm ret}$, i.e. the ratio between the stellar mass initially formed to the stellar mass that survives until the galaxy is observed. G13 assumed a retention factor of $0.6$ which we adopt throughout this paper.

2.2. The Pacifici et al. Inferred SFHs

P16 derived SFHs by fitting the multi-band photometry of 845 quiescent galaxies with spectral energy distributions (SEDs) computed from a large library of theoretical SFHs. The SED library consists of 500,000 simulated galaxy SEDs, created by applying the semi-analytic model of De Lucia & Blaizot (2007) to the merger trees from the Millennium Simulation (Springel et al., 2005, see also Pacifici et al. 2012), and predicting the emission using the stellar population synthesis model of Bruzual & Charlot (2003). This library is flexible enough not to introduce a particular shape of the SFH a priori. P16 restricted themselves to quiescent galaxies, defined to have specific star formation rates

[TABLE]

where the sSFR is in units of $\rm yr^{-1}$ and $t_{\rm obs}$ is the cosmic time at observation in gigayears. For the comparisons with the Illustris and G13 samples, we set $z_{\rm obs}=0$. We consider the median SFHs in six bins in redshift and six bins in mass given by P16 (see their Figure 5).

2.3. The Illustris Simulation

The Illustris simulation (Vogelsberger et al., 2014a) follows a comoving cosmological volume of 106.5 cubic Mpc. The cosmological parameters were set according to the WMAP9 cosmology of Hinshaw et al. (2013), and the same cosmology was adapted for all calculations in this paper. The simulation was run using the moving-mesh code Arepo (Springel, 2010) which includes physical models for gas cooling, star formation, metal enrichment, black hole growth, as well as feedback from stellar winds, supernovae, and AGNs (Vogelsberger et al., 2013, 2014b, 2014a; Genel et al., 2014; Torrey et al., 2014; Sijacki et al., 2015). For the purposes of this paper, the star formation prescription is of particular interest. In Illustris, stars are formed according to a sub-grid model, which stochastically places “star” particles in gas cells that exceed a threshold number density of $0.13\ \rm cm^{-3}$ (Springel & Hernquist, 2003; Vogelsberger et al., 2013). Each star particle represents a population of stars born with a Chabrier initial mass function (Chabrier, 2003). The timescale of star formation depends on the inverse root of the density, leading to a relation between surface density and star formation that roughly follows the observed Kennicutt-Schmidt relation (e.g., Kennicutt, 1998).

As they age, simulation star particles return mass to the interstellar medium via stellar winds and supernovae (Vogelsberger et al., 2013). Star particles therefore have masses that evolve with time to values lower than their initial birth mass. For any given galaxy, this decrease causes the integrated SFR to be higher than the sum of stellar particle masses. Throughout the paper, we distinguish these two stellar mass definitions carefully. The mass decrease in a stellar population can be expressed as a retention factor, which we find to vary between $0.54$ for the oldest and $0.62$ for the youngest stellar populations in Illustris, with a mean of $0.57$. We note that this is very close to the value of $0.6$ assumed by G13, meaning that we can directly compare the cumulative SFRs of the observed and simulated galaxy samples without incurring a large error.

Dark matter halos, and the galaxies that inhabit them, are identified using the Subfind algorithm (Davis et al., 1985; Springel et al., 2001; Dolag et al., 2009). The various properties of galaxies are computed based on the bound matter within twice the the stellar half-mass radius of a subhalo as identified by Subfind (Vogelsberger et al., 2014a). We use a number of galaxy properties listed in the Subfind catalogs, including various mass and size definitions, metallicity, black hole mass, and properties of the parent group if applicable. Furthermore, we use the stellar assembly data provided by Rodriguez-Gomez et al. (2016), for example, the fraction of a galaxy’s stellar mass that was formed in situ, ex situ, and in various types of mergers.

We consider all galaxies with $M_{*}\geq 10^{9}\>{M_{\odot}}$, a sample of \num29203 galaxies (\num19375 of which are centrals at $z=0$, and $9828$ satellites). For the comparisons with the observational galaxy sample described in Section 2.1, we use the sub-set of high-mass galaxies with $M_{*}\geq 10^{10}\>{M_{\odot}}$, a sample of $6947$ galaxies ($4643$ centrals and $2304$ satellites). In both samples, the satellite fraction is almost exactly one-third. Given that the initial mass of star particles in Illustris is about $1.6\times 10^{6}\>{M_{\odot}}$, all galaxies are resolved by $600$ or more particles, and the high-mass sample by $6000$ or more particles. Sparre et al. (2015) tested the resolution dependence of star formation in Illustris using a lower-resolution run (Vogelsberger et al., 2014b), and found good convergence in key properties such as the star formation main sequence, the SFHs, and their variance.

2.4. SFHs and Log-normal Fitting

We extract the SFHs of Illustris galaxies using the same procedure as Sparre et al. (2015). In particular, we consider all stellar particles in the corresponding friends-of-friends subgroup and bin them by their birth times, weighted by their initial mass, in $100$ equal time bins between the first and last snapshots of the simulation (at $t=54$ Myr and $t=13.75$ Gyr, respectively). Defined in this way, the SFH includes the star formation in all progenitors that have merged with a galaxy during its lifetime. This definition matches observations of stellar mass and SFR at $z=0$, which do not distinguish between stars formed in situ (in the galaxy in question) or ex situ (in a smaller galaxy that merged). However, it does not match the observational definitions at higher redshift if a galaxy is to accrete a significant amount of stars at later times, though the effect of such mergers is significant only for high-mass galaxies (Conroy & Wechsler, 2009; Rodriguez-Gomez et al., 2016). In some figures, we use the SFR at $z=0$, which was computed from the SFR in all gas cells in a simulated galaxy. We have verified that this instantaneous SFR matches the SFH averaged over $20$ Myr very well.

We now investigate whether the tabulated SFHs are well fit by log-normal functions. We note that the term “fit” takes on a slightly unusual meaning in this context because the log-normal form is not intended to describe the details of an SFH. First, the SFHs in cosmological simulations suffer from shot noise due to the limited number of stellar particles formed in a given time bin. This noise grows as the bin size decreases, making it difficult to define a meaningful goodness-of-fit statistic such as $\chi^{2}$. Second, the spatial resolution (or force softening length) in Illustris is such that giant molecular clouds are not resolved, which is why the sub-grid model described in Section 2.3 is employed. The short-term structure of an SFH is thus significantly influenced by the characteristics of the sub-grid ISM model and resolution, as shown by high-resolution zoom-in simulations (e.g., Hayward et al., 2011; Torrey et al., 2012; Sparre et al., 2017). We note that other cosmological hydrodynamical simulations use similarly motivated sub-grid star formation models (Schaye et al., 2015; Davé et al., 2016), meaning their SFHs may be subject to comparable uncertainties on short timescales.

For these reasons, we disregard the short-term behavior of the simulated SFHs and instead focus on the overall, cumulative evolution of a galaxy’s stellar mass. We fit the cumulative SFR,

[TABLE]

with the integral of the log-normal SFR, which turns out to be given by a complementary error function (Appendix A.1). We use the Levenberg–Marquardt algorithm to minimize the square residuals between this functional form and the cumulative SFR in $100$ time bins. The values in these bins are highly correlated, which a standard $\chi^{2}$ statistic would not account for. We ignore this issue in the fit itself, but quantify the quality of the fit by computing the Kolmogorov–Smirnov statistic $D$, equivalent to the maximum difference of the cumulative functions at any time. We normalize this number by the total stellar mass created in a galaxy such that

[TABLE]

where $t_{0}$ is the age of the universe today. We use $D$ to quantify the quality of our fits in Section 3.1.

Finally, we note an inconvenient feature of the log-normal functional form: when only data from times well before the peak is available, the peak time and width are not well constrained, resulting in a strong degeneracy between $T_{0}$ and $\tau$. This issue manifests itself in a small fraction of Illustris galaxies that experience rising SFRs $z=0$ and are sometimes assigned unphysically high values of $T_{0}$ and $\tau$. We shall later return to the question of whether such galaxies exist in the real universe, but in a fitting sense, they cause inconvenient artifacts such as SFHs that peak hundreds of gigayears in the future.

Thus, we impose a prior on $t_{\rm peak}$ such that the fit is unaffected at $t_{\rm peak}<t_{0}$ (the time at $z=0$). Later peak times are penalized by multiplying the residuals by the term $(1+(\log_{10}t_{\rm peak}-\log_{10}t_{0})^{2})$ which effectively cuts off the distribution at 50 Gyr, corresponding to a scale factor $a_{\rm peak,max}=10$ in the Illustris cosmology. A similar, less informative prior penalizes extreme widths of $\sigma_{\rm SFR}>10^{2.5}$. The vast majority, over 90% of the galaxies in our sample, are unaffected by these priors. Galaxies in the high-$t_{\rm peak}$ tail experience drastically different best-fit parameters, namely changes of 50% or greater for 4% of the galaxies. However, due to the strong degeneracy of $T_{0}$ and $\tau$, the actual fitted SFHs of these galaxies are almost identical to those fit without a prior. Thus, the fit quality $D$ is barely affected: only 3% of the SFHs experience a change in fit quality of 20% or greater. The prior is even less important for the high-mass sample. For the inferences and comparisons shown in this paper, the tails in the $T_{0}$ and $\tau$ distributions are not important.

2.5. Halo Mass Accretion Histories

Besides the SFH, we also consider the mass accretion history (MAH) of a galaxy’s dark matter halo, which we extract from the Sublink merger trees provided by Rodriguez-Gomez et al. (2015). For each halo at $z=0$, we follow its most massive (or main branch) progenitor back in time and record its mass. In order to reduce the complex shape of the MAHs to characteristics such as a halo’s formation time, we fit the MAHs with the exponential form of Wechsler et al. (2002),

[TABLE]

where $M_{0}$ is the halo mass at $z=0$. This function was designed to fit the MAHs of isolated halos and is not equipped to capture the effects of subhalo accretion onto a larger host. Thus, we restrict the fits to epochs when the halos were not subhalos, which improves the correlation between halo formation redshift, and the SFH of their respective galaxies (Section 3.3). Following Wechsler et al. (2002), we define a formation redshift where the halo has formed half its mass, $z_{\rm wechsler}=\ln(2)/\alpha$. We have also experimented with the more flexible, three-parameter function of Tasitsiomi et al. (2004) and McBride et al. (2009), but find that the resulting formation redshift correlates less strongly with the SFHs.

Finally, we attempt to draw a more direct connection to the SFH fits by fitting the MAHs with a log-normal in time. We find that the log-normal is not particularly adapted to the shape of MAHs. The exponential rise at early times is generally well fit, but the additional free parameter compared to the Wechsler et al. (2002) function rarely improves the fit significantly. If the MAH flattens or even decreases, the power-law term in the Tasitsiomi et al. (2004) formula is a better description, and if the MAH is well described by a pure exponential, the third free parameter is unconstrained. For this reason, about a quarter of the fits fail to converge, even with the priors described in Section 2.4. We return to the relation between SFH and MAH log-normal fits in Section 4.5.

3. Results

In this section, we consider the quality of log-normal fits to Illustris and observationally inferred SFHs, and compare the best-fit parameters from the observational and simulated galaxy samples.

3.1. The Quality of Log-normal Fits

Figure 2 shows an array of example log-normal fits to a range of types of SFH, namely those of early and late-forming galaxies, centrals and satellites, slow and fast decliners, and galaxies that underwent a major merger. Despite the extremely different shapes of the corresponding SFHs, the log-normal captures the main features in the cumulative distribution while glossing over the noise and spikes in the SFR. Figure 2 highlights a large range of fit qualities, including the worst fit to any SFH in the high-mass sample ($M_{*}>10^{10}\>{M_{\odot}}$). In the lower-mass sample ($M_{*}<10^{10}\>{M_{\odot}}$), one particularly difficult case for the log-normal function is a sharp cut-off in star formation, for example the rapid quenching of satellites after infall into their host (see Figure 2 for an example, or Mistani et al. 2016 for an analysis of this effect in Illustris).

Figure 3 shows the distribution of $D$, the maximum fractional residual between the fit and the cumulative SFR. While $D$ ranges from almost zero to $34\%$, the vast majority of galaxies, namely $85\%$ ($92\%$ of centrals and $73\%$ of satellites), are fit to $5\%$ or better at all times. Only $1\%$ of the fits are off by $10\%$ or more at any time ($0.4\%$ of centrals and $2.7\%$ of satellites). The statistics are even better in the high-mass sample, where only $0.3\%$ of galaxies ($0.1\%$ of centrals and $0.8\%$ of satellites) exhibit a $D$ of $10\%$ or greater. Figure 3 also demonstrates that the tail toward poor fit qualities consists almost entirely of satellites. Visual inspection of these objects confirms that most of them are low-mass satellites that experienced rapid quenching after being accreted into a larger halo.

The generally excellent quality of the log-normal fits is by no means guaranteed, as we are fitting $100$ time bins with a function of only three free parameters. For comparison, we also consider the delayed-$\tau$ model (e.g., Gavazzi et al., 2002) and the double power law (e.g., Behroozi et al., 2013b) in Appendix A. We show that, on average, the log-normal function performs slightly better than the delayed-$\tau$ model, which has the same number of free parameters. In contrast, the double power law has four free parameters (normalization, peak time, and rising and falling slope), leading to a significantly more accurate fit, particularly for rapidly quenched satellites.

We also apply our fitting procedure to the inferred median SFHs of P16 (see Appendix B.1 and Figure 15). We find fit qualities very similar to those of the Illustris SFHs, with $D$ between 1.4% and 5%, and a median of 2.8% (compared to 3.1% for the Illustris sample).

In summary, the log-normal functional form provides an excellent description of the majority of the simulated and observationally inferred cumulative SFHs we considered. The key prediction of the log-normal is a steep rise and slow decline in linear time, and that the two timescales are coupled. This fundamental characteristic is shared with the delayed-$\tau$ model, meaning that many of the following results could probably be obtained based on delayed-$\tau$ fits as well. We conclude that a steep rise and slow decline are key features of the SFHs considered here, except for galaxies whose star formation is suddenly truncated by external processes such as satellite quenching. Models where the rise and decline timescales are independent (such as the double power law) can achieve a better fit for some galaxies, but have less predictive power: by fitting to only the rising or declining part of the SFH, one cannot infer the timescale of the other part. We further discuss this issue in Appendix B.2.

3.2. The Log-normal Parameter Space

We now consider the distribution of best-fit parameters in the G13, P16, and Illustris samples. Before attempting such a comparison, we need to check that the galaxy populations in Illustris and in the G13 sample are at least roughly comparable. Figure 4 shows the abundance of galaxies as a function of stellar mass and SFR in Illustris (blue) and the G13 sample (orange). The distributions of $M_{*}$ are relatively similar (left panel), but the SFRs (right panel) differ significantly. The G13 sample contains a much larger fraction of entirely quiescent, non-star forming galaxies, whereas Illustris contains an excess of galaxies with SFRs between $1$ and $10\>{M_{\odot}}/\rm yr$. This disagreement is not surprising because the G13 sample is more or less representative of the quiescent and star-forming fraction in the $z=0$ universe (Abramson et al., 2015), whereas an excess of blue, star-forming galaxies in Illustris at $z=0$ has been noted before (Vogelsberger et al., 2014b). Furthermore, Sparre et al. (2015) detected a paucity of starburst galaxies in Illustris compared to observations. Conversely, we remind the reader that the log-normal parameters of G13 are inferred from a semi-global fit rather than measured from resolved SFHs. With these caveats in mind, we emphasize that our goal is not to test the galaxy population of Illustris or the fitting procedure of G13, but rather to compare the regions of log-normal parameter space occupied by the two samples, and to investigate the physical properties of galaxies that occupy those regions. For those purposes, the populations shown in Figure 4 are sufficiently similar.

We compare the best-fit parameters of the Illustris and G13 samples in Figure 5. The top row shows histograms for each parameter. While the distributions of final stellar masses roughly agree, we notice that Illustris galaxies peak slightly later (higher $t_{\rm peak}$) and have significantly larger widths (higher $\sigma_{\rm SFR}$). We note that $M_{\rm final}$ is a reflection of the stellar mass of galaxies, modulo their current age: for early forming galaxies, $M_{\rm final}\approx M_{*}(z=0)$, whereas late-forming galaxies have yet to form some of the stars that contribute to $M_{\rm final}$.

The differences in the distributions are partially due to the excess of blue, star-forming galaxies in Illustris compared to the G13 sample (Figure 4). The fits to the most extreme of those cases tend to be poorly constrained and are influenced by our prior, which discourages $t_{\rm peak}>t_{0}$ (Section 2.4), causing tails of a few percent of the population toward very high values in all three parameters (M_{\rm final}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}10^{12}, t_{\rm peak}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}t_{0}, \sigma_{\rm SFR}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}30). The high values of $M_{\rm final}$ are an artifact: although the corresponding galaxies have relatively low stellar masses at $z=0$, their fitted log-normal SFHs predict a large amount of future star formation, increasing their $M_{\rm final}$.

The bottom two rows of Figure 5 show the correlations between each pair of parameters. The $M_{\rm final}$-$t_{\rm peak}$ (left column) and $M_{\rm final}$-$\tau$ (right column) planes show the expected trends with stellar mass: larger galaxies form earlier (have lower $t_{\rm peak}$) and decline relatively fast (have lower $\sigma_{\rm SFR}$), in agreement with the downsizing picture. The population of extremely late-forming galaxies in Illustris manifests itself as tails in the top right hand corners of both panels. At low stellar masses, the scatter in peak time increases. Interestingly, these trends are less significant in the G13 sample, especially the relation between mass and peak time.

We now consider the relation between peak time and width (center column of Figure 5). Once again ignoring the Illustris population with a peak time later than $z=0$, the G13 and Illustris samples occupy a similar region of parameter space. In both samples, earlier-forming galaxies form faster, but the correlation is more significant in Illustris. Furthermore, at fixed $t_{\rm peak}$ Illustris galaxies have broader widths, i.e., they form stars for a longer time than inferred by G13.

We note another potentially revealing difference between the G13 and Illustris samples: the absence of young galaxies or “late bloomers” (G13; Oemler et al., 2013a; Dressler et al., 2016) in Illustris. Those are galaxies that form the majority of their stars relatively quickly after $z\approx 1$. The corresponding parameter space of $t_{\rm peak}$ between $z=1$ and $z=0$ and \sigma_{\rm SFR}\lower 3.01385pt\hbox{;\buildrel<\over{\sim};}5 Gyr is visibly populated in G13, but deserted in Illustris (though a few such galaxies exist in the lower-mass sample, Section 3.3). This difference is another manifestation of the relatively wide SFHs in Illustris: if galaxies peak as late as $z=1$, they are most likely still active today (at least according to the log-normal fits).

As an independent, observational check, we compare the G13 and Illustris samples to the log-normal parameters of the P16 median SFHs in Figure 6. The color scale of the P16 points indicates their observation redshift. The G13 and P16 samples match very well and occupy the same $t_{\rm peak}$-$\sigma_{\rm SFR}$ relation. The lower scatter in the P16 sample is expected since those parameters were derived from median SFHs rather than individual galaxies. However, P16 considered only quiescent galaxies as defined in Equation 3. The corresponding limits at each redshift are shown with lines in Figure 6, and points of a given color should be compared to the Illustris and G13 galaxies below the corresponding line. We note that the P16 points themselves lie slightly above the lines at each redshift, indicating that the log-normal fits overestimate the sSFR at the redshift of observation (see Appendix B.1 and Figure 15). Nevertheless, the agreement between the observational samples is reassuring, and lends credibility to the inferred parameters of G13.

3.3. Which Galaxy Properties Determine the SFH?

Having established that our simulated and observed galaxy samples occupy similar areas of the log-normal parameter space, we now ask what galaxy formation physics causes particular shapes of the SFH: are they pre-determined by initial conditions such as the density environment, influenced by the history of a galaxy’s halo, or do they arise from the complicated physics of galaxy formation? For this investigation, we consider the entire Illustris sample with $M_{*}>10^{9}\>{M_{\odot}}$.

3.3.1 The Peak Time–Width Relation

Figure 7 reproduces the $t_{\rm peak}$–$\sigma_{\rm SFR}$ relation for the entire Illustris sample. The correlations with $M_{\rm final}$ are very similar to those in the high-mass sample shown in Figure 5. Figure 7 demonstrates that the $t_{\rm peak}$–$\sigma_{\rm SFR}$ relation is well approximated by a power law,

[TABLE]

with a standard deviation of 0.18 dex (about 50%). The line fit and scatter for the $M_{*}>10^{10}\>{M_{\odot}}$ sample are virtually identical. According to this relation, the SFHs in Illustris follow a scaling relation between their peak time and width, where early forming galaxies also form quickly. We notice a few outliers, for example, relatively late, short starburst SFHs (see Figure 2 for an example), but overall the relation is surprisingly tight. P16 found a similar relation between SFH width and the age of the universe, but they referred to the time of observation rather than peak time (their Figure 6).

Given how well the peak time of an SFH predicts its width, the question of which galaxy properties determine the SFH can be divided into two separate questions: which properties change peak time and width along their degeneracy, and which properties (if any) influence width at fixed peak time?

3.3.2 Centrals vs. Satellites

Part of the answer to these questions might be found in the differences between central and satellite galaxies. Figure 8 shows the distributions of their parameters compared to the overall sample. The distributions do not strongly depend on mass, except at very high masses, where most galaxies are early forming.

Figure 8 shows that even though satellites are, on average, slightly less massive than centrals (first panel), they form earlier (second panel), in contrast with the downsizing expectation that larger galaxies form earlier. Moreover, satellites form significantly faster than centrals (third panel), both because they are an intrinsically older population and because their star formation may be quenched after infall into their hosts (e.g., van den Bosch et al., 2008; Hearin & Watson, 2013; Wetzel et al., 2013).

However, the differences between the central and satellite populations do not appear to be responsible for the scatter in the $t_{\rm peak}$-$\sigma_{\rm SFR}$ relation. The right panel of Figure 8 shows the deviation from this relation, which is more or less symmetric for both centrals and satellites, though centrals have a slight tendency to lie above the relation (i.e., exhibit relatively wide SFHs at fixed $t_{\rm peak}$) and satellites below. The scatter around the relation is comparable for the two sub-samples, 0.16 dex for centrals and 0.2 dex for satellites.

3.3.3 The Galaxy–Halo Connection

As discussed in Section 1, we expect the MAH of a halo to have significant influence on the evolutionary history of its galaxy. Thus, we compare the log-normal parameters to a number of halo properties in Figure 9. The first impression from the gray histograms is that all correlations are subject to large scatter, particularly at low mass. This scatter implies that galaxy SFHs are diverse and not determined by any one parameter (see Dressler et al. 2016 for an observational investigation). Nevertheless, some trends are well-defined in an average sense. Many halo and galaxy properties correlate with halo and stellar mass, introducing a trivial correlation with peak time and width. Thus, we split the overall sample into three stellar mass bins. We have also verified that the central and satellite samples are similarly distributed for the correlations shown.

The left column of Figure 9 shows the correlations with the maximum mass a halo has obtained over its history (for subhalos, the mass today can be significantly lower than this maximum). The corresponding relations with stellar mass are very similar. For large, early forming galaxies, $M_{\rm final}$ is roughly equal to $M_{*}$ today, which manifests itself in a well-defined stellar mass–halo mass relation at the high-mass end. As expected, galaxies in more massive galaxies form earlier. At low halo or stellar masses, however, the correlation disappears (in agreement with Dressler et al., 2016). At first sight, this correlation seems to be in conflict with the expectation that large dark matter halos form late in hierarchical structure formation. However, while the halos of massive clusters are still growing today, their central galaxies are growing by mergers rather than star formation, meaning that their stellar populations whose ages we investigate here were largely formed in other halos (Rodriguez-Gomez et al., 2016). Thus, large halos form late whereas large central galaxies form early (e.g., Neistein et al., 2006).

The next variable we consider is the environment at $z=3$, i.e. whether a galaxy was born in an overdense or underdense region of the universe. We quantify density by the distance to the fifth nearest neighbor, but using a smoothed galaxy overdensity gives similar results. The trend that galaxies born in overdense environments form earlier holds for all stellar masses, and is a manifestation of the earlier collapse times of halos in overdense regions, known as “assembly bias” (Gao et al., 2005; Wechsler et al., 2006; Croton et al., 2007; Dalal et al., 2008; Zentner et al., 2014; Lin et al., 2016; Miyatake et al., 2016; Tinker et al., 2016, see also the arguments in Dressler 1980, Abramson et al. 2016, and Kelson et al. 2016). The environment at earlier and later redshifts shows similar relations with $t_{\rm peak}$, but the correlation weakens at late ($z<1$) and very early ($z>4$) redshifts.

We expect that the age of a galaxy should be strongly influenced by the age of its halo. We have measured halo age in a number of ways: as the redshift where half (or a quarter) of today’s mass has been accreted, by fitting with the MAH models of Wechsler et al. (2002) and Tasitsiomi et al. (2004), and as the starvation redshift $z_{\rm starve}$ suggested by Hearin & Watson (2013) which is defined as the maximum of $z_{\rm wechsler}$, the redshift where the halo first grew to $10^{12}\>{M_{\odot}}$, and the redshift when it was accreted if it is a satellite. The correlations with $z_{\rm wechsler}$ and $z_{\rm starve}$ are shown in the right columns of Figure 9. For all stellar masses, earlier halo formation does correspond to lower $t_{\rm peak}$ as expected, though the correlation exhibits large scatter. This scatter is not surprising given that the stellar mass–halo mass relation is subject to a scatter of about 0.2 dex (e.g., More et al., 2009; Behroozi et al., 2013b; Reddick et al., 2013; Gu et al., 2016).

Interestingly, $z_{\rm starve}$ exhibits a more pronounced (and more or less mass-independent) relation with $t_{\rm peak}$. By construction, $z_{\rm starve}$ is supposed to correspond to the redshift when a halo’s galaxy is expected to quench (Hearin & Watson, 2013). The strong correlation with $t_{\rm peak}$ means that the Illustris galaxy population supports the age matching picture in which the SFR (or color) of a galaxy is directly connected to its halo’s accretion history. However, this connection is non-trivial: massive, early forming galaxies tend to live in halos with a high $z_{\rm starve}$ (because the halo reached $10^{12}\>{M_{\odot}}$ early), whereas $z_{\rm wechsler}$ may indicate a late formation time for the same halo, leading to the large disagreement between the trends for high-mass galaxies (purple lines in Figure 9). All other definitions of formation redshift are less correlated with $t_{\rm peak}$ than $z_{\rm starve}$ and $z_{\rm wechsler}$.

So far, we have discussed how $t_{\rm peak}$ depends on stellar mass, environment, and halo age. However, what influences the duration of a galaxy’s star formation at fixed peak time? From the bottom panels in Figure 9, we conclude that mass and environment have little effect on whether an SFH lies above or below the peak time–width relation. It appears that very massive galaxies ($M_{*}>10^{11}\>{M_{\odot}}$) generally lie slightly below the relation, but this trend is not very strong. However, $z_{\rm wechsler}$ does correlate significantly, and for all stellar masses, with the deviation from the relation. In particular, at fixed peak time, galaxies in early forming halos tend to have relatively broad SFHs, while late-forming halos have relatively narrow SFHs. In Section 4.5, we demonstrate that this trend is generically expected if the the SFR is proportional to the halo mass accretion rate (MAR). We cannot detect a comparable correlation with $z_{\rm starve}$, indicating that $\sigma_{\rm SFR}$ is set by the MAH at times after $z_{\rm starve}$ to which $z_{\rm wechsler}$ is sensitive.

We note that these results do not allow us to directly assess the relative importance of mass quenching and environment quenching (e.g., Peng et al., 2010; Vogelsberger et al., 2014b). The strong correlation with $z_{\rm starve}$ (which relies on a physical halo mass scale) indicates that halo quenching plays an important role, but the formation history of halos is tied to their environment through assembly bias, making it hard to disentangle the two effects.

We conclude that the accretion history of a galaxy’s halo plays a key role in determining $t_{\rm peak}$ and $\sigma_{\rm SFR}$, in agreement with many theoretical models (Section 1). We discuss the galaxy–halo connection further in Section 4.5.

3.3.4 Ancillary Correlations

We now turn to physical properties of galaxies that we do not expect to fundamentally influence the SFH, but that exhibit interesting correlations with the SFH. To that end, the left column of Figure 10 shows correlations of the best-fit parameters with black hole mass at $z=0$. Even at fixed stellar mass, Illustris predicts a strong relation between black hole mass and $t_{\rm peak}$, where earlier-forming galaxies host more massive black holes. The population at high $M_{\rm final}$ and low black hole mass is due to to extremely late-forming galaxies, which are assigned an artificially high $M_{\rm final}$ despite their moderate stellar mass at $z=0$ (see the discussion in Section 3.2).

The accretion onto black holes cannot be resolved in a cosmological simulation, meaning that some of the correlation may be a manifestation of the seeding of black holes and their co-evolution with galaxies in the Illustris model. However, given that Illustris reproduces the observed black hole mass function (Sijacki et al., 2015), Figure 10 makes a clear prediction that earlier-forming galaxies host larger black holes (see also Sijacki et al., 2015; Bluck et al., 2016).

Another property of galaxies that shows a strong correlation with $t_{\rm peak}$ is their size, defined as the stellar half-mass radius at $z=0$ (right column of Figure 10). At a given mass, earlier-forming galaxies are significantly smaller today than their late-forming counterparts (Oesch et al., 2010; Valentinuzzi et al., 2010; Poggianti et al., 2013; Carollo et al., 2013; van der Wel et al., 2014; Fagioli et al., 2016). This correlation becomes insignificant for the largest Illustris galaxies with $M_{*}>10^{11}\>{M_{\odot}}$, presumably because their structure is significantly influenced by their merger history (Rodriguez-Gomez et al., 2017). We have also investigated a number of other galaxy properties such as metallicity, host halo mass (for satellites), merger history, and the fraction of stars formed in situ. Many of those properties correlate strongly with stellar mass, which drives their correlations with $t_{\rm peak}$.

4. Discussion

We have demonstrated that the log-normal functional form describes the majority of simulated and observationally inferred SFHs in the Illustris and P16 samples. We have shown that the log-normal parameter space of normalization, peak time, and width is useful for comparing the stellar histories of simulated and observed galaxies, and to infer the physical parameters that influence star formation in galaxies. In this section, we test the performance of the log-normal framework using more advanced metrics and discuss the theoretical motivation for log-normal SFHs.

4.1. The Global SFRD

Figure 11 compares various estimates of the global SFRD of the universe. The dashed lines show the observational estimates of G13 and Madau & Dickinson (2014). The SFRD of our Illustris sample (dark blue) underestimates the observed normalization due to contributions from galaxies with $M_{*}<10^{9}\>{M_{\odot}}$, which we do not account for. The Illustris SFRD is well matched by the sum of all log-normal fits (light blue). The SFRD of the G13 log-normal fits (orange line in Figure 11) matches their log-normal fit to the data almost perfectly. This match is expected because the SFRD was used as a constraint in their fits (any deviations are due to the re-sampling of the G13 galaxies, see Section 2.1).

We extrapolate the fitted SFRDs into the future (dotted lines), but keep two important caveats in mind. First, log-normal SFHs with significant future star formation cannot react to quenching events, meaning their limiting stellar mass is likely an overestimate of what one would find if Illustris were run into the future. Second, the prediction is affected by our prior on peak time: if left unconstrained, the most extreme cases of late-forming, high-width SFHs dominate the prediction. Even with very strict priors, however, the future SFRD of the log-normal fits does not decrease significantly compared to that shown in Figure 11, implying that the high future SFRD is not solely due to unconstrained fits. In fact, the log-normal extrapolation is not too surprising given that the slope of the Illustris SFRD is already shallower than the observations at $z\leq 1$. This deviation is not an artifact of the log-normal fit, and has been noted before (see Figure 8 of Vogelsberger et al. 2014b and Figure 6 of Torrey et al. 2014, or Figure 2 of Sparre et al. 2015).

4.2. Evolution along and across the Main Sequence

While we have shown that log-normals faithfully capture the overall, cumulative evolution of SFHs, we have yet to investigate metrics that are more sensitive to the instantaneous SFR. For example, the distinction between star-forming and quiescent galaxies is often drawn based on their position in $M_{*}$–sSFR space, where star-forming galaxies fall onto a relatively tight, evolving relation called the main sequence (e.g., Noeske et al., 2007). In this section, we investigate whether the log-normal fits correctly reproduce the $M_{*}$–sSFR evolution of Illustris galaxies.

Figure 12 summarizes various aspects of the $M_{*}$–sSFR trajectories. The gray background histogram in each panel shows the $z=0$ distribution of all Illustris galaxies in our sample, where the line at the bottom contains galaxies with $\log_{10}$ SFR below $-13$. The main sequence in Illustris is somewhat flatter than observations suggest (Sparre et al., 2015), preventing us from comparing the evolution of the SFHs relative to some main sequence fitting function (e.g., Speagle et al., 2014). Instead, we consider trajectories in $M_{*}$–sSFR space itself.

Panel (a) of Figure 12 shows the trajectories of four individual galaxies, which were randomly chosen from four mass bins indicated by different colors. The log-normal fits result in trajectories very similar to those of the SFHs themselves. Panel (b) addresses the same question in a more quantitative manner, by showing the median trajectories of all SFHs and fits. Not only do the fits reproduce the medians very well, the 68% scatter in the fits is also very similar to that in the SFHs. We do, however, note that the scatter of the fits in the high-mass bins ($M_{*}>11$) exceeds the scatter in the SFHs at early times. This difference can be understood by considering that the stellar mass of the galaxies is growing more rapidly than at late times, meaning that fit errors will manifest themselves in a larger deviation from the median.

In Section 3.1, we found that log-normals sometimes have trouble fitting the SFHs of satellites that cease their star formation abruptly. Thus, we separately investigate the median SFHs and fits of satellites in panel (c). Their trajectories do, as expected, reach lower sSFRs at late times, on average, but the fit quality is as good as for the entire sample. Even if the SFR decreases too slowly in a log-normal fit, this decrease will lead to a steep trajectory away from the main sequence because $M_{*}$ evolves very little. This means that the trajectories are essentially insensitive to the duration of quenching. We thus conclude that the good match of our fits in this parameter space is reassuring, but that quenching trajectories are easy to match regardless of how good the SFH fit is at late times.

In the same vein, we find that the $M_{*}$–sSFR trajectories of Illustris and G13 galaxies agree better than one might expect from the somewhat different best-fit parameter distributions (panel (d) of Figure 12). This agreement makes sense considering that Abramson et al. (2016) showed that the G13 sSFR evolution is compatible with the observed evolution of the main sequence normalization and slope, and that Illustris also roughly matches the observed evolution (Sparre et al., 2015).

Finally, the bottom panels of Figure 12 provide an intuition for where the log-normal trajectories in the different mass bins typically end up (dots indicate the end points at $z=0$, dotted lines the future development). The sSFR of log-normals decreases monotonically, and Figure 12 reaffirms that all log-normal trajectories end in a steep drop to low sSFRs eventually. In this picture, the scatter in the main sequence is naturally explained as an age gradient at fixed stellar mass (Abramson et al., 2016; Oemler et al., 2016).

4.3. Can Log-normals Capture Rapid Quenching?

In the previous section, we concluded that $M_{*}$–sSFR trajectories are very well reproduced by log-normal fits, but that these trajectories are insensitive to the timing of the sSFR evolution. In this section, we test this timing explicitly by considering the timescale over which star formation ceases. A galaxy is generally defined to be quiescent (or have quenched) if its sSFR that has fallen below a particular value well below the star formation main sequence, with the so-called green valley as an intermediate space between star-forming and quiescent galaxies. To assess the time evolution of log-normal fits, we measure their green valley crossing time, which we define to be the time an SFH spends between logarithmic sSFRs of $-10.5$ and $-11.5$ (see panel (a) in Figure 12). For this test, we consider only SFHs whose sSFRs have fallen below $-11.5$ in the last five time bins to ensure that they have reliably crossed the green valley. Furthermore, we require that the log-normal fit has also fallen below an sSFR of $-11.5$ (though this criterion makes relatively little difference to the selection). According to these criteria, about 3% of the Illustris galaxies have quenched today. While this fraction is lower than observed (Wetzel et al., 2013; Bluck et al., 2016), the actual number of galaxies does not matter for our comparison. We compute the time an Illustris SFH spent in the green valley by linearly interpolating between time bins. For the fits, we calculate the time spent in the green valley exactly.

Figure 13 shows the results of this exercise. The left panel shows the distribution of the actual green valley crossing times as measured from the SFHs. Satellites quench faster, on average, with a significant population that crosses the green valley in less than half a gigayear. The different timescales for centrals and satellites are expected and hint at different quenching mechanisms (e.g., Hahn et al., 2016). The right panel shows the difference between the fitted and actual crossing times. For centrals, the distribution is relatively symmetric, but systematically offset by $0.85$ Gyr. For satellites, this difference increases to $1.5$ Gyr, largely because the log-normals fail to reproduce the population of rapidly quenching SFHs. On the other hand, Abramson et al. (2016) showed that the G13 log-normals do reproduce the general increase of the average green valley crossing time with increasing age of the universe.

The overestimated quenching times represent the most serious limitation of log-normal fits that we have discovered. However, how important is this limitation? Some insight can be gained from abundance matching studies such as Wetzel et al. 2013 (see also van den Bosch et al. 2008 and Hahn et al. 2016). They find that above $M_{*}=10^{10}\>{M_{\odot}}$, the dominant avenue by which galaxies cease to form stars is intrinsic rather than satellite quenching. This result explains why only a small fraction of SFHs in this mass range were poorly fitted by the log-normal form. Below $10^{10}\>{M_{\odot}}$, satellite quenching accounts for a larger and larger fraction of quiescent galaxies (Wetzel et al., 2013). We did not investigate galaxies less massive than $10^{9}\>{M_{\odot}}$, but the importance of satellite quenching would presumably continue to increase with decreasing mass. For example, Geha et al. (2012) find that dwarf galaxies are more or less only quenched if they are satellites, indicating that the slower, intrinsic cessation of star formation becomes less important at low masses. Moreover, the model of Wetzel et al. (2013) infers very short quenching timescales for satellites, raising the possibility that quenching could occur even more rapidly than in Illustris. Unfortunately, their conclusions are hard to compare to our data because they rely critically on the observed fractions of green valley and quenched galaxies, which are not very well matched by Illustris (e.g., Bluck et al. 2016; see also, however, the recent observational discussion regarding the density of galaxies in the green valley, e.g. Chang et al. 2015 or Eales et al. 2017).

In summary, log-normal fits reproduce very well how galaxies evolve in $M_{*}$–sSFR space, but in some cases do not accurately capture how fast this evolution happens. As expected, this failure is most pronounced in rapidly quenched satellites, which represent a small fraction of the population at $M_{*}>10^{10}\>{M_{\odot}}$, but are more frequent at lower masses.

4.4. The SFH of the Milky Way

One possible application of the log-normal framework is to infer the stellar history of a particular type of galaxy. For example, Figure 14 shows the SFR and SFH of Milky-Way-like galaxies (hereafter MW analogs) as inferred from the log-normal fits to the Illustris (dark blue) and G13 (orange) samples. We have selected MW analogs to have $4<M_{*}<6\times 10^{10}\>{M_{\odot}}$ and $0.5<{\rm SFR}<2\>{M_{\odot}}/{\rm yr}$, resulting in samples of $136$ Illustris and $76$ G13 galaxies. The Illustris SFH is, of course, almost identical to what one would obtain by just summing the SFHs of the galaxies in the sub-sample without log-normal fitting. However, the log-normal fits to the G13 sample allow us to compare their prediction to the simulation.

Interestingly, the G13 Milky Way analogs are inferred to form later, with $t_{\rm peak}=6.1$ Gyr compared to $4.0$ Gyr for the Illustris sample. We further compare these predictions to the observational result of van Dokkum et al. (2013). They assumed a fixed number density of galaxies in order to infer which high-$z$ galaxies correspond to MW-like galaxies today (red dashed lines in Figure 14). While the assumption of constant number density does not hold, in general (Behroozi et al., 2013a; Torrey et al., 2015, 2017; Wellons & Torrey, 2017), tests with synthetic SFHs created using the Guo et al. (2011) semi-analytic model indicate that the systematic issues of the method, namely merging galaxies and scatter, lead to moderate errors on the inferred SFHs (Leja et al., 2013). The van Dokkum et al. (2013) result is in good agreement with the Illustris sample, with an even slightly earlier peak time of $3.4$ Gyr. We compare their measurement to the SFH prediction of Leja et al. (2013, red points).

We also show the SFH inferred by Patel et al. (2013) who used fits to the observed star formation main sequence and its evolution in order to determine MW progenitors. This selection is systematically different from the assumption of constant number density because the main sequence only takes star-forming galaxies into account. They investigated galaxies with a final mass of about $0.2$ dex lower than that of the MW, an offset we have added to their results in Figure 14. Interestingly, Patel et al. (2013) suggest a slightly later-forming MW than van Dokkum et al. (2013), in agreement with the G13 sample.

4.5. The Galaxy–Halo Connection Revisited

In Section 3.3.3, we established that the halo MAH is the most important factor determining the peak time and width of an SFH. As an extension of this result, we investigate whether the peak time–width relation arises simply as a consequence of halo MAHs in a $\Lambda$CDM universe.

The simplest possible connection we can assume is a steady-state model, where the SFR is proportional to the halo MAR (e.g., Bouché et al., 2010; Lilly et al., 2013; Dekel et al., 2013; Becker, 2015; Tacchella et al., 2016a). For concreteness, we quantify the MAH using the function of Wechsler et al. (2002). At z\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}1, we can neglect the influence of dark energy and use $a\propto t^{2/3}$ to convert Equation 6 into an expression for the MAR (e.g., Dekel et al., 2013),

[TABLE]

Because this expression has only two free parameters, it predicts a particular relation between the peak time and width of the MAH (where peak time now refers to the time when the MAR was highest). Interestingly, this relation is a power law with a slope of $1$, shallower than the value of $3/2$ measured for the Illustris SFHs. At $z<1$, we compute the peak time and width numerically, and find slopes even shallower than one. The results of this calculation are shown with gray dashed lines in Figure 9. The SFR $\propto$ MAR assumption significantly underpredicts the peak time at all formation redshifts, meaning that there is an offset between the formation times of halos and galaxies. However, it qualitatively explains the trend with respect to the mean peak time–width relation: at early formation times, the width of the halo MAH is larger than that of the SFH, but due to the shallow slope of the halo peak time–width relation this difference is reversed at late formation times.

We note that the shape of the Wechsler et al. (2002) expression for the halo MAR is similar to a log-normal for certain parameters. Thus, we attempt to fit the Illustris SFHs with log-normals. As mentioned in Section 2.5, a significant fraction of these fits fail. As before, the successful fits exhibit earlier peak times than the SFHs and obey a power-law relation between the peak time and width of the MAH. In the case of log-normal MAH fits, however, the slope is $1.35$, closer to the slope of $3/2$ measured for the SFHs, implying that the correlation between MAHs and SFHs might be even tighter than indicated by the best-fit parameters of the Wechsler et al. (2002) function.

In summary, both the shape of the SFHs and the peak time–width relation are encoded in halo MAHs to some degree, but the simple assumption that the SFR is proportional to the MAR cannot explain the SFHs quantitatively. It seems that baryonic processes act to delay the peak time of the SFHs, and to increase the width at later formation times. For example, a decrease in star formation is expected to follow a decrease in gas supply with some time delay, which might not be constant (e.g., Wetzel et al., 2013; Schawinski et al., 2014, see also the simulation results of Feldmann et al. (2016a, b)). In general, baryonic processes are expected to decouple galaxy and halo growth to some extent in lower-mass systems (e.g., Brook et al., 2012; Übler et al., 2014). Furthermore, the strength of the galaxy–halo connection in Illustris could be artificially enhanced because several physical processes in the Illustris model (such as black hole seeding and wind velocities) are tied to the dark matter properties of galaxies.

5. Conclusions

Inspired by the work of Gladders et al. (2013), we have investigated the log-normal function as a description of the SFHs of individual galaxies in both simulations and observations, and shown it to be a good fit in the majority of cases. We have discussed the strengths and failures of the log-normal picture, and connected the log-normal parameter space to the physical properties of galaxies and their halos. Our main conclusions are as follows.

Log-normals are an excellent fit to the SFHs of the vast majority of Illustris galaxies. While log-normals ignore short-term fluctuations in the SFR, the cumulative SFH is fit to within 5% deviation from the total mass formed for 85% of Illustris galaxies, and to within 10% for 99% of galaxies. The most poorly fit SFHs belong to rapidly quenched satellite galaxies. We find similarly good fit qualities for the observationally inferred median SFHs of P16. 2. 2.

Comparing different fitting functions, we find that the log-normal performs slightly better than the popular delayed-$\tau$ model, though both models predict steeply rising and slowly declining SFHs. While models with more than three free parameters (such as the double power law) result in even better fits, they are less predictive and more difficult to constrain in observational analyses like G13. 3. 3.

The Illustris, G13, and P16 galaxy samples occupy similar regions in the log-normal parameter space of total stellar mass, peak time, and full width at half maximum, though Illustris predicts more extended SFHs than the observations. Illustris galaxies exhibit a tight correlation between peak time and width, $t_{\rm peak}\propto\sigma_{\rm SFR}^{3/2}$, with only 0.18 dex scatter. 4. 4.

Log-normal fits correctly reproduce many global aspects of star formation, such as the SFRD and the evolution of galaxies along and across the star formation main sequence. The timing of this evolution, however, can be poorly matched for rapidly quenched satellite SFHs, leading to significantly overestimated quenching timescales. 5. 5.

Considering the physical properties of galaxies, we find that the formation history of the halo has the strongest influence on peak time and width, and that the starvation redshift of Hearin & Watson (2013) is a particularly good predictor of peak time. Generally, the trends are as expected: at fixed stellar mass, earlier-forming galaxies live in early forming halos, inhabit dense environments, host massive black holes, and have small sizes. Assuming that the SFR is proportional to the dark matter accretion rate can qualitatively explain the peak time–width relation, but does not quantitatively predict the correct peak times or widths. In general, SFHs are a diverse population with large scatter in all parameters, meaning that they are not uniquely determined by any one factor.

We have left a number of possible investigations for future work. For example, we did not evaluate the fidelity of the G13 fitting procedure on an object-to-object basis, a test that could be performed using a simulated mock dataset. Conversely, we did not investigate whether the log-normal assumption can help in recovering masses and SFRs from photometry or spectra (e.g., Simha et al., 2014). Finally, any analysis based on a single simulation is subject to the systematic uncertainty in the galaxy formation models used. Our analysis should be repeated using future cosmological simulations to ascertain the extent to which our results are specific to Illustris.

We are deeply grateful to all those who shared their data with us, namely the authors of G13 (Mike Gladders, Augustus Oemler, Alan Dressler, Bianca Poggianti and Benedetta Vulcani), Camilla Pacifici for sharing her SFHs, and the Illustris team for making their simulation public (Mark Vogelsberger, Shy Genel, Volker Springel, Debora Sijacki, Dandan Xu, Greg Snyder, Dylan Nelson, Lars Hernquist, and, in particular, Vicente Rodriguez-Gomez for his help with the merger trees). We thank Peter Behroozi, Charlie Conroy, Andrew Hearin, Lars Hernquist, Joel Leja, Camilla Pacifici, Ana-Roxana Pop, Laura Sales, Josh Speagle, and Sandro Tacchella for fruitful discussions and/or comments on a draft of this paper. B.D. gratefully acknowledges the financial support of an Institute for Theory and Computation Fellowship.

Appendix A Mathematical Details of SFH Fitting Functions

In this Appendix, we compute certain properties of a number of commonly used SFH fitting functions, such as their peak time, peak SFR, and cumulative SFR.

A.1. The Log-normal Function

The log-normal functional form is given in Equation 1. Counter-intuitively, $T_{0}$ is not the logarithm of its peak time, which is given by

[TABLE]

leading to a peak SFR of

[TABLE]

By integrating the log-normal, we find its cumulative SFH,

[TABLE]

where erf denotes the error function,

[TABLE]

The simplicity of this expression is extremely convenient and allows an intuitive interpretation of $T_{0}$ as the logarithm of the time when half the stars have formed, while $\tau$ determines how wide the error function extends before and after $T_{0}$ in logarithmic time units.

In addition to peak time, we derive a more intuitive description of the width of an SFH in time, namely the interval between the times when the SFH first and last reaches a fraction $1/\beta$ of its peak value (Equation A2). The interval between these times in logarithmic space is

[TABLE]

Translating back into linear space,

[TABLE]

In this paper, we use the full width at half maximum in linear time, or $\beta=2$, to determine the width of the SFH,

[TABLE]

A.2. The Delayed-$\tau$ Model

One of the oldest suggestions for the form of SFHs is an exponentially declining SFR (e.g., Tinsley, 1972; Bruzual A., 1983), physically motivated by the idea of a fixed gas supply that is gradually being exhausted. However, it is clear that this function will not fit simulated SFHs well because it corresponds to an infinitely sharp rise of the SFR from zero to a peak value. An extension of this model is the delayed-$\tau$ model (e.g., Gavazzi et al., 2002; Lee et al., 2010; Simha et al., 2014),

[TABLE]

which represents a linear increase in the SFR at early times and an exponential decrease at late times. Before the initial time $t_{\rm i}$, no stars are formed. Here, we have normalized the function such that, as for the log-normal, the total stellar mass formed is $10^{9}\times A$. The SFR peaks at a time $t_{\rm peak}=t_{\rm i}+\tau$ with ${\rm SFR}_{\rm peak}=A/\tau/e$. The cumulative SFR is

[TABLE]

The delayed-$\tau$ model can be simplified further by fixing $t_{\rm i}$, for example, to 1 Gyr as suggested by Simha et al. (2014). However, in Section B.2, we will show that even the delayed-$\tau$ model with three free parameters has trouble capturing the rising SFHs in a number of cases, meaning that a fixed $t_{\rm i}$ would lead to unacceptable fits.

A.3. The Double-power-law Model

The double power law (e.g., Behroozi et al., 2013b) models the SFH as

[TABLE]

Due to its four free parameters, this function is more flexible than the log-normal or delayed-$\tau$ models. Unfortunately, some of the following mathematical expressions are somewhat cumbersome. We find a peak time of

[TABLE]

and a resulting peak SFR of

[TABLE]

The most inconvenient expression is that for the cSFR. If we define $x\equiv t/\tau$,

[TABLE]

where

[TABLE]

and

[TABLE]

where ${}_{2}F_{1}$ is the hypergeometric function. While relatively straight-forward to implement numerically, this function is not exactly intuitive. Moreover, we could not derive a simple expression for the total stellar mass formed, making the interpretation of the normalization $A$ more complicated than in the cases of the log-normal or delayed-$\tau$ models.

Appendix B More Details on the Log-normal Fits

In this appendix, we discuss certain aspects of our SFH fits in detail, including the fits to the P16 inferred SFHs, the relative fit qualities of different fitting functions, and an alternative parameter space for the log-normal functional form.

B.1. The Pacifici et al. Sample

Figure 15 shows log-normal fits to the median SFHs inferred by P16 in six redshift and six mass bins. Some of the bins are omitted because they contained too few objects. Visually, the quality of the fits appears satisfactory: the log-normals tend to rise slightly faster than the inferred median SFHs, but the fits are well within the 50% confidence curves at all times. The distribution of the fit quality indicator $D$ is similar to the Illustris sample, with $D$ between 1.4% and 5%, and a median of 2.8% (compared to 3.1% for the Illustris sample).

B.2. Comparison of Fit Qualities

We have fit all galaxies in the Illustris sample with the function forms discussed in Appendix A. Except for a very small number of failed delayed-$\tau$ fits, best-fit parameters could be obtained for each SFH and fitting model. Figure 16 shows a few representative cases. For the majority of SFHs, all three fitting functions give similar results, though with subtle differences. For example, the delayed-$\tau$ model generally struggles to capture the earliest times because of the sudden onset of star formation in that model (second SFH in Figure 16). Due to its free declining slope, the double power law significantly outperforms the other models when sudden quenching leads to a sharp drop in the SFH (third SFH in Figure 16). However, even the double power law cannot capture all types of quenching as shown in the rightmost panels of Figure 16.

We quantify the goodness of fit for each model in Figure 17. The fit quality of the log-normals is systematically better than for the delayed-$\tau$ model, but the log-normal is outperformed by the double power law. This difference is expected because the double power law has an extra free parameter: while the slope of the rise and decline in star formation are linked in the log-normal and delayed-$\tau$ models, they are independent in the double power law. This freedom allows for more accurate fits, but also reduces the predictive power of the model. For example, if data is available for only the early, rising part of an SFH, the declining slope of a double-power-law fit is unconstrained, meaning the function makes no prediction for the declining time scale (e.g., for very late-forming galaxies in Illustris). This issue makes the double power law somewhat less attractive for observational investigations such as G13 where limited data is available for each galaxy.

B.3. The Original Log-normal Parameter space

Throughout the paper, we have used a log-normal parameter space of peak time and full width at half maximum. G13, however, use the original parameter space given by Equation 1, i.e. the half-mass time $T_{0}$ and the logarithmic width $\tau$. In order to allow a direct comparison between the Illustris results and the figures in G13, Figure 18 shows the same comparison as in Figure 5, but in $T_{0}$–$\tau$ space. Our normalization $M_{\rm final}$ only differs from G13’s $A$ by a multiplicative factor, meaning the correlations with $A$ are very similar to those shown in Figure 5. In $T_{0}$–$\tau$ space, the peak time–width correlation appears somewhat weaker, supporting the notion that the $t_{\rm peak}$–$\sigma_{\rm SFR}$ parameter space is easier to interpret.