B Magnelli
email: magnelli@astro.uni-bonn.de

P Popesso

S Berta

F Pozzi

D Elbaz

D Lutz

M Dickinson

B Altieri

P Andreani

H Aussel

M Béthermin

A Bongiovanni

J Cepa

V Charmandaris

R.-R Chary

A Cimatti

E Daddi

N M Förster Schreiber

R Genzel

C Gruppioni

M Harwit

H S Hwang

R J Ivison

G Magdis

R Maiolino

E Murphy

R Nordon

M Pannella

A Pérez García

A Poglitsch

D Rosario

M Sanchez-Portal

P Santini

D Scott

E Sturm

L J T Acconi

I Valtchanov

The deepest Herschel-PACS far

Keywords: galaxies, evolution -infrared, galaxies -galaxies, starburst -galaxies, statistics

published or not. The documents may come

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. 2010Berta et al. , 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. 2010Gruppioni et al. , 2013;;Casey et al. 2012) and the examination of their spectral energy distributions (e.g., Hwang et al. 2010;Elbaz et al. 2010Elbaz et al. , 2011;;Nordon et al. 2010Nordon et al. , 2012;;Magnelli et al. 2010Magnelli et al. , 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. 2010Gruppioni et al. , 2013)), and their physical properties differ significantly from those of their local counterparts (e.g., Elbaz et al. 2011;[START_REF] Wuyts | [END_REF]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 (PEP1 ; Lutz et al. 2011) guaranteed time key programme; and the GOODS-Herschel (GOODS-H2 ; 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 catalogues3 . 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 starformation 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 0 = 71 km s -1 Mpc -1 , Ω Λ = 0.73 and Ω M = 0.27.

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 .U s i n g 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 µm4 . 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). Notes. (a) Prior to be combined, observations of the ECDFS were trimmed to match the GOODS-S layout of the PEP observations. (b) These observations were taken during the Herschel performance verification (PV) phase with preliminary instrument settings and thus were appropriately underweighted for the map creation (Sect. 3).

The entire 11 ′ × 17 ′ GOODS-S field was observed by PEP with PACS at 70, 100 and 160 µm5 . 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).

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,a10 ′ × 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)andElbaz et al. (2011), respectively.

Map creation

Observations were reduced using the standard PACS photometer pipeline [START_REF] Wieprecht | Astronomical Data Analysis Software and Systems XVIII[END_REF]) 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)andPopesso 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 µma s t r o m etry, and therefore our PACS maps, match the GOODS ACS version 2 coordinate system6 . 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 24 > 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 highpass 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).

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.

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.

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 µmm a p s ( Magnelli et al. 2011) 7 and then use the 24 µm-detected sources (i.e., with S 24 > 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 .B a s e do nt h er e l ative 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 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. 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.4andz ∼ 1.3where the 18 µma n d9 .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 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 L IR = 10 11.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. 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.

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.

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.

Monte Carlo simulations

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(Berta et al. , 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 24 = 20 µJy) sources, i.e., S min λ = min(S 24 ) × mean(S λ /S 24 ), where mean(S λ /S 24 ) 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.6mJy 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 A132, page 7 of 22 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.

Content of the PEP/GOODS-H released package

The PEP/GOODS-H released package8 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. Notes. (a) Numbers listed in this table are only suitable for sources situated in the regions of the field with exposure time higher than that reported in this column.To convert the PACS coverage maps provided in the released package into hrs/pixel, users should multiply them by 1.06 × 10 -2 .

(b) The 1σ flux density levels have been computed from 10 000 random extractions on residual maps. (c) Number of sources above a given S/N threshold. (d) Completeness is defined as the fraction of simulated sources with S in ∼ 3(5)σ and extracted with a flux accuracy better than 50%, i.e., -0.5 < (S out -S in )/S out < 0.5. (e) Contamination is defined as the fraction of simulated sources introduced with S in < 2σ but extracted with S out > 3σ. ( f ) Flux densities above which our catalogues are 80% complete.

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)andHwang 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

clean_index = Neib24 + Neib100 × 10 + Neib160 × 100, (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 ′′ ;s e eElbaz 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 release9 (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.

Confusion noise

Following the formalism of [START_REF] Dole | [END_REF], 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)a n d 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.1h/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 = 5w a s 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 σ T = σ 2 I + σ 2 c . 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 (σ 100 I = 0.11 mJy) and it is fully dominated by the photometric confusion noise at 160 µm( σ 160 I = 0.21 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 µmand2.0 mJy at 160 µmha ve to be treated with caution (see Sect. 4.4).

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 µm, respectively). 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 2 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 [START_REF] Dole | [END_REF]. The density limit of 16.7 beams/source of [START_REF] Dole | [END_REF] 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 Fig. 6. PACS 100 and 160 µmdifferential 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). 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.

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 µmd i fferential number counts are derived in Berta et al. 2011). For that purpose, we applied the method described in Berta et al. (2010Berta et al. ( , 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 subsample containing sources in the GOODS-S-ultradeep field; and (ii) a deep sub-sample containing sources in the GOODS-Sdeep and GOODS-N fields (these two fields having the same depth). The ultradeep sub-sample covers an effective area of 47 arcmin 2 , while the deep sub-sample covers an effective area of 327 arcmin 2 .

The PACS 100 and 160 µmd i fferential 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.5dex 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 . , 1 2 .61 +8.31 -1.74 and 13.63 +3.53 -0.85 nW m -2 sr -1 at 100 and 160 µm, respectively, by power-law fitting of PEP data), we find that in the GOODS-S-ultradeep field our PACS observations resolve ∼75 +12

-30 %and∼75 +7 -15 % 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 100 [mJy] < 1.5) have a median redshift of z = 1.37 +0.58 -0.52 (errors give the interquartile range), while brighter sources (i.e., S 100 [mJy] > 1.5) have a median redshift of z = 0.85 +0.41 -0.33 . Similarly, faint PACS 160 µm sources (i.e., S 160 [mJy] < 2.5) have a median redshift of z = 1.22 0.68 -0.41 while brighter sources (i.e., S 160 [mJy] > 2.5) have a median redshift of z = 0.94 +0.52 -0.38 . At 100 µm (160 µm), sources at 0 < z < 0.5, 0.5 < z < 1.0, 1.0 < z < 2.0andz > 2.0 contribute 24 (17)%, 36 (33)%, 26 (28)% and 14 (22)% of the CIB resolved by our PACS observations. 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 A132, page 11 of 22 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)a n dN i e m ie ta l .(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.

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 midinfrared 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 2 ), 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 multiwavelength 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 10 (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 2 and 184 arcmin 2 , 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 farinfrared 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 L PACS+MIPS IR /L PACS IR 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)a nd because the MIPS-24 µm flux density might be affected by emission from an AGN, we decided not to use the MIPS-24 µmflux densities to derive L IR . We note that using the SED template library of Chary & Elbaz (2001) instead of that of Dale & Helou (2002)t ode r i v eL IR (again leaving the normalisation as a free parameter), has no impact on our results. Indeed, the L DH IR /L CE01 IR 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. The blue line on a white background shows the L IR 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 L IR .

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,P A C Sa n d SPIRE observations) of a large sample of star-forming galaxies at 0 < z < 2.

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 z max -V z min 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

φ = 1 ∆L 1 V max,i × w i , (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., w i = max(w 70 i ,w 100 i ,w 160 i ). 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 T dust -L IR 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 T dust -L IR relation found by Chapman et al. (2003), linearly extrapolated to 10 13 L ⊙ . Dot-dashed lines show the T dust -L IR 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.

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; z max > z lower bin + (z upper binz lower bin )/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)t ode r i v et heL 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., z max > z lower bin +(z upper bin -z lower bin )/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, Fig. 10. Infrared luminosity functions estimated in six redshift bins with the 1/V max 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 L knee 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. (2009Magnelli et al. ( , 2011) ) using deep MIPS-24 µm observations. Magnelli et al. (2009Magnelli et al. ( , 2011) ) used the same double power-law function that is used here. To illustrate the infrared luminosity range constrained by Magnelli et al. (2009Magnelli et al. ( , 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.

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 11 ; 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 ± 4and 34 ± 5Katz = 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 11 To be consistent with Chapman et al. (2003), the T dust -L IR relation of Symeonidis et al. (2013) is obtained from their f 60/ f 100-L IR relation converted using the f 60/ f 100 and assigned dust temperature of each Dale & Helou template. infrared luminosity errors. For the plot, they were defined as the quadratic sum of the Poissonian errors (∝1/ √ 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. (2009Magnelli et al. ( , 2011)).

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. (2009Magnelli et al. ( , 2011) ) using deep MIPS-24 µm observations and Herschel-based estimates from Gruppioni et al. (2013)thatuse shallower PACS observations covering a wider effective area. At z < 1.8, we find good agreement with the Spitzer analysis of Magnelli et al. (2009Magnelli et al. ( , 2011)). In particular, we observe 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-L IR conversion factors (black triangles) and above-MS-based 24 µm-to-L IR 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.

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. (2009Magnelli et al. ( , 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.T h e most significant difference between the present results and those of Magnelli et al. (2009Magnelli et al. ( , 2011) is observed at z ∼ 2a n da t the highest infrared luminosities (L IR > 10 12 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. (2009Magnelli et al. ( , 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.0keV] > 3 × 10 42 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. (2009Magnelli et al. ( , 2011))areaffected by large pre-Herschel uncertainties on the 24 µm-to-L IR conversion factors (Nordon et al. 2010;Elbaz et al. 2010Elbaz et al. , 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 funct i o nu s e dt ofi tt h ed a t a . Gruppioni et al. (2013)u s e daf a i n t 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.0andφ 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.0t h a ni nMagnelli et al. (2011).

With our ultradeep PACS data we can only constrain the IR LFs of galaxies down to L IR = 10 11 L ⊙ at z ∼ 1a n d L IR = 1012 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 µmluminosity (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 andTable A.5) of "faint" galaxies (i.e., 10 7 L ⊙ < L IR < 10 11 L ⊙ ), LIRGs (i.e., 10 11 L ⊙ < L IR < 10 12 L ⊙ )a n d ULIRGs (i.e., L IR > 10 12 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 11 L ⊙ at z ∼ 1a n d L IR = 10 12 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 ∼ 0v a l u e .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;[START_REF] Wuyts | [END_REF]Nordon et al. 2012;Magnelli et al. 2012).

We find that the number density of ULIRGs strongly evolves from z = 0t oz ∼ 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 = 0a n dz ∼ 1, while relying on extrapolations we find that it is multiplied by a factor ∼30 between z = 0a n dz ∼ 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 = 0a n dz ∼ 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 = 0andz ∼ 2, with most of the evolution happening before z ∼ 1.2d r i v e nb y 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 a res omedi fferences. 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)andSpitzer (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)atz > 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 SPIREbased 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 12 : Green circles show the number density of galaxies which are above the PACS detection limit of our ultradeep sub-sample (i.e., >L min IR , 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 dotdashed 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., >L min IR ). 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 × L IR [L ⊙ ]foraChabrier (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 L IR = 10 11 L ⊙ at z ∼ 1a n dL IR = 10 12 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. We note that these extrapolations are nevertheless corroboratedby direct constraints based on deep MIPS-24 µm observations (see Fig. 11).

SFR [M ⊙ yr -1 ] = 10 -10 × L IR [L ⊙ ], (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 = 0toz ∼ 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 ∼ 1a n dz ∼ 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)andKarim et al. (2011, not shown in Fig. 12).

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 11 L ⊙ at z ∼ 1a n dL IR = 10 12 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 = 0t oz ∼ 1a n dt h e n 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 . ,5 0± 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.

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 crossidentification 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 µmand 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. [log (Mpc -3 dex -1 )] z ∼ 0 -0.60 -2.20 10.48 ± 0.02 -2.52 ± 0.03 0.1 < z < 0.4 -0.60 -2.20 10.84 ± 0.06 -2.85 ± 0.04 0.4 < z < 0.7 -0.60 -2.20 11.28 ± 0.12 -2.82 ± 0.14 0.7 < z < 1.0 -0.60 -2.20 11.53 ± 0.15 -2.96 ± 0.18 1.0 < z < 1.3 -0.60 -2.20 11.71 ± 0.14 -3.01 ± 0.20 1.3 < z < 1.8 -0.60 -2.20 12.00 ± 0.15 -3.29 ± 0.19 1.8 < z < 2.3 -0.60 -2.20 12.35 ± 0.19 -3.47 ± 0.23

Notes. (a) Fixed slopes of the infrared LF. Notes. The local reference is taken from Sanders et al. (2003). Comoving infrared energy densities constrained using our PACS observations are highlighted in bold. Others values rely upon extrapolations based on φ knee and L knee and a faint-end slope of the IR LF fixed at its z ∼ 0value.

× L ⊙ Mpc -3 ][ 1 0 6 × L ⊙ Mpc -3 ][ 1 0 6 × L ⊙ Mpc -3 ][ 1 0 6 × L ⊙ Mpc -3 ] 0.

Fig. 1 .

 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.

Fig. 2 .

 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.

Fig. 4 .

 4 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 (S out -S in )/S out distribution in each flux bin (after 3σ clipping). Red lines show the mean value of the (S out -S in )/S out 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).

Fig. 5 .

 5 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 inBerta et al. (2011) and illustrate the smaller exposure time range probed in their study.

Fig. 7 .

 7 Fig. 7. (Top panel) Uncertainties in determining L IR from monochromatic observations (i.e., MIPS-24 µm, PACS-70 µm, PACS-100 µmor PACS-160 µm bands) and the Dale & Helou (2002) SED library. These uncertainties are derived by taking the standard deviation of the log(L IR ) distribution provided when normalizing all Dale & Helou (2002)S E D 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 µmor70+ 100 + 160 µm detections) as a function of redshift.

Fig. 8 .

 8 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 L IR 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 L IR .

Fig. 12 .

 12 Fig. 12. (Top left panel) Evolution of φ knee and L knee as a function of redshift. (Top right panel)Evolution of the comoving number density of "faint" galaxies (i.e., 10 7 L ⊙ < L IR < 10 11 L ⊙ ; black triangles), LIRGs (10 11 L ⊙ < L IR < 10 12 L ⊙ ; orange diamonds) and ULIRGs (L IR > 10 12 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., >L min IR , defined in Fig.9). The local reference is taken fromSanders 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 dotdashed 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., >L min IR ). Black circles represent the total comoving IR energy density inferred byGruppioni 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 × L IR [L ⊙ ]foraChabrier (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 fromHopkins & Beacom (2006). Green circles and squares are the total SFR density of the Universe and the contribution of ULIRGs estimated byMurphy et al. (2011) using deep MIPS-24 and -70 µm observations. Dark blue diamonds, stars and triangles are estimates fromCasey 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 L IR = 10 11 L ⊙ at z ∼ 1a n dL IR = 10 12 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. We note that these extrapolations are nevertheless corroboratedby direct constraints based on deep MIPS-24 µm observations (see Fig.11).

Fig. A. 1 .

 1 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.

Table 1 .

 1 Properties of the Herschel observations combined for the PEP/GOODS-H data release.

	Field	Survey	RA	Dec	Wavelengths	Size	Time
			[Degree, J2000]	[µm]	[arcmin]	[h]
	GOODS-N GOODS-N GOODS-S ECDFS a (GOODS-S) GOODS-S GOODS-S b	PEP GOODS-H 189.22862 189.22862 PEP 53.12654 PEP 53.10417 GOODS-H 53.12654 PV 53.12654	62.23867 62.23867 -27.80467 70, 100, 160 11 ′ × 17 ′ 100, 160 11 ′ × 17 ′ 100, 160 11 ′ × 17 ′ 124.0 25.8 226 -27.81389 100, 160 32.8 30 ′ × 30 ′ -27.80467 100, 160 10 ′ × 10 ′ 206.3 -27.80467 100, 160 11 ′ × 17 ′ 7.9

Table 2 .

 2 Statistics of released catalogues, for regions above the specified exposure time.

	Field	t exposure	a > S (1σ) b	Num. c	Num. c	Compl. d	Contam. e	Compl. d	Contam. e	S (80%) f
	and bands	hrs/pixel	mJy	S /N > 3 S /N > 5% a t 3 σ	%at3σ	%at5σ	%at5σ	mJy
					Prior source extraction
	GOODS-N 100	1.30 .32	914	574	36	14	72	5	1.97
	GOODS-N 160	1.30 .70	838	546	24	32	62	5	4.56
	GOODS-S 70	1.00 .30	476	245	21	24	77	2	1.39
	GOODS-S 100 ultradeep	6.40 .18	472	346	21	27	65	5	1.22
	GOODS-S 160 ultradeep	8.30 .43	445	313	18	38	45	14	3.63
					Blind source extraction
	GOODS-N 100	1.30 .32	865	596	13	37	56	7	2.04
	GOODS-N 160	1.30 .70	867	521	16	48	53	17	5.88
	GOODS-S 70	1.00 .33	396	205	19	19	72	11	1.51
	GOODS-S 100 ultradeep	6.40 .18	513	377	21	42	74	9	1.02
	GOODS-S 160 ultradeep	8.30 .43	485	368	15	47	45	8	4.89

Table A .

 A 1. PACS 100 and 160 µm number counts, normalized to the Euclidean slope.

			PACS 100 µm					PACS 160 µm
		GOODS-S ultradeep	GOODS-N/S deep		GOODS-S ultradeep	GOODS-N/S deep
	S centre Counts	Err.	Counts	Err.	S centre Counts	Err.	Counts	Err.
	[mJy]	[10 4 deg -2 mJy 1.5 ][ 1 0 4 deg -2 mJy 1.5 ][ m J y ]	[ 1 0 4 deg -2 mJy 1.5 ][ 1 0 4 deg -2 mJy 1.5 ]
	0.56 0.71 0.89 1.13	0.59 0.91 1.34 1.78	0.10 0.14 0.20 0.28	------1.29 0.10	1.42 1.79 2.25 2.84	3.59 4.91 6.49 7.56	0.49 0.72 0.91 1.15	----4.92 0.31 6.48 0.41
	1.42	1.95	0.36	1.90	0.13	3.57	8.62	1.54	8.80	0.56

Table A .

 A 4. Parameter values of the infrared LF.

	Redshift	α 1	a	α 2

a log (L knee )l o g ( φ knee ) [log (L ⊙ )]

Table A .5.

 A Evolution of the comoving IR energy density and the relative contribution of "faint" galaxies (10 7 L ⊙ < L IR < 10 11 L ⊙ ), LIRGs (10 11 L ⊙ < L IR < 10 12 L ⊙ ) and ULIRGs (L IR > 10 12 L ⊙ ).

	Redshift	Total IR LD	"faint" galaxies	LIRGs	ULIRGs
		[10 6

http://www.mpe.mpg.de/ir/Research/PEP

http://hedam.oamp.fr/GOODS-Herschel

http://www.mpe.mpg.de/ir/Research/PEP/public_data_ releases.php

GOODS-H has also covered the entire GOODS-N field with SPIRE at 250, 350 and

µm. SPIRE observations and catalogues are described inElbaz et al. (2011) and are publicly available at http:// hedam.oamp.fr/GOODS-Herschel

SPIRE-250/350/500 µm observations of the GOODS-S field have been performed by the HerMES survey(Oliver et al. 2012). These observations are publicly available through the Herschel Database in Marseille (HeDaM) at http://hedam.oamp.fr/HerMES

See the GOODS ACS data release at MAST: http://archive. stsci.edu/prepds/goods/

The MIPS-24 µm catalogues used here are slightly different from those released byMagnelli et al. (2011). See Sect. 4.5 for more details.

http://www.mpe.mpg.de/ir/Research/PEP/ public_data_releases.php

http://hedam.oamp.fr/GOODS-Herschel

Publicly available at http://www.mpe.mpg.de/ir/Research/ PEP/

Kennicutt (1998) adopts a Salpeter initial mass function so that we divide his normalization by 1.72.

Pages 19 to 22 are available in the electronic edition of the journal at http://www.aanda.org

Acknowledgements. 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

⋆ Based on observations carried out by the Herschel Space Observatory. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. at http://www.aanda.org