Astronomy&Astrophysicsmanuscriptno.Bethermin˙model (cid:13)c ESO2011 February1,2011 Modeling the evolution of infrared galaxies: A Parametric backwards evolution model M.Be´thermin1,2,H.Dole1,2,G.Lagache1,2,D.LeBorgne3,andA.Penin1,2 1UnivParis-Sud,LaboratoireIAS,UMR8617,Orsay,F-91405 e-mail:[email protected] 2CNRS,Orsay,F-91405 3Institutd’AstrophysiquedeParis(IAP),UMR7095CNRS,UPMC,98bisboulevardArago,F-75014Paris,France Submitted30September2010/Accepted21January2011 1 1 Abstract 0 2 Aims.Weaimatmodelingtheinfraredgalaxyevolutioninawayassimpleaspossibleandreproducestatisticalpropertiesamong n whichthenumbercountsbetween15µmand1.1mm,theluminosityfunctions,andtheredshiftdistributions.Wethenaimatusing a thismodeltointerprettherecentobservations(Spitzer,Akari,BLAST,LABOCA,AzTEC,SPTandHerschel),andmakepredictions J forPlanckandfutureexperimentslikeCCATorSPICA. 1 Methods. Thismodel usesanevolution indensity andluminosityof theluminosityfunctionparametrized bybroken power-laws 2 withtwobreaksatredshift∼0.9and2,andcontainsthetwopopulationsoftheLagacheetal.(2004) model:normalandstarburst galaxies.Wealsotakeintoaccounttheeffectofthestronglensingofhigh-redshiftsub-millimetergalaxies.Thiseffectissignificant ] inthesub-mmandmmrangenear50mJy.Ithas13freeparametersand8additionalcalibrationparameters.Wefittheparametersto O theIRAS,Spitzer,HerschelandAzTECmeasurementswithaMonte-CarloMarkovchain. C Results. The model ajusted on deep counts at key wavelengths reproduces the counts from the mid-infrared to the millimeter wavelengths,aswellasthemid-infraredluminosityfunctions.Wediscussthecontributiontothecosmicinfraredbackground(CIB) . h andtotheinfraredluminositydensityofthedifferentpopulations.Wealsoestimatetheeffectofthelensingonthenumbercounts, p anddiscuss therecent discoveryby theSouthPoleTelescope(SPT)of avery bright population lyingathigh-redshift. Wepredict - thecontributionofthelensedsourcestothePlancknumbercounts,theconfusionlevelforfuturemissionsusingaP(D)formalism, o and the Universe opacity to TeV photons due to the CIB. Material of the model (software, tables and predictions) is available at r http://www.ias.u-psud.fr/irgalaxies/. t s a Keywords.Cosmology:diffuseradiation-Galaxies:statistics-Galaxies:evolution-Galaxies:starformation-Infrared:galaxies- [ Submillimeter:galaxies 2 v 0 1. Introduction The physical models (such as Laceyetal. (2010); 5 Wilmanetal. (2010); Younger&Hopkins (2010) for the 1 The extragalactic backgroundlight (EBL) is the relic emission latest)useaphysicalapproachbasedonsemi-analyticalrecipes 1 due to galaxy formation and accretion processes since the anddarkmatternumericalsimulations.Theyusealimitedsetof . recombination.Theinfrared(8µm < λ < 1000µm)partofthis physicalparameters,buttheynowadayspoorlyreproducesome 0 1 emissioncalledcosmicinfraredbackground(CIB)wasdetected basic observationalconstraintslike the infraredgalaxy number 0 for the first time by Pugetetal. (1996) and containsabouthalf counts(Oliveretal.2010). 1 of the energy of the EBL (Doleetal. 2006; Be´therminetal. : 2010a). Nevertheless, in the local universe, the optical/UV v emissions are 3 times larger than infrared/sub-millimeter ones i The backwards evolution models (like Lagacheetal. X (Soifer&Neugebauer 1991; Driveretal. 2008). This pseudo- (2004); Franceschinietal. (2010); Rowan-Robinson (2009); paradoxis explainedby a strong evolution of the propertiesof r Valianteetal. (2009)) use an evolutionthe luminosity function a theinfraredgalaxies. (LF)ofthegalaxiestoreproduceempiricallythegalaxycounts, andotherconstraints.Thesemodelsmakeonlyadescriptionof The infrared luminosity density is dominated by theevolutionandcontainlittlephysics.Theparametersofthese normal galaxies (L < 1011L ) in the lo- models were tuned manually to fit observational constraints. IR,bolometric ⊙ cal Universe (Saundersetal. 1990). At higher redshift, LeBorgneetal. (2009) used an other approach and performed it is dominated by luminous infrared galaxies (LIRG, a non-parametric inversion of the counts to determine the 1011L < L < 1012L ) at z=1 (LeFloc’hetal. LF. Nevertheless, this approach is complex, uses only one ⊙ IR,bolometric ⊙ 2005) and by ultra-luminous infrared galaxies (ULIRG, population of galaxy, and does not manage to reproduce the 1012L < L < 1013L ) at z=2 (Caputietal. 2007). 160µm numbercounts.Anotherfully-empiricalapproachwas ⊙ IR,bolometric ⊙ The infrared luminosity of these galaxies is correlated to the used by Dominguezetal. (2010). They fitted the SED from star formationrate (Kennicutt1998). Thus modelingthis rapid UV to mid-infrared of detected galaxies and extrapolated the evolutionoftheinfraredgalaxiesisveryimportanttounderstand far-infrared spectral energy distribution of these galaxies and thehistoryofthestarformation. the contributionof faint populations.Nevertheless, their model aimsonlytoreproducetheCIB;howeveritsabilitytoreproduce 1 Be´therminetal.:Parametricbackwardsevolutionmodel otherconstraintslikethenumbercountswasnottested. on a very large library of templates and claim that the contri- bution of the AGNs and the dispersion of the dust temperature TheBalloon-borneLarge-ApertureSubmillimeterTelescope of the galaxies must be taken into account to reproduce the (BLAST) experiment (Pascaleetal. 2008; Devlinetal. 2009) observationalconstraint. Our approach is to keep the modelas and the spectral and photometric imaging receiver (SPIRE) simpleaspossible,buttouseadvancedmethodstoconstrainits instrument (Griffin 2010) onboard the Herchel space telescope free parameters. This new parametric modelcan be used as an (Pilbratt&al. 2010) performed recently new observations inputforhalomodelingorP(D)analysisforinstance. in the sub-mm at 250, 350 and 500 µm. In their current version, most of the models fail to reproduce the number Asitwillbeshown,wedidnotneedAGNcontributionand counts measured at these wavelengths (Patanchonetal. 2009; temperature dispersion to reproduce the current observational Be´therminetal. 2010b; Clementsetal. 2010; Oliveretal. constraints. In fact, in the local Universe, the AGNs only 2010). The Valianteetal. (2009) model gives the best results, dominate the ULIRG regime (Imanishi 2009). Alexanderetal. using a Monte Carlo approach (sources are randomly taken in (2005) estimate an AGN contribution of 8% for the submil- libraries) to simulate the temperature scatter and the hetero- limeter galaxies (SMG). Recently Faddaetal. (2010) showed geneityofthepopulationsofactivegalacticnucleus(AGN),but that the proportion of AGN-dominated sources is rather small this model strongly disagrees with the recent measurementsof for LIRGs at z∼1 (5%) and ULIRGs around z∼2 (12%). theredshiftdistributionoftheCIBbyJauzacetal.(2010).Itis Jauzacetal.(2010)showedthatAGNcontributiontotheCIBis thusnecessarytodevelopnewmodelsthatreproducetherecent lessthan10%atz<1.5.Thesecategoryofluminositydominates far-infraredandsub-mmobservations. the infrared output at their redshift. The low contribution of AGN in these categories explains why the AGNs are not The discoveryofverybrightandhigh-redshiftdustygalax- necessary to reproduce the mean statistical properties of the ies by Vieiraetal. (2009) with the south pole telescope (SPT) galaxies. Nevertheless, despite their small contribution to the suggeststhatthecontributionofhigh-redshiftgalaxiesstrongly infrared output, the AGNs play a central role in the physics of lensed by dark matter halos of massive low-redshift galaxies galaxies. on the bright sub-millimeter and millimeter counts is non negligible. This contribution was discussed by Negrelloetal. Ourmodeltakesintoaccountthethestrong-lensingofhigh (2007) and an observationalevidence of this phenomenonwas redshift galaxies by dark matter halos of elliptical galaxies. found very recently by Negrelloetal. (2010). We can also cite According to the results of Sect. 7.3, the effect of the lensing thesimplifiedapproachofLimaetal.(2010)whoreproducethe on the counts we fitted is smaller than 10%. The model of AzTEC and SPT counts using a single population of galaxies lensing does not have free parameters. It is based on WMAP- with a Schechter LF at a single redshift and a lensing model. 7-years-best-fit cosmology and on some parameters taken at We can also cite the very recent work of Hezaveh&Holder valuesgivenbythelitterature.Thelensingisthusnotusefulto (2010)ontheeffectofthelensingontheSPTcounts,basedon reproduce the current observations, but is necessary to make anadvancedlensingmodel. predictionsatbrightfluxes(>100mJy)inthesub-mmandmm range,wheretheeffectsofthelensingarelarge. We present a new simple and parametric model based on Lagacheetal. (2004) SED libraries, which reproduces the new observational constraints. The parameters of this model 3. Descriptionofthemodel (13 free parameters and 8 calibration parameters) were fitted from a large set of recent observations using a Monte-Carlo 3.1.Basicformulas Markovchain(MCMC)method,allowingtostudydegeneracies ThefluxdensityS atafrequencyνofasourcelyingataredshift between the parameters. This model also includes the effects ν zis(Hogg1999)is ofthe stronglensingon theobservations.We makepredictions on the confusion limit for future missions, on the high-energy S = (1+z)L(1+z)ν (1) opacityofthe Universeandontheeffectsofthestronglensing ν 4πD2(z) on the counts. This modelis plugged to a halo model to study L the spatial distribution of the infrared galaxies in a companion where z is the redshift, DL is the luminosity distance of the paper (Peninetal. 2010). Note that an other study using also source,andL(1+z)ν istheluminosityatafrequency(1+z)ν.The MCMCmethodswasperformedbyMarsdenetal.(2010)atthe comoving volume corresponding to a redshift slice between z sametimethanours. andz+dzandaunitsolidangleis dV (1+z)2D2 WeusetheWMAP7yearbest-fitΛCDMcosmologyinthis = D A (2) H paper(Larsonetal.2010).WethushaveH0=71km.s−1.Mpc−1, dz ΩΛ+(1+z)3Ωm ΩΛ=0.734andΩm =0.266. where D pistheHubbledistance(D = c/H ), D theangular H H 0 A distancetotheredshiftz.Ω andΩ arethenormalizedenergy m Λ densityofthematterandofthecosmologicalconstant. 2. Approach The backward evolution models are not built on physical 3.2.Bolometricluminosityfunctionanditsevolution parameters. Each model uses different evolving populations to reproduce the observational constraints. Some recent models Weassumethattheluminosityfunction(LF)isaclassicaldouble (like Franceschinietal. (2010); Rowan-Robinson (2009)) use exponentialfunction(Saundersetal.1990) 4 galaxy populations evolving separately to reproduce the L 1 L observations.Valianteetal.(2009)takerandomlygalaxySEDs Φ(LIR)=Φ⋆ × LI⋆R)1−α × exp −2σ2log210(1+ LI⋆R)) (3) (cid:16) h i 2 Be´therminetal.:Parametricbackwardsevolutionmodel where th is the hyperbolic tangent function. L is the lumi- pop nosityatwhichthenumberofnormalandstarburstgalaxiesare equal.σ characterizesthewidthofthetransitionbetweenthe pop twopopulations.AtL = L ,thefractionofstarburstis50%. IR pop There are 88%of starburstat LIR = Lpop ×10σpop, and 12%at LIR = Lpop ×10−σpop. The contributionof the differentpopula- tionstothelocalinfraredbolometricLFisshowninFig.1. 3.4.Observables Thenumbercountsatdifferentwavelengthsareanessentialcon- straintforourmodel.Thesourceextractionbiasesareingeneral accuratelycorrectedfortheseobservables.Thecountsarecom- putedwiththefollowingformula dN (S )= (5) dS dΩ ν ν ∞ dN dL dV f (L ) IR dz (6) Figure1. Solidline:Localinfraredbolometricluminosityfunc- XpopZ0 pop IR dLIRdV(cid:12)(cid:12)LIR(Sν,z,pop) dSν dzdΩ tionfromourbest-fitmodel.Reddottedline:contributionofthe (cid:12) ∞ dN = dz (7) normalgalaxies.Bluedashedline:contributionofthe starburst XpopZ0 dSνdzdΩ galaxies.Blackverticallongdashedline:luminosityofthetran- sitionbetweenthetwopopulation(Lpop). where dN/dS /dΩ is the number of source per flux unit and ν per solid angle. f (L ) is the fraction of the sources of a pop IR givengalaxypopulationcomputedwiththeEq.4.dN/dL /dV where Φ(L ) is the number of sources per logarithm of lumi- IR IR iscomputedfromtheEq.3 nosityandpercomovingvolumeunitforaninfraredbolometric luminosityL .Φ isthenormalizationconstantcharacterizing IR ⋆ dN dN Φ(L ) the density of sources. L⋆ is the characteristic luminosity at = = IR . (8) the break. 1−α and 1−α−1/σ2/ln2(10) are the slope of the dLIRdV dlog10(LIR)LIRlog(10)dV LIRlog(10) asymptotic power-law behavior at respectively low and high L (S ,z,pop) and dL /dS were computed on a grid in S luminosity. IR ν IR ν ν and z from the cosmologyand the SED templates.These grids donotdependontheevolutionoftheLFnoronthepopulation Weassumeacontinuousevolutioninluminosityandinden- mixing parameters. These grids are thus generated only once sity of the luminosity function with the redshift of the form andsavedtoacceleratethecomputationofthecounts.Notethat L⋆ ∝(1+z)rL andΦ⋆ ∝(1+z)rΦ whererLandrφarecoefficients with this method, it is very easy to change the SED templates drivingtheevolutioninluminosityanddensity,respectively.Itis and/oraddotherpopulations. impossibletoreproducetheevolutionoftheLFwithconstantr L andr .Weconsequentlyauthorizetheirvaluetochangeatsome φ Other measurements help to constraint our model. For ex- specificredshifts.Thepositionofthesebreaksarethesamefor ample,themonochromaticluminosityΦ functionatagiven both r and r . The position of the first redshiftbreak is a free mono L φ redshiftis parameterandconvergetothesamefinalvalueforinitialvalues between0and2.Toavoidadivergenceathighredshift,wealso dlog (L ) addasecondbreakfixedatz=2. Φmono = fpop(LIR(νLν))φ(LIR(νLν)) d(1ν0L )IR . (9) Xpop ν 3.3.Spectralenergydistribution(SED)ofthegalaxies We do not use the bolometric LFs, because they are biased by thechoiceoftheassumedSEDofthesources. We use the Lagacheetal. (2004) SED library. This library contains two populations: a starburst one and a normal one. WecanalsocomputetheredshiftdistributionN(z)forase- This library is parametrized only by the infrared bolometric lectioninfluxS >S with luminosity (L ). There is no evolution of the SED with the ν ν,cut IR redshift.Thenormalpopulationhasaspectrumtypicalofspiral ∞ dN galaxy.TheSEDofthispopulationdoesnotevolvewithLIR.On N(z,Scut)= dSν. (10) the contrary, the starburst SED evolves with LIR. The brighter ZSν,cut dSνdz thestarburstgalaxy,thehotterthedust. Theextragalacticbackgroundduetothegalaxiesatagiven Thenormalgalaxiesaredominantatlowluminosityandthe wavelengthis starburstathighluminosity.Wethuschosearbitrarythefollow- ∞ ∞ dN ∞ dN ingsmoothfunctiontodescribethefractionofstarburstgalaxies B = S dS dz= S dS (11) asafunctionofthebolometricluminosityLIR: ν Zz=0ZSν=0 νdSνdzdΩ ν ZSν=0 νdSνdΩ ν Φstarburst = 1+thhlog10(LIR/Lpop)/σpopi (4) andcanbecomparedtothemeasurementsoftheCIB. Φ 2 3 Be´therminetal.:Parametricbackwardsevolutionmodel Parameter Description Value α Faint-endslopeoftheinfraredbolometricLF 1.223±0.044 σ Parameterdrivingthebright-endslopeoftheLF 0.406±0.019 L (z=0)(×1010 L ) LocalcharacteristicluminosityoftheLF 2.377±0.363 ⋆ ⊙ φ (z=0)(×10−3gal/dex/Mpc3) LocalcharacteristicdensityoftheLF 3.234±0.266 ⋆ r Evolutionofthecharacteristicluminositybetween0andz 2.931±0.119 L⋆,lz break,1 r Evolutionofthecharacteristicdensitybetween0andz 0.774±0.196 phi⋆,lz break,1 z Redshiftofthefirstbreak 0.879±0.052 break,1 r Evolutionofthecharacteristicluminositybetweenz andz 4.737±0.301 L⋆,mz break,1 break,2 r Evolutionofthecharacteristicdensityofbetweenz andz -6.246±0.458 phi⋆,mz break,1 break,2 z Redshiftofthesecondbreak 2.000 (fixed) break,2 r Evolutionofthecharacteristicluminosityforz>z 0.145±0.460 L⋆,hz break,2 r Evolutionofthecharacteristicdensityforz>z -0.919±0.651 phi⋆,hz break,2 L (×1010 L ) Luminosityofthetransitionbetweennormalandstarbursttemplates 23.677±2.704 pop ⊙ σ Widthofthetransitionbetweennormalandstarbursttemplates 0.572±0.056 pop Table1. Parametersofourmodelfittedtoourselectionofinfraredobservations.TheerrorsarederivedfromtheMCMCanalysis. where c is the speed of light and σ the velocity dispersion in v thehalo,whichdependsonthecosmology,theredshiftandthe massofthehalo. TheprobabilityP(µ ,z )forasourceataredshiftz tobe min s s magnifiedbyafactorgreaterthanµ is min P(µ>µ ,z )= (15) min s (1+zs)2 zs ∞ dN σ(µ>µ )dVdMdz (16) 4πDc(zs)Z0 Z0 d(log10(M))dV min dz where z is the redshift of the source, D the comoving radial s c distance, dN is the halo mass function, and dV is the d(log10(M))dV dz comovingvolumeassociatedtotheredshiftslicedz.Weusethe halomassfunctionofReedetal.(2007). Figure2. Effect of the lensing on the number counts at 850 Thecountsderivedby ourmodeltake into accountthe fact microns. The contribution of lensed source is multiply by 10 that a small fraction of the sources are gravitationally magni- to underlinethe effect of the lensing on the counts. Red dotted fied.Theobservednumbercountstakingintoaccountthelensing line:countsofnon-lensedsources.Greendashedline:countsof (dN/dSν/dΩ)lensedarecomputedfrominitialcountsdSν/dz/dΩ lensedsources.Blacksolidline:totalcounts. with dN (S )= (17) The level of the non-correlated fluctuations (shot-noise) of dS dΩ lensed ν (cid:16) ν (cid:17) tthioenC: IBcanbeeasilycomputedfromourmodelwiththeequa- ∞ µmax dP(z)1 dN Sν,z dµdz. (18) Sν,cut dN Z0 Zµmin dµ µdSνdzdΩ µ ! P = S2 dS (12) Practically, this operation is performed multiplying the vector SN Z0 νdSνdΩ ν containingthe countsfor a givenredshift slice by a matrix de- where P is the level of the non-correlated fluctuations and scribing the effect of lensing. This lensing matrix has diagonal SN S thefluxlimitforthecleaningoftheresolvedsources. coefficients values around 1, and small (< 10−3) non-diagonal ν,cut terms. These non-diagonalterms describe how magnified faint sources affect the counts at brighter flux. The effect on the monochromaticluminosityfunctionwas computedin the same 3.5.Effectofthestronglensingonthecounts way. We chose µ = 2 which corresponds to the limit of min We use a simple strong lensing model based on Perrottaetal. the validity of the strong-lensing hypothesis (Perrottaetal. (2001,2002).Itsupposesthatthedarkmatterhalosaresingular 2001). The spatial extension of the lensed galaxies limits the isothermalspheres.Thecross-sectionσofahaloforamagnifi- maximum magnification. According to Perrottaetal. (2002), cationµlargerthanµmin is µmax is in the 10-30 range. We chose to use µmax=20 in this paper.Negrelloetal.(2007)usedµ =2andµ =15. min max 4πα2D σ(µ>µ )= A,ls (13) min µ2 Fig. 2 illustrates how number counts are affected by lens- ing.Thisfigureisbasedonthenumbercountspredictedbythe where D is the angular-diameter distance between the lens A,ls modelat850µmwithaprobabilityofmagnificationmultiplied andthesourceandαisgivenby by a factor 10to bettershow this effect. The greendashedline σ2 is contributionofthe lensed sources.Due to the magnification, α=4π v (14) thepeakofthiscontributionisathigherfluxthanfornon-lensed c2 4 Be´therminetal.:Parametricbackwardsevolutionmodel sources,andduetothesmallprobabilityoflensing,thepeakis error bars, which could be largely underestimated. For in- lowerthanfornon-lensedsources.Thiseffectofthemagnifica- stance, Be´therminetal. (2010a) estimate that the Poissonian tiononthecountsbecomenonnegligiblewhentheslopeofthe uncertainties underestimate the real sample uncertainties by countsisverysteep,likeinthesub-mmandmmdomain. afactor3forcountsaround100mJyat160µmina10deg2field. We do not fit the 850 µm observation because of the large 4. Fittingthemodelparametersonthedata discrepanciesbetweenthesubmillimetercommon-userbolome- ter array (SCUBA) observations (Coppinetal. 2006) and the Ourmodelhasseveralfreeparameters.Wetriedtohavethemin- large APEX bolometer Camera (LABOCA) ones (Weißetal. imumnumberofparameters.Wedeterminedthembyfittingthe 2009).WediscussthisproblemintheSect.5.4. model to published measurements of the counts and LFs. We usedaMonte-CarloMarkovchain(MCMC)tofindthebestpa- rameters,theiruncertainties,anddegeneracies.Wedonotfitthe We fitted the AzTEC measurements at 1.1 mm of measuredredshiftdistributions,becausethecosmicvarianceand Austermannetal. (2010) and Scottetal. (2010). The area theselectioneffectsarecurrentlynotenoughaccuratelyknown. coveredbyAzTECis smallcomparedto Spitzer andHerschel. We used two independentmeasurementsof the AzTEC counts toincreasetheweightofmmobservationsinourfit. 4.1.Data:extragalacticnumbercounts 4.1.1. Datausedforthefit Wehavechosentofitthefollowingdata: 4.2.Data:monochromaticluminosityfunctions – Spitzer MIPS counts of Be´therminetal. (2010a) at 24, 70 4.2.1. Datausedforthefit and160µm, Wehavechosentofitthefollowingdata: – HerschelSPIREOliveretal.(2010)countsat250,350and 500µm, – IRAS localluminosity functionat 60 µm of Saundersetal. – AzTEC counts of Austermannetal. (2010) and Scottetal. (1990) (2010)at1.1mm. – Spitzer local luminosity function at 24 µm of Rodighieroetal.(2009), 4.1.2. Justificationofourchoice – Spitzer luminosity function at 15 µm at z=0.6 of Rodighieroetal.(2009), We fit only the differential number counts since the integral – Spitzer luminosity function at 12 µm at z=1 of countsarehighlycorrelatedandthecorrelationmatrixisrarely Rodighieroetal.(2009), estimated. – Spitzer luminosity function at 8 µm at z=2 of Caputietal. (2007), The number counts were measured at numerous bands between 15 µm and 1.1 mm. We have chosen a collection of points. We were guided by the reliability of the measurements 4.2.2. Justificationofourchoice andtheirerrorbars. Wefittedsomemonochromaticluminosityfunctions.Wechose Numbercountsat15µmbasedontheinfraredspaceobser- only wavelengthsand redshifts for which no K-correctionsare vatory(ISO)data(Elbazetal.1999;Gruppionietal.2002)and needed. These observations strongly constrain the parameters on the Akari data (Pearsonetal. 2010; Hopwoodetal. 2010) drivingtheredshiftevolutionsofourmodel. exhibitadiscrepancybyafactorofabout2,andtheirerrorsdo not include cosmic variance. The results of these papers were From the Rodighieroetal. (2009) LFs measured with the notfitted.Nevertheless,wecomparedaposterioritoourresults Spitzer data at 24 µm, we computed 3 non K-corrected LFs at tocheckconsistencyinSect.5.4. z=0,0.6and1.WeusedtheirlocalLFat24µm.Atz=0.6and 1,insteadofusingdirectlytheirresultsintheirredshiftbins,we We fitted the Spitzer MIPS counts of Be´therminetal. combined their 15 µm LF at z=0.6 (respectively 12 µm LF at (2010a)at24,70and160µm.Thesepointswerebuiltfromthe z=1)inthe0.45 < z < 0.6and0.6 < z < 0.8bins(respectively data of FIDEL, COSMOS and SWIRE legacy programs. The 0.8 < z < 1.0and1.0 < z < 1.4)to obtain15µm LFatz=0.6 errorsbarstakeintoaccountthecosmicvariance.Thesecounts (respectively 12 µm LF at z=1). The error on a point is the agree with previous Spitzer measurements of Papovichetal. maximumofthecombinationofthestatisticalerrorsofthetwo (2004);Shupeetal.(2008);LeFloc’hetal.(2009);Frayeretal. binsandofthedifferencebetweenthemeasuresinthetwobins. (2009) and Herschel measurements of Bertaetal. (2010) (in The second value is often larger due to the quick evolution of whichthedifferentfieldsarenotcombined). the LF and the cosmic variance. We fitted only the points that do not suffer incompletenessto avoid possible biases. We also At250µm,350µmand500µm,wefittedHerschelSPIRE fittedthe8µmatz=2ofCaputietal.(2007). Oliveretal. (2010) counts which take into account the cosmic variance and the deboosting uncertainty. These counts agree We also fitted the local LF at 60 µm determined from the with the BLAST measurements of Patanchonetal. (2009) infrared astronomical satellite (IRAS) data (Saundersetal. and Be´therminetal. (2010b) and the Herschel measurements 1990) to better constrain the faint-end slope of the local LF. of Clementsetal. (2010). We chosen Oliveretal. (2010) Due to the strong AGN contaminationat 60 µm in the ULIRG counts because Herschel data are better than BLAST ones regime,wedidnotfitthepointsbrighterthan1011.5L at60µm. ⊙ and because Clementsetal. (2010) countsuse only Poissonian 5 Be´therminetal.:Parametricbackwardsevolutionmodel Instrument Calibrationparameter(γ ) Calib.uncertainty step n is kept with a probability of 1 if L(θ ) > L(θ ) or b n n−1 MIPS24µm 1.00±0.03 4% with a probability L(θ )/L(θ ) else. The distribution of the n n−1 MIPS70µm 1.06±0.04 7% realizationofthechainisasymptoticallythesameastheunder- MIPS160µm 0.96±0.03 12% lyingprobabilitydensity.Thispropertyis thusveryconvenient SPIRE250µm 0.88±0.05 15% todeterminetheconfidenceareafortheparametersofthemodel. SPIRE350µm 0.97±0.07 15% SPIRE500µm 1.17±0.1 15% We used the Fisher matrixformalism to determine the pro- AzTEC1.1mm 0.98±0.09 9% posaldensityofthechainfrominitialparametersvaluessetman- Table 2. Calibration parameters and 1-σ marginalized errors ually.TheassociatedFishermatrixis: from our MCMC fit compared with calibration uncertainties givenbytheinstrumentalteams. Npoints ∂m ∂m 1 1 F (θ)= model,k model,k + (20) ij ∂θ ∂θ 2σ2  2σ2  Xk=1 i j m  calib,b where θ is a vector containing the model and calibration (γ ) 4.3.Data:CIB b parameters.Theterm in bracketsappearsonly forthe diagonal The bulk of the CIB is not resolved at SPIRE wavelengths. terms corresponding to a calibration parameter. We ran a first We thus used the absolute measurement of the CIB level shortchain(10000steps)andcomputedanewproposaldensity in SPIRE bands as a constraint for our model. We used with the covariance matrix of the results. We then ran a sec- the Lagacheetal. (1999) measurement performed on the far- ond long chain of 300 000 steps. The final chain satisfies the infrared absolute spectrophotometer (FIRAS) data: 11.7±2.9 Dunkleyetal.(2005)criteria(j⋆ >20andr <0.01). nW.m2.sr−1 at 250 µm, 6.4±1.6 nW.m2.sr−1 at 350 µm and 2.7±0.7nW.m2.sr−1 at500µm.WeassumethattheCIBisonly duetogalaxiesandthusneglectapossibleextragalacticdiffuse 5. Resultsofthefit emission. 5.1.Qualityofthefit Our final best-fit model have a χ2 (χ2 = −2log(L) because all 4.4.Calibrationuncertainties errors are assumed to be Gaussian) of 177 for 113 degrees of freedom.Ourfitisthusreasonablygood.Theparametersfound The calibration uncertainty is responsible for correlated un- withthefitaregiveninTable1(theuncertaintiesarecomputed certainties between points measured at a given wavelength fromtheMCMC).Thecalibrationfactorarecompatiblewiththe with thesame instrument.A changein thecalibrationmodifies calibrationuncertaintiesgivenbytheinstrumentalteamswitha globally the number counts and the LF. Assuming the ”good” χ2 of 2.89for 7 points(see Table 2). The resultsare plotted in calibration is obtained in multiplying the fluxes by a factor γ, the”good”normalizedcountsareobtainedwithS = γS and Fig.3. new (S2.5 dN/dS ) = γ1.5(S2.5dN/dS). The effect on the LF in good good dex per volume unit is more simple. We just have to shift the 5.2.Comparisonbetweenthemodelandtheobserved luminosityinabscissabyafactorγ. countsusedinthefit TheBe´therminetal.(2010a)pointsfitgloballywell,withsome We added to our free parameters a calibration parameter exceptions. Our model is lower by about 15% than two points for each fitted band (see Table 2). We took into account the around300µJyat24µm.Thesetwopointsarebuiltcombining uncertainties on the calibration estimated by the instrumen- theFIDEL,COSMOSandSWIREfields.TheSWIREfieldsare tal team in our fit (Stansberryetal. 2007; Gordonetal. 2007; shallowfieldsandthecountscouldbeaffectedbytheEddington Engelbrachtetal.2007;Swinyardetal.2010;Scottetal.2010). bias. We also observe a slight under-prediction of the bright (S >50mJy)countsat70µm.WealsoplottedtheBertaetal. 70 4.5.Fittingmethod (2010) counts at 160 µm measured using the photodetector arraycameraandspectrometer(PACS)ontheHerschelsatellite. Tofitourpoints,weassumedthattheuncertaintiesonthemea- ThesecountsagreewithBe´therminetal.(2010a)andourmodel. surements and on the calibrations are Gaussian and not corre- lated.Thelog-likelihoodisthen Our model fits globally well the Oliveretal. (2010) and Be´therminetal. (2010b) counts, excepting a slight under- Npoints (m −m (θ))2 Nband (γ −1)2 −log(L(θ))= k model,k + b (19) prediction of the counts between 30 mJy and 100 mJy at 2σ2 2σ2 500 µm. There is a mild disagreementwith the Clementsetal. Xk=1 m Xb=1 calib,b (2010)counts,buttheirerrorsbarsdonottakeintoaccountthe where L is the likelihood, θ the parameters of the model, m a cosmic variance and are thus underestimated. We also plotted k measurement,m thepredictionofthe modelforthe same the results of the P(D) analysis of Glennetal. (2010). These model,k measurement, σ the measurement uncertainty on it, γ the points and especially the error bars must be interpreted with m b calibration parameter of the band b and σ the calibration caution(seethecompletediscussioninGlennetal.(2010)).We calib,b uncertaintyforthisband. have plotted the knots of the smooth and power-law models. Theygloballyagreewithourmodel. WeusedaMonteCarloMarkovchain(MCMC)Metropolis- Hastings algorithm (Chib&Greenberg 1995; Dunkleyetal. Our model agrees very well with the AzTEC counts of 2005) to fit our model. The method consists in a random walk Austermannetal. (2010) and Scottetal. (2010). The contribu- in the parameter space. At each step, a random shift of the tionofthe stronglensingobjectstothe AzTECcountsisweak parameters is done using a given fixed proposal density. A (<10%,seeSect.7.3). 6 Be´therminetal.:Parametricbackwardsevolutionmodel Figure3. (a)to(f):Differentialextragalacticnumbercountsusedforthefit.(h):MonochromaticLFsatdifferentwavelengthsand redshifts.(a)to(h):Thefittedpointsarethicker.Blacksolidline:ourbest-fitmodel.Blackdashedline:1-σrangeofthemodel.(a) to(c):Reddiamonds:Be´therminetal.(2010a)Spitzerlegacynumbercounts.(c):Greentriangles:Bertaetal.(2010)Herschel/PEP numbercounts.(d)to(f):Reddiamonds:Oliveretal.(2010)Herschel/Hermesnumbercounts.Greentriangles:Glennetal.(2010) Herschel/Hermes P(D) analysis. Clementsetal. (2010) Herschel/ATLAS number counts. Purple cross: Be´therminetal. (2010b) BLAST number counts. (g): Green triangles: Scottetal. (2010) AzTEC number counts in the CDFS field. Green triangles: Austermannetal. (2010) AzTEC number counts in the SHADES field. (h): Red plus: Rodighieroetal. (2009) local 24 µm LF (notfittedpointsingrey).Greendiamonds:Saundersetal.(1990)local60µmLF(shiftedbya factor10onthey-axis;notfitted pointsingrey);Bluetriangles:Rodighieroetal.(2009)15µmLFatz=0.6(shiftedbyafactor100onthey-axis;notfittedpointsin grey).Purplesquares:Rodighieroetal.(2009)12µmLFatz=1(shiftedbyafactor1000onthey-axis;notfittedpointsingrey). Cyancrosses:Caputietal.(2007)8µmLFatz=2(shiftedbyafactor10000onthey-axis). 7 Be´therminetal.:Parametricbackwardsevolutionmodel 5.3.Comparisonbetweenthemodelandtheobserved monochromaticLFs OurmodelfitswellourcollectionofLFs(Saundersetal.1990; Caputietal.2007;Rodighieroetal.2009),exceptingthebright- est point of Caputietal. (2007). In Fig. 3, we arbitrary shifted thedifferentLFsonthey-axistoobtainaclearerplot.Ourmodel underestimatesthe 60 µm local LF in the ULIRG regime. It is expectedbecauseourmodeldonotcontainAGNsandconfirms ourchoiceofnotfittingthesepoints(Sect.4.2). 5.4.Comparisonbetweenthemodelandtheobserved countsnotusedinthefit We also compared our results with the counts at other wave- lengths. They are plotted on Fig. 4 and 5. The 1-σ region of the model includes the γ uncertainty of Akari at 15 µm (4%, b Ishiharaetal.(2010)),PACSat110µm(about10%,Bertaetal. (2010)) and LABOCA at 850 µm (8.5%, Weißetal. (2009)). The uncertainty on γ is about the same for LABOCA and b SCUBA (∼10%, Scottetal. (2006)). The uncertainties on the model are larger at these non-fitted wavelengths because the correlations between the model and the calibration parameters arenottakenintoaccountbythefit. At 15 µm, the Elbazetal. (1999) counts from different fields are not compatible between them, but our counts pass in the cloud of points. The Gruppionietal. (2002) counts are significantly lower than our model and other works. We marginally agree with the Pearsonetal. (2010) counts. The Hopwoodetal. (2010) counts measured with Akari in a field around Abell 2218 are lower than our model by about 25%. Nevertheless,theirfieldisverynarrowandtheirestimationmay sufferfromcosmicvariance.Finally,wewellagreewiththevery recent Teplitzetal. (2010) measurements performed with the infraredspectrograph(IRS)onboardtheSpitzerspacetelescope. We compare our counts to Hacking&Houck (1987), Lonsdaleetal. (1990), Rowan-Robinsonetal. (1990), Saundersetal. (1990), Gregorichetal. (1995) and Bertinetal. (1997) at 60 µm from IRAS data. There are disagreements betweenthe differentobservationsandsome errorbarsmay be underestimated,butourmodelgloballyagreeswiththecloudof points. Figure4. (a) to (c): Differential extragalactic number counts We can also compare the prediction of our model with not used for the fit. Black solid line: our best-fit model. Black Bertaetal.(2010)countsat110µm.Ourmodelgloballyagrees dashed line: 1-σ range of the model. (a): Red diamonds: withtheirwork.Nevertheless,ourmodeltendstobehigherthan Elbazetal.(1999)ISOcounts.Greentriangle:Gruppionietal. theirmeasurementnear100mJy.Observationsonseverallarger (2002) ISO counts. Blue squares: Pearsonetal. (2010) Akari fieldswillhelptoseeifthiseffectisanartifactornot. counts. Purple cross: Hopwoodetal. (2010) Akari (lensed) counts. Cyan plus: Teplitzetal. (2010) Spitzer/IRS counts. At 850 µm, we very well agree with the P(D) analysis (b): Red diamonds: Hacking&Houck (1987), Lonsdaleetal. of the LABOCA data of Weißetal. (2009) (see Fig. 5). (1990), Rowan-Robinsonetal. (1990), Saundersetal. (1990), But, the measurements performed with SCUBA (Borysetal. Gregorichetal.(1995)andBertinetal.(1997)IRAScounts.(c): 2003; Scottetal. 2006; Coppinetal. 2006) and LABOCA Reddiamonds:Bertaetal.(2010)Herschel/PEPcounts. (Beelenetal.2008)aresignificantlyhigherthanourmodelat6 and 8 mJy. At low flux (<2 mJy), our model agrees very well with the measurementperformed in lensed region (Smailetal. 2002;Knudsenetal.2008;Zemcovetal.2010). sources can be separated from dusty galaxies consideringtheir spectrum.Wethuscompareourresultswiththecountsofdusty We alsocompareourmodelpredictionswithSPTmeasure- sources. Vieiraetal. (2009) measured counts for all the dusty ments at 1.38 mm (Vieiraetal. 2009). At this wavelength, the sources and for the dusty sources without IRAS 60 µm coun- contribution of the synchrotron emission of the local radio- terpart.Ourmodelagreeswith thesetwo measurements.Fig. 6 galaxies to the counts is not negligible. Nevertheless, these showsthecountsofthenon-IRASdustysources.The7.2%cali- 8 Be´therminetal.:Parametricbackwardsevolutionmodel Figure5. Integralnumber counts at 850 µm. Black solid line: our best-fit model. Black dashed line: 1-σ range of the model. Grey plus: Zemcovetal. (2010) combined SCUBA lensed counts. Blue dashed line: Weißetal. (2009) LABOCA P(D) (Schechtermodel).Reddiamonds:Coppinetal.(2006)SCUBA SHADES counts. Cyan square: (Beelenetal. 2008) LABOCA counts around the J2142-4423 Lyα protocluster. Black plus: Knudsenetal. (2008) combined SCUBA lensed counts; Green triangles: Scottetal. (2006) revisited SCUBA counts. Purple cross: Borysetal. (2003) SCUBA HDFN counts.Yellow aster- isks:Smailetal.(2002)lensedcounts. Figure6. Number counts at 1.38 mm of dusty sources with- out IRAS 60 µm counterpart . Black diamonds: Vieiraetal. (2009) south pole telescope (SPT) measurements. Black solid line: Total contribution of S < 0.2 Jy sources. Green dot- 60 dashedline:Contributionofthenon-lensedsources.Reddashed line: Contribution of the strongly-lensed sources. Dotted lines 1-σcontours. brationuncertaintyofSPTistakenintoaccountinthe1-σregion ofthemodel. 5.5.Comparisonwiththeobservedredshiftdistributions Figure7. Redshift distribution of the S >80 µJy (a), 24 In Fig. 7, we compare our model predictions with observed S >300 µJy (b), S >40 mJy (c), S >30 mJy (d), 24 250 350 redshift distributions. At 24 µm, our model over-predicts S >20 mJy (e), S >3 mJy (f), S >3 mJy (g) sources. 500 850 1100 by about 20% the number of sources below z=1 compared These measurements are not fitted. Black solid line: our best- to LeFloc’hetal. (2009) observations for the selection fit model. Black dotted line: 1-σ range of the model. Grey S24 > 80µJy. Nevertheless, they exclude i+AB<20 galaxies and solid line: our best-fit model convolved by a gaussian of σz = their numberof sources at low redshift is thus underestimated. 0.125z. Purple three dot-dashed line: LeBorgneetal. (2009) model.Greendashedline:Valianteetal.(2009)model.Redas- terisks:LeFloc’hetal.(2009)(a,b),Chapinetal.(2010)(c,d, 9 e), Chapmanetal. (2005) (f) andChapinetal. (2009) (g)mea- surements. Blue diamonds: Rodighieroetal. (2009) measure- ments(a,b).Cyansquares:Desaietal.(2008)measurements(b). Be´therminetal.:Parametricbackwardsevolutionmodel Our model also underpredicts the number of sources at z>3. But, the redshifts of the z>2 sources are only moderately accurate (σ ≈ 0.25 for i+ >25 at z∼2). Because of the z AB strong slope of the redshift distribution, a significant number of sourcesmeasurednear z=3.5 could be sources lying around z=3 with overestimated redshift. If we convolve our model with a gaussianerrorwithσ = 0.125zto simulatethe redshift z uncertainties,the modeland the measurementsagrees (Fig. 7). The LeBorgneetal. (2009) model agrees very well with the measurements,exceptingatz<0.5andz>2.5.TheValianteetal. (2009) model poorly reproduces this observation. The same observables was measured by Rodighieroetal. (2009). Their resultsarein agreementwith LeFloc’hetal. (2009), excepting atz>3,wheretheyarehigher.Itcouldbeexplainbyalargerσ z athighredshift. We also compare the model with the redshift distribution of S > 300µJy sources measured by LeFloc’hetal. (2009), 24 Rodighieroetal.(2009)andDesaietal.(2008).Thesedifferent measurements exhibit disagreements below z=0.5. This differ- encecouldbeexplainedbytheremovingofthebrightestoptical sources(see previousparagraph).Our modeloverestimatesthe number of sources at z<0.5 by a factor 2. There is a rather good agreements between the models and the measurements between z=0.5 and z=2.5, except a small overestimation by Valianteetal. (2009) near z=2. At higher redshifts, the mea- surements are significantly higher than the models. It could be explainedbytworeasons:aneffectoftheredshiftuncertainties andtheabsenceofAGNinourmodel. We comparewith the Chapinetal. (2010) redshiftdistribu- tions of the BLAST isolated sources at 250 µm, 350 µm and Figure8. DifferentialcontributionoftheS > 80µJysources 500 µm. This selection of isolated sources does not allow to 24 to the CIB as a function of the redshift at 70 µm (upper know the effective size of the field. We thus normalized our panel) and 160 µm (lower panel). Red asterisks: measurement model and the measured counts to have dN/dzdz = 1. Our by stacking in the COSMOS field (Jauzacetal. 2010). Black predicted redshift distribution globally fitRs the measurements, solid line: Our model (1-σ limit in black dotted line). Purple except at low z at 250 µm and 350 µm. This difference could three dot-dashed line: LeBorgneetal. (2009) model. Green be explained by the selection of isolated sources, which could dashedline:Valianteetal.(2009)model.Bluedot-dashedline: miss sources in structures at low redshift. The other models Franceschinietal.(2010)model. (LeBorgneetal. 2009; Valianteetal. 2009) underpredicts the number of sources at low z. Valianteetal. (2009) also slightly overpredictsthenumberofsourcesatz∼1.5. cept near z=0.5 (see Fig. 8), where their low data pointscould We compared the redshift distribution of the SCUBA come from a large-scale underdensityin the COSMOS field at sources at 850 µm with the prediction of our model. We use thisredshift.TheLeBorgneetal.(2009)modeloverpredictsthe theselection-correctedmeasurementsofChapmanetal. (2005) contribution of the 24 µm sources at z>1. The Valianteetal. used by Marsdenetal. (2010). All the models agrees with this (2009) model does not reproduce the trend of these measure- measurement. ments.Franceschinietal.(2010)underestimatethecontribution of the local sources and overestimate the contribution of z∼1 We also compared the prediction of our model with the sources. redshift distribution of the sources detected at 1.1 mm by AzTEC (Chapinetal. 2009). A significant part of the sources 5.6.ComparisonwiththemeasuredPoissonfluctuationsof detected atthis wavelength(10 over28) are notidentified,and theCIB the selection is not performed in flux, but in signal-to-noise ratio. Consequently, the normalization of the redshift distribu- Table 3 summarizes the recent measurements of the non- tion is not known. We thus use the same normalization than correlated fluctuations of the CIB (P ) and the predictions SN fortheBLASTredshiftdistributions( dN/dzdz = 1).Thebe- of our model. Note that P depends strongly on the S , SN ν,cut haviorpredictedbyourmodelagreeswRellwiththeobservations. the flux density at which the resolved sources are cleaned. We agreewiththemeasurementsofMiville-Descheˆnesetal.(2002) Recently,Jauzacetal.(2010)hasmeasuredthecontribution at 60 µm and 100 µm, Lagacheetal. (2007) at 160 µm and oftheS > 80µJy totheCIBat70and160µm asa function Vieroetal. (2009) at 250 µm and 350 µm. We found a value 24 oftheredshift.Theirstackinganalysisallowstocheckthetotal 35%lowerthanVieroetal.(2009)at500µm.Thisisconsistent far-infrared(FIR)emissionsofthefaintsourcesnotresolvedat with the slight under-estimation of the counts at 500 µm by thesewavelengths.Ourmodelagreeswellwiththeirresults,ex- our model. Our model is also about 40% lower than the SPT 10