© ESO, 2011

1. Introduction

With the exception of the cosmic microwave background, which represents the relic of the Big Bang, emitted at the last scattering surface, the extragalactic background light (EBL) is the integral of the energy radiated by all galaxies, from γ-rays to radio frequencies, across all cosmic epochs. Its energy density distribution is characterized by two primary peaks: the first at λ ~ 1 μm produced by starlight, and a second peak at  ~100–200 μm mainly due to starlight absorbed and reprocessed by dust in galaxies (e.g., Hauser & Dwek 2001). Contribution from other sources, such as active galactic nuclei, are expected to be only 5–20% to the total EBL in the mid- and far-IR (e.g., Matute et al. 2006; Jauzac et al. 2011; Draper & Ballantyne 2011). The optical and infrared backgrounds dominate the EBL by several orders of magnitude with respect to all other spectral domains.

The cosmic infrared background (CIB) was detected and measured for the first time in the mid nineties, analyzing the data obtained with the DIRBE and FIRAS instruments aboard the COsmic Background Explorer (COBE) satellite (Puget et al. 1996; Hauser et al. 1998; Fixsen et al. 1998; Lagache et al. 1999, 2000). In the global budget, the CIB provides roughly half of the total EBL (e.g., Dole et al. 2006; Hauser & Dwek 2001). Since in the local Universe the infrared output of galaxies is only a third of the emission at optical wavelengths (e.g., Soifer & Neugebauer 1991), this implies a strong evolution of infrared galaxy populations, towards an enhanced far-IR output in the past, in order to account for the total measured CIB.

The discovery of large numbers of distant sources emitting a substantial amount of their energy in the IR (e.g., Smail et al. 1997; Aussel et al. 1999; Elbaz et al. 1999; Papovich et al. 2004, among many others) demonstrated that, although locally rare, powerful IR galaxies are indeed numerous at high redshift. The deep extragalactic campaigns carried out in the nineties and 2000’s with the Infrared Space Observatory (ISO, see Genzel & Cesarsky 2000, for a summary), and the Spitzer Space Telescope (Werner et al. 2004; e.g., Papovich et al. 2004) were very efficient in identifying and characterizing large samples of mid-IR sources. In contrast, at wavelengths near the CIB peak (100–200 μm), the performance of these telescopes was strongly limited by their small apertures ( ≤ 85 cm diameter), the prohibitive confusion limits, and detectors’ sensitivity. Spitzer surveys at 70 and 160 μm produced limited samples of distant far-IR objects (e.g., Frayer et al. 2009): in the 160 μm Spitzer/MIPS band only  ~7% of the CIB was resolved into individually detected objects (Dole et al. 2004, see Lagache et al. 2005, for a review). Performing stacking of 24 μm sources, Dole et al. (2006) retrieved more than 70% of the far-IR background at 70 and 160 μm. Béthermin et al. (2010a) extended this analysis by extrapolating to very faint flux densities using a power-law fit to stacked number counts, and produced an estimate of the CIB surface brightness in agreement with direct measurements.

Launched in May 2009, Herschel (Pilbratt et al. 2010) is providing stunning results: its large 3.5 m mirror, and the high sensitivity Photodetector Array Camera and Spectrometer (PACS, performing imaging at 70, 100, 160 μm; Poglitsch et al. 2010) are tailored to overtake the confusion and blending of sources that were hampering the detection of faint far-IR sources in previous space missions. In Berta et al. (2010, hereafter called “Paper I”), we exploited data from the PACS Evolutionary Probe (PEP) survey, covering the GOODS-N, Lockman Hole and COSMOS (to partial depth only) fields, and we resolved into individual sources 45% and 52% of the CIB at 100 and 160 μm. At longer wavelengths, in the spectral domain covered by the Spectral and Photometric Imaging REceiver (SPIRE, Griffin et al. 2010), Oliver et al. (2010) resolved 15%, 10% and 6% of the CIB at 250, 350 and 500 μm. These fractions increased to 64%, 60% and 43%, respectively, using a P(D) analysis (Glenn et al. 2010). Finally, Béthermin et al. (2010b) retrieved roughly 50% of the CIB applying stacking of 24 μm sources to BLAST 250, 350 and 500 μm maps.

The PACS Evolutionary Probe (PEP¹) is one of the major Herschel Guaranteed Time (GT) extragalactic projects. It is structured as a “wedding cake” survey, based on four different layers in order to cover wide shallow areas and deep pencil-beam fields. PEP includes the most popular and widely studied extragalactic blank fields: COSMOS (2 deg²), Lockman Hole, EGS and ECDFS (450–700 arcmin²), GOODS-N and GOODS-S (~200 arcmin²). In addition, the fourth tier of the “cake” consists of ten nearby lensing clusters, offering the chance to break the PACS confusion limit thanks to gravitational lensing. An in depth description of PEP fields and survey properties is presented by Lutz et al. (2011).

Here we extend the analysis carried out in Paper I to the deepest field observed by PEP, GOODS-S, reaching a 3σ depth of 1.2 mJy at 100 μm, and including also the 70 μm PACS bandpass. The surface brightness of the CIB and its spectral energy distribution (SED) has become known with increasing detail, from its initial discovery with COBE to the Spitzer era. The real conundrum is now represented by how the CIB evolves as a function of cosmic time. This information is not only important to constrain models of galaxy evolution, but also to shed light on the intrinsic spectra of TeV sources, such as BLAZARs, whose γ-ray photons interact with the CIB (e.g., Domínguez et al. 2011; Franceschini et al. 2008; Mazin & Raue 2007). Thanks to Herschel capabilities, in Paper I we have initiated the study of the redshift build-up of the local far-IR CIB, by associating to each object detected by PACS in GOODS-N a complete UV-to-FIR SED and splitting number counts and CIB into redshift slices. Jauzac et al. (2011) performed a similar study based on stacking on 70 and 160 μm Spitzer/MIPS maps of COSMOS, and Le Floc’h et al. (2009) performed a similar analysis at 24 μm in the same field. Here we take advantage of deeper PACS observations in GOODS-S and full coverage in COSMOS to further explore the history of CIB, with the aim to understand when it was primarily emitted.

In this paper, the cosmic far-IR background is reconstructed “brick by brick” in three main “storeys”. In Sect. 2 the resolved component of number counts is presented, including slicing in redshift and a comparison to available evolutionary models. Section 3 deals with stacking of 24 μm sources onto PACS maps. The third step (Sect. 4) extends the analysis to even fainter fluxes using the “probability of deflection”, P(D), technique. After a brief digression about source confusion in Sect. 5, we discuss the observed properties of CIB in Sect. 6. Finally, Sect. 7 draws a summary of our findings.

2. Level 1: detected sources

In Paper I we studied number counts based on science demonstration phase (SDP) and early routine phase data, i.e. based on the GOODS-N, Lockman Hole and COSMOS fields. At the time of SDP, COSMOS was available only to 2/3 of its nominal depth. Altieri et al. (2010) exploited gravitational lensing in the Abell 2218 cluster of galaxies to estimate the lens amplification of the fluxes of background galaxies. In this way, they were able to push the analysis down to 1 mJy at 100 and 160 μm, thus breaking the confusion limit for blank field observations (see Sect. 5).

Here we complete the Herschel/PACS view of far-IR number counts and CIB, making use of the deepest PEP field, GOODS-S, and the full coverage of COSMOS, in addition to Paper I results. Data reduction, catalogs construction and simulations, aimed at deriving completeness, fraction of spurious sources and photometric reliability, are described in Lutz et al. (2011) and Berta et al. (2010). Table 1 summarizes the main properties of the fields taken into account here.

We applied the method described by Chary et al. (2004) and Smail et al. (1995) to correct number counts for incompleteness, on the basis of simulations. The distribution of input and output fluxes in simulations is organized in a matrix P_ij, so that i represents the input flux and j the output flux. In other words, the ijth element of the matrix gives the number of sources with ith input flux and jth output flux. The way P_ij is built implies that  ∑ _jP_ij ≤ 1 represents the completeness correction factor at the ith input flux in simulations. In order to correct the observed counts, P_ij is re-normalized such that the  ∑ _iP_ij equals the number of real sources detected in the jth flux bin. Finally, the completeness-corrected counts in the given ith bin are given by .

Figure 1 presents PACS number counts at 70, 100 and 160 μm, normalized to the Euclidean slope (i.e. the slope expected for a uniform distribution of galaxies in Euclidean space, dN/dS ∝ S^-2.5). Error bars include Poisson statistics, flux calibration uncertainties, and photometric errors. The latter have been propagated into number counts via 10⁴ realizations of random Gaussian flux errors applied to each PACS source, using a dispersion equal to the local measured noise. It is worth to note that in most cases the faint end of counts derived in shallow fields is consistent with results from deeper fields, thus confirming the validity of completeness corrections. Number counts, their corresponding uncertainties and completeness values are reported in Tables 3–5.

The left panels in Fig. 1 show the counts for each field separately. Filled symbols represent data above 80% completeness, while open symbols extend down to the 3σ detection threshold. The grey shaded areas in Fig. 1 belong to pre-Herschel number counts. We include the Rodighiero & Franceschini (2004) and Héraudeau et al. (2004) 95 μm ISO data, and the 70 μm and 160 μm Spitzer number counts by Béthermin et al. (2010a), based on GOODS/FIDEL, COSMOS and SWIRE fields. Both individual detection and stacked Spitzer counts (hatched areas) are shown. PACS counts and previous results are in good agreement, over the flux range in common. The PEP deepest field, GOODS-S, extends the knowledge on far-IR number counts one order of magnitude deeper in flux than Spitzer individual detections at 160 μm and roughly 5 times deeper at 70 μm. The 3σ limit in GOODS-S (1.2 mJy and 2.4 mJy at 100 and 160 μm, respectively) is very close to the effective depth reached by Altieri et al. (2010) in the Abell 2218 cluster when studying lensed background galaxies, but our improved source statistics provide much tighter uncertainties.

In order to provide a single reference counts description, the number counts belonging to the four PEP fields studied have been combined via a simple average, weighted by their respective uncertainties in each flux bin. Results are included in Tables 3–5, and are shown in the right panels of Fig. 1, as compared to a collection of model predictions (see Sect. 2.4).

The depth reached by PACS/PEP allows us to accurately probe the faint-end of counts. The peak in the normalized counts is well sampled and turns out to lie at  ~4 mJy at 70 μm,  ~10 mJy at 100 μm, and  ~20–30 mJy at 160 μm. The differential counts are reproduced by a broken power law (dN/dS ∝ S^α), characterized by a break at flux S_break and two distinct slopes at the at faint/bright sides of the break. A weighted least squares fit was performed on the data, and the results are presented in Table 2, for different fields and flux ranges. Breaks happen at  ~3.5 mJy at 70 μm,  ~5.0 mJy at 100 μm, and  ~9.0 at 160 μm. Uncertainties at the bright end are dominated by Poisson statistics, and nearly-Euclidean slopes are allowed.

2.1. Field to field variations

The large area probed by COSMOS (~2 deg²) can be used to test the effect of field-to-field density variations on the bright-end of number counts. We adopt a fully empirical method, describing the variance in number counts coming from the inferred variations, while a full clustering quantification goes beyond the scope of this paper. To this aim, we split the COSMOS field into a number of fully independent tiles, probing different angular scales: 28 tiles with size  ~200 arcmin², similar to the size of GOODS fields, 12 Lockman Hole like tiles (~500 arcmin²), and 4 tiles of  ~1500 arcmin² each.

Number counts in each tile were computed as described before and are plotted in Fig. 2 as thin solid lines. Black symbols represent the average counts in each flux bin, and the associated error bars are computed in the usual manner accounting for Poisson statistics, systematics and photometric uncertainties.

The properties of field-to-field standard deviation are studied in the bottom diagrams of Fig. 2. As expected, the σ(var) in number counts increases as a function of flux (upper panels), both because more luminous sources are rarer in the sky and because lower-redshift objects are dominating at these fluxes, hence probing a smaller volume. Over the whole flux range covered by this analysis, the field-to-field deviation is comparable to the uncertainty σ_P obtained combining Poisson statistics, photometric errors and systematics (bottom panels). Overall, σ(var)/σ_P is slightly larger than unity within the errors. This effect can be explained considering that neighboring sub-tiles are not fully independent: because of clustering, a non-negligible correlation term contributes to σ(var). At larger scales (dotted green lines), σ(var)/σ_P is more noisy because of the limited number of tiles available.

2.2. Ancillary data and multi-wavelength catalogs

The fields observed with PACS as part of the PEP survey benefit from a plethora of ancillary data, spanning from the X-rays to radio frequencies. We took advantage of these data to build reliable multi-wavelength catalogs and associate a full spectral energy distribution (SED) and a redshift estimate to each PACS-detected object.

As described in Paper I, a PSF-matched catalog was created in GOODS-N, including photometry from GALEX far-UV to Spitzer IRAC and MIPS 24 μm. The Southern GOODS field is rich in coverage as well. Here we adopt the PSF-matched catalog built by Grazian et al. (2006), to which we add the 24 μm photometry by Magnelli et al. (2009) and a collection of spectroscopic redshifts for more than 3000 sources (Balestra et al. 2010; Popesso et al. 2009; Santini et al. 2009; Vanzella et al. 2008; Le Fèvre et al. 2005; Mignoli et al. 2005; Doherty et al. 2005; Szokoly et al. 2004; Dickinson et al. 2004; van der Wel et al. 2004; Stanway et al. 2004b,a; Strolger et al. 2004; Bunker et al. 2003; Croom et al. 2001; Cristiani et al. 2000). Finally, we browsed the COSMOS public database² and combined the U-to-K broad- and intermediate-band photometry (Capak et al. 2007; Ilbert et al. 2009, containing 2 017 800 sources), the public IRAC catalogs, the 24 μm data (Le Floc’h et al. 2009) and the available photometric (Ilbert et al. 2009) and spectroscopic (Lilly et al. 2009; Trump et al. 2009) redshifts. As far as the Lockman Hole is concerned, an extensive ancillary catalog is currently on the make by Fotopoulou et al. (in prep.), but is not available at the time of this analysis, hence this field will not be used in this piece of analysis.

PACS catalogs were matched to the ancillary source lists by means of a multi-band maximum likelihood procedure (Sutherland & Saunders 1992), starting from the longest wavelength available (160 μm, PACS) and progressively matching 100 μm (PACS), 70 μm (PACS, GOODS-S only) and 24 μm (Spitzer/MIPS) data.

When no spectroscopic redshift is available, a photometric estimate is necessary. In GOODS-S we make use of the available photometric redshifts by Grazian et al. (2006), while new photo-z’s were produced in GOODS-N, exploiting the wealth of multi-wavelength data collected, and adopting the EAZY (Brammer et al. 2008) code. Up to 14 photometric bands were used, depending on the data available. The top panel in Fig. 3 presents the comparison between photometric and the available spectroscopic redshifts. The fraction of outliers, defined as objects having Δ(z)/(1 + z_spec) ≥ 0.2, is  ~6% over the whole sample of spectroscopic redshifts, and decreases to  ~2% for sources with a PACS detection. Most of these outliers are sources with few photometric points available, or SEDs hardly reproduced by the available templates. The median absolute deviation (MAD³) of the Δ(z) distribution⁴ is 0.040 for the whole catalog, and 0.038 for PACS-detected objects with spec-z available.

In the area covered by ancillary data, roughly 60–65% of GOODS-N sources detected by PACS in either band have a spectroscopic redshift estimate. In the GOODS-S MUSIC area (Grazian et al. 2006) this fraction is as high as  ~80%. Roughly 95% of these spectroscopic redshifts in GOODS-S lie at z < 2.0, with an almost complete coverage. For the remaining PACS sources we adopt photometric redshifts, obtaining a 100% redshift completeness above the 80% photometric completeness threshold (see Table 1). As far as COSMOS is concerned, the public photometric redshift catalog (Ilbert et al. 2009) is limited to a magnitude I ≤ 25, thus producing a redshift incompleteness in the PACS-selected sample. Only  ~75% of PACS sources have a photometric redshift estimate in the Ilbert et al. (2009) catalog. This incompleteness is independent of PACS fluxes. We derived new photometric redshifts using the EAZY code and exploiting the public COSMOS datasets. The bottom panel of Fig. 3 shows the results. The fraction of outliers is  ~1% in this case, and MAD(Δ(z)) ≃ 0.01. This result is similar to that by Ilbert et al. (2009) for the objects in common.

Given the non-null fractions of outliers in Fig. 3, it is possible that the number of sources at high redshift (e.g., z ≥ 2) in PEP catalogs is in part contaminated by “catastrophic photometric-redshift failures”, mainly represented by low-z objects with wrong, high photo-z. Because of the paucity of z ≥ 2 PACS sources benefiting from a spectroscopic follow-up, the only viable approach to test this effect is to assume that the distribution of Δ(z)/(1 + z_spec) of PACS-detected objects is similar to that of the general galaxy population. For this reason, the resulting contamination fraction has to be considered as an upper limit only.

The fraction of potential contaminants at z ≥ 2 is computed in two steps. First the fraction of z < 2 objects having z_spec < 2, but Δ(z)/(1 + z_spec) ≥ 0.2 and z_phot > 2 is derived; then this is re-scaled to the ratio of z ≷ 2 PACS sources. It turns out that up to  ~25% of GOODS-S z ≥ 2 PACS objects might have an improperly attributed high redshift. Similar results were obtained in GOODS-N, while this effect cannot be properly tested in COSMOS, because public z ≥ 2 spectroscopic redshifts are currently still lacking.

We also note that the opposite phenomenon – namely high-redshift sources with wrong photo-z potentially being shifted at low-z – does not play a significant role here, because the vast majority of z < 2 PEP/PACS objects benefits from a spectroscopic redshift measurement.

2.3. Contribution to the counts from different epochs

We exploit the rich ancillary information described in the previous Section to perform a detailed study of counts across cosmic time. The redshift derivative dN/dz of PACS sources above the 80% photometric completeness limits, normalized to 1 deg², is shown in Fig. 4 for the three fields and bands considered here. The covered flux ranges are quoted on each panel. In the COSMOS area, we sample the bright end of PACS counts. The distribution in this field peaks at z ≤ 0.5 in both 100 and 160 μm bands, with 70–80% of all sources lying below z = 1 and  ~20% between z = 1 and z = 2. For the deeper GOODS-N data, the peak of the redshift distribution shifts to z ~ 0.7. Although the small sampled area limits source statistics, objects at high redshift (up to z ≃ 5) start to pop up. Our deepest field, GOODS-S, covered in all three PACS bands, displays some remarkable features. Overall dN/dz is now peaked around z ~ 1. It is possible to recognize two well known structures at z ≃ 0.7 and z ≃ 1.1 (e.g., Gilli et al. 2003; Vanzella et al. 2005), which produce narrow and intense spikes in the distribution at 100 and 160 μm. On the other hand, at 70 μm these structures are barely seen. At higher redshift, a broad “bump” is detected, between z = 2 − 3. This peak cannot be identified in the shallower 70 μm data, but is outstanding at the other two PACS wavelengths. Cutting the GOODS-S catalog at the depth reached by GOODS-N (5.5 mJy and 11.0 mJy at 100 and 160 μm, respectively), the high-z feature disappears and the redshift distribution resembles that of GOODS-N, with only a few sources left above z ≥ 2. Similarly, when cutting GOODS-S at the COSMOS 80% depth, we retrieve a distribution peaked at z = 0−0.5, obviously with much poorer statistics than in the COSMOS field itself. An extensive analysis of PACS GOODS-S large scale structure at z = 2−3 and of a z = 2.2 filamentary over density is being presented by Magliocchetti et al. (2011). Table 6 reports dN/dz  [deg^-2], as derived in these three PEP fields. The median redshifts of the sources detected in GOODS-S are z ~ 0.6,  ~ 0.7,  ~ 0.8 in the three PACS bands, but would shift if the known spikes in dN/dz were not present.

Finally, it is possible to split number counts in GOODS-S, GOODS-N and COSMOS into redshift bins, similarly to what was done in Paper I for GOODS-N. Counts are then constructed in four broad bins, in order to allow for a sufficient number of sources in the small fields. In every redshift bin we apply the same completeness correction derived for the total number counts at the given flux. Similarly to what was done for the total counts, we build combined number counts through a weighted mean in each flux bin. The average number counts, sliced in redshift intervals are reported in Tables 7–9, and are shown in Fig. 5. PACS observations in the three fields nicely complement each other: the southern GOODS field reaches deep flux densities, not probed previously, and also traces the high redshift population that was missing in GOODS-N. The bright end of counts is sampled by observations in the COSMOS area, which not only contribute to the low-redshift counts, but also nicely trace the bright component up to z = 2 and beyond.

2.4. Comparison to model predictions for counts across time

Number counts and redshift distributions encode the evolution of galaxies: the upturn at intermediate fluxes and the over-Euclidean slope at the bright side of the peak are usually interpreted as signatures of strong evolution in the properties of the underlying galaxy populations, and have stimulated a plethora of model interpretations.

We are comparing here our results to examples of two basic classes of models attempting to reproduce these observables. Backward evolutionary models transform the statistical properties of far-IR galaxies observed at known redshift (e.g., local luminosity functions) into observables at any redshift (e.g., galaxy counts, CIB brightness, etc.) assuming a library of template SEDs and parametric laws of luminosity and/or density evolution. These models do not implement fundamental physical information but simply attempt to describe the evolution of galaxy populations (e.g., Lagache et al. 2004; Franceschini et al. 2010; Rowan-Robinson 2009; Le Borgne et al. 2009; Valiante et al. 2009; Béthermin et al. 2011; Marsden et al. 2010; Gruppioni et al. 2011).

Very simple in their principles, backward models are thus parameterizations embedding the knowledge obtained from previous observations. The different adopted recipes were optimized to reproduce available observables such as mid-IR ISO and Spitzer number counts and far-IR Spitzer counts based on detections and – in some cases – stacking, sub-mm counts, as well as redshift distributions. A couple of models took advantage of early Herschel results, namely SDP PACS and SPIRE number counts (Béthermin et al. 2011) and PACS luminosity functions up to z ~ 3 (Gruppioni et al. 2011). Generally speaking, confronting backward models to the new, detailed Herschel data is a test to their flexibility in producing reliable predictions at far-IR wavelengths.

At the opposite side, forward evolution models simulate the physics of galaxy formation and evolution forward in time, from the Big Bang to present days. Current implementations are based on semi-analytic recipes (SAM) to describe the dissipative and non-dissipative processes influencing galaxy evolution, framed into Λ-CDM dark matter numerical simulations (e.g., Lacey et al. 2010; Niemi et al., in prep.). Galaxy radiation, including far-IR emission, is computed with spectrophotometric synthesis and radiation transfer dust reprocessing, given the fundamental properties of the galaxies in the model.

Figures 1 and 5 show a collection of models overlaid to the observed PACS counts. Most of the models reproduce fairly well the observed 100 and 160 μm normalized counts, on average the most successful being Franceschini et al. (2010), Marsden et al. (2010), Rowan-Robinson (2009), Valiante et al. (2009) and Gruppioni et al. (2011). Very different assumptions produce relatively similar number counts predictions. The Franceschini et al. (2010) model adopts four different galaxy populations, including normal galaxies, luminous infrared galaxies, AGNs and a class of strongly-evolving ULIRGs, which dominate above z ≥ 1.5, but is negligible at later epochs, resembling high-z sub-mm galaxies. Also Rowan-Robinson (2009) uses four galaxy populations (cirrus-dominated quiescent galaxies, M82-like starbursts, Arp220-like extreme starbursts, AGN dust torii), but employs analytic evolutionary functions without discontinuities. Valiante et al. (2009) describe galaxy population taking into account the observed local dispersion in dust temperature, the local observed distribution of AGN contribution to L_TIR as a function of luminosity, and an evolution of the local luminosity-temperature relation for IR galaxies. Marsden et al. (2010) base their SEDs on Draine & Li (2007) prescriptions and tune them to reproduce the local color-dependent far-IR luminosity functions; they include luminosity, density and color evolutions to fit sub-mm counts, redshift distributions and EBL. Finally, Gruppioni et al. (2011) include a significant Seyfert-2 population, based on a fit to Herschel LFs (see also Gruppioni et al. 2010).

Semi-analytical approaches surely represent a more complete view of galaxy evolution, including a wide variety of physics in a single coherent model. They cover a wide range of observational data: UV, optical, near-IR luminosity functions, galaxy sizes, metallicity, etc. Moreover, at far-IR wavelengths, even in the case that global properties such as star formation and AGN activity are correctly modeled, a further complication arises from the assumptions about dust content and structure, which need to be invoked. As a result, the large number of parameters involved goes at the expense of inference precision: the performance of SAM models with respect to PACS observables needs still substantial tuning.

As described above, it seems that pre-Herschel models, with only a few exceptions, are quite successful at reproducing total number counts despite the range and diversity of the employed solutions, thus showing that the discriminatory power of this observable is rather limited.

The right answer comes from the redshift information available in the selected PEP fields. Figure 4 presents the redshift distribution dN/dz of PACS galaxies, and Fig. 5 shows number counts split in redshift slices (d ²N/dS/dz). The combination of the two is a real discriminant: model predictions are now dramatically different.

None of the available models provides a convincing coherent prediction of the whole new set of observables (both dN/dz and d ²N/dS/dz), over the flux range covered. One common source of discrepancy seems to be the well-know degeneracy between dust temperature (and hence luminosity and total dust mass) and redshift (e.g., Blain et al. 2003): red objects can be either cool local galaxies or warmer distant galaxies. What is observed is a mis-prediction of the redshift distribution, reflected for example in an overestimation of high redshift counts and an underestimation at later epochs – or vice versa. Several models seem to systematically over predict the number of galaxies above z ≃ 1 in the deep regime. A variety of SED libraries are adopted in these model recipes, based on different assumptions and templates. These differences in SED shapes and their evolution – as well as the implementation of luminosity and density evolution for the adopted galaxy populations – indeed produce significantly different results.

Among all, the Béthermin et al. (2011) distribution is probably the closest to observations, although it presents a significant underestimation of dN/dz at 70 μm. This model was optimized taking into account differential number counts between 15 μm and 1.1 mm, including PACS (Berta et al. 2010) and SPIRE (Oliver et al. 2010) early results, mid-IR luminosity functions (LF) up to z = 2, the far-IR local LF, and CIB measurements. This success demonstrates the need to include in model tuning (ideally by means of a proper automated fit) not only Herschel data, but also the detailed redshift information (e.g., evolving LF, number counts split in redshift bins, redshift distributions, etc.).

In summary, while several pre-Herschel backward evolutionary models already provide reasonable descriptions of the total counts in at least some of the PACS bands, they tend to fail in the synopsis of all counts and in particular redshift distributions. The new Herschel dataset is a treasure box, allowing to explore the evolution of counts and CIB with unprecedented detail and much deeper than previous data, on which models were calibrated. Significant modifications on model assumptions for SEDs and/or evolution will be needed for a satisfactory fit to this new quality of data.

3. Level 2: stacking of 24 μm sources

Limiting the number counts analysis to individually-detected sources only, we miss a significant fraction of the information stored in PACS maps. It is possible to recover part of this information by performing stacking of sources from deeper data (typically at shorter wavelengths, e.g., Dole et al. 2006; Marsden et al. 2009; Béthermin et al. 2010a). The 24 μm catalog in the GOODS-S field, extending down to  ~20 μJy (Magnelli et al. 2009), provides the ideal priors to perform stacking of faint sources on PACS maps.

The aim here is to transform 24 μm faint number counts into PACS counts, given the average PACS/24μm colors of galaxies. The overall procedure consists in building such colors as a function of 24 μm flux, including both sources detected by PACS and un-detected ones. Then the derived colors are used to transform 24 μm counts into PACS counts, probing the faint end of the distribution, below the PACS detection limit.

To reach this goal, we select 24 μm sources not detected in PACS maps, and bin them by their 24 μm flux. Following the standard technique described by Dole et al. (2006) and Béthermin et al. (2010a), in each S(24 μm) flux bin, we pile-up postage stamps at the position of these sources and produce stacked frames at 70, 100 and 160 μm. We measure the flux density of the stacked sources by performing PSF-fitting in the same way as for the individually-detected objects (see Berta et al. 2010). Uncertainties on the stacked fluxes are computed through a simple bootstrap procedure. Stacking was performed on PACS residual images, i.e. after removal of individually-detected sources, and the fluxes of the latter were then added back to stacking results. Finally, we checked that no significant signal was detected when stacking at random positions, and that stacks of PACS-detected sources retrieved their actual total summed flux. Flux corrections due to losses during the high-pass filtering process were tested via simulations and turned out to have a negligible effect on our results (see Lutz et al. 2011 for a description).

We derived average flux densities simply dividing by the total number of sources in each 24 μm flux bin, and then we built PACS/24 μm average colors. Figure 6 shows the results. Stacking of mid-IR sources actually allows us to retrieve missed PACS fluxes down to faint regimes, thus recovering the actual average colors, that would not be possible to derive otherwise. On the bright end, average stacked+detected fluxes and individual detections progressively converge, and become consistent within errors when  ≥ 55% of 24 μm sources in the given flux bin are detected by PACS.

Stacked points were fit with a polynomial function and then used to transform the observed 24 μm differential counts into PACS counts, following the recipe suggested by Béthermin et al. (2010a): (1)where the right terms are computed in the given 24 μm flux bin and the left term is computed at the corresponding S_FIR = f(S₂₄), as defined by the average color functions. The derivative dS₂₄/dS_FIR is estimated numerically.

Although this method has been widely and successfully used by several authors, it is worth to note that stacking results are affected by completeness limitations. In this specific case, a flux cut in the 24 μm priors (e.g., the Magnelli et al. 2009 catalog reaches 20 μJy at the 3σ detection threshold) induces incompleteness in the PACS stacked fluxes, because red objects (faint at 24 μm) are not included in the stack. Consequently, stacking of 24 μm sources on the PACS maps provides only a lower limit to the PACS number counts below a given PACS flux. The bottom panel of Fig. 6 is based on Le Borgne et al. (2009) mock catalogs, and shows that, for a 24 μm cut of 20 μJy, the PACS stacked counts are 90% complete only above  ~0.35,  ~0.7 and  ~1.3 mJy at 70, 100 and 160 μm, respectively. To this effect, one should add the intrinsic properties of the 24 μm parent catalog, which reaches  ~90% completeness at S(24) ~ 35 μJy.

Figure 7 shows the differential number counts dN/dS and includes the results of stacking above the 80% stacking completeness (big blue crosses), in good agreement with the faint-end of the resolved counts presented before (open symbols in Figs. 7 and 1). The result of a power-law fit is reported in Table 2: at 70 and 100 μm these slopes are slightly flatter, while at 160 μm slightly steeper, than those of resolved counts.

4. Level 3: P(D) analysis

It is possible to extend the study of number counts and CIB to even deeper flux regimes, beyond the confusion limit, exploiting the so-called “probability of deflection” statistics, or P(D) distribution (e.g., Scheuer & Ryle 1957; Condon 1974; Franceschini et al. 1989, 2010; Oliver et al. 1997; Patanchon et al. 2009; Glenn et al. 2010). In a simplified notation, this is a representation of the distribution of pixel values in a map: its shape and width are mainly driven by three components: number counts, the instrumental spatial response function, and the instrumental noise. Once the three are known, it is possible to reproduce the observed pixel flux probability distribution, or – vice versa – the observed P(D) can be used to constrain the underlying dN/dS number counts.

The first P(D) applications were carried out at radio frequencies (e.g., Scheuer & Ryle 1957; Condon 1974; Franceschini et al. 1989), but have subsequently been applied to different wavelengths, spanning across the whole electromagnetic spectrum (e.g., X-ray, Scheuer et al. 1974, Barcons et al. 1994; sub-mm, Hughes et al. 1998, Weiß et al. 2009; Scott et al. 2010 Patanchon et al. 2009; infrared, Oliver et al. 1997). Since the Herschel satellite was launched, a number of analyses have made use of P(D) techniques to describe the properties of real maps (Glenn et al. 2010; Oliver et al. 2010), as well as produce model predictions to be compared to actual far-IR observations (Béthermin et al. 2011; Franceschini et al. 2010).

We performed the P(D) analysis on GOODS-S PACS maps, modeling differential counts dN/dS as a broken power-law with three sections and two nodes. To this aim, we developed a new, numeric method to estimate the P(D), given the counts model, the observed PSF and the observed noise. The position of nodes and amplitude of counts at these positions are free parameters. The best combination of parameters reproducing the observed P(D) is sought by minimizing the difference between model prediction and data, via a Monte Carlo Markov Chain (MCMC) engine. Appendix A describes the details of this method.

Figure 8 shows the fit to the observed P(D) in the three PACS bands; Fig. 7 and Table 10 report on the number counts model producing the best fit to the observed P(D). Note that no constraint was set to number counts, but only the observed P(D) distribution was fit. It is remarkable to note the very good agreement between P(D) results and the resolved number counts at the bright end, and over the whole flux range covered by the latter. At the very faint end, the number count model derived from P(D) tends to diverge, because of degeneracies. This is dramatically evident at 70 μm, where the faintest section in dN/dS is totally undefined. For this reason, we limit the 70 μm result to  ~0.35 mJy, when computing CIB surface brightnesses in the next Sections. At 100 and 160 μm, the P(D) approach allows the knowledge of PACS number counts to be extended down to  ~0.2 mJy, one order of magnitude deeper than individual detections. The resulting slopes are consistent with resolved and stacked counts, over the flux range in common.

5. Confusion noise

The extreme depth reached with the P(D) analysis, and the detailed view of number counts coming from individually-detected sources, make worthwhile to take a digression and study the properties of confusion noise in PACS maps.

Following the Dole et al. (2003) formalism, two different approaches can be adopted to describe the confusion noise due to extragalactic sources, i.e. the effect of fluctuations due to the presence of point sources in the beam (see also Frayer et al. 2006; Condon 1974; Franceschini et al. 1989, among others): the photometric and the density criteria.

Given a detection threshold S_lim, the former is derived from the signal fluctuations due to sources below S_lim, while the latter consists in deriving the density of sources detected above S_lim and limiting the probability that objects are missed because blended with neighbors to a small fraction (e.g., 10%).

Here we apply both criteria to PACS/PEP data and directly measure the confusion noise and limits of Herschel deep extragalactic imaging.

5.1. Photometric confusion

Photometric confusion is defined by setting the desired signal-to-noise ratio q between the limiting flux S_lim and the noise σ_c describing beam-to-beam fluctuations by sources fainter than S_lim, i.e. q = S_lim/σ_c. We adopt a value q = 5 for our analysis.

Generally speaking, the noise measurable on PACS images is given by the combination of the instrumental noise σ_I (including photon noise, detector noise and data-processing noise) and σ_c. The instrumental noise was estimated empirically, similarly to Frayer et al. (2006). We build partial-depth maps in the GOODS-S field, exploiting only a fraction of the available PACS observations. By progressively increasing total exposure time, we drew the diagrams in Fig. 9. Sources were detected above S_lim in the usual way (see Berta et al. 2010, and Lutz et al. 2011) and subtracted prior to measuring the noise on the maps. For short exposure times, noise is dominated by instrumental effects, and is proportional to t^-0.5. This allows to estimate σ_I and extrapolate it to deeper regimes. When a long effective exposure time is used, the measured total noise σ_T deviates from the t^-0.5 trend, testifying that confusion noise plays an increasingly significant role. Since both σ_T (after source extraction) and σ_I have nearly Gaussian distributions, the confusion noise can also be approximated by . We iterated between source extraction at different depths and noise measurements, until convergence at q = 5 was reached.

No deviation from the instrumental σ_I ∝ t^-0.5 is detected at 70 μm, at the depth of PEP GOODS-S observations, while we derive σ_c = 0.27 and 0.92 mJy at 100 μm and 160 μm, respectively (with q = 5). Table 11 summarizes the derived values.

5.2. Density criterion for confusion

The second criterion to determine the influence of confusion on extragalactic observations is defined by requiring a minimum completeness level in the detection of sources brighter than S_lim, driven by the fraction of objects missed because blended to their (detected) neighbors. Source number counts provide a neat one-to-one relationship between source density and flux, so that the flux at which the source density confusion (SDC) limit is reached is univocally defined, given dN/dS.

Adopting the Lagache et al. (2003) definition of beam (i.e., ), and allowing a probability P = 10% that sources are tightly blended and thus cannot be extracted, Dole et al. (2003) derived a SDC limit of 16.7 beams/source. In what follows we use this value in order to facilitate comparisons with previous works and future exploitations of this PACS confusion limit estimate. Nevertheless, it is important to recall that such threshold depends on the steepness of source counts. Following the general treatment of counts and confusion by Franceschini (1982), the density of sources with flux S > S_lim can be expressed as , where α is the slope of number counts, and q was defined in Sect. 5.1. The canonical density of 16.7 beams/source is retrieved for α =  −1.8, which lies within the range of slopes estimated by power-law fitting of PEP data (see Table 2), but the SDC limit can significantly vary between 8.3 and 25 beams/source for α =  −1.5 to  − 2.0.

The depth reached in PACS maps, especially in the GOODS fields, is such that the high density of detected sources hinders the extraction of fainter objects. In Paper I we have already treated the case of GOODS-N, showing that the observed catalog is already hitting the density confusion limit at 160 μm.

At 70 μm, despite deep observations, GOODS-S is far from being confused: even at the 3σ limit of 1.1 mJy, PACS catalogs have a density of 32 beams/source including completeness correction. The SED shape in the far-IR is such that, while a large number of sources is detected at 100 and 160 μm at the PEP effective exposure time, positive k-correction strongly suppresses the number of detections at 70 μm. The typical 70/160 color of galaxies is red enough to require a deeper threshold in order to detect a comparable number of sources in the two bands, but PACS sensitivity at 70 μm is not sufficient to compensate this effect. Based on P(D) results, the 16.7 beams/source SDC limit would be reached at  ~0.4 mJy in this band.

At longer wavelengths, GOODS-N and GOODS-S completeness-corrected number counts and P(D) modeling agree in indicating that the SDC limit is reached between 1.5 and 2.0 mJy at 100 μm and around 8.0 mJy at 160 μm (Table 11). At the GOODS-S 3σ levels (1.1 mJy at 100 μm and 2.4 mJy at 160 μm) the density of sources in our catalogs is 10 and 5 beams/source, after completeness correction. It is worth to note that the Dole et al. (2003) SDC limit corresponds to a blending-completeness of 90%, while sources can be reliably extracted at lower completeness also below the 16.7 beams/source threshold. Furthermore, photometric completeness further flattens the detected counts, thus allowing for sporadic very faint objects to be extracted. All values quoted above here are nevertheless based on completeness-corrected number counts. For comparison, the depth reached by P(D) at 100 and 160 μm digs deep into density confusion, down to levels corresponding to 2–4 beams/source.

6. Discussion: the cosmic infrared background

Since its discovery by COBE (Puget et al. 1996; Hauser et al. 1998), several attempts have been carried out to derive the surface brightness of the CIB, based on direct measurements, integration of number counts, statistical analyses, and γ − γ opacity to TeV photons. The PACS data exploited so far can now be used to derive new, stringent lower limits to the actual far-IR background.

6.1. CIB measurements: state of the art

The integral of number counts provides a lower limit to the cosmic IR background, to be compared to direct measurements from COBE maps. Here we adopt the DIRBE measurements by Dole et al. (2006): 14.4 ± 6.3, 12.0 ± 6.9, and 12.3 ± 2.5 [nW m^-2 sr^-1], at 100, 140, and 240 μm, respectively. Performing the geometrical average between the upper and lower limits to the CIB given by the Fixsen et al. (1998) fit to the FIRES spectrum, one obtains a 160 μm CIB brightness estimate of 12.8 ± 6.4 [nW m^-2 sr^-1], consistent with Dole et al. (2006) measurements at 140 and 240 μm. Dole et al. (2006) thoroughly describe the calibration and uncertainty details for CIB direct measurements. The most significant problems arise from the effective brightness of zodiacal light, currently known only within a factor of two accuracy.

Odegard et al. (2007) used the Wisconsin Hα Mapper (WHAM) Northern Sky Survey as a tracer of the ionized medium, in order to study the effects of the foreground interplanetary and Galactic dust on DIRBE CIB measurements. The authors find CIB surface brightness values in agreement with those by Dole et al. (2006), once renormalized to the FIRAS photometric scale. Juvela et al. (2009) derive the CIB surface brightness at 170 μm from ISO maps, obtaining a value  ~1.5 times higher than COBE results, with 30% systematics and 30% statistical uncertainties. Based on recent observations of the AKARI Deep Field South (ADF-S), carried out with this satellite at 65, 90, 140 and 160 μm, Matsuura et al. (2010) detected and measured the absolute brightness and spatial fluctuations of the CIB, deriving values consistent to our Dole et al. (2006) 160 μm reference, but detecting a significant excess at 90–140 μm.

Finally, the most recent estimate of the far-IR CIB comes from the fluctuation analysis carried out by Penin et al. (2011) on Spitzer 160 μm maps of the SWIRE ELAIS-N1 field and reprocessed IRAS 60 and 100 μm data (Miville-Deschênes et al. 2002). The newly-derived CIB brightness is 6.6 ± 2.7 and 14.4 ± 2.3 [nW m^-2 sr^-1], at 100 and 160 μm, respectively. The 160 μm measurement is consistent with the previous COBE interpolation, but the 100 μm value is almost halved.

No direct measurement is available at 70 μm, except for the Miville-Deschênes et al. (2002) fluctuation analysis on the IRAS 60 μm map (νI_ν = 9.0 [nW m^-2 sr^-1]) and the AKARI estimate at 65 μm (νI_ν = 12.4 ± 1.4 ± 9.2 [nW m^-2 sr^-1], including statistical and zodiacal-light uncertainties, Matsuura et al. 2010).

Further constraints on the actual value of the CIB come from the cosmic photon-photon opacity: very high-energy photons suffer opacity effects by γ − γ interactions with local radiation backgrounds, producing particle pairs (e.g., Franceschini et al. 2008; Nikishov 1962). Because of the large photon density of the cosmic microwave background (CMB), any photon with energy ϵ > 100 TeV has a very short free path, and extragalactic sources are undetectable above this energy. Less energetic photons suffer from the opacity induced by extragalactic backgrounds, other than the CMB, including the CIB. Observations of absorption features and cut-offs in the high-energy spectra of BLAZARs have been successfully used to pose upper limits to the intensity of the EBL. Mazin & Raue (2007), among many others, assembled a compilation of 11 TeV BLAZARs (all those known at the time of their analysis) and combined them obtaining constraints to the EBL over the whole 0.44–80 μm wavelength range.

6.2. The total CIB in PEP

Integration of number counts was performed over as wide a flux range as possible, i.e. using the combined counts presented in Sect. 2, including all the four fields analyzed so far: GOODS-S, GOODS-N, Lockman Hole and COSMOS. To this aim, GOODS-S data were extended down to the 3σ detection threshold, thus reaching 1.2 mJy at 70 and 100 μm and 2.0 mJy at 160 μm. At the bright side, COSMOS allows the integration to be carried out to 140 mJy at 100 μm and 360 mJy at 160 μm. Including completeness corrections, the derived CIB locked in resolved number counts is νI_ν = 7.82 ± 0.94 and 9.17 ± 0.59 [nW m^-2 sr^-1], corresponding to 54 ± 7% and 72 ± 5% of the COBE direct measurements at 100 and 160 μm, respectively. If the Penin et al. (2011) 100 μm value were taken as reference, then PACS sources would produce  ~100% of the CIB. The resolved CIB at 70 μm is νI_ν = 3.61 ± 1.12 [nW m^-2 sr^-1].

The PEP fields are missing the very bright end of number counts, that can be probed only over areas even larger than COSMOS. This lack of information is particularly relevant at 70 μm, because this band was employed only for GOODS-S observations. We therefore add previous bright counts measurements to our data, in order to derive the CIB brightness emitted above PEP flux limits, integrated all the way to  + ∞. The available data taken into account here are the 70 and 160 μm Spitzer counts by Béthermin et al. (2010a), built over a more than of  ~50 deg² and extending to  ~1 Jy, and 100 μm ISO and IRAS number counts (Rodighiero & Franceschini 2004; Héraudeau et al. 2004; Rowan-Robinson et al. 2004; Efstathiou et al. 2000; Oliver et al. 1992; Bertin et al. 1997), extending to  ~60 Jy. Beyond the flux range covered by these past surveys, we extrapolate the number counts with an Euclidean law, dN/dS ∝ S^-2.5, normalized so to match the very bright end of observed counts, and extended to infinity.

Figure 10 shows the cumulative CIB surface brightness as a function of flux, as derived from the combination of PACS, Spitzer, ISO and IRAS data. The curve describing the resolved CIB (red circles) rapidly converges to the COBE measurements at 100 and 160 μm, over the flux range covered by PACS/PEP. Our data extend the knowledge of counts and CIB by over an order of magnitude in flux, with respect to previous estimates based on individual detections. The depth reached in GOODS-S is similar to the one obtained in Abell 2218 through gravitational lensing, and the corresponding CIB surface brightnesses are consistent within the uncertainties. The total CIB emitted above PEP 3σ flux limits (1.1 mJy at 70 μm, 1.2 mJy at 100 μm, 2.0 mJy at 160 μm), is νI_ν = 4.52 ± 1.18, 8.35 ± 0.95 and 9.49 ± 0.59 [nW m^-2 sr^-1] at 70, 100 and 160 μm, respectively. Table 12 lists all the values obtained.

Roughly 10% of the CIB down to the 70 μm adopted threshold is emitted by bright galaxies, out of reach for our survey, due to the limited volume sampled. On the other hand, only a few percent is produced at fluxes brighter than PEP COSMOS upper limit at 100 and 160 μm. Including 100 μm ELAIS and IRAS data (up to 60 Jy), only a negligible fraction (<0.05%) is missed at the bright end, while beyond Spitzer counts (reaching  ~1 Jy at 70 and 160 μm) the Euclidean extrapolation provides roughly 1–2% of the total CIB.

The power-law fit to resolved counts (see Sect. 2 and Table 2) provides a new estimate of the total expected CIB. We extrapolate power laws down to 1.0 μJy and obtain νI_ν ≥ 11.09, and  [nW m^-2 sr^-1] in the three PACS bands. At 70 μm our data provide only a lower limit, because the curve obtained from our power-law fit is not fully converging at 1 μJy (see Fig. 10). The uncertainty on the 100 μm CIB is still large because discordant slopes were found at the faint end in the different PEP fields. The results are fully consistent with COBE data, but PEP and PACS pinpoint the total CIB values with unprecedented precision, thanks to the high quality of observations, survey strategy, maps, and – last but not least – thanks to Herschel capabilities (grey shaded areas in Fig. 10 denote the 3σ uncertainty on power-law fits).

Stacking results are fully consistent with resolved number counts, and the P(D) analysis extends down by another decade in flux with the exception of 70 μm, for which we truncated the computation of CIB because of a strong divergence in P(D) uncertainties (see Sect. 4). Consequently, the lower limit set to the CIB by PACS data is further improved: above the flux level reached by P(D) statistics, the contributions to the CIB are 4.98, 9.32 and 11.31 [nW m^-2 sr^-1] in the three bands, thus recovering  ~ 65% and 89% of the total Dole et al. (2006) CIB at 100 and 160 μm, respectively. When referred to the total values obtained through power-law extrapolations, the P(D) analysis recovers 45%, 74% and 83% of the total background at 70, 100 and 160 μm.

6.3. CIB contributions from different cosmic epochs

Thanks to the rich ancillary datasets available in PEP fields, it is possible to estimate the amount of CIB emitted at different epochs. In order to reach this goal, we integrate the number counts split into redshift bins. The same completeness correction was applied in each redshift bin, as derived from total counts. Table 13 reports the results for GOODS-S and for the combination of GOODS-S, GOODS-N and COSMOS.

Figure 11 shows the fraction of resolved CIB emitted at different cosmic epochs in the GOODS-S and COSMOS fields. In this case, we keep the two fields separate, in order to study the details of different flux regimes and aid the comparison to model predictions. As stated above, deep observations resolve most of the CIB at 100 and 160 μm, and thus give a fairly complete census of the redshift dependence of the background. At the PEP 3σ detection threshold, more than half of the resolved CIB is emitted by objects lying at z ≤ 1. The new results are in line with our Paper I analysis, based on shallower GOODS-N detections and stacking of 24 μm sources. As expected, at the depth of GOODS-S, galaxies in the highest redshift bin come into play and provide a non negligible contribution to the CIB, as high as  ~15% of the resolved amount in GOODS-S.

For comparison, based on SED template extrapolations from 15 μm, Elbaz et al. (2002) estimated that the contribution of ISO 15 μm sources to the 140 μm CIB was of the order of 65%, with a median redshift of z = 0.6, and mostly emitted below z = 1. Exploiting Spitzer 24 μm observations in the COSMOS field, Le Floc’h et al. (2009) showed that  ~50% of the 24 μm background intensity originates at z ≃ 1.

The comoving (IR) luminosity (or SFR) density dependence on redshift (Madau et al. 1998, see Gruppioni et al. 2010 for a recent determination based on Herschel data) is known to peak between redshift z = 1.5 and 3.0. Harwit (1999) showed that when transforming the resulting energy generation rate into unit redshift interval, the high redshift component is strongly suppressed, by a factor (1 + z)^5/2 in a q₀ = 0.5 cosmology. Moreover, the effective energy reception rate, thus the background observed locally, is further suppressed by a factor (1 + z), because bandwidths are reduced. Based on these and other simple considerations, it is easy to demonstrate that the distribution of the integrated CIB should be dominated by z < 1 galaxies (Harwit 1999).

The relative fraction among the four bins moves from low to high redshift, as wavelength increases: Fig. 12 shows the cumulative CIB fraction as a function of redshift, computed on the combination of GOODS-S and COSMOS, and including completeness correction. Half of the CIB detected over the flux ranges quoted on each panel was emitted below redshift z = 0.58, 0.67 and 0.73 at 70, 100, 160 μm, respectively.

While at 70 μm  ~ 80% of the resolved CIB was emitted at z ≤ 1.0, at 160 μm this fraction decreases to  ~55%. This trend is mainly due to the positive k-correction of far-IR SEDs short ward of their peak, and will become even more relevant in SPIRE and sub-mm bands, where negative k-correction comes into play (Oliver et al. 2010; Marsden et al. 2009). Note also that including the contribution from 70 μm unresolved galaxies would shift the distribution to slightly higher redshifts (see also discussion about 70 μm depth and confusion in Sect. 5.2). On the other hand, adding the bright end of number counts, sampled by the observations in the COSMOS field, does not significantly influence the high redshift bins in the CIB, because most of the contribution goes into filling the lower redshift slices (see Table 13). The 90th percentiles of the redshift-cumulative CIB fall at z = 1.38, 2.03 and 2.20. Excluding the bright end covered by COSMOS (only at 100 and 160 μm), we witness only a marginal shift of the 90%-light redshift by Δ(z) ≤ 0.05.

Performing SED fitting of each individual object detected in GOODS-S in the PACS bands, Rodighiero et al. (2010) derived IR luminosities (8–1000 μm) of the sources contributing to far-IR counts and CIB. Combining the CIB distribution found above and the notion of Malmquist bias, it is no surprise that the resolved background is dominated by normal galaxies (L < 10¹¹   L_⊙) at low redshift and more luminous sources dominate the high redshift bins. Quantitatively, 95% of the resolved 160 μm CIB at z ≤ 0.5 is produced by normal galaxies,  > 90% of the contribution at 0.5 < z ≤ 1.0 arises from luminous infrared galaxies (LIRGs, 10¹¹ ≤ L < 10¹²   L_⊙) and ultra-luminous infrared galaxies (ULIRGs, 10¹² ≤ L < 10¹³   L_⊙) provide 50% of the background at 1.0 < z ≤ 2.0 and 88% above. Globally, roughly 50% of the CIB resolved in the three PACS bands in GOODS-S was produced by LIRGs. At 160 μm, the remainder is equally distributed between normal IR galaxies and ULIRGs, while at shorter wavelengths the former prevail.

COSMOS shallow observations probe the bright end of PACS number counts, and hence are limited to lower redshifts, but the large area grants a greater detail in the distribution of CIB, allowing a much finer binning. At the COSMOS 80% completeness limits (9.0 and 20.0 mJy at 100 and 160 μm, respectively), the resolved CIB peaks at z ≃ 0.3 and exhibits a weak secondary “bump” at z ≃ 0.7 in both bands. Recently, Jauzac et al. (2011) obtained similar results based on stacking of 24 μm sources on Spitzer MIPS 70 and 160 μm. These authors find a dip at z ≃ 0.5 in the differential CIB brightness in COSMOS, i.e. the same secondary peak a z ≃ 0.7 found in PACS data. This features are certainly due to the z = 0.73 structure known in the COSMOS field. Only a small fraction (≤10%) of the total CIB is locked in z < 0.5 bright objects, which are dominated by the non-evolutionary component of number counts (see for example Franceschini et al. 2010; Berta et al. 2010).

The bottom panels in Fig. 11 present the redshift derivative of the background surface brightness, and compare it to a set of model predictions. The derivative d(νI_ν)/dz peaks in the 0.5 < z ≤ 1.0 bin and then drops at higher redshifts. Le Floc’h et al. (2009) find a similar behavior at 24 μm, based on deep COSMOS data. The decrement is very steep at 70 μm, while it becomes milder at longer wavelengths. While over broad redshift binning the five models considered are overall consistent to the data, at the low redshift bright-end, where detailed information is available, they significantly differ, especially at longer wavelengths. Again, this might be a sign of possible differences in SED color- and temperature-luminosity assumptions and evolution.

The fine details of the comparison between models and data in the COSMOS field are affected by the known z ~ 0.7 overdensity. In order to study the impact of these features, the COSMOS and GOODS-S catalogs have been cut to a common (overlapping) flux range. In this case, the d(νI_ν)/dz derived in the two areas are consistent to each other within the errors, although the fine structure of GOODS-S cannot be studied. If then the same broad redshift binning is adopted, the features due to COSMOS large scale structure are washed out, and the trends seen in the two fields fully resemble each other.

Figure 13 summarizes our findings and compares the PEP CIB estimates to direct measurements and upper limits from γ-rays opacity. Our estimate of the background surface brightness emitted above the PEP detection and P(D) limits is fully enclosed in the COBE 1σ error bars. The uncertainties on the CIB values derived from PACS resolved counts, P(D), and power-law extrapolation are overall significantly smaller than what achieved previously by COBE/DIRBE. Far-IR number counts are currently in the position to provide a more accurate estimate of the extragalactic background light than direct photometric measurements.

7. Summary and conclusions

Using data belonging to the PACS Evolutionary Probe Herschel survey (Lutz et al. 2011), we derived far-IR number counts in the GOODS-S, GOODS-N, Lockman Hole and COSMOS fields. By employing stacking and P(D) analyses, it was possible to exploit the whole information in PACS maps and extend the study to very deep fluxes. The wealth of ancillary data in the GOODS and COSMOS fields allowed counts and CIB to be split into redshift bins. Through integration of number counts, we computed an unprecedented estimate of the contribution of IR galaxies to the extragalactic background light, across cosmic epochs. Among the many analyses carried out here, the main results that were presented are:

• resolved far-IR number counts at 70, 100, 160 μmextend from a few mJy to  ~200 mJy.In Euclidean-normalized units, our combination of four fieldstightly defines the bright side of counts, the peak at intermediatefluxes, the downturn and the slope at the faint-end. Stacking and P(D)analyses push the knowledge of counts by a decade fainter in flux,down to  ~0.2 mJy in GOODS-S;

• based on the observed counts, we derive a photometric confusion noise estimate of σ_c = 0.27 and 0.92 mJy (obtained with q = 5) at 100 and 160 μm respectively, and 16.7 beams/source density confusion limits of 0.4,  ~2.0 and 8.0 mJy in the three PACS bands;

• the fraction of CIB emitted at fluxes brighter that PEP 3σ thresholds is 58% (74%) of the Dole et al. (2006) reference value at 100 μm (160 μm). Using P(D) statistics, this fraction increases to  ~65% (~89%); this estimate lies well within the 1σ uncertainty of the direct CIB measurements, with a 3σ confidence interval smaller than COBE’s by a factor of 2 at 160 μm;

• exploiting power-law extrapolation of number counts at the faint side, and including past IRAS, ISO and Spitzer results and Euclidean extrapolations at the bright end, we derive new expectation values for the total CIB at 70, 100 and 160 μm: νI_ν ≥ 11.09, and [nW m^-2 sr^-1];

• at the current 3σ detection threshold, more than half of the resolved CIB was emitted at redshift z ≤ 1. The half-light redshift lies at z = 0.58, 0.67 and 0.73 at 70 (GOODS-S), 100 and 160 μm (GOODS-S and COSMOS combined), respectively. The balance moves towards higher redshift at longer wavelengths: while at 70 μm roughly 80% of the resolved CIB is emitted at z ≤ 1.0, at 160 μm the CIB budget is almost equi-partitioned below/above z = 1.0. Roughly 50% of the CIB observed at z = 0 and resolved in the three PACS bands in GOODS-S is emitted by LIRGs;

• most of the available evolutionary models fairly reproduce PACS total counts, despite the large variety of adopted assumptions. On the other hand, the detailed high-quality PACS data highlight dramatic differences between models and severe failures when compared to redshift distributions and CIB derivatives. Only models actually tuned taking into account redshift information and early Herschel data are successful in reproducing the new PEP dataset, pointing out the need to include these constraints in future modeling attempts.

Herschel and PACS set a new reference of the brightness of the cosmic far-IR background. Expectantly, next analyses and missions, rather than trying to recover the missing fraction of CIB, envisage the reconstruction of the bolometric background from the mid-IR to sub-mm frequencies. The contribution of galaxies at different redshift to the CIB is being investigated to an increasing details by current surveys. A further study of the redshift-dependent CIB, e.g., based on luminosity functions and properly including the effects of k-correction, goes beyond the purpose of this paper and is deferred to future works. This information is not only important in the fine tuning of galaxy evolution models, but is also of fundamental importance to constrain the intrinsic spectra of distant TeV sources, such as BLAZARs, whose γ-ray photons interact with the CIB. So far, such constraints have been based solely on model predictions, to be confirmed with direct estimates of the CIB as observed from different epochs. Herschel source statistics have demonstrated to provide a precise estimate of the total CIB. Future studies will necessitate to improve direct measurements and thus the study of foreground contaminants, in order to complete the understanding of the total background and pinpoint its measurement.

---

Acknowledgments

We wish to thank M. Bethermin, A. Franceschini, C. Gruppioni, G. Marsden, S. Niemi, M. Rowan-Robinson, and E. Valiante for providing new, extended model outputs for us, and the anonymous referee for useful comments. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain).

References

Online material

Appendix A: A cheap P(D) estimator

In addition to directly detected sources and to stacking of 24 μm dense catalogs, further information about the shape of PACS number counts comes from the statistical properties of observations, probing the counts at even fainter fluxes. The “probability of deflection”, or P(D) distribution, is basically the distribution of pixel values in a map (see Sect. 4).

For a large density of sources, the contribution to the P(D) is a Poisson distribution with a large mean value, convolved with the instrumental noise (typically assumed to be Gaussian). Generally speaking, for sufficiently steep number counts at faint flux densities, the depth reached with the P(D) analysis is significantly higher than what is provided by individually-detected sources, and is a powerful tool to probe counts slopes in cases dominated by confusion.

Given differential number counts dN/dS, and a beam function (e.g., PSF) f(θ,φ), describing the point-source spatial-response at position x, the probability of deflection P(D) can be written as (e.g., Patanchon et al. 2009): (A.1)or, alternatively, isolating the real part only (e.g., Franceschini et al. 2010): (A.2)where the mean number of source responses of intensity x in the beam is given by: (A.3)We defer to Patanchon et al. (2009), Franceschini et al. (2010), Glenn et al. (2010) for a full treatment and description of the P(D) formalism. Only in the case of very simple dN/dS functional forms and trivial beams (e.g., Gaussian), the above equations can be solved analytically. For an effective beam that is not strictly positive, or is not azimuthally symmetric, it is necessary to use the full 2D beam map.

For example, the Herschel/PACS beam is characterized by 3 prominent lobes, which cannot be described analytically. In such a case, if one likes to account for the real beam function, a numerical treatment of R(x), and possibly P(D) is necessary.

Finally, in case of real observations, the instrumental noise contribution to pixel flux densities must be included.

A.1. The numerical P(D) approach

Given an instrumental PSF f(θ,φ) and a number counts model dN/dS, it is possible to predict the pixel flux density distribution of a map with M pixels, simply rolling random numbers. We describe here the principles of the method we developed in the frame of the PEP survey. This approach turns out to be of simple implementation, avoiding integrations of oscillating functions or Fourier Transforms, and it is accessible to any user with a basic knowledge of random number generators. Furthermore, it avoids any analytic simplification of the beam function, thus allowing to employ the real instrumental PSF, regardless of how complex its bi-dimensional profile is.

For a given dN/dS, we describe how a synthetic P(D) is computed. The adopted form of number counts and fitting technique are described in the next Section. Parameters are varied by means of a MCMC engine and for each realization a new P(D) is computed and compared to the observed data.

The input information needed by the P(D) numerical algorithm are:

• description of number counts, dN/dS;

• beam function, e.g., in form of an observed PSF 2D map;

• value of instrumental noise.

Further secondary inputs, derived from the three primaries above, are:

• the maximum radius r_max in pixels from a given source within which apixel is affected by the source itself, e.g., the radius at which 99%of the observed PSF flux is enclosed; if the pixel scale of the map is αarcsec/pix, the corresponding angular radius is αr_max;

• the total expected number, N₀, of sources within the solid angle Ω defined by r_max, simply given by the integral of number counts:

• the cumulative probability that a random source has a flux lower than S, given by: .

As mentioned above, the procedure relies on several dice rolls, using random numbers to determine the actual number of sources expected to contribute to each pixel of the map, the distance of the given source to the pixel, the fluxes of sources.

• Step 1.

  We assume that the distribution of sources in the skyis randomly uniform, i.e. we ignore the clustering prop-erties of sources in the sky. Clustering would mainly af-fect the P(D) over scales larger than the size of the GOODS-S field under exam here (e.g., Glennet al. 2010; Vieroet al. 2009; Lagacheet al. 2007). Given N₀, the number of sources affecting the given pixel in the current realization is a random Poisson deviate for expected N₀. This number is recomputed for each pixel and each dN/dS realization.

• Step 2.

  For each source affecting the given pixel, it is necessary to determine the fraction f_i of flux affecting the pixel. Roll N random, uniformly distributed positions, inside r_max (defined either by r,φ or by x,y coordinates). Based on the position of the given source, and the beam function, derive the fraction of flux affecting the pixel for each source.

• Step 3.

  Finally assign a random flux S_i to each of the N sources affecting the given pixel, according to the input dN/dS. To this aim, the probability P(S) (defined above) is used to transform a flat random number distribution in the needed flux values.

Assuming that any further background has been subtracted from the map, the flux of the given pixel is finally given by , where D₀ is the term describing pure noise (i.e. random flux, in case of no sources at all) and is given by a Gaussian deviate as wide as the noise measured on actual PEP rms maps (see Table 1).

A.2. Number counts model and minimization

Number counts are modeled similarly to Patanchon et al. (2009) and Glenn et al. (2010). We adopt a simple, parametric model, defined by the differential counts dN/dS at a set of flux density knots, or nodes. The model is defined within a S_min and S_max lower and upper flux boundaries; outside of this flux range, we simply set dN/dS to zero. Between nodes, number counts are described as power-laws, connecting at knots.

We used a Markov-Chain Monte Carlo (MCMC) sampling of the parameter space, to explore the likelihood function of the model to reproduce the observed P(D). We use a Metropolis-Hastings algorithm as our MCMC sampler. In order for this sampler to converge efficiently, we draw new steps from univariate Gaussian distributions, whose widths were given by dummy MCMC runs.

The position of knots and the corresponding amplitude of counts are free parameters, and are varied at each MCMC map realization. The upper boundary of the covered flux range is fixed to 100 mJy, a value driven by the limits of observed, resolved counts in GOODS-S. In order to minimize the number of free parameters, and further degeneracies, we limit the description of number counts to a broken power-law with three sections, for a total of seven free parameters (3 positions and 4 amplitudes).

In the past, a standard approach was to subtract bright individually detected sources from the P(D) analysis, and work only at the faint side, but it has been shown that a more robust estimate of the source counts is obtained if the P(D) is fitted across the whole available flux density range (e.g., Patanchon et al. 2009). Therefore, to derive the observed P(D) to be reproduced, we use the observed science maps without subtracting any object.

Section 4, Figs. 8, 7 and Table 10 describe the results of the whole P(D) analysis.

A.3. Effects of map-making high-pass filtering

It is necessary to recall that PEP maps have been obtained by applying a high-pass filtering process along the data time-line, in order to get rid of 1/f noise. This procedure consist in a running-box median filter, and induces two different side effects onto PACS maps: 1) the background is naturally removed from the maps; 2) the PSF profile is eroded and wings are suppressed.

At very faint fluxes, and high source densities, close to and beyond the confusion limit, the contribution of sources to the model P(D) basically generates a “background” sheet of sources, which shifts the P(D) histogram to a non-zero peak and median. In principle, the source component of the P(D) is always positive, and the addition of noise generates the negative nearly-Gaussian behaviors at the faint side of P(D). The adopted high-pass filtering shifts the whole flux pixel distribution function to a zero peak/median. In order to take this into account, with the current form of PEP maps, we shift the model P(D) to the same zero peak, prior of fitting.

The second consequence of applying an high-pass filter is a modification of object profiles. Detected objects are masked and excluded from the median filter derivation. Tests were done adding simulated sources to the time-lines before masking and before high-pass filtering (see Popesso et al., in prep.; Lutz et al. 2011). The result is that the filtering modifies the fluxes of masked sources below the detection limit – by a factor  ≤ 16%. In order to understand the impact of this on the P(D) analysis, we built simulated maps 20 times deeper than the GOODS-S 3σ threshold, once using an un-filtered (i.e. masked) PSF at all fluxes, and once using an un-filtered PSF above 3σ, but a filtered (unmasked) PSF below. The P(D) obtained in the two cases are well consistent to each other, with unit median ratio and a  < 5% scatter across all fluxes; no systematic shift, nor skewness are observed.

All Tables

All Figures