HERSCHEL-ATLAS GALAXY COUNTS AND HIGH-REDSHIFT LUMINOSITY FUNCTIONS: THE FORMATION OF MASSIVE EARLY-TYPE GALAXIES^*

The star formation history in galaxies is one of the key issues we have to understand in order to reconstruct how the universe evolved from small matter perturbations at the recombination epoch to the present richness of structures.

Star formation proceeds at a different pace for different galaxies, depending on the physical conditions of the available gas. Early-type galaxies (ETGs) and massive bulges of Sa galaxies are composed of relatively old stellar populations with mass-weighted ages of ≳ 8–9 Gyr (corresponding to formation redshifts z ≳ 1–1.5), while the disk components of spiral and irregular galaxies are characterized by significantly younger stellar populations. For instance, the luminosity-weighted age for most of Sb or later-type spirals is ≲ 7 Gyr (cf. Bernardi et al. 2010, their Figure 10), corresponding to a formation redshift z ≲ 1. In general, the old stellar populations feature low specific angular momentum as opposed to the larger specific angular momentum of the younger ones.

How can these facts be interpreted in the framework of the hierarchical evolution of the dark matter (DM) galaxy halos that have proven to be remarkably successful in accounting for the power spectrum (or the spatial correlation function) of the large-scale matter distribution (e.g., Springel et al. 2006)? A widely held view is that the merging of halos is also the principal mechanism driving the evolution of the visible part of galaxies (see Benson 2010 for a recent review). However, several bodies of evidence are difficult to reconcile with this scenario.

• 1.  
  ETGs are characterized by old and homogeneous stellar populations. Correlations tight enough to allow little room for random processes such as a sequence of mergers (apart from small mass additions through minor mergers at late epochs; see Kaviraj et al. 2008) and sensitivity to the environment (color−luminosity; fundamental plane relations; dynamical mass−luminosity) have been known for a long time and have been recently confirmed with very large samples, and shown to persist up to substantial redshifts (Renzini 2006; Clemens et al. 2009; Thomas et al. 2010; Rogers et al. 2010; Peebles & Nusser 2010, and references therein). More recently, a remarkably tight luminosity−size correlation has been reported (Nair et al. 2010). In addition, ETGs were found to host supermassive black holes whose mass is proportional to the bulge and to the halo mass of the host galaxy (see Magorrian et al. 1998; also Ferrarese & Ford 2005 for a review). All that indicates that the formation and evolution of ETGs is almost independent of environment, and driven mainly by self-regulation processes and intrinsic galaxy properties such as mass.
• 2.  
  There are rather tight, albeit not inescapable, observational constraints on the star formation timescale in the most massive ETGs. An upper limit comes from the observed α-enhancement or, more properly, iron underabundance compared to α elements. Depending on the slope of the assumed initial mass function (IMF), the observed α/Fe element ratios require star formation timescales ≲ 10⁹ yr (e.g., Matteucci 1994; Thomas et al. 1999). But in merger-driven galaxy formation models star formation in ellipticals typically does not truncate after 1 Gyr (Thomas & Kauffmann 1999; however, see Arrigoni et al. 2010; Khochfar & Silk 2011). If a standard IMF is assumed a lower limit comes from (sub)millimeter counts, implying that several percent of massive galaxies are forming stars at rates of thousands M_☉ yr⁻¹ at z ∼ 2–3 (see Chapman et al. 2003, 2005). This requires that this star formation rate (SFR) is sustained for ≳ 0.5 Gyr, much longer than the timescale of a merger-induced starburst, which is of the order of the dynamical time (∼0.1 Gyr for the massive ETGs of interest here; see, e.g., Benson 2010; Hopkins 2011). In other words, a single starburst episode is too short to account for the space density of (sub)millimeter bright galaxies as well as for their present-day stellar masses. And indeed models envisage a sequence of mergers, and associated starbursts, throughout the galaxy lifetime, i.e., over a time span much longer than the upper limit set by the α-enhancement. The problem of accounting for the counts of submillimeter galaxies may be eased assuming a top-heavy IMF (Baugh et al. 2005). Indeed, recent observational evidences (Gunawardhana et al. 2011; Dunne et al. 2011) indicate that highly star-forming galaxies have IMFs dN/dm∝m^−x with flatter high-mass power-law slopes x than galaxies with low SFRs. Gunawardhana et al. (2011), using a sample of galaxies from the Galaxy And Mass Assembly survey (Driver et al. 2009) covering the redshift range 0 < z < 0.35, find a dependence of x on the SFR that, extrapolated to an SFR = 1000 M_☉ yr⁻¹, would give x ≈ 1.5, substantially flatter than the Salpeter (1955) slope (x_S = 2.35). However, it is not clear that these results apply to the high-z proto-spheroidal galaxies considered in this paper, since the IMF may depend on other parameters, such as age and metallicity. Moreover, the most recent study of the evolution of galaxies in the far-infrared/submillimeter based on starbursts triggered by mergers (Lacey et al. 2010) resorts to an even flatter high-mass IMF (x = 1) and has still serious problems with reproducing the Herschel counts (see Section 7).
• 3.  
  Integral-field near-IR spectroscopy of galaxies with less extreme SFRs (of few 10² M_☉ yr⁻¹) at z ∼ 2, that appear to be very productive star formers in the universe (Dekel et al. 2009), has shown that in many cases they have ordered, rotating velocity fields with no kinematic evidence for ongoing merging (Genzel et al. 2006; Förster-Schreiber et al. 2009). Still, they harbor several star-forming clumps: a complex morphology is not necessarily a symptom of merging. These galaxies show tight SFR−mass correlations, with small dispersions (Daddi et al. 2007; Pannella et al. 2009; Dunne et al. 2009a; Rodighiero et al. 2010; Maraston et al. 2010). This is not easily reconciled with a scenario in which star formation proceeds through a series of short starbursts interleaved by long periods of reduced activity and these galaxies have been caught in a special, starburst moment of their existence. The data are more easily accounted for if the high SFRs are sustained for some 1–2 Gyr, much longer than a dynamical time typical of starbursts. Although the duration of the star formation phase for these objects is longer than that of the more extreme objects considered in the previous bullet by a factor of 2–3, their final stellar mass is several times lower because the SFR is about an order of magnitude smaller. For these galaxies a weak α-enhancement is predicted, consistent with observations.
• 4.  
  A comparison of the stellar mass functions at z ≳ 1.5 with the local one shows that little additional growth can have occurred for z ≲ 1.5 through minor mergers (Mancone et al. 2010; Fan et al. 2010; Kaviraj et al. 2009).

It is clear from the above that the reconstruction of the star formation history through cosmic time for galaxies of different masses provides a critical test for galaxy formation and evolution. Since the star formation occurs within dusty molecular clouds and is deeply obscured at ultraviolet and optical wavelengths, data at far-IR/(sub)millimeter wavelengths, where the absorbed radiation is re-emitted, are essential to provide a complete picture of it.

In this paper we focus on high-redshift (z ≳ 1) galaxies. We exploit the far-IR/submillimeter data collected by the Herschel Space Observatory (Pilbratt et al. 2010) during the Science Demonstration Phase (SDP) of the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010b). The H-ATLAS is an open-time key program that will survey ∼550 deg² with the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) and the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) in five bands, from 100 to 500 μm.

The observed SDP field covers an area of ≈3.8 × 3.8 deg² centered on (α, δ) ≈ (09^h 05^m, +0° 30'). Complete descriptions on reduction of PACS and SPIRE data are given in Ibar et al. (2010) and Pascale et al. (2011), respectively. Source extraction and flux density estimation are described in Rigby et al. (2011). The 5σ detection limits, including confusion noise, are 33.5, 37.7, and 44.0 mJy beam⁻¹ in the SPIRE bands at 250, 350, and 500 μm, respectively; in the PACS bands they are 132 mJy and 121 mJy at 100 and 160 μm, respectively.

The plan of the paper is the following. The selection of the sample is described in Section 2. In Section 3, we discuss the far-IR spectral energy distribution (SED) of high-z star-forming galaxies, a fundamental ingredient for our photometric redshift estimates, presented in Section 4. In Section 5, we estimate the galaxy luminosity functions (LFs) at different redshift in the range 1–4. In Section 6, we discuss some clues on star formation timescales for massive galaxies. In Section 7, model predictions for source counts from 250 μm to 2 mm are compared with the data. Our main conclusions are summarized in Section 8.

Throughout the work we adopt a standard, flat ΛCDM cosmology (see Komatsu et al. 2011) with matter density parameter Ω_M = 0.27 and Hubble constant H₀ = 70 km s⁻¹ Mpc⁻¹. We adopt a Bardeen et al. (1986) cold DM power spectrum with primordial index n_s = 1 and cosmic mass variance σ₈ = 0.81. Stellar masses and luminosities of galaxies are evaluated assuming the Chabrier's (2003) IMF; these can be converted to a Salpeter (1955) IMF on multiplying by a factor ≈1.6.

The H-ATLAS sources comprise both a low-z galaxy population, identified through matching to the Sloan Digital Sky Survey (SDSS; York et al. 2000) data (Smith et al. 2011a), and a high-z population (median redshift ∼2) identified through their far-IR colors (Amblard et al. 2010). Low-z galaxies are generally normal/star-forming late-type galaxies with moderate opacity (Smith et al. 2011b; Dunne et al. 2011). Through analyses of clustering these two populations are found to be very different; the low-z population (z < 0.3) clusters like star-forming blue galaxies (Guo et al. 2011; van Kampen et al. 2011; Maddox et al. 2010), while the high-z population clusters much more strongly, suggesting that the high-z sources reside in more massive halos (Maddox et al. 2010).

In this work, we investigate the evolution with cosmic time of high-redshift galaxies with intense star formation activity, interpreted as massive proto-spheroidal galaxies in the process of forming most of their stellar mass (Granato et al. 2004; Lapi et al. 2006). Since these objects are observed to be in passive evolution at z ≲ 1–1.5, we confine ourselves to z > 1. At these redshifts, sources above the H-ATLAS detection limits have dust luminosities ≳ 10¹² L_☉ and SFRs ≳ 100 M_☉ yr⁻¹. They are therefore ultraluminous infrared galaxies (ULIRGs). Their dust heating mostly comes from young massive stars within molecular clouds, implying, on one side, that their far-IR SEDs are generally (albeit not always, e.g., Hwang et al. 2010) warmer than those of low-z dusty galaxies which have higher contributions from cooler interstellar dust heated by old stars, and, on the other side, that their optical emission is strongly attenuated. The last point means not only that spectroscopic redshifts are available just for a tiny fraction of sources, but also that we do not have at our disposal the multi-frequency optical/near-IR photometry that allowed photometric redshift estimates at lower-z (Dye et al. 2010; Vaccari et al. 2010; Eales et al. 2010a).

In fact, for most H-ATLAS z ≳ 1 galaxies, the only available data is the Herschel photometry, primarily in the three SPIRE bands (250, 350, and 500 μm), plus mostly upper limits in the PACS 100 and 160 μm bands. A key issue is then whether the redshift estimates that can be obtained from such data are sufficient to obtain meaningful estimates of the LFs at least over a limited redshift range. At first sight one would be inclined to answer "no," but a closer investigation can suggest a more optimistic conclusion.

The dust re-radiation in starburst galaxies is expected to come from at least three astrophysical settings (e.g., Silva et al. 1998): molecular clouds, diffuse low-density clouds (cirrus), and circumnuclear regions, heated by active galactic nuclei (AGNs). The AGN dust emission peaks in the mid-IR (Granato & Danese 1994; Andreani et al. 2010; Lutz et al. 2010; Hatziminaoglou et al. 2010) and can be safely ignored in the SPIRE wavelength range. Molecular clouds are the preferential sites of star formation implying that they are endowed with intense radiation fields and relatively warm dust temperatures. The cirrus component is exposed to the less intense general radiation field due to older stellar populations that have come out from their native molecular clouds and is therefore characterized by lower dust temperatures.

In the nearby universe, molecular clouds and cirrus give comparable contributions to the far-IR emission from "normal" late-type galaxies, with relatively low SFRs (Rowan-Robinson et al. 2005). The colder cirrus contribution is especially important in less optically obscured IR galaxies (Hwang et al. 2010) while the warmer molecular cloud (sometimes referred to as "starburst") component becomes increasingly important for higher and higher SFRs (Rowan-Robinson et al. 2010). This argument also highlights a possible degeneracy: a "cold" observed SED may be associated either to a low-z cirrus-dominated galaxy or to a redshifted warm galaxy. If the redshift is estimated using a warm SED, cold low-z galaxies would be erroneously assigned high redshifts. This problem can be overcome, however, because cold, low-z galaxies are only moderately obscured by dust (the cirrus optical depth cannot be very large), and are therefore relatively bright in the optical bands. This is illustrated by Figure 1, which shows that the SEDs of optically identified z < 0.5 SDP galaxies studied by Smith et al. (2011b) imply r-band magnitudes brighter than the SDSS DR7 limit (r = 22.4) at all redshifts even if their 250 μm flux density is at the detection limit, while z > 1 ULIRGs, with SED like that of Arp220 (and even more, younger high-z galaxies with SED like that of SMM J2135-0102, "The Cosmic Eyelash"; Ivison et al. 2010b, Swinbank et al. 2010) are fainter than this magnitude limit for the same 250 μm flux density.²¹

Zoom In Zoom Out Reset image size

Therefore, we may weed out cold low-z galaxies by dropping SDP galaxies with SDSS counterparts (Smith et al. 2011a), except those with optical (spectroscopic or photometric) redshifts in the range of interest here (z > 1). This operation, however, has potential drawbacks. First, it leaves in low-z ULIRGs. This is not a big problem, since these objects have warm SEDs and therefore, as discussed below, their redshifts can be estimated with sufficient accuracy. Second, the reliability of SDSS counterparts can never be 100%. For example, strongly lensed galaxies generally have an apparently reliable counterpart which is most likely the foreground lens. These objects can however be recovered since most frequently the lenses are ellipticals, whose optical colors are incompatible with a large dust emission. The fraction of false identifications among the optical counterparts to H-ATLAS SDP sources with reliability R > 0.8 is estimated to be ≈5.8% (Smith et al. 2011a). Although this fraction is reassuringly small, we need to keep in mind that some truly high-z sources can be missed by our procedure.

Most importantly, not all the true r < 22.4 counterparts can be identified with >80% confidence, due to incompleteness of the SDSS catalog, positional uncertainties, close secondaries, and the random probability of finding a background source within that search radius (Dunne et al. 2011). According to Smith et al. (2011a), about 60% of the 6621 with 250 μm flux density >32 mJy have counterparts brighter than r = 22.4 mag in the SDSS (and are therefore, with few exceptions, at z < 1). Of these, 2423 could be identified with a reliability R > 0.8, implying that another ≈1550 sources, i.e., about 23% of the total sample are really at z < 1 but are missing a reliable identification. From the redshift distribution of reliable identifications (Figure 6 of Smith et al. 2011a) we estimate that about 20% of these sources are at z > 0.5. If the same proportion applies to unidentified sources, the identification incompleteness of z < 0.5 H-ATLAS SDP galaxies is ≈19%, which we conservatively round to 20%. The effect of this incompleteness on our LF estimates is discussed in Section 5.

Starting from the catalog of SDP sources by Rigby et al. (2011), which contains 6876 sources, we drop the galaxies for which Smith et al. (2011a) have identified reliable counterparts. We further require detection at ⩾3σ at 350 μm. In this way we get a sample (that will be taken as our reference sample) defined by the following criteria: (1) S_250 μm ⩾ 35 mJy; (2) no optical identification with R > 0.8 (Smith et al. 2011a); and (3) detection at ⩾3σ at 350 μm. The resulting sample is made of 3469 sources. The redshift estimates presented in Section 4 indicate that some of these (376) are z < 1 ULIRGS with SEDs akin to our templates; they will be excluded from the subsequent analysis.

Most (2763, i.e., ≈80%) galaxies satisfying the above criteria are detected at ⩾4σ at 350 μm, and this obviously helps with the photometric redshift estimates. This is because low-z galaxies have low $S_{350\,\mu \rm m}/S_{250\,\mu \rm m}$ ratios and are therefore lost when we raise the 350 μm flux limit (see also Amblard et al. 2010).

For the objects of interest here, with SFR ≳ 100 M_☉ yr⁻¹, only the "warm" (starburst) component is relevant. It is important to take into account, however, that its SED is much broader than a single-temperature graybody (see, e.g., Silva et al. 1998). Over the 50–500 μm range (in the rest frame) such an SED can be approximated, to better than 10%–20%, by a sum of two graybodies, each described by

$$\begin{equation} S_{\nu }\propto \frac{\nu ^{3+\beta }}{e^ {h\nu /kT_d}-1}, \end{equation} \tag{ 1 }$$

with temperatures T_d ≈ 30 K and ≈60 K, and dust emissivity indices β = 1.7 and 2, respectively (see also Dunne & Eales 2001). A fit to the observed SEDs of standard template starburst galaxies with quite different SFRs (M82, Arp220), and of the z ≈ 2.3 strongly lensed galaxy SMM J2135-0102, confirms the validity of this approximation and shows that the relative normalization of the two graybodies varies by factors of several. Mid/far-IR data on starburst galaxies emphasize the warmer component, while submillimeter data emphasize the colder one; if a single-temperature dust is used, data in different wavelength ranges may yield substantially different temperatures. If the source redshift is in the range 1 ≲ z ≲ 3.5, PACS and SPIRE wavelengths extend across the dust emission peak (which is typically at a rest-frame wavelength λ ≈ 90–100 μm). This may allow reasonably accurate redshift estimates using only Herschel data (Negrello et al. 2010). However, only a minor fraction of the H-ATLAS SDP galaxies have at least one PACS detection, usually at 160 μm, and these are mostly at low redshifts; we must also beware of the flux boosting by confusion, increasing at longer wavelengths (see Section 4).

Arp220 and SMM J2135-0102 are of particular interest because their SEDs are well determined and have SFRs quite typical of the galaxies considered here. We have modeled their SEDs from extreme-UV to radio frequencies through the spectrophotometric code GRASIL,²² which includes a sophisticated treatment of dust reprocessing (Silva et al. 1998, 2011; Schurer et al. 2009). The results are illustrated in Figure 2; it must be stressed that the dominant contribution for wavelengths λ ≳ 30 μm is provided by molecular clouds (more details in A. Bressan et al. 2011, in preparation). Assuming a Chabrier (2003) IMF, we obtain an SFR of ≈140 M_☉ yr⁻¹ for Arp220 and of ≈270 M_☉ yr⁻¹ for SMM J2135-0102 (a Salpeter IMF would yield ≈224 M_☉ yr⁻¹ and ≈430 M_☉ yr⁻¹, respectively). Then assuming the usual linear scaling between the SFR and the continuum far-IR luminosity L_FIR integrated between 8 μm and 1000 μm, i.e.,

$$\begin{equation} {\mathrm{SFR}\over M_{\odot }\,{\rm yr}^{-1}}=k\times 10^{-44}\,{L_{\rm FIR}\over {\rm erg\, s}^{-1}}, \end{equation} \tag{ 2 }$$

we find k_Arp = 2.5 and k_SMM = 3.4. The difference between the two coefficients is well within the expected range. Kennicutt (1998) pointed out that variations of ≈30% can be due to differences in the star formation history, implying different effective ages of the stellar populations. In Figure 2, we also sketch the far-IR SED of G15.141 (H-ATLAS J142413.9+022304), a strongly lensed submillimeter galaxy at z ≈ 4.24 with estimated SFR of several hundreds M_☉ yr⁻¹ (Cox et al. 2011), modeled as the sum of two graybodies with T₁ = 32 K, T₂ = 60 K, β = 2, and a ratio of 0.02 between the coefficients of the warm and of the cold components. We choose the SEDs of these three galaxies (Arp220, SMM J2135-0102, and G15.141) as our references to quantify the effect of different choices for the SED on estimates of the LFs. SMM J2135-0102 has the smallest $S_{60\,\mu \rm m}/S_{100\,\mu \rm m}$ ratio and G15.141 exhibits the steepest decrease at long wavelengths.

Zoom In Zoom Out Reset image size

A major source of concern is the boosting of fluxes in the Rigby et al. (2011) catalog due to confusion by faint sources. Since the effect is larger at longer wavelengths, because of the poorer angular resolution, it tends to bias the high redshift estimates based on Herschel photometry. According to the simulations carried out by Rigby et al. (2011), 56.5% of sources detected at ⩾5σ at 500 μm show a flux boosting by a factor >1.5 and 27.3% by a factor >2, and the effect becomes negligible for ⩾10σ sources.

To investigate the quantitative impact of the flux boosting, we estimated the photometric redshifts of 39 H-ATLAS galaxies at z > 0.5 for which spectroscopic redshifts are available. The results are presented in Figure 3 for our three reference SEDs. The average SED of low-z H-ATLAS SDP galaxies, determined by Smith et al. (2011b), was also used, for comparison. There is no indication that photometric redshifts of the high-z sources are systematically overestimated when we use the SED of SMM J2135-0102 as a template. The median value of Δz/(1 + z) ≡ (z_phot − z_spec)/(1 + z_spec) is 0.01 with a dispersion of 0.21. The situation is only slightly worse in the case of Arp220: the median value of Δz/(1 + z) is 0.06 with a dispersion of 0.27. The median offset between photometric and spectroscopic redshifts increases to 0.18, with a dispersion of 0.28, if we use the cooler SED of G15.141. In the redshift range (z ⩾ 1) we are most interested in that there are 24 objects. The median values of Δz/(1 + z) are −0.002 (dispersion 0.12) for SMM J2135-0102, 0.07 (dispersion 0.16) for Arp220, and 0.16 (dispersion 0.13) for G15.141. There is no statistically significant difference between the 6 strongly lensed objects, having >10σ detections at 500 μm, and the other 33 sources, whose flux densities are representative of those in our sample. The median values of Δz/(1 + z) and the dispersions rapidly increase as we go down in redshift, as expected since our templates do not match those of low-z galaxies whose SEDs include large cirrus contributions. For the 1096 H-ATLAS SDP sources with z < 0.5 and spectroscopic redshift in the Smith et al. (2011a) catalog, we find median Δz/(1 + z) of 0.38 (dispersion 0.45), 0.42 (dispersion 0.54), and 0.63 (dispersion 0.55) for SMM J2135-0102, Arp220, and G15.141, respectively. Not surprisingly, the redshift estimates based on the mean low-z SED of Smith et al. (2011a) go in the opposite direction: the redshifts are systematically underestimated. The mean offset is small (Δz/(1 + z) = −0.08) at z < 0.5 and increases (in absolute value) to −0.21 for z > 0.5 and to −0.34 for z > 1. The dispersions are 0.31, 0.29, and 0.23 for z < 0.5, z > 0.5, and z > 1, respectively. This confirms that the low-z, optically bright, H-ATLAS galaxy population has far-IR properties different from those of high-z galaxies. Hence, the SEDs determined at low-z are not applicable at high-z.

Zoom In Zoom Out Reset image size

This test suggests that, at least for sources at z ≳ 1, errors on photometric redshift estimates are more related to the choice of the SED template than to the signal-to-noise ratio at 500 μm, and hence to flux boosting. This is mostly due to the fact that the relative errors on flux densities are larger at longer wavelengths so that the data points more liable to flux boosting weight less in the minimum χ² fit to the template SED. We do see in most cases that the best-fit SED is below the 500 μm data, as expected if the latter is overestimated. We have checked this by lowering the 500 μm flux densities first by 20%, and this left the derived redshift distribution almost unchanged, and then by a factor of two. In the latter case, we got a median value of (z_deboost − z)/(1 + z) ≈ −0.1, with only 5% of cases at < − 0.2 (here z is the redshift estimated using the fluxes tabulated by Rigby et al. 2011).

A second concern is the effect of the SED variety of active star-forming galaxies. To quantify the corresponding errors in the redshift estimates, we generated a catalog of 9 × 10³ galaxies with a distribution of flux densities reflecting that of our sample and redshifts randomly selected in the range 1 ≲ z ≲ 3.5. Each redshift was randomly assigned to 1 of the 19 well observationally determined SEDs of local (ultra)luminous IR galaxies with SFRs ≳ 20 M_☉ yr⁻¹ and with contributions from an active nucleus to the far-IR flux f_AGN  ≲  10% studied by Vega et al. (2008). To the 250, 350, and 500 μm flux densities we associated errors extracted randomly from the distribution of the errors for real observations. We have then estimated the redshifts of simulated galaxies, based on flux densities at SPIRE wavelengths, using each of our three template SEDs.

In ≈95% of the cases we have |z_phot − z_true|/(1 + z) ≲ 0.2. Note that this is likely an upper limit since the previous test and the analysis of multi-frequency source counts (see Section 7) strongly indicate that the SEDs of high-z galaxies with intense star formation are more uniform, and closer to that of SMM J2135-0102, than those of star-forming low-z galaxies. This can be expected since these galaxies have generally lower SFRs (of order of tens M_☉ yr⁻¹, to be compared with >100 M_☉ yr⁻¹ for high-z galaxies) and higher contributions to dust heating from old stellar populations. As the SFR increases, the ratio between the contributions to the far-IR/sub-mm SED from the molecular cloud component (associated to star formation) and from the cirrus component (heated by the general interstellar radiation field) increases causing a systematic variation of the SED shape with the specific SFR, as indeed observed (da Cunha et al. 2008; Smith et al. 2011b). This trend stops when the former component dominates, as in the case of sources considered in this paper, and the SEDs become more uniform.

We have estimated the redshifts²³ of objects in the samples defined in Section 2. Sources undetected at ≳ 3 σ at 500 μm were attributed a 3σ upper limit of 27 mJy. The calculations were repeated using 5σ upper limits (45 mJy): the derived redshift distribution did not change appreciably. The results are similar, but somewhat more sensitive to the effect of boosting (based on the tests described above), if we use the cataloged flux densities and errors also for sources with S_500 μm < 3σ. The redshift estimate is the result of a minimum χ² fit of each of the three templates to the SPIRE and PACS data (including 3σ upper limits²⁴). Figure 4 shows that, after correcting for the offsets highlighted by Figure 3, the derived redshift distributions are only moderately affected by the choice of the template SED; they have broad maxima in the range 1.5 ≲ z ≲ 2.5 and a tail extending up to z ≈ 3.5, consistent with earlier estimates for BLAST (Ivison et al. 2010a) and Herschel (Amblard et al. 2010; Eales et al. 2010a) samples, when the different selections, and in particular the fact that we have dropped galaxies with R > 0.8 optical counterparts, are taken into account. Sources with estimated z < 1 will be excluded from the subsequent analysis. With our preferred SED, that of SMM J2135-0102, our reference sample (S_250 μm > 35 mJy, S_350 μm > 3σ, and no R > 0.8 optical identifications) contains 3093 galaxies (i.e., ≈45% of the 6876 H-ATLAS SDP galaxies) at z > 1, consistent with the fraction of sources expected to be below the SDSS limit (∼40%) and therefore at high-z. Only 33 of them have ⩾5σ detections in at least one PACS channel; these include the strongly lensed galaxies found by Negrello et al. (2010).

Zoom In Zoom Out Reset image size

We have computed the LFs at rest-frame wavelengths of 100 and 250 μm in four redshift intervals (1.2 ⩽ z < 1.6, 1.6 ⩽ z < 2, 2 ⩽ z < 2.4, and 2.4 ⩽ z < 4), exploiting the classical Schmidt (1968) 1/V_max estimator, our redshift estimates, and K-corrections computed from the SEDs used for the redshift estimates. The V_max was computed, for each galaxy, taking into account the redshift boundaries of the bin and the maximum accessible redshifts implied by the 250 and 350 μm flux density limits. In the case of our preferred SED (SMM J2135-0102) the numbers of sources in these redshift bins are 900, 891, 821, and 502, respectively.

The error on z introduces both a statistical and a systematic effect. The latter is related to the Eddington bias, that can be quantified with a Bayesian approach. The probability that the true redshift of a source is z when its estimated value is z_e reads

$$\begin{equation} p(z|z_e) \propto P(z_e|z)\,n(z), \end{equation} \tag{ 3 }$$

where P(z_e|z) is the probability that the estimated redshift of a source is z_e when its true value is z and n(z) is the true redshift distribution. We take P(z_e|z) to be a Gaussian with mean z and dispersion equal to the rms difference between spectroscopic and photometric redshifts for the galaxies in Figure 3 and n(z_e) as a proxy for n(z). The maximum likelihood estimate of z is the value for which p(z|z_e) peaks. All the LFs presented in this paper are corrected for this bias. We caution, however, that this correction does not entirely account for the Eddington bias due to the errors on estimated luminosities. The effect is expected to be smaller than that of the error on z, which is the main source of the error on luminosity. A full treatment would require extensive numerical simulations that are beyond the scope of this paper.

As mentioned in Section 2, we expect a residual contamination of our reference sample by cold low-z galaxies with genuine r < 22.4 counterparts that could not be identified with R > 0.8. To quantify the effect of this contamination on our LF estimates we have exploited the H-ATLAS SDP sources with spectroscopic redshift in the Smith et al. (2011a) catalog. Out of a total of 1096 galaxies with z_spec < 0.5, 682 comply with our selection criteria (S_250 μm > 35 mJy and S_350 μm > 3σ). For these sources the use of the SMM J2135-0102 template for redshift estimates is inappropriate and yields a substantial positive offset and a rather broad dispersion of z_phot − z_spec (see Section 4). In fact, after the Bayesian correction described above, this template yields z_phot ⩾ 1.2 for 160 of them. Of these, 88 fall in the first redshift bin, 45 in the second, 21 in the third, and 6 in the fourth. If the sample of SDP galaxies with z_spec < 0.5 is representative and about 20% of the 5021 galaxies in our sample prior to the removal of R > 0.8 objects have unrecognized genuine SDSS counterparts, the low-z contaminants would amount to 14% of sources in the first redshift bin, 7.4% in the second, 3.8% in the third, and 1.8% in the fourth. On the other hand, it is estimated (Smith et al. 2011a) that ≈5.8% of the 2418 R > 0.8 identifications are spurious. If these spurious identifications have the same redshift distribution as the 3469 sources not identified with R > 0.8, we would be missing about 4% of sources in each redshift bin. Finally, we have to take into account the incompleteness of the SDP catalog. This can be done using the flux density dependent correction factors given in Table 2 of Rigby et al. (2011).

The estimated LFs, taking into account all the corrections described above, are presented in Figures 5–8 and in Tables 1 and 2. The upper scale in these figures displays the SFR corresponding to the 100 or 250 μm luminosity for the SMM J2135-0102 calibration giving

$$\begin{equation} {L_{100\,\mu \rm m}\over {\rm W\, Hz}^{-1}}=5.9\times 10^{23}\,{\mathrm{SFR}\over M_{\odot }\,{\rm yr}^{-1}}, \end{equation} \tag{ 4 }$$

while for L_250 μm the coefficient is 1.4 × 10²³. Since for galaxies with intense star formation the rest-frame dust emission peaks in the range λ ≈ 90–100 μm, the 100 μm luminosity is a good estimator of the SFR.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Note that the flattening of the LFs at the lowest luminosities may be, at least in part, due to the overestimate of the accessible volume yielded by the 1/V_max estimator for objects near the flux limit, as pointed out by Page & Carrera (2000; see also Eales et al. 2009). On the other hand, the Page & Carrera (2000) estimator holds under the assumption that the LF varies little within the luminosity bin, while in our case it is very steep over most of the luminosity range.

Figure 6 illustrates the effect of varying the template SED, after correcting for the median offsets in the redshift estimates highlighted by Figure 3. The differences are quite limited. Such a stability of the LF estimates follows from some favorable circumstances. First, the large numbers of sources in each redshift bin smooth out the effect of errors on redshift estimates. Second, for the redshift range considered here (1 ≲ z ≲ 4) one of the SPIRE wavelengths is always sampling directly the rest-frame SED in the range 100 μm ≲ λ ≲ 125 μm, implying that K-corrections (and the related uncertainties) are minimal. Third, the LFs are only moderately sensitive to the uncertainty on redshift estimates. This can be shown as follows. The luminosity at a rest-frame frequency ν is related to the flux density at the observed frequency ν_obs by

$$\begin{equation} L_{\nu }=f_{\nu _{\rm obs}}\,{4\pi \, D^2_L(z)\over (1+z)\times l}, \end{equation} \tag{ 5 }$$

where

$$\begin{equation} l={L_{\nu _{\rm obs}\, (1+z)}\over L_{\nu }}. \end{equation} \tag{ 6 }$$

Since, in our case, ν ≈ ν_obs (1 + z), l ≈ 1. Then it is easily checked that, in the redshift range of interest here, L_ν∝(1 + z)³, so that

$$\begin{equation} \Delta \log L\approx {3\over \ln (10)}\,\left({\Delta z\over 1+z}\right). \end{equation} \tag{ 7 }$$

The LF, ϕ, is computed by weighting each galaxy by the inverse of the maximum accessible volume (Schmidt 1968), i.e., ϕ∝1/V. At z ≳ 1 we have, roughly, V∝(1 + z)^2.8, whence

$$\begin{equation} \Delta \log \phi \approx -1.2\, \left({\Delta z\over 1+z}\right). \end{equation} \tag{ 8 }$$

Since Δz/(1 + z) ≈ 0.2, the error on the LF determination due to uncertainty in the redshift estimate amounts to around 0.25 dex. The overall uncertainty is the sum in quadrature of the statistical error and of the error coming from the uncertainty on the redshift estimate. As shown in Figure 8, the correction for the Eddington bias is rather small compared to the overall uncertainty of the LF.

In Figures 5 and 7, our LF estimates are compared with those by Gruppioni et al. (2010) and Eales et al. (2010a), respectively. The latter authors adopted a graybody SED with T = 26 K and a dust emissivity index β = 1.5. The assumed dust temperature is possibly more appropriate for relatively low SFR/far-IR luminosity galaxies while typical dust temperatures of SPIRE-detected z ≈ 2 galaxies are somewhat higher (Chapman et al. 2010; Amblard et al. 2010). The different choice for the SED has a substantial impact on the K-correction and therefore on the estimated rest-frame luminosity. For example, for a galaxy at z ≈ 1.5 and a given observed 250 μm flux density, the Eales et al. SED yields a rest-frame 250 μm luminosity a factor ≈1.5 higher than that obtained using either the SMM J2135-0102 or the Arp220 SED. Once we replace their K-correction with ours, the LF estimates by Eales et al. (2010a) are found to be in good agreement with ours (see Figure 7). The two estimates are to some extent complementary: the deeper samples used by Eales et al. have allowed them to reach lower luminosities, while the larger H-ATLAS SDP area has allowed us to reach higher luminosities at high redshifts.

Gruppioni et al. (2010) exploited the deep PACS data at the 100 and 160 μm in the GOODS-N field (∼150 arcmin²), obtained as part of the PACS Evolutionary Probe (PEP) survey, to estimate the 60 and 90 μm rest-frame LFs up to z ∼ 3. They have used all the available data to derive the SEDs of their sources.

About 31% of the Eales et al. sources have spectroscopic redshifts and almost all the others have photometric redshifts typically based on nine optical and near-IR bands. Gruppioni et al. (2010) have spectroscopic redshifts for ∼70% of their sources; for the remaining ∼30% they have obtained photometric redshifts from multi-frequency data. The agreement between our LFs and theirs is an additional confirmation that photometric redshifts derived from submillimeter photometry are good enough for the present purpose.

In the sampled luminosity range, the LFs exhibit an exponential falloff and a substantial luminosity evolution at least up to z ≈ 2.5, while a weaker evolution at higher-z is indicated (see Figure 8). The 250 μm luminosity corresponding to log (ϕ/Mpc⁻³) = −5 is log (L₂₅₀/W Hz⁻¹) = 26.02, 26.16, 26.32, and 26.38 for z ≈ 1.4, 1.8, 2.2, and 3.2, respectively. At 100 μm the corresponding luminosities are log (L₁₀₀/W Hz⁻¹) = 26.65, 26.78, 26.94, and 27.01. These results are consistent with those based on PEP survey data (Gruppioni et al. 2010), which have poorer statistics at high-z. We remark that our data supplemented with those by Gruppioni et al. (2010) and Eales et al. (2010a) allow the determination of the LFs over only one order of magnitude in luminosity for 1 ≲ z ≲ 2, and over only a factor of about three for z ≳ 2.

We now discuss how the LFs of high-redshift star-forming galaxies at submillimeter wavelengths and the counts in submillimeter and millimeter bands concur with their clustering properties in probing the process of star formation in the progenitors of massive ETGs.

6.1. Clustering Properties and Host Halo Masses of Submillimeter Galaxies

6.2. Timescales of Submillimeter and UV Bright Phases

The bright end of the LFs provides information on the average duty cycle of the star formation in the large halos indicated by the clustering analysis. This can be illustrated using a simple model (hereafter referred to as the "toy model") relying on the following assumptions.

• 1.  
  Most of the star formation occurs soon after the fast collapse phase of the halo, as identified by Zhao et al. (2003) and by several subsequent works (e.g., Diemand et al. 2007; Genel et al. 2010; Wang et al. 2011).
• 2.  
  The halo formation rate is given by the positive term of the derivative with respect to the cosmic time of the Sheth & Tormen (1999, 2002) cosmological mass function (Equation (1) of Lapi et al. (2006)). This was shown to be a good approximation at the redshifts (z ≳ 1) and for the halo masses (M_H ≳ 3 × 10¹² M_☉) of interest here (e.g., Haehnelt & Rees 1993; Sasaki 1994; see also the discussion by Lapi et al. 2006).
• 3.  
  During the main episode of star formation, the SFR is proportional to the halo mass, i.e., SFR∝M_H, with a small scatter.
• 4.  
  The duration Δt_sf of the main star formation episode, before quenching by the AGN feedback, is roughly constant for the considered range of halo masses (3 × 10¹² M_☉ ≲ M_H ≲ 3 × 10¹³ M_☉) and shorter than the expansion timescale at the source redshift.

This very simple model provides a good fit to the data (Figures 6 and 7) if Δt_sf ≈ 7 × 10⁸ yr and

$$\begin{equation} \mathrm{SFR}=35\left({M_{\rm H}\over 10^{12}\, M_{\odot }}\right)\left({1+z\over 2.5}\right)^{2.1}\, M_{\odot }\,{\rm yr}^{-1}, \end{equation} \tag{ 9 }$$

with a rather small (±35%) scatter. We have used Equation (2) with the calibration appropriate for SMM J2135-0102 to pass from SFR to FIR luminosity, and the SED of this object to go from the FIR luminosity to monochromatic luminosities. The fit would significantly worsen if the halo mass range associated with individual galaxies extends above M_H ≳ 3 × 10¹³ M_☉. This suggests that the star formation becomes very inefficient in the most massive halos, possibly due to the cooling time becoming longer than the expansion time. The bright end of the LFs is quite sensitive to the value of σ₈; e.g., raising σ₈ from our fiducial value of 0.81–0.9 would increase by a factor ≈2 the 250 μm LF at z ≈ 2.5 and L_250 μm ≈ 3 × 10²⁶ W Hz⁻¹. Interestingly, the duration of the main star formation episode yielded by the toy model matches that inferred from the observed α-enhancement of massive local ETGs (see Section 1).

A rough estimate of the mass in stars at the end of the main star formation episode, for given halo mass and redshift, can be obtained by multiplying the appropriate SFR by its duration Δt_sf ≈ 7 × 10⁸ yr, and correcting for the fraction of stellar mass returned to the interstellar medium (ISM). For a Chabrier IMF the latter amounts to ≈30% after 0.1 Gyr from a burst of star formation, and increases to ≈40% after 1 Gyr and to 50% after several Gyrs. Combining with the halo mass function at that redshift we get an estimate of the stellar mass function of ETGs. A comparison with observational determinations at different redshifts is shown in Figure 9. The agreement is remarkably good, considering the crudeness of the approach. This suggests that the toy model captures the basic aspects of the star formation in massive galaxy halos at high redshift.

Zoom In Zoom Out Reset image size

A comparison of our LFs with the UV ones at z ≈ 2–3 (Sawicki & Thompson 2005; Reddy & Steidel 2009) for galaxies with comparable SFRs shows that the UV space densities are a factor ≳ 100 lower, implying that the UV bright, dust-free phase is much shorter than the far-IR bright phase. The far steeper slope of the UV LF, compared with the submillimeter ones, for SFR ∼ hundreds M_☉ yr⁻¹ implies that, for the corresponding range of galaxy masses, the UV-bright phase is shorter for larger SFRs, i.e., for more massive galaxies, as predicted by Mao et al. (2007).

If, as implied by most current semi-analytical models (Granato et al. 2001, 2004; Croton et al. 2006; Hopkins et al. 2008; Somerville et al. 2008), the black hole growth is linked to star formation, most of the present-day black hole mass should have been accreted by the end of the star-forming phase (see Bonfield et al. 2011). Therefore, the local ratio M_•/M_⋆ ≈ 2 × 10⁻³ (Marconi & Hunt 2003) between the black hole mass M_• and the stellar mass M_⋆ of the host ETG should be already in place then (but see also McLure et al. 2006 and Peng et al. 2006). If so, the stellar mass function translates immediately into a black hole mass function which, in turn, can be translated into a bolometric luminosity function assuming that the black holes, immediately before being switched off, emit at the Eddington limit for a time t_vis. We get

$$\begin{equation} {L_{\rm BH,bol} \over L_\odot } \approx 79\, {M_\star \over M_\odot }. \end{equation} \tag{ 10 }$$

The corresponding B-band luminosity is L_B = 0.1(10/k_bol)L_{BH, bol}, where k_bol is the bolometric correction (Marconi et al. 2004; Hopkins et al. 2007). We can then compare the comoving density of galaxies with that of the associated QSOs. For example, at z ≈ 2 there are ≈2.7 × 10⁻⁵ submillimeter galaxies per Mpc³ with M_⋆ ≳ 2 × 10¹¹ M_☉, while the density of QSOs brighter than the corresponding luminosity of L_B ≈ 1.6 × 10¹² L_☉ (k_bol ≈ 10) is ≈8 × 10⁻⁷ Mpc⁻³ (Croom et al. 2004). This implies that the optical visibility time of the QSOs is about a factor of 30 shorter than the duration of the submillimeter bright phase, in agreement with independent estimates (Shankar et al. 2004; Marconi et al. 2004; see also Lapi et al. 2006).

The physical model by Granato et al. (2001, 2004), further elaborated by Lapi et al. (2006), complies with the indications coming from the previous analysis. In this model the star formation in massive galaxies occurs on a timescale of several 10⁸ yr, is very soon obscured by dust, and is stopped by quasar feedback. The star formation is triggered by the rapid cooling of the gas within a region with an approximate size ≈1/3–1/4 of the halo virial radius, i.e., 20–30 kpc, encompassing about 40% of the total mass (DM plus baryons), and is regulated by the energy feedback from supernovae (SNe) and AGNs (the latter being relevant especially in the most massive galaxies). As a result, the star formation is a very inefficient process, as proved by the fact that on the average only about 5%–10% of the baryons in the universe are converted into stars (Persic & Salucci 1992; Fukugita & Peebles 2004), with a possible maximum of 15%–20% for halo masses ∼10¹² M_☉ (Shankar et al. 2006; Moster et al. 2010). It is worth noticing that the Negrello et al. (2007) predictions for the lensed submillimeter bright galaxies follow directly from the LF and redshift distribution yielded by the Lapi et al. (2006) model. We have updated this model by replacing the Arp220 GRASIL SED used, e.g., by Negrello et al. (2007) with the SED of SMM J2135-0102. The resulting LFs are shown in Figures 5 and 7, labeled as "full model." We stress that the LFs from the full model are not fits to the data, but have been computed with the same parameters used by Lapi et al. (2006; see their Table 1). This model has proven to be successful in reproducing a wealth of observables, including high-redshift quasar LFs.

Note that, in the present framework, sources dominating the high-z LF are proto-spheroidal galaxies that are in passive evolution at z ≲ 1.5–1 and therefore do not contribute to the low-z FIR/submillimeter LFs computed, e.g., by Dye et al. (2010), Vaccari et al. (2010), and Dunne et al. (2011).

An alternative model whereby intense star formation at z ≈ 2–3 is not supported by mergers was proposed by Dekel et al. (2009). In this model, star formation is driven by steady cold streams supplying gas at an approximate rate of

$$\begin{equation} \dot{M}\approx 52 \left({M_{\rm H}\over 10^{12}\, M_{\odot }}\right)^{1.15}\,\left({1+z\over 2.5}\right)^{2.25} f_{0.165}\,M_{\odot } {\rm yr}^{-1}, \end{equation} \tag{ 11 }$$

where f_0.165 is the baryonic fraction in the halos in units of the cosmological value, f_b = 0.165. The dependence of the cold gas inflow rate, $\dot{M}$, on M_H and z are similar to those implied by our toy model (Equation (9)), implying that the SFR must be roughly proportional to $\dot{M}$. To account for the SFR indicated by the Herschel data, the star formation efficiency must be very high, ∼70% for M_H ≈ 10¹² M_☉ and z ≈ 1.5. For comparison, the mean cosmic star formation efficiency is Ω_star/Ω_baryon ≈ (5.8 ± 1.1)% (Fukugita & Peebles 2004; Komatsu et al. 2011). In halos more massive than ≈10¹³ M_☉ cold streams are suppressed by shock heating.

Figures 10–12 compare the predictions of the "full model" with the observed counts in Herschel/SPIRE bands. In these figures the contributions to the counts of proto-spheroidal galaxies have been complemented with those of normal late-type and starburst galaxies computed by Silva et al. (2005) and Negrello et al. (2007). As for massive proto-spheroidal galaxies, the main difference with the counts in the latter paper comes from having replaced the Arp220 SED yielded by GRASIL with that of SMM J2135-0102. The contribution from lensed proto-spheroidal galaxies has been estimated using the amplification distribution of Perrotta et al. (2003) and Negrello et al. (2007).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An interesting prediction of the model is that massive proto-spheroidal galaxies dominate the (sub)millimeter counts over a limited flux density range (about a decade). At 250 μm the Euclidean-normalized differential counts of these objects peak at ≈30 mJy, i.e., roughly at the detection limit of the H-ATLAS survey, and sink down rapidly at fainter fluxes where the contribution of starburst galaxies takes over and accounts for the results of the P(D) analysis by Glenn et al. (2010). Since the proto-spheroidal galaxies are mostly at z ≳ 1.5 while late-type/starburst galaxies are mostly at z ≲ 1.5, the model implies that the redshift distribution drifts toward lower redshifts as we go fainter. This does not contradict the finding of Bourne et al. (2011) that the contribution of z < 0.28 galaxies to the cosmic infrared background is <5% in the SPIRE bands.

The model accurately fits the source counts from 250 μm to 850 μm (Figures 10–13). The observed counts at 1.1 mm (Figure 14) span the peak of the Euclidean-normalized counts, while those at 1.4 and 2 mm (Figures 15 and 16) only cover the brightest tail of the counts, dominated by strongly lensed sources (apart from radio sources). The model somewhat overestimates the counts of strongly lensed galaxies at 1.4 and 2 mm, and the most recent counts at 1.1 mm. This may suggest that higher-z galaxies, that yield larger and larger contributions to the bright counts at increasing millimeter wavelengths, have SEDs slightly colder than SMM J2135-0102 and closer to that of G15.141.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The flux density range over which massive proto-spheroidal galaxies dominate the counts increases somewhat with increasing wavelength (Figures 10–16), reflecting the increase of the redshift range through which the strongly negative K-correction makes the flux corresponding to a given luminosity essentially independent of distance. The broadening of the peak in the Euclidean-normalized counts of proto-spheroidal galaxies is not very large, however, because massive halos become increasingly rare at high-z. At all wavelengths the bright counts of proto-spheroidal galaxies drop down very steeply, reflecting the exponential decline of the halo mass function. This rapid falloff, borne out by the data, can hardly be accounted for by phenomenological models that evolve the LF of local populations of dusty galaxies backward in time: spheroidal galaxies are essentially in passive evolution since z ≈ 1.5 and are therefore not represented in local LFs at far-IR to millimeter wavelengths.

The physical model of Lacey et al. (2010) works pretty well at 850 μm but does not correctly predict the rapid falloff of the counts above ≈30 mJy at 250 μm, 350 μm, and 500 μm (Figures 10–13). The decline of the counts is most easily understood if the duration of the most active star formation phase in massive halos is relatively short, as also required by the observed α-enhancement (see Section 1). However, in merger-driven evolutionary models, such as the one by Lacey et al. (2010), the star formation does not truncate after ≲ 1 Gyr. At 850 μm, consistency with the data can be recovered if submillimeter bright galaxies contain large amounts of cold dust, but this leads to 250 μm, 350 μm, and 500 μm counts less steep than observed.

The fact that using, for high-z galaxies, the SMM J2135-0102 SED we are able to simultaneously reproduce the source counts over a broad wavelength range indicates that this SED is indeed typical for submillimeter bright galaxies. In fact the results on the counts vary substantially if we use different SEDs (see Figure 17). Almost by construction, in all cases the 250 μm (the main selection wavelength) counts are accurately reproduced. But if, for example, we use the SED of G15.141 instead of the SMM J2135-0102 one, the longer wavelength counts are somewhat underestimated. The situation is better, but still not as good as for our preferred SED, if we use the Arp220 SED. The quite cold (T_d = 26 K) SED used by Eales et al. (2010a) leads to counts substantially in excess of the observed ones. Thus, the "colors" of the counts constrain the SEDs of high-redshift star forming and in particular on the abundance of galaxies with cold dust temperatures. As pointed out by Hwang et al. (2010) apparently very cold values of T_d could be caused by an overestimation of the submillimeter fluxes due to blending problems or inappropriate single-temperature SED fitting.

Zoom In Zoom Out Reset image size

Figure 18 compares the redshift distributions of z > 1 H-ATLAS SDP sources detected at S > 35, 40, 48 mJy at 250, 350, and 500 μm, respectively, with those predicted by the model. The effective flux correction factors tabulated by Clements et al. (2010) have been applied for this comparison. Since the redshift distributions are a key ingredient for the successful prediction of the abundance of strongly lensed galaxies and of their redshift range (Negrello et al. 2010), they cannot be badly wrong. Thus, the agreement between our estimates and model predictions is a confirmation of the global consistency of our results.

Zoom In Zoom Out Reset image size

The shift in the redshift distribution of submillimeter bright galaxies to higher redshifts with increasing selection wavelength happens because the strongly negative K-correction extends to larger and larger redshifts and the lensing probability increases with z. For example, we expect that, after removing local objects that are found to be only ≈1/3 of the total (Vieira et al. 2010), the redshift distribution of sources with S_1.4 ≳ 10 mJy has a broad peak around z ≈ 4.5. Using the SMM J2135-0102 SED we also find that the observed 1.4–2 mm spectral indices are correlated with z: for z ≲ 4.5, α²_1.4 ≳ 3.2, while for z ≳ 4.5, α²_1.4 ≲ 2.6. Thus, the value of this spectral index may be a rough redshift indicator.

In Figure 19, we illustrate the contribution of proto-spheroids to the extragalactic infrared background intensity (Lagache et al. 1999; Hauser & Dwek 2001) according to our physical model. At λ ≳ 850 μm the background is almost entirely accounted for by high-redshift proto-spheroids, while these contribute about 65% and about 50% of the background at λ ≈ 500 μm and λ ≈ 250 μm, respectively.

Zoom In Zoom Out Reset image size

We have exploited the H-ATLAS SDP survey data to investigate the evolution of the 100 and 250 μm LFs of bright star-forming galaxies (SFR ≳ 10² M_☉ yr⁻¹) at z ≳ 1. Redshifts have been estimated using three SED templates representative of the range of well-measured SEDs of galaxies with such a high star formation. The rms uncertainties on redshift estimates have been assessed comparing our redshift estimates with spectroscopic redshift measurements for 39 H-ATLAS galaxies at z > 0.5 as well as by means of simulations. Both methods yield rms values of Δz/(1 + z) ≈ 0.2 or smaller. The LFs derived using the redshift estimates based on each of the three SED templates are very similar to each other. The uncertainties due to the spread of redshift estimates are added in quadrature to Poisson errors to compute the global uncertainties on the LFs. Our LF estimates are in close agreement (in the common redshift and luminosity range, after applying the same K-corrections) with those at 250 μm by Eales et al. (2010a), based on a substantial fraction of spectroscopic redshifts complemented with photometric redshifts of optical/near-IR counterparts, as well with those at 90 μm by Gruppioni et al. (2010), based on an even higher fraction of spectroscopic redshifts. This is a further confirmation that our redshift estimates are sufficiently accurate to allow reliable estimates of the LFs.

The SED of SMM J2135-0102 was found to perform significantly better than the others in the tests we made, at least up to z ≈ 3, and was therefore adopted as our reference. Remarkably, this SED allowed us to simultaneously fit the counts over a broad wavelength range.

We find (see Figure 8) a significant luminosity evolution at least up to z ≈ 2.5 while the LF shows a modest variation between the last two redshift bins, centered at z ≈ 2.2 and z ≈ 3.2, respectively, even though the corresponding time interval is ≈1 Gyr, substantially larger than the time interval (≈0.65 Gyr) between the central redshifts of the second and third bin (z ≈ 1.8 and ≈2.2). This is consistent with the results based on PEP survey data (Gruppioni et al. 2010). We show that the evolution of the LF reflects that of the halo formation rate if, for the very massive galaxies represented in our sample, the SFR obeys a simple relationship with halo mass (estimated from the clustering properties) and redshift (Equation (9)), and the lifetime of the main star formation phase is Δt_sf ≈ 7 × 10⁸ yr, consistent with the constraint coming from the α-enhancement observed in the most massive ETGs.

The stellar mass function resulting from the derived SFRs and timescale Δt_sf nicely fits the observed stellar mass function of passively evolving ETGs at z ≈ 1–2.

A comparison of submillimeter and UV LFs indicates that the UV visibility time of massive galaxies is much (≳ 100 times) shorter than the duration of the submillimeter bright phase (Δt_sf ≈ 7 × 10⁸ yr). This implies that dust forms very rapidly after the onset of the main episode of star formation, either in the ISM (Draine 2009; Dunne et al. 2011) or due to the effect of Type II SNe (Dunne et al. 2003, 2009b; Matsuura et al. 2011).

As discussed by Mao et al. (2007), a longer UV bright phase is expected for less massive galaxies. This prediction is supported by the flatter slope of the far-IR LF, compared to the UV one, for SFR ∼ hundreds M_☉ yr⁻¹, implying a decrease with increasing galaxy mass of the ratio of UV to far-IR space densities, i.e., of the ratio between UV and far-IR lifetimes.

In the same vein we find that the duration of the optically bright QSO phase is of order 1/30 of Δt_sf, i.e., of about (2–3) × 10⁷ yr, consistent with independent evidences. On the other hand, the bright end of the QSO LF is somewhat less steep than that of submillimeter galaxies. This is not surprising because, although the evolution of the SFR and of accretion into the central black hole may be well linked, the relationship is mediated by several steps that introduce a substantial scatter that, coupled with the curvature of the LF, flattens its bright tail.

The 100 and 250 μm LFs at different redshifts are quite well reproduced by the physical model of ETG formation and evolution by Granato et al. (2001, 2004), further elaborated by Lapi et al. (2006), without any adjustment of the parameters. As discussed in these papers, the model is built in the framework of the standard hierarchical clustering scenario. Many simulations (e.g., Zhao et al. 2003; Diemand et al. 2007; Genel et al. 2010; Wang et al. 2011) have shown that the growth of a halo occurs in two different phases: a first regime of fast accretion in which the potential well is built up by the sudden mergers of many clumps with comparable masses; and a second regime of slow accretion in which mass is added in the outskirts of the halo, only occasionally affecting the central region where the galactic structure resides. According to the model, the fast accretion phase triggers a burst of star formation that, in massive halos at z ≳ 1, starts an evolutionary sequence that can be summarized as follows. There is an early, short phase of (almost) dust-free star formation when the galaxy shines as a bright UV source. It is followed by a dust-enshrouded star formation phase when the galaxy shines in the far-IR/submillimeter range. The duration of both the UV bright phase and the far-IR/submillimeter bright phase is shorter for the most massive galaxies, with the highest SFRs; for these objects the UV phase lasts ≲ 10⁷ yr and the far-IR/submillimeter phase lasts <10⁹ yr. Then there is a phase, lasting several 10⁷ yr, when the nucleus shines as a bright QSO after having swept away most of the interstellar gas and dust. Finally, passive evolution of the stellar populations follows, and the galaxy evolves into a local ETG.

According to this model, star-forming proto-spheroidal galaxies account for a substantial fraction of the cosmic infrared background (see Figure 19) and dominate the cosmic SFR at z > 1.5, while at lower-z the SFR is dominated by late-type (normal and starburst) galaxies. This model was the basis for the successful predictions of the submillimeter counts of strongly lensed galaxies by Perrotta et al. (2003) and Negrello et al. (2007). It also accurately reproduced the epoch-dependent galaxy LFs in different spectral bands, as well as a variety of relationships among photometric, dynamical, and chemical properties, as shown in previous papers (see Table 2 of Lapi et al. 2006 and additional results, especially on the galaxy chemical evolution, in Mao et al. 2007 and Lapi et al. 2008).

The Herschel-ATLAS is a project with Herschel, which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. The H-ATLAS Web site is http://www.h-atlas.org/

The work has been supported in part by ASI/INAF agreement No. I/009/10/0 and by INAF through the PRIN 2009 "New light on the early Universe with sub-mm spectroscopy." We thank the referee for helpful comments and suggestions, and Cedric Lacey who provided in tabular form the submillimeter counts yielded by the Lacey et al. (2010) model. A. Lapi acknowledges useful discussions with A. Cavaliere, G. L. Granato, P. Salucci, L. Silva, and F. Shankar, and thanks SISSA and INAF-OATS for warm hospitality.