THE IMACS CLUSTER BUILDING SURVEY. IV. THE LOG-NORMAL STAR FORMATION HISTORY OF GALAXIES

Current data (see, for example, the compilation in Cucciati et al. 2012) clearly show both a rise and fall in the SFRD over cosmic time. More specifically, the rise in the SFRD is fast at early times, with a long continuing but declining tail of star formation to the present epoch. Ignoring mergers, the shape of the SFRD should represent the shape of the SFH of the mean galaxy in the universe. Of course, the actual SFH of any recognized single galaxy at z = 0 will likely not fit this description, not least because a typical z = 0 galaxy is thought to represent the merger of a number of smaller systems over cosmic time.

Putting aside this complexity for the moment, we note that the presence of both a rise and fall in the SFRD also indicates that simple SFH models which decline from initially high star formation rates (i.e., exponentially declining models—referred to as τ models—commonly used in modeling the spectrophotometric properties of galaxies at later times) cannot easily describe the overall SFRD evolution; a delay time must be allowed—such as in the delayed exponential SFH of Sandage (1986)—in order to even approximate current data. The rise in the SFRD from high redshift can be produced by an ensemble of simple τ models if a second parameter—the time at which such a SFH "turns on"—is allowed and is appropriately distributed at early times.

Recent analysis of the photometric properties of distant galaxies also argues for a SFH in individual galaxies that is smoothly rising at early times (Maraston et al. 2010). They found that an SFH of ${\rm SFR}=Ae^{-(t-t_o)/\tau }$, a so-called inverted-τ model, produced a better fit to the photometric properties of a sample of z ∼ 2 actively star forming galaxies than did the more typically used declining-τ models. A similar result, for galaxies with a fixed comoving number density of 2 × 10⁻⁴ Mpc⁻³, was found to hold from z = 8 z to z = 3 by Papovich et al. (2011); Reddy et al. (2012) find that rising SFHs at 7 > z > 2 are the preferred history for a sample of galaxies observed at z ∼ 2. Maraston et al. also note that the results of Cimatti et al. (2008) from a study of z ∼ 1.6 elliptical galaxies also argue for an early SFH that is rising in time. Additionally, some studies suggest that the sSFR of galaxies is flat beyond the peak of the SFRD at z ∼ 2 out to at least z ∼ 7 (e.g., Stark et al. 2009; González et al. 2010), though see de Barros et al. (2012) for an alternate view; one simple interpretation of such a flattening is that individual galaxies at these redshifts have SFHs with a star formation rate smoothly increasing to later times. This overall picture of rising SFHs in individual galaxies at high redshift is also found in recent cosmological hydrodynamic simulations (Finlator et al. 2011; Jaacks et al. 2012).

Maraston et al. point out that the remarkably ubiquitous use of declining τ models in the literature over the past few decades is not well justified, and has become thoroughly divorced from the original aim of these models of measuring the age of early-type galaxies at z ∼ 0. It thus seems reasonable to seek a new simple model which captures both the rising early-time and declining late-time behavior of galaxy SFHs.

In seeking a functional form to describe the SFRD over cosmic time, we have explored a number of possible descriptions, and have found that a remarkably simple one—a log-normal in time—works extremely well. We describe the scale free log-normal distribution of the star formation rate as

$$\begin{equation} {\rm SFR}(t,t_0,\tau)=\frac{1}{t\sqrt{2\pi \tau ^2}}e^{-\frac{(\ln t-t_0)^2}{2\tau ^2}}, \end{equation} \tag{ 1 }$$

where t is the elapsed time since the big bang, t₀ is the logarithmic delay time, and τ sets the rise and decay timescale. One great advantage of such a form over some other models—such as the inverted-τ model of Maraston et al. (2010) is that it naturally subsumes both a rising and falling SFH at different times; simple τ models—inverted or otherwise—are aphysical, whereas a log-normal in time could arguably describe the SFH of a galaxy (or at least a galaxy's components) at all times. A further advantage is that the delay time is de-coupled from the width of the SFH, unlike the delayed exponential model of Sandage (1986). Behroozi et al. (2013) have recently suggested two other functional forms—a double power law, or a hybrid exponential+powerlaw—which we will explore further below.

Figure 1 shows a fit to the SFRD with a single log-normal. This simple functional form produces an reasonable fit, with a χ² of 2.1. It is worth re-emphasizing that the log-normal, besides providing an good fit, is a functional form that emerges again and again in the analysis of distributions in natural systems. From the failure rate of electronic components (Salemi et al. 2008) to the latency periods of infectious diseases (Kondo 1977) and many other systems, the log-normal distribution appears as the preferred rate model. Log-normal distributions occur when multiplicative effects dominate; that the cosmic SFRD is log-normal in time suggests a deeper meaning than that it is simply a good fit.

Zoom In Zoom Out Reset image size

Figure 1 also shows fits to several other functional forms suggested in the literature:

• 1.  
  the best fitting delayed exponential SFH, given by the SFH∝t/τ²exp (− t²/2τ²); the fit is obviously poor—with a χ² of 10.2. This result is not unexpected, given the coupling of delay time and width in this model with a single shape parameter;
• 2.  
  the best fitting exponential + power-law (Behroozi et al. 2013), given by SFH∝t^Aexp (− t/τ); this model, with the same number of parameters as the single log-normal, produces a poorer fit, with a χ² of 4.3;
• 3.  
  the best fitting double power-law (Behroozi et al. 2013), given by SFH∝((t/τ)^A + (t/τ)^−B)⁻¹; this three-parameter model fits the distribution extremely well, with a χ² of 1.7.

The apparent agreement between different datasets seen in Figure 1 suggests that the uncertainties on the data as reported are overestimated—or include a significant correlated systematic component. We caution that a reduced χ² is thus not trivial to compute; however the relative χ² of models with the same number of parameters—i.e., the single log-normal and the exponential + power-law is still informative. In this comparison the log-normal is clearly preferred.

Figure 1 also shows the result of a double log-normal fit (where each log-normal component has equal weight). This fit resolves any lingering tension between the peak height and late-time SFRD present in the single log-normal fit. The χ² of this fit is 1.6. We include the double log-normal fit here primarily because it yields an interesting physical interpretation; one of the fitted SFHs is characterized as an early onset with a fast decline and the other as a somewhat later onset with a slower decline. These two basic SFHs could be interpreted as corresponding to early- and late-type galaxies, or even the bulge and disk components of the typical galaxy. We also use the double log-normal fit as a smooth description of the cosmic SFRD.

The simple model fits above have been computed using a downhill simplex, but in general in this paper we fit ensembles of parameterized SFH models to data using a simulated annealing algorithm; the choice of simulated annealing to thoroughly explore the parameter space for large aggregate samples of SFHs will become apparent in the next section.

The z = 0 sSFR of a given galaxy is the current star formation rate divided by total mass, and the total mass is the integral of the SFH for that galaxy, adjusted for stellar mass loss. The z = 0 sSFR distribution is thus a reasonably orthogonal measure to the SFRD over cosmic time, and a potentially useful constraint on a SFH model. We illustrate this point in Figure 2, which shows the τ, t₀ parameter plane for a log-normal SFH, with lines of constant z = 0 sSFR and lines of constant time of peak star formation shown. The displayed range of sSFRs brackets values seen in several samples of galaxies discussed below. Over much of that range, lines of constant sSFR and lines of constant peak time are approximately orthogonal—i.e., imposing both a history and a current sSFR selects one particular curve for most galaxies. Only galaxies with large sSFRs at the present epoch are poorly constrained in this way.

Zoom In Zoom Out Reset image size

As discussed in Paper III, an evolution in the rate of starbursts in galaxies is insufficient to explain the observed evolution in the sSFR distributions from low to high redshift. The presence of a main sequence of star formation at both low and high redshifts has led numerous authors (e.g., Noeske et al. 2007; Rodighiero et al. 2011) to conclude that starbursts are subdominant in affecting the observed evolution in star formation in galaxies, as is also seen in simulations (Di Matteo et al. 2008). While it is apparent that secondary starbursts do happen, our aim in this paper is to explore the ability—or lack thereof—of smooth SFH models to reproduce the observed trends in sSFR, completely absent any starbursts.

Motivated by Figure 2, we consider the z ∼ 0 sSFR distribution for galaxies, for which we use the local galaxy sample described in Paper III. Two sub-samples are included. The first, taken from the PG2MC survey (Calvi et al. 2011), covers a larger volume, but has a higher mass limit of 4 × 10¹⁰ M_☉, with minimum and maximum redshifts of 0.03 and 0.11 and a median redshift of 0.0918. The second, from the Sloan Digital Sky Survey (SDSS) observations of the northern galactic cap, is more restricted in redshift, with minimum and maximum redshifts of 0.035 and 0.045 and a median redshift of 0.0401. This second sub-sample has a lower mass limit of 1 × 10¹⁰ M_☉, and is cut at the upper end at the lower limit of the first sub-sample. Galaxies are weighted to bring the two sub-samples to a common volume. As described in Oemler et al. (2013a; Paper I of this sequence) masses for this second sub-sample have been computed using a variant of the technique in Bell & de Jong (2001) with delayed exponential models (see Section 2.1) as the underlying form of the SFH.

The total sample is 2094 galaxies, distributed in sSFR versus mass as shown in Figure 3, with a mean weighted redshift of 0.0678. The sSFR values are computed from Hα fluxes. The detectability of star formation depends on Hα equivalent widths and as a result the sSFR limit is not trivial to describe, as it relies both on the flux of the Hα line in emission, as well as the continuum strength, including the depth of the Hα line in absorption. The resulting incompleteness does not appear exactly fixed in either star formation rate (this would be the expected result if only the Hα line flux were relevant, for example) or sSFR (which would be expected if only the Hα equivalent width were relevant). Figure 3 also shows the fraction of galaxies which are measured as having a sSFR identically zero as a function of mass; as expected early-type systems with no measurable Hα emission are proportionately more common at the high mass end. For the purposes of this paper, we estimate the threshold sSFR—i.e., the allowed upper limit of the actual value of the sSFR for a galaxy of a given mass measured in the data in Paper III to have sSFR = 0—as simply a fixed star formation rate of 0.05 M_☉ yr⁻¹. Note that the exact choice of limits does not significantly affect any of the results which follow.

Zoom In Zoom Out Reset image size

To create a sample of SFHs that match both the z ∼ 0 sSFR distribution and the cosmic SFRD, we proceed as follows. We consider a simulated sample of 2094 galaxies with masses identical to the sSFR data discussed above. The SFH for each galaxy is described by two parameters, τ and t₀. We jointly solve for these parameters for each galaxy in the ensemble using simulated annealing, requiring at the same time that the mass- and sample-weighted sum of the individual SFHs match the shape of the double log-normal model fit to the cosmic SFRD as detailed in Figures 1 and 2. We do not use the raw SFRD data detailed in Figure 1, but the smooth fit to these data, since this modeling process produces an under-constrained realization rather than a unique best-fit model. Galaxies with a sSFR measured to be zero in the data are allowed to take any value up to the mass-dependent threshold shown in Figure 3.

The sSFR value for each simulated galaxy at its measured redshift is computed simply as the ratio of the star formation rate divided by the mass, where the latter is the integral of the former from early times to the time of observation, modified downward by stellar mass loss that occurs over the galaxy's history. We take the functional form of this mass loss from Jungwiert et al. (2001) and as in Maraston et al. (2010) scale the mass loss so that the total loss at 10 Gyr is 40%. We do not attempt to compute the mass loss individually for each galaxy, convolved across its SFH; rather, for computational simplicity and efficiency we simply compute the mass loss from the peak of the SFRD at z ∼ 2 to the epoch of observation for each galaxy. The differences between this approach and a more refined computation are in general small because much of the mass loss that occurs at early times, and by z ∼ 0 the bulk of the star formation in galaxies is well in the past.

The resulting distribution of sSFR values, and the distribution of τ and t₀ parameters, is shown in Figure 4. The fit is reasonable; however, Figure 4 reveals a certain amount of tension between the SFRD constraints and the z = 0 sSFR constraints when fitting each galaxy with a log-normal SFH. Specifically, the SFRD shape prefers a somewhat higher normalization of the z = 0 sSFR distribution than the actual measurements suggest. This is apparent in the fit shown in Figure 4 in two ways: (1) the final fitted values of the sSFR for each galaxy are on average slightly high, and (2) the galaxies with sSFRs that are nominally zero in the data are fit with sSFR values that pile up at the adopted upper threshold for the allowed sSFR values.

Zoom In Zoom Out Reset image size

The source of this tension is, effectively, the mass limitation of the input z ∼ 0 dataset, coupled to the changing contribution of galaxies of various luminosities to the cosmic SFRD as a function of redshift. Cucciati et al. (2012) explore the contribution of galaxies of varying far-ultraviolet luminosities and find that fainter galaxies contribute more significantly at lower redshifts. This is the expected result from an overall downsizing in the star formation in galaxies, with more massive (and generally more luminous) galaxies forming their stars earlier, as first noted by Cowie et al. (1996). The consequence of this effect in our analysis is that the contribution of the least massive and faintest galaxies that are not included in our z ∼ 0 sample varies as a function of redshift, and hence the nominal shape of the cosmic SFRD as reported in Cucciati et al. (2012) is not completely appropriate to the particular set of z ∼ 0 galaxies we are considering here.

However, we might reasonably expect that any subset of galaxies drawn from the overall population would itself have a cosmic SFRD that is also well described by a log-normal function in time. Indeed the SFRDs in Cucciati et al. (2012) decomposed by luminosity appear to make a homologous set. Figure 5 demonstrates a model very similar to Figure 4, except that we have allowed for an additional log-normal component designed to mimic the contribution of the faintest galaxies to the cosmic SFRD. We include this as a component whose mass is set at 5 × 10⁹ M_☉ (a factor of two below the lower mass threshold of the data) and with an sSFR set from the mean mass–sSFR relation apparent in Figure 3. The total weight of this component (i.e., equivalent to setting the number of such galaxies) is such that it produces 20% of the total star formation in the SRFD diagram. This extra component—designed to mimic the contribution of the faintest and least massive galaxies—is fit as part of the modeling and behaves exactly as expected; i.e., interpreted as individual galaxies these are somewhat later-forming systems with a peak of star formation at z = 1.0, which contributes most significantly to the SFRD at z = 0. This component is shown in Figure 5. The addition of this component resolves the tension in the model noted above.

Zoom In Zoom Out Reset image size

Finally, note that measured precisely from the model in Figure 5 at the weighted mean redshift of the input sSFR sample, this component produces 43% of the star formation at that redshift. Figure 6 shows the mass-sorted cumulative star formation in the input data, as well as the nominal asymptote that this additional component suggests. This extra amount of star formation appears perfectly reasonable. Note that this agreement is at least somewhat arbitrary, as our choice of weight for this extra component does affect (and was informed by) Figure 6.

Zoom In Zoom Out Reset image size

That the realization of a galaxy population described above matches both the SFRD as well as the z ∼ 0 sSFR distribution does not in itself clearly validate the choice of a log-normal SFH for individual galaxies. For example, one could in principle assemble the SFRD from a large ensemble of galaxies and from a large range of functional forms (simple tophat or Gaussian SFHs, for example), and one can imagine that at least some of these may very well reproduce both the SFRD and the z ∼ 0 sSFR distributions.

To provide some insight on this point, we next compare the results of the above log-normal realization to the data from Paper III on galaxy populations at intermediate redshifts; specifically we compare the predicted sSFR distributions to the measured values. In this first analysis we simply compare the predictions of the model to that data; a model using the higher-z data as a constraint is presented in Section 3.2. Data are drawn from two sources, as described in Paper III: the field sample from the IMACS Cluster Building Survey (ICBS), and the All-wavelength Extended Groth Strip International Survey (AEGIS; Davis et al. 2007). SFRs in the ICBS are computed from a hierarchy of methods calibrated against 24 μm fluxes, with Hα used where available, then Hβ, or [O ii] λ3727 when not. 24 μm fluxes are used for the brightest such objects (see Paper I and III for further details). As described in Paper III, the SFRs in the AEGIS sample have been recalibrated to our measurements and our calibration by reconsidering their 24 μm fluxes and also measuring and calibrating emission lines in a set of DEEP2 spectra (upon which the AEGIS sample is based) using our own tools. All SFRs at all redshift are thus measured, as best possible, on the same calibrated system.

We compare the model predictions to data in four redshift bins: 0.2 < z < 0.4 and 0.4 < z < 0.6 from the ICBS field galaxy sample, and 0.6 < z < 0.8 and 0.8 < z < 1.0 from AEGIS. The mass limits, the approximate sSFR limits for 100% completeness, the total number of galaxies, and the fraction above the sSFR 100% completeness limit in each redshift bin are given in Table 1. Masses for the ICBS samples are as described in Paper I, computed in a manner akin to Bell & de Jong (2001) but with delayed exponential SFHs as the underlying star formation model. Masses for the AEGIS subsamples are taken from Davis et al. (2007) but converted to a Saltpeter initial mass function to match the other samples.

Before comparing the log-normal modeled sSFR distributions to the galaxy data two further corrections are necessary. First, note that the modeling as constructed does not explicitly account for merging; the individual SFHs tagged to individual galaxies at z ∼ 0 must be unmerged as one goes to higher redshift, since local galaxies are built up, to some extent, from smaller components over cosmic time. We will compare the model realization to the actual data by integrating the model over mass down to the appropriate mass limit, so a simple treatment of the merger history of the model galaxies is sufficient. Specifically we wish to know whether any individual galaxy in the z ∼ 0 sample remains above the mass limit in a higher redshift sample, or has "un-merged" and one or more of its un-merged components have been removed from the higher redshift sample by virtue of now falling individually below the mass limit. We draw a description of the mass-dependent rate for major mergers from Xu et al. (2012) and from this compute the merger rate over redshift for each z ∼ 0 galaxy from the present epoch back to a redshift assigned at random from the redshifts measured in the higher redshift bin. There are no strong trends in mass or sSFR versus redshift in these higher redshift datasets, so this simple process should adequately model the higher redshift data, and intra-bin variations in merger rates and sSFRs, without overly complicating this computation of a best-fit realization. Using this merger rate from Xu et al. we construct many realizations of the merger history of each galaxy, allowing for mergers, if they occur, with mass ratios up to a 3:1 split, and allowing mergers up to two steps deep in the merger tree for each z ∼ 0 galaxy. From this ensemble of mock histories we compute the likelihood that each z ∼ 0 galaxy (or a portion thereof) remains in the sample in the higher redshift bin when cut at the appropriate mass limit, and use these likelihoods as weights when computing the final sSFR distributions for the model. Note the correction due to accounting for merging in this way is not dramatic; in the highest redshift bin the variation in weights across the entire mass range is less than a factor of two.

The second correction is straightforward; when assessing a given SFH at higher redshifts (closer in time to the peak of the star formation in a given galaxy) somewhat less stellar mass has been lost. We adjust the galaxy masses and sSFRs appropriately.

Figure 7 shows the comparison of the sSFR distributions measured from the data reported in Table 1 and for the model realization shown in Figure 5. Modeled galaxies are assessed at the same redshift in the higher redshift bin used to compute the merger histories. The agreement is remarkably good. That the overall scaling matches—i.e., that the model tracks the movement of the sSFR distributions to higher values at earlier times—is not in itself surprising, as this is implicit in the fact that the model matches the cosmic SFRD. However, the fact that the shapes of the distributions at sSFRs higher than the 100% completeness limit agree rather well with the measured distributions—in all four redshift bins, comprising some 60% of the total timeline of star formation in the universe—is striking. The agreement is not perfect—for example in the highest redshift bin the model somewhat overpredicts the abundance of the highest sSFR galaxies—but nevertheless it suggests quite strongly that the choice of a log-normal SFH for each galaxy is well motivated. This agreement need not have happened; it is not explicit in the modeling process, as we will explore further in Section 5 below.

Zoom In Zoom Out Reset image size

The next step in our analysis is to create a model realization that uses the higher redshift sSFR data as an explicit constraint, in addition to the z ∼ 0 sSFRs. In this realization we retain the mass distribution given by the z ∼ 0 galaxies. An explicit galaxy-by-galaxy comparison to the higher redshift data is not possible, since the different datasets described in Paper III sample different co-moving volumes to different mass and sSFR limits in each redshift bin, and moreover a methodology to account in detail for the merging of galaxies across redshift bins—beyond the simple calculation sketched above—is not apparent. Thus any constraint must, to some extent, integrate across mass to ameliorate differences in the galaxy mass functions between the various datasets, and then compare distributions of sSFR values or related quantities, rather than making a galaxy-by-galaxy comparison. Indeed the comparison of model outputs to sSFR distributions in Section 3.1 does exactly that, integrating across mass directly, down to the mass limit of each redshift bin, and then simply comparing histograms of the logarithm of the sSFRs. This is likely not an optimal approach for constraining the model, however.

A closer look at the sSFR versus mass distribution for each redshift bin of the data, shown in Figure 8, illustrates the model-to-data comparison we have chosen to impose as a constraint on the model. Specifically, note that the sSFR versus mass values form a linear sequence in log–log space: this is the star-formation main sequence described by many authors at many redshifts (see Rodighiero et al. 2011 for a recent summary). This sequence appears to have an approximately constant width with mass. Over the redshift interval considered here, there is also no strong evidence for an evolution in the slope of this linear sequence; the zeropoint clearly evolves to higher values of sSFR at higher redshift but the slope remains approximately constant. To construct a distribution of sSFR-like values to constrain the model we thus consider the cumulative distribution of the difference between log(sSFR) at various redshifts and this relation at z = 0, over the appropriate mass range in each redshift bin. The resulting cumulative distributions for the data are given in Figure 8.

Zoom In Zoom Out Reset image size

As in Section 3.1 above, the effect of mergers is computed as a weight for each z ∼ 0 galaxy, and appropriately used when computing the cumulative distribution of residual sSFRs in the model. This final realization of the model, now additionally using the sSFR data in the higher redshift bins as a constraint, is shown in Figure 9.

Zoom In Zoom Out Reset image size

We remind the reader that all galaxies in these final model realizations have only smooth SFHs described by log-normal functions. No starburst activity is included. Starbursts might induce rapid—albeit temporary—changes in a given galaxy's sSFR, and we might expect to see this reflected as a tension between the sSFR distributions at low and high redshift, when attempting to connect the low and high redshift datasets using a smoothly varying SFH. However, the model in Figure 9 successfully incorporates all of the sSFR distributions to z = 1 without including any starbursts, a point discussed at length in the next section.