A backward evolution model for infrared surveys: the role of AGN- and Color-L_TIR distributions
aa r X i v : . [ a s t r o - ph . C O ] J un Accepted by ApJ
Preprint typeset using L A TEX style emulateapj v. 08/22/09
A BACKWARD EVOLUTION MODEL FOR INFRARED SURVEYS: THE ROLE OF AGN- AND COLOR- L TIR
DISTRIBUTIONS
E. Valiante † , D. Lutz , E. Sturm , R. Genzel , E. Chapin Accepted by ApJ
ABSTRACTEmpirical “backward” galaxy evolution models for infrared bright galaxies are constrained usingmulti-band infrared surveys. We developed a new Monte-Carlo algorithm for this task, implementingluminosity dependent distribution functions for the galaxies’ infrared spectral energy distributions(SEDs) and for the AGN contribution, allowing for evolution of these quantities. The adopted SEDstake into account the contributions of both starbursts and AGN to the infrared emission, for the firsttime in a coherent treatment rather than invoking separate AGN and star-forming populations.In the first part of the paper we consider the quantification of the AGN contribution for localuniverse galaxies, as a function of total infrared luminosity. It is made using a large sample of LIRGsand ULIRGs for which mid-infrared spectra are available in the
Spitzer archive. We find the ratioof AGN 6 µ m luminosity and total infrared luminosity to rise with L . over the infrared luminosityrange 10 to 10 L ⊙ and estimate its spread. Judging from the modest number of distant sourceswith Spitzer spectroscopy, the relation changes at high z .In the second part we present the model. Our best-fit model adopts very strong luminosity evolution, L = L (1 + z ) . , up to z = 2 .
3, and density evolution, ρ = ρ (1 + z ) , up to z = 1, for the populationof infrared galaxies. At higher z , the evolution rates drop as (1 + z ) − and (1 + z ) − . respectively. Toreproduce mid-infrared to submillimeter number counts and redshift distributions, it is necessary tointroduce both an evolution in the AGN contribution and an evolution in the luminosity-temperaturerelation. At a given total infrared luminosity, high redshift infrared galaxies have typically smallerAGN contributions to the rest frame mid-infrared, and colder far-infrared dust temperatures thanlocally. We also suggest an extension of the local infrared galaxy population towards lower dusttemperatures.Our models are in plausible agreement with current photometry-based estimates of the typical AGNcontribution as a function of mid-infrared flux, and well placed to be compared to upcoming Spitzer spectroscopic results. As an example of future applications, we use our best-fitting model to makepredictions for surveys with
Herschel . Subject headings: cosmology: source counts, cosmology: redshift distribution, galaxies: active, galax-ies: starburst, galaxies: evolution, infrared: galaxies INTRODUCTION
The discovery of the cosmic infrared background(CIB) (see Hauser & Dwek 2001 for a review), to-gether with recent deep cosmological surveys in theinfrared (IR) and submillimeter bands, have offerednew perspectives on our understanding of galaxy for-mation and evolution. The surprisingly large amountof energy contained in the CIB showed that it is cru-cial to probe its contributing galaxies to understandwhen and how the bulk of stars formed in the Uni-verse. Thanks to the deep cosmological surveys car-ried out by ISO (Kessler et al. 1996, e.g. Aussel et al.1999; Oliver et al. 2000), SCUBA (Holland et al. 1999,e.g. Hughes et al. 1998), MAMBO (Kreysa et al. 1998,e.g. Bertoldi et al. 2000; Dannerbauer et al. 2004) and
Spitzer (Werner et al. 2004, e.g. Papovich et al. 2004;Frayer et al. 2006a; Dole et al. 2006) it is now pos-sible, to various degrees, to resolve the CIB intodiscrete sources. The source counts are high Max-Planck-Institut f¨ur extraterrestrische Physik, Postfach1312, 85741 Garching, Germany Dept. of Physics and Astronomy, Univ. of British Columbia,6224 Agricultural Road, Vancouver, B.C., V6T 1Z1, Canada † [email protected] , when compared to no-evolution or moderate-evolutionmodels for infrared galaxies (Guiderdoni et al. 1998;Franceschini et al. 1998). The striking results of thesesurveys concerning the evolution of the infrared and sub-millimeter galaxy population, and the constraints fromthe measurements of the CIB, require the developmentof new models to explain the high rate of evolution ofinfrared galaxies.The problem of describing the number distribution ofgalaxies in the Universe is usually tackled via one of twomethods. In the first method, known as forward evolu-tion , the calculation assumes some initial conditions andphysical processes for chemical and photometric evolu-tion of the stellar populations that heat the dust. The so-called semi-analytical approach combines these assump-tions about the chemical/photometric evolution of galax-ies with models for the dissipative and non-dissipativeprocesses affecting galaxy formation within dark matterhalos (Baugh et al. 2006) and has provided a reasonablefit to the source counts in the infrared (Guiderdoni et al.1998; Lacey et al. 2008). Arguably, the forward evolution approach has the advantage of being based on a morefundamental set of assumptions. However the obviousdisadvantage is the large number of free/unknown pa- Valiante et al.rameters assumed in these models. Compared to otherwavelengths, there is one further complication for for-ward modelling of infrared surveys. Even if global prop-erties of the evolving galaxies like star formation ratesor AGN activity are correctly modelled, additional as-sumptions about the dust content and structure of thegalaxy have to be invoked to convert them into a predic-tion of the luminosity and spectral energy distribution ofthe re-emitted infrared radiation.The alternative method, often called the backward evo-lution approach, takes the observed, present day ( z = 0)luminosity function (LF) and evolves it in luminosityand/or density out to higher redshifts assuming someparametrization of the evolution (e.g. Franceschini at al.1988; Pearson & Rowan-Robinson 1996; Xu et al.2001; Franceschini et al. 2001; Rowan-Robinson 2001;Lagache et al. 2003). This method has the advantage ofbeing both direct and relatively simple to implement.The disadvantage in the past was that the informationon which the measurements of evolution were oftenmade with IRAS or ISO data which extended out onlyto relatively low redshifts. Because of their simplicity,backward evolution models have traditionally played astrong role in the planning of new infrared surveys. Abackward evolution model fitting the main constraintsprovided by previous missions can be easily modified topredict the results of new missions and help in the firststeps of interpreting their results. Thanks to the latestdeep and comprehensive observations cited above, thetools to constrain the backward evolution methodologyto significantly higher redshifts are now available.Very recently, several empirical approaches have beenproposed to model the high rate of evolution of infraredgalaxies, in particular to reproduce source counts of themid-IR surveys made with Spitzer (e.g. Lagache et al.2003, 2004; Gruppioni et al. 2005; Rowan-Robinson2009).The “classical” backward evolution model starts fromdifferent populations of galaxies, typically cirrus, star-burst, active galactic nuclei (AGN) and ultraluminousinfrared galaxies (ULIRGs) (e.g. Rowan-Robinson 2001),or a subset of those. Each population is assigned a spec-tral energy distribution (SED) and local luminosity func-tion, and evolved independently, assuming it has a propercosmic evolution rate. These hypotheses are not alwayssatisfied. First of all, the SEDs adopted for starburstgalaxies are usually either represented by a single SED orcome from a set of templates where a unique relation be-tween temperature and luminosity is assumed. This as-pect was already discussed by Chapman et al. (2003) andChapin et al. (2009): analyzing a sample of local infraredgalaxies, they found that the luminosity-temperature( L − T ) relation presents a significant spread. The lat-ter, in particular, showed evidence for some evolutionof the relation with redshift, comparing the local trendwith high-redshift data from SHARC-II 350 µ m obser-vations of SMGs (e.g. Kov´acs et al. 2006; Coppin et al.2008). Second, and perhaps more important, differentpopulations of galaxies are not as distinct as assumed inthe models: very often AGN and starbursts co-exist inthe same object, can be predominant at different timesdepending of the evolutionary stage of the object itselfand can influence each other (e.g. AGN feedback on starformation). Backward evolution models based on single SEDs orthe simple SED family approach described above havebeen quite successful in fitting the number counts fromthe IRAS and ISO missions, the first submillimetercounts, the global CIB level, and for making predictionsfor Spitzer and
Herschel . The number and quality of ob-served constraints is increasing, however, and already thefirst
Spitzer results have led to on-the-spot modifications(c.f. the modified SEDs adopted by Lagache et al. 2004)that may either genuinely represent improved knowledge,or reflect limitations of the simple assumptions made pre-viously. Furthermore, questions gaining increased impor-tance, like the co-existence of AGN and star-formationin infrared galaxies, cannot be addressed by this gener-ation of models, not even in the simple sense of fittingexisting data and extrapolating to new observations. Thegoal of this work is hence to take the next step and de-velop a backward evolution model that considers realisticspreads in far-infrared SEDs, in AGN contributions atdifferent luminosities, and their possible evolution withredshift. To that end, we first have to describe the lo-cal situation that is the starting point of the backwardevolution scenario.Luminous (LIRGs: 10 < L TIR ≡ L − µ m < L ⊙ ) and ultraluminous (ULIRGs: L TIR > L ⊙ )infrared galaxies (Sanders et al. 1988) have been studiedextensively in the local Universe, for instance with the Infrared Astronomical Satellite (IRAS; e.g. Soifer et al.1987; Saunders et al. 1990), the
Infrared Space Obser-vatory (ISO; e.g. Lutz et al. 1998; Genzel & Cesarsky2000; Tran et al. 2001), and more recently, with the
In-frared Spectrograph (IRS; Houck et al. 2004) on boardof
Spitzer (e.g. Weedman et al. 2005; Brandl et al.2006; Armus et al. 2007; Desai et al. 2007; Veilleux et al.2009). These galaxies exhibit a large range of propertiesin the mid-IR, some showing strong PAH emission fea-tures characteristic of powerful (up to ≈ M ⊙ yr − )star formation rates (e.g. Brandl et al. 2006; Smith et al.2007), and exhibiting a large range in 9 . µ m silicate ab-sorption or emission strengths (e.g. Weedman et al. 2005;Desai et al. 2007; Imanishi et al. 2007). Spitzer
IRS isnow enabling the study of mid-IR spectra of infraredgalaxies to much higher redshifts ( z &
2; e.g. Houck et al.2005; Yan et al. 2007; Men´endez-Delmestre et al. 2007;Valiante et al. 2007; Pope et al. 2008). Although lo-cally rare, infrared galaxies are an important popula-tion at high redshifts and account for an increasingfraction of the star-formation activity in the universe(Le Floc’h et al. 2005). By studying their infrared prop-erties, we are just starting to estimate the extent to whichAGN and star-formation contribute to their infrared lu-minosities, and therefore determine a correct census ofstarbursts and AGN at epochs in the Universe when theirluminosity density was at its maximum.The local luminosity function of IR-bright galaxiescan be described considering two different populations:normal galaxies, dominating the “low-luminosity” partof LF, and starburst galaxies, dominant in the “high-luminosity” part. The AGN contribution appears domi-nant only at very high luminosity ( L TIR & × L ⊙ ).Nevertheless, a small fraction of the total infrared lu-minosity can be due to the presence of an active galac-tic nucleus even in star formation dominated cases, and backward evolution model for infrared surveys 3vice versa (Genzel et al. 1998). There are several stud-ies of the AGN content for luminous galaxies at ULIRGand HYLIRG (hyperluminous infrared galaxies, L TIR > L ⊙ ) levels, both using X-ray (Franceschini et al.2003; Teng et al. 2005) and mid-IR (e.g. Genzel et al.1998; Lutz et al. 1998; Tran et al. 2001) emission, andit is generally believed that AGN are typically less im-portant at lower infrared luminosities. Still missingis a comprehensive study including lower luminosities( < L ⊙ ), with the aim of quantifying the distributionof AGN luminosity respect to the infrared luminosity ofthe host and the fraction of infrared luminosity due toaccretion. Such a study is made in the following section(see § § z . In the secondpart ( § THE AGN CONTRIBUTION IN LOCAL ULIRGS ANDLIRGS
In order to be useful for the backward modellingproject outlined above, a local calibration of the dis-tribution of AGN contributions at given infrared lumi-nosity has to span the highest luminosities down to atleast the low end of the LIRG regime (10 L ⊙ ). Giventhe evidence for the major contribution of objects abovethis luminosity threshold to the CIB and to the cos-mic star formation rate at z & > L ⊙ )regime and even its upper end where the AGN con-tribution is highest, while studies at lower luminositiesmostly focussed on individual interesting objects. Whatis needed is a quantification of AGN content for an unbi-ased far-IR selected sample reaching down to 10 L ⊙ . Ofthe two principle routes towards this observational goal,X-ray observations and mid-IR spectroscopy, we makeuse of the second for two reasons. First, there is not yeta full unbiased and deep X-ray dataset for such a sam-ple available. There are XMM and Chandra data of thedepth required for good X-ray spectral analysis for manylocal LIRGs/ULIRGs, but still with a tendency to tar-get known AGN. Second, since the implicit goal of ourmodelling effort is to characterize the AGN effect on theinfrared SED, mid-IR data can provide useful constraintsquite directly even when adopting simple analysis meth-ods, while evidence from other wavelengths would have to go through a conversion to the mid-IR range, or evenin two steps via the AGN bolometric power, a processthat will be subject to the considerable scatter of AGNSEDs.Our basic approach is to determine the distribution ofAGN contributions for local LIRGs and ULIRGs fromthe IRAS Revised Bright Galaxy Sample (Sanders et al.2003) for which archival Spitzer spectroscopic coverageis essentially complete. We adopt a simple quantificationof the AGN contribution from the rest frame 6 µ m con-tinuum which, while not making use of the full details ofthe best S/N spectra, is applicable to the entire sampleand allows for straightforward use in the later modelling.In mid-IR spectra of AGN, the 6 µ m flux after subtrac-tion of the starbursts contribution ( f ), is related tothe nuclear activity, as indicated by the intrinsic hardX-ray flux (Lutz et al. 2004). Hard X-rays, unless ex-tremely obscured in fully Compton-thick objects, canprovide a direct view to the central activity, while theinfrared continuum is due to AGN emission reprocessedby dust, either in the torus or on somewhat larger scales.The mid-IR continuum has the advantage of showing nosignificant differences between type 1 and type 2 AGNand of being a good tracer of nuclear activity even inthose cases where hard X-rays are strongly absorbed.Low resolution mid-IR spectra of galaxies can be de-composed into three components (Laurent et al. 2000): acomponent dominated by the aromatic “PAH” featuresarising from photodissociation regions or from the dif-fuse interstellar medium of the host, a continuum risingsteeply toward wavelengths beyond 10 µ m due to HII re-gions and a typically flatter thermal AGN dust contin-uum present in active galaxies.In the range covered by the low resolution modules ofIRS (SL1 and SL2, from ∼ ∼ µ m as described inHouck et al. 2004), the AGN emission is most easily iso-lated shortwards of the complex of aromatic emission fea-tures (Laurent et al. 2000). We determine the continuumat 6 µ m rest wavelength and eliminate non-AGN emissionby subtracting a star formation template scaled with thestrength of the aromatic “PAH” features arising from thehost or from circumnuclear star formation. This methoddoes not require spatially resolving the AGN from thehost and, thanks to the high sensitivity of IRS, allowsthe detection of a weak AGN in the presence of strongstar formation in many nearby galaxies.For the purposes of this work, the use of mid-IR spec-tra, in particular of f , is the best way to quantify theAGN fraction for several reasons: (1) in order to derivea local luminosity function, we need a large sample: thenumber of objects for which mid-IR spectra are availableis much greater than the number of sources with X-raydata of sufficient quality; (2) we are interested in thecontribution of AGN to the infrared emission: L ismore strongly correlated with the AGN’s infrared emis-sion than, for example, L −
10 KeV ; (3) f is the onlyconsistent way to constrain the AGN contribution in ourcomparison sample at high redshifts, due to the limi-tations (in wavelength coverage and S/N) of the high-zspectra; (4) as explained above, f is also sensitive toCompton-thick AGN, not detected in the X-rays: thesekinds of galaxies become increasingly important at high- z , as shown both in mid-IR (Daddi et al. 2007) and X-ray(Gilli et al. 2007) surveys. Valiante et al. Sample selection and data reduction
The starting point for this study is the IRAS Re-vised Bright Galaxy Sample (RBGS), a complete flux-limited survey of all extragalactic objects with 60 µ mflux f > .
24 Jy (Sanders et al. 2003). The galax-ies with L TIR > L ⊙ were selected from this sam-ple and the mid-IR IRS data were obtained from the Spitzer archive for all of the them, with the exception of6 sources that were not yet observed or for which therewere no observation programs at all. Most of our ob-servations have been obtained within the program PID30323 (PI L.Armus). A comprehensive study of LIRGsand ULIRGs is being carried on by the Great Observa-tory All-sky LIRG Survey (GOALS ) and a more com-plete approach to derive AGN contributions from mid-IRspectra will be presented in forthcoming studies by theGOALS team. Additional data have been retrieved fromother Spitzer programs .Because of the limit in flux, this sample misses themost luminous infrared galaxies that are rare in the lo-cal universe. In fact, it does not include any galaxy with L > . L ⊙ , excluding all HYLIRGs, for example, fromthe analysis. In order to enlarge our range in luminos-ity, in particular including representatives of the mostluminous tail of the local galaxy population, we added Spitzer spectra for 14 galaxies from the sample of 16 byTran et al. (2001). These latter galaxies were selectedfrom surveys at different limits in f and with a prefer-ence for high-luminosity targets. They do not introduceany bias in mid-IR spectral properties or AGN content.In particular, no AGN-related IRAS color criteria, likethe f /f ratio, were applied in the selection.The total sample includes 211 sources spanning a lu-minosity range 10 . ≤ L TIR ≤ . L ⊙ (see Fig. 1 andTab. ?? ).We reduced the data as follows. We subtracted, foreach cycle, the two nod positions of the basic calibrateddata (BCD) frames. In the difference just calculated,we replaced deviant pixels by values representative oftheir spectral neighbourhoods. We subtracted residualwavelength dependent backgrounds, measured in source-free regions of the two dimensional difference spectra.In averaging all the cycles of the 2-dimensional sub-tracted frames, we excluded values more than three timesthe local noise away from the mean. The calibrated 1-dimensional spectra for the positive and the negativebeams were extracted and the two 1-dimensional spec-tra averaged in order to obtain the final spectrum, usingthe SMART analysis package (Higdon et al. 2004).We derive f by a decomposition over a range from5 . µ m to 6 . µ m rest wavelength. We fit the spectrumby the superposition of a star formation component dom-inated by the 6 . µ m PAH feature (M82 spectrum bySturm et al. 2000) and a simple linear approximation ofthe AGN continuum. Fig. 2 shows the decomposition ofthe spectrum for sources with different nature: NGC1614is a pure starburst; NGC1365 is classified as a Seyfert1.8 from its optical spectrum, has a Compton-thin AGNobserved in the X-ray emission (Risaliti et al. 2000) but http://goals.ipac.caltech.edu/ PIDs 14, 105, 666 (PI J.R.Houck), PID 1096 (PI H.Spoon),PID 3237 (PI E.Sturm), PID 3247 (PI C.Struck), PID 20549 (PIR.Joseph), PID 20589 (PI C.Leitherer) shows also PAH features in the mid-IR spectrum; MCG-03-34-064 has a typical AGN-dominated mid-IR spec-trum; the galaxy pair NGC3690/IC694 (Arp299) is clas-sified as Compton-thick AGN from its X-ray spectrum(Della Ceca et al. 2002) and its mid-IR spectrum revealsthat the nuclear activity is quite high but coincident withstrong star formation.The results of the fitting procedure are ˜ f , repre-sented in Fig. 2 by the thick vertical line, and ˜ f . ,an average flux density of the PAH component over therest wavelength range 6 . . µ m (see Tab. ?? ). Amore comprehensive and in-depth study of the PAH fea-tures of this sample will be presented in future papers bythe GOALS team. The spectrum of M82 is technicallywell suited for this decomposition but represents justa single object that may not be representative for starforming objects in general. We deal with this using theresults and approach of Lutz et al. (2004) (their section2.1). Specifically, we first use the M82 SWS spectrumfor decomposition and then correct for systematic differ-ences between M82 and a sample of 11 starbursts. Totake into account the dispersion among star-forming ob-jects in relative importance of all the components in themodel of Laurent et al. (2000), Lutz et al. (2004) havedecomposed this small sample of 11 star-forming objectsusing the same approach and same M82 template. Theyfind a mean residual of 0 . × ˜ f . with a disper-sion of 0.085. This corresponds to a a 6 µ m continuum inM82 that is fainter than in other starbursts. This findingagrees with the results of Brandl et al. (2006) who findmore 6 µ m continuum in their Spitzer starburst templatecomposed of 13 sources than in the M82 SWS spectrum.Decomposition of the Brandl et al. (2006) template us-ing our approach and correction for M82 systematics inthe way described results in an insignificant residual con-tinuum of less than 0.5% of the 6.2 µ m PAH peak, show-ing that also quantitatively the agreement between theirtemplate and our approach is good. Tab. ?? implementsour approach in the following way: the values f used for all further analysis are obtained after subtract-ing 0 . × ˜ f . from the direct fit result ˜ f .They are thus corrected for the slight systematic differ-ence between M82 and other starbursts. Error estimatesfor f are the quadratic sum of two components. Thefirst is a measurement error based on individual pixelnoise derived from the dispersion in the difference of ob-servation and fit. The second is 0 . × ˜ f . , thusconsidering the dispersion in the properties of the com-parison star forming galaxies.In summary, AGN 6 µ m continua obtained by our de-composition, are thus consistent with those obtained us-ing Brandl et al. (2006) starburst template, and in addi-tion include an estimate of the uncertainty arising fromthe spread of starburst properties. Distribution of the AGN contribution
Next we study the AGN contribution to L TIR inULIRGs and LIRGs as a function of luminosity. Ourresults will be used in the modeling of § µ m AGN continuum luminosity and the total infrared backward evolution model for infrared surveys 5luminosity ( νL /L TIR ) for five different luminositybins, spanning 10 . to 10 . L ⊙ . These distributions donot show the intrinsic AGN luminosity, but the fractionof the IR emission due to the contribution of an AGN,because of the normalization used.As shown from the distributions of Fig. 3 (solid andcross-hatched histograms) and also in the diagram ofFig. 7 (open squares), there is no evident trend with lu-minosity in the values of the means of the detections . Inthe bins at lower L TIR there are good statistics, but thepercentage of upper limits to L is rather high (up-per limits are represented by cross-hatched histograms).In the bins at higher L TIR , where the f detectionrates increase in a significant way (solid histograms rep-resent detections), the numbers of objects are more mod-est. In both cases, the averages obtained from detectionsonly (open squares in Fig. 7) have to be considered up-per limits to the real mean values of the distributions of νL /L TIR .The clearest trend in our measurements is the detec-tion rate as a function of L TIR . The decreasing frac-tion of upper limits clearly demonstrates that the distri-bution is changing across the luminosity bins and thatthere must be a trend also in the intrinsic distributionof νL /L TIR . With the assumption of a normal dis-tribution, it is possible to estimate the main parameters,mean and standard deviation, of the intrinsic distribu-tion of νL /L TIR by Monte-Carlo simulations, tak-ing into account the flux limit defining the RBGS sampleas well as the
Spitzer -IRS detection limit.
Simulating the distribution of the AGNcontribution
In order to reproduce our measurements, in terms ofdetection rate, mean and standard deviation of the detec-tions (hereafter d.r., m det , σ det ), using Monte-Carlo sim-ulations, we need to make assumptions about the intrin-sic distribution of νL /L TIR we are looking for andthe quantities that affect our upper limits. νL /L TIR is tightly correlated with two different quantities. L depends in a trivial way on f , i.e. on the sensitivityof the instrument. L TIR is related to the 60 µ m flux by re-lations well known and tested on IRAS data (Helou et al.1985; Dale et al. 2001).Our RBGS-based sample is selected at 60 µ m. There-fore, the first step of the simulation consists of generat-ing a population of sources having a 60 µ m flux densityspanning a range from 5 .
24 to 200 Jy and a distribution N ( f ν ) ∝ f − . ν (Sanders et al. 2003), expected for a com-plete sample of objects in a non-evolving Euclidean uni-verse that is a reasonable approximation for the relativelysmall redshift range covered by our sample. We derive L TIR for each source generated starting from its 60 µ mflux density, adopting L TIR ∼ × νL µ m as appropri-ate for our sample. Assuming that the distribution of νL /L TIR is gaussian with a certain mean and stan-dard deviation ( m , σ ), we calculate the expected value of f . At this point we apply a detection limit to f and calculate d.r., m det and σ det . The values of m and σ that best reproduce the measured d.r., m det and σ det foreach luminosity bin are assumed to be the parameters ofthe νL /L TIR intrinsic distribution.Our IRS data come from the public
Spitzer archive,consisting of several observing campaigns to varying depths. To assign a detection limit to f , we con-sider the values of 3 σ for the detections and the values ofthe upper limits for the non-detections in each luminositybin. We assume the median values of the distributions ofthese quantities as the detection limits of our measure-ments. These values span ∼ − × objects. This number was chosen suchthat the sampled distributions converge. Each diagramshows the best gaussians that, given the detection limit,reproduce the observed values of d.r., m det and σ det . Asummary of the distribution parameters of measurementsand simulations is presented in Tab. 1. The values foundfor m and σ are also plotted in Fig. 7 (green stars). Here,a correlation between the contribution to the infraredluminosity due to the AGN and the total infrared lu-minosity is clearly visible, with a dispersion decreasingwith the luminosity. The best fit, νL /L TIR ∝ L α TIR ,gives α = 1 . ± . νL /L TIR distribution tobe gaussian for each luminosity bin. This is a reason-able first hypothesis and it reproduces the observationsin terms of detection rate, mean and dispersion of thedetections with a small numbers of parameters. Moreattention is needed when applying it to our Monte Carlosimulations of the infrared sky. The low AGN activ-ity end of the distribution, corresponding mostly to IRSnon-detections, is certainly poorly constrained, but suchdifferences between already weak and totally insignifi-cant AGN will not have a major effect on the predictionsof our model. Of course, using these distributions tostudy physics of low luminosity AGN in these systemswould be misleading. Differences at the high AGN activ-ity end can be more important by introducing or miss-ing very luminous mid-IR sources dominated by AGN.Comparing for example the log( νL /L TIR ) distribu-tion for the third luminosity bin in Fig. 3 and Fig. 4(11 . < log( L TIR ) < . νL /L TIR ) ∼ νL /L TIR ) end. We will return to this issuein § L TIR ∼ . L ⊙ ) and then the relation be-tween νL /L TIR and L TIR becomes flat (Eq. 1; seealso Fig. 7, green dotted line).
Evolution with redshift
It is interesting to check if and how the relation justfound between νL /L TIR and L TIR evolves with red-shift. Previous studies already showed that, even if inthe local universe ULIRGs tend to be AGN-dominatedat L TIR > . L ⊙ (Lutz et al. 1998; Veilleux et al.1999; Tran et al. 2001), this trend does not extend tohigher redshifts. For example, submillimeter galaxiesthat, with their typical infrared luminosity of ∼ L ⊙ ,are among the most luminous objects known, are mainlystarburst dominated. Mid-IR spectroscopy in fact showsthat the contribution to the total luminosity due toAGN continuum is small in most of the sources ob- Valiante et al.served with Spitzer
IRS at high redshifts (Lutz et al.2005a; Men´endez-Delmestre et al. 2007; Valiante et al.2007; Pope et al. 2008). These results are consis-tent with X-ray observations of the same population(Alexander et al. 2005).A systematic study at high redshifts is difficult with ex-isting data. Ideally, to be consistent with the locally de-rived IRAS luminosity function, a high- z sample shouldbe selected at 60 × (1 + z ) µ m, at the peak of the dustemission, but this cannot be done with our pre- Herschel observations.Our high redshift sample includes all the sourcesobserved by Men´endez-Delmestre et al. (2007),Valiante et al. (2007), Pope et al. (2008) andBrand et al. (2008) for which the observations cover5 . − . µ m in the rest frame (see Tab. 1). Objectsfrom the first three works are submillimeter galaxies,some of them radio pre-selected. These galaxies aremainly starburst-dominated, but can present continuumemission in their mid-IR spectra, even if AGN is notthe dominant power source. Galaxies from Brand et al.(2008) are optically faint objects selected at 70 µ m:their mid-IR spectra show that they can be either PAHor absorption dominated. All the IRS data have beenretrieved from the Spitzer archive and reduced using thesame techniques described in § L in a consistent way for all the objects. Theresults of the fit process are listed in Tab. 1. The valuesof L TIR have been taken from the works cited above.Fig. 5 shows the distributions of the νL /L TIR .Unfortunately, the statistics are rather modest for thebins at lowest and highest luminosities. The redshiftrange spanned by this sample is very broad (0 . 1, while themost luminous bins include the most distant sources. Wetreat the whole sample as a unique population, becausethe inclusion of the z . νL /L TIR distribution in any of the lu-minosity bins.In order to characterize the intrinsic distribution of νL /L TIR , we run Monte-Carlo simulations with thesame techniques used for the local galaxies explained in § µ mflux density spanning 40 to 1100 mJy. We calculate theexpected value of the 60 µ m rest-frame flux density fromthe 850 µ m and 70 µ m fluxes, respectively for the sub-millimeter and the 70 µ m selected sources, assuming agrey-body spectrum with emissivity index β = 1 . T d = 32K. The values obtained give us the range of the60 µ m population to simulate. Because of the low countsin the bin at highest luminosity, we are able to repro-duce only the distributions of the first two bins, shownin Fig. 6.At high redshift, AGN contributions to the total in-frared luminosity differ from local objects. For L TIR ∼ L ⊙ , the AGN contribution is consistent with the lo-cal universe distribution given the small number statis-tics. For the better populated 12 . < log( L TIR ) < . νL /L TIR relation. Forour subsequent simulations, we simply assume that therelation is the same as derived locally, with a changein the behaviour of the flat part: while in the localuniverse the increase of the AGN contribution stops at L TIR ∼ . L ⊙ , for distant galaxies it occurs at lowerluminosities ( L TIR ∼ L ⊙ ). It can be summarized as: νL µ m /L TIR ∝ (cid:26) L . ( L TIR < L f ) L . ( L TIR ≥ L f ) ; (1) L f = (cid:26) . L ⊙ ( z < . L ⊙ ( z ≥ . Discussion The different role of AGN at low and high redshiftfor objects in the same L TIR class should not be a sur-prise. Fig. 7 can be read in two (equivalent) ways:(1) the contribution of AGN to L TIR is lower at high- z with respect to the local universe for sources with L TIR ∼ L ⊙ , or (2) the AGN content of high- z sourceswith L TIR ∼ L ⊙ is similar to that of local galaxieswith L TIR ∼ L ⊙ . This latter conclusion points outan analogy in the properties of SMGs and local ULIRGs.Such an analogy has been already discussed byTacconi et al. (2006): from CO observations of a sam-ple of 14 SMGs, they conclude that the density of themolecular gas, the luminosity surface density and thetemperature of the dust are the same in the two pop-ulations and SMGs are similar to local ULIRGs merg-ers, suitably scaled for their larger masses, luminositiesand star formation rates, as well as their greater gasfractions. In spatially resolved observations, SMGs andULIRGs show considerable similarities in their kinemat-ics (Tacconi et al. 2008): the SMGs presenting multiplecomponents are interacting systems and are similar to lo-cal double-nucleus ULIRGs, while SMGs showing char-acteristics of a rotational star-forming gas disk presentthe same surface and volume densities, for example, ofthe extremely compact ULIRG Arp220. All SMGs stud-ied with sub-arcsecond millimeter interferometer so far,appear to be major mergers in different stages, similarto local ULIRGs.Our analysis, nevertheless, is model-oriented: its pur-pose is not to find a theoretical explanation for the dif-ferent populations observed, but to develop an evolution-ary model based on observations capable of reproducingavailable data and to help plan forthcoming surveys withnew instruments. THE MODEL: NUMBER COUNTS AND REDSHIFTDISTRIBUTIONS OF INFRARED SOURCES The strong evolution of infrared galaxies:observational evidence There has been, in the past years, strong observationalevidence indicating high rates of evolution for infraredgalaxies.First, galaxy evolution can be observed through its im-print on the far-IR extragalactic background. Weaklyconstrained even as recently as the end of the 1990s, var-ious observations now measure or give upper/lower lim-its on the background from the ultraviolet (UV) to the backward evolution model for infrared surveys 7millimeter waveband (e.g. Hauser & Dwek 2001). Datashow the existence of a minimum between 3 and 10 µ mseparating direct stellar radiation from the infrared partdue to radiation re-emitted by dust. This re-emitteddust radiation contains a comparable integrated poweras the optical/near-IR: this amount is much larger thanwhat is measured locally ( ∼ 30 per cent). The CIB isthus likely to be dominated by a population of stronglyevolving redshifted infrared galaxies. Since the long-wavelength spectrum of the background is significantlyflatter than the spectrum of local star-forming galaxies,it constrains the far-IR radiation production rate history(Gispert et al. 2000). The energy density must increaseby a factor larger than 10 between the present time andredshift z ∼ − z ∼ µ m have resolved a frac-tion of the CIB into discrete sources. For all surveys,number counts indicate a very strong cosmological evo-lution of infrared galaxies, not only in the total powerradiated but also in the shape of the LF. This is par-ticularly obvious at submillimeter wavelengths, wherethe background is dominated by high-luminosity galaxies(SCUBA and MAMBO sources). The high rates of evolu-tion exceed those measured in other wavelength regimesas well as those observed for quasars and active galacticnuclei (AGN).Finally, high rates of evolution are suggestedby the detection of Poissonian fluctuations of theCIB at a high level at 60 and 100 µ m withIRAS (Miville-Deschˆenes et al. 2002), and 170 µ m withISOPHOT (Lagache et al. 2000; Matsuhara et al. 2000).For example, the constraints given by Matsuhara et al.(2000) on the galaxy number counts indicate the exis-tence of a strong evolution in the counts. SEDs Most existing backward evolution models ( § 1) sepa-rate the infrared sources into different populations. Forexample, Franceschini at al. (1988) consider three differ-ent populations: (1) normal late-type galaxies, (2) in-teracting/starburst galaxies and (3) galaxies with AGN.IRAS studies showed that galaxies of different nature inthe local Universe have different spectral energy distri-butions (SEDs), e.g different values of the f /f ra-tio, so usually each population is associated with a par-ticular SED family. The different populations of themodel can evolve as a single population (Xu et al. 1998;Xu 2000) or at different rates, assuming a local LF foreach component (Xu et al. 2001). Nevertheless, it isknown that infrared galaxies, in particular ULIRGs, canhost both starburst activity and an AGN (Genzel et al.1998; Lutz et al. 1998) and that galaxies classified asSeyfert show starburst tracers at infrared wavelengths(Lutz et al. 2004), while most of the far-IR emission inlocal QSOs is powered by starbursts (Schweitzer et al.2006).The SEDs can either be constructed empiri-cally, considering the individual components respon-sible for the overall emission from the galaxies(e.g. Rowan-Robinson 1992), or using radiative transfermodels (e.g. Efstathiou et al. 2000). While a complicatedarray of dust properties contributes to the SED of each galaxy, studies of IRAS galaxies have typically reducedthe description to a best-fit single dust temperature, T d ,with a one-to-one mapping to the f /f flux ratio(hereafter referred to as R (60 , R (60 , R (60 , R (60 , Spitzer data (e.g. Lagache et al. 2004). There is al-ready evidence for a change in SED properties with red-shift, with the intrinsic SED shapes of high z SMGs re-sembling lower luminosity local objects rather than simi-larly luminous local HYLIRGs, but a deep study of SEDsand temperatures of high-redshift galaxies is still missing.The main problems are the degeneracy of dust tempera-ture with redshift (Blain et al. 2002) , the limited sampleof sub-millimeter galaxies with measurements of spectro-scopic redshifts (e.g. Chapman et al. 2005; Valiante et al.2007; Pope et al. 2008), and the selection effects on dusttemperature induced by current methods to select highredshift infrared galaxies.In building the SEDs for our model, we will take intoaccount the coexistence of starburst and AGN in thesame galaxies, their varying contribution with L TIR , andthe spread in the luminosity-temperature ( L − T ) rela-tion. AGN are hosted by infrared galaxies of differentluminosities and their contribution to the total infraredemission can be either negligible or predominant. Westudied and discussed the problem of the AGN contribu-tion in § 2. We found that, on average, the contributionto L TIR due to an AGN in IRAS galaxies is proportionalto L . of the host. Even if the dispersion in the relationallows the presence of low luminosity AGN-dominatedobjects (e.g. NGC1068, log( L/L ⊙ ) = 11 . L ) = 12 . L ⊙ ), this correlation holds over two-orders-of-magnitude, log( L ) ∼ − L ⊙ . In § 2, weestimated the AGN contribution measuring the emissionat 6 µ m, after the subtraction of a starburst-like spec-trum; this method does not require to spatially resolvethe AGN within the host.In the simulation, we go in the opposite direction.Each source is assigned a redshift and a total infraredluminosity, ˜ L TIR , consisting of both the starburst ( ˜ L SB )and AGN ( ˜ L AGN ) contributions ( ˜ L TIR = ˜ L SB + ˜ L AGN ).Starting from ˜ L TIR , we calculate the expected value h νL µ m /L TIR i for that luminosity, using Eq. 1 (see alsoFig. 7).Then, we generate a random value of νL µ m /L TIR assuming its distribution is a gaussian with mean h νL µ m /L TIR i and sigma corresponding to the luminos-ity bin where ˜ L TIR lies (see Tab. 1). From this randomvalue we derive ˜ f for the redshift of the source gener-ated. Last, we scale an AGN template in order to achieve Valiante et al.the desired flux density ˜ f at 6 µ m. The AGN tem-plate used is the model calculated by Efstathiou et al.(1995) for the nuclear infrared continuum spectrum ofthe Seyfert galaxy NGC1068. The luminosity of thescaled template, ˜ L AGN , is the contribution of the AGNto ˜ L TIR .Spectral templates for starbursts are taken from theDale et al. (2001) (see also Dale & Helou 2002) cata-log, divided into 64 classes from R (60 , . 29 to1.64, corresponding roughly to single-component dust-temperature models of 23 − 45 K. SEDs are normal-ized by integrating each spectrum over the 8 − µ mrange and scaling to the expected value of L TIR forits R (60 , R (60 , C ∗ × (cid:18) L ′ L SB (cid:19) − δ × (cid:18) L SB L ′ (cid:19) γ (2)with γ = 0 . δ = − . C ∗ = − . 50 and L ′ = 5 . × L ⊙ . In order to fit this model we follow an identicalmethodology to Chapin et al. (2009), calculating L TIR by integrating the same Dale & Helou (2002) SED cat-alogue, normalized to IRAS 60 and 100 µ m data. In thisway we obtain an SED family in the range 10 − . Hzspanning luminosities log( L ) = 8 − . L ⊙ . During thesimulation, in Eq. 2 we use the starburst contribution˜ L SB to the total infrared luminosity only, which we ob-tain by subtraction ˜ L AGN from ˜ L TIR . From Eq. 2 weobtain the corresponding mean h R (60 , i for the ˜ L SB of the generated source. Then, we calculate the spreadin the L − T relation and generate a random value of R (60 , h R (60 , i (see Fig. 8). Finally, we rescaleto the correct ˜ L SB the Dale & Helou (2002) library SEDassociated to the R (60 , L − T relation.To obtain the final SED, we add the two spectra de-rived for the AGN and starbursts parts. Fig. 9 showsexamples of SEDs with different L TIR and different AGNcontributions. Because the model takes into accountboth the average relation and its dispersion, it permits awide range of combinations of L TIR and AGN emission. Model parameters In this section we describe the algorithm to coherentlymodel the number counts in different bands. A backwardevolutionary model requires some basic “ingredients”: • the spectral energy distributions, already describedin § L TIR , R (60 , • a cosmological model describing the accessible vol-ume in a survey solid angle as a function of redshift.This is now well constrained from WMAP results(Spergel et al. 2007). A concordance cosmology of H = 75 km s − Mpc − ,Ω M = 0 . Λ = 0 . • a luminosity function at z = 0, e.g. the most recentderivation of infrared luminosity function for IRASgalaxies by Chapin et al. (2009) • the evolution functions for density and luminosityevolution.The solid angle observed (and therefore the numberdensity of the sources) is geometrically defined by the as-sumed Universe model. The predicted number of sourcesin a given redshift interval ( z − . z, z + 0 . z ) and in agiven infrared luminosity interval ( L − . L, L + 0 . L )is given byd N ( L, z ) = g ( z ) φ ( L/f ( z )) d V d z d L d z (3)where φ is the local luminosity function and f ( z ) and g ( z )are respectively the luminosity evolution function andthe density evolution function and d V / d z is the comovingvolume element.The luminosity function assumed at z = 0 is a powerlaw parametrization of IRAS galaxies. As explained byChapin et al. (2009), the LF calculated accounts for atemperature bias in the underlying 60 µ m flux-limitedsample, as well as luminosity evolution. The LF wascomputed using the 1 /V max method (Schmidt 1968) witha modified formalism. The parametric form of the LF forthe L TIR is given by: φ ( L ) = ln(10) Lρ ∗ × (cid:18) LL ∗ (cid:19) − α × (cid:18) LL ∗ (cid:19) − β (4)with α = 2 . β = 2 . ρ ∗ = 6 . × − Mpc − L − ⊙ and L ∗ = 6 . × L ⊙ (see Fig. 10). As in Eq. 2 wehave re-calculated the fit for L TIR , since the results inChapin et al. (2009) are for narrower-bandwidth FIR lu-minosities.There exist observational evidences for strong evolu-tion in the infrared galaxy population (see § µ m by SCUBA (Archibald at al. 2001).Our model assumes power law luminosity, f ( z ), anddensity, g ( z ), evolutions. The following functional forms backward evolution model for infrared surveys 9are adopted: f ( z ) = (cid:26) (1 + z ) n ( z < z )(1 + z ) n ( z ≥ z ) (5) g ( z ) = (cid:26) (1 + z ) m ( z < z )(1 + z ) m ( z ≥ z ) (6)Both functions are assumed to be continuous at z and z .The algorithm calculates the number of sources pre-dicted in each redshift bin (see Eq. 3), introducing apoissonian scattering to the expected value. Then, avalue ˜ L TIR is associated to each source, by Monte-Carloextraction, according to the luminosity function calcu-lated for the redshift of the source, φ ( L/f ( z )). ˜ L TIR is then split into two parts: the fraction due to AGNcontribution ( ˜ L AGN ), calculated following Eq. 1, and theremaining part associated to starbursts ( ˜ L SB ). The SEDis thus “built” (see § L AGN with the starburst tem-plate corresponding to ˜ L SB , after the application of ascatter in the L − T relation, and scaled to the appro-priate luminosity distance. The sources’ flux densitiesin different bands are then calculated by convolving theredshifted SED with the bandpasses of the filters. Inthis way, we effectively simulate a virtual sky for a givenevolution model. For each source, we know not only theflux densities for several bands and the L TIR , but alsothe temperature of the dust and the L SB /L AGN ratio.The latter information, in particular, allows us to followthe co-evolution of starbursts and accretion, exploringphotometric indicators to select different populations ofinfrared sources. RESULTS AND COMPARISONS WITH AVAILABLESURVEYS Backward evolution models for the infrared spectralrange include a significant number of tunable parame-ters, and using the proper observational constraints toadjust these parameters is the key to the successful use ofsuch models. Traditionally, the first quantities to fit arethe total IR/submm background as measured by COBE(CIB), and the number counts at mid-IR to submillime-ter wavelengths. With improving identification and fol-low up of SCUBA, ISO and Spitzer sources, redshift dis-tributions of different IR-selected populations now pro-vide further powerful constraints. We are now about toenter the next level of physical characterization includ-ing determination of full rest frame far-IR SEDs, a fieldthat will be given a big boost with BLAST and Herschel ,and better disentangling the role of star formation andAGN in high redshift infrared galaxies. The model pre-sented here is designed to allow comparison with this newtype of constraints on SED and AGN content. By implic-itly coupling AGN evolution to infrared galaxy evolution,this model is also amenable to comparison with AGN sur-veys at other wavelengths, if making use of assumptionsabout AGN SEDs and obscuration.The most relevant results of the last years are relatedto surveys carried out by SCUBA and Spitzer . Thus,we will compare our results with the data available forthese instruments. Hereafter, we will use number counts from Papovich et al. (2004) and Shupe et al. (2008) forthe 24 µ m sources, Frayer et al. (2006a,b) for the 70 µ msources, Frayer et al. (2006a) for the 160 µ m sources,Coppin et al. (2006) for the 850 µ m sources. Red-shift distributions are from Wuyts et al. (2008) (24 µ msources) and Chapman et al. (2005) (850 µ m sources).Some of the models simulated are also compared with theCIB measurements (Fixsen et al. 1998; Lagache et al.2000; Renault et al. 2001; Miville-Deschˆenes et al. 2002;Lagache et al. 2003). All the simulations are run ona field of view of 10 deg , in order to be comparablewith the main infrared surveys and sufficiently samplethe high luminosity end. The redshift range covered is0 ≤ z ≤ 8. In the plots of the redshift distributions(see Fig. 11 − Fig. 14), the available data have been ap-propriately scaled in order to be comparable with thesimulations.In order to find the best set of parameters( n , n , z , m , m , z ) fitting available data we have runsimulations on small fields of view, changing slightly eachparameter “by hand”. Changing the parameters one byone allows us to understand how each parameter influ-ences the final results: we have used this approach tohave full control of the simulation and to understandthe role and weight of each parameter. During the opti-mization process, we have considered differently the sev-eral constraints given from current surveys, depending ontheir reliability: first, we have worked to reproduce themost accurate measurements and only later we have triedto fit the less constrained results. For example, SHADES(Coppin et al. 2006) has made accurate measurementsof the number counts at 850 µ m down to faint fluxes( ∼ Spitzer has provided very sensitive mea-surements with minimum uncertainties of the numbercounts at 24 µ m, either with observations on large scales(49 deg , Shupe et al. 2008) or very deep, down to fluxes ∼ . z . In fact, the radio-selectedsample of Chapman et al. (2005) is biased towards z . z ∼ µ m sources are obtained from photometric redshiftson a K S -band selected galaxies in the GOODS-S field(Wuyts et al. 2008): cosmic variance can have an impor-tant role on the redshift distribution in a field so small, aswell as the selection and the choice of the SEDs for pho-tometric redshifts. All these effects could overestimatethe number of low redshift objects. Finally, the surveysat 160 µ m (Frayer et al. 2006a) are still rather shallow.Below, we present a summary of the optimization pro-cess, showing the main steps followed to obtain the bestparameter set ( n , n , z , m , m , z ). A summary of theparameters for all the models discussed are shown inTab. 1. Pure starbursts SEDs As a starting point, we compare the observations witha model considering only the starburst part of the SED,including its scatter in the L − T relation (model M1).The submillimeter galaxies are starburst dominated, sowe look for the best parameter set reproducing at least0 Valiante et al.the 850 µ m number counts. We assume a strong evolu-tion in luminosity ( n = 3, n = 0) and a weaker evolu-tion in density ( m = 1, m = 0), with z = z = 2. Theresults are shown in Fig. 11. The bright 850 µ m countsare well reproduced. However, all the other counts areunderestimated: using this evolutionary model, there isstill room to introduce an SED taking into account theAGN contribution. The redshift distribution at 24 µ mis quite well reproduced, while at 850 µ m the modeledredshift distribution shows a strong excess at z & . Adding the AGN contribution In this second step, we build SEDs as described in § § 2, we assume agaussian distribution in each infrared luminosity bin (seeTab. 1). The results, assuming the same evolution ap-plied previously, are shown in Fig. 12 (M2). Comparedto model M1, the 24 µ m number counts are now betterreproduced, even if the shape of the peak is slightly dif-ferent. The redshift distribution at 24 µ m has becomeworse: the AGN contribution has introduced a large ex-cess at z ∼ 2. The parameter set used up to now is nomore acceptable and we have to explore the parameterspace further, in order to reproduce all the observations.Comparing different kinds of evolution, we learn thatwe need a strong luminosity evolution in order to repro-duce the number counts at 24 µ m and the 850 µ m, butthis always introduces an overprediction of the 24 µ m z ∼ µ m thatis very similar to the observed one, but we completelylose the peak in the 24 µ m number counts. We note,moreover, that a shift of the peak of the density evolu-tion, z , towards lower redshifts, causes a change in theshape of the redshift distribution at 24 µ m, reducing thenumber of the z ∼ µ mnumber counts. Last, if we want to reduce the numberof high redshift 850 µ m sources, we need a luminosityevolution peaking at z and then dropping down ( n hasto be negative). Fig. 13 shows number counts and red-shift distributions for a model with n = 3 . n = − z = 2 . m = 1, m = − . z = 1 (M3). Dataat 24 µ m are well reproduced, in particular the predictedpeak of the number counts has now the same shape of themeasurements even if it is slightly lower, while at 850 µ mthe simulation underpredicts the integral number countsbut reproduces the redshift distribution. Evolving the L − T relation In this simulation, the starburst part of the SED stilldoes not vary from the local to the distant universe. Thislast point is in contradiction to recent works alreadymentioned in § L − T relation, accounting for the cooler temperatures ofhigh redshift infrared galaxies. We modify Eq. 2 in thefollowing way: R (60 , C ∗ × (cid:16) L ′ L SB / (1+ z ) . (cid:17) − δ × (cid:16) L SB / (1+ z ) . L ′ (cid:17) γ (7)with all the parameters as in Eq. 2. This evolution isconsistent with the observations cited above. Fig. 14shows number counts and redshift distributions for thisnew model with n = 3 . n = − z = 2 . m = 1, m = − . z = 1 (M4). All the submillimeter mea-surements are now well reproduced, with perhaps someexcess in the number counts at strong fluxes, togetherwith the mid-IR data, whose number counts have im-proved and redshift distribution has not felt the effectof the new SEDs. The number counts in the far-IR, inparticular at 160 µ m, are nevertheless underestimated.Even if the peak of the 24 µ m number counts is well re-produced by the model M4, there is still some excess inthe number counts at strong fluxes, where AGN dom-inate, but the uncertainties are large given the smallnumber statistics. We compare the number of AGN pre-dicted by our model with the number of sources expectedadopting QSO luminosity functions. We select only thebrightest 24 µ m sources ( f ≥ L AGN /L SB ≥ L BOL ∼ νL ν at mid-IR rest wavelengths, we can calculate the ex-pected number of QSOs for the same L BOL and z rangesof our subsample. We integrate the luminosity functionof Hopkins et al. (2007) (“full” model, with pure lumi-nosity evolution and bright- and faint-end slopes evolv-ing with redshift) for 10 . ≤ L BOL ≤ . L ⊙ and0 . ≤ z ≤ . 3. The number of the sources selectedfrom our simulation does not exceed the expected num-ber given by the bolometric luminosity function. Never-theless, we already pointed out that the high part of the νL (AGN) / L TIR distribution is not well constrained (see § § Evidence for a local “cold” population? Our models were no able to reproduce the num-ber counts at 160 µ m, underestimating them by a fac-tor of ∼ 5. The ISOPHOT Serendipity and FIR-BACK surveys have revealed a population of nearby coldgalaxies (Stickel et al. 1998, 2000; Chapman et al. 2002;Patris et al. 2003; Dennefeld et al. 2005; Sajina et al.2006), under-represented in the 60 µ m IRAS sample.These objects are often associated with bright opticalspiral galaxies, and their far-IR colours ( f /f ∼ . µ m, similar to backward evolution model for infrared surveys 11that seen for example in the Milky Way galactic ridge(Serra et al. 1978).Lagache et al. (2003) implemented this class of objects(hereafter “cold galaxies”) in their evolutionary model:cold galaxies dominate the z = 0 luminosity functionat low luminosities and become less important at higherredshifts, because their evolution is passive and short( n = 1 , n = 0 , z = 0 . L TIR ∼ L ⊙ , while at larger lu-minosities it is quite small. Lagache et al. (2003) showthat the fraction of cold galaxies can contribute to thetotal number counts up to ∼ 50% at 170 µ m in the fluxrange where a comparison with measurements is possible(0 . − . µ m (Lagache et al. 2004, their Fig. 4).In order to introduce a similar population in our model,we have to assume that probably the scatter in the L − T relation calculated by Chapin et al. (2009) and imple-mented in our simulations needs some corrections in or-der to take into account the cold galaxies population.The lack of these corrections can explain why we un-derestimate the 160 µ m number counts. Moreover, anychanges of the adopted L − T relation and its scat-ter has to be constrained by the 850 µ m counts. The160 µ m counts need to be increased without creating ex-cess 850 µ m counts. This may support a mostly localchange. Moreover, colder galaxies with a higher 160 µ mflux have typically lower 70 µ m flux: any change intro-duced cannot be too much strong. Chapin et al. (2009)discussed two SMGs that are outliers respect to their L − T relation: they explain the exceptions assuming amismatching in the optical counterparts. Nevertheless,they do not exclude that both the galaxies are really atlow redshift and a population of cold galaxies do exist at z < L − T relation as-suming an asymmetric gaussian distribution. The orig-inal distribution calculated by Chapin et al. (2009) is astandard gaussian in log( R (60 , z < n = 3 . n = − z = 2 . m = 2, m = − . z = 1 (M5) are shown in Fig. 15. Finally, both thefar-IR (160 µ m and 70 µ m wavebands) number countsare reproduced within the uncertainties as well as thesubmm measurements. We note that a qualitatively sim-ilar need to add local cold sources arose when using inour models both the earlier IRAS-based determination oftemperature spread by Chapman et al. (2003) and therecent work of Chapin et al. (2009). We can considerM5 our best-fit model: all the number counts are wellreproduced (within a factor . Herschel surveys is given in § bottom ).A further constraint on our evolutionary model comesfrom the comparison with the current CIB measure- ments. Our model predictions for the CIB are derivedby summing up the flux densities of all sources for thebands in question. Fig 16 shows the predicted CIB inten-sity at specific wavelengths together with the comparisonwith present observations for the model M5: it agreeswith all the available measurements and limits. We alsocompare the IR luminosity function calculated at dif-ferent redshifts with measurements. Fig. 17 shows theluminosity function of our best model M5 fitting all datafrom low redshifts z = 0 . z =0 . , . , . , . , . z = 1 , z ∼ . STAR FORMATION HISTORY There has been much interest in the star formationrate and metal production history of the Universe sinceMadau et al. (1996) showed that they could be derivedfrom deep UV/optical cosmological surveys. Since thenthe so-called Madau diagram has been revised manytimes through various improvements, including (1) thecorrections for the effect of dust extinction on UV/opticalluminosities (e.g. Steidel et al. 1999), (2) results from lessextinction-sensitive Balmer line surveys (e.g. Yan et al.1999), and (3) results from submillimeter SCUBA sur-veys (e.g. Chapman et al. 2005).Adopting the conversion factor of Kennicutt (1998),SFR [M ⊙ yr − ] = 4 . · − × L TIR [erg s − ] , (8)we can convert the cosmic luminosity density evolutionof our model to a SFR curve. In the calculation, we takeinto account only the fraction on infrared luminosity dueto star formation ( L SB ). In Fig. 18 the results of the sim-ulation using model M5 are compared with the surveydata already shown by Chapman et al. (2005) (see theirFig. 12). The model is in very good agreement with theobservations for the full redshift range. In particular, for z & 4, the model results are similar to the results fromSCUBA surveys, while it is slightly lower than the mea-surements obtained from the UV/optical surveys. Thisdoes not surprising, since the model itself is constrainedby the infrared/submillimeter data. AGN CONTRIBUTION TO THE INFRARED EMISSION As discussed earlier, LIRGs and ULIRGs, with theirhuge infrared luminosities, correspond to an extremelyactive phase of dust enshrouded star formation and/orAGN activity. They become an increasingly significantpopulation at high redshift, representing an importantphase in the buildup of massive galaxy bulges and inthe growth of their central supermassive black holes. Tounderstand these processes, it is essential to separate therelative contributions of AGN and starburst activity tothe infrared luminosity of LIRGs and ULIRGs.Recent work has shown how the IRAC color-color di-agram and the MIPS 24 to 8 µ m color can be usedto identify AGN-dominated sources (Lacy et al. 2004;Sajina et al. 2005; Stern at al. 2005; Yan et al. 2004;Brand et al. 2006). Brand et al. (2006), in particular,have demonstrated how the 24 to 8 µ m flux ratio ( ζ ≡ log[ νf nu (24 µ m) /νf ν (8 µ m)]) can be used to disentanglethe contribution of AGN and starbursts to the total re-processed mid-IR ( ∼ − µ m) emission as a function2 Valiante et al.of the 24 µ m flux. We use this study as a first exam-ple of how to compare our model with observations thatconstrain the AGN content of high − z sources.Nevertheless, a photometric study of this type is sub-ject to several caveats. Many broad emission- andabsorption-line features are known to be present in themid-IR spectrum of ULIRGs (e.g. Houck et al. 2005;Yan et al. 2005), and this may affect ζ as function ofredshift. For example, at 1 . . z . . 7, the 24 µ m ob-served emission can be strongly attenuated by the sili-cate absorption feature at 9 . µ m. It is also possible thatan AGN could be heavily embedded in large amountsof cooler dust and could remain undetected in the 8 µ mband even though it dominates the bolometric infraredemission. In practice: on one hand, the galaxies, which inBrand et al. (2006) are classified as “AGN-dominated”,do host an AGN, but it is not necessarily the dominantsource powering the infrared emission; on the other hand,objects that are known to be dominated by AGN in themid-IR, like NGC1068, do not present the color ( ζ ∼ µ m selected samples of infrared sources will be impor-tant in putting more constraints on the AGN contribu-tion to the total infrared emission and in characteriza-tion of the 24 µ m selected population. Large Spitzer -IRSobserving programs covering sources with 24 µ m fluxesaround and below 1 mJy are under analysis and will bewell suited for this task. In combination with IRS sur-veys at brighter fluxes, they will trace the evolution ofthe AGN/SB ratio, strength of PAH emission and mid-IRopacities as a function of L TIR and z . In the meantime,results of photometric studies like the one of Brand et al.(2006) can be compared to our models (see Fig. 19 a ).The models developed here allow us to predict the frac-tion of AGN-dominated sources as a function of the mid-IR flux. In particular, it is possible to distinguish be-tween sources dominated by AGN at 24 µ m and sourcesAGN-dominated with respect to the total infrared lumi-nosity. This distinction is important when comparing theresults of the simulation with the observations: the infor-mation about f will be useful in comparison with theupcoming results from the IRS observations cited above,while the predictions related to L TIR will be confrontedwith Herschel results.At the moment, as already explained, the shape of theAGN distribution is not well constrained at the highestratios of AGN to total infrared luminosity. As an exper-iment, in order to reduce the number of luminous AGN,we change the assumed gaussian distribution introducinga cut for large values of νL µ m /L TIR : we reject all val-ues of νL µ m /L TIR more than 3 σ away from its expectedvalue. This model, again using the same parameter set( n = 3 . n = − z = 2 . m = 1, m = − . z = 1, model M6), does not significantly change thenumber counts and redshift distributions with respect tothe best model M5, so we do not show them here.Differences between the two models are evident whendoing different kinds of analysis. Fig. 19 a shows thefraction of all sources whose mid-IR emission is dom-inated by AGN ( f /f > 1) as a function ofthe 24 µ m flux density for model M5 ( solid histogram )and M6 ( dashed histogram ). The trend and the values obtained are similar to those measured by Brand et al.(2006) ( open squares ), even if the measurements have tobe compared carefully for the reasons explained above.Our predictions will be easily compared instead with theupcoming IRS results, that will help us to put more con-straints on the modeled AGN contribution. In fact, theAGN fraction of model M6, where the AGN contribu-tion has been slightly modified, is a bit lower comparedto model M5: only more constraints coming from mid-IRspectroscopic observations can discriminate between thetwo models M5 and M6. Fig. 19 b shows the fraction ofall sources whose total infrared emission is dominated byAGN ( L AGN /L SB > 1) as a function of the 24 µ m fluxdensity for the two models M5 and M6. As expected, thevalues of the fraction are lower, because the AGN startsto dominate the mid-IR part of the SED before beingdominant in the entire infrared range (see e.g. Fig. 9). OPEN ISSUES Our best-fit model M5 is able to reproduce most ofthe measurements available. Nevertheless, there are stillsome issues that need a deeper analysis and discussionand will be dealt with in future studies.The first point concerns the so called “local cold pop-ulation”. Even if there are several clues about the exis-tence of a population of colder IR galaxies not detectedby IRAS, only forthcoming surveys at longer wavelengthscarried out with Herschel will be able to reveal and char-acterize this kind of objects.The second point regards the distribution of the AGNcontribution and the fraction of AGN-dominated sources.We already pointed out that the models M5 and M6 pre-dict similar 24 µ m number counts, but different amount ofAGN-dominated sources. Galaxies at z ∼ K − . µ m colours than nor-mal galaxies: this is evidence for an AGN contribution tothe mid-IR continuum due to warm dust. The presence ofCompton-thick AGN is confirmed by the stacked Chan-dra X-ray data (Daddi et al. 2007). The sky density andthe volume density of this population of obscured AGNagree reasonably well with those predicted by the back-ground synthesis models of Gilli et al. (2007). In orderto discriminate between our best models, M5 and M6, weneed to compare the AGN-dominated sources predictedwith upcoming results with IRS. This will introduce afurther constraint, in addition with the bounds alreadyused to the total number counts, and it will enable thefinal model to reproduce the co-evolution of AGN andstarbursts in more detail.The third point regards the evolutions of the AGN − and Color − L TIR relations. It is clear that some evolu-tion is needed in order to reproduce existing observations,nevertheless the evolution implemented in this model isquite simple. More information regarding the AGN con-tent of high − z far-IR selected sources and the shape ofthe SED of SMGs will enable a better characterization ofthese details. Understanding evolution is intimately re-lated to the black hole growth and to the star-formationrate history of massive galaxies. PREDICTIONS FOR MULTIBAND HERSCHEL SURVEYS In this section, we show predictions using the best-fitmodel M5. We concentrate on future surveys with Her-schel , which will be launched in early 2009, in particular backward evolution model for infrared surveys 13on the PEP survey . All three PACS bands (70, 100,160 µ m) are considered. Predicted number counts at 70and 160 µ m are already shown in Fig. 14 and comparedwith the currently available data in those wavebands.In Fig. 20 top the predicted integral number counts at100, 250, 350, 500 µ m are plotted. While at 70 µ m theAGN-dominated sources can give a non-negligible con-tribution to the total number counts, in particular for f & 50 mJy, in the other Herschel bands the numbercounts are exclusively dominated by the starburst com-ponent.Fig. 20 middle shows predicted redshift distributions forthe deep ( f > . f > ′ × ′ ) and the large (85 ′ × ′ )fields are centred respectively on the GOODS-S and onthe COSMOS fields, in order to have the maximum avail-ability of multi-wavelength data and follow up opportuni-ties. The chosen fields are in fact fully covered in severalof the following: deep X-ray surveys, UV/optical/near-IR/IRAC imaging, HST imaging, Spitzer mid-IR sur-veys, submillimeter surveys and radio mapping. The pre-dicted distributions show a prominent peak at z ∼ z ∼ z ∼ 3, allowing a consistentstep forward in the study of the cosmic evolution of dustystar formation and of the infrared luminosity function, aswell as in the determination of the overall SEDs of activegalaxies.In Fig.20 bottom the predicted redshift and dust tem-perature distributions of galaxies selected in the PACSbands (100 and 160 µ m) for models M4 (no local coldpopulation) and M5 (with local cold population) areplotted. The selection is assumed to be the same asexpected in the CDFS field ( ∼ , f > 27 mJy, f > . z . ∼ 3. Dust temper-ature distributions are completely different: M4 predictsa mean temperature of T d ∼ 38 K, while M5 presents askewed distribution with a prominent peak at T d ∼ 23 Kand a tail up to T d ∼ 60 K. Multiband Herschel surveyswill measure these distributions. CONCLUSIONS The aim of this work has been to find a “backwardevolution” model that can viably fit the galaxy sourcecounts and redshift distributions from the mid-infraredto submillimeter wavelengths whilst not violating theconstraints set by the Cosmic Infrared Background mea-surements. Increasingly detailed observations are char-acterizing high redshift galaxies simultaneously at manyinfrared wavelengths, and are starting to unravel therole of AGN within them, suggesting a refinement ofprevious models. We adopted a Monte-Carlo-based ap-proach that considers the evolution of a coherent popu-lation of infrared galaxies with distribution functions de-scribing the contribution of AGN and the spread aroundthe luminosity-temperature relation for the far-infraredemission.As a necessary ingredient of this model, we have char- acterized the local distribution of AGN contributions asa function of total infrared luminosity. We have appliedspectral decomposition to a large sample of Spitzer -IRSspectra of ULIRGs and LIRGs with L TIR > L ⊙ to isolate the AGN 6 µ m continua. The distribution of L /L TIR changes with L TIR . We compare the AGNdetections and limits to simple Monte-Carlo simulations,making the assumption of an intrinsically gaussian dis-tribution, to quantify this increase with luminosity of theAGN contribution to the infrared luminosity. The bestfit, νL µ m /L TIR ∝ L α TIR , gives α = 1 . ± . 6. The rela-tion does not hold any more at high redshifts. The mostluminous high redshift infrared galaxies, L TIR ∼ L ⊙ ,show a small contribution from AGN, being mainly star-burst powered.The local start of our backward evolution modelscombines this result with the luminosity function andtemperature-luminosity relation of Chapin et al. (2009).We have explored several models considering evolution-ary variation of increasing subsets of the parameterswithin this framework. We find that we can satisfactorilyreproduce number counts and redshift constraints fromthe literature with such a population of infrared galaxiesevolving in luminosity and density. The contribution ofstarbursts and AGN varies with redshift and luminosity,as does the far-infrared dust temperature. The luminos-ity evolution is mainly constrained by number counts at24 and 850 µ m, while the density evolution is stronglyinfluenced by the 24 µ m redshift distribution. While thisevolution in luminosity, density, AGN contribution anddust temperature is clearly suggested by this sequenceof models, the detailed functional forms and parametersare necessarily uncertain with current data.As in previous work (e.g. Lagache et al. 2003), we haverealized a need to invoke a population of local cold galax-ies in order to reproduce the number counts at 160 µ m.Even with such clues about the existence of this popu-lation, its observational characterization remains incom-plete.The coherent consideration of both AGN and starburstin our model allows for direct comparison to observationsconstraining the role of AGN in high redshift infraredgalaxies. Trends with modelled mid-infrared flux fromour best fit model are consistent with photometric es-timates (Brand et al. 2006). Forthcoming IRS observa-tions of large flux-limited samples will both reduce theobservational uncertainties in characterizing the role ofAGN and allow to better constrain the role of the mostluminous AGN in our models.Characterization of the SEDs of z ∼ . − Herschel -PACS.4 Valiante et al. backward evolution model for infrared surveys 15 Steidel, C.C., Adelberger, K.L., Giavalisco, Dickinson, M., M.,Pettini, M. 1999, ApJ, 519, 1Stickel, M., et al. 1998, A&A, 336, 116Stickel, M., et al. 2000, A&A, 359, 865Sturm, E., Lutz, D., Tran, D., Feuchtgruber, H., Genzel, R.,Kunze, D., Moorwood, A.F.M., Thornley, M.D. 2000, A&A,358, 481Tacconi, L.J., et al. 2006, ApJ, 640, 228Tacconi, L.J., et al. 2008, ApJ, 680, 246Teng, Stacy H., Wilson, A.S., Veilleux, S., Young, A.J., Sanders,D.B., Nagar, N.M. 2005, ApJ, 633, 664Tran, Q.D., et al. 2001, ApJ, 552, 527Valiante, E., Lutz, D., Sturm, E., Genzel, R., Tacconi, L.J.,Lehnert, A., Baker, A.J. 2007, ApJ, 660, 1060Veilleux, S., Kim, D.-C., Sanders, D.B. 1999, ApJ, 522, 113Veilleux, S., et al. 2009, ApJ, in press, (arXiv:0905.1577) Weedman, D.W., et al. 2005, ApJ, 633, 706Werner, M.W. et al. 2004, ApJS, 154, 1Wuyts, S., Labbe, I., F¨orster Schreiber, N.M., Franx, M.,Rudnick, G., Brammer, G.B., van Dokkum, P.G. 2008, ApJ,682, 985Xu, C. et al. 1998, ApJ, 508, 576Xu, C. 2000, ApJ, 541, 134Xu, C., Lonsdale, C.J., Shupe, D.L., O’Linger, J., Masci, F. 2001,ApJ, 562, 179Yan, L., McCarthy, P.J., Freudling, W., Teplitz, H.I., Malumuth,E.M., Weymann, R.J., Malkan, M.A. 1999, ApJ, 519, L47Yan L. et al. 2004, ApJS, 154, 60Yan, L., et al. 2005, ApJ, 628, 604Yan, L., et al. 2007, ApJ, 658, 778 TABLE 1Full sample of infrared galaxies: results of the fit andcorrected f fluxes Name z log( L TIR ) ˜ f f L J mJy mJy Sources from RBGS catalog (Sanders et al. 2003) NGC0023 0.015 11.05 13.9 ± < ± < ± < ± ± ± < ± < ± < ± ± ± < ± < ± ± ± < ± < ± < ± < ± < ± < ± < ± < ± < ± ± ± < ± ± ± < ± ± ± < ± ± ± < ± ± ± ± ± < ± ± ± < ± < ± ± ± < ± ± ± < ± < ± < ± < ± ± ± < ± < ± ± ± < ± < ± < ± < ± ± ± ± ± ± ± < ± < ± < ± < ± ± ± < ± < ± < ± < ± < ± ± ± < ± < ± ± ± < ± ± ± < ± < ± ± backward evolution model for infrared surveys 17 TABLE 1 continued – Full sample of infrared galaxies: results of thefit and corrected f fluxes Name z log( L TIR ) ˜ f f L J mJy mJyESO432-IG006 0.016 11.02 18.7 ± < ± ± ± ± ± ± ± < ± < ± ± ± < ± < ± < ± ± ± ± ± < ± < ± ± ± < ± < ± < ± ± ± < ± ± ± ± ± < ± < ± < ± < ± < ± ± ± ± ± ± ± ± ± ± ± < ± < ± < ± < ± ± ± < ± ± ± < ± < ± < ± ± ± < ± < ± < ± < ± < ± ± ± < ± < ± < ± < ± < ± < ± ± ± < ± < ± < ± < ± < ± ± ± < ± ± ± ± ± < ± < ± < ± < ± ± ± < TABLE 1 continued – Full sample of infrared galaxies: results of thefit and corrected f fluxes Name z log( L TIR ) ˜ f f L J mJy mJyIRASF16399-0937 0.027 11.56 13.5 ± ± ± < ± ± ± < ± < ± < ± < ± < ± < ± < ± < ± < ± < ± < ± < ± < ± < ± < ± < ± < ± < ± < ± ± ± ± ± < ± < ± ± ± < ± ± ± ± ± < ± < ± < ± ± ± < ± < ± < ± < ± < ± ± ± < ± ± ± < ± ± ± < ± < ± < ± ± ± < ± ± ± < ± ± ± < ± < ± < Sources from Tran et al. (2001) IRASF00183-7111 0.327 12.77 45.1 ± ± ± ± ± ± ± ± ± < ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± < ± ± backward evolution model for infrared surveys 19 TABLE 2Distributions parameters of measurements and simulations log( L TIR ) Measurements Detect.limit Simulations L ⊙ m det σ det d.r.[%] mJy m det σ det d.r.[%] m σ ÷ ÷ ÷ ÷ ÷ TABLE 3High redshift sample of infrared galaxies: results of the fit andcorrected f fluxes Name z log( L TIR ) ˜ f ˜ f . f L J mJy mJy mJy Sources from Brand et al. (2008) ± ± ± ± ± < ± ± < ± ± < ± ± ± ± ± ± ± ± < ± ± ± ± ± < ± ± ± ± ± ± Sources from Pope et al. (2008) C3 1.88 12.78 0.15 ± ± ± ± ± < ± ± < ± ± < ± ± ± ± ± < ± ± ± ± ± < ± ± ± Sources from Men´endez-Delmestre et al. (2007) SMMJ030228+000654 1.41 13.44 0.07 ± ± < ± ± ± ± ± ± Sources from Valiante et al. (2007) SMMJ00266+1708 2.73 12.70 0.15 ± ± < ± ± ± ± ± < ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± TABLE 4Evolution models Model Description n n z m m z M1 SB 3.6 0.0 2.0 1.0 0.0 2.0M2 SB+AGN 3.6 0.0 2.0 1.0 0.0 2.0M3 SB+AGN 3.4 -1.0 2.3 2.0 -1.5 1.0M4 SB+AGN+ L − T evolution 3.4 -1.0 2.3 2.0 -1.5 1.0M5 SB+AGN+ L − T evolution+“cold spread” 3.4 -1.0 2.3 2.0 -1.5 1.0M6 SB+“cut” AGN+ L − T evolution+“cold spread” 3.4 -1.0 2.3 2.0 -1.5 1.0 backward evolution model for infrared surveys 21 Fig. 1.— Total infrared luminosity distribution of the galaxies of our sample, including galaxies from RBGS catalog (Sanders et al. 2003)and sources from Tran et al. (2001). Fig. 2.— Examples for the decomposition used to isolate the AGN continuum. Top panels : Spitzer IRS low resolution spectra of fourrepresentative galaxies. Bottom panels : Cutout of the region around the 6 . µ m PAH feature. Top continuous line = observed spectrum.Top dotted line = fit by the sum of the M82 spectrum and a linear AGN continuum. Bottom dotted line = fitted AGN continuum. Bottomcontinuous line = difference of observed spectrum and fitted PAH component. The thick vertical line indicates the 6 µ m AGN continuumflux density ˜ f . backward evolution model for infrared surveys 23 Fig. 3.— Distribution of νL /L TIR , for different luminosity bins. Solid histograms show the measurements, while the cross-hatchedhistograms indicate objects with upper limits. As indicated in each diagram, the mean values of the detections do not vary significantlybetween the different bins. What noticeably varies is the detection rate, spanning values from 23.3% to 100%. Fig. 4.— Simulations of νL /L TIR , assumed to be normally distributed, reproducing the same detection rates (d.r.), means ( m det )and sigma ( σ det ) of the detections as the data. The detection limits assumed for f in each bin is reported in Tab. 1. Symbols havethe same meaning as in Fig. 3. Each simulation is made for a population of 2 × objects. Each diagram reproduces the values observedin the corresponding panel of Fig. 3 and shows also mean ( m ) and sigma ( σ ) of the adopted distribution. backward evolution model for infrared surveys 25 Fig. 5.— Distribution of νL /L TIR for different luminosity bins. The sample lies at redshift 0 . < z < . 35. Symbols have thesame meaning as in Fig. 3. Each diagram also shows the mean, the dispersion and the detection rate of the distribution. Fig. 6.— Simulations of νL /L TIR normal distributions reproducing the same detection rates (d.r.), means ( m det ) and sigma ( σ det )of the detections obtained from the data of the high redshift sample, assuming a detection limit of 0 . 13 mJy. Each diagram reproduces thevalues observed in the corresponding panel of Fig. 5, while the last luminosity bin has not been simulated because the statistic is too lowfor the fit to converge. Symbols have the same meaning as in Fig. 3. Each diagram also shows mean ( m ) and sigma ( σ ) of the adopteddistribution. backward evolution model for infrared surveys 27 Fig. 7.— νL /L TIR , vs. L TIR . Dots and small arrows represent measurements and upper limits of all the galaxies in our local sample.The large open squares are the means of the detections ( m det ) in bins of luminosity (see Tab. 1 and Fig.3). Due to the large number ofnon-detections, they represent only an upper-limit on the true log νL /L TIR for almost all the bins. The green stars and the associatedbars show means ( m ) and sigma ( σ ) of the gaussian distributions adopted in the simulations in order to reproduce detection rates (d.t.),means ( m det ) and sigma ( σ det ) of the detections from the measurements (see Tab. 1 and Fig.4). The best fit, νL /L TIR ∝ L α TIR ,gives α = 1 . ± . L TIR . The red stars and their bars showmeans and sigma of the gaussian distributions adopted to reproduce measurements of the high redshifts sample (see Fig. 5 and Fig. 6).The relation found for the local galaxies no longer holds. Fig. 8.— Relation between the color R (60 , L TIR (solid and dotted lines give the mean and 1 − σ envelope of Eq. 2 respectively). Fig. 9.— Some examples of adopted SEDs (solid lines) obtained adding a star-forming galaxy spectrum (dashed lines) (Dale et al. 2001;Dale & Helou 2002) and an AGN template (dotted lines) (Efstathiou et al. 1995). The flux density is expressed in Jansky and is scaledfor a distance of 100 Mpc. In each diagram we indicate L TIR and the ratio between the luminosity of the AGN and the luminosity ofthe starburst spectra. For each luminosity class, two examples are shown: an AGN-dominated SED ( L AGN /L SB ≥ 1, left column) and astarburst-dominated SED ( L AGN /L SB < 1, right column). backward evolution model for infrared surveys 29 Fig. 10.— The ‘double power-law” parametrization of the TIR luminosity function. Fig. 11.— Number counts at 24, 70, 160, 850 µ m and redshift distributions at 24 and 850 µ m, considering evolutions with n = 3 . n = 0, m = 1, m = 0, z = z = 2, without AGN contribution (M1, solid lines). Data are from Papovich et al. (2004) ( filled circles ) andShupe et al. (2008) ( open squares ) for the 24 µ m number counts, Frayer et al. (2006a,b) (70 µ m and 160 µ m counts), Coppin et al. (2006)(850 µ m counts), Wuyts et al. (2008) (24 µ m redshifts) and Chapman et al. (2005) (850 µ m redshifts). Euclidean normalized differentialcounts are shown at 24, 70, 160 µ m, while integrated counts are shown at 850 µ m following the style of the papers in which they wereoriginally presented. Redshift data ( open histograms ) have been scaled for comparison with the simulations ( filled histograms ). backward evolution model for infrared surveys 31 Fig. 12.— Same as Fig. 11, but now with AGN contribution (M2). AGN-dominated galaxies ( L AGN /L SB ≥ 1) are represented withdashed line, while starburst dominated galaxies ( L AGN /L SB < 1) are in dotted line. Fig. 13.— Number counts at 24, 70, 160, 850 µ m and redshift distributions at 24 and 850 µ m, considering evolutions with n = 3 . n = − z = 2 . m = 2, m = − . z = 1, with AGN contribution (M3). Symbols are like in Fig. 12. Data are like in Fig. 11. backward evolution model for infrared surveys 33 Fig. 14.— Same as Fig. 13, assuming an evolution in the L − T relation (M4). Fig. 15.— Same as Fig. 13, assuming an evolution and a spread towards lower temperatures in the L − T relation (M5). backward evolution model for infrared surveys 35 Fig. 16.— Cosmic background from mid-IR to millimeter wavelengths for model M5. The CIB derived from the models at 24, 70,100, 160 and 850 µ m is shown by big open squares. The analytic form of the CIB at the FIRAS wavelengths is from Fixsen et al. (1998)( solid and dashed lines ). Measurements are at DIRBE wavebands (140 and 240 µ m, plus signs , Lagache et al. 2000), at 100 µ m ( triangles ,Lagache et al. 2000; Renault et al. 2001) and at 60 µ m ( triangles , Miville-Deschˆenes et al. 2002). See Lagache et al. (2003) for referencesabout the limits at shorter wavelengths. Fig. 17.— Luminosity function of model M5 calculated at different redshifts, compared with available measurements. Data arefrom Huang et al. (2007) ( squares ), Le Floc’h et al. (2005) ( triangles ), Magnelli et al. (2009) ( circles ), Caputi et al. (2007) ( diamonds ),Chapman et al. (2005) ( stars ). backward evolution model for infrared surveys 37 Fig. 18.— Histogram showing the SFR evolution of the model M5 compared with the observational results shown by Chapman et al.(2005) (see their Fig. 12). Yellow triangle : radio surveys; red circles : mid-IR surveys; green stars UV and optical surveys, corrected fordust extinction; violet squares : SMGs corrected for completeness. See Chapman et al. (2005) for a full list of references. Fig. 19.— Fraction of all sources whose mid-IR emission ( a ) and total infrared emission ( b ) are dominated by AGN as a function of f for model M5 ( solid histogram ) and M6 ( dashed histogram ). Redshift range of the sources is 0 ≤ z ≤ Open squares represent the resultsfrom Brand et al. (2006). backward evolution model for infrared surveys 39 Fig. 20.— Predictions using the best-fit model M5 for the upcoming Herschel surveys. Top : integral number counts at 100 µ m ( solidline ), 250 µ m ( dotted line ), 350 µ m ( dashed line ), 500 µ m ( dot-dashed line ). Middle : redshift distributions of PEP surveys in the deep(10 ′ × ′ , f > . left ) and large (85 ′ × ′ , f > right ) fields. Bottom : expected distributions for the survey in the CDFSfield for model M4 ( dotted line ) and M5 ( solid line ) of redshifts ( left ) and temperatures ( rightright