Are star formation rates of galaxies bimodal?

TL;DR

This paper challenges the traditional view of bimodal star formation rate distributions in galaxies, proposing a unified model using zero-inflated negative binomial distributions that naturally incorporate zero SFR galaxies and asymmetric scatter.

Contribution

It introduces a novel, unified statistical model for galaxy SFR distributions that accounts for zero SFR galaxies and avoids the need for bimodal classification.

Findings

Zero-inflated negative binomial distributions fit galaxy SFR data well.

Excluding 'dead' galaxies, the SFR distribution is unimodal.

Measurement uncertainties can mimic bimodality.

Abstract

Star formation rate (SFR) distributions of galaxies are often assumed to be bimodal with modes corresponding to star-forming and quiescent galaxies, respectively. Both classes of galaxies are typically studied separately and SFR distributions of star-forming galaxies are commonly modelled as lognormals. Using both observational data and results from numerical simulations, I argue that this division into star-forming and quiescent galaxies is unnecessary from a theoretical point of view and that the SFR distributions of the whole population can be well fit by zero-inflated negative binomial distributions. This family of distributions has 3 parameters that determine the average SFR of the galaxies in the sample, the scatter relative to the star-forming sequence, and the fraction of galaxies with zero SFRs, respectively. The proposed distributions naturally account for (i) the discrete nature of star formation, (ii) the presence of 'dead' galaxies with zero SFRs, and (iii) asymmetric scatter. Excluding 'dead' galaxies, the distribution of log SFR is unimodal with a peak at the star forming sequence and an extended tail towards low SFRs. However, uncertainties and biases in the SFR measurements can create the appearance of a bimodal distribution.

Equations

$$\mathrm{SFR}=\mathrm{SFC}{m}_{\mathrm{SFC}}/t_{\mathrm{av}}.$$

$$P_{\mathrm{NB}}(k;\mu ,\theta )=\left(\genfrac{}{}{0ex}{}{\theta +k-1}{\theta -1}\right){\left(\frac{\mu }{\theta +\mu }\right)}^k{\left(\frac{\theta }{\theta +\mu }\right)}^{\theta }.$$

$$P_{\mathrm{zNB}}(k;\pi ,\mu ,\theta )=\pi \delta {}_{0k}+(1-\pi )P_{\mathrm{NB}}(k;\mu ,\theta ),$$

Full text

1 Introduction

The bimodal colour distribution of nearby galaxies leads to a natural classification into ‘blue’ and ‘red’ galaxies (e.g., Baldry et al. 2004; Brammer et al. 2009). On average, blue (red) galaxies in the local Universe have high (low) star formation rates (SFRs) although dust extinction complicates this basic picture (Whitaker et al., 2012; Taylor et al., 2014). Hence, galaxies are often divided based on their level of star formation activity into star-forming and quiescent galaxies (e.g., Balogh et al. 2004; Moustakas et al. 2013). However, whether the distribution of SFRs is also bimodal remains an open question (e.g., Elbaz et al. 2007; McGee et al. 2011).

The SFRs of star-forming galaxies strongly correlate with their stellar masses resulting in a well-defined ‘star-forming sequence’ (Noeske et al., 2007). In contrast, quiescent galaxies have generally very low (or vanishing levels) of SFRs whose exact amount is challenging to infer observationally (e.g., Brinchmann et al. 2004; Utomo et al. 2014; Chang et al. 2015). Hence, most observational studies focus on the SFRs of ‘star-forming galaxies’ alone and model their distribution at fixed stellar mass with a lognormal distribution (e.g., Noeske et al. 2007) or with the sum of two lognormal distributions (Sargent et al., 2012). The intrinsic scatter around the star-forming sequence is found to be about $\sim{}0.3-0.4$ dex, essentially independent of redshift ($z\sim{}0-6$; e.g., Chang et al. 2015; Salmon et al. 2015; Schreiber et al. 2015; Shivaei et al. 2015). The approach of approximating the distribution of SFRs of star-forming galaxies with a lognormal is also adopted in most analyses of galaxy formation simulations (e.g., Schaye et al. 2015; Sparre et al. 2015; Dave et al. 2016; Feldmann et al. 2016).

However, this approach has a number of shortcomings. First, a clear separation of galaxies into star-forming and quiescent is challenging in practice. Classifications based on colour–magnitude diagrams suffer from a large population of dust-obscured star-forming galaxies with colours intermediate between blue and red (Salim et al., 2009; Taylor et al., 2014; Chang et al., 2015). While colour–colour diagrams offer a more robust alternative (Wuyts et al. 2007), the mapping from colours to SFRs can be biased by relatively small amounts of recent star formation, by dust, and by the presence of evolved stellar populations (e.g., Salim et al. 2009; Wuyts et al. 2011; Fumagalli et al. 2014). For instance, while high-redshift galaxies classified as quiescent based on colour–colour diagrams have typically significantly reduced SFRs (Man et al., 2016; Straatman et al., 2016), some of them show non-negligible levels of SFR and dust extinction (e.g., Brammer et al. 2009; Spitler et al. 2014; Mancini et al. 2015). It is thus legitimate to ask whether quiescent and star-forming galaxies are actually two separate populations or whether galaxies simply form a continuum from low to high specific SFRs (sSFRs) without a natural dividing point.

Secondly, a lognormal distribution predicts a symmetric scatter around the star-forming sequence, in contrast to the predictions of many observational and simulations studies (Brinchmann et al., 2004; Chang et al., 2015; Dave et al., 2016). Often this difference is attributed to an ‘imperfect separation’ into star-forming and quiescent galaxies (e.g., Chang et al. 2015). Instead, I argue that the level of asymmetry of the scatter contains important information about galactic star formation and should not be ignored.

Thirdly, star formation is correlated in space and time as stars are typically born in clusters (Lada & Lada, 2003). These correlations introduce a level of discreteness into the star formation process. For instance, star clusters in the Milky Way have masses ranging from about $\sim{}10^{2}-10^{3}$ $M_{\odot}$ for open clusters to $\sim{}10^{5}$ $M_{\odot}$ for young massive clusters (Portegies Zwart et al., 2010). Star clusters can be even more massive in star-bursting galaxies (Zhang & Fall, 1999) and at high redshift, when galaxies are generally more gas rich (e.g., Swinbank et al. 2010; Kruijssen 2012). The discreteness effect is further exacerbated in observationally inferred SFRs as those trace energy injections from rare, massive stars (Kennicutt, 1998). A lognormal distribution, which is a continuous distribution, does not account for this discrete mode of star formation.

Finally, the star formation activity, especially for higher redshift and/or lower mass galaxies, is thought to be bursty (e.g., Dominguez et al. 2015; Sparre et al. 2016). Consequently, many galaxies may experience intermittent episodes of low or vanishing SFR (Feldmann et al., 2017). However, galaxies with fully suppressed SFRs cannot be modelled by a lognormal distribution. Thus, in practice, galaxies with SFRs below the detection limit are excluded from (non-stacked) analyses of the star-forming sequence even if such galaxies are star forming according to their colours (e.g., Whitaker et al. 2014).

As I argue in this Letter, these shortcomings can be mitigated by dropping the assumption of a lognormal SFR distribution and by not dividing galaxies into star-forming and quiescent galaxies in the first place. In particular, I propose to replace the lognormal ansatz with (zero-inflated) negative binomial distributions (NBDs). This family of distributions found wide applicability in particle physics (e.g., multiplicity distributions of charged particles in hadronic collisions; Alner et al. 1985), astrophysics (e.g., the number of globular clusters in galaxies and event rates of fast radio bursts; De Souza et al. 2015; Wiel et al. 2016), and cosmology (e.g., modelling count-in-cell distributions and void probability functions; Carruthers & Duong-Van 1983; Gaztanaga 1992) but, to my knowledge, has not been used to model distributions of SFR or sSFRs.

Choosing an appropriate model for the distributions of SFRs is not self-evident as the origin of the scatter around the star-forming sequence is not well understood. As gas accretion on to galaxies and star formation within galaxies are likely linked (e.g. Bouché et al. 2010; Sánchez Almeida et al. 2014), the scatter may be related to variations in the gas accretion rates (Dutton et al., 2010; Forbes et al., 2014) or halo growth rates (Feldmann & Mayer, 2015; Feldmann et al., 2016; Rodríguez-Puebla et al., 2016). The scatter may also arise from random stochasticity (Kelson, 2014), gas fraction variations (e.g., Magdis et al. 2012; Saintonge et al. 2012; Scoville et al. 2016), changes in the efficiency of star formation (e.g., Genzel et al. 2010; Saintonge et al. 2012), a natural diversity in star formation histories (Gladders et al., 2013; Dressler et al., 2016), or combinations of some of these processes (e.g., Tacchella et al. 2016; Feldmann et al. 2017). Hence, it appears justified to explore empirically how well SFRs follow various basic distributions.

With the recent availability of large numbers of reliable SFR measurements these basic distributions can be compared against observations and numerical simulations. I use spectral energy distribution (SED) based SFR estimates of galaxies in the local Universe based on the Sloan Digital Sky Survey (SDSS) (Chang et al., 2015), ultraviolet (UV) and infrared (IR) based SFR measurements of galaxies at $z\sim{}2$ from 3D-HST (Brammer et al., 2012; Skelton et al., 2014), as well as SFRs measured in cosmological simulations (Feldmann et al., 2016) that are part of the Feedback in Realistic Environments (FIRE) Project (Hopkins et al., 2014).

This Letter is organized as follows: Section §2 introduces probability distributions to model SFRs in galaxy samples. Section §3 shows that zero-inflated NBDs (zNBDs) provide adequate approximations to SFR distributions in observations and numerical simulations. I discuss the implications for a possible SFR bimodality in the final section.

2 Modelling the distribution of SFRs

At fixed stellar mass, SFRs of star-forming galaxies are typically assumed to obey lognormal distributions (e.g., Chang et al. 2015), i.e., the logarithm111In the following, $\log$, $\lg$, and $\ln$ denote the logarithm to an arbitrary base, to base 10, and to base ${\rm e}$, respectively. of the ${\rm SFR}$ is assumed to be a continuous variable that is normally distributed with standard deviation $\sigma\equiv\sigma_{\ln{\rm SFR}}=\sigma_{\lg{\rm SFR}}\ln{}10$. The mean, median, and most probable value of $\lg{\rm SFR}$ coincide, and they define the position of the star-forming sequence for the given stellar mass.

Instead, I propose to model SFRs of galaxies with NBDs. As NBDs describe count data, I assume that the star formation activity over time $t_{\rm av}$ consists of individual star formation events, each adding mass $m_{\rm SFC}$, i.e.,

[TABLE]

${\rm SFC}$ (‘star formation count’) is a non-negative integer-valued random variable with a [potentially zero-inflated (discussed later)] NBD.

The probability mass function (PMF) of an NBD is specified by two parameters, e.g., the expected count $\mu$ and a shape parameter $\theta$, both positive real numbers. The probability of the outcome ${\rm SFC}=k\in\mathbb{N}$ is (e.g., Cameron & Trivedi 2013)

[TABLE]

If $\theta$ is a positive integer, eq. (2) describes the probability distribution of the number of Bernoulli trials (each with a success probability $\theta/[\theta+\mu]$) undertaken before the $\theta$th trial is successful. Cases with $\theta=1$ are known as geometric distributions. There exist many other characterizations of NBDs (Boswell & Patil, 1970). For the remainder of this Letter, I will assume $\theta=O(1)$ and $\mu\gg{}\theta$.

NBDs are often used to model count data that is ‘overdispersed’ relative to a Poisson distribution222A Poisson-distributed random variable has a variance equal to its mean value $\mu$. In contrast, the variance of a random variable with an NBD is equal to $\mu+\mu^{2}/\theta$, i.e., the variance at given $\mu$ can be adjusted by changing $\theta$.. NBDs are Poisson–Gamma mixtures which suggests a simple physical interpretation: The overdispersion arises from galaxy-by-galaxy variations of the expected number of SFCs even among galaxies of the same stellar mass. This is plausible as the expected number of SFCs in a galaxy should vary depending on the state of its interstellar medium.

Zero-inflated models (Mullahy, 1986; Lambert, 1992) increase the probability of obtaining zero counts, i.e., of producing galaxies with zero SFRs. For zNBDs,

[TABLE]

where $\pi\in{}[0,1]$ parametrizes the excess probability to obtain zero counts: $P_{\rm zNB}(0;\pi,\mu,\theta)-P_{\rm NB}(0;\mu,\theta)=\pi\,[1-P_{\rm NB}(0;\mu,\theta)]$. Candidate processes responsible for $\pi>0$ include galactic outflows powered by starbursts or active galactic nuclei (e.g. King & Pounds 2015; Somerville & Davé 2015, and references therein) and strong environmental effects (Gunn & Gott, 1972).

Examples of lognormal, negative binomial, zero-inflated negative binomial, and Poisson distributions are shown in Fig. 1. For illustrative purposes I convert between PMFs and probability density distributions (PDFs) by approximating probability point masses with intervals of uniform probability density. As the figure shows, a count of zero is the most probable outcome for random variables with (z)NBDs if $\theta\leq{}1$. Interestingly, while the probability of obtaining a certain count value is unimodal and decreases monotonically with the count value for (z)NBDs, the distribution of the logarithm of the count variable shows a well-defined peak near333Provided $\mu\gg{}\theta$, the peak of $\lg$ count is near $\lg(\mu+\theta/2)\sim{}\lg\mu$. In contrast, for a lognormally distributed random variable $X$, the distributions of $\lg{}X$ and $X$ peak at $\lg\mu-\frac{\sigma^{2}}{2\ln{}10}$ and $\mu{}/{\rm e}^{1.5\sigma^{2}}$, respectively.

the logarithm of $\mu$. In addition, a conversion to $\log$ counts requires that zero and non-zero counts are treated as separate components. As I argue in §4, this split into two components may lead to an apparent ‘bimodality’ of $\log$ SFR distributions.

Fig. 1 also illustrates that the distribution of $\log$ counts around the peak is asymmetric for (z)NBDs. There is a significant tail towards lower values, resulting in a non-negligible probability of a zero count outcome if $\theta<1$. Furthermore, the degree of asymmetry and the width of the $\log$ count distribution increases with decreasing $\theta$ (see the right-hand panel). There, I plot how the upward scatter ($\Delta{}_{+}$) and the downward ($\Delta{}_{-}$) scatter444$\Delta_{+}$ and $\Delta_{-}$ are defined as follows. Let $\lg{}c_{*}$ denote the position where the PDF of $\lg\,{\rm counts}$ reaches its maximum value $p_{*}$. Increasing (decreasing) the counts by $\Delta_{+}$($\Delta_{-}$) dex relative to $c_{*}$ results in a decrease of the PDF by $\chi={\rm e}^{-1/2}$ relative to $p_{*}$. The factor ${\rm e}^{-1/2}$ is chosen such that $\Delta_{+}=\Delta_{-}=1$ standard deviation for a normal distributed PDF of $\lg\,{\rm counts}$. For $\mu\gg{}\theta$, $\Delta_{\pm}=\lg(-W_{\mp}(-\chi^{1/\theta}/{\rm e}))$ where $W_{+}$ and $W_{-}$ are the principal and the -1 branch of the Lambert $W$ function, respectively. relative to the peak of the PDF scale with $\theta$. For $\theta<1$, upward scatter and downward scatter differ significantly from each other, while for $\theta\gtrsim{}1$, $\Delta{}_{+}\sim{}\Delta{}_{-}\sim{}\theta^{-1/2}/\ln{}10$.

3 Comparison with observations and simulations

I use three samples of galaxies at different redshifts, with different stellar masses, and from different sources to test whether SFR distributions at fixed stellar mass can be approximated by zNBDs.

The first sample comprises 38246 nearby galaxies from SDSS. SFR and stellar mass estimates are based on multiwavelength (UV to mid-IR) SED modelling (Chang et al., 2015). I select a mass-complete sample of galaxies555Excluding galaxies with flag $\neq$ 1, i.e., those without reliable aperture corrections, WISE photometry, or SED fits. with stellar masses in the range of $10^{10}-10^{11}$ $M_{\odot}$ and with redshifts $z<0.05456$. I set the SFR of a galaxy to zero if the best-fitting sSFR is very low ($<3\times{}10^{-12}$ yr*-1*) and the modelling error is at least twice the spread of the star-forming sequence.

The second sample contains 2317 galaxies at $z\sim{}2$ from the 3D-HST survey (Brammer et al. 2012; Skelton et al. 2014; catalogue v4.1). SFRs in the catalogue are based on UV$+24\mu{}$m luminosities (SFRs of non-detected galaxies are set to zero) while stellar masses are derived from SEDs fits. I select all galaxies666Excluding galaxies with star_flag = 1, near_star = 1, or use_phot $\neq$ 1. from the five available fields with stellar masses in the range of $10^{10}-10^{11}$ $M_{\odot}$ and with redshifts $z=1.5-2$. For both observational samples, I convert SFR into SFC with the help of eq. (1) and by adopting a conversion factor777This choice corresponds to, e.g., $m_{\rm SFC}=2\times{}10^{4}$ $M_{\odot}$ and $t_{\rm av}=20$ Myr. However, provided $k\gg{}1$ and $\mu\gg{}\theta$, (z)NBDs are well approximated by (zero-inflated) gamma distributions and the PDF of $\lg\,{\rm SFR}$ for SFR¿0 approaches $\ln(10)\,x\,P_{\rm Gamma}(x;\alpha=\theta,\beta=1)$ with $x=\theta\,{\rm SFR}/\langle{}{\rm SFR}\rangle{}$, i.e., it is independent of the conversion factor. of $m_{\rm SFC}/t_{\rm av}=10^{-3}M_{\odot}/{\rm yr}$.

The last sample combines SFRs and stellar masses from cosmological galaxy simulations. Specifically, I use 1648 $z=6$ galaxies with $M_{\rm star}=10^{7}-10^{9}$ $M_{\odot}$ from the MassiveFIRE simulation suite (Feldmann et al., 2016, 2017). Gas and stellar components of the simulated galaxies are represented by gas and star particles with masses $m_{\rm b}=3.3\times{}10^{4}$ $M_{\odot}$. Star formation occurs probabilistically based on the local conditions of the interstellar medium. Each individual star formation event results in the formation of a star particle of mass $m_{\rm b}$. I refer the reader to Hopkins et al. (2014) for a background on the simulation methodology. Stellar masses and SFRs are measured within radii of 0.1$$R_{\rm vir} of the primary dark matter haloes hosting a given galaxy, excluding satellites. SFRs are averaged over the past 20 Myr and converted into SFCs via eq. (1) with $m_{\rm SFC}=m_{\rm b}$ and $t_{\rm av}=20$ Myr.

The data sets contain galaxies of a range of stellar masses to increase the sample size. To combine the SFR distribution of galaxies with different stellar masses, I convert SFRs into SFCs and perform generalized linear regressions of the SFCs as a function of $\ln{}M_{\rm star}$ with a $\log$ link function888Statistical analyses are carried out with R (https://www.r-project.org) using the standard glm function to fit Poisson and geometric distributions, the glm.nb function from the MASS package to fit NBDs, and the zeroinfl function from the pscl package to fit zNBDs. An example regression script is provided at http://www.ics.uzh.ch/~feldmann/resources.html.. I thus simultaneously fit for the position of the star-forming sequence as a function of stellar mass and constrain the parameters for the assumed distribution of SFCs at fixed stellar mass.

The Akaike information criterion (AIC; Akaike 1974) measures how well, relative to each other, different statistical models describe a given data set. Among the SFC distributions I tested (negative binomial, geometric, lognormal, Poisson distributions as well as their zero-inflated versions), zNBDs performed best under the AIC metric. Fig. 2 clarifies why zNBDs work so well. The $\lg\,{\rm sSFR}$ distribution is highly asymmetric with a tail towards low values. Furthermore, there is a significant fraction of galaxies with vanishing SFRs or SFRs below the detection limit. These properties are captured by zNBDs but not by, e.g., lognormal distributions.

I tested the sensitivity of the fitted model parameters (overdispersion parameter $\theta$, stellar mass scalings of the average SFR and of the excess probability $\pi$) by analysing additional MassiveFIRE-based samples for different redshifts and mass resolutions. The main findings are as follows: The slope and normalization of the star-forming sequence show significant changes with redshift (as expected) and only slight changes with particle resolution. The excess probability $\pi$ is not affected by redshift or mass resolution. The overdispersion $\theta$ does not depend on redshift but varies mildly with mass resolution in a not obviously systematic way.

Finally, I also explored the effect of varying the time interval over which SFRs are averaged. As expected, reducing $t_{\rm av}$ results in a smaller number of SFCs and a larger fraction of galaxies with zero SFRs. Specifically, $\pi$ increases from 20% to 52% as $t_{\rm av}$ is reduced from 100 Myr to 5 Myr. Moreover, $\theta$ increases, i.e., the scatter of the star-forming sequence decreases, with increasing averaging time. However, the change in $\theta$ is relatively modest (a factor of 1.9 when $t_{\rm av}$ increases from 20 to 100 Myr). Furthermore, $\theta$ remains unchanged if $t_{\rm av}$ is lowered to 5 Myr. This suggests that the star formation activity of MassiveFIRE galaxies is strongly correlated on $\lesssim{}20$ Myr time-scales.

4 Discussion and Conclusions

(Zero-inflated) NBDs predict significant numbers of galaxies with vanishing SFRs. As discussed in §2, these ‘dead’ galaxies form a separate component upon log-transforming SFRs, while the remaining galaxies have a unimodal log SFR distribution with an extensive tail towards low SFRs. The reader may ask whether these findings are consistent with claims of an SFR bimodality (e.g., Elbaz et al. 2007; McGee et al. 2011).

I address this question in Fig. 3 in a schematic way. A detailed analysis is left for future work. Using the best fitting zNBD of the SDSS-based sample shown in the top panel of Fig. 2, I create a mock sample consisting of the SFRs of 40000 galaxies with $M_{\rm star}=10^{10.5}$ $M_{\odot}$ (black histogram). I then subject this sample to non-linear biases and measurement uncertainties999For the bias, I assume that sSFR below ${\rm sSFR}_{0}=3\times{}10^{-10}$ yr*-1* are overestimated: ${\rm sSFR}_{\rm obs}={\rm sSFR}_{0}\,({\rm sSFR}/{\rm sSFR}_{0})^{0.463}$, following Utomo et al. (2014). To mimic measurement uncertainties, each SFR is replaced by $\max(0,{\rm SFR}+X)$, where $X$ is a normally distributed random variable with zero mean and standard deviation of $0.01$ $M_{\odot}$ yr*-1*.. The latter moves galaxies out of the ‘dead’ pile, thereby introducing in a second peak in the $\log$ SFR distribution (red histogram), while biases distort the apparent shape of the star-forming sequence (blue histogram). The combination of measurement noise and bias can create the illusion of a bimodal SFR distribution101010 Similarly, colour bimodalities can arise if sSFRs are non-linearly mapped to colours, e.g., if ${\rm colour}=\tanh(3\times{}{\rm sSFR}/{\rm sSFR}_{\rm MS}-1)$. with low numbers of galaxies at intermediate-to-low SFRs.

Fitting the distribution of SFRs with (z)NBDs is straightforward (see footnote 8), and it offers substantial benefits compared with the current standard approach of fitting only the star-forming sub-sample with lognormals. I recommend its use for the modelling of SFR distributions both in observations [e.g., in the simplified form of (zero-inflated) gamma distributions, see footnote 7] and in simulations (where the conversion parameters $m_{\rm SFC}$ and $t_{\rm av}$ are given).

Acknowledgments

RF thanks the referee for valuable comments. RF acknowledges financial support from the Swiss National Science Foundation (grant no 157591). 3D-HST observations were taken by the 3D-HST Treasury Program (GO 12177 and 12328) under NASA contract NAS5-26555. Simulations were run under allocations SMD-14-5492, SMD-14-5189, SMD-15-5950 (NASA HEC) and AST120025, AST150045 (NSF XSEDE).