© ESO, 2013

1. Introduction

The detection of a cosmic infrared background, as energetic as the optical/near-infrared background (Puget et al. 1996; Hauser et al. 1998), has revealed the importance of the energy absorbed by the dust in galaxies and re-emitted at mid- to far-infrared wavelengths. From this finding it became clear that a complete census on the formation and evolution of galaxies could not be obtained without accounting for this dust emission. Since then, many studies have confirmed the importance of dust emission using individual detections of galaxies from infrared facilities such as the Infrared Space Observatory (ISO), or the Spitzer Space Telescope. However, because of their relatively small onboard optics, the far-infrared capabilities of these observatories were strongly limited by source confusion. At the time, far-infrared studies were restricted to the analysis of local galaxies or to the analysis of rare high-redshift very luminous galaxies. As a consequence, only a small fraction of the cosmic far-infrared (i.e., λ > 40  μm) background was resolved into individual objects, so our knowledge of the high-redshift Universe at far-infrared wavelengths was very incomplete.

Thanks to the advent of the PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) instruments onboard the Herschel Space Observatory (Pilbratt et al. 2010), this limitation has been largely overcome. Indeed, using the relatively high spatial resolution (provided by a 3.5  m mirror) and sensitivity of Herschel, deep extragalactic surveys can be pursued, thereby resolving a large fraction of the cosmic far-infrared (i.e., ~58% and ~74% with PACS at 100 and 160 μm, respectively; Berta et al. 2010, 2011) and submillimetre (i.e., ~15% with SPIRE at 250 μm; Oliver et al. 2010) background into individually detected galaxies. From these observations, one can study the origin and nature of the cosmic infrared background through, e.g., the determination of the infrared luminosity functions of galaxies (e.g., Gruppioni et al. 2010, 2013; Casey et al. 2012) and the examination of their spectral energy distributions (e.g., Hwang et al. 2010; Elbaz et al. 2010, 2011; Nordon et al. 2010, 2012; Magnelli et al. 2010, 2012; Symeonidis et al. 2013; Berta et al. 2013). All these studies point towards the diversity and redshift evolution, both in term of numbers and properties, of the infrared luminous galaxy population. Relatively rare in the local Universe, the infrared luminous galaxies dominate the cosmic star-formation history at z > 1 (e.g., Gruppioni et al. 2010, 2013), and their physical properties differ significantly from those of their local counterparts (e.g., Elbaz et al. 2011; Wuyts et al. 2011; Nordon et al. 2012).

In this context, we study here the infrared luminous galaxy population further by deriving the PACS numbers counts and infrared luminosity functions down to unprecedented depths using the combination of the two main extragalactic surveys designed to take advantage of the full PACS capabilities: the PACS Evolutionary Probe (PEP¹; Lutz et al. 2011) guaranteed time key programme; and the GOODS-Herschel (GOODS-H²; Elbaz et al. 2011) open time key programme. The PEP survey is structured as a “wedding cake” (i.e., with large area shallow images and smaller deep images) and includes many widely studied blank and lensed extragalactic fields, such as the Great Observatories Origins Deep Surveys North (GOODS-N) and South (GOODS-S) fields and the cosmological evolution survey (COSMOS). The GOODS-H survey only focuses on the GOODS fields, but using very deep observations, close to the Herschel confusion limit. The extensive observations of the GOODS fields made by the PEP and GOODS-H surveys is explained by the availability of a deep multi-wavelength database including X-ray (Alexander et al. 2003; Xue et al. 2011), optical (Giavalisco et al. 2004), near-infrared (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011), mid-infrared (GOODS-Spitzer; PI: M. Dickinson; see Magnelli et al. 2011), (sub)mm (e.g., Oliver et al. 2012; Borys et al. 2003; Weiß et al. 2009) and radio (Morrison et al. 2010; Miller et al. 2008) observations.

We present here the combination of the PEP and GOODS-H observations of the GOODS fields. This combination provides the deepest observations of the GOODS fields obtained by PACS. In particular, in the GOODS-South field the combination of these observations is not limited by the exposure time but by confusion. In this field, we thus obtain the deepest blank field observations achievable with PACS onboard the Herschel Space Observatory. In this paper we present in detail the data analysis method used to produce the publicly available PEP/GOODS-H maps and catalogues³. Then we use these deep observations to constrain the PACS-100 μm and –160 μm confusion noises and number counts, and to study the evolution of the infrared luminosity function and of the star-formation rate history of the Universe up to z ~ 2.

The paper is structured as follows. Observations are described in Sect. 2. Section 3 presents the method used to produce the PACS maps. In Sect. 4, we present our source extraction methods and contents of the released package. In Sect. 5, we estimate the PACS-100 and –160 μm confusion noise. Number counts are presented in Sect. 6 while in Sect. 7 we derive the infrared luminosity functions of galaxies as well as the star-formation rate history of the Universe up to z ~ 2. Finally, we summarise our results in Sect. 8. Throughout the paper we use a cosmology with H₀ = 71  km  s^-1  Mpc^-1,Ω_Λ = 0.73 and Ω_M = 0.27.

2. Observations

The PACS maps and catalogues used and released in this paper are obtained from the combination of the PEP and GOODS-H observations of the GOODS-N and GOODS-S fields (see Table 1). Both programmes have observed the GOODS fields using the scan mode of the PACS photometer on board Herschel. This mode consists of slewing the spacecraft back and forth along parallel lines at a constant speed of 20′′ s^-1. Using this scan mode, astronomical observing requests (AORs) were designed to observe the GOODS fields in both nominal and orthogonal directions. To reach the desired depth many AORs per field were required. The central position of each AOR was dithered by ~8′′ in order to improve the spatial redundancy of the data.

Using this strategy, the PEP and GOODS-H surveys covered the entire 11′ × 17′ GOODS-N field with PACS at 100 and 160 μm⁴. The total observing time (i.e., including overheads) in GOODS-N was 25.8 and 124 h for PEP and GOODS-H, respectively (Table 1).

The entire 11′ × 17′ GOODS-S field was observed by PEP with PACS at 70, 100 and 160 μm⁵. The total observing times were 113, 113 and 226 h at 70, 100 and 160 μm, respectively, including 10 h guaranteed time contributed by Herschel mission scientist Martin Harwit. The Extended Chandra Deep Field South (ECDFS), containing in its centre the GOODS-S field, has also been observed by PEP with a total observing time of 32.8 h over a 30′ × 30′ region. Observations of ECDFS are included in our combined maps, trimmed to match the GOODS-S layout of the PEP observations. GOODS-H has observed a 10′ × 10′ region centred on the GOODS-S field with PACS at 100 and 160 μm. The choice of covering only a fraction of the GOODS-S field was made in order to obtain, in a reasonable amount of time, a PACS-100 μm map close to the confusion limit of Herschel. The GOODS-H observations of the GOODS-S field lasted for a total of 206.3 h. Finally in our combined maps we also include 7.9 h of observations of the GOODS-S field taken with preliminary instrument settings during the Herschel performance verification phase (Table 1).

Fig. 1 PACS-100  μm coverage map of the GOODS-S field. Due to the different layout of the PEP and GOODS-H observations, the 10′ × 10′ region centred on the GOODS-S field has much deeper PACS-100 μm observations than the outskirts. We observe the same pattern in the PACS-160  μm coverage map of the GOODS-S field.

Because PEP and GOODS-H observations were executed using the same observing mode, a combination of these data sets can easily be performed. Of course, due to the different layout of PEP and GOODS-H observations in GOODS-S, this field is not homogeneously covered at 100 and 160 μm. As illustrated by Fig. 1, a 10′ × 10′ region centred on the GOODS-S field has higher PACS-100 and 160 μm coverage than the outskirts. In contrast, the GOODS-S 70 μm coverage is uniform across the field, with the exception of fall-off at the edges. In the rest of the paper, the centred region of the GOODS-S field with ultradeep observations is referred as “GOODS-S-ultradeep”, while the outskirts are referred to as “GOODS-S-deep”. This combined data set provides the deepest Herschel far-infrared observations of the GOODS-N and -S fields. In the rest of the paper this combined data set is referred to as the “PEP/GOODS-H” observations.

Our combined PACS-100 μm maps reach a total observing time per sky position of ~2.6 h, ~2.6 h and ~10.0 h in GOODS-N, GOODS-S-deep and GOODS-S-ultradeep, respectively. The PACS-160 μm maps reach a total observing time per sky position of ~2.6 h, ~4.7 h and ~12.1 h in GOODS-N, GOODS-S-deep and GOODS-S-ultradeep, respectively. In the GOODS-S field, the PACS-70 μm map reaches a total observing time per sky position of ~2.1 h.

We note that a detailed description of the observational strategies adopted by PEP and GOODS-H is provided in Lutz et al. (2011) and Elbaz et al. (2011), respectively.

3. Map creation

Observations were reduced using the standard PACS photometer pipeline (Wieprecht et al. 2009) and some custom procedures, all implemented within the HIPE environment in the Herschel common science system (HCSS). This data reduction process is described in more detail in Lutz et al. (2011) and Popesso et al. (2012). Here we summarise the main steps.

The data reduction process starts at the AOR level and is based on the scanmap script of the PACS photometer pipeline. First, the pipeline flags bad or saturated pixels, converts detector signals from digital units to volts, finds pixels affected by short glitches and replaces their values using a standard interpolation method, and finally it applies a recentring correction to the pointing product of Herschel using reference positions of 24 μm sources with accurate astrometry (see Lutz et al. 2011) from deep observations with the Multiband Imaging Photometer (MIPS; Rieke et al. 2004) onboard the Spitzer Space Telescope. Thanks to these recentring corrections, positions of PACS sources measured in a “blind” catalogue (see Sect. 4.2) are offset by less than 0.2″ and have a rms difference of ~1″ with respect to the position of their MIPS-24 μm counterparts (see also Lutz et al. 2011). We note that the MIPS-24 μm astrometry, and therefore our PACS maps, match the GOODS ACS version 2 coordinate system⁶. Second, the pipeline removes from the timeline of each AOR the “1/f” noise which is the main source of instrumental noise in PACS data. For deep cosmological surveys such as ours, the PACS pipeline does so by using a high-pass filtering method which subtracts, from each timeline, a version of the timeline filtered by a running box median of a given radius (expressed in readouts, i.e., in numbers of points of the timeline). The presence of sources in the timeline affects the high-pass filtering method by artificially boosting this running box median, i.e., leading to the subtraction of part of the source flux (Popesso et al. 2012). Consequently, sources have to be “masked” from the timelines. Using results from Popesso et al. (2012), we choose a masking strategy based on circular patches at prior positions. This method reduces the amount of flux loss due to the high-pass filter and more importantly leads to flux losses which are independent of the PACS flux densities (Popesso et al. 2012). We note that this improved masking strategy differs from that adopted in Lutz et al. (2011). Timelines are masked at the position of 24 μm sources with S₂₄ > 60  μJy using circular patches with radius of 4″, 4″ and 6″ at 70, 100 and 160 μm, respectively. This strategy allows us to mask almost all PACS detections and corresponds to the masking of about ~12% of the timeline. A small number of resolved sources (mainly in the GOODS-S field) are masked using larger patches which are adjusted visually. After being masked, timelines are high-pass filtered using a running box median with radius of 12 readouts (24′′), 12 readouts (24′′) and 20 readouts (40′′) at 70, 100 and 160 μm, respectively. These radii are chosen based on results from Popesso et al. (2012) . Although we use an optimized high-pass filtering strategy, PACS flux densities still have to be corrected for flux losses. The corrections are provided in Popesso et al. (2012) for our high-pass filter setting, i.e., 13%, 12% and 11% at 70, 100 and 160 μm, respectively (see Sect. 4).

After the removal of the “1/f” noise, AORs are flat-fielded and flux calibrated. Then, the map of each AOR is created using a “drizzle” method (Fruchter & Hook 2002) implemented in the HCSS. Because all AOR maps are projected on the same world coordinate system, they are coadded into a final map using, as weight, the effective exposure time of each pixel, appropriately underweighting the performance verification phase data which were taken with preliminary instrument settings. Taking advantage of the large number of AORs, we create uncertainty maps from the standard deviation of this weighted mean in each pixel. Although, a “drizzle” method is applied for the creation of each AOR map, the final map still contains some correlated noise. Thanks to the high redundancy of our data at each sky position, we are able to calculate the mean correlated noise across the map and within our point spread functions (PSFs; see Lutz et al. 2011). These mean correlation corrections were taken into account while measuring uncertainties from the uncertainty maps (i.e., uncertainties have to be corrected upward by a factor ~1.5; Sect. 4.2).

Fig. 2 (Top left) PEP/GOODS-H map of the GOODS-north field (11′ × 17′) at 100 μm. (Top right) PEP/GOODS-H map of the GOODS-south field (11′ × 17′) at 100 μm. Contours correspond to exposure times greater than 0.5 (at the edge), 3, 6 and 9 (centre) hrs/pix. (Bottom left) Colour composite image of the GOODS-north field at 24 μm (blue), 100 μm (green) and 160 μm (red). (Bottom right) Colour composite image of the GOODS-south field at 24 μm (blue), 100 μm (green) and 160 μm (red). The 24 μm images (PI: M. Dickinson) were obtained by the Spitzer Space Telescope, while the 100 and 160 μm images were obtained by the Herschel Space Observatory. In the colour composite images, sources with contribution from the 24 + 100  μm, 24 + 160  μm, 100 + 160  μm and 24 + 100 + 160  μm bands would correspond to a cyan, magenta, yellow and white colours, respectively. The relatively noisy edges of the colour composite images have been trimmed.

Figure 2 shows the combined PEP/GOODS-H PACS-100  μm maps of the GOODS-N (top left panel) and GOODS-S (top right panel) fields. Three-colour composite images of the GOODS-N (bottom left panel) and GOODS-S (bottom right panel) fields at 24, 100, 160 μm are also shown in this figure.

Fig. 3 PACS-to-24 μm flux density ratio as a function of the PACS flux density inferred from the “blind” source catalogue cross-matched with the MIPS-24 μm catalogue. Shaded regions show the space density distribution of galaxies at 70 μm (top left), 100 μm (top right) and 160 μm (bottom) in the GOODS-S field. In this field, the MIPS 24 μm catalogue reaches a 3σ limit of 20 μJy, while the PACS-70, 100 and 160 μm data reach 3σ limits of 0.9, 0.6 and 1.3 mJy (vertical dotted black lines). The parameter spaces reachable by our prior source extraction method using the MIPS catalogue are located below the dashed black lines: PACS sources located above these lines will not be detected at 24 μm. Orange triangles, red circles, light blue squares and dark blue stars show the PACS-to-24 μm flux density ratios of galaxies with LIR = 1011.5 L⊙ as predicted at different redshifts using the spectral energy distributions (SEDs) of normal (MS), starburst (SB) galaxies (Elbaz et al. 2011), IRAS 22491 (Berta 2005) and Arp220 (Silva et al. 1998), respectively. The right point of each track corresponds to z = 0.5, while other points correspond to increasing redshifts, with intervals of Δz = 0.5 (arrows indicate the path followed for increasing redshift). Predictions for galaxies with lower or higher infrared luminosities will correspond to a global vertical shift of these tracks towards lower or higher PACS flux densities, respectively. Most of the PACS population lies well within the limits of the parameter space reachable by our MIPS 24 μm catalogues, i.e., below the dashed lines.

4. Source extraction

At the resolution of PACS most of the sources in our fields are point sources (i.e. FWHM ~ 4.7″, 6.7″ and 11″ at 70, 100 and 160 μm, respectively). Therefore, we use PSF fitting to derive their flux densities. Two catalogues are derived using complementary approaches. First, we construct a catalogue using as priors the source positions expected on the basis of a deep 24 μm catalogue. This provides good deblending of neighbouring sources, but will miss a few Herschel sources that are not detected at 24 μm (i.e., <1% and <4% of the PACS sources in the GOODS-N and GOODS-S fields, respectively). Second, we provide a “blind” catalogue using PSF fitting without positional priors.

4.1. Prior source extraction

Starting from the positions of IRAC-3.6 μm sources (Infrared Array Camera; FWHM ~ 1.6′′) from the GOODS Spitzer Legacy Program (PI: M. Dickinson), we extract sources in the MIPS-24 μm maps (Magnelli et al. 2011)⁷ and then use the 24 μm-detected sources (i.e., with S₂₄ > 3σ ~ 20  μJy) as priors for the source extraction in the PACS maps. The main advantage of this approach is that it deals with a large part of the blending issues encountered in dense fields and provides a straightforward association between IRAC, MIPS and PACS sources (Magnelli et al. 2011). The disadvantage of this method is that we must assume that all sources present in our PACS images have already been detected at 24 μm. Based on the relative depth of the mid- and far-infrared images of the GOODS fields, Magdis et al. (2011) have investigated this assumption. They found that in the GOODS-H observations of the GOODS-S field, less than 2% of the PACS sources are missed in the MIPS-24 μm catalogue. In our combined dataset, PACS observations of the GOODS-N field are shallower than those analysed in Magdis et al. (2011), while MIPS-24 μm observations are equivalently deep (i.e., 3σ ~ 20  μJy in both GOODS-N and GOODS-S). Consequently, in the GOODS-N field, the fraction of PACS sources with no MIPS-24 μm counterparts should be significantly lower than 2%. This result is confirmed by cross-matching our GOODS-N “blind” PACS and MIPS-24 μm catalogues using a matching radius of 4′′ and 6′′ at 100 and 160 μm, respectively: less than 1% of the PACS “blind” sources are missed in the MIPS-24 μm catalogue. In GOODS-S, the situation is somehow different as in this field our PACS observations are deeper than those analysed in Magdis et al. (2011). By cross-matching our GOODS-S “blind” PACS and MIPS-24 μm catalogues, we find that 4.1%, 4.2% and 3.5% of the PACS sources have no MIPS-24 μm counterparts within 3′′, 4′′ and 6′′ at 70, 100 and 160 μm, respectively. These fractions should be treated as upper limits, given that most of these sources have faint PACS flux densities (i.e., 80% with S/N < 5) and thus could be spurious (see Sect. 4.4). To understand the spectral properties of PACS sources with no MIPS-24 μm counterparts, we study the typical PACS-to-MIPS flux ratio observed in the GOODS-S field (Fig. 3). While the bulk of the PACS-100 (–160) μm population lies well within the parameter space reachable by our MIPS-24 μm catalogue (i.e., below the dashed lines of Fig. 3), there is, at faint flux densities, a slight truncation of the high-end of the dispersion of the PACS-to-MIPS flux ratio. This truncation will translate into faint PACS-100 and –160 μm sources with no MIPS-24 μm counterparts. In addition, galaxies with similar spectral properties to Arp220 (i.e., with high PACS-to-MIPS flux ratio and strong silicate absorptions) will also be missed in MIPS-24 μm catalogues, especially at z ~ 0.4 and z ~ 1.3 where the 18 μm and 9.4  μm silicate absorption features are shifted into the MIPS-24 μm passband (see also Magdis et al. 2011). The existence of PACS sources with no MIPS-24 μm counterparts will naturally introduce slight incompleteness (i.e., <1% and <4% in the GOODS-N and GOODS-S fields, respectively) in our prior catalogues. However, because this incompleteness mostly appears at faint flux densities (i.e., S/N < 5), it should not be much larger than that introduced by prior-free source extraction methods at such low signal-to-noise ratio (S/N). Indeed, at faint flux densities, there are only small discrepancies, in terms of number of sources, between our blind and prior catalogues (see Appendix A). We note that we cannot use IRAC priors to reduce this incompleteness because the high IRAC source density (i.e., 2–4 priors per PACS beam) would force each far-infrared source to be deblended into several unrealistic counterparts.

From the expected positions of sources, we fit empirical PSFs (see Sect. 4.3) to our PACS maps. PACS flux densities are then defined as the intensity of the scaled PSFs. The simultaneous fit of nearby sources optimizes the deblending of their flux densities. As for a standard aperture measurement, the photometry of each source is corrected to account for the finite size of our empirical PSFs, i.e., they do not include the power contained in the wings of the real PSFs (see Sect. 4.3). Finally, because of flux losses from the high-pass filtering, additional corrections are applied to our flux density measurements (see Sect. 3).

Flux uncertainties are estimated on residual maps. These uncertainties are defined as the pixel dispersion, around a given source, of the residual map convolved with the appropriate PSF. This method has the advantage of taking into account, at the same time, the rms of the map (including correlated noise and inhomogeneities in the exposure time) and the quality of our fitting procedure. Because our far-infrared observations have been designed to reach the confusion limit of Herschel, our flux uncertainties also contain a part of the confusion noise left in the residual maps (i.e., confusion noise, due to sources fainter than PACS detection limits). In this context of complex combination of instrumental and confusion noise, we test the accuracy of our flux uncertainties using extensive Monte Carlo (MC) simulations (see Sect. 4.4). We find that the flux uncertainties inferred from the residual maps and the photometric accuracies inferred from the MC simulations are in good agreement.

4.2. Blind source extraction

PACS flux densities are also extracted using a “blind” PSF-fitting analysis, i.e., without taking, as prior information, the expected positions of sources. This blind source extraction is performed with Starfinder (Diolaiti et al. 2000a,b). This code is appropriate for the extraction of the unrevolved PACS sources, because it was especially designed to obtain high precision astrometry and photometry of point sources in crowded fields. Basically, Starfinder carries out source extraction using three steps: (i) detection of candidate point sources; (ii) verification of the likelihood of the candidate point sources; (iii) fits of empirical PSFs to the map, at positions of the candidate point sources.

Input parameters of Starfinder are fine-tuned using MC simulations (see Sect. 4.4). Starfinder performed its PSF-fitting analysis with the same empirical PSFs as those used for the prior source extraction (see Sect. 4.3). Flux densities inferred by Starfinder are corrected for the finite size of the empirical PSFs and for the effect of high-pass filtering (see Sects. 3 and 4.3). From the MC simulations we observe that, despite the fine-tuning, Starfinder still tends to overestimate the flux density of sources at low S/N levels (about 10% at S/N ~ 3). This behaviour is corrected using factors derived from the MC simulations.

Starfinder being exclusively designed for point source extraction, the presence of a few extended sources (~8 sources per field) in the PACS maps affects its detection process: extended sources are split into sub-components. To fix this problem, we run Sextractor (Bertin & Arnouts 1996) on these extended sources using Kron elliptical apertures. Extended sources are identified using an empirical method exploiting the isophotal area vs flux parameter space. In this parameter space, extended sources have large isophotal areas compared to their flux densities. PACS flux densities of the sub-components of the extended sources are erased from our blind source catalogue and replaced by flux densities inferred by Sextractor. We note that the accuracy of the Sextractor flux densities was verified using the Sextractor residual maps.

Flux uncertainties are estimated with Starfinder using our PACS uncertainty maps (see Sect. 3) and empirical PSFs. These flux uncertainties are corrected for the finite size of the PSFs and for the effect of high-pass filtering. In addition, because our PACS uncertainty maps do not account for correlated noise, these flux uncertainties have to be corrected upwards using mean correlation corrections (Sect. 3). We note that after correction, the flux uncertainties inferred by Starfinder are in good agreement with the photometric accuracies inferred from the MC simulations.

4.3. PSF reconstruction

Empirical PSFs used for the prior and blind source extraction are derived directly from the PACS maps using Starfinder. A number of isolated point-like sources present in the maps are stacked and then normalized to unit total flux. Because of the limited number of isolated point-like sources, the wings of these empirical PSFs have limited S/N. Therefore, they have to be truncated to smaller radii (9.6″, 7.2″ and 12″ at 70, 100 and 160 μm, respectively) and then renormalized to a unit total flux.

Because of their finite extent, the empirical PSFs do not include all the power contained in the wings of the real PSFs. Consequently, flux densities inferred from these PSFs have to be corrected, as would be done for any standard aperture flux measurement. Aperture corrections are derived by comparing the empirical PSFs with the reference in-flight PSF for our observing mode, obtained from observations of the asteroid Vesta. First, the Vesta observations are manipulated to account for variations of the position angle (PA) between our AORs, i.e., the Vesta observations are rotated to the PA of each AOR and then stacked. Second, Vesta observations are smoothed to match the slightly broader FWHM of the empirical PSFs (broader by a factor 1.04, 1.02 and 1.01 at 70, 100 and 160 μm, respectively). These broader FWHM can be explained by residual pointing uncertainties left by our recentring procedure (see Sect. 3). Finally, aperture corrections are defined as the fraction of the total energy contained in these manipulated Vesta PSFs within the radii of the empirical PSFs. At 70, 100 and 160 μm, the empirical PSFs contain 77, 68 and 69% of the total energy of the manipulated Vesta PSFs, respectively.

4.4. Monte Carlo simulations

Fig. 4 Results of MC simulations in the GOODS-North and -South fields using the prior source extraction method. These MC simulations correspond to the central deepest regions of each field, i.e., avoiding the rather noisy edges of the GOODS-N field and concentrating on the GOODS-S-ultradeep part of the GOODS-S field. Blue lines represent the average photometric accuracy defined as the standard deviation of the (Sout − Sin)/Sout distribution in each flux bin (after 3σ clipping). Red lines show the mean value of the (Sout − Sin)/Sout distribution in each flux bin. Inset plots show the fraction of artificial sources detected in the image (i.e., completeness) as a function of input flux (orange plain histogram) and the fraction of spurious sources (i.e., contamination) as a function of flux density (striped black histogram).

The PEP/GOODS-H observations have been designed to reach the confusion limit of the Herschel Space Observatory at 100 and 160 μm. Flux uncertainties are therefore a complex combination of instrumental and confusion noise. In order to estimate these complex flux uncertainties and to characterize the quality of our catalogues we perform extensive MC simulations.

We add simulated sources to our PACS 70, 100 and 160 μm images with a flux distribution, approximately matching the measured number counts (Berta et al. 2010, 2011). The flux densities of the faintest simulated sources are defined as the PACS flux densities (i.e., S_λ) expected for the faintest MIPS-24 μm (i.e., S₂₄ = 20  μJy) sources, i.e., ${\mathit{S}}_{\mathit{\lambda }}^{\min }$ = min(S₂₄) × mean(S_λ/S₂₄), where mean(S_λ/S₂₄) is 10, 20 and 30 at λ = 70, 100 and 160 μm, respectively (see Fig. 3). To preserve the original statistics of the image (especially its source number density), the number of simulated objects added at one time is kept small (i.e., 20 sources). To recreate the clustering properties of PACS sources, simulated objects are positioned, with respect to the MIPS-24 μm sources, to reproduce the distance to the closest neighbour distribution observed in the MIPS-24  μm prior catalogues. Simulated sources are created using the manipulated Vesta PSFs which contain the wing of the real PSFs (see Sect. 4.3). On these simulated images, we then perform both our blind and prior source extraction using the empirical PSFs and compare the resulting flux densities to the input values. To improve the statistics, this process is repeated a large number of times using each time different positions and fluxes for the simulated sources. For each field and at each wavelength, a total of 20 000 artificial sources are used. Figure 4 shows, as an example, results from the MC simulations performed in the GOODS fields using the prior source extraction method. In both GOODS fields, these MC simulations correspond to the central deepest regions (i.e., avoiding the rather noisy edges of the GOODS-N field and concentrating on the GOODS-S-ultradeep part of the GOODS-S field). MC simulations on the outskirts of the GOODS-S field (i.e., GOODS-S-deep) have been run, but are not shown here. In any case, results from these MC simulations are very similar to those for the GOODS-N field, as expected from the similarity between the GOODS-S-deep and GOODS-N average exposure times.

From these MC simulations we derive three important quantities: the completeness; the contamination; and the photometric accuracy of our catalogues as a function of flux density. Completeness is defined as the fraction of simulated sources extracted with a flux accuracy better than 50%. The contamination is defined as the fraction of simulated sources introduced with S < 2σ but extracted with S > 3σ. The photometric accuracy is defined as the standard deviation of the (S_out − S_in)/S_out distribution as a function of S_out (blue lines in Fig. 4). These photometric accuracies have the advantage of taking into account simultaneously nearly all sources of noise, i.e., including confusion. Table 2 summarises these quantities for both the blind and prior source extraction approaches.

From these MC simulations we conclude that the prior and blind PSF-fitting methods perform accurate extraction of PACS sources. The prior and blind source catalogues are characterized by high completeness and low contamination levels, as well as by good photometric accuracy (i.e., <33%). We observe in Fig. 4 that, using our flux uncertainties, a S/N > 3 cut almost always translates into a photometric accuracy better than 33% (i.e., as expected from sources with S/N > 3). We can thus conclude that our flux uncertainties are fairly accurate. However, we also note that in the GOODS-S field the photometric accuracy is worse than 33% for flux densities below 0.6  mJy and 2.0  mJy at 100 and 160 μm, respectively. At such faint flux densities our PACS 100 and 160 μm maps are likely to be affected by confusion noise (i.e., σ_c(100  μm) ~ 0.15  mJy and σ_c(160  μm) ~ 0.68  mJy; see Sect. 5) that is not fully accounted in flux uncertainty estimated by our source extraction methods. This conclusion is further confirmed by the analysis of the (S_in − S_out)/σ_s distribution (σ_s being the flux uncertainty inferred by our source extraction methods). Indeed, in all but the GOODS-S 100 and 160 μm fields, the (S_in − S_out)/σ_s distribution follows, as expected, a Gaussian distribution with a dispersion almost equal to 1. Instead, in the GOODS-S 100 and 160 μm fields, we find that the (S_in − S_out)/σ_s distribution follow a Gaussian distribution with a dispersion of ~ 1.2, indicating that our flux uncertainties are underestimated and do not fully account for confusion noise. Therefore, we recommend caution in using flux densities lower than 0.6  mJy and 2.0  mJy at 100 and 160 μm, respectively.

From this analysis we conclude that in the GOODS-N field the combined PEP/GOODS-H data are 3.0 (1.16) and 2.8 (1.16) times deeper than the PEP (GOODS-H) data at 100 and 160 μm, respectively. Similarly, in the GOODS-S-ultradeep subarea covered by both projects, the combined PEP/GOODS-H data are 2.3 (1.5) and 2.0 (1.5) times deeper than the PEP (GOODS-H) data at 100 and 160 μm, respectively.

4.5. Content of the PEP/GOODS-H released package

The PEP/GOODS-H released package⁸ contains the scientific, uncertainty and coverage PACS maps. We reiterate that the GOODS-S 100/160 μm coverages and noise levels are inhomogeneous due to the combination of observations with different layouts (see Sect. 2). In contrast, the coverages and noise levels for GOODS-N 100/160 μm and GOODS-S 70 μm are fairly homogeneous, except for degradation near the edges.

The released package contains PACS blind and prior source catalogues, down to 3σ significance. Because our source extraction methods might be inaccurate on the noisy edge of the PACS maps, we crop these regions from our catalogues. In order to track coverage variations across the fields, we provide for each source its exposure time (or equivalent) in each passband. The completeness and contamination levels inferred from the MC simulations are part of the released package. To use them, users should restrict their sample to sources with exposure time greater than that quoted in the second column of Table 2. Finally, we remind users that because our flux uncertainties do not fully account for confusion noise, sources with flux densities below 0.6  mJy at 100 μm and 2.0  mJy at 160 μm have to be treated with caution.

Calibration factors used to generate the final PACS maps are derived assuming an in-band SED of νS_ν =   constant (see the Herschel data handbook for more details). Thus, some moderated colour corrections (i.e., <7%) might need to be applied to our flux density measurements, as in-band SEDs of distant galaxies could be different from those assumed here. However, because these colour corrections depend on the redshift of the source, we decided not to apply any colour correction to the released catalogues.

PACS blind catalogues have been cross-matched to our MIPS-24 μm catalogues using a maximum likelihood analysis (Ciliegi et al. 2001; Sutherland & Saunders 1992). This method simultaneously accounts for fluxes and positions of MIPS-24 μm counterparts as well as positional errors in both the PACS and MIPS samples. The cross-identification of PACS and MIPS sources is included in the released package. In Appendix A, we compare the PACS blind and prior source catalogues using this cross-identification.

Following the results of Elbaz et al. (2011) and Hwang et al. (2010), we provide for each source of the blind and prior catalogues, its “clean index”. Because this “clean index” is a measure of the number of bright neighbours for a given source in all passbands, it supplies information on the potential contamination of its flux densities. Here, “bright” neighbours are defined as sources brighter than half of the flux density of the source of interest, and closer than 20′′, 6.7′′ and 11′′ at 24, 100 and 160, respectively. Having counted the number of neighbours of a given source (i.e., Neib24, Neib100 and Neib160), its “clean index” is given by $\mathrm{clean\dot{}index}\mathrm{=}\mathrm{Neib}24\mathrm{+}\mathrm{Neib}100\mathrm{\times }\mathrm{10}\mathrm{+}\mathrm{Neib}160\mathrm{\times }\mathrm{100}\mathit{,}$(1)We note that the radius of 20′′ at 24 μm was defined in order to provide a “clean region” to the SPIRE-250 μm galaxy population (FWHM ~ 18″; see Elbaz et al. 2011). Therefore a “clean index” corresponding to Neib24≤ 1 would provide a very conservative selection for accurate PACS flux density estimates. Thus, we recommend the use of this conservative criterion (i.e., Neib24≤ 1; 36% and 33% of the PACS sources in the GOODS-N and GOODS-S fields, respectively) not to select PACS sources with accurate flux densities but only to test whether any scientific results obtained with the full catalogue do not change when this criterion is applied.

We stress that the MIPS-24 μm catalogues released here are slightly different from those released by Magnelli et al. (2011). These new MIPS-24 μm catalogues are generated using, as prior information, newer versions of the IRAC catalogues which are publicly available as part of the GOODS-H data release. These MIPS-24 μm catalogues are identical to those used for the GOODS-H data release⁹ (Elbaz et al. 2011).

Finally, we note that the PEP/GOODS-H released package also contains ancillary products, i.e., the uncertainty and coverage maps, the empirical PSFs used by our source extraction methods, the Vesta PSFs used to measure aperture corrections, the results of the MC simulations and the PACS residuals maps.

5. Confusion noise

Following the formalism of Dole et al. (2003), deep PACS observations might be affected by two different types of confusion due to extragalactic sources. First, the photometric confusion noise due to sources below the detection limit S_lim which produce signal fluctuations (i.e., σ_c) within the beam of the PACS sources, i.e., S_lim/σ_c = q with q = 3 or 5. Second, the density confusion due to the high density of sources above S_lim which increases the probability to miss objects that are blended with bright neighbours. Because our GOODS-S 100 and 160 μm maps are the deepest blank field observations obtained by the Herschel Space Observatory, they are the most suitable datasets to estimate the PACS-100 and -160 μm photometric and density confusion noise. These estimates are presented in this section while the PACS-70 μm confusion noise is derived and discussed in Berta et al. (2011, they find that the GOODS-S PACS-70 μm observations are too shallow to be affected by either type of confusion.

The photometric confusion noise is estimated empirically following the procedure described in Frayer et al. (2006) and Berta et al. (2011). First, we built several GOODS-S maps using only a fraction of the Herschel observations available. These partial-depth maps regularly probe the logarithmic exposure time parameter space from ~ 0.1  h/pixel to ~ 10  h/pixel. Secondly, we performed prior source extractions on these partial-depth maps and produced the corresponding residual maps removing sources with S > S_lim, with S_lim/σ_c = q using q = 5. Thirdly, we measured the total noise σ_T in these residual maps using the procedure described in Sect. 4.1, i.e., σ_T is defined as the pixel dispersion of the residual map convolved with the Vesta PSF. Finally, we estimated σ_c by analyzing the variation of σ_T as a function of exposure time (see Fig. 5). This procedure was iterated until convergence at S_lim/σ_c = 5 was reached.

Noise in partial-depth maps with short exposure time is dominated by the instrumental noise σ_I and therefore decreases as t^-0.5. In contrast, noise in partial-depth maps with long exposure time departs from the t^-0.5 trend and is a combination of instrumental noise and confusion noise. Assuming that both components follow a Gaussian distribution, σ_T can then be approximated by ${\mathit{\sigma }}_{\mathrm{T}}\mathrm{=}\sqrt{{\mathit{\sigma }}_{\mathit{I}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{c}}^{\mathrm{2}}}$. Because σ_c does not depend on the exposure time, one can thus estimate σ_c by fitting σ_T with a two components function, i.e., σ_c and σ_I ∝ t^-0.5.

We find a photometric confusion noise σ_c of 0.15 and 0.68  mJy in the PACS-100 and -160 μm passbands, respectively. The total noise of the PEP/GOODS-H observations is thus nearly dominated by the photometric confusion noise at 100 μm (${\mathit{\sigma }}_{\mathrm{I}}^{\mathrm{100}}\mathrm{=}\mathrm{0.11}\hspace{0.17em}$mJy) and it is fully dominated by the photometric confusion noise at 160 μm (${\mathit{\sigma }}_{\mathrm{I}}^{\mathrm{160}}\mathrm{=}\mathrm{0.21}\hspace{0.17em}$mJy). In both bands, this significant contribution of the confusion noise to the total noise affects our flux uncertainty estimates: our flux uncertainties (1σ ~ 0.18 and 0.43 mJy at 100 and 160 μm, respectively) do not fully account for confusion noise and sources with flux densities below 0.6  mJy at 100 μm and 2.0  mJy at 160 μm have to be treated with caution (see Sect. 4.4).

Fig. 5 Noise in the PACS-100 μm (top panel) and PACS-160 μm (bottom panel) maps as a function of exposure time. Empty circles represent the total noise σT of the maps. The total noise is fitted with two noise components added in quadrature: an instrumental noise σI component following a t-0.5 trend (dashed line and empty squares) and an constant confusion noise σc component. The dotted lines present the two components fitted to σT. The dot-dashed lines present the two components fitted to σT made in Berta et al. (2011) and illustrate the smaller exposure time range probed in their study.

Fig. 6 PACS 100 and 160 μm differential number counts, normalized to the Euclidean slope (dN/dS ∝ S-2.5). Filled and open symbols show flux density bins above and below the 80% completeness limit, respectively. Grey shaded areas present estimates obtained using pre-Herschel observations. Blue shaded areas present estimates obtained using PACS observations (Berta et al. 2011). Lines represent predictions from backwards or forwards evolutionary models (Lagache et al. 2004; Rowan-Robinson 2009; Valiante et al. 2009; Le Borgne et al. 2009; Franceschini et al. 2010; Gruppioni et al. 2010; Lacey et al. 2010; Marsden et al. 2011; Rahmati & van der Werf 2011; Niemi et al. 2012; Béthermin et al. 2012).

We note that our photometric confusion noise estimates are lower than those of Berta et al. (2011, σ_c ~ 0.27 and 0.92  mJy at 100 and 160 μ. These discrepancies likely come from the fact that the observations used in Berta et al. (2011) were not deep enough to be dominated by the photometric confusion noise and therefore led to more uncertain estimates. This limitation is illustrated by the dot-dashed lines in Fig. 5 which represent the range of exposure time probed in Berta et al. (2011).

The confusion due to the high density of bright sources is usually defined as the flux limit S_lim at which the source density corresponds to 16.7 beams/source, i.e., the density at which 10% of the sources are separated by less than 0.8 × FWHM and thus “blended” (Dole et al. 2003, where Ω = 1.14 × FWHM² is the area of a beam). Using the PACS number counts, Berta et al. (2011) found that this density criterion corresponds to S_lim ~ 2.0 and ~8 mJy at 100 and 160 μm, respectively. These limits are thus greater than the detection limits derived from our MC simulations. This discrepancy can be explained by the very conservative approach adopted by Dole et al. (2003). The density limit of 16.7 beams/source of Dole et al. (2003) corresponds to 10% of blended sources (i.e., separated by less than 0.8 × FWHM). This requirement translates into a blending-completeness of 90% while sources can be accurately extracted at lower completeness: at the 3σ level a catalogue has a typical completeness value of 40–60% (i.e., 60–40% of blended sources). Using this more realistic assumption (i.e., measuring the density for which 40% of the sources are separated by less than 0.8 × FWHM), we infer a confusion density limit of ~3.5 beams/source. In our GOODS-S-ultradeep catalogues, the density of sources corresponds to 8 and 4 beams/source at 100 and 160 μm, respectively. In these catalogues, the density of sources is thus much higher than the density confusion limit defined by Dole et al. (2003, i.e., 16.7 beams/source) but lower than those defined using more realistic assumptions (i.e., 3.5 beams/source). We stress that while an analytical estimate of the confusion density limit is useful, it should be used with caution as it does not account for the fact that our ability to separate pairs of sources depends on their S/N. In contrast, empirical estimates through MC simulations take this effect into account.

6. Number counts

Using the blind catalogues (i.e., without applying any “clean index” selection) we inferred the PACS 100 and 160 μm differential number counts down to an unprecedented depths (PACS-70 μm differential number counts are derived in Berta et al. 2011). For that purpose, we applied the method described in Berta et al. (2010, 2011), which accounts for the incompleteness and contamination of our catalogues using results from the MC simulations (see also Chary et al. 2004; Smail et al. 1995). In this method, observations with different depths cannot be treated simultaneously. Therefore we divided our PACS observations into two sub-samples: (i) a ultradeep sub-sample containing sources in the GOODS-S-ultradeep field; and (ii) a deep sub-sample containing sources in the GOODS-S-deep and GOODS-N fields (these two fields having the same depth). The ultradeep sub-sample covers an effective area of 47  arcmin², while the deep sub-sample covers an effective area of 327  arcmin².

The PACS 100 and 160 μm differential number counts inferred from these two sub-samples are presented in Fig. 6 and provided in Table A.1. These differential number counts are normalized to the Euclidean slope (i.e., dN/dS ∝ S^-2.5), expected for a uniform distribution of galaxies in Euclidean space. Error bars include Poisson statistics, flux calibration uncertainties, and photometric uncertainties, as described in Berta et al. (2011). Filled and open symbols show flux density bins above and below the 80% completeness limit, respectively. Our differential number counts are compared with various estimates from the literature based on Spitzer or Herschel observations (for more details see Berta et al. 2011). There is good agreement between all these estimates over the flux density range in common. However, thanks to the use of deeper PACS observations, our differential number counts extend to fainter flux densities than any previous estimates. This extension of ~ 0.5  dex and ~ 0.2  dex at 100 and 160 μm, respectively, allows us to resolve into individual galaxies an even larger fraction of the cosmic infrared background (CIB). Using Fig. 10 and CIB estimates of Berta et al. (2011,i.e., $\mathrm{12.6}{\mathrm{1}}_{-1.74}^{\mathrm{+}\mathrm{8.31}}$ and $\mathrm{13.6}{\mathrm{3}}_{-0.85}^{\mathrm{+}\mathrm{3.53}}$ nW m^-2 sr^-1 at 100 and 160 μ, we find that in the GOODS-S-ultradeep field our PACS observations resolve $\mathrm{~}\mathrm{7}{\mathrm{5}}_{-30}^{\mathrm{+}\mathrm{12}}\mathrm{\%}$ and $\mathrm{~}\mathrm{7}{\mathrm{5}}_{-15}^{\mathrm{+}\mathrm{7}}\mathrm{\%}$ of the CIB at 100 and 160 μm, respectively. Furthermore, from the wealth of ancillary data available for the GOODS-S and -N fields (see Sect. 7), we can also study the redshift distribution of the PACS sources. Faint PACS 100 μm sources (i.e., S₁₀₀  [mJy]  < 1.5) have a median redshift of $\mathit{z}\mathrm{=}\mathrm{1.3}{\mathrm{7}}_{-0.52}^{\mathrm{+}\mathrm{0.58}}$ (errors give the interquartile range), while brighter sources (i.e., S₁₀₀  [mJy]  > 1.5) have a median redshift of $\mathit{z}\mathrm{=}\mathrm{0.8}{\mathrm{5}}_{-0.33}^{\mathrm{+}\mathrm{0.41}}$. Similarly, faint PACS 160 μm sources (i.e., S₁₆₀  [mJy]  < 2.5) have a median redshift of $\mathit{z}\mathrm{=}\mathrm{1.2}{\mathrm{2}}_{-0.41}^{\mathrm{0.68}}$ while brighter sources (i.e., S₁₆₀  [mJy]  > 2.5) have a median redshift of $\mathit{z}\mathrm{=}\mathrm{0.9}{\mathrm{4}}_{-0.38}^{\mathrm{+}\mathrm{0.52}}$. At 100 μm (160 μm), sources at 0 < z < 0.5, 0.5 < z < 1.0, 1.0 < z < 2.0 and z > 2.0 contribute 24 (17)%, 36 (33)%, 26 (28)% and 14 (22)% of the CIB resolved by our PACS observations.

Fig. 7 (Top panel) Uncertainties in determining LIR from monochromatic observations (i.e., MIPS-24 μm, PACS-70 μm, PACS-100 μm or PACS-160 μm bands) and the Dale & Helou (2002) SED library. These uncertainties are derived by taking the standard deviation of the log(LIR) distribution provided when normalizing all Dale & Helou (2002) SED templates to the same monochromatic flux density (± 20% to account for typical photometric uncertainties, i.e., S/N ~ 5) at a given observed wavelength (i.e., this observed wavelength depends on the band and redshift of interest). (Bottom panel) Fraction of PACS sources detected (i.e., S/N > 3) in only one of our PACS passbands (i.e., not with a 70 + 100  μm, 70 + 160  μm, 100 + 160  μm or 70 + 100 + 160  μm detections) as a function of redshift.

The PACS 100 and 160 μm counts are finally compared to predictions from backwards or forwards evolutionary models (Lagache et al. 2004; Rowan-Robinson 2009; Valiante et al. 2009; Le Borgne et al. 2009; Franceschini et al. 2010; Gruppioni et al. 2010; Lacey et al. 2010; Marsden et al. 2011; Rahmati & van der Werf 2011; Niemi et al. 2012; Béthermin et al. 2012). Although most of the models reproduce the observed PACS 100 and 160 μm number counts fairly well, some of them can be ruled out because of significant discrepancies with our estimates. In particular, we note that the models of Marsden et al. (2011), Valiante et al. (2009) and Niemi et al. (2012) cannot reproduce the steep faint-end slope of the observed 100 and 160 μm counts. From this cursory comparison, it is clear that ultradeep PACS number counts allow for a better refinement of the models and thus for better constraints on the evolution of star-forming galaxies.

7. The infrared luminosity function

Deep Herschel observations give us the opportunity to determine the infrared luminosity function (LF) of galaxies with an unprecedented accuracy. Indeed, far-infrared observations provide more accurate infrared luminosity estimates than mid-infrared observations from Spitzer: for sources with multiple far-infrared detections (i.e., ~70% of sources in the PEP/GOODS-H fields; see Fig. 7), the SED-shape-L_IR degeneracy is broken and infrared luminosity estimates are only limited by photometric uncertainties; for sources with only one far-infrared detection, uncertainties on the monochromatic-to-L_IR conversion are significantly reduced compared to those provided by single mid-infrared detection (see Fig. 7). The importance of far-infrared data increases further in cases where the mid-infrared may be contaminated by the emission from active galaxy nuclei (AGN).

Taking advantage of PACS and SPIRE far-infrared observations for several multi-wavelength fields, Gruppioni et al. (2013) derived the infrared LF of galaxies. Thanks to the large area covered by their observations (~2.5 deg²), they were able to robustly constrain the intermediate and bright-end part of the infrared LF up to z ~ 2 and the bright-end of the infrared LF up to z ~ 4. Here, we extend such study down to unprecedented depths using deeper PACS observations (~2 deeper in term of flux density than those of Gruppioni et al.). From these deeper observations, we are able to better constrain the faint-end and intermediate part of the infrared LF up to z ~ 2 and therefore to obtain better constraints on its redshift evolution.

The GOODS-N and -S fields benefit from an extensive multi-wavelength coverage necessary to obtain redshift information for the PACS sources. In the GOODS-N field, we use a z + K bands selected PSF-matched catalogue created for the PEP survey¹⁰ (Berta et al. 2010, 2011), with photometry in 16 bands and a collection of spectroscopic redshifts (mainly from Barger et al. 2008). In the GOODS-S field, we use the GOODS-MUSIC z + K bands selected PSF-matched catalogue (Grazian et al. 2006; Santini et al. 2009), with photometry in 15 bands and a collection of spectroscopic redshifts. These multi-wavelength catalogues also include photometric redshift estimates computed using all optical and near-infrared data available (see Berta et al. 2011; Grazian et al. 2006; Santini et al. 2009). The quality of these photometric redshifts was tested by comparing them with the redshifts of spectroscopically confirmed galaxies. The high quality of these photometric redshifts is characterized by a relatively small scatter in Δz/(1 + z) of 0.04 and 0.06 in the GOODS-N and -S fields, respectively (Berta et al. 2011; Santini et al. 2009).

Multi-wavelength catalogues were cross-matched with our MIPS-PACS catalogues using their IRAC positions and a matching radius of 0.8″ (i.e., approximately the HWHM of the IRAC-3.6 μm observations). In case of multiple associations (~ 10%), we select the closest optical counterparts. In the GOODS-N and GOODS-S fields, the common area covered by these catalogues is 164  arcmin² and 184  arcmin², respectively. In these regions, 97% and 96% of the PACS sources have a multi-wavelength counterpart. Among those sources, 64% and 61% have a spectroscopic redshift in the GOODS-N and -S fields, respectively. The rest of the sources has photometric redshift estimates.

The total infrared luminosities (8–1000 μm) of PACS sources with redshift estimates were inferred by fitting their far-infrared flux densities (i.e., 70, 100 and 160 μm) with the SED template library of Dale & Helou (2002), i.e., leaving the normalization of each SED template as a free parameter. For sources with only one far-infrared detection, infrared luminosities were defined as the geometric mean across the range of infrared luminosities given by all SED templates. As shown on Fig. 7, even in this case of only one far-infrared detection, the uncertainties in the inferred infrared luminosities are small, i.e., better than ~0.2 dex. To use the MIPS-24 μm flux density of the PACS sources during our fitting procedure does not change our results. Indeed, the ${\mathit{L}}_{\mathrm{IR}}^{\mathrm{PACS}\mathrm{+}\mathrm{MIPS}}\mathit{/}{\mathit{L}}_{\mathrm{IR}}^{\mathrm{PACS}}$ distribution has a mean value of 1 and a dispersion of 3%. Naturally, for sources with only one far-infrared detection, the addition of the MIPS-24 μm flux densities allow us to break the SED-shape-L_IR degeneracy and thus to reduce uncertainties on our infrared luminosity estimates. However, because the fraction of PACS sources with single far-infrared detection is low (i.e., ~30%; see Fig. 7) and because the MIPS-24 μm flux density might be affected by emission from an AGN, we decided not to use the MIPS-24 μm flux densities to derive L_IR. We note that using the SED template library of Chary & Elbaz (2001) instead of that of Dale & Helou (2002) to derive L_IR (again leaving the normalisation as a free parameter), has no impact on our results. Indeed, the ${\mathit{L}}_{\mathrm{IR}}^{\mathrm{DH}}\mathit{/}{\mathit{L}}_{\mathrm{IR}}^{\mathrm{CE}\mathrm{01}}$ distribution has a mean value of 1 and a dispersion of 13%. Figure 8 illustrates the detection limits of our PACS samples in term of L_IR as a function of redshift.

Uncertainties in determining L_IR from monochromatic observations (i.e., Fig. 7) are inferred assuming that the Dale & Helou (2002) SED library reproduces both the full range of models appropriated for star-forming galaxies and the correct distribution of SEDs within this population. Because neither of these assumptions are necessarily true, the absolute values of these uncertainties should be taken with caution. However, even with these limitations, uncertainties derived here are fully consistent with those inferred in Elbaz et al. (2011) by analyzing the mid-to-far-infrared SEDs (i.e., based on Spitzer, PACS and SPIRE observations) of a large sample of star-forming galaxies at 0 < z < 2.

Fig. 8 Infrared luminosities as a function of redshift for Herschel sources situated in the GOODS-S-ultradeep field (red squares) and situated in the GOODS-S-deep and GOODS-N fields (open circles). The blue line on a white background shows the LIR limits above which the GOODS-S-ultradeep sample could be considered as a unbiased sample of the star-forming galaxy population at this redshift. These limits are inferred in Fig. 9, and “steps” correspond to the redshift bins used for the LF analysis. For the GOODS-S-deep and GOODS-N fields, these limits are shifted by ~0.2  dex towards higher LIR.

The infrared (IR) LFs were derived using the standard 1/V_max method (Schmidt 1968). The comoving volume of a given source is defined as V_max = V_zmax − V_zmin where z_min is the lower limit of the redshift bin being used, and z_max is the maximum redshift at which the source could be seen given the flux density limits of our observations, with a maximum value corresponding to the upper limit of the redshift bin. Here z_max was defined by redshifting the Dale & Helou (2002) template fitted to the far-infrared flux densities of our sources until it fell below the detection limits of our PACS observations, or until z_max is greater than the upper limit of our redshift bin. For each luminosity bin, the LF is then given by $\mathit{\phi }\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Delta }\mathit{L}}\begin{array}{c}\sum \end{array}\frac{\mathrm{1}}{{\mathit{V}}_{\max \mathit{,i}}\mathrm{\times }{\mathit{w}}_{\mathit{i}}}\mathit{,}$(2)where V_max,i is the comoving volume over which the ith galaxy of the luminosity bin could be observed, ΔL is the size of the luminosity bin, and w_i is the completeness correction factor of the ith galaxy. The value of w_i is given by the MC simulations (Sect. 4.4) and depends on the flux densities of each source: w_i equals 1 for bright PACS sources and decreases at faint flux densities due to the incompleteness of our PACS catalogues. For sources with multiple PACS detections it is defined as the maximum of the three passbands, i.e., ${\mathit{w}}_{\mathit{i}}\mathrm{=}\max \mathrm{\left(}{\mathit{w}}_{\mathit{i}}^{\mathrm{70}}\mathit{,}{\mathit{w}}_{\mathit{i}}^{\mathrm{100}}\mathit{,}{\mathit{w}}_{\mathit{i}}^{\mathrm{160}}\mathrm{\right)}$. Because we limit our LFs to infrared luminosities where w_i > 0.5, none of our results strongly depend on these corrections.

Fig. 9 Selection limits introduced in the Tdust−LIR parameter space by our deepest PACS observations, i.e., in GOODS-S-ultradeep. Red continuous lines are these selection limits at different redshifts using our PACS 5σ detection limits. Each line corresponds to the central redshift of our redshift bins, i.e., z = 0.25 (further left), z = 0.55, z = 0.85, z = 1.15, z = 1.55 and z = 2.05 (further right). The shaded area shows the local Tdust−LIR relation found by Chapman et al. (2003), linearly extrapolated to 1013 L⊙. Dot-dashed lines show the Tdust−LIR relation inferred by Symeonidis et al. (2013) using a sample of high-redshift (i.e., 0.2 < z < 1.2) Herschel-detected galaxies. Dashed black lines show, for each redshift, the lowest infrared luminosities probed by our ultradeep PACS observations without any dust temperature biases and yet populated by star-forming galaxies.

Fig. 10 Infrared luminosity functions estimated in six redshift bins with the 1/Vmax method. Red squares and black triangles show results from our ultradeep (i.e., GOODS-S-ultradeep) and deep (i.e., GOODS-S-deep and GOODS-N) samples, respectively. Dashed lines represent the best fits to these data points with a double power-law function with fixed slopes. The shaded areas span all the solutions of these fits which are compatible, within 1σ, with our data points: the dark shaded parts of these areas highlight the luminosity ranges directly constrained by our PACS observations, while the light shaded parts highlight the luminosity ranges where our constraints rely upon extrapolations based on φknee and Lknee and a faint-end slope fixed at its z ~ 0 value. Asterisks show the local reference, taken from Sanders et al. (2003), and the dotted line is the best fitted to these data points with our double power-law function with fixed slopes. Red dot-dashed lines are results from Magnelli et al. (2009, 2011) using deep MIPS-24 μm observations. Magnelli et al. (2009, 2011) used the same double power-law function that is used here. To illustrate the infrared luminosity range constrained by Magnelli et al. (2009, 2011), we show, as open circles, their lowest infrared luminosity bins. Blue triple-dot-dashed lines present results from Gruppioni et al. (2013) who used a different analytical function to fit their data points, in particular, a much shallower faint-end slope.

The minimum infrared luminosities that can be probed by the IR LFs in each of our redshift bins (i.e., 0.1 < z < 0.4, 0.4 < z < 0.7, 0.7 < z < 1.0, 1.0 < z < 1.3, 1.3 < z < 1.8 and 1.8 < z < 2.3) depend on the depth of our observations: at a given infrared luminosity, a large fraction of the galaxies has to be observable (i.e., small completeness correction; w_i > 0.5) over at least half of our redshift bin (i.e., small V_max correction; ${\mathit{z}}_{\max }\mathit{>}{\mathit{z}}_{\mathrm{bin}}^{\mathrm{lower}}\mathrm{+}\mathrm{(}{\mathit{z}}_{\mathrm{bin}}^{\mathrm{upper}}\mathrm{-}{\mathit{z}}_{\mathrm{bin}}^{\mathrm{lower}}\mathrm{)}\mathit{/}\mathrm{2}$). From this definition, it is clear that observations with different depths cannot be treated simultaneously. Therefore, as for the number counts, we first divided our PACS sources into two sub-samples: (i) a ultradeep sub-sample from the GOODS-S-ultradeep field; and (ii) a deep sub-sample from the GOODS-S-deep and GOODS-N fields. Then, for each of our redshift bins, we estimated the lowest infrared luminosities probed by these sub-samples.

At a given redshift, the minimum infrared luminosity observable by PACS depends on the dust colour temperature of galaxies: at a given L_IR, galaxies with warmer dust have brighter PACS flux densities. Therefore, we estimated, for each point of the T_dust−L_IR parameter space its detectability by our PACS observations using the SED templates of Dale & Helou (2002), i.e., a dust temperature was assigned to each Dale & Helou template in a manner that is consistent with procedures used in Chapman et al. (2003) to derive the L_IR − T_dust relation. Then, assuming that the local T_dust−L_IR correlation of Chapman et al. (2003) remains the same at high-redshift (see also Chapin et al. 2009), we defined the minimum L_IR of our IR LFs as the minimum L_IR observable by PACS without any T_dust biases and populated by star-forming galaxies (i.e., within the T_dust−L_IR correlation). In this analysis we used the 5σ limits of our PACS observations (i.e., where w_i > 0.5) and the central redshift of our redshift bins (i.e., ${\mathit{z}}_{\max }\mathit{>}{\mathit{z}}_{\mathrm{bin}}^{\mathrm{lower}}\mathrm{+}\mathrm{(}{\mathit{z}}_{\mathrm{bin}}^{\mathrm{upper}}\mathrm{-}{\mathit{z}}_{\mathrm{bin}}^{\mathrm{lower}}\mathrm{)}\mathit{/}\mathrm{2}$). Figure 9 presents the results of this analysis for our ultradeep sub-sample. The minimum L_IR observable by PACS strongly increases with increasing redshift, and sources with hotter dust temperatures can be detected down to fainter infrared luminosities (red lines of Fig. 9). Combined with the expected positions of galaxies in the T_dust−L_IR parameter space, we can define the minimum L_IR of our IR LFs (dashed lines of Fig. 9). The same analysis was performed for our deep sub-sample. Results are very similar but systematically shifted by ~ 0.2  dex towards higher L_IR. We note that the modest evolution of the T_dust−L_IR relation with redshift (see dot-dashed lines of Fig. 9¹¹; Symeonidis et al. 2013; Hwang et al. 2010; Magnelli et al., in prep.) only has a minor effect on these limits, i.e., they are shifted by 0.05 − 0.1  dex towards higher luminosities. We also note that using the dust temperature of galaxies situated on the main sequence (MS) of the SFR-M_∗ plane, we would derive similar limits to those obtained using the T_dust−L_IR relation. Indeed, MS galaxies have T_dust of 27 ± 3, 28 ± 3, 29 ± 3, 30 ± 4, 32 ± 4 and 34 ± 5  K at z = 0.25, z = 0.55, z = 0.85, z = 1.15, z = 1.55 and z = 2.05 (Magnelli et al., in prep.; see also Magdis et al. 2012).

Uncertainties in the IR LFs depend on the number of sources per luminosity bin, on the photometric redshift errors and on the infrared luminosity errors. For the plot, they were defined as the quadratic sum of the Poissonian errors (∝$\mathrm{1}\mathit{/}\sqrt{\mathit{N}}$) and errors computed from MC simulations which account for both photometric redshift and infrared luminosity uncertainties. The methodology of these MC simulations is described in Magnelli et al. (2009, 2011).

Fig. 11 Comparisons between the infrared luminosity functions derived using PACS observations and those derived using MIPS observations. PACS-based IR LFs are shown as dark/light shaded areas and dashed lines (see also Fig. 10). MIPS-based IR LFs are derived using MS-based 24 μm-to-LIR conversion factors (black triangles) and above-MS-based 24 μm-to-LIR conversion factors (red squares; see text for more details). The dotted lines show the local reference taken from Sanders et al. (2003). Blue triple-dot-dashed lines present results from Gruppioni et al. (2013). At low luminosities, the agreement between PACS-based and MIPS-based IR LFs confirms that PACS extrapolations towards lower infrared luminosities (i.e., light shaded areas) are reliable, at least down to the infrared luminosities probed by Spitzer.

Figure 10 represents the IR LFs derived in six redshift bins using our ultradeep and deep PACS observations (Tables A.2 and A.3). We fit the IR LFs with a double power-law function similar to that used to fit the local IR LF which we also plot for reference (Sanders et al. 2003, φ ∝ L^-0.6 for log (L/L_⊙) < L_knee and φ ∝ L^-2.2 for log (L/L_⊙) > L_knee). In this fitting procedure, the normalization (i.e., φ_knee) and the transition luminosity (i.e., L_knee) of the double power-law function are left as free parameters. The shaded areas of Fig. 10 present the solutions compatible with the data within 1σ. The evolution of φ_knee and L_knee with redshift is presented in the upper left panel of Fig. 12 and given in Table A.4.

We compare our IR LFs with estimates made by Magnelli et al. (2009, 2011) using deep MIPS-24 μm observations and Herschel-based estimates from Gruppioni et al. (2013) that use shallower PACS observations covering a wider effective area. At z < 1.8, we find good agreement with the Spitzer analysis of Magnelli et al. (2009, 2011). In particular, we observe good agreement, within the uncertainties, between our low luminosity extrapolations (based on φ_knee and L_knee and a faint-end slope fixed at its z ~ 0 value) and direct constraints obtained by Magnelli et al. (2009, 2011, empty circles of Fig. 10. This agreement shows that even though our deepest PACS data do not allow us to probe luminosities far below the “knee” of the IR LFs, we obtain accurate low luminosity extrapolations, at least down to the infrared luminosities probed by Spitzer. The most significant difference between the present results and those of Magnelli et al. (2009, 2011) is observed at z ~ 2 and at the highest infrared luminosities (L_IR > 10¹² L_⊙). There, our new IR LF has a higher normalization. Because, the MIPS-24 μm band can be affected by a significant contribution from an AGN, Magnelli et al. (2009, 2011) excluded from their sample all X-ray AGNs, i.e., sources detected in the GOODS-S/N Chandra catalogues (Alexander et al. 2003; Lehmer et al. 2005) with either L_X [0.5 − 8.0 keV]  > 3 × 10⁴²  erg s^-1 or a hardness ratio greater than 0.8 (Bauer et al. 2004). We excluded in the same manner all X-ray AGNs from our sample and found that these exclusions could not reconcile our IR LFs at z ~ 2. The observed difference is thus likely due to the fact that at z ~ 2, estimates from Magnelli et al. (2009, 2011) are affected by large pre-Herschel uncertainties on the 24 μm-to-L_IR conversion factors (Nordon et al. 2010; Elbaz et al. 2010, 2011; Nordon et al. 2012,and reference therein). At this redshift, PACS estimates are more accurate.

At all redshifts we observe very good agreement between our estimates and those of Gruppioni et al. (2013). At high luminosities this agreement is very encouraging, because the data analysed by Gruppioni et al. (2013) sample a larger volume at shallower flux limits, and thus can more accurately measure the number densities for rare, bright sources. The only disagreement between our results and those of Gruppioni et al. (2013) appear at very low infrared luminosities, i.e., at luminosities not probed by both studies and which depend on the analytic function used to fit the data. Gruppioni et al. (2013) used a faint end slope shallower than that used here, i.e., − 0.2 compared with − 0.6.

We find an evolution of φ_knee (φ_knee = 10^{−2.57 ± 0.12} × (1 + z)^{−1.5 ± 0.7} for z < 1.0 and φ_knee = 10^{−2.03 ± 0.72} × (1 + z)^{−3.0 ± 1.8} for z > 1.0) and L_knee (L_knee = 10^{10.48 ± 0.10} × (1 + z)^3.8 ± 0.6 for z < 1.0 and L_knee = 10^{10.31 ± 0.47} × (1 + z)^4.2 ± 1.2 for z > 1.0) in broad agreement with what was observed in Magnelli et al. (2011) and Gruppioni et al. (2013). However, we note that in our study we find a stronger evolution of L_knee at z > 1.0 than in Magnelli et al. (2011).

With our ultradeep PACS data we can only constrain the IR LFs of galaxies down to L_IR = 10¹¹ L_⊙ at z ~ 1 and L_IR = 10¹² L_⊙ at z ~ 2. Consequently, at lower infrared luminosities, our constraints only rely on extrapolations based on φ_knee and L_knee and a faint-end slope fixed at its z ~ 0 value. To test these extrapolations, we take advantage of deep MIPS-24 μm observations using appropriate 24 μm-to-L_IR conversion factors. Thanks to recent Herschel studies, we know that the 24 μm-to-L_IR conversion factors depend not only on the MIPS-24 μm luminosity (as, e.g., in the Chary & Elbaz SED library) but also on the localization of galaxies with respect to the “main sequence” of star formation (MS; log (SFR)= α × log  (M_∗) + C(z), where 0.5 < α < 1.0; Elbaz et al. 2011; Nordon et al. 2012): galaxies situated on the MS have different 24 μm-to-L_IR conversion factors than galaxies situated above the MS (Elbaz et al. 2011; Nordon et al. 2012). Consequently, we derive MIPS-based IR LFs using, on one hand, MS-based 24 μm-to-L_IR conversion factors and on the other hand, above-MS-based 24 μm-to-L_IR conversion factors (Elbaz et al. 2011, Fig. 11). The real IR LFs should be situated between these two estimates, but with a faint-end closer to the MS-based estimates (i.e., in a luminosity range dominated by MS galaxies) and a bright-end closer to the above-MS-based estimates (i.e., in a luminosity range dominated by above-MS galaxies). Up to z ~ 1.3, MS-based and above-MS-based estimates are in agreement: in this redshift range, MS-based and above-MS-based 24 μm-to-L_IR conversion factors are not significantly different. In contrast, at z > 1.3, above-MS-based IR LFs have higher normalization than MS-based estimates: at z > 1.3, above-MS-based 24 μm-to-L_IR conversion factors are much larger than those for MS galaxies (Elbaz et al. 2011; Nordon et al. 2012). At z ~ 2 and at high infrared luminosities, the agreement found between PACS-based and above-MS-based estimates confirms the assumption that the bright-end of the IR LF is dominated by above-MS galaxies. At all redshifts and at low infrared luminosities (i.e., where MS galaxies should dominate the IR LFs), MS-based estimates are in perfect agreement with PACS extrapolations towards lower infrared luminosities (i.e., dashed lines and light shaded areas). We note that these PACS extrapolations also agree, within the error bars, with the above-MS-based estimates. All these agreements confirm that our PACS extrapolations are reliable, at least down to the infrared luminosities probed by Spitzer.

By integrating our PACS-based IR LFs we derive the evolution of the comoving number density (top right panel of Fig. 12) and comoving infrared luminosity density (IR LD, bottom panels of Fig. 12 and Table A.5) of “faint” galaxies (i.e., 10⁷ L_⊙ < L_IR < 10¹¹ L_⊙), LIRGs (i.e., 10¹¹ L_⊙ < L_IR < 10¹² L_⊙) and ULIRGs (i.e., L_IR > 10¹² L_⊙). Here, one has to keep in mind that even with our ultradeep PACS data we can only constrain the IR LFs of galaxies down to L_IR = 10¹¹ L_⊙ at z ~ 1 and L_IR = 10¹² L_⊙ at z ~ 2, respectively. Therefore, number densities and IR LDs of galaxies below these limits rely upon extrapolations based on φ_knee and L_knee and a faint-end slope fixed at its z ~ 0 value. The redshift and luminosity ranges for which the inferred number densities and IR LDs rely upon these extrapolations are highlighted by dashed lines and striped regions, respectively. Although these extrapolations seem to be corroborated by direct constraints from Spitzer (see Fig. 11), we recommend caution when interpreting values not directly constrained by our PACS observations. We also emphasize that the LIRG and ULIRG designations are used here strictly to segregate the luminosity bins, but not to imply physical properties. Indeed, Herschel studies have unambiguously revealed that high-redshift (U)LIRGs do not have the same properties as their local counterparts (e.g., Elbaz et al. 2011; Wuyts et al. 2011; Nordon et al. 2012; Magnelli et al. 2012).

Fig. 12 (Top left panel) Evolution of φknee and Lknee as a function of redshift. (Top right panel) Evolution of the comoving number density of “faint” galaxies (i.e., 107 L⊙ < LIR < 1011 L⊙; black triangles), LIRGs (1011 L⊙ < LIR < 1012 L⊙; orange diamonds) and ULIRGs (LIR > 1012 L⊙; red stars). Green circles show the number density of galaxies which are above the PACS detection limit of our ultradeep sub-sample (i.e., >${\mathit{L}}_{\mathrm{IR}}^{\min }$, defined in Fig. 9). The local reference is taken from Sanders et al. (2003). (Bottom left panel) Evolution of the total comoving IR energy density (black area) and the relative contribution of the “faint” galaxies (light blue area), LIRGs (orange area) and ULIRGs (red area). The black dot-dashed line shows the best fit of the total comoving IR energy density, i.e., defined by integrating the IR LFs best-fitting our data points (see dashed lines in Fig. 10). The green dashed line shows the best fit of the comoving IR energy density of galaxies which are above the PACS detection limit of our ultradeep sub-sample (i.e., >${\mathit{L}}_{\mathrm{IR}}^{\min }$). Black circles represent the total comoving IR energy density inferred by Gruppioni et al. (2013). The right axis of the panel shows the evolution of the obscured SFR density assuming that SFR   [M⊙  yr-1]  = 10-10 × LIR  [L⊙] for a Chabrier (2003) IMF. (Bottom right panel) Evolution of the comoving IR energy density. Areas are the same as in the bottom left panel. The dotted lines represents the unobscured SFR density of the Universe (i.e., not corrected for extinction; Cucciati et al. 2012). The red dashed line on a white background shows the total SFR density of the Universe, defined as the sum of the obscured and unobscured SFR densities. Yellow diamonds are a compilation of SFR density estimates from Hopkins & Beacom (2006). Green circles and squares are the total SFR density of the Universe and the contribution of ULIRGs estimated by Murphy et al. (2011) using deep MIPS-24 and -70 μm observations. Dark blue diamonds, stars and triangles are estimates from Casey et al. (2012) for all galaxies, LIRGs and ULIRGs, respectively, using SPIRE (i.e., submillimetre) observations. One has to keep in mind that even with our ultradeep PACS data we only constrain the IR LFs of galaxies down to LIR = 1011 L⊙ at z ~ 1 and LIR = 1012 L⊙ at z ~ 2, respectively. Therefore, number densities and IR LDs of galaxies below these limits rely upon extrapolations based on φknee and Lknee and a faint-end slope fixed at its z ~ 0 value. The redshift and luminosity ranges for which the inferred number densities and IR LDs rely upon these extrapolations are highlighted by dashed lines and striped regions, respectively. We note that these extrapolations are nevertheless corroborated by direct constraints based on deep MIPS-24 μm observations (see Fig. 11).

We find that the number density of ULIRGs strongly evolves from z = 0 to z ~ 2, as it is multiplied by a factor ~800. Similarly, the number density of LIRGs significantly evolves with redshift: using direct constraints from PACS, we find that the LIRGs number density is multiplied by a factor ~25 between z = 0 and z ~ 1, while relying on extrapolations we find that it is multiplied by a factor ~30 between z = 0 and z ~ 2. Naturally, the redshift evolution of the LIRGs and ULIRGs number densities translate into the redshift evolution of their IR LDs. The IR LDs of LIRGs and ULIRGs are multiplied by a factor ~130 (~1000) and ~40 (~45; based on extrapolations) between z = 0 and z ~ 1 (2), respectively. The redshift evolution of the IR LD of ULIRGs is consistent with that found by Murphy et al. (2011) using deep MIPS-24/70 μm observations (green squares in the bottom right panel of Fig. 12). In contrast, at z ~ 2, Magnelli et al. (2011) found an IR LD of ULIRGs lower by a factor ~3. This inconsistency reflects the discrepancies observed at z ~ 2 between their and our IR LF (see Fig. 10).

Relying upon extrapolations to low infrared luminosities, we find that the total IR LD strongly increases from z = 0 and z ~ 2, with most of the evolution happening before z ~ 1.2 driven by the strong increase of the IR LD of LIRGs. At z ~ 1, LIRGs account for 50 ± 26% of the total IR LD while ULIRGs contribute only 10 ± 6% of it. In contrast, at z ~ 2 ULIRGs contribute 50 ± 24% while LIRGs contribute only 30 ± 20% of the total IR LD.

Although the total IR LDs derived here are consistent, within the uncertainties, with those from Magnelli et al. (2011), there are some differences. At z ~ 1, the total IR LD derived here is lower by a factor ~1.3 than that of Magnelli et al. (2011) while at z ~ 2 it is higher by a factor ~1.3. More importantly, at z ~ 2, the ULIRGs contribution to the total IR LD found here (i.e., 50% ± 24) is much higher than that found in Magnelli et al. (2011, ~17%. Naturally, this finding reflects the disagreement at z ~ 2 between our two IR LFs.

The total IR LDs inferred in Gruppioni et al. (2013, black circles in the bottom left panel of Fig. 12 agree with our findings. Noticeably, their and our estimates are obtained by integrating IR LFs which significantly differ at infrared luminosities lower than those probed by Spitzer, i.e., Gruppioni et al. (2013) used a much shallower faint-end slope (see Fig. 11). This demonstrates that extrapolations towards very low infrared luminosities do not affect much the estimates of the total IR LDs. Indeed, unless the faint-end slope of the IR LFs is significantly steeper than that used here, most of the total comoving IR LD is emitted by galaxies with luminosities reachable with PACS (i.e., ~50%; see black dot-dashed and green dashed lines in the bottom left panel of Fig. 12) and Spitzer (i.e., ~75%).

While the IR LDs derived here for all galaxies and for LIRGs are consistent with SPIRE-based results from Casey et al. (2012), the IR LDs of ULIRGs at z > 0.8 derived here are lower than those inferred by Casey et al. (2012). This disagreement might be due to the large spectroscopic redshift incompleteness corrections applied to the IR LFs of Casey et al. (2012) at z > 0.8 (i.e., corrections >50%). This seems to be confirmed by the good agreement, at z > 0.8, between the IR LFs of Gruppioni et al. (2013, i.e., also in agreement with our IR LFs) and the SPIRE-based IR LFs of Vaccari et al. (in prep.). Finally we note that the total IR LD derived here, as well as the significant contribution of ULIRGs at z ~ 2, are consistent with results from Murphy et al. (2011).

Assuming that the IR LD is totally produced by star formation (i.e., without any contribution from AGNs), it can be converted into the obscured SFR density of the Universe using the relation of Kennicutt (1998), scaled to a Chabrier (2003) initial mass function¹² : $\mathit{SFR}\hspace{0.17em}\mathrm{[}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathrm{yr}}^{-1}\mathrm{]}\mathrm{=}{\mathrm{10}}^{-10}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{L}}_{\mathrm{IR}}\hspace{0.17em}\mathrm{\left[}{\mathit{L}}_{\mathrm{\odot }}\mathrm{\right]}\mathit{,}$(3)About 10–15% of our PACS sources are associated with X-ray AGNs that might significantly contribute to the IR LD derived here. However, Herschel studies of X-ray AGNs (Rosario et al. 2013; Shao et al. 2010; Mullaney et al. 2012; Rosario et al. 2012) have demonstrated that in the vast majority of cases (i.e., >94%) the PACS flux densities are dominated by emission from the host galaxy and thus provide an uncontaminated view of their star-formation activities. Consequently, we assume that the IR LD derived here has no significant contribution from AGNs and can be converted into the obscured SFR density of the Universe (right axis of the bottom panels of Fig. 12). Combined with the unobscured SFR density of the Universe derived by Cucciati et al. (2012) using rest-frame UV observations, we can then infer the total SFR density and its evolution up to z ~ 2 (dashed red line on a white background in the bottom right panel of Fig. 12). The cosmic star-formation history strongly evolves from z = 0 to z ~ 1 and flattens at z > 1. The unobscured SFR density accounts for about ~25%, ~12% and ~17% of the total SFR density of the Universe at z ~ 0, z ~ 1 and z ~ 2, respectively. The contribution of the unobscured SFR density and its evolution with redshift is consistent with the redshift evolution of the mean rest-frame UV dust attenuation (e.g., Cucciati et al. 2012; Tresse et al. 2007). Finally, we note that the cosmic star-formation history derived here is fully consistent with the combination of indicators, either obscured or corrected for dust extinction, as compiled by Hopkins & Beacom (2006), Seymour et al. (2008) and Karim et al. (2011, not shown in Fig. 12.

8. Summary

By combining observations of the GOODS fields from the PEP and GOODS-Herschel key programmes, we obtain the deepest PACS far-infrared blank field extragalactic survey carried out by the Herschel Space Observatory. In particular, in the GOODS-S field the combination of these observations is not limited by the exposure time but by confusion. PACS flux densities are extracted from the maps using two complementary PSF-fitting approaches. Firstly, we extract PACS flux densities using, as prior information, the expected positions of the sources on the basis of deep 24 μm catalogues. Secondly, PACS flux densities are extracted “blindly”, i.e., without positional priors. The accuracy of both approaches is tested through MC simulations. In the deepest parts of the GOODS-S field, these catalogues reach 3σ depths of 0.9, 0.6 and 1.3 mJy at 70, 100 and 160 μm, respectively. From these catalogues we derive number counts down to these unprecedented depths, and determine the infrared luminosity function down to L_IR = 10¹¹ L_⊙ at z ~ 1 and L_IR = 10¹² L_⊙ at z ~ 2, respectively. By integrating these infrared luminosity functions, we estimate the evolution of the SFR density of the Universe up to z ~ 2. We find that the obscured SFR density of the Universe strongly increases from z = 0 to z ~ 1 and then increases more moderately up to z ~ 2.3. The obscured SFR density of the Universe is dominated by the contribution of the LIRGs at z ~ 1 (i.e., 50 ± 26%) while it is dominated by the contribution of ULIRGs (50 ± 24%) and LIRGs (30 ± 20%) at z ~ 2.

Maps and source catalogues (>3σ) are now publicly released. Combined with the large wealth of multi-wavelength data available for the GOODS fields, they provide a powerful new tool for studying galaxy evolution over a broad range of redshifts.

Online material

Appendix A: Comparison between blind and prior catalogues

Figure A.1 presents the comparison between the blind and prior catalogues in the GOODS-N and GOODS-S fields. These comparisons are restricted to the region of the field with exposure time higher than those mentioned in the second column of Table 2. For the GOODS-S 100 and 160 μm maps, this restriction corresponds to the central deepest region of the field, i.e., GOODS-S-ultradeep. The cross-identification of PACS blind and prior sources is based on their MIPS-24 μm positions, i.e., for the blind catalogues we use the PACS/MIPS-24 μm cross-identification performed using a maximum likelihood analysis (see Sect. 4.5). For each field and passband we present the direct comparison of flux densities for sources in common to both catalogues; the distribution of unmatched sources in absolute number; and the distribution of unmatched sources in fraction relative to total in the given flux density bin.

Flux densities extracted with these two independent source extraction methods are consistent with each other. At faint flux densities, there are however few outliers, with higher flux densities in the blind catalogues than in the prior catalogues. Such outliers are expected since blind source extraction is more severely affected by blending issues leading to the overestimation of PACS flux densities.

In the GOODS-N field the absolute number of unmatched sources in the blind and prior source catalogues are very

similar. At faint flux densities, the increase of the fraction of unmatched sources is likely explained by the increase of the incompleteness of our catalogues. Indeed, because our two source extraction methods are independent, at faint flux densities their incomplete samples might not fully overlap.

The blind GOODS-S 100 and 160  μm source catalogues contain a larger number of unmatched sources than the prior source catalogues. Some sources in “excess” in the blind catalogues should correspond to sources effectively missed by our prior extraction due to the lack of MIPS-24 μm counterparts (see Sect. 4.1). However, a significant fraction of the sources in “excess” in the blind catalogue should correspond to spurious detections, as in this field the contamination of the blind catalogues is supposed to be higher than that of the prior catalogues (see Table 2). The GOODS-N 100 μm, GOODS-S 100 μm and GOODS-S 160 μm blind catalogues also contain each one a bright source (i.e., > 8  mJy) missed by the prior extraction. Examining these sources, it turns out that they are not missed because of a lack of MIPS-24 μm priors, but likely correspond to “spurious” bright sources created by the blind extraction by over-deblending very bright and crowded regions of the field. Finally, we observe that the prior GOODS-S 70 μm source catalogue contains a larger number of unmatched sources as it reaches a slightly lower 3σ detection limit than the blind catalogue.

Fig. A.1 Comparison between blind and prior source catalogues. (Left panel) Direct comparison of flux densities for sources in common to both catalogues. (Centre panel) Distribution of unmatched sources in absolute number. (Right panel) Distribution of unmatched sources in fraction relative to total in the given flux density bin.

Acknowledgments

We thank the anonymous referee for suggestions which greatly enhanced this work. 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); and 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). Support for this work was provided by the NASA through an award issued by JPL/Caltech

References

All Tables

All Figures

	Fig. 1 PACS-100  μm coverage map of the GOODS-S field. Due to the different layout of the PEP and GOODS-H observations, the 10′ × 10′ region centred on the GOODS-S field has much deeper PACS-100 μm observations than the outskirts. We observe the same pattern in the PACS-160  μm coverage map of the GOODS-S field.
In the text

Fig. 2 (Top left) PEP/GOODS-H map of the GOODS-north field (11′ × 17′) at 100 μm. (Top right) PEP/GOODS-H map of the GOODS-south field (11′ × 17′) at 100 μm. Contours correspond to exposure times greater than 0.5 (at the edge), 3, 6 and 9 (centre) hrs/pix. (Bottom left) Colour composite image of the GOODS-north field at 24 μm (blue), 100 μm (green) and 160 μm (red). (Bottom right) Colour composite image of the GOODS-south field at 24 μm (blue), 100 μm (green) and 160 μm (red). The 24 μm images (PI: M. Dickinson) were obtained by the Spitzer Space Telescope, while the 100 and 160 μm images were obtained by the Herschel Space Observatory. In the colour composite images, sources with contribution from the 24 + 100  μm, 24 + 160  μm, 100 + 160  μm and 24 + 100 + 160  μm bands would correspond to a cyan, magenta, yellow and white colours, respectively. The relatively noisy edges of the colour composite images have been trimmed.

In the text

Fig. 3 PACS-to-24 μm flux density ratio as a function of the PACS flux density inferred from the “blind” source catalogue cross-matched with the MIPS-24 μm catalogue. Shaded regions show the space density distribution of galaxies at 70 μm (top left), 100 μm (top right) and 160 μm (bottom) in the GOODS-S field. In this field, the MIPS 24 μm catalogue reaches a 3σ limit of 20 μJy, while the PACS-70, 100 and 160 μm data reach 3σ limits of 0.9, 0.6 and 1.3 mJy (vertical dotted black lines). The parameter spaces reachable by our prior source extraction method using the MIPS catalogue are located below the dashed black lines: PACS sources located above these lines will not be detected at 24 μm. Orange triangles, red circles, light blue squares and dark blue stars show the PACS-to-24 μm flux density ratios of galaxies with LIR = 1011.5 L⊙ as predicted at different redshifts using the spectral energy distributions (SEDs) of normal (MS), starburst (SB) galaxies (Elbaz et al. 2011), IRAS 22491 (Berta 2005) and Arp220 (Silva et al. 1998), respectively. The right point of each track corresponds to z = 0.5, while other points correspond to increasing redshifts, with intervals of Δz = 0.5 (arrows indicate the path followed for increasing redshift). Predictions for galaxies with lower or higher infrared luminosities will correspond to a global vertical shift of these tracks towards lower or higher PACS flux densities, respectively. Most of the PACS population lies well within the limits of the parameter space reachable by our MIPS 24 μm catalogues, i.e., below the dashed lines.

In the text

Fig. 4 Results of MC simulations in the GOODS-North and -South fields using the prior source extraction method. These MC simulations correspond to the central deepest regions of each field, i.e., avoiding the rather noisy edges of the GOODS-N field and concentrating on the GOODS-S-ultradeep part of the GOODS-S field. Blue lines represent the average photometric accuracy defined as the standard deviation of the (Sout − Sin)/Sout distribution in each flux bin (after 3σ clipping). Red lines show the mean value of the (Sout − Sin)/Sout distribution in each flux bin. Inset plots show the fraction of artificial sources detected in the image (i.e., completeness) as a function of input flux (orange plain histogram) and the fraction of spurious sources (i.e., contamination) as a function of flux density (striped black histogram).

In the text

Fig. 5 Noise in the PACS-100 μm (top panel) and PACS-160 μm (bottom panel) maps as a function of exposure time. Empty circles represent the total noise σT of the maps. The total noise is fitted with two noise components added in quadrature: an instrumental noise σI component following a t-0.5 trend (dashed line and empty squares) and an constant confusion noise σc component. The dotted lines present the two components fitted to σT. The dot-dashed lines present the two components fitted to σT made in Berta et al. (2011) and illustrate the smaller exposure time range probed in their study.

In the text

Fig. 6 PACS 100 and 160 μm differential number counts, normalized to the Euclidean slope (dN/dS ∝ S-2.5). Filled and open symbols show flux density bins above and below the 80% completeness limit, respectively. Grey shaded areas present estimates obtained using pre-Herschel observations. Blue shaded areas present estimates obtained using PACS observations (Berta et al. 2011). Lines represent predictions from backwards or forwards evolutionary models (Lagache et al. 2004; Rowan-Robinson 2009; Valiante et al. 2009; Le Borgne et al. 2009; Franceschini et al. 2010; Gruppioni et al. 2010; Lacey et al. 2010; Marsden et al. 2011; Rahmati & van der Werf 2011; Niemi et al. 2012; Béthermin et al. 2012).

In the text

Fig. 7 (Top panel) Uncertainties in determining LIR from monochromatic observations (i.e., MIPS-24 μm, PACS-70 μm, PACS-100 μm or PACS-160 μm bands) and the Dale & Helou (2002) SED library. These uncertainties are derived by taking the standard deviation of the log(LIR) distribution provided when normalizing all Dale & Helou (2002) SED templates to the same monochromatic flux density (± 20% to account for typical photometric uncertainties, i.e., S/N ~ 5) at a given observed wavelength (i.e., this observed wavelength depends on the band and redshift of interest). (Bottom panel) Fraction of PACS sources detected (i.e., S/N > 3) in only one of our PACS passbands (i.e., not with a 70 + 100  μm, 70 + 160  μm, 100 + 160  μm or 70 + 100 + 160  μm detections) as a function of redshift.

In the text

Fig. 8 Infrared luminosities as a function of redshift for Herschel sources situated in the GOODS-S-ultradeep field (red squares) and situated in the GOODS-S-deep and GOODS-N fields (open circles). The blue line on a white background shows the LIR limits above which the GOODS-S-ultradeep sample could be considered as a unbiased sample of the star-forming galaxy population at this redshift. These limits are inferred in Fig. 9, and “steps” correspond to the redshift bins used for the LF analysis. For the GOODS-S-deep and GOODS-N fields, these limits are shifted by ~0.2  dex towards higher LIR.

In the text

Fig. 9 Selection limits introduced in the Tdust−LIR parameter space by our deepest PACS observations, i.e., in GOODS-S-ultradeep. Red continuous lines are these selection limits at different redshifts using our PACS 5σ detection limits. Each line corresponds to the central redshift of our redshift bins, i.e., z = 0.25 (further left), z = 0.55, z = 0.85, z = 1.15, z = 1.55 and z = 2.05 (further right). The shaded area shows the local Tdust−LIR relation found by Chapman et al. (2003), linearly extrapolated to 1013 L⊙. Dot-dashed lines show the Tdust−LIR relation inferred by Symeonidis et al. (2013) using a sample of high-redshift (i.e., 0.2 < z < 1.2) Herschel-detected galaxies. Dashed black lines show, for each redshift, the lowest infrared luminosities probed by our ultradeep PACS observations without any dust temperature biases and yet populated by star-forming galaxies.

In the text

Fig. 10 Infrared luminosity functions estimated in six redshift bins with the 1/Vmax method. Red squares and black triangles show results from our ultradeep (i.e., GOODS-S-ultradeep) and deep (i.e., GOODS-S-deep and GOODS-N) samples, respectively. Dashed lines represent the best fits to these data points with a double power-law function with fixed slopes. The shaded areas span all the solutions of these fits which are compatible, within 1σ, with our data points: the dark shaded parts of these areas highlight the luminosity ranges directly constrained by our PACS observations, while the light shaded parts highlight the luminosity ranges where our constraints rely upon extrapolations based on φknee and Lknee and a faint-end slope fixed at its z ~ 0 value. Asterisks show the local reference, taken from Sanders et al. (2003), and the dotted line is the best fitted to these data points with our double power-law function with fixed slopes. Red dot-dashed lines are results from Magnelli et al. (2009, 2011) using deep MIPS-24 μm observations. Magnelli et al. (2009, 2011) used the same double power-law function that is used here. To illustrate the infrared luminosity range constrained by Magnelli et al. (2009, 2011), we show, as open circles, their lowest infrared luminosity bins. Blue triple-dot-dashed lines present results from Gruppioni et al. (2013) who used a different analytical function to fit their data points, in particular, a much shallower faint-end slope.

In the text

Fig. 11 Comparisons between the infrared luminosity functions derived using PACS observations and those derived using MIPS observations. PACS-based IR LFs are shown as dark/light shaded areas and dashed lines (see also Fig. 10). MIPS-based IR LFs are derived using MS-based 24 μm-to-LIR conversion factors (black triangles) and above-MS-based 24 μm-to-LIR conversion factors (red squares; see text for more details). The dotted lines show the local reference taken from Sanders et al. (2003). Blue triple-dot-dashed lines present results from Gruppioni et al. (2013). At low luminosities, the agreement between PACS-based and MIPS-based IR LFs confirms that PACS extrapolations towards lower infrared luminosities (i.e., light shaded areas) are reliable, at least down to the infrared luminosities probed by Spitzer.

In the text

Fig. 12 (Top left panel) Evolution of φknee and Lknee as a function of redshift. (Top right panel) Evolution of the comoving number density of “faint” galaxies (i.e., 107 L⊙ < LIR < 1011 L⊙; black triangles), LIRGs (1011 L⊙ < LIR < 1012 L⊙; orange diamonds) and ULIRGs (LIR > 1012 L⊙; red stars). Green circles show the number density of galaxies which are above the PACS detection limit of our ultradeep sub-sample (i.e., >${\mathit{L}}_{\mathrm{IR}}^{\min }$, defined in Fig. 9). The local reference is taken from Sanders et al. (2003). (Bottom left panel) Evolution of the total comoving IR energy density (black area) and the relative contribution of the “faint” galaxies (light blue area), LIRGs (orange area) and ULIRGs (red area). The black dot-dashed line shows the best fit of the total comoving IR energy density, i.e., defined by integrating the IR LFs best-fitting our data points (see dashed lines in Fig. 10). The green dashed line shows the best fit of the comoving IR energy density of galaxies which are above the PACS detection limit of our ultradeep sub-sample (i.e., >${\mathit{L}}_{\mathrm{IR}}^{\min }$). Black circles represent the total comoving IR energy density inferred by Gruppioni et al. (2013). The right axis of the panel shows the evolution of the obscured SFR density assuming that SFR   [M⊙  yr-1]  = 10-10 × LIR  [L⊙] for a Chabrier (2003) IMF. (Bottom right panel) Evolution of the comoving IR energy density. Areas are the same as in the bottom left panel. The dotted lines represents the unobscured SFR density of the Universe (i.e., not corrected for extinction; Cucciati et al. 2012). The red dashed line on a white background shows the total SFR density of the Universe, defined as the sum of the obscured and unobscured SFR densities. Yellow diamonds are a compilation of SFR density estimates from Hopkins & Beacom (2006). Green circles and squares are the total SFR density of the Universe and the contribution of ULIRGs estimated by Murphy et al. (2011) using deep MIPS-24 and -70 μm observations. Dark blue diamonds, stars and triangles are estimates from Casey et al. (2012) for all galaxies, LIRGs and ULIRGs, respectively, using SPIRE (i.e., submillimetre) observations. One has to keep in mind that even with our ultradeep PACS data we only constrain the IR LFs of galaxies down to LIR = 1011 L⊙ at z ~ 1 and LIR = 1012 L⊙ at z ~ 2, respectively. Therefore, number densities and IR LDs of galaxies below these limits rely upon extrapolations based on φknee and Lknee and a faint-end slope fixed at its z ~ 0 value. The redshift and luminosity ranges for which the inferred number densities and IR LDs rely upon these extrapolations are highlighted by dashed lines and striped regions, respectively. We note that these extrapolations are nevertheless corroborated by direct constraints based on deep MIPS-24 μm observations (see Fig. 11).

In the text

	Fig. A.1 Comparison between blind and prior source catalogues. (Left panel) Direct comparison of flux densities for sources in common to both catalogues. (Centre panel) Distribution of unmatched sources in absolute number. (Right panel) Distribution of unmatched sources in fraction relative to total in the given flux density bin.
In the text