A HYBRID MODEL FOR THE EVOLUTION OF GALAXIES AND ACTIVE GALACTIC NUCLEI IN THE INFRARED

The huge amount of infrared (IR) to millimeter-wave data that has been accumulating over the last several years has not yet led to a fully coherent, established picture of the cosmic star formation history, of the IR evolution of active galactic nuclei (AGNs), and of the interrelations between star formation and nuclear activity.

Many increasingly sophisticated phenomenological models for the cosmological evolution of the galaxy and AGN luminosity functions (LFs) over a broad wavelength range have been worked out (e.g., Béthermin et al. 2012a, 2011; Gruppioni et al. 2011; Rahmati & van der Werf 2011; Marsden et al. 2011; Franceschini et al. 2010; Valiante et al. 2009; Le Borgne et al. 2009; Rowan-Robinson 2009). These models generally include multiple galaxy populations, with different spectral energy distributions (SEDs) and different evolutionary properties, described by simple analytic formulae. In some cases AGNs are also taken into account. All of them, however, admittedly have limitations.

The complex combination of source properties (both in terms of the mixture of SEDs and of evolutionary properties), called for by the richness of data, results in a large number of parameters, implying substantial degeneracies that hamper the interpretation of the results. The lack of constraints coming from the understanding of the astrophysical processes controlling the evolution and the SEDs limits the predictive capabilities of these models. In fact, predictions of pre-Herschel phenomenological models, matching the data then available, yielded predictions for Herschel counts quite discrepant from each other and with the data.

The final goal is a physical model linking the galaxy and AGN formation and evolution to primordial density perturbations. In this paper we make a step in this direction presenting a comprehensive "hybrid" approach, combining a physical, forward model for spheroidal galaxies and the early evolution of the associated AGNs with a phenomenological backward model for late-type galaxies and for the later AGN evolution. We start from the consideration of the observed dichotomy in the ages of stellar populations of early-type galaxies on one side and late-type galaxies on the other. Early-type galaxies 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. Thus proto-spheroidal galaxies are the dominant star-forming population at z ⩾ 1.5, while IR galaxies at z < 1.5 are mostly late-type "cold" (normal) and "warm" (starburst) galaxies.

Fuller hierarchical galaxy formation models, whereby the mass assembly of galaxies is related to structure formation in the dark matter and the star formation and merger histories of galaxies of all morphological types are calculated based on physical prescriptions have been recently presented by several groups (Lacey et al. 2008; Fontanot et al. 2009; Narayanan et al. 2010; Shimizu et al. 2012). However, the predictions for the IR evolution of galaxies are limited to a small set of wavelengths and frequently highlight serious difficulties with accounting for observational data (Lacey et al. 2010; Niemi et al. 2012; Hayward et al. 2013).

While the evolution of dark matter halos in the framework of the "concordance" ΛCDM cosmology is reasonably well understood thanks to N-body simulations such as the Millennium, the Millennium-XXL, and the Bolshoi simulations (Springel et al. 2005; Boylan-Kolchin et al. 2009; Angulo et al. 2012; Klypin et al. 2011), establishing a clear connection between dark matter halos and visible objects proved to be quite challenging, especially at (sub-)millimeter wavelengths. The early predictions of the currently favored scenario, whereby both the star formation and the nuclear activity are driven by mergers, were more than one order of magnitude below the observed SCUBA 850 μm counts (Kaviani et al. 2003; Baugh et al. 2005). The basic problem is that the duration of the star formation activity triggered by mergers is too short, requiring non-standard assumptions either on the initial mass function (IMF) or on dust properties to account for the measured source counts. The problem is more clearly illustrated in terms of redshift-dependent far-IR/submillimeter LF, estimated on the basis of Herschel data (Eales et al. 2010; Gruppioni et al. 2010; Lapi et al. 2011). These estimates consistently show that z ≃ 2 galaxies with star formation rates (SFRs) SFR ≃ 300 M_☉ yr⁻¹ have comoving densities Φ₃₀₀ ∼ 10⁻⁴ Mpc⁻³ dex⁻¹. The comoving density of the corresponding halos is n(M_vir) ∼ Φ₃₀₀(t_exp/τ_SFR), where M_vir is the total virial mass (mostly dark matter), τ_SFR is the lifetime of the star-forming phase, and t_exp is the expansion timescale. For the fiducial lifetime τ_SFR ≃ 0.7 Gyr advocated by Lapi et al. (2011), log (M_vir/M_☉) ≃ 12.92, while for τ_SFR ≃ 0.1 Gyr, typical of a merger-driven starburst, log (M_vir/M_☉) ≃ 12.12. Thus while the Lapi et al. (2011) model implies an SFR/M_vir ratio easily accounted for on the basis of standard IMFs and dust properties, the latter scenario requires an SFR/M_vir ratio more than a factor of six higher.

To reach the required values of SFR/M_vir or, equivalently, of L_IR/M_vir, Baugh et al. (2005) resorted to a top-heavy IMF while Kaviani et al. (2003) assumed that the bulk of the submillimeter emission comes from a huge amount of cool dust. But even tweaking with the IMF and with dust properties, fits of the submillimeter counts obtained within the merger-driven scenario (Lacey et al. 2010; Niemi et al. 2012) are generally unsatisfactory. Further constraints on physical models come from the clustering properties of submillimeter galaxies that are determined by their effective halo masses. As shown by Xia et al. (2012), both the angular correlation function of detected submillimeter galaxies and the power spectrum of fluctuations of the cosmic infrared background (CIB) indicate halo masses larger than implied by the major mergers plus top-heavy initial stellar mass function scenario (Kim et al. 2012) and smaller than implied by cold flow models but consistent with the self-regulated baryon collapse scenario (Granato et al. 2004; Lapi et al. 2006, 2011).

As is well known, the strongly negative K-correction emphasizes high-z sources at (sub-)millimeter wavelengths. The data show that the steeply rising portion of the (sub-)millimeter counts is indeed dominated by ultra-luminous star-forming galaxies with a redshift distribution peaking at z ≃ 2.5 (Chapman et al. 2005; Aretxaga et al. 2007; Yun et al. 2012; Smolčić et al. 2012). As shown by Lapi et al. (2011), the self-regulated baryon collapse scenario provides a good fit of the (sub-)millimeter data (counts, redshift-dependent LFs) as well as of the stellar mass functions at different redshifts. Moreover, the counts of strongly lensed galaxies were predicted with remarkable accuracy (Negrello et al. 2007, 2010; Lapi et al. 2012; González-Nuevo et al. 2012). Further considering that this scenario accounts for the clustering properties of submillimeter galaxies (Xia et al. 2012), we conclude that it is well grounded, and we adopt it for the present analysis. However, we upgrade this model in two respects. First, while on one side, the model envisages a co-evolution of spheroidal galaxies and active nuclei at their centers, the emissions of the two components have so far been treated independently of each other. This is not a problem in the wavelength ranges where one of the two components dominates, as in the (sub-)millimeter region where the emission is dominated by star formation, but is no longer adequate at mid-IR wavelengths, where the AGN contribution may be substantial. In this paper, we present and exploit a consistent treatment of proto-spheroidal galaxies including both components. Second, while the steeply rising portion of (sub-)millimeter counts is fully accounted for by proto-spheroidal galaxies, late-type (normal and starburst) galaxies dominate both at brighter and fainter flux densities and over broad flux density ranges at mid-IR wavelengths. At these wavelengths, AGNs not associated with proto-spheroidal galaxies but either to evolved early-type galaxies or to late-type galaxies are also important. Since we do not have a physical evolutionary model for late-type galaxies and the associated AGNs, these source populations have been dealt with adopting a phenomenological approach.

Another distinctive feature of the present model is that we have attempted to simultaneously fit the data over a broad wavelength range, from mid-IR to millimeter waves. As mentioned in several papers, this presents several challenges. First, the data come from different instruments and the relative calibration is sometimes problematic (see the discussion in Béthermin et al. 2011). Systematic calibration offsets may hinder simultaneous fits of different data sets. For example, Marsden et al. (2011) pointed out that there is considerable tension between the SCUBA 850 μm counts and the AzTEC counts at 1.1 mm, and indeed the 850 μm and millimeter-wave counts have been repeatedly corrected (generally downward) as biases were discovered and better data were acquired. Also, the very complex SEDs in the mid-IR, where strong polycyclic aromatic hydrocarbon (PAH) emission features show up, make the counts exceedingly sensitive to the details of the spectral response function of the specific instrument and introduce large uncertainties in the conversion from broadband measurements to monochromatic flux densities giving rise to strong discrepancies among data sets nominally referring to the same wavelength. In fact, large discrepancies are present among different determinations of 15 μm and 60 μm source counts.

The plan of the work is as follows. In Section 2, we describe the physical model for the evolution of proto-spheroidal galaxies and of the associated AGNs and the SEDs adopted for these sources. Section 3 deals with the evolutionary model for late-type galaxies and z ⩽ 1.5 AGNs. In Section 4, we present the formalism to compute the source counts of unlensed and lensed sources, the cumulative flux density as a function of redshift, and the contributions to the CIB. In Section 5, we report on the determination of the best-fit values of the model parameters. In Section 6, the model results are compared with data on multi-frequency LFs at various redshifts and on source counts, both total and per redshift slices. The multi-frequency power spectra of CIB fluctuations implied by the model are discussed in Section 7. Finally, Section 8 contains a summary of the paper and our main conclusions.

Tabulations of multi-frequency model counts, redshift distributions, SEDs, redshift-dependent LFs at several wavelengths, and a large set of figures comparing model predictions with the data are available at the Web site http://people.sissa.it/~zcai/galaxy_agn/.

Throughout this paper we adopt a flat cosmology with present-day matter and baryon density in units of the critical density, Ω_{m, 0} = 0.27 and Ω_{b, 0} = 0.044; Hubble constant ${h}=\hbox{{H}}_0/100=0.71$; spectrum of primordial density perturbations with slope n = 1 and normalization on a scale of 8 h⁻¹ Mpc σ₈ = 0.81.

2.1. Overview of the Model

2.2. Luminosity Functions

The bolometric LF of proto-spheroids is obtained convolving the halo formation rate dN_ST(M_vir, z)/dt with the galaxy luminosity distribution, P(L, z; M_vir). The halo formation rate is well approximated, for z ≳ 1.5, by the positive term of the cosmic time derivative of the halo mass function N_ST. For the latter, giving the average comoving number density of halos of given mass, M_vir, we adopt the Sheth & Tormen (1999) analytical expression

$$\begin{equation} N_{\rm ST}(M_{\rm vir}, z) dM_{\rm vir} = \frac{\bar{\rho }_{\rm m,0}}{M_{\rm vir}^2} f_{\rm ST}(\nu) \frac{d\ln \nu }{d\ln M_{\rm vir}}dM_{\rm vir}, \end{equation} \tag{ 1 }$$

where $\bar{\rho }_{\rm m,0}=\Omega _{\rm m,0} \rho _{\rm c,0}$ is the present-day mean comoving matter density of the universe and ν ≡ [δ_c(z)/σ(M_vir)]², with δ_c(z) = δ₀(z)D(0)/D(z). The critical value of the initial overdensity that is required for spherical collapse at z, δ₀(z), is (Nakamura & Suto 1997)

$$\begin{eqnarray*} \delta _0(z)&= &\frac{3(12\pi)^{2/3}}{20} [1+0.0123\log \Omega _{\rm m}(z)] \nonumber \\ &\simeq & 1.6865[1+0.0123\log \Omega _{\rm m}(z)]. \end{eqnarray*}$$

The linear growth factor can be approximated as (Lahav et al. 1991; Carroll et al. 1992)

$$\begin{eqnarray*} &&D(z)=\frac{5\Omega _{\rm m}(z)}{2(1+z)}/\left[\frac{1}{70} +\frac{209}{140}\Omega _{\rm m}(z) - \frac{1}{140}\Omega _{\rm m}^2(z) + \Omega _{\rm m}^{4/7}(z)\right]. \end{eqnarray*}$$

The mass variance σ(M_vir) of the primordial perturbation field smoothed on a scale containing a mass M_vir with a top-hat window function was computed using the Bardeen et al. (1986) power spectrum with correction for baryons (Sugiyama 1995), for our choice of cosmological parameters (see Section 1). The results are accurately approximated (error <1% over a broad range of M_vir, 10⁶ < M_vir/M_☉ < 10¹⁶) by

$$\begin{eqnarray} \sigma (M_{\rm vir}) &=& \frac{0.8}{0.84} [ 14.110393-1.1605397x-0.0022104939x^2\nonumber \\ && + 0.0013317476x^3 -2.1049631\times 10^{-6}x^4 ] \end{eqnarray} \tag{ 2 }$$

where x ≡ log (M_vir/M_☉). Furthermore

$$\begin{eqnarray*} &&f_{\rm ST}(\nu) = A [1+(a\nu)^{-p}] \left(\frac{a\nu }{2}\right)^{1/2} \frac{e^{-a\nu /2}}{\pi ^{1/2}}, \end{eqnarray*}$$

where A = 0.322, p = 0.3, and a = 0.707.

The halo formation rate is then

$$\begin{eqnarray} &&{{dN_{\rm ST}(M_{\rm vir}, z)\over dt}} = N_{\rm ST}(M_{\rm vir}, z){{d\ln f_{\rm ST}(\nu)\over dt}}\nonumber \\\ms && \quad = - N_{\rm ST}(M_{\rm vir}, z){\left[{a\delta _{\rm c}\over \sigma ^2} + {2p\over \delta _{\rm c}} {\sigma ^{2p}\over \sigma ^{2p}+a^p \delta _{\rm c}^{2p}} - {1\over \delta _{\rm c}} \right]{d\delta _{\rm c}\over dz} {dz\over dt}} \nonumber \\\ms && \quad \simeq N_{\rm ST}(M_{\rm vir}, z){\left[{a\nu \over 2} + {p\over 1+(a\nu)^{p}} \right]{d\ln \nu \over dz} |{dz\over dt}|}, \end{eqnarray} \tag{ 3 }$$

where dz/dt = −H₀(1 + z)E(z) with $E(z) \equiv \sqrt{\Omega _{\Lambda ,0} + \Omega _{\rm m,0}(1+z)^3}$.

The comoving differential LF Φ(log L, z), i.e., the number density of galaxies per unit log L interval at redshift z, is given by

$$\begin{eqnarray} \Phi (\log L, z) &=&\!\!\! \int ^{M^{\rm max}_{\rm vir}}_{M^{\rm min}_{\rm vir}}\! dM_{\rm vir} \int ^{z^{\rm max}_{\rm vir}}_z\! dz_{\rm vir} \left |\frac{dt_{\rm vir}}{dz_{\rm vir}}\right | \frac{dN_{\rm ST}}{dt_{\rm vir}}(M_{\rm vir}, z_{\rm vir}) \nonumber \\\ms && \cdot P(\log L,z; M_{\rm vir}, z_{\rm vir}), \end{eqnarray} \tag{ 4 }$$

where P(log L, z; M_vir, z_vir) is the luminosity distribution of galaxies at redshift z inside a halo of mass M_vir virialized at redshift z_vir. We set $z^{\rm min}_{\rm vir}=1.5$ and $z^{\rm max}_{\rm vir}=12$.

As mentioned in Section 2.1, the total luminosity of a galaxy is the sum of those of the stellar component and of the active nucleus. For each component, we assume a log-normal luminosity distribution

$$\begin{eqnarray} &&P[\log L|\log \bar{L}] d\log L = \frac{\exp [-\log ^2(L/\bar{L})/2\sigma ^2]}{\sqrt{2\pi \sigma ^2}} d\log L, \nonumber \\ \end{eqnarray} \tag{ 5 }$$

with dispersion σ_* = 0.10 around the mean stellar luminosity $\bar{L}_*(z; M_{\rm vir}, z_{\rm vir})$ and σ_• = 0.35 around the mean AGN luminosity $\bar{L}_\bullet (z; M_{\rm vir}, z_{\rm vir})$. The mean luminosities are computed solving the equations detailed in the Appendix. The higher luminosity dispersion for the AGN component reflects its less direct relationship, compared to the stellar component, with M_vir and z_vir. The distribution of the total luminosity, L_tot = L_* + L_•, is then (Dufresne 2004)

$$\begin{eqnarray} && P[\log L_{\rm tot}|\log \bar{L}_*, \log \bar{L}_\bullet ] d\log L_{\rm tot} = d\log L_{\rm tot} \nonumber \\ && \quad \times {\int ^{\log L_{\rm tot}}_{-\infty } \frac{dx}{2\pi \sigma _* \sigma _\bullet } \frac{L_{\rm tot}}{L_{\rm tot} - 10^x}} \exp \lbrace -(x - \log \bar{L}_*)^2/2\sigma ^2_* \rbrace \nonumber \\ && \quad \times \exp \lbrace -[\log (L_{\rm tot} - 10^x) - \log \bar{L}_\bullet ]^2/2\sigma ^2_\bullet \rbrace. \end{eqnarray} \tag{ 6 }$$

In the upper left panel of Figure 1 we show, as an example, the bolometric LFs at z = 1.5 of the stellar and of the AGN components, as well as the LF of the objects as a whole. As shown in Equation (6), the latter is different from the sum of the first two, although in this case the difference is difficult to perceive. The bright end is dominated by QSOs shining unobstructed after having swept away the interstellar medium (ISM) of the host galaxy. In this phase the QSOs reach their maximum luminosity. Around log (L_bol/L_☉) ≃ 13 the AGNs and the starbursts give similar contributions to the bolometric LF. The inflection at log (L_IR/L_☉) ≃ 11.7 corresponds to the transition from the regime where the feedback is dominated by SNe (lower halo masses) to the regime where it is dominated by AGNs. While the star formation in massive halos is abruptly stopped by the AGN feedback after 0.5–0.7 Gyr, it lasts much longer in smaller galaxies, implying a fast increase of their number density.

Zoom In Zoom Out Reset image size

The upper right panel of the same figure illustrates the evolution with cosmic time of the global LF. The cooling and free-fall timescales shorten with increasing redshift because of the increase of the matter density and this drives a positive luminosity evolution, thwarted by the decrease in the comoving density of massive halos. The two competing factors result, for both the starburst and the AGN component (see the lower panels of the figure), in a positive evolution up to z ≃ 2.5 followed by a decline at higher z, consistent with the observational determinations by Gruppioni et al. (2010) and Lapi et al. (2011) for the starburst component and by Assef et al. (2011) and Brown et al. (2006) for AGNs. The decrease of the LF at low luminosities, more clearly visible at the higher redshifts, is an artifact due to the adopted lower limit for the considered halo masses. This part of the LF, however, does not contribute significantly to the observed statistics and therefore is essentially irrelevant here. Below the minimum virialization redshift, $z^{\rm min}_{\rm vir}=1.5$, the bolometric LF of proto-spheroidal galaxies rapidly declines as they evolve toward the "passive" phase. The decline is faster at the bright end (above log (L_bol/L_☉) ≃ 12) since the switching off of the star formation for the more massive halos occurs on a shorter timescale.

The monochromatic LFs of each component or of objects as a whole can obviously be computed using the same formalism, given the respective SEDs. We define $\bar{\cal L}_{*, \nu } \equiv \nu \bar{L}_{*, \nu } = \nu f_*(\nu) \bar{L}_{*, {\rm IR}}$, $\bar{\cal L}_{\bullet , \nu } \equiv \nu \bar{L}_{\bullet , \nu } = \nu f_\bullet (\nu) \bar{L}_{\bullet , {\rm bol}}$, and $\bar{\cal L}_{\nu } \equiv \bar{\cal L}_{*, \nu } + \bar{\cal L}_{\bullet , \nu }$, where f(ν) is the normalized SED (∫dν f(ν) = 1).

Since the model cannot follow in detail the evolution of the AGN SEDs during the short phase when they shine unobstructed by the ISM of the host galaxy, the distinction between obscured and unobscured AGNs in the model is made in two ways. First, following Lapi et al. (2006), we choose a fixed optical (B-band) "visibility time," Δt_vis = 5 × 10⁷ yr, consistent with current estimates of the optically bright QSO phase. Alternatively, we set the beginning of the optical bright phase at the moment when the gas mass fraction is low enough to yield a low optical depth. We estimate that this corresponds to a gas fraction within the dark matter potential well f_gas = M_gas/M_vir ≲ f_{gas, crit} = 0.03. The two approaches give very similar results and we have chosen the criterion f_gas ≲ f_{gas, crit} to compute the LFs at optical wavelengths.

2.3. Spectral Energy Distributions

Although there is evidence that the galaxy SEDs vary with luminosity (e.g., Smith et al. 2012), Lapi et al. (2011) have shown that the submillimeter data can be accurately reproduced using a single SED for proto-spheroidal galaxies, i.e., the SED of the strongly lensed z ≃ 2.3 galaxy SMM J2135-0102 (Swinbank et al. 2010; Ivison et al. 2010), modeled using GRASIL (Silva et al. 1998). The basic reason for the higher uniformity of the SEDs of high-z active star-forming galaxies compared to galaxies at low-z is that the far-IR emission of the former objects comes almost entirely from dust in molecular clouds, heated by newly formed stars, while in low-z galaxies there are important additional contributions from colder "cirrus" heated by older stellar populations.

This SED worked very well at submillimeter wavelengths but yielded millimeter-wave counts in excess of the observed ones. To overcome this problem the submillimeter slope of the SED has been made somewhat steeper, preserving the consistency with the photometric data on SMM J2135-0102 (see Figure 2). Moreover, the SED used by Lapi et al. has a ratio between the total (8–1000 μm) IR and the 8 μm luminosity (IR8 = L_IR/L₈) of ≃ 30, far higher than the mean value for z ≃ 2 galaxies (IR8 ≃ 9; Reddy et al. 2012). We have therefore modified the near- and mid-IR portions of the SED adopting a shape similar to that of Arp 220. The contribution of the passive evolution phase of early-type galaxies is small in the frequency range of interest here and will be neglected.

Zoom In Zoom Out Reset image size

As mentioned in Section 2.1, the model follows the AGN evolution through two phases (a third phase, reactivation, will be considered in Section 3.2). For the first phase, when the black hole growth is enshrouded by the abundant, dusty ISM of the host galaxy, we adopt the SED of a heavily absorbed AGN taken from the AGN SED library by Granato & Danese (1994). Note that these objects differ from the classical type 2 AGNs because they are not obscured by a circumnuclear torus but by the more widely distributed dust in the host galaxy. They will be referred to as type 3 AGNs. In the second phase the AGN shines after having swept out the galaxy ISM. For this phase, we adopted the mean QSO SED by Richards et al. (2006) extended to submillimeter wavelengths assuming a graybody emission with dust temperature T_dust = 80 K and emissivity index β = 1.8. These SEDs imply that the IR (8–1000 μm) band comprises 92% of the bolometric luminosity of obscured AGNs and 19% of that of the unobscured ones. As illustrated by Figure 2, except in the rare cases in which the AGN bolometric luminosity is much larger than that of the starburst, the AGN contribution is small at (sub-)millimeter wavelengths, while it is important and may be dominant in the mid-IR. This implies that the statistics discussed here are insensitive to the parameters describing the extrapolation of the Richards et al. SED to (sub-)millimeter wavelengths.

Figure 3 shows the global SEDs and the contributions of the stellar and AGN components for two galaxy ages and three host halo masses virialized at z_vir = 3. The shorter evolution timescale of the AGNs is clearly visible. It is worth noting that the effect of feedback as a function of halo mass on the SFR is very different from that on accretion onto the supermassive black hole. In the less massive halos the AGN feedback has only a moderate effect on the evolution of the SFR and of the accretion rate, which are mostly controlled by the SN feedback. With reference to the figure, for log (M_vir/M_☉) = 11.4, the star formation continues at an almost constant rate for a few Gyr. On the other hand, the accretion rate onto the central black hole is at the Eddington limit only up to an age of ≃ 0.3 Gyr and afterward drops to a strongly sub-Eddington regime. This is because the growth rate of the reservoir is approximately proportional to the SFR (and therefore slowly varying for few Gyr), while the accretion rate grows exponentially until the mass contained in the reservoir is exhausted. From this moment on the accretion rate is essentially equal to the (strongly sub-Eddington) inflow rate. For more massive halos the quenching of both the SFR and of the accretion occurs more or less simultaneously at ages of ≃ 0.5–0.6 Gyr, but while the SFR stops very rapidly, the AGN activity continues until the flow of the matter accumulated in the reservoir runs out. At ages ≳ 0.6 Gyr, the more massive galaxies are in passive evolution and therefore very weak in the far-IR, while star formation and the dust emission are still present in lower-mass galaxies.

Zoom In Zoom Out Reset image size

3.1. Late-type and Starburst Galaxies

We consider two z ≲ 1.5 galaxy populations: "warm" starburst galaxies and "cold" (normal) late-type galaxies. For the IR LF of both populations, we adopt the functional form advocated by Saunders et al. (1990):

$$\begin{eqnarray} && \Phi (\log L_{\rm IR}, z) d\log L_{\rm IR} = \Phi ^* \left (\frac{L_{\rm IR}}{L^*}\right)^{1-\alpha } \nonumber \\ && \quad \times \exp \left [- \frac{\log ^2 (1+L_{\rm IR}/L^*)}{2\sigma ^2} \right ] d\log L_{\rm IR}, \end{eqnarray} \tag{ 7 }$$

where the characteristic density Φ* and luminosity L*, the low-luminosity slope α and the dispersion σ of each population are, in principle, free parameters. However, the low-luminosity portion of the LF is dominated by "cold" late-type galaxies and, as a consequence, the value of α of the warm population is largely unconstrained; we have fixed it at α_warm = 0.01. In turn, the "warm" population dominates at high luminosities so that the data only imply an upper limit to σ_cold. We have set σ_cold = 0.3.

For the "warm" population, we have assumed power-law density and luminosity evolution [$\Phi ^*(z) = \Phi ^*_0(1+z)^{\alpha _\Phi }$; $L^*(z)=L^*_0(1+z)^{\alpha _{\rm L}}$] up to z_break = 1, $\alpha _\Phi$ and α_L being free parameters. The "cold" population comprises normal disk galaxies for which chemo/spectrophotometric evolution models (Mazzei et al. 1992; Colavitti et al. 2008) indicate a mild (a factor ≃ 2 from z = 0 to z = 1) increase in the SFR, hence of IR luminosity, with look-back time. Based on these results we take, for this population, α_L = 1 and no density evolution. At z > z_break both Φ*(z) and L*(z) are kept to the values at z_break multiplied by the smooth cutoff function {1 − erf[(z − z_cutoff)/Δz)]}/2, with z_cutoff = 2 and Δz = 0.5. The choice of the redshift cutoff for both populations of late-type galaxies is motivated by the fact that the disk component of spirals and the irregular galaxies are characterized by relatively young stellar populations (formation redshift z ≲ 1–1.5). Above z = 1.5, proto-spheroidal galaxies (including bulges of disk galaxies) dominate the contribution to the LF, at least in the observationally constrained luminosity range. The other parameters are determined by minimum χ² fits to selected data sets, as described in Section 5. Their best-fit values and the associated uncertainties are listed in Table 1.

Table 1. Parameters for Low-z AGNs and for "Warm" and "Cold" Galaxy Populations

	AGN 1	AGN 2	Warm	Cold
	(12 μm)	(12 μm)	(IR)	(IR)
$\log (\Phi ^*_0/{\rm Mpc}^{-3})$	−5.409 ± 0.098	−4.770 ± 0.122	−2.538 ± 0.051	−1.929 ± 0.112
$\log (L^*_0/L_\odot)$	9.561 ± 0.084	10.013 ± 0.093	10.002 ± 0.076	9.825 ± 0.087
α	1.1	1.5	0.01	1.372 ± 0.121
σ	0.627 ± 0.017	0.568 ± 0.021	0.328 ± 0.014	0.3
$\alpha _\Phi$	2.014 ± 0.400	4.499 ± 0.317	0.060 ± 0.200	0.0
αL	2.829 ± 0.297	0.0	3.625 ± 0.097	1.0
zbreak	1.0	1.0	1.0	1.0
zcutoff	2.0	2.0	2.0	2.0

Notes. The parameters of the AGN luminosity functions refer to 12 μm (νL_ν), while those for galaxies refer to IR (8–1000 μm) luminosities. Values without error were kept fixed.

Download table as:  ASCIITypeset image

Although there is clear evidence of systematic variations of the IR SEDs of low-z galaxies with luminosity (e.g., Smith et al. 2012), we tried to fit the data with just two SEDs, one for the "warm" and one for the "cold" population. These SEDs were generated by combining those of Dale & Helou (2002), which are best determined at mid-IR wavelengths, with those of Smith et al. (2012), primarily based on Herschel data in the range 100–500 μm. Dale & Helou (2002) give SED templates for several values of the 60–100 μm flux density ratio, log [f_ν(60 μm)/f_ν(100 μm)]. Using the relation between this ratio and the 3–1100 μm luminosity, L_TIR, given by Chapman et al. (2003), we established a correspondence between their SEDs and those by Smith et al., labeled by the values of log (L_IR/L_☉). The combined SEDs are based on Smith et al. above 100 μm and on Dale & Helou at shorter wavelengths. By trial and error, we found that the best fit to the data is obtained using for the "cold" population the SED corresponding to log (L_IR/L_☉) = 9.75 (actually the SEDs change very slightly for log (L_IR/L_☉) ≲ 9.75) and for the "warm" population the SED corresponding to log (L_IR/L_☉) = 11.25. These two SEDs are displayed in Figure 4.

Zoom In Zoom Out Reset image size

3.2. Reactivated AGNs

In the framework of our reference galaxy and AGN evolutionary scenario, most of the growth of supermassive black holes is associated with the star-forming phase of spheroidal components of galaxies at z ≳ 1.5 when the great abundance of ISM favors high accretion rates, at, or even slightly above, the Eddington limit. At later cosmic times the nuclei can be reactivated by, e.g., interactions, mergers, or dynamical instabilities. The accretion rates are generally strongly sub-Eddington. Our evolutionary scenario cannot predict their amplitudes and duty cycles. We therefore also adopted for these objects a phenomenological backward evolution model analogous to that used for the "warm" galaxy population, i.e., LFs of the same form of Equation (7) and power-law density and luminosity evolution with the same break and cutoff redshifts. However, the parameters of the LFs refer to 12 μm (see Section 3.1 and Table 1). The data do not allow a determination of the slopes, α, of the faint portions of the LFs. We have set α = 1.1 for type 1 AGNs and α = 1.5 for type 2. The steeper slope for type 2 was chosen on account of the fact that these dominate over type 1 at low luminosities. As in the case of normal late-type and of starburst galaxies, the other parameters are obtained by minimum χ² fits, as detailed in Section 5, and the best-fit values are listed, with their uncertainties, in Table 1. For type 2 AGNs, pure density evolution was found to be sufficient to account for the data.

For type 1 AGNs we adopted the mean QSO SED by Richards et al. (2006), extended to millimeter wavelengths as described in Section 2.3, while for type 2 AGNs we adopted the SED of the local AGN-dominated ULIRG Mrk 231, taken from the SWIRE library (Polletta et al. 2007). These SEDs are shown in Figure 5, where the SED of type 3 AGNs associated with dusty star-forming proto-spheroidal galaxies is also plotted for comparison. The SED of type 3 AGNs is the most obscured at optical/near-IR wavelengths due to the effect of the dense, dusty ISM of the high-z host galaxies. This means that the counts at optical/near-IR wavelengths are dominated by type 1 AGNs with type 2 AGNs becoming increasingly important in the mid-IR. The three AGN populations have approximately the same ratio between the rest-frame 12 μm and the bolometric luminosity, as first pointed out by Spinoglio & Malkan (1989).

Zoom In Zoom Out Reset image size

The type 1/type 2 space density ratio yielded by the model increases with luminosity, consistent with observations (e.g., Burlon et al. 2011) and with the receding torus model (Lawrence 1991). In the framework of the standard unified model of AGNs type 1 and type 2 AGNs differ only in terms of the angle which the observers' line of sight makes with the axis of a dusty torus. If the line of sight to the central region is blocked by the torus, the AGN is seen as a type 2. According to the receding torus model, the opening angle of the torus (measured from the torus axis to the equatorial plane) is larger in more luminous objects, implying that obscuration is less common in more luminous AGNs. Since our model implies that type 1 AGNs (but not type 2's) are evolving in luminosity, they become increasingly dominant with increasing redshift.

The surface density of sources per unit flux density and redshift interval is

$$\begin{equation} \frac{d^3N(S_\nu , z)}{dS_\nu dz d\Omega }= \frac{\Phi (\log {L}_{\nu ^{\prime }}, z)}{L_{\nu ^{\prime }}\ln 10} \frac{dL_{\nu ^{\prime }}}{dS_\nu } \frac{d^2V}{dzd\Omega }, \end{equation} \tag{ 8 }$$

where ν' = ν(1 + z),

$$\begin{equation} S_\nu = \frac{(1+z)L_{\nu ^{\prime }}}{4\pi D^2_{\rm L}(z)}, \end{equation} \tag{ 9 }$$

the comoving volume per unit solid angle is

$$\begin{equation} \frac{d^2V}{dzd\Omega } = \frac{c}{H_0} \frac{(1+z)^2D^2_{\rm A}(z)}{E(z)}, \end{equation} \tag{ 10 }$$

and the luminosity distance D_L and the angular diameter distance D_A are related, in a flat universe, by

$$\begin{equation} \frac{D_{\rm L}}{1+z} = (1+z) D_{\rm A} = \frac{c}{H_0} \int ^z_0 \frac{dz^{\prime }}{E(z^{\prime })}. \end{equation} \tag{ 11 }$$

The differential number counts, i.e., the number of galaxies with flux density in the interval S_ν ± dS_ν/2 at an observed frequency ν per unit solid angle, are then

$$\begin{equation} \frac{d^2N}{dS_\nu d\Omega }(S_\nu) = \int ^{z_{\rm max}}_{z_{\rm min}} dz \frac{\Phi (\log {L}_{\nu ^{\prime }}, z)}{L_{\nu ^{\prime }}\ln 10} \frac{dL_{\nu ^{\prime }}}{dS_\nu } \frac{d^2V}{dzd\Omega }. \end{equation} \tag{ 12 }$$

The integral number counts, i.e., the number of galaxies with flux density S_ν > S_{ν, inf} at frequency ν per unit solid angle, are given by

$$\begin{eqnarray} \frac{dN}{d\Omega }(S_\nu > S_{\nu ,\rm inf}) &=& \int ^{z_{\rm max}}_{z_{\rm min}} dz \frac{d^2V}{dzd\Omega } \int ^{\infty }_{\log { L}_{\nu ^{\prime },\rm inf}} \nonumber \\ && \times \Phi (\log {L}_{\nu ^{\prime }}, z) d\log { L}_{\nu ^{\prime }}, \end{eqnarray} \tag{ 13 }$$

where ν' = (1 + z)ν and $L_{\nu ^{\prime },\rm inf}$ is the monochromatic luminosity of a source at the redshift z observed to have a flux density S_{ν, inf}. Counts (per steradian) dominated by local objects (z ≪ 1) can be approximated as

$$\begin{eqnarray} &&S^{2.5}_{\nu } \frac{d^2N}{dS_\nu d\Omega } \simeq \frac{1}{4\pi } \frac{1}{4\sqrt{\pi }} \int ^{\infty }_{0} \Phi (\log {L}_\nu , z \simeq 0) L_\nu ^{3/2} d\log {L}_\nu. \nonumber\\ \end{eqnarray} \tag{ 14 }$$

The redshift distribution, i.e., the surface density of sources with observed flux densities greater than a chosen limit S_{ν, inf} per unit redshift interval, is

$$\begin{equation} \frac{d^2N}{dzd\Omega }(z, S_\nu > S_{\nu , \rm inf}) = \int ^{\infty }_{S_{\nu , \rm inf}} \frac{d^3N}{dS^\prime _{\nu }dzd\Omega } dS^\prime _{\nu }. \end{equation} \tag{ 15 }$$

The steepness of the (sub-)millimeter counts of proto-spheroidal galaxies and their substantial redshifts imply that their counts are strongly affected by the magnification bias due to gravitational lensing (Blain 1996; Perrotta et al. 2002, 2003; Negrello et al. 2007):

$$\begin{equation} \frac{d^3N_{\rm lensed}(S_\nu , z)}{d\log S_{\nu } dz d\Omega } = \int _{\mu } d\mu \frac{d^3N(S_\nu /\mu , z)}{d\log S_\nu dz d\Omega } \frac{dP(\mu |z)}{d\mu }, \end{equation} \tag{ 16 }$$

where dP/dμ is the amplification distribution that describes the probability for a source at redshift z to be amplified by factor μ. Here, we have approximated to unity the factor 1/〈μ〉 that would have appeared on the right-hand side, as appropriate for large-area surveys (see Jain & Lima 2011).

We have computed dP/dμ using the SISSA model worked out by Lapi et al. (2012). The differential counts including the effect of lensing can be computed integrating Equation (16) over z. The effect of lensing on counts of other source populations and on proto-spheroidal counts at shorter wavelengths is small and will be neglected in the following.

Interesting constraints on the halo masses of proto-spheroidal galaxies come from the auto- and cross-correlation functions of intensity fluctuations. A key quantity in this respect is the flux function, d²S_ν/dzdΩ, i.e., the redshift distribution of the cumulative flux density of sources below the detection limit S_{ν, lim}

$$\begin{equation} \frac{d^2S_\nu }{dzd\Omega } = \int ^{S_{\nu , \rm lim}}_{0} \frac{d^3N}{dS_\nu dzd\Omega } S_\nu dS_\nu. \end{equation} \tag{ 17 }$$

The contribution of a source population to the extragalactic background at the frequency ν is

$$\begin{equation} I_\nu = \int ^\infty _0 S_\nu \frac{d^2N(S_\nu)}{dS_\nu d\Omega } dS_\nu. \end{equation} \tag{ 18 }$$

A minimum χ² approach for estimating the optimum values of the parameters of the physical model for proto-spheroidal galaxies and associated AGNs is unfeasible because of the lengthy calculations required. Some small adjustments compared to earlier versions (Granato et al. 2004; Lapi et al. 2006; Mao et al. 2007) were made, by trial and error, to improve the agreement with observational estimates of LFs at z > 1.5. An outline of the model, including the definition of the relevant parameters, is presented in the Appendix. The chosen values are listed in Table 2. Discussions of physically plausible ranges can be found in Granato et al. (2004), Cirasuolo et al. (2005), Lapi et al. (2006), Shankar et al. (2006), Cook et al. (2009), and Fan et al. (2010).

Table 2. Parameters of the Physical Model for the Evolution of Proto-spheroidal Galaxies and Associated AGNs

Parameter	Value	Plausible Range	Description
$\tau ^0_{\rm RD}$	3.0	1–10a	Normalization of optical depth of gas cloud (Equation (A9))
	0.10	0.06–0.42	Black hole accretion radiative efficiency (Equation (A16))
λEdd	1–4	≲a fewb	Redshift-dependent Eddington ratio (Equation (A14))
QSO	3.0	1–10a	Strength of QSO feedback (Equation (A19))
k⋆, IR	3.1	2–4c	Conversion factor from SFR to IR luminosity (Equation (A7))
σ*	0.10	≲ 0.5	Dispersion of mean stellar luminosity (Equation (6))
σ•	0.35	≲ 0.5b	Dispersion of mean AGN luminosity (Equation (6))
fgas, crit	0.03	≲ 0.165	Gas mass fraction at transition
			from obscured to unobscured AGNs (Section  2.2))
SN	0.05	0.01–0.1d	Strength of SN feedback (Equation (A6))
αRD	2.5	1–10e	Strength of radiation drag (Equation (A8))

Notes. The values of the first eight parameters used here are somewhat different from those used in previous papers, but are still well within the plausible ranges listed in Column 3 and discussed in the references given in the footnotes. ^aGranato et al. (2004). ^bLapi et al. (2006). ^cLapi et al. (2011). ^dShankar et al. (2006). ^eA. Lapi et al. (in preparation).

Download table as:  ASCIITypeset image

On the contrary, the minimum χ² approach was applied to late-type/starburst galaxies and to reactivated AGNs. The χ² minimization was performed using the routine MPFIT⁹ exploiting the Levenberg–Marquardt least-squares method (Moré 1978; Markwardt 2009).

The huge amount of observational data in the frequency range of interest here and the large number of parameters coming into play forced us to deal with subsets of parameters at a time using specific data for each subset. The parameters of the evolving AGN LFs were obtained using:

• 1.  
  the B-band local QSO LF of Hartwick & Schade (1990);
• 2.  
  the g-band QSO LFs at z = 0.55 and 0.85 of Croom et al. (2009);
• 3.  
  the z ≲ 0.75, 1.24 μm AGN LFs of Assef et al. (2011);
• 4.  
  the bright end [log (L₆₀/L_☉) ⩾ 12] of the local 60 μm LF of Takeuchi et al. (2003);
• 5.  
  the Spitzer AGN counts at 8 and 24 μm of Treister et al. (2006).

The B- and g-band LFs were used to constrain the parameters of type 1 AGNs (type 2 being important only at the low-luminosity end), while the 1.24 μm LFs were regarded as made by a combination of type 1 and type 2 AGNs, the latter being dominant at low luminosities.

As for the evolving LFs of "warm" and "cold" galaxy populations we used the following data sets:

• 1.  
  the IRAS 60 μm local LF of Soifer & Neugebauer (1991);
• 2.  
  the Planck 350, 550, and 850 μm local LFs of Negrello et al. (2013);
• 3.  
  the Spitzer MIPS counts at 24, 70, and 160 μm of Béthermin et al. (2010);
• 4.  
  the Herschel PACS counts at 160 μm of Berta et al. (2011);
• 5.  
  the Herschel SPIRE counts at 250, 350, and 500 μm (Béthermin et al. 2012b).

The fits of the counts were made after having subtracted the contributions of proto-spheroidal galaxies, which are only important at wavelengths ⩾160 μm. The best-fit values of the parameters are listed in Table 1, where values without errors denote parameters that were kept fixed, as mentioned in Section 3.

In comparing model results with observational data, the instrumental spectral responses were taken into account. This is especially important in the mid-IR because of the complexity of the SEDs due to PAH emission lines. The monochromatic luminosity at the effective frequency ν_eff in the observer's frame is given by

$$\begin{equation} L(\nu _{\rm eff}) \equiv \int T(\nu ^{\prime }) L_{\nu ^{\prime }(1+z)} d\nu ^{\prime } / \int T(\nu ^{\prime }) d\nu ^{\prime }, \end{equation} \tag{ 19 }$$

where T(ν) is spectral response function and the integration is carried out over the instrumental bandpass. When the model is compared with LF data at frequency ν_i (in the source frame) coming from different instruments for sources at redshift z, we use the response function of the instrument for which ν_eff is closest to ν_i/(1 + z). In the case of source counts, we use the response function appropriate for the most accurate data.

6.1. Model versus Observed Luminosity Functions and Redshift Distributions

The most direct predictions of the physical model for proto-spheroidal galaxies are the redshift-dependent SFRs and accretion rates onto the supermassive black holes as a function of halo mass. During the dust enshrouded evolutionary phase, the SFRs can be immediately translated into the IR (8–1000 μm) LFs of galaxies. As mentioned above, according to our model, the transition from the dust obscured to the passive evolution phase is almost instantaneous and we neglect the contribution of passive galaxies to the IR LFs. In turn, the accretion rates translate into bolometric luminosities of AGNs given the mass-to-light conversion efficiency for which we adopt the standard value = 0.1. The SEDs then allow us to compute the galaxy and AGN LFs at any wavelength.

In contrast, the phenomenological model for late-type/starburst galaxies yields directly the redshift-dependent IR LFs and that for reactivated AGNs yields the 12 μm LFs. Again these can be translated to any wavelength using the SEDs described in the previous sections.

In Figure 6, the model IR LFs are compared with observation-based determinations at different redshifts. At z > 1.5, the dominant contributions come from the stellar and AGN components of proto-spheroidal galaxies. These contributions fade at lower redshifts and essentially disappear at z < 1. The model implies that AGNs associated with proto-spheroidal galaxies are important only at luminosities higher than those covered by the Lapi et al. (2011) LFs which therefore have been converted to bolometric LFs using their galaxy SED, i.e., neglecting the AGN contribution, so that log (L_IR/L_☉) = log (L₁₀₀/L_☉) + 0.21 and log (L_IR/L_☉) = log (L₂₅₀/L_☉) + 1.24. At z ⩽ 1.5 "warm" and "cold" star-forming galaxies take over, "cold" galaxies being important only at low luminosities. Type 2 AGNs (long-dashed pink lines) may dominate at the highest IR luminosities, while type 1 AGNs (long-dashed light-blue lines) are always sub-dominant (in the IR).

Zoom In Zoom Out Reset image size

The scale on the top x-axis in Figure 6 gives the SFRs corresponding to the IR luminosities

$$\begin{equation} \log \left (\frac{L_{\rm IR}}{L_\odot }\right) = \log \left (\frac{{\rm SFR}}{M_\odot \ {\rm yr}^{-1}}\right) + 9.892, \end{equation} \tag{ 20 }$$

and is therefore meaningful only to the extent that the AGN contribution is negligible. Moreover, the normalization constant applies to high-z proto-spheroidal galaxies whose IR luminosity comes almost entirely from star-forming regions. For more evolved galaxies, older stellar populations can contribute significantly to the dust heating (da Cunha et al. 2012); therefore L_IR is no longer a direct measure of the SFR and therefore the upper scale has to be taken as purely indicative.

Observational determinations of LFs are available in many wave bands and for many cosmic epochs. The comparison between the model and the observed g-band (0.467 μm) AGN LFs at several redshifts is presented in Figure 7, while the comparison in the J band (1.24 μm) is shown in Figure 8. The conversion from monochromatic absolute AB magnitude M_{λ, AB} to the corresponding monochromatic luminosity νL_ν(λ) is given by $\log (\nu L_\nu /[L_\odot ]) = -0.4 M_{\lambda , \rm AB} - \log (\lambda /[{\rm \mathring{A}}]) + 5.530$. The contribution of type 2 AGNs at z < 1.5 strongly increases from the g to the J band. Apart from the low-luminosity portion of the J-band LF, very likely affected by incompleteness, the agreement between the model and the data is remarkably good.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The comparisons between the global (stellar plus AGN components) luminosity functions yielded by the model and those observationally determined at several redshifts and wavelengths are shown in Figures 9–11. In Figure 12, we compare model and observed redshift distributions at various wavelengths and flux density limits. The comparisons for all the other wavelengths for which estimates of the LF are available can be found in the Web site http://people.sissa.it/~zcai/galaxy_agn/.

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 a substantial fraction of sources have only photometric redshifts. For example, the fraction of photometric redshifts is 91% for the VVDS-SWIRE survey with S_24 μm > 0.4 mJy (Rodighiero et al. 2010), 67.5% for the GOODS-N, and 36% for the GOODS-S samples with S_24 μm > 0.08 mJy (Rodighiero et al. 2010). Only few sources at z > 2 have spectroscopic redshifts (Berta et al. 2011). Note that photometric redshift errors tend to moderate the decline of the distributions at high-z. The effect is analogous to the Eddington bias on source counts: errors move more objects from the more populated lower z bins to the less populated higher z bins than in the other way. Thus, the observed distributions may be overestimated at the highest redshifts. In addition, optical identifications are not always complete. On the whole, observational estimates of LFs and of redshift distributions may be affected by systematic effects difficult to quantify and the true uncertainties may be larger than the nominal values.

6.2. Model versus Observed Source Counts and Contributions to the CIB

Model and observed source counts at wavelengths from 15 μm to 1.38 mm are compared in Figure 13. At wavelengths ⩾350 μm, where, in the present framework, proto-spheroidal galaxies are most important, the model provides a simple physical explanation of the steeply rising portion of the counts, which proved to be very hard to account for by other both physical (Hayward et al. 2013; Niemi et al. 2012; Lacey et al. 2010) and phenomenological (e.g., Béthermin et al. 2012a; Gruppioni et al. 2011) models.

Zoom In Zoom Out Reset image size

Table 3. References for Data on Number Counts (See Figure 13)

Wavelength	Instrument	Field	Reference
(μm)
15, 24	Akari/IRC	NEP-deep	Takagi et al. (2012)
15	Akari/IRC	De-lensed A2218	Hopwood et al. (2010)
15	Akari/IRC	NEP-deep+wide	Pearson et al. (2010)
15	Akari/IRC	CDFS	Burgarella et al. (2009)
15	Akari/IRC	NEP-deep	Wada et al. (2007)
15	ISO/ISOCAM	ELAIS-S	Gruppioni et al. (2002)
15	ISO/ISOCAM	ISOCAM deep surveys	Elbaz et al. (1999)
16	Spitzer/IRS	GOODS-N+S	Teplitz et al. (2011)
24, 70	Spitzer/MIPS	ADF-S	Clements et al. (2011)
24, 70	Spitzer/MIPS	Spitzer legacy fields	Béthermin et al. (2010)
24	Spitzer/MIPS	SWIRE fields	Shupe et al. (2008)
24	Spitzer/MIPS	NDWFS Bootes	Brown et al. (2006)
24	Spitzer/MIPS	GOODS-N	Treister et al. (2006)
24	Spitzer/MIPS	Deep Spitzer fields	Papovich et al. (2004)
70, 100	Herschel/PACS	GOODS, LH, COSMOS	Berta et al. (2011)
70	Spitzer/MIPS	xFLS	Frayer et al. (2006)
70	Spitzer/MIPS	Bootes, Marano, CDF-S	Dole et al. (2004)
100	Herschel/PACS	A2218	Altieri et al. (2010)
250, 500	Herschel/SPIRE	HerMES	Béthermin et al. (2012b)
250, 500	Herschel/SPIRE	H-ATLAS	Clements et al. (2010)
250, 500	Herschel/SPIRE	HerMES	Oliver et al. (2010)
250, 500	Herschel/SPIRE	HerMES P(D)	Glenn et al. (2010)
250, 500	BLAST	BGS P(D)	Patanchon et al. (2009)
500	Herschel/SPIRE	H-ATLAS	Lapi et al. (2012)
550, 850	Planck	Planck all-sky survey	Planck Collaboration (2013)
850	SCUBA	Clusters	Noble et al. (2012)
850	SCUBA	A370	Chen et al. (2011)
850	SCUBA	Clusters	Zemcov et al. (2010)
850	SCUBA	Clusters and NTT-DF	Knudsen et al. (2008)
850	SCUBA	SHADES	Coppin et al. (2006)
850	SCUBA	Clusters	Smail et al. (2002)
870	APEX/LABOCA	Clusters	Johansson et al. (2011)
1100	ASTE/AzTEC	AzTEC blank-field survey	Scott et al. (2012)
1100	ASTE/AzTEC	COSMOS	Aretxaga et al. (2011)
1100	ASTE/AzTEC	ADF-S, SXDF, and SSA22	Hatsukade et al. (2011)
1100	ASTE/AzTEC	GOODS-S	Scott et al. (2010)
1100	JCMT/AzTEC	SHADES	Austermann et al. (2010)
1100	JCMT/AzTEC	COSMOS	Austermann et al. (2009)
1400	SPT	SPT survey	Vieira et al. (2010)

Download table as:  ASCIITypeset image

In our model the sudden steepening of the (sub-)millimeter counts is due to the appearance of proto-spheroidal galaxies that show up primarily at z ≳ 1.5, being mostly in passive evolution at lower redshifts. Their counts are extremely steep because, due to the strongly negative K-correction, the submillimeter flux densities corresponding to a given luminosity are only weakly dependent on the source redshift. Then, since the far-IR luminosity is roughly proportional to the halo mass, the counts reflect the high-z LF whose bright end reflects, to some extent, the exponential decline of halo mass function at high masses. This situation results in a very strong magnification bias due to gravitational lensing (Blain 1996; Perrotta et al. 2002, 2003; Negrello et al. 2007). The counts of strongly lensed galaxies depend on the redshift distribution of the unlensed ones. Thus, the good agreement between the model and the observed counts of strongly lensed galaxies (see the 350 μm, 500 μm, and 1380 μm panels of Figure 13) indicates that the model passes this test on the redshift distribution.

Low-z "warm" and "cold" star-forming galaxy populations become increasingly important with decreasing wavelength. At λ ⩾ 160 μm, proto-spheroidal galaxies yield only a minor contribution to the counts. The AGN (mostly type 2) contribution implied by the model is always sub-dominant. We find a maximum contribution in the mid-IR. At 15 μm it is ≃ 8% up to 1 mJy and then rapidly increases up to ≃ 20% above 10 mJy, while at 24 μm it is ≃ 7%–8% up to 0.5 mJy and increases up to ≃ 20% above 2 mJy, in fair agreement with the observational estimates (Treister et al. 2006; Teplitz et al. 2011).

Another test on the redshift distribution is provided by the estimated counts in different redshift slices (Figure 14), although we caution that the true uncertainties may be larger than the nominal ones since the observational estimates are partly based on photometric redshifts and on stacking. The consistency between the model and the data is reasonably good.

Zoom In Zoom Out Reset image size

Figure 15 shows the contributions of the different populations to the CIB. The model accounts for the full CIB intensity over the whole wavelength range. Only at λ  ⩽  10 μm other galaxy populations, such as passively evolving galaxies, become important. According to the model, for λ  ⩾  350 μm the main contribution to the CIB comes from proto-spheroidal galaxies and the fraction contributed by these objects increases with increasing wavelengths. Below λ  =  350 μm lower z "warm" galaxies take over, with "cold" galaxies adding a minor contribution. AGNs are always sub-dominant. The model gives a total (type 1+type 2+type 3) AGN contribution of 8.6% at 16 μm and of 8.1% at 24 μm. For comparison, Teplitz et al. (2011) estimate a contribution of ∼15% at 16 μm; Treister et al. (2006) and Ballantyne & Papovich (2007) find a contribution of ∼10% at 24 μm. It must be noted that these observational estimates are endowed with substantial uncertainties: on one side they may be too low because strongly obscured AGNs may be missed, and on the other side they may be too high because a significant fraction of the observed emission may come from the host galaxy.

Zoom In Zoom Out Reset image size

An important test of our physical model for the evolution of dusty proto-spheroidal galaxies is provided by their clustering properties that are informative on their halo masses. A specific prediction of our model is that proto-spheroidal galaxies are the main contributors to the CIB at (sub-)millimeter wavelengths with "warm" starburst galaxies becoming increasingly important with decreasing wavelength. Since, in our model, proto-spheroidal galaxies are much more strongly clustered than starburst galaxies, the variation in the mixture with wavelength translates in quantitative predictions on the frequency dependence of the amplitude of the CIB power spectra and on the level of correlations among the maps at different frequencies.

We have updated the analysis by Xia et al. (2012) taking into account the new auto- and cross-frequency power spectra obtained by Viero et al. (2012) from Herschel/SPIRE measurements and the power spectrum at 100 μm derived by Pénin et al. (2012). The latter authors actually also give an estimate of the power spectrum at 160 μm. However, the amplitude of the latter is anomalously large. As an example, for the wavenumber k_θ = 0.03 arcmin⁻¹ we find that the amplitude normalized to the CIB intensity,

$$\begin{equation} {\delta I\over I} = {\left[2\pi k_\theta ^2 P(k_\theta)\right]^{1/2}\over I_{\rm CIB}} \end{equation} \tag{ 21 }$$

(Equation (13) of Viero et al. 2012), is ≃ 0.08–0.09 at 100, 250, 350, and 500 μm but jumps to ≃ 0.2 at 160 μm. Since such a jump over a small wavelength range looks odd, we decided not to use the 160 μm power spectrum.

All the relevant details on the formalism used are given by Xia et al. (2012). Briefly, the power spectrum of the galaxy distribution is parameterized as the sum of the one-halo term that dominates on small scales and depends on the distribution of galaxies within the same halo, and the two-halo term that dominates on large scales and is related to correlations among different halos. The halo occupation distribution, which is a statistical description of how dark matter halos are populated with galaxies, is modeled using a central-satellite formalism (see, e.g., Zheng et al. 2005). This assumes that the first galaxy to be hosted by a halo lies at its center, while any remaining galaxies are classified as satellites and are distributed in proportion to the halo mass profile. The mean halo occupation function of satellite galaxies is parameterized as: $\langle {N_{\rm sat}}\rangle \propto (M_{\rm vir}/M_{\rm sat})^{\alpha _{\rm sat}}$, where M_vir is the halo mass and the power-law index α_sat is a free parameter. The key parameter in the two-halo term is the minimum halo mass, M_{vir, min}, that determines the amplitude of the effective bias function b_eff(z).

In the Xia et al. (2012) paper the only free parameters are the minimum halo mass, M_{min, protosph}, and the power-law index of the mean occupation function of satellites, α_{sat, protosph}, of proto-spheroidal galaxies. This is because the contribution of late-type galaxies to the power spectra at λ ⩾ 250 μm is always sub-dominant and therefore the parameters characterizing their clustering properties were poorly constrained. This is no longer true if we add the 100 μm power spectrum, which, however, still provides only weak constraints on α_{sat, late-type}. We therefore fixed that parameter to α_{sat, late-type} = 1. The fits to the Herschel/SPIRE power spectra determined by Viero et al. (2012) give log (M_{min, protosph}/M_☉) = 12.15 ± 0.04 and α_{sat, protosph} = 1.55 ± 0.05 (1σ errors), close to the values found by Xia et al. (2012). The 100 μm data do not constrain these parameters further but yield log (M_{min, late-type}/M_☉) = 11.0 ± 0.06. The nominal errors on each parameter have been computed marginalizing on the other and correspond to Δχ² = 1. We caution that the true uncertainties are likely substantially higher than the nominal values, both because the model relies on simplifying assumptions that may make it too rigid and because of possible systematics affecting the data. Our value of M_{min, protosph} implies an effective halo mass (Equation (17) of Xia et al. 2012) at z ≃ 2 of proto-spheroidal galaxies, making up most of the CIB, M_eff ≃ 4.5 × 10¹² M_☉. This value is close to the estimated halo mass of the most effective star formers in the universe. Tacconi et al. (2008) estimated their mean comoving density at z ∼ 2 to be ∼2 × 10⁻⁴ Mpc⁻³. For the standard ΛCDM cosmology, this implies that they are hosted by dark matter halos of ∼3.5 × 10¹² M_☉.

The best-fit model power spectra are plotted in Figure 16 where the one- and two-halo contributions of proto-spheroidal and late-type galaxies are also shown. The relative contribution of the latter galaxy population increases with decreasing wavelength and becomes dominant at 100 μm. This trend implies a decrease of the level of correlations among the maps with increasing separation in wavelength. As illustrated by Figure 17, the model is in very good agreement with the cross-wavelength correlations measured by Viero et al. (2012) and defined by (Equation (14) of Viero et al. 2012)

$$\begin{equation} C_{A\times B} = {P^{A\times B}_{k_\theta } \over \sqrt{P^A_{k_\theta }\cdot P^B_{k_\theta }}}. \end{equation} \tag{ 22 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Studies of galaxy properties as a function of morphological type (e.g., Bernardi et al. 2010) have highlighted a dichotomy between the luminosity-weighted ages of early- and late-type galaxies. The former are mostly older than 8 Gyr while most of Sb or later-type spirals are younger than 7 Gyr, corresponding to a formation redshift z ≲ 1–1.5. Building on this datum, we have worked out a model whereby the proto-spheroidal galaxies, in the process of forming the bulk of their stars, are the dominant population in the IR at z ≳ 1.5 while late-type galaxies dominate at lower redshifts. The model is "hybrid" in the sense that it combines a physical, forward model for spheroidal galaxies and the early evolution of the associated AGNs with a phenomenological backward model for late-type galaxies and for the later AGN evolution.

To describe the cosmological evolution of proto-spheroidal galaxies and of the associated AGNs, we adopted the physical model by Granato et al. (2004), upgraded working out, for the first time, the epoch-dependent LFs of sources as a whole (stellar plus AGN component), taking into account in a self-consistent way the variation with galactic age of the global SED. With only minor adjustments of the parameters, the model accurately reproduces the observed LFs at all redshifts (z ≳ 1.5) and IR wavelengths at which they have been determined. The model naturally accounts for the observed positive evolution of both galaxies and AGNs up to z ≃ 2.5 and for the negative evolution at higher redshifts. This is the result of the combination of two competing effects. On one side cooling and free-fall timescales shorten with increasing redshift because of the increase of the matter density and this yields higher SFRs, i.e., higher galaxy luminosities at given halo mass. The higher gas densities are also responsible for a delay of the AGN switch-off time by feedback implying positive luminosity and density evolution of these objects. These effects are thwarted by the decrease in the comoving density of massive halos that prevails above z ≃ 2.5 causing a decline of the bolometric LFs of both galaxies and AGNs.

The model also provides a simple physical explanation of the steeply rising portion of the (sub-)millimeter counts, which proved to be very hard to account for by other physical and phenomenological models. The sharp steepening is due to the sudden appearance of proto-spheroidal galaxies that do not have, in this spectral band, an evolutionary connection with nearby galaxies because their descendants are in passive evolution at z ≲ 1.5. Their (sub-)millimeter counts are extremely steep because, due to the strongly negative K-correction, the flux densities corresponding to a given luminosity are only weakly dependent on the source redshift. Then, since the far-IR luminosity is roughly proportional to the halo mass, the counts reflect, to some extent, the exponential decline of halo mass function at high masses.

The steepness of the counts implies a strong magnification bias due to gravitational lensing. The counts of strongly lensed sources depend on the redshift distribution that determines the distribution of lensing optical depths. In fact, this model was the only one that correctly predicted (Negrello et al. 2007) the strongly lensed counts at 500 μm and the correct redshift distribution of bright (S_500 μm ⩾ 100 mJy) submillimeter sources (Negrello et al. 2010; González-Nuevo et al. 2012).

The epoch-dependent LF of late-type galaxies has been modeled in terms of two populations, "warm" and "cold" galaxies with different SEDs and different evolution properties. Simple truncated power-law models have been adopted for the evolution of these populations. "Cold" (normal) late-type galaxies evolve (weakly) only in luminosity, while "warm" (starburst) galaxies evolve both in luminosity and in density.

Below z = 1.5, the far-IR emission of proto-spheroidal galaxies and the associated AGNs fade out rather rapidly. The AGNs, however, can be reactivated, e.g., by interactions. This later phase of AGN emission has been described by a phenomenological model analogous to that used for late-type galaxies, distinguishing between type 1 and type 2 AGNs.

In this framework, there is a systematic variation with wavelength of the populations dominating the counts and the contributions to the extragalactic background intensity. Above 350 μm the main contributors to the CIB are proto-spheroidal galaxies. In this wavelength range, late-type galaxies dominate the counts only at the brightest (where normal "cold" star-forming galaxies prevail) and at the faintest flux densities (where "warm" starburst galaxies outnumber the proto-spheroids). But these galaxies become increasingly important with decreasing wavelength. Proto-spheroids are always sub-dominant below 250 μm. This strong variation with wavelength in the composition of IR sources implies specific predictions for the auto- and cross-power spectra of the source distribution, which may help discriminate between different models. Essentially, all the alternative models have all source populations present over the full relevant redshift range. This implies a high correlation between the CIB intensity fluctuations at different frequencies. On the contrary, the present model predicts a high (close to unity) cross-correlation only at the longest wavelengths (⩾500 μm). At shorter wavelengths the cross-correlation progressively weakens and we expect little cross-correlation between CIB fluctuations at, say, 100 and 500 μm. No observational determination is available for correlations among these wavelengths, but in the Herschel/SPIRE wavelength range, where cross-correlations have been measured, the model results are in good agreement with observations.

According to our model, the AGN contribution to the CIB is always sub-dominant. It is maximal in the mid-IR, where it reaches 8.6% at 16 μm and 8.1% at 24 μm. These contributions are close to, but somewhat lower than, most observation-based estimates which, however, are complicated by the difficulty of separating the AGN emission from that of the host galaxy. The AGN contribution to the counts is also always sub-dominant. We find a maximum contribution in the mid-IR where the model gives AGN fractions in fair agreement with the observational estimates (Treister et al. 2006; Teplitz et al. 2011).

We are indebted to Matthieu Béthermin for several useful clarifications on flux calibration and color correction issues, to Roberto Assef for having sent his IR luminosity functions of AGNs in tabular form, and to Aurelie Pénin for having provided a tabulation of CIB power spectra at 100 and 160 μm and clarifications on error estimates. We also benefited from useful comments from an anonymous referee. Z.Y.C. acknowledges support from the joint PhD project between XMU and SISSA. A.L. thanks SISSA for warm hospitality. The work has been supported in part by ASI/INAF agreement No. I/072/09/0 and by INAF through the PRIN 2009 "New light on the early universe with submillimeter spectroscopy." J.Q.X. is supported by the National Youth Thousand Talents Program and grants No. Y25155E0U1 and No. Y3291740S3.

The gas initially associated with a galactic halo of mass M_vir, with a cosmological mass fraction f_b = M_gas/M_vir = 0.165 is heated to the virial temperature at the virialization redshift, z_vir. Its subsequent evolution partitions it in three phases: a hot diffuse medium with mass M_inf infalling and/or cooling toward the center; cold gas with mass M_cold condensing into stars; low angular momentum gas with mass M_res stored in a reservoir around the central supermassive black hole, and eventually viscously accreting onto it. In addition, two condensed phases appear and grow, namely, stars with a total mass M_⋆ and the black hole with mass M_•. As mentioned in Section 2.1, we restrict ourselves to the ranges 11.3 ≲ log (M_vir/M_☉) ≲ 13.3 and z_vir ≳ 1.5.

The evolution of the three gas phases is governed by the following equations:

$$\begin{eqnarray} \dot{M}_{\rm inf} &=& - \dot{M}_{\rm cond} - \dot{M}_{\rm inf}^{\rm QSO},\nonumber \\ \dot{M}_{\rm cold} &=& \dot{M}_{\rm cond} - [1-{\cal {R}}(t)]\dot{M}_{\star } - \dot{M}_{\rm cold}^{\rm SN} - \dot{M}_{\rm cold}^{\rm QSO},\\ \dot{M}_{\rm res} &=& \dot{M}_{\rm inflow} - \dot{M}_{\rm BH},\nonumber \end{eqnarray} \tag{ A1 }$$

which link the mass infall rate, $\dot{M}_{\rm inf}$, the variation of the cold gas mass, $\dot{M}_{\rm cold}$, and the variation of the reservoir mass, $\dot{M}_{\rm res}$, to the condensation rate of the cold gas, $\dot{M}_{\rm cond}$, to the SFR $\dot{M}_\star$, to the cold gas removal by SN and AGN feedback, $\dot{M}_{\rm cold}^{\rm SN}$ and $\dot{M}_{\rm cold}^{\rm QSO}$, respectively, to the fraction of gas restituted to the cold component by the evolved stars, $\mathcal {R}(t)$, to the inflow rate of cold gas into the reservoir around the central supermassive black hole, $\dot{M}_{\rm inflow}$, and to the back hole accretion rate, $\dot{M}_{\rm BH}$.

The hot gas cools and flows toward the central region at a rate

$$\begin{equation} \dot{M}_{\rm cond} \simeq \frac{M_{\rm inf}}{t_{\rm cond}}, \end{equation} \tag{ A2 }$$

with $M_{\rm inf}^{0}=f_{\rm b} M_{\rm vir}$ and

$$\begin{equation} t_{\rm cond} \simeq 8\times 10^8 \left (\frac{1+z}{4}\right)^{-1.5} \left (\frac{M_{\rm vir}}{10^{12}M_{\odot }}\right)^{0.2}\ {\rm yr}, \end{equation} \tag{ A3 }$$

where the coefficient is 10% smaller than the value used by Fan et al. (2010). Note that the cooling and inflowing gas we are dealing with is the one already present within the halo at virialization. In this respect it is useful to keep in mind that the virial radius of the halo (R_vir ≃ 220(M_vir/10¹³ M_☉)^1/3[3/(1 + z_vir)] kpc) is more than 30 times larger than the size of the luminous galaxy, and that only a minor fraction of the gas within the halo condenses into stars. Indeed, we need strong feedback processes, capable of removing most of the halo gas, to avoid an overproduction of stars. This implies that any gas infalling from outside the halo must also be swept out by feedback; it could, however, become important for the formation of a disk-like structure surrounding the preformed spheroid once it enters the passive evolution phase, with little feedback (Cook et al. 2009). As mentioned in Section 2.1, the additional material (stars, gas, dark matter) infalling after the fast collapse phase that creates the potential well, i.e., during the slow-accretion phase, mostly produces a growth of the halo outskirts, and has little effect on the inner part where the visible galaxy resides.

The SFR is given by

$$\begin{equation} \dot{M}_\star \simeq \frac{M_{\rm cold}}{t_\star }, \end{equation} \tag{ A4 }$$

where the star formation timescale is t_⋆ ≃ t_cond/s with s ≃ 5. For a Chabrier (2003) IMF of the form ϕ(m) = m^−x with x = 1.4 for 0.1 ⩽ m ⩽ 1M_☉ and x = 2.35 for m > 1M_☉, we find ${\cal R} \simeq 0.54$ under the instantaneous recycling approximation.

The gas mass loss due to the SN feedback is

$$\begin{equation} \dot{M}_{\rm cold}^{\rm SN} = \beta _{\rm SN} \dot{M}_\star , \end{equation} \tag{ A5 }$$

with

$$\begin{eqnarray} \beta _{\rm SN} = &{\frac{N_{\rm SN} \epsilon _{\rm SN} E_{\rm SN}}{E_{\rm bind}}} \simeq 0.6{\left(\frac{N_{\rm SN}}{8\times 10^{-3}/M_\odot }\right)}{\left(\frac{\epsilon _{\rm SN}}{0.05}\right)} \nonumber \\ &\times {\left(\frac{E_{\rm SN}}{10^{51}\ {\rm erg}}\right)}{\left(\frac{M_{\rm vir}}{10^{12}\ M_\odot }\right)^{-2/3}}{\left(\frac{1+z}{4}\right)^{-1}}. \end{eqnarray} \tag{ A6 }$$

We adopt the following values: number of SNe per unit solar mass of condensed stars N_SN ≃ 1.4 × 10⁻²/M_☉; fraction of the released energy used to heat the gas _SN = 0.05; kinetic energy released per SN E_SN ≃ 10⁵¹ erg; halo binding energy E_bind ≃ 3.2 × 10¹⁴(M_vir/10¹² M_☉)^2/3([(1 + z)/4] cm² s⁻² (Mo & Mao 2004).

The infrared luminosity (8–1000 μm) associated with dust enshrouded star formation is

$$\begin{equation} L_{\star ,\rm IR}(t) = k_{\star ,\rm IR} \times 10^{43} \left (\frac{\dot{M}_\star }{M_\odot \ {\rm yr}^{-1}}\right)\ {\rm erg\ s}^{-1}, \end{equation} \tag{ A7 }$$

where the coefficient k_{⋆, IR} depends on the SED. We adopt k_{⋆, IR} ∼ 3 (Lapi et al. 2011; Kennicutt 1998).

The cold gas inflow rate into the reservoir around the supermassive black hole, driven by radiation drag, is given by

$$\begin{equation} \dot{M}_{\rm inflow} \simeq \frac{L_\star }{c^2}(1-e^{-\tau _{\rm RD}}) \simeq \alpha _{\rm RD} \times 10^{-3} \dot{M}_\star (1-e^{-\tau _{\rm RD}}), \end{equation} \tag{ A8 }$$

with

$$\begin{equation} \tau _{\rm RD}(t) = \tau _{\rm RD}^0 \left (\frac{Z_{\rm cold}(t)}{Z_\odot }\right) \left (\frac{M_{\rm cold}}{10^{12}\ M_\odot }\right) \left (\frac{M_{\rm vir}}{10^{13}\ M_\odot }\right)^{-2/3}. \end{equation} \tag{ A9 }$$

For the strength of the radiation drag, we adopt α_RD = 2.5 and set $\tau ^0_{\rm RD} = 3.0$. The model also follows the evolution of the cold gas metallicity, Z_cold(t). An approximate solution of the equations governing the chemical evolution is (A. Lapi et al., in preparation)

$$\begin{eqnarray} && Z_{\rm cold} (t) = Z^0_{\rm inf} + \frac{s}{s \gamma - 1}{\cal E}_Z(t) - \frac{s t / t_{\rm cond}}{e^{(s \gamma -1)t/t_{\rm cond}} -1} \nonumber \\ && \quad \cdot \left \lbrace {\cal E}_Z(t) + B_Z \sum ^{\infty }_{i=2} \frac{1}{i \cdot i!} \left [ (s \gamma -1) \frac{\min (t, t_Z)}{t_{\rm cond}} \right ]^{i-1} \right \rbrace , \nonumber \\ \end{eqnarray} \tag{ A10 }$$

where $\gamma = 1 - {\cal R} - \beta _{\rm SN}$, the metallicity of the primordial infalling gas is $Z^0_{\rm inf} = 10^{-5}$, and the mass fraction of newly formed metals ejected from stars, ${\cal E}_Z(t)$ is given by

$$\begin{equation} \simeq A_Z + B_Z \ln \left [ \frac{\min (t, t_{\rm saturation})}{t_Z} \right ], \end{equation} \tag{ A11 }$$

with A_Z = 0.03, B_Z = 0.02, t_Z = 20 Myr, and t_saturation = 40 Myr for the Chabrier's IMF (Z_☉ ≃ 0.02). Equation (A11) accounts for the fact that, soon after the onset of star formation, the metal yield, mainly contributed by stars with large masses (⩾20 M_☉) and short lifetimes (t_Z ⩽ 20 Myr), is a relatively large fraction of the initial stellar mass (f_metal ⩾ 0.12) while, as the star formation proceeds, it progressively lowers to f_metal ∼ 0.06 as the main contribution shifts to stars with intermediate masses ∼9–20 M_☉ and lifetimes t_Z ∼ 20–40 Myr, and finally saturates to values f_metal ∼ 0.013 as stars with masses ⩽9 M_☉ and long lifetimes (t_saturation ⩾ 40 Myr) take over (Bressan et al. 1998). The two parameters A_Z and B_Z depends mainly on the IMF.

The accretion rate into the central black hole obeys the equation

$$\begin{equation} \dot{M}_{\rm BH} = \min \left(\dot{M}_{\rm BH}^{\rm visc}, \lambda _{\rm Edd} \dot{M}_{\rm Edd}\right), \end{equation} \tag{ A12 }$$

where $\dot{M}_{\rm BH}^{\rm visc}$ is the accretion rate allowed by the viscous dissipation of the angular momentum of the gas in the reservoir

$$\begin{eqnarray} \dot{M}_{\rm BH}^{\rm visc} &=& \frac{M_{\rm res}}{\tau _{\rm visc}} = \kappa _{\rm accr}\ 5\times 10^3{\left(\frac{V_{\rm vir}}{500\ {\rm km\ s}^{-1}}\right)^3} \nonumber \\ && \times {\left(\frac{M_{\rm res}}{M_{\bullet }}\right)^{3/2}}{\left(1+\frac{M_{\bullet }}{M_{\rm res}}\right)^{1/2}}, \end{eqnarray} \tag{ A13 }$$

with κ_accr ≃ 10⁻² and $V_{\rm vir}^2 = G M_{\rm vir}^{2/3} [4\pi \Delta _{\rm vir}(z) \bar{\rho }_{\rm m}(z)/3]^{1/3}$, Δ_vir being the overdensity of a virialized halo at redshift z_vir within its virial radius r_vir. $\dot{M}_{\rm Edd} \equiv M_\bullet /\epsilon \, t_{\rm Edd}$ is the accretion rate corresponding to the Eddington luminosity given the mass-to-light conversion efficiency (we set = 0.1 so that the Salpeter time  t_Edd = 4.5 × 10⁷ yr) and λ_Edd(z) is the Eddington ratio that we assume to slightly increase with redshift for z ≳ 1.5

$$\begin{equation} \lambda _{\rm Edd}(z) \simeq 0.1 (z-1.5)^2 + 1.0 \end{equation} \tag{ A14 }$$

up to a maximum value λ_{Edd, max} = 4. The growth rate of the black hole mass is

$$\begin{equation} \dot{M}_{\bullet }(t) = (1-\epsilon) \dot{M}_{\rm BH} \end{equation} \tag{ A15 }$$

starting from a seed mass $M_{\bullet }^{\rm seed}=10^2\,M_\odot$. The bolometric AGN luminosity is

$$\begin{eqnarray} &&L_{\bullet } = \epsilon \dot{M}_{\rm BH} c^2 = 5.67\times 10^{45} \left (\frac{\epsilon }{0.1}\right) \left (\frac{\dot{M}_{\rm BH}}{M_\odot \ \rm yr^{-1}}\right)\ {\rm erg\ s^{-1}}. \nonumber \\ \end{eqnarray} \tag{ A16 }$$

A minor fraction of it couples with the ISM of the host galaxy giving rise to an outflow at a rate

$$\begin{equation} \dot{M}_{\rm inf,cold}^{\rm QSO} = \dot{M}_{\rm wind} \frac{M_{\rm inf,cold}}{M_{\rm inf}+M_{\rm cold}}, \end{equation} \tag{ A17 }$$

with

$$\begin{equation} \dot{M}_{\rm wind} = \frac{L_{\rm QSO}^{\rm ISM}}{E_{\rm bind}}, \end{equation} \tag{ A18 }$$

and

$$\begin{equation} L_{\rm QSO}^{\rm ISM} \simeq 2\times 10^{44}\epsilon _{\rm QSO}\left (\frac{\dot{M}_{\rm BH}}{M_\odot \ {\rm yr}^{-1}}\right)^{3/2}\ {\rm erg\ s}^{-1}. \end{equation} \tag{ A19 }$$

$L_{\rm QSO}^{\rm ISM}$ is the mechanical AGN luminosity, used to unbind the gas. The coefficient quantifying the strength of the QSO feedback is chosen to be _QSO = 3. The ratio of the mechanical to the total AGN luminosity

$$\begin{equation} L_{\rm QSO}^{\rm ISM}/L_{\bullet } \simeq 3.5\times 10^{-3}\ \frac{\epsilon _{\rm QSO}}{\epsilon } \left(\frac{\dot{M}_{\rm BH}}{M_\odot \ {\rm yr}^{-1}}\right)^{1/2}, \end{equation} \tag{ A20 }$$

is constrained to be in the range 0.006–0.15.

Examples of the resulting evolution with galactic age of properties of the stellar and of the AGN component are shown in Figure 18 for three values of the virial mass and z_vir = 3.

Zoom In Zoom Out Reset image size

As mentioned in Section 5, to improve the fits of the data we have modified, by trial and error, the values of some model parameters used in previous papers, still within their plausible ranges (see Table 2). The impact of these parameters on the derived LFs can be more easily understood with reference to the time lag between the halo virialization and the peak in black hole accretion rate, Δt_peak (Lapi et al. 2006). The duration of star formation is Δt_SF ≲ Δt_peak (see Figure 18) due to the drastic effect of QSO feedback in massive halos which dominate the bright end of the LFs. Note that longer Δt_peak (or Δt_SF) imply higher bright tails of the LFs. The final black hole mass increases with increasing the coefficient, $\tau ^0_{\rm RD}$, of the optical depth of gas clouds (Equation (A9)) because it implies a higher efficiency of the radiation drag driving the gas into the reservoir. There is a degeneracy, to some extent, between $\tau ^0_{\rm RD}$ and the gas metallicity Z_cold, implying that $\tau ^0_{\rm RD}$ cannot be tightly constrained (see Granato et al. 2004). The value of Δt_peak grows substantially in response to a small increase of the radiative efficiency that yields a slower growth of the black hole mass and a weaker QSO feedback. Higher values of the Eddington ratio, λ_Edd, result in lower values of both Δt_peak and of the final black hole mass. A rise of λ_Edd at high-z is required to account for the observed space density of very luminous QSOs (see the high-z data in Figure 7 and Figure 8; Lapi et al. 2006). A higher QSO feedback efficiency (higher _QSO) shortens the duration of star formation, Δt_SF, but has a minor effect on Δt_peak and on the final black hole mass. Finally, the coefficient relating the SFR to the IR luminosity, k_{⋆, IR}, varies with age mix of stellar populations, chemical composition, and IMF. Increasing it, we shift the LFs toward higher luminosities.