## Abstract

We combine a semi-analytic model of galaxy formation with simple analytic recipes describing the absorption and re-emission of starlight by dust in the interstellar medium of galaxies. We use the resulting models to predict galaxy counts and luminosity functions from the far-ultraviolet (FUV) to the submillimetre, from redshift five to the present, and compare with an extensive compilation of observations. We find that in order to reproduce the rest-UV and optical luminosity functions at high redshift, we must assume an evolving normalization in the dust-to-metal ratio, implying that galaxies of a given bolometric luminosity (or metal column density) must be less extinguished than their local counterparts. In our best-fitting model, we find remarkably good agreement with observations from rest ∼1500 Å to m. At longer wavelengths, most dramatically in the submillimetre, our models underpredict the number of bright galaxies by a large factor. The models reproduce the observed total IR luminosity function fairly well. We show the results of varying several ingredients of the models, including various aspects of the dust attenuation recipe, the dust emission templates and the cosmology. We use our models to predict the integrated extragalactic background light, and compare with an observationally motivated extragalactic background light EBL model and with other available observational constraints.

## 1 INTRODUCTION

New observational facilities have greatly extended the range of the electromagnetic spectrum over which emission from galaxies can be measured, while simultaneously expanding the range of cosmic history that can be probed. For example, in recent years, survey observations over significant fractions of the sky have been carried out in the far-ultraviolet (FUV) and near-ultraviolet (NUV) by Galaxy Evolution Explorer ( GALEX ), in the optical by the Sloan Digital Sky Survey (SDSS), in the near-infrared (NIR) by Two Micron All-Sky Survey (2MASS), in the NIR and mid-IR by the Spitzer Space Telescope and, most recently, in the far-IR (FIR) by the Herschel telescope. In addition, multiple ∼0.5–1 deg 2 sized fields have now been deeply imaged in the X-ray with Chandra , in the optical with the Advanced Camera for Surveys (ACS) on the Hubble Space Telescope ( HST ) as well as ground-based facilities, and in the NIR with United Kingdom Infrared Telescope (UKIRT) and Spitzer , allowing large samples of high-redshift galaxies to be identified and studied. The recently installed Wide Field Camera 3 (WFC3) on HST is in the process of observing many of these fields at high resolution in the NIR. Crucial to the extraction of physical quantities and scientific insight from these deep surveys has also been the availability of accurate redshift information for large numbers of galaxies from multi-wavelength medium-band surveys (COMBO-17, COSMOS, MUSYC, NEWFIRM) and multi-object spectroscopy (DEEP, VIMOS).

However, perhaps some of the most surprising and poorly understood observational results of the past two decades have come from long-wavelength observations in the mid- to FIR. The IRAS satellite revealed that ∼30 per cent of the bolometric luminosity of nearby galaxies, mostly normal spirals, is reprocessed by dust in the IR ( Soifer & Neugebauer 1991 ), and discovered a population of heavily obscured luminous and ultra-luminous infrared galaxies (LIRGs and ULIRGS; Sanders & Mirabel 1996 ). IRAS already provided hints of the very strong evolution of this IR-bright population, the number density of which seems to have been enormously larger in the past. This was confirmed and quantified first by Infrared Space Observatory ( ISO ; e.g. Elbaz et al. 1999 , 2002 ), and then by submillimetre Common-user Bolometer Array (SCUBA; Smail, Ivison & Blain 1997 ; Hughes et al. 1998 ; Chapman et al. 2005 ) and Spitzer (e.g. Le Floc’h et al. 2005 ; Babbedge et al. 2006 ). The physical interpretation of these high-redshift LIRGS and ULIRGS, however, has until recently been hampered by the fact that in many cases observations existed only in the mid-IR or submillimetre (sub-mm), relatively far from the m peak of the dust emission (e.g. 15 m in the case of ISO , 24 m for Spitzer , and 450 and 850 m in the case of SCUBA). This situation should improve greatly in the next few years, as observations bracketing the 100 m peak are taken by the PACS (57 to 210 m) and SPIRE (250, 350, 500 m) instruments on the recently launched Herschel telescope.

This rainbow of observations presents a challenge to theoretical models of galaxy formation. To date, most studies have focused on making predictions for rest-optical or intrinsic [e.g. stellar mass, star formation rate (SFR)] properties of galaxies, mainly because of the difficulty of modelling the absorption and emission of light by dust in the interstellar medium (ISM) of galaxies. However, observational estimates of intrinsic physical properties from multi-wavelength photometry suffer from poorly constrained biases ( Lee et al. 2009 ), and quantities such as the SFR can differ greatly depending on which observational tracer is used to estimate them. Therefore, in order to interpret the zoo of galaxies selected at different wavelengths [e.g. LBGs, EROs, DRGs, dust obscured galaxies (DOGs), BzKs, BM/BX, SMGs, etc.], it is important to develop models that can accurately predict observable quantities over the full range of wavelengths probed by modern panchromatic surveys.

Important theoretical advances have been made in the past few years, as well, with the development of radiative transfer (RT) codes that, coupled with a model of the distribution of stars, gas and dust in a galaxy, can produce detailed panchromatic predictions of the galaxy’s appearance and photometric properties ( Silva et al. 1998 ; Jonsson 2006 ; Jonsson et al. 2006 ; Jonsson, Groves & Cox 2010 ). One approach is to use idealized galaxy models (e.g. spheroid plus disc), coupled with an RT code, as in Silva et al. (1998) . Another is to use hydrodynamic simulations to provide more detailed spatial and morphological information, as in Jonsson et al. (2006) .

While powerful, these tools are still computationally quite expensive. Producing panchromatic predictions for statistical samples of galaxies in a cosmological context remains beyond reach, without some shortcuts. Semi-analytic models (SAMs) of galaxy formation, which apply simple but physically motivated recipes for the physical processes that shape galaxy formation, within the framework of structure formation predicted by Λ cold dark matter (ΛCDM), can provide predictions of bulk galaxy properties (such as star formation and chemical enrichment history, radial size, total stellar mass or luminosity, ratio of spheroid to disc, etc.) for very large numbers of galaxies. SAMs have been shown to reproduce many observed properties of galaxies (e.g. Kauffmann, White & Guiderdoni 1993 ; Cole et al. 1994 ; Kauffmann et al. 1998 ; Somerville & Primack 1999 ; Cole et al. 2000 ; Somerville, Primack & Faber 2001 ), and to agree reasonably well with the results of numerical hydrodynamic simulations in their predictions for basic quantities such as the rate of accretion of cold gas and galaxy mergers ( Yoshida et al. 2002 ; Cattaneo et al. 2007 ). In particular, recent models that include ‘radio mode’ feedback from active galactic nuclei (AGN) reproduce quite well the global properties of massive galaxies over a broad range of cosmic history (e.g. S08; Bower et al. 2006 ; Croton et al. 2006 ; Kang, Jing & Silk 2006 ; Menci et al. 2006 ; Monaco, Fontanot & Taffoni 2007 ), although reproducing the properties of low-mass galaxies remains a challenge ( Fontanot et al. 2009a ; Guo et al. 2010 ).

Using a set of recipes for gas accretion and cooling, merging, star formation, stellar feedback, chemical enrichment, and optionally black hole (BH) growth and AGN feedback, an SAM outputs a distribution of ages and metallicities for all the stars within the spheroid and disc components of a galaxy (more details are given in Section 2 ). This information is convolved with ‘simple stellar population’ (SSP) models (e.g. Bruzual & Charlot 2003 ), which specify, for a given stellar initial mass function (IMF), the luminosity as a function of wavelength for a stellar population of a given age and metallicity, in order to predict the unattenuated spectral energy distribution (SED) of starlight in the galaxy. The predictions of stellar population models from different groups have largely converged in their predictions, particularly in the UV and optical, making this component of the modelling relatively robust. 1

For the more difficult step of predicting how this starlight is absorbed and re-radiated by dust, one possible approach is to couple the predictions of a SAM directly with an RT code (e.g. Granato et al. 2000 ; Baugh et al. 2005 ; Fontanot et al. 2007 ). Because SAMs are not able to track the detailed internal structure or morphology of galaxies, this requires the assumption of an idealized geometry such as a spheroid plus disc (where the sizes and masses of the components are specified by the SAM). However, even this approach is prohibitively expensive for large numbers of galaxies, and also has the disadvantage that the simplified geometries may not be representative of the diversity of galaxy types, particularly for LIRG or ULIRG-like objects, many of which are known to be merging systems. Moreover, the dust models contain a large number of free parameters, which must be tuned by fitting a chosen set of observations, and may or may not be constant from galaxy to galaxy or over cosmic time.

An alternative approach is to develop an analytic or SAM to estimate the fraction of starlight that is absorbed by dust in a given galaxy, based on its geometry, metal content and stellar populations. Early SAMs ( Guiderdoni & Rocca-Volmerange 1987 ; Lacey et al. 1993 ; Guiderdoni et al. 1998 ; Kauffmann et al. 1998 ; Somerville & Primack 1999 ; Devriendt & Guiderdoni 2000 ) approached this by assuming that the face-on B- or V -band optical depth of the dust in the disc is proportional to the column density of metals in the gas phase, that the inclination dependence is that predicted by a simple ‘slab’ model in which the dust and stars are uniformly mixed, and that the wavelength dependence is given by a fixed ‘attenuation law’, such as a Galactic extinction law or the starburst attenuation law of Calzetti et al. (2000) . Charlot & Fall (2000) proposed a two-component model that separately accounts for the extinction due to diffuse ‘cirrus’ in the disc and that due to the dense ‘birth clouds’ surrounding newly born stars. De Lucia & Blaizot (2007) combined the two approaches, using a ‘slab’ model to treat the cirrus component and adopting the Charlot & Fall (2000) model to treat the young stars (≲10 7 yr). Fontanot et al. (2009b) tested a wide range of such simple analytic approaches from the literature against full RT, applied within the MORGANA ( Monaco et al. 2007 ) SAM. They concluded that bulk properties [such as UV, optical and NIR luminosity functions (LFs)] predicted by the SAMs using analytic dust recipes agreed quite well with the results of the full RT, at a fraction of the computational cost.

With an estimate for the total energy absorbed by dust in hand, and assuming that all of this energy is re-radiated in the IR, one can use observationally derived or observationally calibrated template SEDs describing the wavelength dependence of the dust emission ( Devriendt, Guiderdoni & Sadat 1999 , hereafter DGS99; Chary & Elbaz 2001 ; Dale & Helou 2002 ; Lagache et al. 2004 ; Rieke et al. 2009 , hereafter R09) or modified Planck functions ( Kaviani, Haehnelt & Kauffmann 2003 ) to compute IR luminosities ( Guiderdoni et al. 1998 ; Devriendt & Guiderdoni 2000 ; Hatton et al. 2003 ; Blaizot et al. 2004 ). Observationally, it is known that the mid- to FIR colours (i.e. the ratio of warm to cool dust) are correlated with the total IR luminosity of the galaxy ( Sanders & Mirabel 1996 ). Accordingly, models based on this approach use an SED library indexed by the total IR luminosity, i.e. the total IR luminosity of the model galaxy is used to select the appropriate FIR template. Fontanot & Somerville (2010) compared this kind of approach, again applied to the MORGANA SAM, with the results of coupling the same SAM with the full RT model GRASIL ( Silva et al. 1998 ). Again, the agreement for statistical quantities such as LFs was quite good.

Our goal in this paper is to develop fully SAMs of galaxy formation that can predict photometric properties of galaxies from the FUV to the FIR with reasonable accuracy. To do this, we adopt a modified version of the Charlot & Fall (2000) model to estimate how much light is absorbed by dust in each galaxy, assume that all of this energy is re-radiated by dust, and employ observationally calibrated templates to estimate luminosities in the mid- to FIR. We first ask how successfully this simple approach can reproduce a compilation of observations spanning the FIR to the FUV and a redshift range of zero to five. We expect such a naïve approach to fail in some respects, and we attempt to identify the physical lessons that we can draw from the points of failure. To aid in this, we vary some of the more uncertain ingredients of our models to study the effect on the observables. The resulting fiducial models will be used to create mock catalogues for pan-chromatic surveys such as CANDELS ( Grogin et al. 2011 ; Koekemoer et al. 2011 ), and to make predictions of the extragalactic background light (EBL) and its build-up over cosmic time. The implications for gamma-ray observations will be explored in a companion paper ( Gilmore et al. 2011 , hereafter GSPD).

This paper is structured as follows. In Section 2 we describe the ingredients of the SAM, including our treatment of dust attenuation and emission. In Section 3 we present predictions for LFs, counts and related quantities from the FUV to the sub-mm for several model variants, and compare them with an extensive compilation of observations. We discuss and summarize our results in Section 4 .

## 2 MODELS

### 2.1 Semi-analytic models

#### 2.1.1 Galaxy formation

The SAMs used here have been described in detail in Somerville & Primack (1999) , Somerville et al. (2001) and especially S08, and we refer the reader to those papers for details. Here we provide a brief summary of the basic ingredients of the SAM, which include the growth of structure in the dark matter component in a hierarchical clustering framework, radiative cooling of gas, star formation, supernova feedback, AGN feedback, galaxy merging within dark matter haloes, metal enrichment of the ISM and intracluster medium (ICM), and the evolution of stellar populations.

We assume a standard ΛCDM universe and a Chabrier IMF ( Chabrier 2003 ). We consider two sets of cosmological parameters here. One is the ‘concordance’ cosmology (C-ΛCDM), = 0.3, Ω Λ = 0.7, H0 = 70.0 and σ 8 = 0.9, which was used by S08 and has also been used by a large number of other studies in the literature. We also present results for an updated set of cosmological parameters that are consistent with the 5 year Wilkinson Microwave Anisotropy Probe ( WMAP ) results ( WMAP 5): = 0.28, Ω Λ = 0.72, H0 = 70.0, σ 8 = 0.81 and ns = 0.96 ( Komatsu et al. 2009 ). We note that these values are generally consistent with those obtained from the analysis of the seven-year WMAP data release ( Komatsu et al. 2010 ). The adopted baryon fraction is 0.1658. We consider WMAP 5 to be our ‘fiducial’ model and present results from C-ΛCDM for ease of comparison with the previous work.

The merging histories (or merger trees) of dark matter haloes are constructed based on the extended Press–Schechter formalism using the method described in Somerville & Kolatt (1999) , with improvements described in S08. These merger trees record the growth of dark matter haloes via merging and accretion, with each ‘branch’ representing a merger of two or more haloes. We construct grids of ‘root haloes’ spanning the range Vvir = 30–1200 km s −1 , where Vvir is the circular velocity at the virial radius, and resolve the merger history of each root halo down to progenitors at least 0.01 times the mass of the root. For root haloes more massive than 10 12 M , we follow the merger histories down to a minimum progenitor mass of 10 10 M . We have checked that our results are robust to the chosen mass resolution of the trees.

Whenever dark matter haloes merge, the central galaxy of the largest progenitor becomes the new central galaxy, and all others become ‘satellites’. Satellite galaxies lose angular momentum due to dynamical friction as they orbit and may eventually merge with the central galaxy. To estimate this merger time-scale we use a variant of the Chandrasekhar formula from Boylan-Kolchin, Ma & Quataert (2008) . Tidal stripping and destruction of satellites are also included as described in S08. We have checked that the resulting mass function and radial distribution of satellites (sub-haloes) agree with the results of high-resolution N -body simulations that explicitly follow sub-structure ( Macciò et al. 2010 ).

Before the Universe is reionized, each halo contains a mass of hot gas equal to the universal baryon fraction times the virial mass of the halo. After reionization, the photoionizing background suppresses the collapse of gas into low-mass haloes. We use the results of Gnedin (2000) and Kravtsov, Gnedin & Klypin (2004) to model the fraction of baryons that can collapse into haloes of a given mass after reionization, assuming that the universe was fully reionized by z = 11.

When a dark matter halo collapses, or experiences a merger that at least doubles the mass of the largest progenitor, the hot gas is shock-heated to the virial temperature of the new halo. This radiating gas then gradually cools and collapses. The cooling rate is estimated using a simple spherically symmetric model, based on the following picture. Assuming that the density profile of the gas decreases monotonically with increasing radius, and the cooling rate is more rapid for dense gas, at any moment we can define the ‘cooling radius’ as the radius within which all the gas will have had time to cool within a time tcool . Then, assuming that the initial density profile of the gas is a singular isothermal sphere (ρ gasr−2 ), the cooling rate is given by

1
where mhot is the mass of the hot halo gas, rvir is the virial radius of the dark matter halo and rcool is the cooling radius. We calculate the cooling radius using the metallicity-dependent atomic cooling curves of Sutherland & Dopita (1993) . Previous studies have used different values for the cooling time tcool (e.g. the Hubble time, the time since the last halo merger, or the halo dynamical time). In the models of S08, and also here, it is assumed to be equal to the halo dynamical time, tdynrvir / Vvir , where Vvir is the virial velocity of the halo.

In some cases the cooling radius can be formally larger than the virial radius. In this case, the cooling rate is limited by the infall rate:

2
The cooling radius limited regime ( rcool < rvir ) is often associated with what have been termed ‘hot flows’ in hydrodynamic simulations, in which gas is shock heated in a diffuse halo and then cools. The infall limited cooling regime ( rcool > rvir ) is associated with ‘cold flows’, in which gas streams into the halo along dense filaments, without ever getting heated ( Birnboim & Dekel 2003 ; Dekel & Birnboim 2006 ; Kereš et al. 2005 ).

In the present SAM, we assume that the cold gas is accreted only by the central galaxy of the halo, but in reality satellite galaxies should also receive some measure of new cold gas. In addition, we assume that all newly cooling gas initially collapses to form a rotationally supported disc. The scale radius of the disc is computed based on the initial angular momentum of the gas and the halo profile, assuming that angular momentum is conserved and that the self-gravity of the collapsing baryons causes contraction of the matter in the inner part of the halo ( Blumenthal et al. 1986 ; Flores et al. 1993 ; Mo, Mao & White 1998 ). This approach has been shown to reproduce the observed size versus stellar mass relation for discs from z ∼ 0 to 2 ( Somerville et al. 2008b ).

Star formation occurs in two modes, a ‘quiescent’ mode in isolated discs, and a merger-driven ‘starburst’ mode. Star formation in isolated discs is modelled using the empirical Schmidt–Kennicutt relation ( Kennicutt 1998 ), assuming that only gas above a fixed critical surface density is eligible to form stars. The efficiency and time-scale of the merger-driven burst mode are a function of the merger mass ratio and the gas fractions of the progenitors, and is based on the results of hydrodynamic simulations ( Robertson et al. 2006a ; Hopkins et al. 2009a ).

Some of the energy from supernovae and massive stars is assumed to be deposited in the ISM, resulting in the driving of a large-scale outflow of cold gas from the galaxy. The mass outflow rate is proportional to the SFR and inversely proportional to the galaxy circular velocity (escape velocity) to the power of α rh , where α rh ≃ 2, as expected for ‘energy-driven’ winds. Some fraction of this ejected gas escapes from the potential of the dark matter halo, while some is deposited in the hot gas reservoir within the halo, where it becomes eligible to cool again. The fraction of gas that is ejected from the disc but retained in the halo versus ejected from the disc and halo is a function of the halo circular velocity (see S08 for details), such that low-mass haloes lose a larger fraction of their gas.

Each generation of stars also produces heavy elements, and chemical enrichment is modelled in a simple manner using the instantaneous recycling approximation. For each parcel of new stars d m* , we also create a mass of metals d MZ = y d m* , which we assume to be instantaneously mixed with the cold gas in the disc. The yield y is assumed to be constant, and is treated as a free parameter. When gas is removed from the disc by supernova-driven winds as described above, a corresponding proportion of metals is also removed and deposited either in the hot gas or outside the halo, following the same proportions as the ejected gas.

Mergers are assumed to remove angular momentum from the disc stars and to build up a spheroid. The efficiency of disc destruction and spheroid growth is a function of progenitor gas fraction and merger mass ratio, and is parametrized based on hydrodynamic simulations of disc–disc mergers ( Hopkins et al. 2009a ). These simulations indicate that more ‘major’ (closer to equal-mass ratio) and more gas-poor mergers are more efficient at removing angular momentum, destroying discs and building spheroids. Note that the treatment of spheroid formation in mergers used here has been updated relative to S08 as described in Hopkins et al. (2009b) . The updated model produces good agreement with the observed fraction of disc versus spheroid dominated galaxies as a function of stellar mass.

In addition, mergers drive gas into galactic nuclei, fuelling BH growth. Every galaxy is born with a small ‘seed’ BH (typically in our standard models). Following a merger, any pre-existing BHs are assumed to merge fairly quickly, and the resulting hole grows at its Eddington rate until the energy being deposited into the ISM in the central region of the galaxy is sufficient to significantly offset and eventually halt accretion via a pressure-driven outflow. This results in self-regulated accretion that leaves behind BHs that naturally obey the observed correlation between BH mass and spheroid mass or velocity dispersion ( Di Matteo, Springel & Hernquist 2005 ; Robertson et al. 2006b ; S08).

There is a second mode of BH growth, termed ‘radio mode’, that is thought to be associated with powerful jets observed at radio frequencies. In contrast to the merger-triggered mode of BH growth described above (sometimes called ‘bright mode’ or ‘quasar mode’), in which the BH accretion is fuelled by cold gas in the nucleus, here, hot halo gas is assumed to be accreted according to the Bondi–Hoyle approximation ( Bondi 1952 ). This leads to accretion rates that are typically only about ≲10 −3 times the Eddington rate, so that most of the BH’s mass is acquired during episodes of ‘bright mode’ accretion. However, the radio jets are assumed to couple very efficiently with the hot halo gas, and to provide a heating term that can partially or completely offset cooling during the ‘hot flow’ mode (we assume that the jets cannot couple efficiently to the cold, dense gas in the infall-limited or cold flow regime.). This ‘radio mode feedback’ appears to be able to quite successfully solve many of the problems experienced by previous generations of ΛCDM-based galaxy formation models; for example, the excess number density and overly high specific SFRs and blue colours of massive galaxies ( Croton et al. 2006 ; Bower et al. 2006 ; S08).

For each galaxy, we store a two-dimensional grid of the mass of stars with a given age and metallicity, separately for the disc and spheroid. We then convolve this distribution with the predictions of the SSP models of Bruzual & Charlot (2003) to obtain the SED of the unattenuated starlight. We use the models based on the Padova (1994) isochrones with a Chabrier (2003) IMF.

#### 2.1.2 Dust attenuation

Our model for dust extinction is based on the model proposed by Charlot & Fall (2000) . As in their model, we consider extinction by two components, one due to the diffuse dust in the disc and another associated with the dense ‘birth clouds’ surrounding young star-forming regions. The V -band, face-on extinction optical depth of the diffuse dust is given by

3
where τ dust, 0 is a free parameter, Zcold is the metallicity of the cold gas, mcold is the mass of the cold gas in the disc and rgas is the radius of the cold gas disc, which is assumed to be a fixed multiple of the stellar scale length (see S08). To compute the actual extinction we assign a random inclination to each disc galaxy and use a standard ‘slab’ model, i.e. the extinction in the V band for a galaxy with inclination i is given by
4

Additionally, stars younger than tBC are enshrouded in a cloud of dust with optical depth τ BC, VBC τ V , 0 , where we treat tBC and μ BC as free parameters. Finally, to extend the extinction estimate to other wavebands, we assume a starburst attenuation curve ( Calzetti et al. 2000 ) for the diffuse dust component and a power-law extinction curve Aλ ∝ (λ/5500 Å) n , with n = 0.7, for the birth clouds ( Charlot & Fall 2000 ).

#### 2.1.3 Dust emission

Using the approach described above, we can compute the total fraction of the energy emitted by stars that is absorbed by dust, over all wavelengths, for each galaxy. We then assume that all of this absorbed energy is re-radiated in the infrared (we neglect scattering), and thereby compute the total IR luminosity of the galaxy LIR . We make use of dust emission templates to determine the SED of the dust emission, based on the hypothesis that the shape of the dust SED is well-correlated with LIR . The underlying physical notion is that the distribution of dust temperatures is set by the intensity of the local radiation field; thus more luminous or actively star-forming galaxies should have a larger proportion of warm dust, as is in fact observed ( Sanders & Mirabel 1996 ).

There are two basic kinds of approaches for constructing these sorts of templates. The first is to use a dust model along with either numerical or analytic solutions to the standard RT equations to create a library of templates, calibrated by comparison with local prototypes. This approach was pioneered by Desert, Boulanger & Puget (1990) , and has been followed by many other workers (e.g. DGS99; Guiderdoni et al. 1998 ; Devriendt & Guiderdoni 2000 ). Desert et al. (1990) posited three main sources of dust emission: polycyclic aromatic hydrocarbons (PAHs), very small grains and big grains. The grains are composed of graphite and silicates, with small and big grains probably dominated by graphite and silicate, respectively. The thermal properties of each species are determined by the size distribution and thermal state. Big grains are assumed to be in near thermal equilibrium, and their emission can be modelled as a modified blackbody spectrum. However, small grains and PAHs are probably in a state that is intermediate between thermal equilibrium and single photon heating. They are therefore subject to temperature fluctuations and their emission spectra are much broader than a modified blackbody spectrum. The detailed size distributions are modelled using free parameters, which are calibrated by requiring the model to fit a set of observational constraints, such as the extinction/attenuation curves, observed IR colours and the IR spectra of local galaxies. Here we make use of the templates derived by DGS99 using a similar approach to Desert et al. (1990) .

The second approach is to make direct use of observed SEDs for a set of prototype galaxies (e.g. Chary & Elbaz 2001 ; Dale & Helou 2002 ; Lagache et al. 2004 ). We also make use of the empirical SED templates recently published by R09. They constructed detailed SEDs from published ISO , IRAS and NICMOS data as well as previously unpublished IRAC, MIPS and IRS observations. They modelled the FIR SEDs assuming a single blackbody with wavelength-dependent emissivity. The R09 library includes 14 SEDs covering the 5.6 × 10 9 L < LIR < 10 13 L range. Examples of the DGS99 and R09 SED templates are shown in Fig. 1 .

Figure 1

Comparison of the dust emission templates of R09 (red) and DGS99 (blue). The four panels show the dust emission templates used in this work for bolometric IR luminosities of 10 10 L , 10 11 L (a LIRG), 10 12 L (a ULIRG) and 10 13 L (an extremely IR-bright ‘Hyper-LIRG’).

Figure 1

Comparison of the dust emission templates of R09 (red) and DGS99 (blue). The four panels show the dust emission templates used in this work for bolometric IR luminosities of 10 10 L , 10 11 L (a LIRG), 10 12 L (a ULIRG) and 10 13 L (an extremely IR-bright ‘Hyper-LIRG’).

The R09 templates have less emission in the PAH and mid-IR regions than those of DGS99, particularly at the brightest luminosities. The R09 templates are also considerably more detailed in their representation of PAH emission. Being observationally based, the shortest wavelengths predicted by the R09 templates are contaminated by emission from direct starlight. We have attempted to remove this component by subtracting from each template the average amount of starlight in the SED for galaxies of that IR luminosity in the local universe. The R09 templates also end abruptly at 5 m, and we have smoothed the transition to the shorter wavelength starlight regime by extrapolating using a power law of slope ∼λ 3 . Our results are not very sensitive to the choice of this power-law slope.

#### 2.1.4 Galaxy formation parameters

The galaxy formation models contain a number of free parameters which are tuned using observational constraints. A full list of the physical parameters in the SAM is given in S08 (table 2). The most important parameters for the quantities presented in this paper are those controlling supernova and AGN feedback; specifically, the efficiency of supernova-driven outflows ( ) and its dependence on galaxy circular velocity (α rh ), and the efficiency of heating by ‘radio mode’ AGN feedback (κ radio ). As in S08, the values of these parameters are adjusted in order to obtain a good match to the observed z ∼ 0 stellar mass function. The low-mass end of the stellar mass function is insensitive to the AGN feedback recipe and is mainly controlled by the supernova feedback recipe, while the reverse is true for the high-mass end. The effective yield used for chemical evolution is fixed by matching the zero-point of the galaxy stellar mass–metallicity relation (see S08). The values of the parameters used in the C-ΛCDM model presented in this paper are identical to those for the C-ΛCDM model in S08, except that we have slightly increased the strength of the supernovae feedback ( and α rh = 2.5), and decreased the efficiency of the radio mode feedback (κ radio = 2.5 × 10 −3 .) For the fiducial WMAP 5 model presented here, we have adjusted the galaxy formation parameters slightly to account for the modified cosmology, but they are quite similar to the parameters for the WMAP 3 model presented in S08.

#### 2.1.5 Dust parameters

We also have three additional parameters that control the dust attenuation in our model: the normalization of the face-on V -band optical depth τ dust, 0 , the opacity of the birthclouds relative to the cirrus component μ BC and the time that newly born stars spend enshrouded in their birthclouds, tBC . We first set τ dust, 0 by matching the normalization of the observed relationship between Ldust / LUV versus bolometric luminosity Lbol , where Ldust is the total luminosity absorbed by dust and re-emitted in the mid- to FIR and LUV is the luminosity in the FUV (∼1500 Å). The predicted median relation and 1σ scatter is shown in Fig. 2 , along with several observational estimates. We also show the predicted relationship between the cold gas metallicity and Ldust / LUV compared with observations in Fig. 2 . We find good agreement with τ dust, 0 = 0.2. With this value, we also obtain good agreement with the observed optical through NIR LFs at z = 0 (see Section 3.1 ).

Figure 2

Left: Galaxy luminosity in the IR relative to the luminosity in the UV versus bolometric luminosity. Open symbols with error bars are observational estimates of this relationship for nearby galaxies from Buat et al. (2007) . The long dashed green line is the observational relation from Bell (2003) , and the dotted line is the observational relation from Xu et al. (2006) . Right: Galaxy luminosity in the IR relative to the UV versus metallicity of the cold gas in the galaxy. The dashed (green) line is the observational relation of Heckman et al. (1998) and the solid green line is the relation of Cortese et al. (2006) . In both panels, solid black squares show the medians for our fiducial model, and dashed black lines show the 16th and 84th percentiles.

Figure 2

Left: Galaxy luminosity in the IR relative to the luminosity in the UV versus bolometric luminosity. Open symbols with error bars are observational estimates of this relationship for nearby galaxies from Buat et al. (2007) . The long dashed green line is the observational relation from Bell (2003) , and the dotted line is the observational relation from Xu et al. (2006) . Right: Galaxy luminosity in the IR relative to the UV versus metallicity of the cold gas in the galaxy. The dashed (green) line is the observational relation of Heckman et al. (1998) and the solid green line is the relation of Cortese et al. (2006) . In both panels, solid black squares show the medians for our fiducial model, and dashed black lines show the 16th and 84th percentiles.

The birthcloud parameters μ BC and tBC mainly control the attenuation of UV light relative to longer wavelengths. At z = 0, the g through K -band LFs are insensitive to the birthcloud parameters, while the FUV through u bands are quite dependent on them. We adjust these parameters in order to match the z = 0 FUV and NUV observed LFs, finding good agreement with μ BC = 4.9 and tBC = 2 × 10 7 yr.

It remains an open question whether the properties of interstellar dust have evolved over cosmic time and, if so, what impact this might have on the appearance of high-redshift galaxies. Recent work has suggested that, at fixed Lbol , galaxies may be less extinguished at high redshift than one would expect if the relation shown in Fig. 2 were constant over all times (e.g. Reddy et al. 2006 , 2010 ). In addition, we find, in agreement with some previous studies ( Lo Faro et al. 2009 ; Guo & White 2009 ), that if we keep the dust parameters fixed at constant values, we underpredict the number of UV-bright galaxies at high redshift. Therefore we introduce an ad hoc redshift dependence into the dust parameters, which we use in our fiducial model and most of its variants: τ dust, 0 ( z ) =τ dust, 0 ( z = 0)/(1 + z ), and both μ BC and tBC scale with z−1 above z = 1. We chose this dependence because it allows us to fit the bright end UV and B -band LFs at all redshifts where they are well constrained observationally. We also show results from a model with dust parameters that do not vary with redshift (‘fixed dust’ model).

#### 2.1.6 Model variants

We vary several of the uncertain ingredients of our models in order to study the sensitivity of our results to these assumptions. As we have already discussed, we consider two cosmologies, the ‘concordance’ (C-ΛCDM) or WMAP 1 cosmology and the currently favoured WMAP 5 cosmology. We consider a model variant in which instead of using the composite dust attenuation model (with cirrus plus birth clouds) based on Charlot & Fall (2000) , we instead use a fixed attenuation curve from Calzetti et al. (2000) . We consider two different sets of empirical dust emission SED templates (R09 and DGS99). As discussed above, we also consider a model in which the dust parameters do not evolve with redshift. This results in five models, which are summarized in Table 1 . This table also specifies the line style that will be used for each model in the plots of the following section. For some quantities, some of the models produce the same predictions, in which case we show a subset of the models.

Table 1

Summary of Models.

 Name Cosmology Dust attenuation Dust emission Dust parameters Designated line style WMAP 5 fiducial WMAP 5 Composite Rieke et al. (2009) Evolving Solid black WMAP 5+fixed dust WMAP 5 Composite Rieke et al. (2009) Fixed Dash–dotted purple WMAP 5+Calzetti WMAP 5 Calzetti Rieke et al. (2009) Fixed Dashed red WMAP 5+DGS WMAP 5 Composite Devriendt et al. (1999) Evolving Long-dashed blue C-ΛCDM C-ΛCDM Composite Rieke et al. (2009) Evolving Dotted black
 Name Cosmology Dust attenuation Dust emission Dust parameters Designated line style WMAP 5 fiducial WMAP 5 Composite Rieke et al. (2009) Evolving Solid black WMAP 5+fixed dust WMAP 5 Composite Rieke et al. (2009) Fixed Dash–dotted purple WMAP 5+Calzetti WMAP 5 Calzetti Rieke et al. (2009) Fixed Dashed red WMAP 5+DGS WMAP 5 Composite Devriendt et al. (1999) Evolving Long-dashed blue C-ΛCDM C-ΛCDM Composite Rieke et al. (2009) Evolving Dotted black

## 3 RESULTS

### 3.1 Luminosity functions and luminosity densities

Although our models are very similar to those presented in S08, the resolution used here is much higher than the resolution used in the simulations presented in S08. In S08, haloes were followed down to a minimum mass of 10 10 M , while here, root haloes at z = 0 are followed down to a minimum progenitor mass of 1.3 × 10 8 M . Therefore, we first investigate whether changing the resolution of the merger trees has any effect on our results. In Fig. 3 we show the two main quantities that we use to normalize our models, the stellar mass function and the gas fraction of disc-dominated galaxies as a function of stellar mass, at z = 0. In S08, we stated that galaxies with stellar masses above ∼10 9 M should be reliably resolved. We see from Fig. 3 that the results above this mass change negligibly when we increase the mass resolution by almost two orders of magnitude. In addition, with no tuning, our model fits the observed stellar mass function extremely well down to the smallest masses where it has been measured. The gas fractions also appear to match observations down to very low mass objects (∼10 7 M ).

Figure 3

Left: stellar mass function at z = 0. The solid line shows the fiducial ( WMAP 5) model; dotted line shows the C-ΛCDM model and the dashed grey line shows the fiducial model at the resolution of the simulations presented in S08. Symbols show the observationally derived stellar mass functions from Baldry, Glazebrook & Driver (2008 , blue squares), Baldry et al. (2011 , purple crosses), Panter et al. (2007 , dark green triangles) and Li & White (2009 , light green hexagons). Note that the quoted values for the observed mass functions at mstar ≲10 8.5 M are likely to be lower limits, due to surface brightness selection effects. Right: gas fraction for central disc-dominated galaxies in the fiducial model (grey shading as in S08). Solid and dashed lines show model median and 16th and 84th percentiles. Large open circles show the observations of Kannappan (2004) . Filled squares show gas fractions from galaxies in the THINGS survey ( Leroy et al. 2008 ), and small filled circles show observations from Baldry et al. (2008) .

Figure 3

Left: stellar mass function at z = 0. The solid line shows the fiducial ( WMAP 5) model; dotted line shows the C-ΛCDM model and the dashed grey line shows the fiducial model at the resolution of the simulations presented in S08. Symbols show the observationally derived stellar mass functions from Baldry, Glazebrook & Driver (2008 , blue squares), Baldry et al. (2011 , purple crosses), Panter et al. (2007 , dark green triangles) and Li & White (2009 , light green hexagons). Note that the quoted values for the observed mass functions at mstar ≲10 8.5 M are likely to be lower limits, due to surface brightness selection effects. Right: gas fraction for central disc-dominated galaxies in the fiducial model (grey shading as in S08). Solid and dashed lines show model median and 16th and 84th percentiles. Large open circles show the observations of Kannappan (2004) . Filled squares show gas fractions from galaxies in the THINGS survey ( Leroy et al. 2008 ), and small filled circles show observations from Baldry et al. (2008) .

Next, as a further orientation to our models, we show in Fig. 4 the global SFR density and stellar mass density across cosmological time. The global star formation histories are similar in the two models below redshift two, but the C-ΛCDM model has more early star formation because of the larger amount of small-scale power. We can see that there is a very small difference between the fiducial model at the resolution of S08 and with the higher resolution. This shows that the resolution adopted by S08 was sufficient to resolve the galaxies that contribute the bulk of the global SFR and stellar mass density from z ∼ 6 to the present.

Figure 4

Left: global star formation history. Right: global stellar mass assembly history. Both panels: solid black lines show the predictions of the WMAP 5 model; dotted lines show the C-ΛCDM model; dark grey dashed lines show the fiducial model at the resolution of the simulations presented by S08. Symbols show a compilation of observational estimates (references given in S08). The long-dashed lines are the observational estimates from Hopkins & Beacom (2006) .

Figure 4

Left: global star formation history. Right: global stellar mass assembly history. Both panels: solid black lines show the predictions of the WMAP 5 model; dotted lines show the C-ΛCDM model; dark grey dashed lines show the fiducial model at the resolution of the simulations presented by S08. Symbols show a compilation of observational estimates (references given in S08). The long-dashed lines are the observational estimates from Hopkins & Beacom (2006) .

As a first test we investigate the LFs at z = 0 from the rest-UV to the IR. In Fig. 5 , we show the FUV and NUV LFs from GALEX , ugriz LFs from SDSS, the K -band LF from 2MASS and the total IR LF (references given in figure caption). Our fiducial model agrees very well with all of these data, other than slightly overpredicting the numbers of the faintest galaxies in the GALEX FUV and NUV bands. We see that the model with a fixed attenuation curve (Calzetti model) is unable to simultaneously reproduce the UV through optical bands, while the composite (modified Charlot & Fall model) does this well. In this model, the young stars that produce most of the UV light are heavily enshrouded in dense birth clouds, resulting in an effectively age-dependent extinction curve.

Figure 5

LFs at z = 0 in the FUV, NUV, SDSS ugriz and K bands, as well as the total IR (8–1000 m). The solid black line is our fiducial WMAP 5 model using the composite dust recipe (Charlot–Fall), and dashed red shows the model with a single component dust model (Calzetti). The dotted black line shows the C-ΛCDM model. The black long-dashed line shows the predictions of the fiducial model without dust attenuation. Data are from Wyder et al. (2005) (blue points, 〈 z 〉= 0.05), Bell et al. (2003) (green points) and Rodighiero et al. (2010) (red points, 〈 z 〉= 0.15). For the total IR panel, the x -axis shows the logarithmic luminosity in solar units, and the y -axis has units of N dex −1 Mpc −3 ; all other axes are as indicated.

Figure 5

LFs at z = 0 in the FUV, NUV, SDSS ugriz and K bands, as well as the total IR (8–1000 m). The solid black line is our fiducial WMAP 5 model using the composite dust recipe (Charlot–Fall), and dashed red shows the model with a single component dust model (Calzetti). The dotted black line shows the C-ΛCDM model. The black long-dashed line shows the predictions of the fiducial model without dust attenuation. Data are from Wyder et al. (2005) (blue points, 〈 z 〉= 0.05), Bell et al. (2003) (green points) and Rodighiero et al. (2010) (red points, 〈 z 〉= 0.15). For the total IR panel, the x -axis shows the logarithmic luminosity in solar units, and the y -axis has units of N dex −1 Mpc −3 ; all other axes are as indicated.

Next we explore the evolution over cosmic time of LFs in several wavebands. All LFs are shown in the rest frame at the indicated redshift. We assume that the observations have been properly completeness corrected and make no attempt to determine whether our model galaxies would in fact obey the secondary selection criteria of any given survey (for example, the colour criteria in Lyman-break selected samples). We defer these sorts of more detailed comparisons to future work. Fig. 6 shows the rest-frame ∼1500 Å LFs from z = 0.5 to z = 5. We see that the models produce more than enough UV luminous galaxies at all redshifts, requiring some extinction even at z ∼ 5. However, the model with fixed dust parameters (normalized at z = 0) systematically underproduces UV-luminous galaxies by a larger and larger factor as redshift increases. This has been found before by other studies using SAMs ( Guo & White 2009 ; Lo Faro et al. 2009 ). Moreover, as already mentioned, there is direct observational evidence that high-redshift galaxies may be less extinguished than their low-redshift counterparts ( Reddy et al. 2010 ). With the simple evolving dust model, we obtain fairly good agreement over all redshifts, although we somewhat overproduce low-luminosity galaxies at high redshift. Obviously, we could have adopted a more complex dust model with luminosity-dependent evolution in the dust parameters to get a better overall fit. Alternatively, this could be an indication that low-mass galaxies are forming too early in the models ( Fontanot et al. 2009a ). Fig. 7 shows a similar comparison in the rest B band. We see a similar discrepancy with our non-evolving dust model to that observed in the UV: we need to adopt lower extinctions in high-redshift galaxies to reproduce the number of luminous galaxies. The discrepancy is smaller, and sets in at higher redshift, but is still pronounced by z ∼ 3. The combination of the rest UV and optical observations is what forced us to adopt evolving parameters in both components of our two-component model (the cirrus and the birthclouds).

Figure 6

LFs in the rest-frame 1500 Å UV band at several non-local redshifts. The solid black line is our fiducial WMAP 5 model using the evolving composite dust recipe (Charlot–Fall), the dash–dotted purple shows the composite dust model with fixed parameters, and dashed red is with the single component dust model (Calzetti). The black long-dashed line shows the predictions of the fiducial model without dust attenuation. Blue squares are data from Arnouts et al. (2005) , violet circles are from Hathi et al. (2010) , red stars are from Reddy et al. (2008) and green circles are from Bouwens et al. (2007) .

Figure 6

LFs in the rest-frame 1500 Å UV band at several non-local redshifts. The solid black line is our fiducial WMAP 5 model using the evolving composite dust recipe (Charlot–Fall), the dash–dotted purple shows the composite dust model with fixed parameters, and dashed red is with the single component dust model (Calzetti). The black long-dashed line shows the predictions of the fiducial model without dust attenuation. Blue squares are data from Arnouts et al. (2005) , violet circles are from Hathi et al. (2010) , red stars are from Reddy et al. (2008) and green circles are from Bouwens et al. (2007) .

Figure 7

Rest-frame LFs in the B band at several redshifts. Line types are the same as in Fig. 6 . Red squares are data from Salimbeni et al. (2008) , and open blue hexes are from the study of Giallongo et al. (2005) . The green shaded regions in the top two panels are the Schechter fits to DEEP survey LFs ( Faber et al. 2007 ), with 1σ errors. Green circles in the two middle panels are from Marchesini et al. (2007) .

Figure 7

Rest-frame LFs in the B band at several redshifts. Line types are the same as in Fig. 6 . Red squares are data from Salimbeni et al. (2008) , and open blue hexes are from the study of Giallongo et al. (2005) . The green shaded regions in the top two panels are the Schechter fits to DEEP survey LFs ( Faber et al. 2007 ), with 1σ errors. Green circles in the two middle panels are from Marchesini et al. (2007) .

Next, in Fig. 8 we show the comparison with the rest-frame K -band LFs in the same redshift bins. Here, interestingly, the models overpredict the number of galaxies at high redshift, particularly below the knee in the LF. This is consistent with the findings of Fontanot et al. (2009a) , who showed that three independently developed SAMs all overproduce low-mass galaxies at high redshift (see also Marchesini et al. 2009 ). The good agreement of our models at the brightest K -band magnitudes also indicates that, contrary to the findings of Henriques et al. (2010) , there is no need to invoke an enhanced contribution to the NIR due to the TP-AGB stellar phase, as in e.g. the stellar population models of Maraston (2005) .

Figure 8

Rest-frame LFs in the K band at several redshifts. Line types are the same as in Fig. 6 . Blue squares are observations from Cirasuolo et al. (2010) , red stars at z = 0.5 are from Pozzetti et al. (2003) , and open green stars at redshifts 1 and 2 are from Caputi et al. (2006) . Note that all observations have been converted to the AB magnitude system.

Figure 8

Rest-frame LFs in the K band at several redshifts. Line types are the same as in Fig. 6 . Blue squares are observations from Cirasuolo et al. (2010) , red stars at z = 0.5 are from Pozzetti et al. (2003) , and open green stars at redshifts 1 and 2 are from Caputi et al. (2006) . Note that all observations have been converted to the AB magnitude system.

We now move from bands that are dominated by light directly emitted by stars to light that has been reprocessed by dust. Fig. 9 shows the LF in the Spitzer IRAC rest 8 m band. Our models produce about the right number of galaxies around the knee in the LF, but underproduce luminous galaxies, especially at high redshift ( z ∼ 1–2). As one can see from Fig. 1 , this portion of the spectrum is very complex, and highly dependent on whether the spectrum is dominated by emission from PAHs, or absorption features such as the deep silicate absorption seen at m. In addition, this part of the spectrum is particularly subject to contamination by obscured AGN. If there are a significant number of buried AGN at z ∼ 2, as has been suggested ( Daddi et al. 2007 ), this could at least partially account for the discrepancy between the observations and our model. However, it would also be unsurprising if our simple unevolving template approach were simply inadequate to accurately model the very complex physics that determine the 8 m flux in galaxies. Indeed, we see a large difference between the predictions using the R09 templates and those using the DGS99 dust emission templates. Very similar remarks apply to the comparison between the Spitzer MIPS 24 m LFs in Fig. 10 , although here the models agree reasonably well with the observations at z = 0.5 and z = 1, and show a somewhat milder discrepancy at z = 2. The differences in the predictions with the two different dust emission templates are also considerably smaller than at 8 m.

Figure 9

Rest-frame LFs in the IRAC 8 m band at several redshifts. Long-dashed blue lines show the predictions of the model with the DGS99 dust emission templates. Other lines show results with the IR templates of R09: solid black is our fiducial WMAP 5 model using our evolving composite dust recipe, dash–dotted purple shows the model with fixed dust parameters, dashed red shows the single component dust (Calzetti) model and dotted black is the C-ΛCDM model. Green stars are data from Dai et al. (2009) . Blue hexagons show the data from Rodighiero et al. (2010) and red stars are measurements from Caputi et al. (2007) .

Figure 9

Rest-frame LFs in the IRAC 8 m band at several redshifts. Long-dashed blue lines show the predictions of the model with the DGS99 dust emission templates. Other lines show results with the IR templates of R09: solid black is our fiducial WMAP 5 model using our evolving composite dust recipe, dash–dotted purple shows the model with fixed dust parameters, dashed red shows the single component dust (Calzetti) model and dotted black is the C-ΛCDM model. Green stars are data from Dai et al. (2009) . Blue hexagons show the data from Rodighiero et al. (2010) and red stars are measurements from Caputi et al. (2007) .

Figure 10

LFs in the rest-frame MIPS 24 m band. Line types are the same as in Fig. 9 . Blue hexagons show the observational data from Rodighiero et al. (2010) , and red stars are from Babbedge et al. (2006) . Green asterisks are from Magnelli et al. (2011) . Orange squares are from Rujopakarn et al. (2010) ; we have interpolated these points from data presented for two different redshift bins.

Figure 10

LFs in the rest-frame MIPS 24 m band. Line types are the same as in Fig. 9 . Blue hexagons show the observational data from Rodighiero et al. (2010) , and red stars are from Babbedge et al. (2006) . Green asterisks are from Magnelli et al. (2011) . Orange squares are from Rujopakarn et al. (2010) ; we have interpolated these points from data presented for two different redshift bins.

Fig. 11 shows the LF at rest 250 m, compared with early results from Herschel. At z = 0.5, the only redshift where the LF at this wavelength has been measured to date, the models overpredict the number of luminous galaxies. The last LF we examine is that in Fig. 12 for the total integrated IR from 8 to m, estimated from multi-wavelength observations. Here we see fairly good agreement between the models and the observational estimates, although with a small deficit of luminous galaxies, particularly at z ∼ 2. However, considering that the true errors on the observations are probably considerably larger than the error bars shown, this level of agreement is encouraging. It is also encouraging that the different model variants produce very similar results for this quantity.

Figure 11

The rest-frame LFs at 250 m. Line types are the same as in Fig. 9 , and results from early Herschel observations ( Dye et al. 2010 ) at z = 0.5 are shown by the blue hexagons.

Figure 11

The rest-frame LFs at 250 m. Line types are the same as in Fig. 9 , and results from early Herschel observations ( Dye et al. 2010 ) at z = 0.5 are shown by the blue hexagons.

Figure 12

LFs for the total IR emission. Line types are the same an in Fig. 9 ; however as the total IR predictions are insensitive to the templates used, the model with the DGS99 templates is not shown. Blue hexagons are the data from Rodighiero et al. (2010) , green stars are from Le Floc’h et al. (2005) and red diamonds are from Caputi et al. (2007) .

Figure 12

LFs for the total IR emission. Line types are the same an in Fig. 9 ; however as the total IR predictions are insensitive to the templates used, the model with the DGS99 templates is not shown. Blue hexagons are the data from Rodighiero et al. (2010) , green stars are from Le Floc’h et al. (2005) and red diamonds are from Caputi et al. (2007) .

Fig. 13 shows the total emissivity from all galaxies as a function of wavelength in our fiducial model, at various redshifts. One can see a slight shift in the wavelength of the peak of the dust emission, as well as the ratio of warm to cold dust with redshift. Because we have assumed non-evolving dust emission templates in this work, this is entirely due to the changing mix of galaxies of different total IR luminosities with time (i.e. the declining contribution of very IR luminous galaxies with decreasing redshift). Fig. 14 shows the integrated luminosity density as a function of redshift in the FUV and NUV, the B , z and K band, and the total IR, as a function of redshift for all the models presented in this work. The model predictions are compared with observational estimates obtained by integrating observed LFs (see figure caption for references). Such observational estimates in the UV and IR bands are often used to estimate the global SFR density (e.g. as shown in Fig. 4 ). It is interesting that while our fiducial model is in good agreement with the observations in the FUV over the whole redshift range, it is somewhat low compared with the total-IR luminosity at z ∼ 1–2. As already discussed, this could be due to the total IR having been overestimated, or could indicate that the SFR in the models is too low. Also interesting to note is that the UV and IR luminosity densities peak at a higher redshift ( z ∼ 3) than the z or K band, because the latter arise from older stellar populations.

Figure 13

Integrated luminosity density as a function of wavelength in our WMAP 5 fiducial model, shown at various slices in redshift. All data shown are measurements in the local universe. Measurements are from GALEX (blue circle), SDSS (red stars; Montero-Dorta & Prada 2009 ), 6dF (cyan squares; Jones et al. 2006 ) and 2MASS (green stars; Cole et al. 2001 ; Bell et al. 2003 ). In the mid- and FIR, the orange squares ( Soifer & Neugebauer 1991 ) and blue stars ( Takeuchi et al. 2001 ) show observations from IRAS and SCUBA.

Figure 13

Integrated luminosity density as a function of wavelength in our WMAP 5 fiducial model, shown at various slices in redshift. All data shown are measurements in the local universe. Measurements are from GALEX (blue circle), SDSS (red stars; Montero-Dorta & Prada 2009 ), 6dF (cyan squares; Jones et al. 2006 ) and 2MASS (green stars; Cole et al. 2001 ; Bell et al. 2003 ). In the mid- and FIR, the orange squares ( Soifer & Neugebauer 1991 ) and blue stars ( Takeuchi et al. 2001 ) show observations from IRAS and SCUBA.

Figure 14

Luminosity density as a function of redshift for various wavelengths, as well as for the total IR luminosity (8–1000 m). Line types are the same as those discussed in previous plots; see captions for Figs 6 and 9 and Table 1 . Additionally, orange lines with alternate short and long dashes are the predictions of the empirical model of Domínguez et al. (2011) , for comparison with our work. Observational data are as follows: 1500 Å: blue squares are from Dahlen et al. (2007) , red stars are from Schiminovich et al. (2005) , green stars are from Bouwens et al. (2007) and orange circles are from Reddy et al. (2008) . The solid purple circle is a local measurement with GALEX by Wyder et al. (2005) . 2800 Å: blue squares and the purple circle are again from Dahlen et al. (2007) and Wyder et al. (2005) , respectively. Red stars are from Gabasch et al. (2006) . B -band: blue squares are from Dahlen et al. (2005) , and are incomplete at higher redshifts. DEEP and COMBO-17 data from Faber et al. (2007) are shown as red stars and open red squares, respectively (these are very similar and difficult to differentiate here). Other data shown are from Wolf et al. (2003) (green star) and Marchesini et al. (2007) (open purple hexes). z -band: local measurements are provided by Montero-Dorta & Prada (2009) (red) and Blanton et al. (2003) (green). Blue squares are from Gabasch et al. (2006) . K band: the local determination is from Kochanek et al. (2001) . High-redshift data are from Barro et al. (2009) (blue squares) and Arnouts et al. (2007) (open red hexagons). Total IR luminosity: observational estimates of the IR emissivity are from Caputi et al. (2007) (open blue pentagons), Reddy et al. (2008) (green circles), Rodighiero et al. (2010) (purple stars) and Le Floc’h et al. (2005) (red crosses).

Figure 14

Luminosity density as a function of redshift for various wavelengths, as well as for the total IR luminosity (8–1000 m). Line types are the same as those discussed in previous plots; see captions for Figs 6 and 9 and Table 1 . Additionally, orange lines with alternate short and long dashes are the predictions of the empirical model of Domínguez et al. (2011) , for comparison with our work. Observational data are as follows: 1500 Å: blue squares are from Dahlen et al. (2007) , red stars are from Schiminovich et al. (2005) , green stars are from Bouwens et al. (2007) and orange circles are from Reddy et al. (2008) . The solid purple circle is a local measurement with GALEX by Wyder et al. (2005) . 2800 Å: blue squares and the purple circle are again from Dahlen et al. (2007) and Wyder et al. (2005) , respectively. Red stars are from Gabasch et al. (2006) . B -band: blue squares are from Dahlen et al. (2005) , and are incomplete at higher redshifts. DEEP and COMBO-17 data from Faber et al. (2007) are shown as red stars and open red squares, respectively (these are very similar and difficult to differentiate here). Other data shown are from Wolf et al. (2003) (green star) and Marchesini et al. (2007) (open purple hexes). z -band: local measurements are provided by Montero-Dorta & Prada (2009) (red) and Blanton et al. (2003) (green). Blue squares are from Gabasch et al. (2006) . K band: the local determination is from Kochanek et al. (2001) . High-redshift data are from Barro et al. (2009) (blue squares) and Arnouts et al. (2007) (open red hexagons). Total IR luminosity: observational estimates of the IR emissivity are from Caputi et al. (2007) (open blue pentagons), Reddy et al. (2008) (green circles), Rodighiero et al. (2010) (purple stars) and Le Floc’h et al. (2005) (red crosses).

### 3.2 Galaxy counts

Figs 15 and 16 show galaxy number counts from the UV through the sub-mm. This is an important cross-check on the LF comparison as many high-redshift samples do not have spectroscopic redshifts available particularly for the fainter objects. Our fiducial model provides quite a good agreement with the data at all UV and optical wavelengths. There is a factor of 2–3 excess of faint galaxies in the models in the optical and NIR, which corresponds to the excess of faint galaxies at z ∼ 1–2 seen in the LF comparison. The agreement at 3.6 and 8 m is also quite good, except below Jy, where the models overpredict faint galaxies. The agreement is also good for MIPS 24 and 70 m, except for 24 m sources between 100 Jy to 1 mJy. The models with the R09 templates underproduce galaxies in this range, while the model with the DGS99 templates performs better. This is because of the much lower fluxes in the 10–20 m range in the R09 templates relative to the DGS99 templates. None of the models reproduces the ‘bump’ seen in the observed 24 m counts between ∼0.1 and 1 mJy. The agreement with the 70 m counts is fairly good, but the observed counts do not go as deep.

Figure 15

Number counts in the GALEX UV bands and the four HST ACS bands. Line types are the same as in Fig. 9 ; note that some models do not deviate significantly from the fiducial WMAP 5+evolving dust model (solid black line) and are therefore not visible. Note that results here have been rescaled to a Euclidean geometry. In the UV bands, data are from GALEX ( Xu et al. 2005 , green squares), STIS on HST ( Gardner, Brown & Ferguson 2000 , purple asterisks) and the balloon-borne FOCA experiment ( Iglesias-Páramo et al. 2004 ; Milliard et al. 1992 , red stars and open pentagons, respectively). The FOCA points have been converted to the GALEX bands using the method described in Xu et al. (2005) . Blue crosses are from HST ACS/SBC observations of multiple fields in GOODS-N and -S ( Voyer et al. 2011 ). In the ACS bands, red, blue and green squares are from the compilation by Dolch & Ferguson (in preparation), which includes data from the Hubble Ultra Deep Field . Additional data in orange from SDSS-DR6 are from Montero-Dorta & Prada (2009) . In the K band, we show data from 6dF (orange crosses, Jones et al. 2006 ), from Keenan et al. (2010 , open red hexagons), from Barro et al. (2009 , blue squares) and McCracken et al. (2010 , green pentagons). All observational data have been converted to AB magnitudes.

Figure 15

Number counts in the GALEX UV bands and the four HST ACS bands. Line types are the same as in Fig. 9 ; note that some models do not deviate significantly from the fiducial WMAP 5+evolving dust model (solid black line) and are therefore not visible. Note that results here have been rescaled to a Euclidean geometry. In the UV bands, data are from GALEX ( Xu et al. 2005 , green squares), STIS on HST ( Gardner, Brown & Ferguson 2000 , purple asterisks) and the balloon-borne FOCA experiment ( Iglesias-Páramo et al. 2004 ; Milliard et al. 1992 , red stars and open pentagons, respectively). The FOCA points have been converted to the GALEX bands using the method described in Xu et al. (2005) . Blue crosses are from HST ACS/SBC observations of multiple fields in GOODS-N and -S ( Voyer et al. 2011 ). In the ACS bands, red, blue and green squares are from the compilation by Dolch & Ferguson (in preparation), which includes data from the Hubble Ultra Deep Field . Additional data in orange from SDSS-DR6 are from Montero-Dorta & Prada (2009) . In the K band, we show data from 6dF (orange crosses, Jones et al. 2006 ), from Keenan et al. (2010 , open red hexagons), from Barro et al. (2009 , blue squares) and McCracken et al. (2010 , green pentagons). All observational data have been converted to AB magnitudes.

Figure 16

Number counts from four Spitzer (IRAC and MIPS) infrared bands, as well as Herschel 250 m and SCUBA 850 m. Line types are the same as in Fig. 9 ; for clarity models similar to the fiducial model are not shown. Results are scaled to a Euclidean geometry. Solid blue circles in the 3.6 IRAC band are from Sanders et al. (2007) ; all other points in the IRAC 3.6 and 8.0 bands are from Fazio et al. (2004) . The MIPS data at 24 m shown here are the S-COSMOS ‘Extragalactic Wide’ points from Sanders et al. (2007) (green hexes), and the Wide and Deep Legacy Survey points from Béthermin et al. (2010) (blue squares). At 70 m data shown are the normal (blue squares) and stacked (cyan squares) measurements from Béthermin et al. (2010) , while red stars are from Frayer et al. (2006) . Herschel data at 250 m are from Clements et al. (2010 , red squares) and Glenn et al. (2010 , blue stars); the latter is from the spline model with FIRAS priors. We show data from the SCUBA SHADES survey ( Coppin et al. 2006 ) at 850 m in the lower-right panel.

Figure 16

Number counts from four Spitzer (IRAC and MIPS) infrared bands, as well as Herschel 250 m and SCUBA 850 m. Line types are the same as in Fig. 9 ; for clarity models similar to the fiducial model are not shown. Results are scaled to a Euclidean geometry. Solid blue circles in the 3.6 IRAC band are from Sanders et al. (2007) ; all other points in the IRAC 3.6 and 8.0 bands are from Fazio et al. (2004) . The MIPS data at 24 m shown here are the S-COSMOS ‘Extragalactic Wide’ points from Sanders et al. (2007) (green hexes), and the Wide and Deep Legacy Survey points from Béthermin et al. (2010) (blue squares). At 70 m data shown are the normal (blue squares) and stacked (cyan squares) measurements from Béthermin et al. (2010) , while red stars are from Frayer et al. (2006) . Herschel data at 250 m are from Clements et al. (2010 , red squares) and Glenn et al. (2010 , blue stars); the latter is from the spline model with FIRAS priors. We show data from the SCUBA SHADES survey ( Coppin et al. 2006 ) at 850 m in the lower-right panel.

Moving towards longer wavelengths, the agreement becomes increasingly poor. Our models underpredict galaxies fainter than Jy at 250 m by a factor of 4–5, and underproduce sub-mm galaxies at 850 m by more than an order of magnitude. It is well known that CDM-based models in general have difficulty reproducing the observed population of sub-mm galaxies (e.g. Devriendt & Guiderdoni 2000 ; Baugh et al. 2005 ). We discuss different possible interpretations of these results in Section 4 , and present more detailed predictions for galaxies in the Herschel PACS and SPIRE wavebands in Niemi et al. (2012) .

### 3.3 The extragalactic background light

One of the major goals of this work is to predict the integrated EBL. The EBL consists of light emitted at all epochs, modified by redshifting and dilution due to the expansion of the universe. Because the EBL consists of light emitted by all types of sources at all cosmological epochs, it forms a unique record of the history of photon production in the universe. A detailed measurement of the EBL flux can potentially inform us about the existence or non-existence of photon sources beyond those that can be resolved in galaxy surveys, and its SED encodes the details of the redshifts and spectral characteristics of these sources. As discussed in our companion paper (GSPD), the EBL also affects the propagation of gamma rays in the GeV and TeV bands. The EBL at a redshift z0 and frequency ν 0 in proper coordinates can be computed in our model as the following integral over emissivities from our predicted galaxy population:

5
where ε(ν, z ) is the galaxy emissivity at redshift z and frequency ν=ν 0 (1 + z )/(1 + z0 ), and d l /d z is the cosmological line element, defined as
6
for a flat ΛCDM universe ( Peebles 1993 ). We assume here that the EBL photons evolve passively after leaving their source galaxies and are not affected by any further interactions except for cosmological redshifting. This is an acceptable approximation for photons at energies below the Rydberg energy of 13.61 eV.

We note that our estimate of the EBL does not include the contribution from light radiated by AGN. However, previous studies have shown that this contribution should be less than ∼10–20 per cent in the mid-IR, and smaller (a few per cent) at other wavelengths (see Domínguez et al. 2011 , and references therein).

Observational estimates of the EBL are obtained either by direct detection or by integration of galaxy counts. Direct detection is complicated by foreground emission from our own galaxy and reflected zodiacal light from our Sun, which are much brighter than the EBL across most of the optical and IR spectrum ( Hauser & Dwek 2001 ). Integrated galaxy counts provide a firm lower limit, but there has been considerable debate about how much light these estimates might be missing because of undetected, low surface brightness galaxies or the faint extended wings of detected galaxies, which can be difficult to disentangle from the background. Fig. 17 shows the predictions of our five model variants along with a compilation of recent observational estimates of the EBL at various wavelengths from both methods. There is a significant gap between the direct detection measurements and the integrated counts, with the former always lying higher than the latter, indicating that there must still be biases or systematic sources of error in one or both methods. Unsurprisingly, since we have already seen that our models reproduce the galaxy counts at most wavelengths, our model predictions all lie close to the estimates from integrated counts. It is noteworthy that although our models, if anything, seem to over -produce faint galaxies, the integrated EBL predictions are nowhere nearly as high as the direct detection estimates. The model with the one-component (Calzetti) attenuation law overproduces the FUV EBL and is low in the mid-IR. This is because in the two-component (modified Charlot–Fall) dust model, young stars are enshrouded in dense birth clouds with higher optical depth. The largest difference between the models with different dust emission templates is seen in the mid-IR, with the DGS99 models predicting a higher EBL in the mid-IR than the R09 models. The two models bracket the error bars on the existing observations. Hopefully new observations will tighten up these constraints. The peak of the EBL at ∼100–250 m is a bit low compared with observational estimates, but is within the errors. In terms of the overall partition of the EBL in the UV–NIR versus NIR–FIR regimes, our models are in good agreement with the observational constraints (see GSPD).

Figure 17

The predicted integrated EBL spectrum, compared with experimental constraints. Line types follow the convention of earlier plots. Data: upwards pointing arrows indicate lower bounds from number counts; other symbols are results from direct detection experiments. Note that some points have been shifted slightly for clarity. Lower limits: the blue–violet triangles are results from Hubble and STIS ( Gardner et al. 2000 ), while the purple open triangles are from GALEX ( Xu et al. 2005 ). The solid green and red triangles are from the Hubble Deep Field ( Madau & Pozzetti 2000 ) and Ultra Deep Field (Dolch & Ferguson, in preparation), respectively, combined with ground based-data, and the solid purple triangle is from a measurement by the Large Binocular Camera ( Grazian et al. 2009 ). In the NIRJ , H and K bands, open violet stars are the limits from Keenan et al. (2010) . Open red triangles are from IRAC on Spitzer ( Fazio et al. 2004 ), and the purple triangle at 15 m is from ISOCAM ( Hopwood et al. 2010 ) on ISO . The lower limits from MIPS at 24, 70 and 160 m on Spitzer are from Béthermin et al. (2010) , Chary et al. (2004) , Frayer et al. (2006) and Dole et al. (2006) (solid blue, solid gold, open gold and open green, respectively). Lower limits from Herschel number counts ( Berta et al. 2010 ) are shown as solid red triangles. In the submillimetre, limits are presented from the BLAST experiment ( Devlin et al. 2009 ). Direct detection: in the optical, orange hexagons are based on data from the Pioneer 10/11 Imaging Photopolarimeter ( Matsuoka et al. 2011 ). The blue star is a determination from Mattila et al. (2011) , and the triangle at 520 nm is an upper limit from the same. In the NIR, the points at 1.25, 2.2 and 3.5 m are based upon DIRBE data with foreground subtraction: Wright (2001 , dark red squares), Cambrésy et al. (2001 , orange crosses), Levenson & Wright (2008 , red diamond), Gorjian, Wright & Chary (2000 , purple open hexes), Wright & Reese (2000 , green square) and Levenson, Wright & Johnson (2007 , red asterisks). In the FIR, direct detection measurements are shown from DIRBE ( Wright 2004 ; Schlegel, Finkbeiner & Davis 1998 , blue stars and solid red circles), and FIRAS ( Fixsen et al. 1998 , purple bars). Blue–violet open squares are from direct photometry with the AKARI satellite ( Matsuura et al. 2010 ).

Figure 17

The predicted integrated EBL spectrum, compared with experimental constraints. Line types follow the convention of earlier plots. Data: upwards pointing arrows indicate lower bounds from number counts; other symbols are results from direct detection experiments. Note that some points have been shifted slightly for clarity. Lower limits: the blue–violet triangles are results from Hubble and STIS ( Gardner et al. 2000 ), while the purple open triangles are from GALEX ( Xu et al. 2005 ). The solid green and red triangles are from the Hubble Deep Field ( Madau & Pozzetti 2000 ) and Ultra Deep Field (Dolch & Ferguson, in preparation), respectively, combined with ground based-data, and the solid purple triangle is from a measurement by the Large Binocular Camera ( Grazian et al. 2009 ). In the NIRJ , H and K bands, open violet stars are the limits from Keenan et al. (2010) . Open red triangles are from IRAC on Spitzer ( Fazio et al. 2004 ), and the purple triangle at 15 m is from ISOCAM ( Hopwood et al. 2010 ) on ISO . The lower limits from MIPS at 24, 70 and 160 m on Spitzer are from Béthermin et al. (2010) , Chary et al. (2004) , Frayer et al. (2006) and Dole et al. (2006) (solid blue, solid gold, open gold and open green, respectively). Lower limits from Herschel number counts ( Berta et al. 2010 ) are shown as solid red triangles. In the submillimetre, limits are presented from the BLAST experiment ( Devlin et al. 2009 ). Direct detection: in the optical, orange hexagons are based on data from the Pioneer 10/11 Imaging Photopolarimeter ( Matsuoka et al. 2011 ). The blue star is a determination from Mattila et al. (2011) , and the triangle at 520 nm is an upper limit from the same. In the NIR, the points at 1.25, 2.2 and 3.5 m are based upon DIRBE data with foreground subtraction: Wright (2001 , dark red squares), Cambrésy et al. (2001 , orange crosses), Levenson & Wright (2008 , red diamond), Gorjian, Wright & Chary (2000 , purple open hexes), Wright & Reese (2000 , green square) and Levenson, Wright & Johnson (2007 , red asterisks). In the FIR, direct detection measurements are shown from DIRBE ( Wright 2004 ; Schlegel, Finkbeiner & Davis 1998 , blue stars and solid red circles), and FIRAS ( Fixsen et al. 1998 , purple bars). Blue–violet open squares are from direct photometry with the AKARI satellite ( Matsuura et al. 2010 ).

## 4 DISCUSSION

### 4.1 Comparison with previous results

Guiderdoni et al. (1998) was one of the earliest studies using ‘forward evolution’, cosmologically motivated CDM models combined with modelling of both dust extinction and emission. This work was extended in Devriendt & Guiderdoni (2000) and subsequent papers by the GALICS collaboration ( Hatton et al. 2003 ; Blaizot et al. 2004 ). They coupled a SAM with analytic recipes for dust attenuation and theoretical templates for dust emission, very much in the same spirit as our work here. The Devriendt & Guiderdoni (2000) models overpredicted the UV and optical counts, possibly because of their use of a single dust attenuation relation rather than a two-component model like our modified Charlot–Fall model. They also found that their standard model failed to reproduce enough sub-mm galaxies at 850 m, but were able to fit the sub-mm counts by introducing by hand a population of heavily dust extinguished starburst galaxies at high redshift.

The Durham group has coupled their SAM of galaxy formation ( Cole et al. 1994 ; Cole et al. 2000 ) with the dust models and RT code GRASIL, developed by Silva et al. (1998) . An important feature of the GRASIL approach is that the dust emission SED is computed self-consistently based on the assumed properties of the dust and the predicted radiation field. This work is presented in an extensive series of papers ( Granato et al. 2000 ; Baugh et al. 2005 ; Swinbank et al. 2008 ; Lacey et al. 2008 , 2010 ). Granato et al. (2000) showed that the fiducial model of Cole et al. (2000) , when coupled with the GRASIL machinery, could reproduce local galaxy LFs from the UV to the FIR (12-100 m). However, Baugh et al. (2005) showed that this model failed to reproduce sufficient numbers of bright sub-mm galaxies by an order of magnitude or more. They therefore made several modifications to their model. They modified the star formation recipes in their model in two ways: (1) they adopted a quiescent star formation recipe with a constant star formation time-scale, rather than a standard Kennicutt–Schmidt relation, leading to larger gas reservoirs in high-redshift galaxies; (2) they adopted an efficient starburst mode in minor as well as major mergers. As shown by Somerville et al. (2001) , this leads to a much larger population of luminous starburst galaxies at high redshift. In addition, they adopted a top-heavy stellar IMF in bursts, in which the slope of the IMF d N /dln mmx is x = 0 over the whole range 0.15 < m < 125 M . In addition to producing much more UV light per unit mass of stars formed, the top-heavy IMF also produces much larger amounts of metals and dust. They found that all of these changes combined were needed to reproduce the sub-mm counts. The same model simultaneously reproduces the rest-UV and optical LFs at high redshift.

Lacey et al. (2008) argue that the top-heavy IMF is needed in order to reproduce the evolution of the mid-IR LF as measured by Spitzer , and the counts in Spitzer bands as well. However, Swinbank et al. (2008) showed that the K band and IRAC 3–8 m luminosities of both SMGs and LBGs in the Lacey–Baugh models are a factor of 10 too low (because of the dearth of low-mass stars). Moreover, the predictions of the models do not seem to be in good agreement with early results from BLAST and Herschel at 350–500 m ( Clements et al. 2010 ; Lacey et al. 2010 ).

Comparison between our results and those of Lacey et al. is complicated because not only the approach for treating dust is different, but also many of the other model ingredients are significantly different. Perhaps most significantly, the published Lacey–Baugh models do not include ‘radio mode’ feedback from AGN, which is now commonly adopted in SAMs in order to prevent the formation of grossly overmassive galaxies (e.g. S08; Bower et al. 2006 ; Croton et al. 2006 ). In an attempt to solve the overcooling problem via other means, Cole et al. (2000) adopted modified hot gas profiles that suppressed early cooling in massive haloes. This also suppresses the early formation of massive galaxies ( Bower et al. 2006 ). Interestingly, Fontanot et al. (2007) use basically the same approach to modelling dust (GRASIL) within a different SAM, and find that they are able to reproduce the K band and sub-mm counts with a Salpeter IMF. They ascribe this success to an improved cooling model. However, their models overproduce massive galaxies at low redshift ( z < 1). Clearly, more work is needed in order to develop a model that reproduces all the available observations.

One should also keep in mind that the GRASIL dust+RT models used in the work described above assume a simplified, regular geometry (spheroid+disc). Particularly when considering populations that are likely to correspond to merger-driven starbursts, such as the SMGs, it is probably important to capture the complex physics and geometry of these systems. Some significant progress has been made recently in this regard, using hydrodynamic simulations of isolated and merging galaxies in combination with the polychromatic Monte Carlo Radiative Transfer code sunrise ( Jonsson 2006 ; Jonsson et al. 2006 , 2010 ; Jonsson & Primack 2010 ). Narayanan et al. (2010a) studied the SEDs of merger simulations combined with sunrise RT calculations and identified objects with properties consistent with both bright and fainter SMGs. Narayanan et al. (2010b) used similar simulations to propose a physical model for the z ∼ 2 DOGs (optically faint galaxies identified at 24 m), which partially overlap with the SMG population. The SEDs from these simulations can be combined with predictions of merger rates from semi-empirical models like those of Hopkins et al. (2010) , or with SAMs, to compute statistics of the population, such as counts (Hayward et al., in preparation).

There is an extensive literature on ‘backwards evolution’ models for galaxy counts and the EBL, which we review in our companion paper GSPD, but do not discuss here. Recently, empirical and ‘semi-empirical’ approaches have been adopted by several authors to predict the EBL. Younger & Hopkins (2011) used the observed stellar mass function at different redshifts, in combination with a semi-empirical model of galaxy evolution and template SEDs from Chary & Elbaz (2001) to predict the mid- to FIR EBL. Domínguez et al. (2011 , D11) made use of empirical template SEDs and observed fractions of 25 different galaxy types to predict the EBL and its evolution. An explicit comparison with the luminosity density evolution estimated by their approach is shown in Fig. 14 . The D11 estimates are anchored to the observed K -band LFs, which our SAMs reproduce fairly well, so the predictions are very similar at NIR wavelengths. At shorter wavelengths (optical and UV), the D11 approach predicts lower luminosity densities at high redshift than our SAMs. These differences are due to the use by D11 of SWIRE templates ( Polletta et al. 2007 ) from the UV to IR, while in our approach we model the star formation history and dust attenuation of each galaxy. In the FIR, D11 estimate a higher and sharper peak at z ∼ 1–2 (again because of the use of different SED templates), in better agreement with observations, and a steeper decline at higher redshift z ≳ 2. Note that the observed K -band LFs that ground their empirical approach are available only up to z ∼ 4, and the results shown at higher redshifts are extrapolations.

### 4.2 Summary and conclusions

We have presented predictions for the luminosity and flux distributions of galaxies from the FUV to the FIR and over the bulk of cosmic history ( z = 0–6). Our predictions are based on SAMs of galaxy formation, set within the hierarchical CDM paradigm of structure formation, and including modelling of gas cooling, star formation, stellar feedback, chemical enrichment and AGN feedback. In addition, crucial to the present study is modelling of the attenuation and re-emission of starlight by dust in the ISM of galaxies. We use a simple but physically motivated analytic approach to estimate the dust attenuation as a function of wavelength. In our fiducial models, based on the approach proposed by Charlot & Fall (2000) , young stars are enshrouded in dense ‘birth clouds’, while older stellar populations are embedded within a more diffuse ‘cirrus’ component. Stars emerge from the dense birth clouds as they age. This two-component dust model results in an effectively age-dependent attenuation relation, such that younger stars are more extinguished. We find that the two-component model gives much better agreement with the UV-optical colours of galaxies than the widely used approach of a fixed attenuation curve.

We then assume that all light absorbed by dust is re-radiated in the IR, and use a set of ‘template’ SEDs to estimate IR luminosities. Our current assumption is that the total IR luminosity of a galaxy determines which template SED to use, based on the empirical correlation between bolometric or total IR luminosity and dust temperature in local LIRGS and ULIRGS ( Sanders & Mirabel 1996 ). However, this is certainly too simplistic, and a goal of our future research is to try to understand and characterize how the physical parameters of galaxies impact the shape of their FIR SEDs. One avenue towards this goal is to use detailed dust and RT models implemented within hydrodynamic simulations ( Jonsson 2006 ; Jonsson et al. 2006 , 2010 ; Narayanan et al. 2010a , b ). Another approach is to use the much richer set of mid- and FIR data that are becoming available, spanning a broader range of cosmic epoch and galaxy type, to develop a more complete set of template SEDs (e.g. Chary & Pope 2010 ).

We found that in order to fit the LFs and counts in the UV and (to a lesser extent) optical, it was necessary to adopt a dust optical depth normalization that varied with redshift. This has the net effect that galaxies of a given bolometric luminosity are less extinguished at high redshift, and could be interpreted as an evolving dust-to-metal ratio or as a geometric effect (e.g. perhaps the distribution of gas and dust relative to the stars is different in high-redshift galaxies). This result has also been found by other groups working with SAMs ( Guo & White 2009 ; Lo Faro et al. 2009 ) and is supported by direct observational evidence ( Reddy et al. 2010 ). However, our current approach is completely ad hoc (we simply tuned the dust normalization parameters to match the observed FUV LFs from z = 0 to 5), and it would be desirable to develop better observational constraints as well as a deeper physical understanding of this effect. With the evolving dust model, we find very good agreement with observed FUV LFs from z ∼ 0 to 5 and B -band LFs from z ∼ 0 to 3, and slightly overpredict galaxies in the rest NIR ( K band) at high redshift ( z ∼ 2–3). In all UV-NIR bands our models predict an excess of low-luminosity galaxies, which confirms the excess found by Fontanot et al. (2009a) and others based on stellar mass function comparisons.

It would have been unsurprising if our simple approach for computing IR luminosities disagreed drastically with more detailed RT calculations or with observations. However, Fontanot et al. (2009b) and Fontanot & Somerville (2010) showed that implementing a recipe similar to the one we adopt here within the MORGANA SAM produced very similar results to the full RT calculations using the GRASIL code of Silva et al. (1998) , at least for global quantities such as LFs and counts. Moreover, we have shown that our models reproduce galaxy counts from the UV to the mid-IR (∼70 m) remarkably well. A growing discrepancy starts to emerge at longer wavelengths: our models underproduce intermediate luminosity (∼30–100 mJy) sources in the SPIRE 250 m band by a factor of 2–5, and S ≳ 5 mJy sources at 850 m by an order of magnitude or more. This problem is far from being unique to our models, as we have discussed above. Clements et al. (2010) show that none of the models that they compare with their SPIRE 250, 350 and 500 m count data agrees with their results very well. Most of these models are empirical ‘backwards evolution’ models, but they also compare with the SAM of Lacey et al. (2008) , which predicts a significant excess of luminous galaxies in all three SPIRE bands.

When comparing with estimates of LFs at z ∼ 0.5–2, we find deficits of luminous galaxies at high redshifts in the rest 8 m and 24 m bands. In the mid-IR (8 and 24 m), it is possible that there could be significant contamination from obscured AGN, particularly at z ∼ 2 (e.g. Daddi et al. 2007 ). The agreement with the estimated total IR LF to z ∼ 2 is not terrible (within the observational errors). It is important to remember that these observational estimates rely on k -correcting from an observed wavelength to the rest frame in a messy part of the SED (particularly in the case of 8 and 24 m), or on estimating a total IR luminosity from a single or a limited number of observed wavelengths. These conversions are themselves highly uncertain and in general rely on SED templates. Therefore, it is encouraging that the agreement between our models and the more directly observed quantity, the galaxy counts, is in general superior to the agreement with the derived quantities (rest-frame or total LFs). In our companion paper (GSPD), we also show the redshift dependence of the build-up of the EBL at 24 m, 70 m and 160 m compared with available observations and find fairly good agreement.

By integrating over all galaxies, accounting for the redshifting and dilution of light, we estimate the integrated EBL predicted by our models. As expected, our EBL predictions lie close to the lower limits from integrated galaxy counts. The largest uncertainties (model-to-model differences) in our predictions are in the mid-IR, due to the limitations of available IR templates. These will improve as more data from multi-wavelength observations are synthesized. In particular, important constraints on this part of the EBL, and correspondingly on the star formation history and dust SEDs of galaxies, may be obtained from observations of GeV and TeV gamma rays. High energy gamma rays are attenuated via electron–positron pair production against the EBL. In principle, the cosmological history of the EBL could be reconstructed by comparing observations of high-energy sources at different redshifts to their known intrinsic spectra. In a companion paper (GSPD), we provide a detailed analysis of the implications of our predictions for current and future gamma-ray observations.

We dedicate this paper to the memory of Donald Lee MacMinn. We warmly thank James Bullock, Julien Devriendt, Bruno Guiderdoni, David Elbaz and Fabio Fontanot for useful discussions. We thank Matthieu Bethermin and Chris Kochanek for pointing out additional observational data to include in our compilations, and we thank the anonymous referee for comments that improved the paper. We thank Elysse Voyer and Timothy Dolch for providing their observational data ahead of publication. RSS and RCG thank UCSC for hospitality on numerous occasions during the gestation of this paper. RCG acknowledges support from a Fermi Guest Investigator grant and a research fellowship from the SISSA Astrophysics Sector. JRP acknowledges support from NASA ATP grant AST-1010033.

## Footnotes

1
The convergence is not as good as in the NIR, where there are still significant uncertainties regarding the importance of contributions from thermally pulsating asymptotic giant branch (TP-AGB) stars ( Maraston 2005 ).

## REFERENCES

Arnouts
S.
et al. ,
2005
,
ApJ
,
619
,
L43
Arnouts
S.
et al. ,
2007
,
A&A
,
476
,
137
Babbedge
T. S. R.
et al. ,
2006
,
MNRAS
,
370
,
1159
Baldry
I. K.
Glazebrook
K.
Driver
S. P.
,
2008
,
MNRAS
,
388
,
945
Baldry
I. K.
et al. ,
2011
, preprint (arXiv e-prints)
Barro
G.
et al. ,
2009
,
A&A
,
494
,
63
Baugh
C. M.
Lacey
C. G.
Frenk
C. S.
Granato
G. L.
Silva
L.
Bressan
A.
Benson
A. J.
Cole
S.
,
2005
,
MNRAS
,
356
,
1191
Bell
E. F.
,
2003
,
ApJ
,
586
,
794
Bell
E. F.
McIntosh
D. H.
Katz
N.
Weinberg
M. D.
,
2003
,
ApJS
,
149
,
289
Berta
S.
et al. ,
2010
,
A&A
,
518
,
L30
Béthermin
M.
Dole
H.
Beelen
A.
Aussel
H.
,
2010
,
A&A
,
512
,
A78
Birnboim
Y.
Dekel
A.
,
2003
,
MNRAS
,
345
,
349
Blaizot
J.
Guiderdoni
B.
Devriendt
J. E. G.
Bouchet
F. R.
Hatton
S. J.
Stoehr
F.
,
2004
,
MNRAS
,
352
,
571
Blanton
M. R.
et al. ,
2003
,
ApJ
,
592
,
819
Blumenthal
G.
Faber
S. M.
Flores
R.
Primack
J. R.
,
1986
,
ApJ
,
301
,
27
Bondi
H.
,
1952
,
MNRAS
,
112
,
195
Bouwens
R. J.
Illingworth
G. D.
Franx
M.
Ford
H.
,
2007
,
ApJ
,
670
,
928
Bower
R. G.
Benson
A. J.
Malbon
R.
Helly
J. C.
Frenk
C. S.
Baugh
C. M.
Cole
S.
Lacey
C. G.
,
2006
,
MNRAS
,
370
,
645
Boylan-Kolchin
M.
Ma
C.-P.
Quataert
E.
,
2008
,
MNRAS
,
383
,
93
Bruzual
G.
Charlot
S.
,
2003
,
MNRAS
,
344
,
1000
Buat
V.
et al. ,
2007
,
ApJS
,
173
,
404
Calzetti
D.
Armus
L.
Bohlin
R. C.
Kinney
A. L.
Koornneef
J.
Storchi-Bergmann
T.
,
2000
,
ApJ
,
533
,
682
Cambrésy
L.
Reach
W. T.
Beichman
C. A.
Jarrett
T. H.
,
2001
,
ApJ
,
555
,
563
Caputi
K. I.
McLure
R. J.
Dunlop
J. S.
Cirasuolo
M.
Schael
A. M.
,
2006
,
MNRAS
,
366
,
609
Caputi
K. I.
et al. ,
2007
,
ApJ
,
660
,
97
Cattaneo
A.
et al. ,
2007
,
MNRAS
,
377
,
63
Chabrier
G.
,
2003
,
PASP
,
115
,
763
Chapman
S. C.
Blain
A. W.
Smail
I.
Ivison
R. J.
,
2005
,
ApJ
,
622
,
772
Charlot
S.
Fall
S. M.
,
2000
,
ApJ
,
539
,
718
Chary
R.
Elbaz
D.
,
2001
,
ApJ
,
556
,
562
Chary
R.
Pope
A.
,
2010
, preprint (arXiv e-prints)
Chary
R.
et al. ,
2004
,
ApJS
,
154
,
80
Cirasuolo
M.
McLure
R. J.
Dunlop
J. S.
Almaini
O.
Foucaud
S.
Simpson
C.
,
2010
,
MNRAS
,
401
,
1166
Clements
D. L.
et al. ,
2010
,
A&A
,
518
,
L8
Cole
S.
Aragón-Salamanca
A.
Frenk
C. S.
Navarro
J. F.
Zepf
S. E.
,
1994
,
MNRAS
,
271
,
781
Cole
S.
Lacey
C. G.
Baugh
C. M.
Frenk
C. S.
,
2000
,
MNRAS
,
319
,
168
Cole
S.
et al. ,
2001
,
MNRAS
,
326
,
255
Coppin
K.
et al. ,
2006
,
MNRAS
,
372
,
1621
Cortese
L.
et al. ,
2006
,
ApJ
,
637
,
242
Croton
D. J.
et al. ,
2006
,
MNRAS
,
365
,
11
E.
et al. ,
2007
,
ApJ
,
670
,
173
Dahlen
T.
Mobasher
B.
Somerville
R. S.
Moustakas
L. A.
Dickinson
M.
Ferguson
H. C.
Giavalisco
M.
,
2005
,
ApJ
,
631
,
126
Dahlen
T.
Mobasher
B.
Dickinson
M.
Ferguson
H. C.
Giavalisco
M.
Kretchmer
C.
Ravindranath
S.
,
2007
,
ApJ
,
654
,
172
Dai
X.
et al. ,
2009
,
ApJ
,
697
,
506
Dale
D. A.
Helou
G.
,
2002
,
ApJ
,
576
,
159
De Lucia
G.
Blaizot
J.
,
2007
,
MNRAS
,
375
,
2
Dekel
A.
Birnboim
Y.
,
2006
,
MNRAS
,
368
,
2
Desert
F.
Boulanger
F.
Puget
J. L.
,
1990
,
A&A
,
237
,
215
Devlin
M. J.
et al. ,
2009
,
Nat
,
458
,
737
Devriendt
J. E. G.
Guiderdoni
B.
,
2000
,
A&A
,
363
,
851
Devriendt
J. E. G.
Guiderdoni
B.
R.
,
1999
,
A&A
,
350
,
381
(DGS99)
Di Matteo
T.
Springel
V.
Hernquist
L.
,
2005
,
Nat
,
433
,
604
Dole
H.
et al. ,
2006
,
A&A
,
451
,
417
Domínguez
A.
et al. ,
2011
,
MNRAS
,
410
,
2556
Dye
S.
et al. ,
2010
,
A&A
,
518
,
L10
Elbaz
D.
et al. ,
1999
,
A&A
,
351
,
L37
Elbaz
D.
Cesarsky
C. J.
Chanial
P.
Aussel
H.
Franceschini
A.
D.
Chary
R. R.
,
2002
,
A&A
,
384
,
848
Faber
S. M.
et al. ,
2007
,
ApJ
,
665
,
265
Fazio
G. G.
et al. ,
2004
,
ApJS
,
154
,
39
Fixsen
D. J.
Dwek
E.
Mather
J. C.
Bennett
C. L.
Shafer
R. A.
,
1998
,
ApJ
,
508
,
123
Flores
R.
Primack
J. R.
Blumenthal
G.
Faber
S. M.
,
1993
,
ApJ
,
412
,
443
Fontanot
F.
Somerville
R. S.
,
2010
, preprint (arXiv e-prints)
Fontanot
F.
Monaco
P.
Silva
L.
Grazian
A.
,
2007
,
MNRAS
,
382
,
903
Fontanot
F.
De Lucia
G.
Monaco
P.
Somerville
R. S.
Santini
P.
,
2009a
,
MNRAS
,
397
,
1776
Fontanot
F.
Somerville
R. S.
Silva
L.
Monaco
P.
Skibba
R.
,
2009b
,
MNRAS
,
392
,
553
Frayer
D. T.
et al. ,
2006
,
ApJ
,
647
,
L9
Gabasch
A.
et al. ,
2006
,
A&A
,
448
,
101
Gardner
J. P.
Brown
T. M.
Ferguson
H. C.
,
2000
,
ApJ
,
542
,
L79
Giallongo
E.
Salimbeni
S.
Menci
N.
Zamorani
G.
Fontana
A.
Dickinson
M.
Cristiani
S.
Pozzetti
L.
,
2005
,
ApJ
,
622
,
116
Gilmore
R.
Somerville
R.
Primack
J.
Domínguez
A.
,
2011
,
MNRAS
, in press (doi:10.1111/j.1365-2966.2012.20841.x)
Glenn
J.
et al. ,
2010
,
MNRAS
,
409
,
109
Gnedin
N. Y.
,
2000
,
ApJ
,
542
,
535
Gorjian
V.
Wright
E. L.
Chary
R. R.
,
2000
,
ApJ
,
536
,
550
Granato
G. L.
Lacey
C. G.
Silva
L.
Bressan
A.
Baugh
C. M.
Cole
S.
Frenk
C. S.
,
2000
,
ApJ
,
542
,
710
Grazian
A.
et al. ,
2009
,
A&A
,
505
,
1041
Grogin
N. A.
et al. ,
2011
, preprint (arXiv e-prints)
Guiderdoni
B.
Rocca-Volmerange
B.
,
1987
,
A&A
,
186
,
1
Guiderdoni
B.
Hivon
E.
Bouchet
F. R.
Maffei
B.
,
1998
,
MNRAS
,
295
,
877
Guo
Q.
White
S. D. M.
,
2009
,
MNRAS
,
396
,
39
Guo
Q.
et al. ,
2010
, preprint (arXiv e-prints)
Hathi
N. P.
et al. ,
2010
,
ApJ
,
720
,
1708
Hatton
S.
Devriendt
J. E. G.
Ninin
S.
Bouchet
F. R.
Guiderdoni
B.
Vibert
D.
,
2003
,
MNRAS
,
343
,
75
Hauser
M. G.
Dwek
E.
,
2001
,
ARA&A
,
39
,
249
Heckman
T. M.
Robert
C.
Leitherer
C.
Garnett
D. R.
van der Rydt
F.
,
1998
,
ApJ
,
503
,
646
Henriques
B.
Maraston
C.
Monaco
P.
Fontanot
F.
Menci
N.
De Lucia
G.
Tonini
C.
,
2010
, preprint (arXiv e-prints)
Hopkins
A. M.
Beacom
J. F.
,
2006
,
ApJ
,
651
,
142
Hopkins
P. F.
Cox
T. J.
Younger
J. D.
Hernquist
L.
,
2009a
,
ApJ
,
691
,
1168
Hopkins
P. F.
et al. ,
2009b
,
MNRAS
,
397
,
802
Hopkins
P. F.
et al. ,
2010
,
ApJ
,
715
,
202
Hopwood
R.
et al. ,
2010
,
ApJ
,
716
,
L45
Hughes
D. H.
et al. ,
1998
,
Nat
,
394
,
241
Iglesias-Páramo
J.
Buat
V.
Donas
J.
Boselli
A.
Milliard
B.
,
2004
,
A&A
,
419
,
109
Jones
D. H.
Peterson
B. A.
Colless
M.
Saunders
W.
,
2006
,
MNRAS
,
369
,
25
Jonsson
P.
,
2006
,
MNRAS
,
372
,
2
Jonsson
P.
Primack
J. R.
,
2010
,
New Astron.
,
15
,
509
Jonsson
P.
Cox
T. J.
Primack
J. R.
Somerville
R. S.
,
2006
,
ApJ
,
637
,
255
Jonsson
P.
Groves
B. A.
Cox
T. J.
,
2010
,
MNRAS
,
403
,
17
Kang
X.
Jing
Y. P.
Silk
J.
,
2006
,
ApJ
,
648
,
820
Kannappan
S. J.
,
2004
,
ApJ
,
611
,
L89
Kauffmann
G.
White
S. D. M.
Guiderdoni
B.
,
1993
,
MNRAS
,
264
,
201
Kauffmann
G.
Colberg
J. M.
Diaferio
A.
White
S. D. M.
,
1998
,
MNRAS
,
303
,
188
Kaviani
A.
Haehnelt
M. G.
Kauffmann
G.
,
2003
,
MNRAS
,
340
,
739
Keenan
R. C.
Barger
A. J.
Cowie
L. L.
Wang
W.
,
2010
,
ApJ
,
723
,
40
Kennicutt
R. C.
,
1998
,
ApJ
,
498
,
181
Kereš
D.
Katz
N.
Weinberg
D. H.
Davé
R.
,
2005
,
MNRAS
,
363
,
2
Kochanek
C. S.
et al. ,
2001
,
ApJ
,
560
,
566
Koekemoer
A. M.
et al. ,
2011
, preprint (arXiv e-prints)
Komatsu
E.
et al. ,
2009
,
ApJS
,
180
,
330
Komatsu
E.
et al. ,
2010
, preprint (arXiv e-prints)
Kravtsov
A. V.
Gnedin
O. Y.
Klypin
A. A.
,
2004
,
ApJ
,
609
,
482
Lacey
C.
Guiderdoni
B.
Rocca-Volmerange
B.
Silk
J.
,
1993
,
ApJ
,
402
,
15
Lacey
C. G.
Baugh
C. M.
Frenk
C. S.
Silva
L.
Granato
G. L.
Bressan
A.
,
2008
,
MNRAS
,
385
,
1155
Lacey
C. G.
Baugh
C. M.
Frenk
C. S.
Benson
A. J.
Orsi
A.
Silva
L.
Granato
G. L.
Bressan
A.
,
2010
,
MNRAS
,
405
,
2
Lagache
G.
et al. ,
2004
,
ApJS
,
154
,
112
Le Floc’h
E.
et al. ,
2005
,
ApJ
,
632
,
169
Lee
S.
Idzi
R.
Ferguson
H. C.
Somerville
R. S.
Wiklind
T.
Giavalisco
M.
,
2009
,
ApJS
,
184
,
100
Leroy
A. K.
Walter
F.
Brinks
E.
Bigiel
F.
de Blok
W. J. G.
B.
Thornley
M. D.
,
2008
,
AJ
,
136
,
2782
Levenson
L. R.
Wright
E. L.
,
2008
,
ApJ
,
683
,
585
Levenson
L. R.
Wright
E. L.
Johnson
B. D.
,
2007
,
ApJ
,
666
,
34
Li
C.
White
S. D. M.
,
2009
,
MNRAS
,
398
,
2177
Lo Faro
B.
Monaco
P.
Vanzella
E.
Fontanot
F.
Silva
L.
Cristiani
S.
,
2009
,
MNRAS
,
399
,
827
McCracken
H. J.
et al. ,
2010
,
ApJ
,
708
,
202
Macciò
A. V.
Kang
X.
Fontanot
F.
Somerville
R. S.
Koposov
S.
Monaco
P.
,
2010
,
MNRAS
,
402
,
1995
P.
Pozzetti
L.
,
2000
,
MNRAS
,
312
,
L9
Magnelli
B.
Elbaz
D.
Chary
R. R.
Dickinson
M.
Le Borgne
D.
Frayer
D. T.
Willmer
C. N. A.
,
2011
, preprint (arXiv e-prints)
Maraston
C.
,
2005
,
MNRAS
,
362
,
799
Marchesini
D.
et al. ,
2007
,
ApJ
,
656
,
42
Marchesini
D.
van Dokkum
P. G.
Förster Schreiber
N. M.
Franx
M.
Labbé
I.
Wuyts
S.
,
2009
,
ApJ
,
701
,
1765
Matsuoka
Y.
Ienaka
N.
Kawara
K.
Oyabu
S.
,
2011
,
ApJ
,
736
,
119
Matsuura
S.
et al. ,
2010
, preprint (arXiv e-prints)
Mattila
K.
Lehtinen
K.
Vaisanen
P.
von Appen-Schnur
G.
Leinert
C.
,
2011
, preprint (arXiv e-prints)
Menci
N.
Fontana
A.
Giallongo
E.
Grazian
A.
Salimbeni
S.
,
2006
,
ApJ
,
647
,
753
Milliard
B.
Donas
J.
Laget
M.
Armand
C.
Vuillemin
A.
,
1992
,
A&A
,
257
,
24
Mo
H. J.
Mao
S.
White
S. D. M.
,
1998
,
MNRAS
,
295
,
319
Monaco
P.
Fontanot
F.
Taffoni
G.
,
2007
,
MNRAS
,
375
,
1189
Montero-Dorta
A. D.
F.
,
2009
,
MNRAS
,
399
,
1106
Narayanan
D.
Hayward
C. C.
Cox
T. J.
Hernquist
L.
Jonsson
P.
Younger
J. D.
Groves
B.
,
2010a
,
MNRAS
,
401
,
1613
Narayanan
D.
et al. ,
2010b
,
MNRAS
,
407
,
1701
Niemi
S.-M.
Somerville
R. S.
Ferguson
H. C.
Huang
K.-H.
Lotz
J.
Koekemoer
A.
,
2012
,
MNRAS
, in press (doi:10.1111/j.1365-2966.2012.20425.x)
Panter
B.
Jimenez
R.
Heavens
A. F.
Charlot
S.
,
2007
,
MNRAS
,
378
,
1550
Peebles
P. J. E.
,
1993
,
Principles of Physical Cosmology
.
Princeton Univ. Press
, Princeton, NJ
Polletta
M.
et al. ,
2007
,
ApJ
,
663
,
81
Pozzetti
L.
et al. ,
2003
,
A&A
,
402
,
837
Reddy
N. A.
Steidel
C. C.
D.
Yan
L.
Pettini
M.
Shapley
A. E.
Erb
D. K.
K. L.
,
2006
,
ApJ
,
644
,
792
Reddy
N. A.
Steidel
C. C.
Pettini
M.
K. L.
Shapley
A. E.
Erb
D. K.
Dickinson
M.
,
2008
,
ApJS
,
175
,
48
Reddy
N. A.
Erb
D. K.
Pettini
M.
Steidel
C. C.
Shapley
A. E.
,
2010
,
ApJ
,
712
,
1070
Rieke
G. H.
Alonso-Herrero
A.
Weiner
B. J.
Pérez-González
P. G.
Blaylock
M.
Donley
J. L.
Marcillac
D.
,
2009
,
ApJ
,
692
,
556
(R09)
Robertson
B.
Cox
T. J.
Hernquist
L.
Franx
M.
Hopkins
P. F.
Martini
P.
Springel
V.
,
2006a
,
ApJ
,
641
,
21
Robertson
B.
Hernquist
L.
Cox
T. J.
Di Matteo
T.
Hopkins
P. F.
Martini
P.
Springel
V.
,
2006b
,
ApJ
,
641
,
90
Rodighiero
G.
et al. ,
2010
,
A&A
,
515
,
A8
Rujopakarn
W.
et al. ,
2010
,
ApJ
,
718
,
1171
Salimbeni
S.
et al. ,
2008
,
A&A
,
477
,
763
Sanders
D. B.
Mirabel
I. F.
,
1996
,
ARA&A
,
34
,
749
Sanders
D. B.
et al. ,
2007
,
ApJS
,
172
,
86
Schiminovich
D.
et al. ,
2005
,
ApJ
,
619
,
L47
Schlegel
D. J.
Finkbeiner
D. P.
Davis
M.
,
1998
,
ApJ
,
500
,
525
Silva
L.
Granato
G. L.
Bressan
A.
Danese
L.
,
1998
,
ApJ
,
509
,
103
Smail
I.
Ivison
R. J.
Blain
A. W.
,
1997
,
ApJ
,
490
,
L5
Soifer
B. T.
Neugebauer
G.
,
1991
,
AJ
,
101
,
354
Somerville
R. S.
Kolatt
T. S.
,
1999
,
MNRAS
,
305
,
1
Somerville
R. S.
Primack
J. R.
,
1999
,
MNRAS
,
310
,
1087
Somerville
R. S.
Primack
J. R.
Faber
S. M.
,
2001
,
MNRAS
,
320
,
504
Somerville
R. S.
Hopkins
P. F.
Cox
T. J.
Robertson
B. E.
Hernquist
L.
,
2008a
,
MNRAS
,
391
,
481
(S08)
Somerville
R. S.
et al. ,
2008b
,
ApJ
,
672
,
776
Sutherland
R.
Dopita
M. A.
,
1993
,
ApJS
,
88
,
253
Swinbank
A. M.
et al. ,
2008
,
MNRAS
,
391
,
420
Takeuchi
T. T.
Ishii
T. T.
Hirashita
H.
Yoshikawa
K.
Matsuhara
H.
Kawara
K.
Okuda
H.
,
2001
,
PASJ
,
53
,
37
Voyer
E.
Gardner
J.
Teplitz
H.
Siana
B.
de Mello
D.
,
2011
,
ApJ
, submitted
Wolf
C.
Meisenheimer
K.
Rix
H.-W.
Borch
A.
Dye
S.
Kleinheinrich
M.
,
2003
,
A&A
,
401
,
73
Wright
E. L.
,
2001
,
ApJ
,
553
,
538
Wright
E. L.
,
2004
,
New Astron. Rev.
,
48
,
465
Wright
E. L.
Reese
E. D.
,
2000
,
ApJ
,
545
,
43
Wyder
T. K.
et al. ,
2005
,
ApJ
,
619
,
L15
Xu
C. K.
et al. ,
2005
,
ApJ
,
619
,
L11
Xu
C. K.
et al. ,
2006
,
ApJ
,
646
,
834
Yoshida
N.
Stoehr
F.
Springel
V.
White
S. D. M.
,
2002
,
MNRAS
,
335
,
762
Younger
J. D.
Hopkins
P. F.
,
2011
,
MNRAS
,
410
,
2180

## Author notes

Visiting researcher at the Santa Cruz Institute for Particle Physics (SCIPP), University of California, Santa Cruz, CA 95064, USA.