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.
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.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 (ρ gas ∝ r−2 ), the cooling rate is given bySutherland & 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, tdyn ∝ rvir / Vvir , where Vvir is the virial velocity of the halo. 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
Additionally, stars younger than tBC are enshrouded in a cloud of dust with optical depth τ BC, V =μ BC τ 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 .
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 ).
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.
|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.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 ⊙ ).
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.
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.
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).
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) .
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.
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.
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.
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.
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: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).
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 m ∝ m − x 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.