The contribution of very massive high-redshift SWIRE galaxies to the stellar mass function
S. Berta, C.J. Lonsdale, M. Polletta, R.S. Savage, A. Franceschini, H. Buttery, A. Cimatti, J. Dias, C. Feruglio, F. Fiore, E.V. Held, F. La Franca, R. Maiolino, A. Marconi, I. Matute, S.J. Oliver, E. Ricciardelli, S. Rubele, N. Sacchi, D. Shupe, J. Surace
aa r X i v : . [ a s t r o - ph ] O c t Astronomy & Astrophysics manuscript no. 7491 c (cid:13)
ESO 2018October 24, 2018
The contribution of very massive high-redshift SWIRE galaxiesto the stellar mass function
Stefano Berta , , ,⋆ , Carol J. Lonsdale , , Maria Polletta , , Richard S. Savage , Alberto Franceschini ,Helen Buttery , Andrea Cimatti , Joao Dias , Chiara Feruglio , Fabrizio Fiore , Enrico V. Held , FabioLa Franca , Roberto Maiolino , Alessandro Marconi , Israel Matute , Seb J. Oliver , Elena Ricciardelli ,Stefano Rubele , Nicola Sacchi , David Shupe , and Jason Surace Dipartimento di Astronomia, Universit`a di Padova, Vicolo dell’Osservatorio 3, 35122 Padova, Italy. Center for Astrophysics and Space Sciences, University of California, San Diego, 9500 Gilman Dr., La Jolla, CA 92093-0424,USA. Max-Planck-Institut f¨ur Extraterrestrische Physik (MPE), Postfach 1312, 85741 Garching, Germany. Infrared Processing & Analysis Center, California Institute of Technology 100-22, Pasadena, CA 91125, USA. Institut d’Astrophysique de Paris, 98bis bld Arago, 75014 Paris, France. Astronomy Centre, CPES, University of Sussex, Falmer, Brighton BN19QJ, UK. INAF – Osservatorio Astronomico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy. Dipartimento di Astronomia, Universit`a di Bologna, Via Ranzani 1, 40127 Bologna, Italy. INAF – Osservatorio Astronomico di Roma, Via Frascati 33, I-00044 Monteporzio Catone, Italy. INAF – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy. Dipartimento di Fisica, Universit`a degli Studi ‘Roma Tre’, Via della Vasca Navale 84, I-00146 Roma, Italy. Dipartimento di Astronomia e Scienza dello Spazio, Universit`a di Firenze, Largo E. Fermi 2, 50125 Firenze, Italy. Spitzer Science Center, California Institute for Technology, 220-6, Pasadena, CA 91125, USA.Received 16 March 2007; accepted 24 September 2007
ABSTRACT
Context.
In the last couple of years a population of very massive ( M ⋆ > M ⊙ ), high-redshift ( z ≥
2) galaxies has beenidentified, but its role in galaxy evolution has not yet been fully understood.
Aims.
It is necessary to perform a systematic study of high-redshift massive galaxies, in order to determine the shape of thevery massive tail of the stellar mass function and determine the epoch of their assembly.
Methods.
We selected high- z massive galaxies at 5.8 µ m, in the SWIRE ELAIS-S1 field (1 deg ). Galaxies with the 1.6 µ m stellarpeak redshifted into the IRAC bands ( z ≃ −
3, called “IR-peakers”) were identified. Stellar masses were derived by means ofspectro-photometric fitting and used to compute the stellar mass function (MF) at z = 1 − −
3. A parametric fit to theMF was performed, based on a Bayesian formalism, and the stellar mass density of massive galaxies above z = 2 determined. Results.
We present the first systematic study of the very-massive tail of the galaxy stellar mass function at high redshift. Atotal of 326 sources were selected. The majority of these galaxies have stellar masses in excess of 10 M ⊙ and lie at z > . M ⋆ /L ratio. The influence of near-IR data and of the chosen stellar library on the SED fitting are alsodiscussed. The z = 2 − and ∼ M ⊙ is probed with unprecedented detail. A significantevolution is found not only for galaxies with M ∼ M ⊙ , but also in the highest mass bins considered. The comovingnumber density of these galaxies was lower by more than a factor of 10 at z = 2 −
3, with respect to the local estimate.SWIRE 5.8 µ m peakers more massive than 1 . × M ⊙ provide 30 −
50% of the total stellar mass density in galaxies at z = 2 − Key words.
Galaxies: evolution - Galaxies: mass function - Galaxies high-redshift - Galaxies: fundamental parameters - Galaxies:statistics - Infrared: galaxies
Send offprint requests to : Stefano Berta, e-mail: [email protected], ste [email protected] ⋆ SB was supported by the Ing. Aldo Gini Foundation
1. Introduction
Tracing the formation of galaxies and understanding theepoch when the bulk of their baryonic mass was assembledrepresents one of the major problems of modern cosmol-
Berta S., et al.: IR-peakers mass function ogy, particularly controversial when dealing with massive( M stars > M ⊙ ) objects.The assembly of massive galaxies is one of the crit-ical questions in the cosmic evolutionary scenario. Theuniform properties of local early-type galaxies and of thefundamental plane have inspired the so called “monolithiccollapse” scenario (Eggen et al. 1962; Chiosi & Carraro2002), in which galaxies formed in the remote past throughhuge events of star formation and subsequently evolvedpassively across cosmic time. On the other hand, inthe more recent “hierarchical” scenario (White & Rees1978; Kauffmann et al. 1996; Kauffmann & Charlot 1998;Somerville & Primack 1999), massive galaxies assemble bymergers of lower-mass units, with the most massive ob-jects being born in the latest stages of evolution, at z ≤ z ≃ −
2, with a subsequent decline of at leastone order of magnitude to the present time (e.g. Hopkins2004; Rudnick et al. 2003; Flores et al. 1999; Madau et al.1996, among others).An alternative approach consists of studying the massalready assembled in galaxies, instead of the amount ofstars being formed. The integral of the past star formationdensity provides the stellar mass density at a given epoch:a complementary constraint on cosmic galaxy evolution.Thus, the build up of the stellar mass across cosmictime has become one of the major topics in observa-tional cosmology, and has overtaken the classic Madau-Lilly diagram as the central tool for studying galaxy evo-lution. The large observational effort dedicated to thissubject has shown that the global stellar mass densityincreases from early epochs to the low-redshift Universe(e.g. Brinchmann & Ellis 2000; Dickinson et al. 2003a;Fontana et al. 2003, 2004, 2006; Rudnick et al. 2003, 2006;Drory et al. 2005). Very deep surveys have been ex-ploited to describe the shape of the stellar mass functionat high redshift (Fontana et al. 2006; Drory et al. 2005;Gwyn & Hartwick 2005), but a clear picture on the roleof very massive galaxies has not yet emerged.Several pieces of evidence exist that fully formed mas-sive galaxies were already in place at redshift z ∼ − z > M = 1 − × M ⊙ ), evolved (ages of1 − . z >
4. Based on FIRESdata, Rudnick et al. (2003) inferred that DRGs contribute ∼
50% of the global stellar mass density at z = 2 −
3. Dickinson et al. (2003a) exploited Hubble Deep FieldNorth NICMOS data to derive the stellar masses of galax-ies up to z = 3. The study of the global stellar mass den-sity highlighted that 50 −
75% of the present-day stellarmass was already in place at z ≃
1, while only 3 − z ≃ . z ≃ I − K ] > z ≃ M ⊙ . These objects contribute roughly 30% of the to-tal stellar mass density of the Universe at that epoch.McCarthy et al. (2004) estimated the age of red galaxiesat z = 1 . − . − −
500 M ⊙ yr − ). These massivegalaxies must have undergone a rapid formation processat z > K -band selected galaxies at 1 . < z < . M ≃ − × [ M ⊙ ]. Combining deep K20 spec-troscopy and HST-ACS imaging, Cimatti et al. (2004) dis-covered four old, fully assembled spheroidal galaxies at1 . < z < .
9: the most distant such objects currentlyknown. The stellar mass of these galaxies turned out tobe in the range 1 − × M ⊙ . Fontana et al. (2004) stud-ied K20 galaxies as well, showing that massive ( M > M ⊙ ) galaxies are easily found up to z ≃
2. These authorsalso report on the stellar mass function: only mild evolu-tion ( ∼ − z = 1, but only ∼
35% ofthe z = 0 stellar mass locked up in massive objects wasassembled by z = 2.At even higher redshifts, Rigopoulou et al. (2006) haverecently studied a population of z ∼ M ⊙ .McLure et al. (2006) identified nine LBGs at z ≥ . A stacking analysis suggests that the typical stellarmass of these sources is > × M ⊙ . Mobasher et al.(2005) analyzed the properties of J-dropouts in the HubbleUltra Deep Field (HUDF, Beckwith et al. 2006), exploit-ing Spitzer photometry. They identified a z ∼ . . × M ⊙ (but see,for example, Yan et al. 2004 for a different interpretation).Drory et al. (2005) studied the stellar mass function ofgalaxies in the FORS Deep Field (FDF, Heidt et al. 2003)and GOODS/CDFS (Giavalisco et al. 2004) field, over atotal area of 90 arcmin and found that the total stellarmass density at z = 1 is 50% of the local value. At z = 2,25% of the local mass density was already assembled, and erta S., et al.: IR-peakers mass function 3 at z = 3 and z = 5, at least 15% and 5% of stellar mass,respectively, was already in place. Massive ( M > M ⊙ ) galaxies existed over the whole redshift range probed,up to z = 5. The number density of these massive galaxiesevolves very similarly to galaxies with M > M ⊙ ,decreasing by 0.4 dex to z = 1, 0.6 dex to z = 2, and 1dex to z = 4.By analyzing the properties of K -selected galaxiesin the ∼
131 arcmin of the GOODS-CDFS survey,Caputi et al. (2006) found that the vast majority (85 − M > . × M ⊙ galaxies appears to bealready in place at z ∼
1. These authors also infer thatroughly 65 −
70% of these galaxies assembled at z = 1 − z = 3 − ∼
49 deg in the seven Spitzer chan-nels, and is therefore ideal to find rare objects. Its vol-ume is large enough to detect ∼
85 DM haloes of mass > M ⊙ in the 2 < z < L bol > L ⊙ ) and most mas-sive ( M ⋆ > several 10 M ⊙ ) galaxies ever to exist.The Infrared Array Camera (IRAC, Fazio et al. 2004),onboard Spitzer (Werner et al. 2004), samples the rest-frame near-IR light of distant galaxies. A near-IR selectionnot only directly probes the low-mass stars dominatingthe baryonic mass of a galaxy, but also is minimally af-fected by dust extinction. Therefore Spitzer is best suitedto the study of the stellar content of galaxies up to z = 3.Moreover Spitzer allows detection of galaxies that wouldbe missed by restframe UV selection, for example distantred galaxies (e.g. Daddi et al. 2004).We take advantage of the shape of near-IR spectralenergy distribution (SEDs) of galaxies to identify high-redshift objects on the basis of IRAC colors. Our selec-tion is based on the detection of the 1.6 µ m stellar peakin galaxies (Sawicki 2002; Simpson & Eisenhardt 1999),redshifted to the IRAC domain. It is also worth notingthat galaxies selected in this way benefit from a negativek-correction in the IRAC bands, because the slope of agalaxy’s SED is negative redward of the 1.6 µ m peak. Inthis way we performed a systematic search for M ∼ > M ⊙ galaxies at z > J , K s ) photometry and optical spectroscopy are available(Berta et al., 2006, Dias et al., in prep., La Franca et al.,in prep.). This area, and the sampled volume, are big-ger than any other previously explored for studying verymassive galaxies at z = 1 −
3, which were limited to verydeep, pencil-beam surveys (e.g. Dickinson et al. 2003a;Drory et al. 2005; Gwyn & Hartwick 2005; Fontana et al.2006). This paper is structured as follows. Section presentsthe data available in the ELAIS-S1 field; in Sect. wepresent our selection criterion; then Sect. discusses thephotometric estimate of redshifts. Section deals withthe estimate of the stellar mass in galaxies, and presents avery detailed analysis of the influence of mid-IR and near-IR constraints on it. Experiments with different stellar li-braries and IMFs are also discussed. Section presents thestellar mass function of our galaxies, including complete-ness correction and a parametric fit based on a Bayesianformalism. Finally, Sects. and discuss results and drawour conclusions.Throughout this work, we adopt a standard H = 71[km s − Mpc − ], Ω m = 0 .
27, Ω Λ = 0 .
73 cosmology, unlessotherwise stated.
2. Available Data
The ELAIS-S1 field (RA = 00 h m s , Dec = − ◦ ′ ′′ , J2000.0) represents the minimum of theMilky Way 100 µ m cirrus emission (Schlegel et al. 1998)in the southern hemisphere, with an average emissivity of0.38 MJy/sterad. This area was targeted by the Spitzer Space Telescope(Werner et al. 2004), as part of the
Spitzer Wide areaInfra-Red Extragalactic survey (SWIRE, Lonsdale et al.2003, 2004). A total of ∼ were observed with theIRAC (Fazio et al. 2004) and MIPS (Rieke et al. 2004)cameras.Data processing is described by Surace et al. (2004),Surace et al. (in prep.), Shupe et al. (in prep.) and AfonsoLuis et al. (in prep.). It consists of Basic Calibrated Data(BCD) by the SSC pipeline plus post-processing aimedat artifact removal, mosaicking and source extraction.Mosaicking was performed with the SSC routine MOPEX,and source extraction with SExtractor (Bertin & Arnouts1996). For unresolved sources (i.e. the case examinedhere), IRAC fluxes were extracted through a 1.9 ′′ diameteraperture and corrected to total fluxes following SSC pre-scriptions; MIPS fluxes were extracted by means of PRFfitting (see Surace et al., and MIPS Data Handbook 2006).The resulting 5 σ depths are 4.1, 8.5, 43, 48, 400 µ Jy at3.6, 4.5, 5.8, 8.0, and 24 µ m respectively. Observations at70 and 160 µ m lead to 26 and 166 mJy 5 σ limits (Suraceet al., Shupe et al., Afonso Luis et al., in prep.).For more details on Spitzer data reduction, calibrationand catalog extraction see the SWIRE delivery documen-tation (Surace et al. 2004, and following releases ). available at the Spitzer Science Center Legacy Programweb page, http://ssc.spitzer.caltech.edu/legacy/ Berta S., et al.: IR-peakers mass function The ELAIS-S1 area is the target of extensive optical fol-low up carried out with ESO telescopes: the
ESO-SpitzerImaging extragalactic Survey (ESIS, Berta et al. 2006).The optical ancillary data cover ∼ in the B , V , R bands, observed with the Wide Field Imager (Baadeet al. 1999, on the 2.2m ESO-MPI telescope in La Silla)and in the I , z bands, obtained with VIMOS (Le F`evreet al. 2002, on the VLT). The B , V , R observations andresults over the central 1.5 deg area were presented inBerta et al. (2006), while the VIMOS data are in an ad-vanced reduction stage (Berta et al., in prep.).The reduction of B , V , R data was performed by us-ing IRAF standard tasks and self-built routines. Sky flat-field frames were acquired during each night and appliedto images obtained on the same date. In order to obtaina uniform background and photometric zeropoint acrossthe field, super-sky-flat frames were built and applied foreach observing night, taking care of masking very brightsources.Astrometric calibration was performed on observationsof Stone et al. (1999) astrometric fields. In order to ac-curately map coordinates, a TNX algorithm was chosen,combining a gnomonic projection and non-linear polyno-mial distortions. The r.m.s coordinate difference betweenthe ESIS and the GSC 2.2 catalogs turned out to be ∼ . BV R color-curves were built, and mag-nitudes were transformed to the Johnson-Cousins stan-dard photometric system. Catalogs were extracted by us-ing SExtractor (Bertin & Arnouts 1996), Kron total mag-nitudes are adopted for extended sources, while aperturemagnitudes, corrected using the observed PSF, are usedfor unresolved objects. The final catalog reaches 95% com-pleteness at
B, V ≃
25 and R ≃ . J and K s imaging The square degree at the center of the ELAIS-S1 areaincludes J and K s observations carried out with the SOFI(Moorwood et al. 1998) camera on the NTT, during fourdifferent periods in 2002 and 2003 (Dias et al., in prep.).Pre-reduction, sky-subtraction and mosaicking werecarried out in the standard way, using the IRAF environ-ment. The astrometric mapping was calibrated by usingthe ESIS R band images (Berta et al. 2006), reaching a0.22 ′′ and 0.16 ′′ r.m.s. uncertainty for RA and Dec, re-spectively. Photometric calibration was computed using aperture photometry of point-like objects. The average J and K s magnitude difference for the sources in commonwith the 2MASS catalog is ∼ .
02 mag in both bands.Catalog extraction was performed using SExtractor(Bertin & Arnouts 1996). In this work we make use oftotal magnitudes. The resulting catalog is 95% completeto J ∼ . Ks = 18 .
73 (Vega).The optical, near-IR and SWIRE/Spitzer catalogswere matched with a simple closest-neighbor algorithm,adopting a 1.0 arcsec matching radius (see Berta et al.2006).
The central area benefits from X-ray observation by theXMM-Newton telescope (Puccetti et al. 2006). The wholeELAIS-S1 field was observed at 1.4 GHz with the ATCAradio telescope by Gruppioni et al. (1999) down to a 80 µ Jy r.m.s. and by Middelberg et al. (submitted) down to30 µ Jy r.m.s. Finally, the ultraviolet Galaxy EvolutionExplorer (GALEX, Martin et al. 2005) Deep ImagingSurvey (DIS) included the ELAIS-S1 field.
The ELAIS-S1 area has been spectroscopically surveyedby different projects, in particular focused on the centralfield.La Franca et al. (2004) performed optical spec-troscopy of ISOCAM 15 µ m counterparts, with the2dF/AAT, ESO-Danish 1.5m, ESO 3.6m and NTT tele-scopes, over the spectral range between 4000 − − K -selected galaxies ( K s < . µ m SWIRE objects. Further spectroscopy of opticallybright ( R <
21) X-ray and mid-IR sources was obtainedat the 3.6m/ESO telescope, during Fall 2005. These ob-servations will be presented in future works by La Francaand collaborators.Overall, at the present time, ∼
3. Selection criteria
In order to select high redshift galaxies dominated by stel-lar emission in the restframe near-IR, we have exploitedthe “IR-peak” technique, described in Lonsdale et al. (inprep.) and Berta et al. (2007).The near-IR emission of a galaxy is characterized by apeak centered at 1.6 µ m, due to the combination of thePlanck spectral peak of low-mass stars (dominated bytype M), a minimum in the H − opacity in stellar atmo-spheres and molecular absorptions in the spectra of coldstars. The Spitzer IRAC camera was partly designed todetect this feature in high redshift galaxies (Sawicki 2002;Simpson & Eisenhardt 1999). The peak is fully sampled erta S., et al.: IR-peakers mass function 5 by the IRAC photometric bands when it falls in the 4.5 or5.8 µ m channel, i.e. when at least one band lies shortwardand one longward of the peak. We thus focus our anal-ysis on 4.5 µ m-peakers and 5.8 µ m-peakers (objects withSEDs peaking in the 4.5 µ m and 5.8 µ m channels). The 1.6 µ m peak is redshifted to the 4.5 and 5.8 µ m channels forsources in the redshift range z = 1 . − h m s Dec = − ◦ ′ ′′ , where allBVR, J , K s , IRAC and MIPS photometric data are avail-able, for a total area of 1 deg . This area includes 48101SWIRE sources, and 29859 SWIRE+ESIS matches. Inthis area, 2848 objects are detected at 24 µ m by the MIPScamera, while the remaining Spitzer sources are detectedby IRAC only.We first apply a 5.8 µ m flux cut at the 3 σ SWIREdepth (25.8 µ Jy), in order to favor sources detected in theIRAC bands, which sample the restframe near-IR emissionand hence the stellar mass (as dominated by low massstars). This sub-sample consists of 5546 objects (over onesquare degree).
The IR-peak selection is well exemplified in the IRACcolor space (Lacy et al. 2004; Stern et al. 2005), as shownin the top panels of Fig. 1. Here red triangles repre-sent 5.8 µ m-peakers: formally sources showing S ν (3 . S ν (8 . µ m-peak objects, characterized by S ν (3 . S ν (5 . µ m band, butfor which the 8.0 µ m 3 σ upper limit is consistent with thedefinition of the IR-peak. The small dots represent thegeneral galaxy population detected in the central ELAIS-S1 SWIRE field, light blue having a 4 band detection, andorange having a 8.0 µ m upper limit only.We overplot a set of template tracks as a function ofredshift: a Seyfert-1 (dotted lines, Mrk 231, Fritz et al.,2006), a Seyfert-2 (dot-dashed, IRAS 19254-7245, Bertaet al., 2003), a starburst galaxy (solid line, M82, Silvaet al., 1998, enhanced with observed PAHs by F¨orster-Schreiber et al. 2001), a spiral (long-dashed, M51, Silvaet al., 1998) and an optically-blue spiral (short-dashed,NGC4490, Silva et al., 1998) galaxy. The tracks are limitedto the range z = 1 −
3, for the sake of clarity, apart forthe starburst one, which is extended down to z ≃ IR-peak-like colors could in principle be produced not onlyby the 1.6 µ m feature redshifted in the IRAC domain, butalso by a strong 3.3 µ m PAH feature at lower redshift, as isdemonstrated by Lonsdale et al. (in prep.) and Berta et al.(2007). In this case, the real 1.6 µ m peak lies shortward of the IRAC channels, resulting in a blue NIR-IRAC ob-served color.As far as 5.8 µ m-peakers are concerned, the require-ment that S ν (3 . < S ν (4 .
5) automatically avoids low-redshift interlopers, because it forces rejection of thosesources for which the 1.6 µ m peak falls at wavelengthsshorter than the IRAC 3.6 µ m channel. On the other hand, in the case of 4.5 µ m-peakers, onlyone IRAC band samples the SED on the blue side of theobserved peak and there is no trivial way to avoid aliaseswith low-redshift objects with a bright 3.3 µ m PAH. Inthis case the K s data are of fundamental importance inbreaking the degeneracy and ruling out interlopers. Themiddle panel in Fig. 1 involving the K s magnitude showshow near-IR data can effectively break this degeneracyand help in selecting z = 1 − K s − . AB < z <
1. Thereforewe adopt this value in order to disentangle low- and high- z sources.In the two top panels of Fig. 1 (discussed above), opencircles represent objects with ( K s − . AB <
0, while filledcircles lie above this threshold. In the top left diagram(Lacy et al. 2004), these low-redshift interlopers seem tobe separated from the high-redshift objects and we empir-ically define the transition line:(4 . − . AB = 1 . × (3 . − . AB + 0 .
7. (1)We then use this line to distinguish between low- and high- z candidates among those sources that are not detected inthe near-IR ( J , K s ) survey.It is also interesting to test whether the optical-IRACcolor space shows any segregation of low- z versus high- z objects (as defined on the basis of the IR-peak techniqueand K s band photometry), in order to define alternativeselection criteria to be used when no near-IR data areavailable. The bottom plots in Fig. 1 show a couple ofdiagrams involving the available optical BVR photometryfrom the ESIS survey (Berta et al. 2006). The majority ofIR-peakers with ( K s − . AB < R C − . AB < R C − . AB = 1 . × ( B J − R C ) AB . (2)In the end, combining the IR-peak criterion and the near-IR and optical constraints, our sample contains 149 and231 sources peaking at 4.5 µ m and 5.8 µ m respectively,over one square degree. Among these, 149 (i.e. ∼ ∼ J or K s , and 109 ( ∼ µ m counterpart. Finally, ∼
77% of the IR-peakers haveno 8.0 µ m detection. Table 1 summarizes these numbers. Berta S., et al.: IR-peakers mass function
Fig. 1.
Selection of the IR-peak galaxies. The two top panels show the position of 4.5 µ m-peakers (solid circles) and5.8 µ m-peakers (triangles) in the IRAC color space (Lacy et al. 2004; Stern et al. 2005). The middle plot highlights theuse of K s band data to recognize low-redshift interlopers aliasing 4.5 µ m-peaker colors. The open circles in the other 4panels belong to sources with ( K s − . < bottom panels present additional informationbased on optical-IRAC colors. Small dots represent the general galaxy population in the SWIRE ELAIS-S1 field.Template tracks are overlaid for a Sy-1, Sy-2, starburst (M82), spiral (M51) and blue spiral (NGC4490) galaxies. Thecovered redshift range is reported aside each track. erta S., et al.: IR-peakers mass function 7 Fig. 2.
Effect of sky-noise on IR-peakers selection. Random fluctuations of IRAC and K s fluxes can cause objects tobe scattered in and out of the sample. Solid lines are the actual distribution of IR-peakers. Shaded areas are the resultof simulations and represent the variation in the number of peakers due to random fluctuations of their colors. Thethree vertical lines refer to the 3, 4, 5 σ thresholds at 5.8 µ m in the SWIRE ELAIS-S1 survey. Description NumberSWIRE sources 48101SWIRE + opt. sources 2985924 µ m sources 2848 S (5 . ≥ . µ Jy (3 σ ) 5546IR-peakers 380IR-peakers 4 σ µ m-peakers 1494.5 µ m-p. + 24 µ m 444.5 µ m-p. + optical 694.5 µ m-p. + NIR 1404.5 µ m-peakers 4 σ µ m-peakers 2315.8 µ m-p. + 24 µ m 655.8 µ m-p. + optical 805.8 µ m-p. + NIR 1555.8 µ m-peakers 4 σ Table 1.
Summary of IR-peaker selection in the centralsquare degree of the SWIRE ELAIS-S1 field. The descrip-tions “+ optical” and “+ NIR” refer to sources with atleast one optical or JK detection.
Random fluctuations of fluxes within the photometric un-certainties in the bands defining the IR-peaker selectioncan cause non-peaker objects to be scattered into the sam-ple and actual peakers to fall out of it.It is worth noting that, since the 3.6 µ m and 4.5 µ mdetection thresholds are much fainter than at 5.8 µ m (seeSection 2), our sources are detected above ∼ σ in the two bluest IRAC channels, and sky-noise affects their fluxes bya factor smaller that ∼ µ mand 8.0 µ m bands can significantly modify the colors ofthese objects and scatter them in and out of the sample.In order to quantify this effect, we have performed abootstrap simulation of IR-peaker selection, starting fromthe multi-wavelength SWIRE catalog in the ELAIS-S1field (48101 sources). The procedure applies a randomfluctuation to the K s -to-8.0 µ fluxes of all SWIRE sources,with a dispersion given by the sky-noise associated to eachsource.This process was looped 10000 times and a new IR-peaker selection was performed at each step. Figure 2 re-ports the result of this simulation, for 4.5 µ m-peakers (leftpanel) and 5.8 µ m-peakers (right). Solid lines represent theactual distribution of IR-peakers, while shaded areas arethe result of simulations. Vertical dashed lines representthe 3, 4, 5 σ detection thresholds in the SWIRE ELAIS-S15.8 µ m survey.This analysis shows that, at the 3 σ detection levelrandom noise fluctuations significantly affect the selectionof IR-peakers, and tend to diverge. On the other hand atthe 4 and 5 σ flux levels the contamination decreases to ∼
20% and ∼
15% respectively.Therefore the subsequent analysis will be carried outonly for those IR-peakers detected above the 4 σ level (34.4 µ Jy) in the 5.8 µ m band. Thus the sample of IR-peakersreduces to 326 objects, 123 of which are 4.5 µ m-peakersand the remaining 203 are 5.8 µ m-peakers. The electronicTable associate to this work reports the main data of thefinal selected sample. Berta S., et al.: IR-peakers mass function
4. Photometric redshifts
The spectroscopic survey in ELAIS-S1 has targetedmainly X-ray sources and generally low-redshift ( z ≤ µ m excesses due to torus warm dust.These objects do not belong to the IR-peaker sample.As a consequence, we need to rely mostly on photo-metric redshift. Nevertheless, recent spectroscopic anal-yses of subsamples of IR-peakers in other SWIRE fields(Lockman Hole, ELAIS-N1, ELAIS-N2) have confirmedthat the adopted selection effectively identifies galaxiesat z ≃ . − .
0. Weedman et al. (2006) performed IRS(Houck et al. 2004) mid-IR spectroscopy of IR-peakers de-tected at 24 µ m and found redshifts between z = 1 . − . z = 1 . Hyper-z code (Bolzonella et al. 2000). We adopt a semi-empirical template library including GRASIL (Silva et al.1998) models of spiral and elliptical galaxies, M82and Arp220 templates (Silva et al.) upgraded withobserved PAH mid-IR features (F¨orster Schreiber et al.2001; Charmandaris et al. 1999), and a ULIRG template(IRAS 19254-7245, Berta et al. 2003).The
Hyper-z performance has been optimized forgalaxies: we have excluded from the analysis all sourceswith a type-1 AGN spectroscopic classification and allthose showing a power-law like IRAC SED. In this wayAGN-dominated objects are avoided. This class is partic-ularly delicate to fit with a template-based procedure andphotometric redshifts are typically mis-interpreted in 50%of the cases (see for example Berta et al. 2007), becausesharp features are missing in their SEDs. Moreover, theIR-peak criterion automatically selects the sample againstpower-law AGNs.The estimate of photometric redshifts was performedincluding the optical, J , K s and IRAC (3 . − µ m) datain the χ computation, but ignoring the MIPS 24 µ m flux.In fact, including the 24 µ m data turned out to produce ahigher degree of degeneracy and aliasing.The choice of templates in the library and the allowedextinction have been tested on galaxy populations in theELAIS-S1 field, taking advantage of the available spec-troscopic redshifts. Figure 3 reports the comparison ofspectroscopic and photometric redshifts, as obtained with Hyper-z , The dashed and dotted lines represent ± ±
20% uncertainties. The difference between the twoestimates has an r.m.s. of 0.19 and s.i.q.r. of 0.076. s.i.q.r. = semi inter-quartile range, defined as ( q − q ) / q and q are the first and third quartiles. The firstquartile is the number below which 25% of the data are foundand the third quartile is the value above which 25% of the dataare found. Fig. 3.
Comparison of photometric and spectroscopic red-shifts for all galaxies with available spectroscopy in theELAIS-S1 field. Objects classified as type-1 AGNs on thebasis of spectroscopy or having a power-law like IRACSED have been excluded from this analysis. Dashed anddotted lines represent ±
10% and ±
20% uncertainty levels.
Fig. 4.
Distribution of photometric redshifts of IR-peakers in ELAIS-S1, as derived with the code
Hyper-z (Bolzonella et al. 2000).The main outliers are sources with few photometricdata available, especially those with a detection in thetwo more sensitive IRAC channels (3.6 µ m and 4.5 µ m) andone optical band (typically R) only. The photometric codecalibrated in this way was then used to derive redshiftsfor the general IR-peak population selected as describedabove. erta S., et al.: IR-peakers mass function 9 The redshift distribution of IR-peakers is shown in Fig.4, where 4.5 µ m-peakers are represented by the dotted lineand 5.8 µ m-peakers by the dashed one.
5. The estimate of stellar mass
The estimate of the stellar mass of the galaxiesin our sample has been obtained with the code
Sim-Phot-Spec (SPS), developed in Padova (Berta et al.2004; Poggianti et al. 2001).This code performs mixed stellar population (MSP)spectro-photometric synthesis. Several phases in the life ofa simple stellar population (SSP) are combined together,adopting a different star formation rate (SFR) for eachage. Each SSP phase is meant to represent a formationepisode of average constant SFR, over a suitable time pe-riod ∆ t . The effective number of SSP ages involved in thefit depends on the redshift of the given source, in ordernot to exceed the age of the Universe at that redshift. Themaximum number of SSP phases considered for our sam-ple is 6, corresponding to a redshift z = 1, and decreases to4 if z = 3. Moreover, during the minimization phase, thecode checks for SSPs that do not effectively contribute tothe emitted light, either in the optical and IR, and deletesthem. In the end the number of SSPs contributing to thefinal fit is typically 3-4.Age-selective extinction is applied, assuming that theoldest stars have abandoned the dusty medium long ago.Keeping in mind that disc populations are on average af-fected by a moderate A V ( < A V = 0 . − . A V ≤ χ minimization of the dif-ferences between the observed photometric data and thesynthetic SED, including photometric uncertainties in allbands. Photometric redshifts, as computed with Hyper-z (Bolzonella et al. 2000, see Sect. 4) are adopted. TheSPS code takes into full account the uncertainty associ-ated with the photometric redshift estimate, while explor-ing parameter space. This uncertainty is of the order of0 . − .
2, depending on how many photometric bands areavailable. When the given source is not detected in one(or more) bands, upper limits are used.For each galaxy in the sample, the SPS code builds alarge number of models ( ∼ ), exploring the parameterspace by means of the “Adaptive Simulated Annealing”algorithm (ASA, Ingber 2001, 1989; Berta et al. 2004).First the parameters are varied with a coarse resolution,and then the algorithm focuses on the found minima. Inorder to avoid “freezing” in local minima, re-annealing isapplied and the code literally “jumps” out of the mini-mum in order to explore a different region of the param-eter space and find the absolute best fit. In this way, thepossible fluctuations of colors, within the measured pho-tometric errors, are accounted for and are automaticallypropagated to the output stellar mass uncertainty. The SPS code is optimized to derive the assembledstellar mass in galaxies and the associated uncertainty dueto degeneracies in star formation history (SFH) space (seeBerta et al. 2004). The χ of each model considered dur-ing the minimization (i.e. each combination of parametervalues) is recorded and 1, 2, 3 σ contours are computed. Inthis way, the resulting uncertainty on the stellar mass esti-mate accounts also for minimal- and maximal-mass mod-els, derived from the projection of the parameters spaceinto “young” and “old” sub-spaces. The resulting best fitmasses and 3 σ mass ranges are reported in an electronicTable.The adopted SSP library is based on thePadova evolutionary sequences of stellar models(Fagotto et al. 1994a,b; Bressan et al. 1993) andisochrones (Bertelli et al. 1994), and was computedby assuming a solar metallicity and a Salpeter initialmass function (IMF) between 0.15 and 120 M ⊙ . TheSSP spectra were built with the Pickles (1998) spectrallibrary, and extended to the near-IR with the Kurucz(1993) atmosphere models. Nebular features were addedthrough the ionization code CLOUDY (Ferland 1996).The spectra thus obtained provide a reliable descriptionof simple stellar generations up to ∼ µ m (restframe).Beyond this wavelength, dust emission is no longernegligible.We assume that the total energy absorbed in the UV-optical domain is processed by dust in the thick molecularclouds embedding young stars and re-emitted in the mid-and far-IR (8-1000 µ m) in the form of a starburst tem-plate. By convolving the template with the 24 µ m pass-band, the MIPS 24 µ m flux expected for the given modelis computed. Finally, this synthetic flux is compared tothe observed 24 µ m data and included in the χ computa-tion. Mid-IR photometry is mainly sensitive to the powerof the ongoing starburst, as well as to the amount of dustobscuring it, therefore it provides a valuable constrainton the amount of dust and the strength of young stellarpopulations (see also Berta et al. 2004).Typically an M82 template is adopted, but a dif-ferent choice would not significantly affect the result-ing stellar mass estimate, as demonstrated in Berta et al.(2004). Increasing observational evidence exists that high- z IR-peakers detected in the mid-IR resemble the M82prototype. Based on Spitzer mid-IR IRS spectroscopy,Weedman et al. (2006) found that z ≃ . µ m absorption. Rowan-Robinson et al. (2005)studied and classified the SEDs of SWIRE sources over6.5 deg in the ELAIS-N1 fields, finding that M82-likestarbursts are 3 times more numerous than colder Arp220-like objects. Finally, millimeter (250 GHz) observations ofSWIRE 24 µ m-selected IR luminous galaxies, performedwith the MAMBO bolometer on the IRAM/30m telescope(Lonsdale et al., in prep.) showed that the 1.2mm/24 µ mflux ratio of these sources resembles that of M82, lowerthat for an Arp220-like population. Fig. 5.
Examples of SED fits of IR-peak sources. The top panels belong to 4.5 µ m peakers, the bottom ones show two5.8 µ m-peak galaxies. The two sources on the left are not detected at 24 µ m, and only an upper limit is available, whilethose on the right have an actual 24 µ m flux above 250 µ Jy (i.e. the 3 σ SWIRE limit in ELAIS-S1). The dashed linerepresents the contribution to the best fit model by young and intermediate-age ( < ≥ µ m (restframe). Thelong-dashed line longward of 5 µ m (restframe) is the adopted starburst template re-processing the UV-optical lightabsorbed by dust to the IR.Figure 5 shows a few examples of SED fits of IR-peaksources. Overplotted on the observed fluxes is the bestfit solution of the MSP synthesis: the dashed line refersto young-intermediate (age < ≥ µ m (restframe), the long-dashed line is the starburst template used to model the IR emission from dust, heatedby young stars. The open triangle represents the 24 µ mflux predicted by the model, which often is beneath theobserved data-point.Figure 6 reports the distribution of stellar masses asa function of redshift. Circles belong to 4.5 µ m-peakers,while triangles represent 5.8 µ m-peak sources. erta S., et al.: IR-peakers mass function 11 Fig. 6.
Distribution of stellar masses as a function of redshift. Circles and triangles represent 4.5 µ m- and 5.8 µ m-peakers respectively. Left panel: filled and open symbols belong to sources dominated (best fit) by young or old stellarpopulations. The mass thresholds expected for the maximal and minimal M ⋆ /L ratios in the sample are overplottedas dashed (5.8 µ m-peakers) and dotted (4.5 µ m-peakers) lines. Histograms are plotted for 4.5 µ m-peakers (dotted) and5.8 µ m-peakers (dashed). The solid histogram is the total distribution of sources. Right panel: filled and open symbolsbelong to sources with or without a 24 µ m detection. Histograms have the same meaning, dotted ones representingobjects detected by MIPS and dashed ones objects not detected.In the left panel, open symbols indicate objects whosebest fit synthetic SEDs are dominated by old (age ≥ yr) stars, while for filled circles the bulk of the luminos-ity is provided by younger populations. The dashed anddotted histograms belong to 5.8 µ m- and 4.5 µ m-peakers.Objects with a resulting lower mass are dominated byyounger stellar populations, as expected on the basis ofthe dependence of the M ⋆ /L ratio on age. The effect ofMalmquist bias is evident not only from the trend seen inthe M ⋆ vs. z plot, but also from the distribution of sourcesin the mass space: 4.5 µ m-peak galaxies have masses lowerthan 5.8 µ m-peakers, on average.Typical error bars are shown: uncertainties on the stel-lar mass can be as high as 0.3 dex for sources that lack two( B and V ) or all optical bands and — in the worst cases— also near-IR photometric data. This situation is repre-sented by the big error bar. When good multi-wavelengthcoverage is available, the uncertainty on the stellar massreduces significantly (small error bar).The mass thresholds expected at the 4 σ flux limitof S (5 . µ m) = 34 . µ Jy, as derived from the maximaland minimal M ⋆ /L ratios in the SED fitting results, areoverlaid on the data. Dotted lines belong to 4.5 µ m-peakersand dashed lines to 5.8 µ m-peakers. See Sect. 6.2 for adiscussion on the completeness of the sample.The right panel of Fig. 6 shows the same plot, butcoded on the basis of 24 µ m detections. Filled and opensymbols refer to sources detected or not detected in theMIPS 24 µ m channel, respectively, at the 3 σ level. Thedashed and dotted histograms have the same meaning. µ m data on the mass The availability of mid-IR observations provides a valuabletool for constraining the amount of dust and the luminos-ity of young stellar populations in the synthesized models,thus reducing the uncertainty on the derived stellar mass(e.g. Berta et al. 2004).In order to understand the effect of the 24 µ m flux onthe mass, i.e. how it effectively influences the star forma-tion history of the best fit solution, we have performed afitting run ignoring the available 24 µ m photometry.The results are shown in Fig. 7. The top-left panelshows the distribution of stellar masses as a function ofredshift (symbols are the same as in Fig. 6), while thetop-right plot compares the new results (solid line, ob-tained excluding the 24 µ m flux from the minimization)to the standard fit mass distribution (dashed histogram).The two bottom panels present the same comparison, forsources dominated by young and old stellar populationsrespectively (more that 50% of the total mass being as-sembled in stars younger or older than 1 Gyr).If no mid-IR detection is available, the amount of dustand young stars are free to vary with no control, andvery extinguished ( A V ≥
4) young populations, not vis-ible in the optical domain, may exist. As a consequence,the amount of young stellar mass is found to be larger, thetotal stellar mass inferred from SED fitting smaller on av-erage, and the mass spread slightly wider. The availabilityof MIPS data provides a valuable constraint on the recenthistory of star formation of galaxies and helps to avoid asignificant number of SED fitting solutions that could notbe a priori ruled out.
Fig. 7.
Results of SED fitting, ignoring the 24 µ m fluxes (solid histogram), and comparison to the fit including 24 µ m(dashed, see Fig. 6). In the top left panel, filled and open symbols represent sources dominated (best fit) by youngor old stellar populations, respectively. The mass distributions for these two sub-classes are shown in the two bottom panels. The near-IR J and K s data turned out to be very usefulat the selection stage, in particular to exclude low-redshiftinterlopers with IRAC colors similar to 4.5 µ m-peakers.It is worthwhile exploring how observations at thesewavelengths (1-2 µ m) constrain the estimate of stellarmasses, by means of SED fitting. We have therefore per-formed another fitting run, without taking into accountthe J , K s photometry. Figure 8 shows a couple of SEDsand the comparison to masses from the standard fit.The top panels include the SEDs of a 4.5 µ m-peaker(left) and a 5.8 µ m-peak galaxy (right). The dashed linerepresents the fit obtained with the standard technique,i.e. accounting for all the available photometric data, whilethe dotted lines are the best fit models obtained by exclud-ing the near-IR ( J and K s ) fluxes from the minimization.In some cases, the J , K s data turn out to be use-ful to constrain the age of the dominant stellar popula-tions hosted by IR-peakers, since they sample the depthof the D4000 break, when combined with optical data.Nevertheless, the J , K s constraint turns out to only beeffective, in this respect, for a small fraction of cases( ≤ µ m-peakers than for objects that peakat longer wavelength (in the 5.8 µ m channel, i.e. at higherredshift). For the latter, the SED slope blueward of the1.6 µ m peak is, in fact, already defined by the 3 . − . − . µ m colors well enough to reasonably constrain the D4000break, and therefore the best fit solutions obtained withor without J , K s do not differ too much.Several authors (e.g. Berta et al. 2004) had shown thatthe introduction of IRAC data in SED fitting providestighter constraints on the stellar mass, reducing its un-certainty by a factor as high as 5, for z > JHK band datacan reduce the uncertainty on the stellar mass of bluegalaxies with ( U − V )rest <
1, while IRAC photometry isreally effective only for redder sources, having a R − IRACcolor significantly redder than 1; this is the case for ourIR-peakers.Overall, the distributions of stellar masses, as obtainedwith and without the near-IR J and K s data (Fig. 8), arenot significantly different, confirming the result that thesedata do not play a fundamental role in deriving the stellarmass of IR-peakers. On the other hand, as we have alreadypointed out, near-IR magnitudes are critical to avoid low- z aliases. We also expect these data to be more effective inconstraining the age and M ⋆ /L ratio of galaxies in other erta S., et al.: IR-peakers mass function 13 Fig. 8.
Comparison of SED fitting with and without the J , K s data. The top panels show an example of 4.5 µ m- ( left )and 5.8 µ m-peaker ( right ). The dotted line shows the best fit obtained ignoring the real J , K s data, while the dashed lineis the standard fit. The three bottom panels show the distribution of masses for all sources and young/old-dominatedgalaxies, for the standard fit (dashed histograms) and the fit obtained ignoring the J , K s data (solid lines).redshift bins, e.g. when the blue side of the 1.6 µ m featureis not fully sampled by IRAC, z ≤ . The choice of an initial mass function (IMF) different fromSalpeter’s can significantly affect the estimate of the stel-lar mass.The Salpeter (1955) IMF has a constant α = 1 .
35 slopethroughout the whole mass range considered, but in factit was never measured down to 0.15 M ⊙ by Salpeter.More recent determinations of the IMF showed that aflatter slope is needed at low masses ( M ≤ ⊙ ), in orderto reproduce the observed Galactic data (Miller & Scalo1979; Kennicutt 1983; Kroupa et al. 1993; Kroupa 2001;Chabrier 2003a, e.g.) Consequently, assuming Salpeter’sIMF results in including too many low-mass stars (whichdominate the stellar mass budget) in the galaxy modeling.Introducing a drop-off at M ≤ ⊙ leads to lower valuesof the stellar mass of the analyzed galaxies.At the bright end, Elmegreen (2006) pointed out thatabove 1 M ⊙ the IMF slope is not steeper than Salpeter’s;Miller & Scalo (1979) and the other IMFs with a drop atthe highest masses were based on galactic disk measure-ments which cannot be safely used to trace the high-massend of the IMF, because of the complicated SFH of theGalaxy. We study here the effect of a different choice in theIMF, by performing a new fitting run with the Chabrier(2003a) IMF instead of the classic Salpeter (1955) one.The Chabrier (2003a) IMF is described by a power-law for M > ⊙ and a lognormal form below this limit.Chabrier (2003a) and Kroupa (2001) are very similar toeach other, but here we prefer to adopt the former be-cause it is physically motivated and provides a better fit tocounts of low-mass stars and brown dwarfs in the Galacticdisc (see for example Chabrier 2001; 2002; 2003b and alsoBruzual & Charlot 2003).The results of this analysis are reported in Fig. 9,where we show the shape of different IMF’s (top leftpanel), and the direct comparison between stellar massesas obtained with the Salpeter and Chabrier IMFs (topright panel). In the bottom panels, the comparison ofthe stellar mass distribution for the two cases is shown.Solid histograms represent stellar masses derived with theChabrier IMF and dashed ones with Salpeter’s.The diagrams show that the Chabrier (2003a) IMFleads to systematically lower masses than the Salpeter(1955), as expected. The difference between the two turnsout to be ∼ < yr) and old stars ( ≥ yr), a few objects migrate from Fig. 9.
Comparison of IR-peakers stellar masses as derived with a Salpeter (1955) and a Chabrier (2003a) IMFs.The top left panel shows different derivations of the IMF: Salpeter (1955, solid line), Miller & Scalo (1979, dotted),Kroupa (2001, dot-dashed) and Chabrier (2003a, dashed, almost superimposed over Kroupa’s). The top-right panelshows the direct comparison of stellar masses in the two examined cases, for 4.5 µ m-peakers (circles) and 5.8 µ m-peakers(triangles). The bottom panels show the stellar mass distribution for the Salpeter (dashed) and Chabrier (solid) fits.The three panels report the results for all sources and for objects dominated by young-intermediate (age < yr) orold ( ≥ yr) SSPs.the former to the latter sub-samples, but the overall dis-tributions are maintained.Despite the fact that the choice of a Chabrier (2003a)or Kroupa (2001) IMF restricts the inclusion of toomany low-mass stars, the majority of literature studieson the stellar mass function of galaxies are based on theSalpeter (1955) description of the IMF (e.g Fontana et al.2006, 2004; Drory et al. 2005, 2004; Franceschini et al.2006; Rudnick et al. 2006; Dickinson et al. 2003a;Papovich et al. 2001; Shapley et al. 2001; Sawicki & Yee1998). For this reason, we will derive the stellar massfunction of IR-peakers from the Salpeter-based SEDfitting, keeping in mind that a different choice (e.g.Chabrier 2003a) would produce a shift of the massfunction to lower masses (e.g. by about 0.3 dex). Recently, Maraston et al. (2006) have performed stellarpopulation synthesis of z = 1 . − . Z ≥ Z ⊙ . As a consequence, the MS lifetimeis longer for the Padova tracks and the RGB phase is de-layed, with respect to the Frascati models.In any case, the key difference of the Maraston (2005)approach is the way the thermally-pulsing asymptotic gi-ant branch (TP-AGB) phase is included in the evolution ofstellar populations, i.e. by means of a semi-empirical fuelconsumptions table, in contrast to “a posteriori” recipesused in isochrone synthesis. TP-AGB stars dominate thenear-IR (e.g. K band) luminosity for SSP ages between0.3 and 2 Gyr (for a Salpeter IMF and solar metallicity),while main-sequence stars dominate in the optical (e.g. V band).At the transition between the early-AGB branch andthe TP-AGB, i.e. when the TP-AGB phase begins, thenear-IR luminosity of stars significantly increases, andtherefore the near-IR mass-to-light ratio ( M ⋆ /L K ) of theSSP drops suddenly by a factor 3 −
5. In models withoutthis phase, M ⋆ /L monotonically increases, independent ofwavelength. At ages between 0.8 and 1 Gyr, the modelsbased on the Padova tracks have M ⋆ /L K larger than the erta S., et al.: IR-peakers mass function 15 Fig. 10.
Comparison of SSPs from the Padova 1994 (e.g.Fagotto et al. 1994a, solid lines) and the Maraston (2005)library (dotted), for a Salpeter IMF and solar metallicity.The critical age range, when the TP-AGB phase is active,is shown. The SEDs are normalized at 5000 ˚A.Maraston (2005) ones, as a consequence of the differenttreatment of the TP-AGB phase. At older ages, M ⋆ /L K is smaller for the Padova 1994 case, because of the coolerRGB phase.Maraston et al. (2006) model the observed SEDs ofhigh- z galaxies, adopting different star formation histories(SFH): an instantaneous burst, a exponentially decliningstar formation (SF), a prolonged burst, and a constantSF. As a result, they find that the Maraston (2005) li-brary leads to systematically younger best fit solutionsfor these systems, and hence to lower stellar masses, withrespect to the Padova 1994 library. The IRAC Spitzer dataturn out to be very useful in constraining the mass, andthe Maraston (2005) models provide a better fit to theobserved photometry.With this in mind, we have performed a similar anal-ysis on our sample of IR-peakers, adapting the Maraston(2005) SSP library to the SPS code and running it in thesame exact way as before. Figure 10 compares the SEDs ofthe Padova and Maraston SSPs, in the age range in whichthe TP-AGB phase is active.The results of this further analysis are reported in Fig.11, where the stellar masses derived with the Marastonlibrary (solid histogram) are compared to those obtainedearlier with our “standard” fit (dashed). The top-left panelreports the M ⋆ /L . ratio for the SSPs in the two librariesand the top-right plot shows the direct comparison of stel-lar masses as derived in the two cases. As far as the low-mass end of the distribution is concerned, stellar massesbased on the Maraston library turn out to be higher than in the Padova-94 case. Conversely, the Padova-94 fit pro-vides slightly lower masses at the bright end.The free-form MSP technique adopted for fitting theIR-peakers produces a different result, compared to thatof Maraston et al. (2006). The new library of models pro-duces best fit solutions characterized by completely differ-ent SFH’s, with respect to the previous ones, based on thePadova tracks. The relative contribution of young and oldstars in the best fits changes significantly, with the frac-tion of objects dominated by old populations increasingin the Maraston case. The two bottom plots show the M ⋆ distribution for sources dominated (in mass) by old andyoung stars: not only the balance between old and youngSSPs changes, but also the distributions peak at lowermasses (for old-dominated sources) and higher masses (foryoung-dominated ones), than in the Padova-94 case. Onthe other hand, the overall stellar mass distribution (topright panel) is not significantly changed, apart for a shiftof M < M ⊙ galaxies to higher masses. First of all, the sample of IR-peakers has been checkedagainst bright X-ray sources, thanks to the availableXMM-Newton survey in the area (Puccetti et al. 2006).Ten sources turn out to host an AGN X-ray component,seven 4.5 µ m-peakers and three 5.8 µ m-peakers. The for-mer are characterized by a bright 8.0 µ m excess, likely dueto the presence of torus warm dust in the mid-IR SED (seealso Berta et al. 2007, Lonsdale et al. 2007, for an exten-sive description of SED shapes). The three 5.8 µ m-peakersemitting in the X-rays show a smooth stellar peak, dilutedby the AGN dust. These sources have been excluded formfurther analyses.The remaining sources show a sharp IR-peak, and, inperforming the SED analysis, we have assumed that onlystellar emission contributes to the observed SEDs, whileno AGN component is present. In fact, the warm dust inthe AGN torus would emit at IRAC-MIPS wavelengths,producing a power-law SED or diluting the 1.6 µ m (rest-frame) stellar peak. Hence a sharp stellar peak shifted tothe 4.5 or 5.8 µ m bands should identify sources whosenear-IR emission is dominated by stars. Nevertheless, thepresence of a possible AGN component cannot be fullyexcluded.Berta et al. (2007) have analyzed UV-optical (rest-frame) spectra of IR-peakers observed with the Keck-I 10m telescope, selected in order to show no evidenceof AGN components on IRAC color-color diagrams (e.g.Stern et al. 2005; Lacy et al. 2004). Because of instrumentlimitations, the sample was limited to the brightest IR-peakers in the SWIRE survey, with r ′ ∼ < . z = 1 . − . µ m-peakers. In order to explain the pres- Fig. 11.
Comparison of stellar mass distributions as obtained with the Maraston (2005) SSP library (solid lines) andthe Padova 1994 (e.g. Fagotto et al. 1994a) models (dashed). The top left panel reports the restframe mass-to-lightratio in the 3.6 µ m IRAC band, as computed for the two SSP libraries, a Salpeter IMF and solar metallicity. The topright panel shows the direct comparison of stellar masses, derived in the two cases, for 4.5 µ m-peakers (circles) and5.8 µ m-peakers (triangles). The bottom panels show the comparison of the stellar mass distribution for all sources andobjects dominated by young/old stars (in the best fit).ence of the 1.6 µ m peak in the 4.5 and 5.8 µ m channels,a multi-component (stars and AGN) fit to the SEDs ofthese sources was performed. The AGN component par-tially dilutes the 1.6 µ m peak, and — more importantly —modifies its shape, resulting in an apparent shift to longerwavelengths.At optical restframe wavelengths, the AGN emissionis overwhelmed by stars in the host galaxy, while it re-emerges in the IRAC domain, especially in the two longestwavelength channels (5.8 and 8.0 µ m): torus dust is notnegligible for these sources, and thus might affect the1.6 µ m-peak selection. All the sources hosting an AGNwere detected at 24 µ m above 250 µ Jy, have moderate24 µ m excess ([3 . −
24] = 1 . − . R C − . AB = 2 − r ′ ∼ < .
8, corresponding to V J ≃ . ∼ V J = 23 . µ m, corresponding to ∼
6% ofthe sample, and finally only ∼ . −
24) color.We will not discuss this problem further, since the frac-tion of sources that might be affected is very small.
6. The mass function
The results obtained from spectro-photometric synthesishave been exploited to derive the contribution of z ≃ − ρ ⋆ .The derived stellar masses are used to build thestellar mass function of IR-peakers in two ways: withthe V a technique and applying the STY (Sandage et al.1979) method. The latter uses a parametric function,Φ ( M, param) and a Monte Carlo Markov Chain algorithmadopting a Bayesian formalism. erta S., et al.: IR-peakers mass function 17 Fig. 12.
Mass incompleteness.
Left and right panels belong to 4.5 µ m-peakers and 5.8 µ m-peakers, respectively. Thedistribution of stellar mass as a function of observed fluxes is shown for different redshift bins. The diagonal linesrepresent the tracks at the bin’s central redshift described by the minimal (dotted) and maximal (dashed) M ⋆ /L ratios in our samples. The vertical dashed lines set the adopted S (5 . µ m) = 34 . µ Jy flux limit. The horizontal linesrepresent the mass threshold above which the samples are fully complete. See text and Fontana et al. (2004) for detailson the technique adopted to correct incompleteness.
In order to compute the comoving number density of IR-peakers at the given redshift, we have adopted the well-known 1 /V a method (e.g. Schmidt 1968). The accessiblevolume , V a , in which each galaxy could be detected in thesurvey is computed by taking into account the selectioncriteria: – the flux cut at S (5 . µ m) = 34 . µ Jy, – the IR-peak color selection.The former provides the classical V max estimate (Schmidt1968), related to the maximum redshift z max at which thegiven galaxy would be observable in our survey. The z max is computed by adopting the best fit SED model (see Sect.5), redshifted and k-corrected until cosmological dimmingfades it below the adopted flux limit (e.g. Hogg 1999): S ( ν ) = L ( ν )4 πd L K ( ν, z ) , (3)where d L is the luminosity distance and K ( ν, z ) is the k -correction through the filter T ( ν ): K ( ν, z ) = (1 + z ) R ν ν L [ ν (1 + z )] T ( ν ) dν R ν ν L ( ν ) T ( ν ) dν . (4)The IR-peak condition defines a spherical shell, de-limited by redshifts z peakmin and z peakmax , where the IR-peakselection is valid, i.e. where the 1.6 µ m restframe peak is detected in the IRAC 4.5 or 5.8 µ m channel. The effectiveaccessible volume for each source in our color-selected sur-vey is thus given by: V a = Ω4 π Z min ( z max ,z peakmax ) z peakmin dVdz dz , (5)where Ω is the surveyed sky area, and the volume element dV depends on the adopted cosmology.It is worth noting that the ( K s − . > z peakmin . The IR-peaker sample was first selected by applying aflux cut at 5.8 µ m. Despite the fact that we are directlyprobing the restframe near-IR emission of these galaxies,which is primarily powered by low-mass stars dominatingtheir stellar mass budget, it is not possible to define asharp mass limit encompassing the whole sample.In fact, the mass-to-light ratios ( M ⋆ /L ) of these galax-ies spans a relatively wide range, between ∼ ∼ µ m (restframe). At a given redshift and flux, thesetranslate into very different mass values, as shown by theminimal and maximal galaxy tracks shown in Fig. 6 (leftpanel). As an example, it would thus be more correctto say that at z = 2 the sample is limited to masses M > × or M > . × M ⊙ for low or high M ⋆ /L objects, respectively.This effect has been thoroughly examined byDickinson et al. (2003b), Fontana et al. (2006) andFontana et al. (2004). The latter perform a very detailedanalysis and provide a valuable recipe to partially correctthe mass incompleteness of flux-limited samples.Figure 12 reports the distribution of stellar masses asa function of the observed 5.8 µ m flux for the galaxiesin our sample, split into different redshift sub-bins. Theleft panels refer to 4.5 µ m-peakers, while on the right side5.8 µ m-peak galaxies are shown. The vertical dashed linerepresents the adopted 5.8 µ m flux 4 σ cut. The diagonaldashed and dotted lines are the tracks described by thesources with minimal and maximal M ⋆ /L ratios in thesample (see also Fig. 6), at the central redshift of eachsub-bin.Since the majority of 4.5 µ m- and 5.8 µ m-peak sourcesare in the 1 . ≤ z < .
75 and 2 . ≤ z < .
50 sub-bins(see the redshift distribution in Fig. 6), we will use theseas reference to derive the threshold completeness mass forthe two different classes.The horizontal lines in the two panels focused on theseredshift sub-bins represent the stellar mass levels abovewhich the samples are definitely complete, because thesethresholds are set by the maximal M ⋆ /L observed values.These masses turn out to be ∼ × and 3 . × M ⊙ for the two populations. Below these values, the samplebecomes increasingly incomplete.Following Fontana et al. (2004), it is possible to com-pute the fraction of sources lost at each given mass, byexploiting the observed distribution of M ⋆ /L ratios andfluxes. We defer the reader to their paper for further de-tails on this technique. We keep those mass bins for whichthe fraction of missing sources turns out to be ∼ < .
5. Inthis way it was possible to extend our analysis down to1 . × and 1 . × M ⊙ for 4.5 µ m- and 5.8 µ m-peakers respectively (see Tab. 2). Note also that belowthese limits the amplitude of random flux fluctuationsdue to sky-noise is large and the contamination of theIR-peaker sample by other classes of sources can not becontrolled (see Sect. 3.4). Within the V a formalism, the comoving number of galaxiesper unit volume, in each redshift bin and in the mass bin∆ M , is obtained by:Φ ( M ) ∆ M = X i V ia ∆ M , (6)where the sum is made over all galaxies in the given massbin.Figure 13 shows the distribution of the comovingnumber of galaxies as a function of stellar mass (i.e.the “observed” stellar mass function) for 4.5 µ m-peakers(dashed line, filled circles) and 5.8 µ m-peakers (solid line, triangles). Table 2 reports the data. The stellar massescome from the spectro-photometric fit obtained with thePadova-94 library and the Salpeter IMF, and accountingfor all the available photometric data and upper limits.Filled symbols represent the mass function corrected forincompleteness, while open symbols show the data as ob-tained before applying any correction.The number density of 4.5 µ m-peak population turnsout to be significantly lower than that provided by 5.8 µ m-peakers. It is worth noting that the selection based on a5.8 µ m flux cut is optimized for the detection of the 1.6 µ mstellar peak in IRAC channel 3, while in order to havea comparable selection for 4.5 µ m-peakers we would haveneeded to apply a 4.5 µ m flux cut. The 5.8 µ m flux cutcorresponds to a restframe H band selection for 5.8 µ m-peakers and to a K band selection for 4.5 µ m-peakers.Assuming that the typical S ( H ) /S ( K ) flux ratio of agalaxy (e.g. a IR-peaker) is ∼ .
4, then the S (5 . µ m) =34 . µ Jy cut corresponds to S (4 . µ m) = 48 . µ Jy for a4.5 µ m-peak galaxy. As a confirmation of this effect, Figure14 shows the distribution of 4.5 µ m and 5.8 µ m fluxes forour sources.However, the 5.8 µ m-peakers sample includes also allthose galaxies not detected at 8.0 µ m, but still consistentwith the 5.8 µ m-peak selection, when using the 8.0 µ m up-per limit. Unfortunately, since the 5.8 µ m band is the leastsensitive among the 3.6, 4.5 and 5.8 µ m SWIRE/IRACchannels, it is not possible to perform a comparable 4.5 µ mselection. In fact, in order to include similar objects tothe 4.5 µ m-peakers sample, all the z = 1 − S (4 . µ m) ≃ . µ Jy, not detected at 5.8um, but stillconsistent with the 4.5um-peak selection, should be takeninto account. The resulting sources are detected only intwo IRAC channels (3.6 µ m and 4.5 µ m) and are barelydetected in J or Ks. Hence aliasing by low- z interlopersor contamination by power-law objects would be hard tocontrol.We conclude that our completeness correction is notable to control this selection deficit in the 4.5 µ m-peakerssub-sample. Therefore we will limit the parametric deriva-tion of the mass function to the 5.8 µ m peakers only.The error bars on the mass function plotted in Fig. 13account only for Poisson statistics. Nevertheless, the un-certainties on the stellar mass estimate, caused by degen-eracies in the SFH space and errors on photometric red-shifts, introduce a non negligible contribution to the un-certainty on the observed stellar mass function and mustbe taken into account.Due to these uncertainties, an object can move fromthe mass bin it formally belongs to and actually contributeto the stellar mass function in another mass regime. Theprobability for each galaxy to contribute to the comov-ing number density in other mass bins is given by the χ distribution obtained during SED fitting with the ASAalgorithm.For sake of clarity we do not plot these additional errorbars in Fig. 13, but they will be fully taken into account erta S., et al.: IR-peakers mass function 19 Fig. 13.
Mass function of 4.5 µ m-peakers (dashed line, circles) and 5.8 µ m-peakers (solid line, triangles) in the ELAIS-S1 area, as obtained with the standard spectro-photometric fit and the Padova-94 library (see Sect. 5). Filled symbolsbelong to the completeness-corrected mass function, while open ones represent the data before any correction wasapplied. Error bars account only for Poisson noise statistics. log(Mass) Φ( M ) 4.5 µ m-p. Φ( M ) 5.8 µ m-p. Φ( M ) 5.8 µ m-p (evo) Missing fraction[h − M ⊙ ] [10 − h Mpc − dex − ] 4.5 µ m-p. 5.8 µ m-p.11.15 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± Table 2.
Observed mass function of SWIRE IR-peakers, as obtained after correction for incompleteness effects, for aSalpeter IMF and Padova-94 stellar tracks. We report the measure values of the galaxy comoving number density andthe fraction of missing sources in each mass bin (see the discussion on completeness in Sect. 6.2). For 5.8 µ m-peakersthe mass function obtained after evolutionary correction is also included (column 4).when fitting the observed data with a parametric massfunction, by using a Bayesian approach (see Sect. 6.6). The IR-peakers mass function obtained with the Padova-94 library is compared to the results from the Maraston(2005) approach in Fig. 15. Data are corrected for incom-pleteness, in the same way described for the Padova-94case.As far as 4.5 µ m-peakers are concerned, the two li-braries lead to very different mass functions, mainly be-cause low-mass objects migrate to higher mass values,when the Maraston (2005) SSPs are used (see also Fig.11 and Sect. 5.4). Nevertheless, this mass function is still affected by a strong incompleteness effect, which cannotbe recovered.In the case of 5.8 µ m-peakers, the difference betweenthe two results is less dramatic and the two libraries leadto very similar mass functions.Note, however, that the comparison of our results toprevious estimates of the stellar mass function at highredshift (e.g. Fontana et al. 2006, 2004; Drory et al. 2005;Franceschini et al. 2006, among others), will make use ofthe results obtained exploiting the Padova-94 library. The mass density of galaxies increases very rapidly be-tween redshift 3 and 2, by a factor of ∼
4, according to
Fig. 16.
Left panel: mass function of 5.8 µ m-peakers split into three redshift sub-bins (same symbols as in Fig. 13). Right panel: evolutionary trend of the 5.8 µ m-peakers stellar mass density as a function of redshift for three differentmass bins. Solid lines connect the observed data, dotted lines represent a power law fit to the evolution. Fig. 14.
Distribution of 4.5 µ m and 5.8 µ m fluxes for IR-peakers. Circles represent 4.5 µ -peakers, triangles belongto 5.8 µ m-peakers. The vertical dashed line corresponds tothe S (5 . µ m) = 34 . µ Jy flux cut. The horizontal dashedline sets the S (4 . µ m) = 48 . µ Jy limit, obtained by as-suming a restframe S ( H ) /S ( K ) = 1 . µ m- Fig. 15.
Comparison of the mass function as obtainedwith the Padova-94 (filled symbols, solid lines) andMaraston (2005, open symbols, dotted lines) libraries, af-ter correction of incompleteness.peakers split for the first time into three redshift sub-bins:2 . ≤ z ≤ .
25, 2 . < z ≤ . . < z ≤ .
0. Notethat the integration boundaries in Eq. 5 are now set bythe used redshift sub-bins.The number of massive galaxies decreases in the higherredshift sub-bins, as the normalization of the MF becomessmaller by a factor 2-4, depending on mass. Despite the erta S., et al.: IR-peakers mass function 21Mass range a b [10 M ⊙ ]1.6-3.2 -0.64 ± ± ± ± ± ± Table 3.
Evolutionary coefficients describing the depen-dence of the stellar mass density on redshift and mass, inthe form of a (1 + z ) power law (see text for more details). Fig. 17.
Effect of evolutionary correction on the stellarmass function of IR-peakers. The open symbols and thedotted line belong to the mass function derived withoutevolutionary correction, while the filled symbols and thesolid line account for it.large error bars, a flattening of the stellar mass functionseems to be detected, as redshift increases. The right panelin Fig. 16 shows the evolution of the stellar mass densityas a function of redshift and for three different mass bins,highlighting that the highest-mass tail ( M ≥ × M ⊙ )of the mass function evolves less rapidly than lower-massbins.For each mass range considered, we have reproducedthe evolutionary trends with a power law: dNdM ( z ) ∝ b ( M ) × (1 + z ) a ( M ) . (7)The resulting evolutionary parameters are reported inTab. 3. Dotted lines in the right panel of Fig. 16 representthe power law fit.In any case, it is worth to point out that above z =2 . µ m sample,and the significance of the derived mass function is ratherpoor. Only by exploiting the whole SWIRE area (49 deg )will it be possible to build a catalog of massive galaxiesat z = 2 − z = 2 − z range considered, but thelow-mass bins are contributed by galaxies populating onlythe low- z part of the whole z = 2 − V a and completeness analyses, on the other hand the intrinsicevolution is not accounted for in this way. Moreover thisevolution strongly depends on mass.We therefore computed the necessary correction fac-tor for each galaxy, by comparing the average stellar massdensity over the z = 2 − µ -peakers corrected for evolution. An independent characterization of the stellar mass func-tion of massive galaxies at high redshift is the paramet-ric approach by Sandage, Tammann, Yahil (1979, STY).We focus the analysis on 5.8 µ m-peakers only, because the4.5 µ m-peaker sample turned out to be affected by incom-pleteness effects that we could not control (see Sect. 6.3).The results of this analysis will then also be used toestimate the contribution of this population to the globalstellar mass density of galaxies.As commonly found in the literature, we choose to re-produce the observed data with the usual parametric de-scription by Schechter (1976):Φ ( M ) = Φ ∗ (cid:18) MM ∗ (cid:19) α e − MM ∗ , (8)where — for the sake of clarity — M here denotes thegalaxy stellar mass, previously called M ⋆ .The three free parameters Φ ∗ , α and M ∗ represent thenormalization, the slope in the low mass regime, and thetransition mass between a power-law and the exponentialdrop-off or the e-folding mass of the latter.Since we are sampling the very high mass tail of thegalaxy stellar mass function, SWIRE IR-peakers are well suited to constrain the M ∗ parameter, while we expectto probe α only in a poor way. We will discuss this issuelater.While performing this analysis, it is very important toaccount for the main sources of uncertainty on the galaxynumber density, as described above (see Sect. 6.3). Thefirst source of uncertainty is represented by Poisson noise,since the sample consists only of 203 5.8 µ m-peakers. Thesecond is that the stellar mass estimate for each individ-ual galaxy has non-negligible uncertainty, which must beaccurately included in the parametric fit.We used a Markov-Chain Monte Carlo (MCMC) sam-pling of the parameter space, to explore the posteriorprobability function of the model, with parameters com-prising both the three Schechter parameters and also ad-ditional hyper-parameters to represent each of the indi-vidual galaxy masses.According to the Bayes’ theorem: P ( θ | d ) = P ( d | θ ) P ( θ ) P ( d ) , (9)i.e. the conditional probability of the Schechter parame-ters, given the data, is equal to the conditional probabilityof the data, given the parameters, times the prior proba-bility of the parameters, divided by a normalization factor.The theorem can be paraphrased as:posterior = likelihood × priorevidence . (10)The model’s posterior probability function (given thedata) is explored using MCMC sampling.The likelihood of the Schechter parameters is deter-mined using the standard STY method (Sandage et al.1979), derived from the luminosity function formalism.The mass function prior probability is assumed to be atop-hat (i.e. weak) prior in the Schechter parameters, withboundaries sufficiently wide as to have negligible effects.The only adopted exceptions are that M ∗ and Φ ∗ are as-sumed to be non-negative (as we are dealing with physi-cal stellar masses), while α <
0. The latter assumption ismade in order to prevent Φ( M ) deflections in the low massdomain, where no constraints are available for SWIRE IR-peakers.The hyper-parameters are constrained solely by theprior knowledge of the stellar mass probability distribu-tion function (PDF, provided by the SED fitting); there-fore their marginalization automatically accounts for themass uncertainties. The prior probability is thus computedasPrior ∝ e − χ / . (11)We use a Metropolis-Hastings algorithm as our MCMCsampler. In order for this sampler to converge efficiently,we use a combined strategy for generating proposal steps.For the Schechter parameters, we draw new steps fromunivariate Gaussian distributions, whose widths weregiven by dummy MCMC runs. The galaxy mass sampling for each galaxy is made by drawing randomly an input-sample from the set for each galaxy. We use rejection sam-pling, based on the ∆ χ of each input sample, relative tothe best fit (von Neumann 1951) to improve the efficiencyof this part of the sampling. We note that this second partof the sampling is independent of the current position inparameter space; this makes the Metropolis-Hastings ac-ceptance criterion easy to assess.Because each sample is quick to evaluate, we are ableto calculate a very large number (10 ) of steps in one hourof 3 GHz CPU time. The results of the Schechter parametric analysis of theobserved comoving number density of 5.8 µ m-peakers areshown in Fig. 18, based on the Padova-94 library.By considering the Schechter sub-space using aBayesian approach, we automatically marginalize over thehyper-parameters, producing the posterior probability dis-tributions for the mass function parameters.The top left panel in Fig. 18 reports the 1, 2, 3 σ contours in the ( M ∗ , α ) space, as computed to includethe 68.3, 95.5 and 99.7 % of the total volume occupied bythe explored samples.The median values obtained for the Schechter charac-teristic mass M ∗ are 1 . × and 1 . × [ h − M ⊙ ], in case the evolutionary correction has or has notbeen considered. The 3 σ uncertainty range is ∼ α parameter. Wefind a median value of − .
30, which in fact is reasonablysimilar to the value α = − .
35 obtained by Fontana et al.(2006) on a GOODS-MUSIC sample ranging from ∼ to ∼ × M ⊙ . If no evolutionary correction is applied,the MCMC space exploration highlights the presence ofa secondary solution at α = − .
16 (see top left panel inFig. 18). The range of suitable models extends between α = − .
12 and − . α to − .
35 andproceed to explore the ( M ∗ , Φ ∗ ) space (upper right panelin Fig. 18). The parameter Φ ∗ has a most probable medianvalue of 0 . Mpc − ], with a 1 σ range smallerthan 0 . M ∗ is similarly distributed as in theformer case with α as a free parameter.The bottom panels in Fig. 18 overlay the Schechter re-sults on the actual observed 5.8 µ m-peaker comoving num-ber density, with and without the correction for evolution-ary effects in the z = 2 − σ (dark) and 3 σ (light) ranges, asderived with our Bayesian approach. The difference be-tween the two cases is very significant, particularly as faras M ∗ and Φ ∗ are concerned. erta S., et al.: IR-peakers mass function 23 Fig. 18.
Results of Schechter parametric analysis of the observed comoving number density of 5.8 µ m-peakers (Padova-94 case). The two top panels report the distribution of models in the M ∗ , α, Φ ∗ space explored by MCMC sampling,in the case of no evolutionary correction. Contours belong to 68.3, 95.5 and 99.7% confidence levels. The bottom plotsshow the possible Schechter mass functions overlaid on the observed comoving number density of 5.8 µ m-peakers, asobtained without ( left ) and with ( right ) correction for evolutionary effects. The shaded areas represent the 1 σ and 3 σ ranges obtained in the fit. Error bars include the probability of galaxies to be shifted from one mass bin to another,due to uncertainties in the stellar mass estimate.The error bars shown in this Figure account also for thepossible shifts of galaxies from one mass bin to another,caused by the uncertainty on stellar mass, as derived fromSED fitting. The error bars are computed taking into ac-count the probability of a galaxy to be shifted, which isgiven by the actual χ distribution of all ∼ solutionsin the exploration of the SED parameter space.The same MCMC Schechter analysis has been per-formed also on the galaxy mass function obtained byadopting the Maraston (2005) SSP library. The median values of the Schechter parameters are very similar tothose obtained with the Padova-94 library.Table 4 summarizes the results, including the medianand the 1 σ , 3 σ ranges for the Schechter parameters, inboth the Padova-94 and the Maraston (2005) cases, andwith or without evolutionary correction.
7. Discussion
We have taken advantage of the wide area offered bythe SWIRE survey and the extensive multiwavelength
Fig. 19.
Comparison between the observed z = 2 − µ m-peakers and literature data, trans-formed to a Salpeter IMF (in the range 0 . −
100 M ⊙ ). Symbols represent the observed comoving galaxy number den-sities from Fontana et al. (2006), Gwyn & Hartwick (2005) and Drory et al. (2005). The two thick, long dashed curvesbelong to the Schechter STY fit to the 5.8 µ m-peakers, as obtained with (lower line) and without (upper line) the evolu-tionary correction. The thick, sort-dashed line represents the Schechter fit to z = 2 − M > M ⊙ ) galaxies at high redshift. Thus the verymassive tail of the z = 2 − z = 2 − Figure 19 compares the observed comoving numberof 5.8 µ m-peakers to z = 2 − . −
100 M ⊙ ). Only Poisson error bars are shownhere, in order to directly compare our results to literatureestimates.Our data are quite consistent with previous deriva-tions of the stellar mass function, at the bright end, butdeviate in the lower mass bins. This is mainly due to factthat no evolutionary correction is usually applied in theliterature. For comparison, we show the IR-peakers data as obtained both with (filled triangles) and without (opentriangles) this correction. The two thick long-dashed linesrepresent the two corresponding Schechter best fits. TheFontana et al. (2006) estimate (short dashes) is very sim-ilar to ours, when no evolution is taken into account.SWIRE has the advantage of probing the highestmasses in an enormous cosmic volume, approaching 10 M ⊙ , a regime where no previous studies have ever suc-ceeded because of the rarity of sources. On the other hand,SWIRE is a shallow survey and not enough information isavailable in the low mass regime. As a consequence, below M ∗ the SWIRE IR-peakers Schechter function is poorlyconstrained.A significant evolution with respect to the local galaxystellar mass function is confirmed in the highest mass bins( M > × M ⊙ ). This evolution is of the same or-der to that derived at lower masses from literature data:the number of galaxies at z = 2 − ∼ >
10 at all mass regimes, with respect to the erta S., et al.: IR-peakers mass function 25Padova-94 (no evo.)median 1 σ σ log( M ∗ ) 11.12 11.07 to 11.13 11.01 to 11.19Φ ∗ α − . − − − − σ σ log( M ∗ ) 11.22 11.18 to 11.24 11.14 to 11.27Φ ∗ α − . − − − − σ σ log( M ∗ ) 11.10 11.04 to 11.12 10.97 to 11.15Φ ∗ α − . − − − − σ σ log( M ∗ ) 11.18 11.16 to 11.23 11.12 to 11.31Φ ∗ α − . − − − − Table 4.
Results of Schechter STY analysis of the comov-ing number density of 5.8 µ m-peakers.current epoch. On the other hand, at lower redshift only aweak evolution is detected for massive galaxies ( M > M ⊙ ), at least to z ∼ . z ∼ M ⋆ ). The results of this simulation are fairly con-sistent with our 5.8 µ m-peaker data, above 2 × M ⊙ .Nevertheless, the slope of their prediction appears to betoo steep, because the number of z = 2 − M = 1 − × M ⊙ is overpredicted. Note alsothat below 10 M ⊙ , the galaxy stellar mass function issignificantly overestimated by the Millennium model, withrespect to the observed literature data.Thin lines and other shaded areas belong to semi-analytic (SAM, Bower et al. 2006; Menci et al. 2006) andhydro-dynamical (Nagamine et al. 2005a,b) models. Boththe two SAMs include AGN feedback on star formation,although in two different ways: in the Durham simulation(Bower et al. 2006) the AGN feedback is related to the gassmooth accretion onto the central black hole, continuinginto the AGN quiescent phase; while in the Menci et al.(2006) case it is produced by shock waves originated dur-ing the short active AGN phase only. Further details on the two models and their differences are found inMenci et al. (2006).The Nagamine et al. (2005a,b) hydro-dynamical sim-ulations include radiative heating and cooling, super-novae feedback, and standard star formation recipes. Thetwo models have been obtained with two different ap-proaches. In Nagamine et al. (2005a) a Smoothed ParticleHydrodynamics (SPH) entropy formulation was adopted,while the other case is based on a Eulerian mesh code witha Total Variation Diminishing (TVD) scheme. The SPHsimulation takes into account also feedback by galacticwinds and a multi-component interstellar medium influ-encing the star formation process.The agreement of these four models with the data is —unfortunately — rather poor. The two hydro-dynamicalsimulations systematically overpredict the number densityof galaxies at z = 2 −
3, at 10 < M < M ⊙ , as shownby comparison to literature data. Moreover the SPH runapplies a cutoff at M = 5 × M ⊙ , while our data alsoshow a non-null number density at these very high masses.The SAM predictions are closer to the actual data,but merely intersect the observed Φ( M ), under-predictingthe comoving number density of very massive objects andoverpredicting that of lower-mass sources. This happensat M ≃ M ⊙ for the Menci et al. (2006) model and at M ≃ × M ⊙ for Bower et al. (2006). The latter showsa very steep slope above M ∗ , too steep to be consistentwith the observational evidence. A complete view on the contribution of high-redshift mas-sive galaxies to the mass budget in the Universe is givenby their global stellar mass density, as obtained by inte-grating the stellar mass function.As far as 5.8 µ m-peakers are concerned, we have de-rived this quantity by exploiting the properties of theSchechter function. Its integral can be expressed as: ρ ⋆ ( M > M inf ) = Z ∞ M inf Φ( M ) dM = M ∗ Φ ∗ × Γ (cid:18) α + 1 , M inf M ∗ (cid:19) (12)where Γ( a, x ) is the incomplete gamma function, esti-mated from a to infinity, and M inf is the lower-mass inte-gration cutoff (set by the completeness limit, in our case).We integrate the stellar mass function over the mass rangefor which we were able to apply a reliable completenesscorrection, i.e. M > . × M ⊙ , taking into accountthe large uncertainties in the Schechter parameters, dueto poor sampling in the low mass regime. The resultingstellar mass density (accounting for evolutionary correc-tion) is ρ ⋆ = 1 . × [h M ⊙ Mpc − ], with a 3 σ rangebetween 6 . × and 2 . × [h M ⊙ Mpc − ].In the case of 4.5 µ m-peakers, instead, we can set onlya lower limit to the actual stellar mass density built in M ≥ M ⊙ . By simply integrating the observed data, Fig. 20.
Integrated stellar mass density as a function of redshift. The two large symbols represent the values of ρ ⋆ derived from our data. The filled triangle represents the contribution to the stellar mass density by 5.8 µ m-peakersabove 1 . × M ⊙ , as obtained by integrating the Schechter solutions of our MCMC fit. The filled circle sets alower limit to the stellar mass density of 4.5 µ m-peakers above 1 . × M ⊙ , obtained by integrating their observednumber density. Literature data from numerous works are shown (see right panel for references), integrated down to10 M ⊙ . Thin lines belong to different models, while the thick solid line simply connects the ρ ⋆ obtained for galaxiesheavier than 10 M ⊙ by Fontana et al. (2006) and Cole et al. (2001). All data and models have been transformed toa Salpeter IMF, extended to the mass range 0 . −
100 M ⊙ .we obtain ρ ⋆ ≥ . × [h M ⊙ Mpc − ]. Table 5 reportsthe stellar mass densities thus obtained. µ m-peak 5.8 µ m-peak ρ ⋆ [h M ⊙ Mpc − ]Observed ( V a , no evo) > . × . ± . × Schechter STY fit (no evo) – 1 . × Schechter 3 σ (no evo) – 0 . − . × Observed ( V a , evo) – 1 . ± . × Schechter STY fit (evo) – 1 . × Schechter 3 σ (evo) – 0 . − . × Table 5.
Stellar mass density for 4.5 µ m-peakers ( z = 1 − M ≥ . × M ⊙ ) and 5.8 µ m-peakers ( z = 2 − M ≥ . × M ⊙ ), as obtained with the V a andparametric analyses (Padova-94 library).Figure 20 compares the ρ ⋆ based on IR-peakers toa collection of data from the literature (Rudnick et al.2006, 2003; Fontana et al. 2006, 2004, 2003; Drory et al.2005, 2004; Franceschini et al. 2006; Bundy et al.2005; Gwyn & Hartwick 2005; Caputi et al. 2006,2005; Dickinson et al. 2003b; Brinchmann & Ellis 2000; Cole et al. 2001), all transformed to a Salpeter IMFextended between 0.1 and 100 M ⊙ . This graph highlightsa dramatic scatter in the current estimate of the stellarmass density, with significant discrepancies between thevarious authors.While the literature data were obtained by integrat-ing the mass function down to 10 M ⊙ , the IR-peakerdata points belong to sources above ∼ M ⊙ only. Thesolid thick line simply connects the ρ ⋆ estimate for galax-ies more massive than 10 M ⊙ by Fontana et al. (2006).Despite the large uncertainty on the actual value, dueto degeneracies in the Schechter fit, our measured ρ ⋆ forthe 5.8 µ m-peakers is fully consistent with Fontana et al.(2006) data. On average between 30% and 50% of thetotal stellar mass in galaxies at z = 2 − M > . × M ⊙ ) 5.8 µ m-peakers. Although very unlikely (on the basis of data byFontana et al. 2006), this value could in fact grow to ahigher fraction if we account for the uncertainties and/orcompare to the lowest estimate of the overall stellar massdensity from the literature.The thin lines in Fig. 20 represent the evolution of theglobal stellar mass density, as predicted by the same semi- erta S., et al.: IR-peakers mass function 27 analytical and hydro-dynamical models shown in Fig.19 (Bower et al. 2006; Menci et al. 2006; Nagamine et al.2005a,b). For all models the evolution of the stellar massdensity becomes steeper at z ≥ . − .
0. The values of ρ ⋆ predicted by different models gradually diverge at redshift z >
1, therefore the observational data could in principleconstrain the actual scenario, but the current wide scatterand large error bars in the stellar mass density makes itvery hard to disentangle the various models. However, weprovide the first tight constraint at the largest masses, for z = 2 − µ m-peakers andbecause they provide Schechter parameters in all redshiftbins. The dotted, dashed and solid lines represent threedifferent mass bins. Triangles represent 5.8 µ m peakers,and are obtained by integrating our Schechter fit in the10 ≤ M < × and 5 × ≤ M < M ⊙ massbins.Despite the poor statistics and the high masses andthe large “noise” previously pointed out in Fig. 20,Fontana et al. (2006) data clearly show that very massivegalaxies evolve more rapidly than lower mass objects andreach the local density at earlier epochs. It is also inter-esting to note that 5.8 µ m-peakers sample a redshift rangewhere ρ ⋆ /ρ ⋆ (0) assumes similar values in all the threemass bins considered. This effect explains why the ratioof the stellar mass function at z ∼ . φ/φ (0) spans values within a factor of 3 − × and 10 M ⊙ .At least 50-60% of the present-day stellar mass in verymassive systems seems to have assembled between z ≃ . z = 1 .
5, i.e. when the Universe was between ∼ . ∼ . z ∼
1, when more than 80% of massive galaxies were alreadyin place. SWIRE 5.8 µ m-peakers ( z = 2 −
3) represent lessthan ∼
10% of the stellar mass in the local population ofmassive galaxies ( M = 10 − M ⊙ ).According to theoretical models, the formation ofgalaxies and large scale structure occurs in the frame ofsome variant of ‘biased’ hierarchical buildup within a Λ-CDM cosmology (e.g. Cole et al. 2000; Hatton et al. 2003;Granato et al. 2004). In these scenarios the most mas-sive objects (e.g. M ⋆ > several 10 M ⊙ ) are predictedto assemble earlier, more quickly and in richer environ-ments than less massive ones (e.g. Somerville et al. 2001;Nagamine et al. 2005b).It is worth noting that no information is currentlyavailable on the environment hosting these very massivegalaxies at high-redshift. The study of the clustering prop-erties of mid-IR powerful-emitting IR-peakers has shownthat their spatial correlation length is r = 14 . ± . Fig. 21.
Top panel: ratio of the stellar mass density ata given redshift and in the local Universe, as a functionof redshift and for different mass bins. Datapoints withcrosses were computed by integrating the Fontana et al.(2006) Schechter fit to the stellar mass function.
Bottom panel: ratio of the mass function at z = 2 − µ m-peakers (triangles) and fromthe literature (other symbols). The Cole et al. (2001) localmass function was adopted.[h − Mpc] at z = 2 − M ≃ M ⊙ dark matter haloes, thus they are excellent candidates forbeing part of protoclusters at high redshifts. If this is thecase, and if less active IR-peakers (i.e. with lower mid-IRluminosities or not detected by MIPS at all) showed anal-ogous properties, they could trace quiescent stages of simi-lar systems. Broadly speaking, the studied 5.8 µ m-peakers,sampling the very massive tail of the stellar mass functionat z = 2 −
3, might represent a different class with respectto the field general population. They could, in fact, beat the center of large dark matter haloes, their compan-ions being too faint to be detected by SWIRE. Differentenvironmental conditions would rule their evolution andpartially explain the mass “downsizing” effect exemplifiedby a faster high- z evolution of massive systems.
8. Summary and conclusions
Sampling the very massive tail of the stellar mass func-tion at high redshift and estimating its contribution tothe global stellar mass density is a critical task in mod-ern cosmology, motivated by the recent evidence that a“downsizing” effect exists in the evolution of stellar massacross cosmic time.We have exploited SWIRE/Spitzer and ancillary datain the ELAIS-S1 area to perform a systematic search forhigh redshift ( z ∼ > .
0) massive (
M > M ⊙ ) galaxies.High redshift systems have been isolated by identifyingthe 1.6 µ m restframe stellar peak shifted to IRAC wave-lengths (3.6 − µ m). The availability of near-IR ( J and K s band) data allowed us to avoid low-redshift interlopersand selection aliasing due to bright 3.3 µ m PAH featuresat z ≃ .
4. A total of 203 5.8 µ m-peakers and 123 4.5 µ m-peakers have been identified over one square degree in theELAIS-S1 field.We have performed an extensive SED analysis, basedon mixed stellar population synthesis, focused on derivingthe stellar masses of the selected sample. The advantage ofnear-IR and mid-IR constraints, as well as the dependenceof results on the choice of the IMF and SSP library havebeen explored in detail. The main results are: – because of the shallow 5.8 µ m flux cut adopted, theSWIRE IR-peaker sample consists of very massivegalaxies, the majority of sources having M ⋆ > M ⊙ . Objects in the range z = 1 − M ⋆ ≃ M ⊙ , while the distribution of 5.8 µ m-peakers ( z = 2 − M ⋆ ≃ × M ⊙ . Typical uncertain-ties in stellar mass estimate (due to degeneracies in theSFH space) range between 0.1 and 0.3 dex, dependingon multiwavelength coverage. The emission of ∼
30% ofthe sources turns out to be dominated by stars olderthan 1 Gyr, and also in the majority of the remain-ing cases old stellar populations do contribute to theobserved SEDs. – the availability of mid-IR data provides a valuable con-straint on the recent star formation history of individ-ual galaxies. If no 24 µ m flux (nor upper limit) wereavailable, the resulting stellar masses could be under-estimated and the spread in mass would be wider. – despite being very useful in the selection process, near-IR ( J and K s ) data turned out to be effective in con-straining the D4000 break only in ∼
15% of cases, andpreferentially for 4.5 µ m-peakers. – the choice of a Chabrier (2003a) IMF, instead of aSalpeter (1955) one, leads to systematically lower stel-lar masses, resulting in a rigid shift of the stellar massdistribution by ∼ . – when including thermally-pulsing AGB stars in theSSP library (Maraston 2005), the fraction of objectsdominated by old stellar populations increases, but theoverall stellar mass distribution does not change signif-icantly, because of the different M ⋆ /L ratios of SSPs. The stellar mass estimates have been used to compute thecomoving number density of galaxies as a function of stel-lar mass (i.e. the observed stellar mass function), adoptingthe accessible volume formalism. Following Fontana et al.(2004), we have corrected the samples for mass incom-pleteness and recovered the mass function down to 1 . × and 1 . × M ⊙ for 4.5 µ m- and 5.8 µ m-peakersrespectively.Unfortunately, the selection, based on a 5.8 µ m flux cutturned out to be only partially sensitive to 4.5 µ m peakers,therefore only a lower limit on the mass function of z =1 − µ m-peakerswas reproduced with a parametric function, using the STY(Sandage et al. 1979) approach. The uncertainties in thestellar masses of individual sources were automatically in-cluded in the analysis by using a Bayesian formalism, andthe best fit was obtained with a MCMC sampling of theparameter space. The stellar mass function was finally in-tegrated to derive the stellar mass density locked in 5.8 µ m-peakers with M > . × M ⊙ , at z = 2 −
3. The resultsof the mass function analysis are: – the wide area surveyed by SWIRE allows the very mas-sive tail of the stellar mass function to be probed withunprecedented detail at z = 2 −
3, extending previousanalyses up to M ⋆ = 7 × [h − M ⊙ ]. – at M < × M ⊙ , a significant intrinsic evolutionhas been detected across the redshift range z = 2 − µ m-peakers number density on (1 + z ) has powersof ∼ − . ∼ − .
65 for M ∼ × and 4 × M ⊙ respectively. In the highest mass bins ( M ≥ × M ⊙ ) the number density of IR-peakers keepsnearly constant. – comparison to literature data for the z = 2 − µ m-peaker comoving number density to that of the K20,MUNICS, and GOODS-MUSIC surveys, in the highermass regime. At lower masses a significant evolution-ary correction should be applied, when the mass func-tion is averaged over a wide redshift range. – a significant evolution of the stellar mass function of M ∼ > M ⊙ galaxies with respect to the local esti-mate was detected: Φ( z = 2 − ≤ . × Φ( z = 0).Combining 5.8 µ m-peakers and literature data (e.g.Fontana et al. 2006), this implies that the bulk of mas-sive galaxies was not yet in place by the time theUniverse was ∼ ∼ . – current hydro-dynamical models significantly overesti-mate the number density of massive galaxies, while thesemi-analytic approach underestimates it. – since SWIRE 5.8 µ m-peakers sample only the very mas-sive tail of the mass function, the Schechter slope α cannot be constrained. Despite its very low signifi-cance, the best fit value α = − .
30 is similar tothat found in the literature. The other best fit pa- erta S., et al.: IR-peakers mass function 29 rameters are: M ∗ = 1 . × [h − M ⊙ ] ± ∗ = 0 . +0 . − . [ h Mpc − ]. – the integrated stellar mass density of 5.8 µ m peakers is ρ ⋆ = 1 . × [h M ⊙ Mpc − ], with a 3 σ range of ± . µ m-peak galaxies: ρ ⋆ ≥ . × [h M ⊙ Mpc − ]. – on average SWIRE massive 5.8 µ m-peakers provide30 −
50% of the total stellar mass density in galaxiesat z = 2 − – µ m-peakers provide less than ∼
10% of the stellarmass locked in massive galaxies ( M = 10 − M ⊙ )in the local Universe.The analysis carried out on SWIRE massive galaxies at z > . . × − [h Mpc − ] and it is necessary to explore large volumes ofUniverse in order to fully characterize them. The natu-ral extension of this analysis is to build the stellar massfunction of high- z galaxies over the whole SWIRE 49 deg area, taking advantage of what we have learned thanks tothe full multi-wavelength coverage in ELAIS-S1.At the time being, not much information is knownabout the environment hosting these massive high- z ob-jects. Further analyses of this population should be carriedout probing also the environmental frame they belong to,in order to correctly interpret their role in the “downsiz-ing” scenario. Acknowledgements.
We also would like to thank the wholeSWIRE team for preparing the ES1 Spitzer data. We are grate-ful to Stephan Charlot for very useful discussions about SEDfitting, IMF, and AGB stars, and Paolo Cassata on the stellarmass function and units. Finally, SB wishes to thank people atthe Sussex Astronomy Center for their warm hospitality.This work made use of observations collected at theEuropean Southern Observatory, Chile: ESO projects No.168.A-0322, 170.A-0143, 073.A-0446, 075.A-0428.The Spitzer Space Telescope is operated by the JetPropulsion Laboratory, California Institute of Technology, un-der contract with NASA. SWIRE was supported by NASAthrough the SIRTF Legacy Program under contract 1407 withthe Jet Propulsion Laboratory.
References
Abraham, R. G., Glazebrook, K., McCarthy, P. J., et al.2004, AJ, 127, 2455Baade, D., Meisenheimer, K., Iwert, O., et al. 1999, TheMessenger, 95, 15Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al.2006, AJ, 132, 1729 Berta, S., Fritz, J., Franceschini, A., Bressan, A., &Lonsdale, C. 2004, A&A, 418, 913Berta, S., Fritz, J., Franceschini, A., Bressan, A., &Pernechele, C. 2003, A&A, 403, 119Berta, S., Lonsdale, C. J., Siana, B., et al. 2007, A&A,467, 565Berta, S., Rubele, S., Franceschini, A., et al. 2006, A&A,451, 881Bertelli, G., Bressan, A., Chiosi, C., Fagotto, F., & Nasi,E. 1994, A&AS, 106, 275Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393Bolzonella, M., Miralles, J.-M., & Pell´o, R. 2000, A&A,363, 476Bower, R. G., Benson, A. J., Malbon, R., et al. 2006,MNRAS, 370, 645Bressan, A., Fagotto, F., Bertelli, G., & Chiosi, C. 1993,A&AS, 100, 647Brinchmann, J. & Ellis, R. S. 2000, ApJ, 536, L77Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000Bundy, K., Ellis, R. S., & Conselice, C. J. 2005, ApJ, 625,621Caputi, K. I., Dunlop, J. S., McLure, R. J., & Roche, N. D.2005, MNRAS, 361, 607Caputi, K. I., McLure, R. J., Dunlop, J. S., Cirasuolo, M.,& Schael, A. M. 2006, MNRAS, 366, 609Cassisi, S., Castellani, M., & Castellani, V. 1997, A&A,317, 108Chabrier, G. 2001, ApJ, 554, 1274Chabrier, G. 2002, ApJ, 567, 304Chabrier, G. 2003a, PASP, 115, 763Chabrier, G. 2003b, ApJ, 586, L133Charmandaris, V., Laurent, O., Mirabel, I. F., et al. 1999,Ap&SS, 266, 99Chiosi, C. & Carraro, G. 2002, MNRAS, 335, 335Cimatti, A., Daddi, E., Mignoli, M., et al. 2002a, A&A,381, L68Cimatti, A., Daddi, E., Renzini, A., et al. 2004, Nature,430, 184Cimatti, A., Mignoli, M., Daddi, E., et al. 2002b, A&A,392, 395Cimatti, A., Pozzetti, L., Mignoli, M., et al. 2002c, A&A,391, L1Cole, S., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2000,MNRAS, 319, 168Cole, S., Norberg, P., Baugh, C. M., et al. 2001, MNRAS,326, 255Daddi, E., Cimatti, A., Renzini, A., et al. 2004, ApJ, 600,L127Dickinson, M., Giavalisco, M., & The Goods Team. 2003a,in The Mass of Galaxies at Low and High Redshift, ed.R. Bender & A. Renzini, 324Dickinson, M., Papovich, C., Ferguson, H. C., & Budav´ari,T. 2003b, ApJ, 587, 25Drory, N., Bender, R., Feulner, G., et al. 2004, ApJ, 608,742Drory, N., Salvato, M., Gabasch, A., et al. 2005, ApJ, 619,L131Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, erta S., et al.: IR-peakers mass function 31erta S., et al.: IR-peakers mass function 31