# The HerMES submillimetre local and low-redshift luminosity functions

L. Marchetti, M. Vaccari, A. Franceschini, V. Arumugam, H. Aussel, M. Béthermin, J. Bock, A. Boselli, V. Buat, D. Burgarella, et al.

► To cite this version:

L. Marchetti, M. Vaccari, A. Franceschini, V. Arumugam, H. Aussel, et al. The HerMES submillimetre local and low-redshift luminosity functions. Monthly Notices of the Royal Astronomical Society, 2016, 456, pp.1999-2023. 10.1093/mnras/stv2717. hal-01433389

HAL Id: hal-01433389

https://hal.science/hal-01433389

Submitted on 27 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The HerMES submillimetre local and low-redshift luminosity functions*

L. Marchetti,<sup>1,2†</sup> M. Vaccari,<sup>2,3,4</sup> A. Franceschini,<sup>2</sup> V. Arumugam,<sup>5</sup> H. Aussel,<sup>6</sup>

M. Béthermin,<sup>6,7,8</sup> J. Bock,<sup>9,10</sup> A. Boselli,<sup>11</sup> V. Buat,<sup>11</sup> D. Burgarella,<sup>11</sup>

D.L. Clements,<sup>12</sup> A. Conley,<sup>13</sup> L. Conversi,<sup>14</sup> A. Cooray,<sup>9,15</sup> C. D. Dowell,<sup>9,10</sup>

D. Farrah,<sup>16</sup> A. Feltre,<sup>17</sup> J. Glenn,<sup>13,18</sup> M. Griffin,<sup>19</sup> E. Hatziminaoglou,<sup>8</sup> S. Heinis,<sup>11</sup>

E. Ibar,<sup>20</sup> R. J. Ivison,<sup>5,21</sup> H. T. Nguyen,<sup>9,10</sup> B. O’Halloran,<sup>12</sup> S. J. Oliver,<sup>22</sup>

M. J. Page,<sup>23</sup> A. Papageorgiou,<sup>19</sup> C. P. Pearson,<sup>1,24</sup> I. Pérez-Fournon,<sup>25,26</sup>

M. Pohlen,<sup>19</sup> D. Rigopoulou,<sup>24,27</sup> I. G. Roseboom,<sup>5,22</sup> M. Rowan-Robinson,<sup>12</sup>

B. Schulz,<sup>9,28</sup> Douglas Scott,<sup>29</sup> N. Seymour,<sup>30</sup> D. L. Shupe,<sup>9,28</sup> A. J. Smith,<sup>22</sup>

M. Symeonidis,<sup>23</sup> I. Valtchanov,<sup>14</sup> M. Viero,<sup>9</sup> L. Wang,<sup>31,32</sup> J. Wardlow,<sup>33</sup>

C. K. Xu<sup>9,28</sup> and M. Zemcov<sup>9,10</sup>

*Affiliations are listed at the end of the paper*

Accepted 2015 November 17. Received 2015 November 16; in original form 2014 December 1

## ABSTRACT

We used wide-area surveys over 39 deg² by the HerMES (Herschel Multi-tiered Extragalactic Survey) collaboration, performed with the Herschel Observatory SPIRE multiwavelength camera, to estimate the low-redshift, 0.02 < z < 0.5, monochromatic luminosity functions (LFs) of galaxies at 250, 350 and 500 µm. Within this redshift interval, we detected 7087 sources in five independent sky areas, ~40 per cent of which have spectroscopic redshifts, while for the remaining objects photometric redshifts were used. The SPIRE LFs in different fields did not show any field-to-field variations beyond the small differences to be expected from cosmic variance. SPIRE flux densities were also combined with Spitzer photometry and multiwavelength archival data to perform a complete spectral energy distribution fitting analysis of SPIRE detected sources to calculate precise k-corrections, as well as the bolometric infrared (IR; 8–1000 µm) LFs and their low-z evolution from a combination of statistical estimators. Integration of the latter prompted us to also compute the local luminosity density and the comoving star formation rate density (SFRD) for our sources, and to compare them with theoretical predictions of galaxy formation models. The LFs show significant and rapid luminosity evolution already at low redshifts, 0.02 < z < 0.2, with $L_{IR}^* \propto (1+z)^{6.0 \pm 0.4}$ and $\Phi_{IR}^* \propto (1+z)^{-2.1 \pm 0.4}$, $L_{250}^* \propto (1+z)^{5.3 \pm 0.2}$ and $\Phi_{250}^* \propto (1+z)^{-0.6 \pm 0.4}$ estimated using the IR bolometric and the 250 µm LFs, respectively. Converting our IR LD estimate into an SFRD assuming a standard Salpeter initial mass function and including the unobscured contribution based on the UV dust-uncorrected emission from local galaxies, we estimate an SFRD scaling of SFRD₀ + 0.08z, where SFRD₀ ≃ (1.9 ± 0.03) × 10⁻² [M☉ Mpc⁻³] is our total SFRD estimate at z ~ 0.02.

**Key words:** galaxies: evolution – galaxies: luminosity function, mass function – galaxies: statistics – submillimetre: galaxies.

## 1 INTRODUCTION

Observations carried out in roughly the past 20 years have revealed a rapid evolution of cosmic sources, both normal, actively star-forming and AGN-dominated galaxies, over the last several billion years. This was mostly achieved from continuum rest-frame UV

<footnote>* Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.</footnote>

<footnote>† E-mail: marchetti.lu@gmail.com</footnote>
photometric imaging in the optical (e.g. Lilly et al. 1995) and Hα or [O II] line spectroscopy (e.g. Gallego et al. 1995), all, however, including very uncertain dust extinction corrections. GALEX has also been used for similar purposes by Martin et al. (2005) and Bothwell et al. (2011), among others.

In the far-IR (FIR), the pioneering exploration by the IRAS satellite revealed a particularly dramatic evolution of the galaxy luminosity functions (LFs; Saunders et al. 1990), illustrating the importance of local studies at infrared (IR) wavelengths. This result was later confirmed up to $z \simeq 1$ by ISO (Pozzi et al. 2004) and Spitzer studies using the Multi-Band Imaging Photometer (MIPS) 24 µm (Le Floc'h et al. 2005; Marleau et al. 2007; Rodighiero et al. 2010) and 70 µm (Frayer et al. 2006; Huynh et al. 2007; Magnelli et al. 2009; Patel et al. 2013) channels. At longer, submillimetre wavelengths, the Balloon-borne Large Aperture Submillimeter Telescope (BLAST) was able to estimate the galaxy LF at low $z$ and map its evolution (Eales et al. 2009), although with limited statistics and uncertain identification of the sources. Finally, surveys in the radio bands have also been exploited, with the necessity to include large bolometric corrections, for LF estimates (Condon 1989; Serjeant, Gruppioni & Oliver 2002).

Interpretations of these fast evolutionary rates are actively debated in the literature, with various processes being claimed as responsible (like gas fuel consumption in galaxies, heating of the gas so as to prevent cooling and collapse, decreasing rates of galaxy interactions with time, etc.). Indeed, galaxy evolution codes have often found it difficult to reproduce these data, and slower evolution seems predicted by the models than it is observed.

However, the estimates of the low-redshift LFs of galaxies, and correspondingly the total star formation and AGN accretion rates, still contain some significant uncertainties. In particular, due to the moderate volumes sampled at low redshift, an essential prerequisite for determining the local luminosity functions (LLFs) is the imaging of large fields, where it is difficult however to achieve the required multiwavelength homogeneous coverage and complete redshift information.

In the very local Universe, at $z < 0.02$, a sample of a few hundred sources from the Early Release Compact Source Catalogue by the Planck all-sky survey (Planck Collaboration VII 2011) have been used by Negrello et al. (2013) to estimate LFs at 350, 500 and 850 µm. Although the authors were very careful to account for various potentially problematic factors, namely the photometric calibration from the large Planck beam, removal of Galactic emission and CO line contribution to the photometry, their estimate might not be completely immune to the effects of large inhomogeneities (like the Virgo cluster) inherent in their very local spatial sampling (see Section 5 for further details).

Vaccari et al. (2010) report a preliminary determination of the local submillimetre LFs of galaxies, exploiting the much improved angular resolution and mapping speed of the SPIRE instrument (Griffin et al. 2010) on the Herschel Space Observatory (Pilbratt et al. 2010). They used data from the Lockman Hole (LH) and Extragalactic First Look Survey (XFLS) fields of the Herschel Multi-tiered Extragalactic Survey (HerMES) programme (http://hermes.sussex.ac.uk; Oliver et al. 2012) over about 15 deg² observed during the Herschel Science Demonstration Phase (SDP), and including a few hundred sources to a flux limit of about 40 mJy in all three SPIRE bands (250, 350, 500 µm). Their published functions were integrated over a wide redshift interval at $0 < z < 0.2$. Because of the limited source statistics, Vaccari et al. (2010) could not take into account any evolutionary corrections, while significant evolution is expected to be present over this large redshift bin.

Still based on the HerMES data base, but using a much larger total area, many more independent sky fields and deeper fluxes, this paper reports on a systematic effort to characterize the local and low-redshift LFs of galaxies in the submillimetre bins. The Herschel survey catalogue has been cross-correlated with existing optical photometry and spectroscopy in the fields, as well as with photometric data in the mid-IR (MIR) and FIR from Spitzer (Werner et al. 2004). By fitting the source-by-source multiwavelength photometry with spectral templates, the bolometric IR luminosities and bolometric LFs can also be estimated. Importantly, the much improved statistics allows us to work in narrow redshift bins, so as to disentangle LF shapes from evolution, and to obtain the most robust and complete statistical characterization over the last few Gyr of galaxy evolution. By combining this long-wavelength information with similar analyses in the optical–UV, we can determine the local bolometric luminosity density and comoving star formation rate (SFR) and their low-$z$ evolution.

The paper is structured as follows. In Section 2, we describe the multiwavelength data set that we use, as well as the selection of the samples, source identification and spectral energy distribution (SED) fitting. In Section 3, we detail the statistical methods used in our data analysis and the various adopted LF estimators, including the Bayesian parametric recipe that we develop. Our results are reported in Section 4, including the multiwavelength LFs, the local luminosity densities (LLDs) and the SFRs. Our results are then discussed in Section 5 and our main conclusions summarized in Section 6.

Throughout the paper, we adopt a standard cosmology with $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$ and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

## 2 THE HERMES WIDE SAMPLE

HerMES is a Herschel Guarantee Time Key Programme (Oliver et al. 2012¹) and the largest single project on Herschel, for a total 900 h of observing time. HerMES was designed to comprise a number of tiers of different depths and areas, and has performed observations with both SPIRE (Griffin et al. 2010) and PACS (Poglitsch et al. 2010), surveying approximately 70 deg² over 12 fields whose sizes range from 0.01 to 20 deg².

To estimate the SPIRE LLF, we use HerMES L5/L6 SPIRE observations (see table 1 in Oliver et al. 2012 for more details on the observations) covering five fields: LH, XFLS, Bootes, ELAIS-N1 (EN1) and XMM-Large Scale Structure (XMM-LSS). In the following, these fields and the SPIRE sample arising from them will collectively be referred to as the HerMES Wide Fields and Sample, respectively. These fields are the widest Herschel HerMES fields where imaging data are available with both Spitzer Infrared Array Camera (IRAC) and MIPS, thus enabling the detailed study of the full IR SED of a significant number of sources in the local Universe.

### 2.1 SPIRE source extraction

Source confusion is the most serious challenge for Herschel and SPIRE source extraction and identification. In particular, confusion is an important driver in determining the optimal survey depth. By making a maximum use of the full spectrum of ancillary data, it is possible to limit the confusion problem at the source detection and identification steps. For this reason, the choice of HerMES survey fields has been largely driven by the availability of extensive multiwavelength ancillary data at both optical and IR wavelengths.

<footnote n="1">¹ http://hermes.sussex.ac.uk/</footnote>
In particular, Roseboom et al. (2010, and Roseboom et al. 2012) developed a new method for SPIRE source extraction, hereafter referred to as XID, which improves upon more conventional techniques (e.g. Smith et al. 2012; Wang et al. 2014) based on existing source extraction algorithms optimized for Herschel purposes.

The XID technique makes use of a combination of linear inversion and model selection techniques to produce reliable cross-identification catalogues based on *Spitzer* MIPS 24 µm source positions. The tiered nature of HerMES is well matched to the variable quality of the *Spitzer* data, in particular the MIPS 24 µm observations. This is confirmed by simulation performed using pre-*Herschel* empirical models (e.g. Fernandez-Conde et al. 2008; Le Borgne et al. 2009; Franceschini et al. 2010) which shared the comparable sensitivities of the 250 and 24 µm source densities. Since the HerMES Wide Fields are homogeneously covered by the *Spitzer* Data Fusion (described in Section 2.3), which provides homogeneous MIPS 24 µm source lists, the SPIRE flux densities used in this paper are obtained with the XID technique using the *Spitzer* Data Fusion MIPS 24 µm positional priors (or, in other words, the MIPS 24 µm positions are used as a prior to guide the SPIRE flux extraction on the SPIRE maps).

As reported in Roseboom et al. (2010), using a prior for the SPIRE source identification based on MIPS 24 µm detections could, in principle, introduce an additional incompleteness related to the relative depth at 24 µm catalogues used in the process and the distribution of intrinsic SED shapes. However, Roseboom et al. (2010) show how incompleteness would affect only the fainter SPIRE sources with the higher 250 µm/24 µm flux density ratios, which are very likely to be ultra-red high-redshift objects. We can therefore be confident that for relatively nearby sources the XID catalogues are complete at the relatively bright flux limits used in this work. This relatively complex procedure of association is reported in the dedicated papers by Roseboom et al. (2010, 2012), to which we refer the reader for further details about this method.

## 2.2 SPIRE sample selection

To define the sample to be used for our LLF determinations, we use SPIRE flux density estimates obtained using the XID method (Roseboom et al. 2010 and Roseboom et al. 2012), applied to SPIRE maps produced by Levenson et al. (2010) and using MIPS 24 µm positional priors based on the *Spitzer* Data Fusion detailed in Sec-

tion 2.3. The SPIRE 250 µm channel is the most sensitive of the three SPIRE bands and thus we select sources on the basis of a SPIRE 250 µm reliability criterion (discussed in Roseboom et al. 2010) defined as $\chi^2_{250} < 5$ and $\text{SNRT}_{250} > 4$, where the first criterion is the $\chi^2$ of the source solution in the neighbourhood of a source (7 pixel radius) and the second is the signal-to-noise ratio (SNR) at a given selection $\lambda$, including confusion, and referred to as 'total' $\text{SNR}_\lambda$ or $\text{SNRT}_\lambda$.

The SPIRE 250 µm catalogues of L5/L6 HerMES observations are highly complete and reliable down to approximately 25/30/35 mJy at 250/350/500 µm, respectively, as shown in Fig. 1 (left). In order to combine the data collected in these different fields, we have to ensure a uniform completeness both in flux and in redshift coverage across fields; thus, due to some minor differences across the fields, we decided to cut our sample at 30 mJy at 250 µm. These minor differences are visible in Fig. 1 (right) where we compare SPIRE 250 µm number counts estimated for the five wide fields and for the COSMOS deep field (the COSMOS sample is from Vaccari et al., in preparation). These discrepancies are consistent with the levels of cosmic variance predicted by theoretical models for fields of this size (Moster et al. 2011), as well as with the slightly different depths of MIPS 24 µm observations available for these fields, which were used to guide HerMES XID source extraction. In any case, differences are on the whole small and have major effects only at low flux densities, well below our selected limit. The greatest discrepancy is shown in XFLS, where the SPIRE 250 µm completeness reflects the slightly brighter flux limit of the XFLS MIPS 24 µm and IRAC catalogues, due to a shorter exposure time in comparison with the other fields.

## 2.3 The *Spitzer* Data Fusion

As previously mentioned, the HerMES fields were chosen so as to have the best multiwavelength data for sky areas of a given size. In particular, the fields used in this work are covered by *Spitzer* seven-band IRAC and MIPS imaging data which enable not only an improved identification process but also the detailed characterization of the IR SEDs of *Herschel* sources.

In this work, we exploit the *Spitzer Data Fusion* (Vaccari et al. 2010 and Vaccari et al., in preparation, http://www.mattiavaccari.net/df). The *Spitzer* Data Fusion combines *Spitzer* MIR and FIR

<figure><img src="image_3.png" /><figcaption>**Figure 1.** SPIRE 250 µm source counts (right) and completeness (left) based on XID catalogues from Roseboom et al. (2012) for the HerMES Wide Fields sample used to estimate the SPIRE LLFs, compared with COSMOS estimates from Vaccari et al. (in preparation). The black solid line signs the flux limit of our selection.</figcaption></figure>
data from the *Spitzer* Wide-area InfraRed Extragalactic (SWIRE; Lonsdale et al. 2003) survey in six fields, the *Spitzer* Deep-Wide Field Survey (SDWFS, PI: Daniel Stern, Spitzer PID 40839), the *Spitzer* XFLS (PI: Tom Soifer, Spitzer PID 26), with photometric data at UV, optical and near-infrared (NIR) wavelengths, as well as optical spectroscopy over about 70 deg<fnref n="2" /> in total. It thus makes full use of public survey data from the Galaxy Evolution Explorer (GALEX), the Sloan Digital Sky Survey (SDSS), the Issac Newton Telescope Wide Field Survey (INT WFS), the 2 Micron All-Sky Survey (2MASS), the UKIRT Infrared Deep Sky Survey (UKIDSS) and the Visible and Infrared Survey Telescope for Astronomy (VISTA) projects, as well as further optical imaging obtained by the SWIRE, SDWFS and XFLS teams. It also provides spectroscopic information variously available from SDSS, NASA/IPAC Extragalactic Database (NED; http://ned.ipac.caltech.edu), recent literature and proprietary follow-up programmes.

The *Spitzer* Data Fusion thus represents an ideal starting point to perform statistical studies of IR galaxy populations, such as detailed SED fitting analyses to estimate photometric redshifts and masses, as well as SFRs; an early version of the data base has already been used to that effect by Rowan-Robinson et al. (2013). It has been used to validate Herschel SDP observations within the HerMES consortium team and to produce current and future public HerMES catalogues.<fnref n="2" /> Since this paper only uses the *Spitzer* Data Fusion to derive SPIRE LLF estimates, we refer the reader to Vaccari et al. (in preparation) for a complete description of the data base and in the following we only summarize its basic properties as they relate to this work.

The *Spitzer* Data Fusion is constructed by combining *Spitzer* IRAC and MIPS source lists, as well as ancillary catalogues, following a positional association procedure. Source extraction of IRAC four-band images and of MIPS 24 µm images is carried out using SExTRACTOR (Bertin & Arnouts 1996), whereas MIPS 70 and 160 µm source extraction is carried out using *APEX* (Makovoz & Marleau 2005). Catalogue selection is determined by a reliable IRAC 3.6 or IRAC 4.5 µm detection. We then associate MIPS 24 µm detections with IRAC detections using a 3 arcsec search radius, while MIPS 70 and 160 µm catalogues are matched against MIPS 24 µm positions using a search radius of 6 and 12 arcsec, respectively. UV, optical and NIR catalogues are then matched against IRAC positions using a 1 arcsec search radius. This multistep approach increases the completeness and reliability of the longer wavelength associations, while better pin-pointing MIPS sources using their IRAC positions.

The HerMES Wide Fields used in this work are part of the *Spitzer* Data Fusion and are all covered both by *Spitzer* seven-band IR imaging and by SDSS five-band optical imaging and optical spectroscopy (Csabai et al. 2007; Abazajian et al. 2009; Carliles et al. 2010; Bolton et al. 2012). They also benefit by a vast quantity of additional homogeneous multiwavelength observations and additional spectroscopic redshifts available from NED, as well as the recent literature, and our own *Spitzer/Herschel* proprietary follow-up programmes. We thus associate a reliable spectroscopic redshift with our sources whenever this is available and otherwise rely on SDSS photometric redshift estimates based on a KD-tree nearest-neighbour search (see Csabai et al. 2007 for more details). In so doing, we follow a commonly adopted photometric reliability criterion for SDSS good photometry, only selecting detections with SDSS *cmodelmag* $r_{AB} < 22.2$, thus avoiding unreliable photomet-

ric redshifts. In Fig. 2, we report SDSS $r_{AB}$ and redshift histograms of the HerMES Wide sample. In order to avoid effects of incompleteness in redshift, we limit our HerMES Wide sample to $z \lesssim 0.5$, below the completeness and reliability limit of SDSS redshift estimates. Moreover, to avoid the possible redshift incompleteness that affects the very bright and nearby galaxies in SDSS data, we cut our sample to the lowest redshift of $z = 0.02$, as suggested by e.g. Montero-Dorta & Prada (2009). As discussed in Roseboom et al. (2010), the SPIRE source extraction works very well for point-like sources, but can underestimate the fluxes of the extended sources; cutting the sample at $z > 0.02$ also avoids this problem since the vast majority of extended sources are located at lower redshifts. The numbers of sources of the HerMES Wide sample are detailed in Table 1.

## 2.4 SED fitting

Thanks to the *Spitzer* Data Fusion, we are able to perform the multi-wavelength SED fitting analysis of our HerMES Wide Fields sample and thus estimate the IR bolometric (8-1000 µm) and monochromatic rest-frame luminosities and relative *k*-corrections. We perform the SED fitting analysis using LE PHARE (Arnouts et al. 1999 and Ilbert et al. 2006). To perform the fit, we use SDSS *ugriz*, 2MASS *JHKs*, IRAC-3.6, IRAC-4.5, IRAC-5.8, IRAC-8.0, MIPS-24, MIPS-70, MIPS-160, SPIRE-250, SPIRE-350 and SPIRE-500 flux densities, which are available over the whole area covered by our sample. As template SEDs we use two different set of empirical templates according to the range of wavelengths we are fitting: in the optical-MIR range (up to 7 µm rest-frame), we use the same templates and extinction laws exploited by the COSMOS team to estimate the COSMOS photometric redshifts as in Ilbert et al. (2009), while to fit the IR/submm range (from 7 µm rest-frame upwards) we use the SWIRE templates of Polletta et al. (2007) and their slightly modified version described in Gruppioni et al. (2010), for a total of 32 and 31 SEDs, respectively; this includes elliptical, spiral, AGN, irregular and starburst spectral types as summarized in Table 2. As an example we report two typical examples of our SED fitting results in Fig. 3. Splitting the overall wavelength coverage into two provides us with a particularly good fit to the FIR bump and a reasonably good fit at all other wavelengths for all sources, with a mean value of the reduced $\chi^2$ of around 0.5. In Fig. 2 (upper panels), we report the $L-z$ distribution for both the $L_{250}$ and the $L_{IR}$ rest-frame luminosities obtained through the SED fitting procedure.

Thanks to this multiwavelength SED fitting, we are able to also investigate the relation between monochromatic rest-frame luminosities at different wavelengths. As an example we report in Fig. 4 a comparison between SPIRE 250 µm and PACS 100 µm monochromatic rest-frame luminosities plotted against the IR bolometric luminosity. Historically, the monochromatic rest-frame luminosity at 60-100 µm has been considered a good indicator of the IR bolometric luminosity, due to a strong correlation between the two [e.g. Patel et al. (2013) used the relation between MIPS 70 µm and $L_{IR}$]. In Fig. 4, we show that we confirm this trend in our SED fitting results while, on the other hand, the SPIRE 250 µm luminosity does not show a strong correlation with the IR bolometric luminosity and thus cannot be used as a reliable indicator of the total IR emission of the galaxy. As also confirmed by other HerMES works that have carefully studied the SED shape of the HerMES sources (e.g. Symeonidis et al. 2013), we find that the SEDs in the FIR regime of our local HerMES sample peak close to the PACS 100 µm band and thus the monochromatic luminosity at this wavelength best traces the total IR bolometric luminosity integrated between 8 and

<footnote n="2">² available at http://hedam.oamp.fr/HerMES/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top: SPIRE 250 µm (expressed as $vL_v$) and IR bolometric luminosity versus redshift. Bottom: SDSS $r_{AB}$ (left) and redshift histograms (right) for the HerMES Wide Fields sample used to estimate the SPIRE LLFs. The $L-z$ plots are colour-coded according to the SED best-fitting class obtained by the SED fitting procedure following the list reported in Section 2.4. The histograms report the relative quantities for the photometric and spectroscopic samples in blue and in green, respectively, with the total sample illustrated in red.</figcaption></figure>

<table><caption>**Table 1.** Number of $0.02 < z \lesssim 0.5$ sources used to estimate the SPIRE LLFs. The number of sources with spectroscopic/photometric redshifts is indicated in parentheses after the total number of sources. The 250 µm sample is cut at $S_{250} > 30$ mJy, according to the SPIRE 250 µm completeness (see the text for details). ‘Set’ refers to table 1 in Oliver et al. (2012) and identifies the HerMES specific observing mode in each field.</caption><thead><tr><th>Field</th><th>250 µm detections</th><th>Area (deg²)</th><th>Set</th></tr></thead><tbody><tr><td>LH</td><td>2336 (942/1394)</td><td>11.29</td><td>34</td></tr><tr><td>XFLS</td><td>801 (427/374)</td><td>4.19</td><td>40</td></tr><tr><td>BOOTES</td><td>1792 (1220/572)</td><td>9.93</td><td>37</td></tr><tr><td>EN1</td><td>693 (246/447)</td><td>3.91</td><td>35</td></tr><tr><td>XMM</td><td>1606 (367/1239)</td><td>9.59</td><td>36</td></tr><tr><td>Total</td><td>7087 (3195/3892)</td><td>38.9</td><td></td></tr></tbody></table>

1000 µm. It is also interesting to notice the very different behaviour of the k-corrections estimated at SPIRE 250 µm and PACS 100 µm (lowest panels of Fig. 4). The differences between these two are remarkable, and this is reflected in the different behaviour of the resulting luminosities.

While a detailed physical analysis of our sample is beyond the scope of this paper, we did exploit our SED fitting analysis and

the IRAC colour–colour criteria by Lacy et al. (2004) and Donley et al. (2012) to search for any potential AGN contamination in our sample. On the whole, the vast majority of our sources show galaxy- or starburst-like best-fitting SEDs with less than 10 per cent of the sample being best fitted by AGN-like SEDs (SED classes between 17 and 25 and between 28 and 31 as reported in Table 2). These numbers do not change significantly even if we fit a single SED template to the whole range of available photometry (from optical to SPIRE bands). Fig. 5 confirms that our objects mostly lie within the starburst-dominated region of the IRAC colour–colour plot, with only a small fraction of the sources (mainly located at $z > 0.25$) sitting in the area usually occupied by AGN-like objects. On the whole, we find that about 20 per cent of our sources sit in the AGN region identified by Lacy et al. (2004), with less than 6 per cent at $z \le 0.2$ and about 30 per cent at $0.2 < z \le 0.5$. These fractions change significantly if we apply the selection reported in Donley et al. (2012) which is able to better discriminate pure bona fide AGNs from samples that are contaminated by low- and high-redshift star- forming galaxies as the one selected by Lacy’s criterion. We find that only 3 per cent of our total sample is identified as AGN-dominated by Donley’s criterion, less than 2 per cent at $z \le 0.2$ and 4 per cent at $0.2 < z \le 0.5$.
<table><caption>**Table 2.** List of the SEDs used to perform the SED fitting analysis in the IR/submm. The ‘Spectral type’ column shows the grouping procedure we implemented in order to collect together those SED classes with similar properties in terms of FIR colours.</caption><thead><tr><th>Index</th><th>SED class</th><th>Spectral type</th><th>Reference</th></tr></thead><tbody><tr><td>01</td><td>Ell13</td><td>Elliptical</td><td>Polletta+07</td></tr><tr><td>02</td><td>Ell5</td><td>Elliptical</td><td>Polletta+07</td></tr><tr><td>03</td><td>Ell2</td><td>Elliptical</td><td>Polletta+07</td></tr><tr><td>04</td><td>S0</td><td>Spiral</td><td>Polletta+07</td></tr><tr><td>05</td><td>Sa</td><td>Spiral</td><td>Polletta+07</td></tr><tr><td>06</td><td>Sb</td><td>Spiral</td><td>Polletta+07</td></tr><tr><td>07</td><td>Sc</td><td>Spiral</td><td>Polletta+07</td></tr><tr><td>08</td><td>Sd</td><td>Spiral</td><td>Polletta+07</td></tr><tr><td>09</td><td>Sdm</td><td>Spiral</td><td>Polletta+07</td></tr><tr><td>10</td><td>Spi4</td><td>Spiral</td><td>Polletta+07</td></tr><tr><td>11</td><td>N6090</td><td>Starburst</td><td>Polletta+07</td></tr><tr><td>12</td><td>M82</td><td>Starburst</td><td>Polletta+07</td></tr><tr><td>13</td><td>Arp220</td><td>Starburst</td><td>Polletta+07</td></tr><tr><td>14</td><td>I20551</td><td>Starburst</td><td>Polletta+07</td></tr><tr><td>15</td><td>I22491</td><td>Starburst</td><td>Polletta+07</td></tr><tr><td>16</td><td>N6240</td><td>Starburst</td><td>Polletta+07</td></tr><tr><td>17</td><td>Sey2</td><td>Obscured AGN</td><td>Polletta+07</td></tr><tr><td>18</td><td>Sey18</td><td>Obscured AGN</td><td>Polletta+07</td></tr><tr><td>19</td><td>I19254</td><td>Obscured AGN</td><td>Polletta+07</td></tr><tr><td>20</td><td>QSO2</td><td>Unobscured AGN</td><td>Polletta+07</td></tr><tr><td>21</td><td>Torus</td><td>Unobscured AGN</td><td>Polletta+07</td></tr><tr><td>22</td><td>Mrk231</td><td>Obscured AGN</td><td>Polletta+07</td></tr><tr><td>23</td><td>QSO1</td><td>Unobscured AGN</td><td>Polletta+07</td></tr><tr><td>24</td><td>BQSO1</td><td>Unobscured AGN</td><td>Polletta+07</td></tr><tr><td>25</td><td>TQSO1</td><td>Unobscured AGN</td><td>Polletta+07</td></tr><tr><td>26</td><td>Sb</td><td>Spiral</td><td>Gruppioni+10</td></tr><tr><td>27</td><td>Sdm</td><td>Spiral</td><td>Gruppioni+10</td></tr><tr><td>28</td><td>Sey2</td><td>Obscured AGN</td><td>Gruppioni+10</td></tr><tr><td>29</td><td>Sey18</td><td>Obscured AGN</td><td>Gruppioni+10</td></tr><tr><td>30</td><td>Mrk231</td><td>Obscured AGN</td><td>Gruppioni+10</td></tr><tr><td>31</td><td>qso_high</td><td>Unobscured AGN</td><td>Gruppioni+10</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Typical LE PHARE SED fits. The two best-fitting SEDs used to fit short- and long-wavelength photometry are shown by the red and magenta solid line, respectively. The black solid circles are the photometric data used to perform the fit. The ID and the redshift of the source are reported on top of each panel.</figcaption></figure>

# 3 STATISTICAL METHODS

Accurately estimating the LF is difficult in observational cosmology since the presence of observational selection effects like flux detection thresholds can make any given galaxy survey incomplete and thus introduce biases into the LF estimate.

Numerous statistical approaches have been developed to overcome this limit, but, even though they all have advantages, it is only by comparing different and complementary methods that we can be confident about the reliability of our results. For this reason, to estimate the LLFs in the SPIRE bands reported in this paper, we exploit different LF estimators: the $1/V_{max}$ approach of Schmidt (1968) and the modified version $\phi_{est}$ of Page & Carrera (2000); the Bayesian parametric maximum likelihood (ML) method of Kelly, Fan & Vestergaard (2008) and Patel et al. (2013); and the semi-parametric approach of Schafer (2007). All these methods are explained in the following sections.

## 3.1 $1/V_{max}$ Estimator

Schmidt (1968) introduced the intuitive and powerful $1/V_{max}$ estimator for LF evaluation. The quantity $V_{max}$ for each object represents the maximum volume of space which is available to such an object to be included in one sample accounting for the survey flux limits and the redshift bin in which the LF is estimated. $V_{max}$ thus depends on the distribution of the objects in space and the way in which detectability depends on distance. Once the $V_{max}$ (or $V_{max}(L_i)$, since it depends on the luminosity of each object) is defined, the LF can be estimated as

$$ \Phi(B_{j-1} < L \le B_j) = \sum_{B_{j-1} < L \le B_j} \frac{1}{V_{max}(L_i)}, \quad (1) $$

in which its value is computed in bins of luminosity, within the boundary luminosities of a defined bin $[B_{j-1}, B_j]$. It is usually expressed in the differential form as

$$ \phi_{1/V_{max}}(L, z) = \frac{1}{\Delta L} \sum_{i=1}^{N} \frac{1}{V_{max,i}}, \quad (2) $$
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top: relation between rest-frame SPIRE 250 µm or PACS 100 µm luminosities and IR bolometric luminosity colour-coded as a function of redshift. Middle: relations between rest-frame SPIRE 250 µm/PACS 100 µm luminosities and IR bolometric luminosity colour-coded according to the SED best-fitting class obtained by the SED fitting procedure following the list reported in Table 2. Bottom: SPIRE 250 µm and PACS 100 µm k-corrections as a function of redshift colour-coded according to the SED best-fitting class.</figcaption></figure>

where N is the number of objects within some volume–luminosity region. Errors in the LF can be evaluated using Poisson statistics:

$$ \sigma_{\phi(L)}^2 = \sum_{B_{j-1} < L \le B_j} \frac{1}{(V_{\max}(L_i))^2} \quad (3) $$

In our case, there are three main selection factors that may constrain the $V_{\max}$ for each object in our sample: the limit in r magnitude that guides the photometric redshift estimates in the SDSS survey, $r_{AB} < 22.2$; the MIPS 24 µm flux limit that guides the SPIRE 250 µm extraction, $S_{24} > 300$ µJy; and finally the flux density limit
<figure><img src="image_7.png" /><figcaption>**Figure 5.** IRAC colour–colour plot as from Lacy et al. (2004) and Donley et al. (2012). In the left- and right-hand panels, respectively, the complete samples of sources in the two redshift ranges 0.02 < z ≤ 0.2 and 0.2 < z ≤ 0.5 are reported. In both panels, the overplotted solid green circles and blue open circles show the AGN-like objects selected using the Lacy et al. (2004) and Donley et al. (2012) criteria, respectively, while the solid red circles represent the rest of the sample in each redshift bin.</figcaption></figure>

in the SPIRE 250 µm band, $S_{250} > 30$ mJy. Moreover, since we estimate the $1/V_{max}$ in a number of redshift bins, the $V_{max}$ value is actually also limited by $z_{min}$ and $z_{max}$ for each z-bin. Taking into account all these considerations, the $V_{max}$ estimator used in equation (2) is described by

$$V_{\max} = \frac{\Omega}{4\pi} \int_{z_{\min}}^{z_{\max}} dz \frac{dV}{dz}, \qquad (4)$$

where $z_{min}$ and $z_{max}$ are the redshift boundaries resulting from taking into account both the redshift bin range and the selection factors

$$z_{k,\min} = z_{\text{bin}_k,\min} \qquad (5)$$

$$z_{k,\max} = \min[z_{0,\max}, \dots, z_{n,\max}, z_{\text{bin}_k,\max}] \qquad (6)$$

for all 0, ..., n selection factors and for each k redshift bin. For instance, in the case of the SPIRE 250 µm LF estimate in z-bin 0.02 < z < 0.1, the conditions just shown become

$$z_{0.02<z<0.1,\min} = 0.02 \\ z_{0.02<z<0.1,\max} = \min[z_{r_{AB},\max}, z_{f24,\max}, z_{f250,\max}, 0.1],$$

where $z_{r_{AB},\max}$, $z_{f24,\max}$ and $z_{f250,\max}$ are the redshift at which a source in the sample reaches the SDSS $r_{AB}$ magnitude limit (= 22), the 24 µm flux limit (= 300 µJy) and the SPIRE 250 µm limit (= 30 mJy), respectively; 0.02 and 0.1 are the minimum and the maximum of the redshift bin.

This method implies binning of the luminosity data, a non-parametric technique, and as such does not need to assume an analytic form. It does however contain the underlying assumption that galaxies have a uniform distribution in space. In principle, this could be tested with the $V/V_{max}$ distribution, but that still remains difficult when there are multiple selection factors limiting the sample.

The simple $V_{max}$ estimator has evolved, being improved and refined over the years to accommodate the many different types of survey that have steadily grown in size and complexity. One of these approaches is the one implemented in Page & Carrera (2000), the so-called $V_{est}$ method, which we also used here to check whether with our $1/V_{max}$ estimates we are ignoring any important incompleteness factor in our sample. Page & Carrera (2000) improved

the method to take into account systematic errors in the $V_{max}$ test introduced for objects close to the flux limit of a survey. This new method defines the value of the LF $\phi(L)$ as $\phi_{est}$, which assumes that $\phi$ does not change significantly over the luminosity and redshift intervals $\Delta L$ and $\Delta z$, respectively, and is defined as

$$\phi_{\text{est}} = \frac{N}{\int_{L_{\min}}^{L_{\max}} \int_{z_{\min}}^{z_{\max}(L)} \frac{dV}{dz} dz dL}, \qquad (7)$$

where N is the number of objects within some volume–luminosity region.

Due to how the methods work in practice, for LFs in most of the redshift intervals, the two will produce the same results, particularly for the highest luminosity bins of any given redshift bin. However, for the lowest luminosity objects in each redshift bin, which are close to the survey limit and occupy a portion of volume–luminosity space much smaller than the rectangular $\Delta L \Delta z$ region, the two methods can produce the most discrepant results. Nevertheless, in our case we do not find any substantial differences between the $1/V_{max}$ and $1/V_{est}$ solutions, as shown in the following sections.

## 3.2 Bayesian parametric maximum likelihood estimator

The maximum likelihood estimator (MLE) has first been applied in studies of observational cosmology by Sandage, Tammann & Yahil (1979), the so-called STY estimator. In ML analysis, one is interested in finding the estimate that maximizes the likelihood function of the data. For a given statistical model, parametrized by $\theta$, the likelihood function, $p(x|\theta)$, is the probability of observing the data, denoted by x, as a function of the parameters $\theta$. In Bayesian analysis, one attempts to estimate the probability distribution of the model parameters, $\theta$, given the observed data x. Bayes' theorem states that the probability distribution of $\theta$ given x is related to the likelihood function as

$$p(\theta|x) \propto p(x|\theta)p(\theta), \qquad (8)$$

where $p(x|\theta)$ is the likelihood function of the data and the term $p(\theta)$ is the prior probability distribution of $\theta$; the result, $p(\theta, x)$, is called the posterior distribution. The prior distribution, $p(\theta)$, should convey information known prior to the analysis. In general, the
prior distribution should be constructed to ensure that the posterior distribution integrates to 1, but does not have a significant effect on the posterior. In particular, the posterior distribution should not be very sensitive to the choice of prior distribution, unless the prior distribution is constructed with the purpose of placing constraints on the posterior distribution that are not conveyed by the data. The contribution of the prior to $p(\theta|x)$ should become negligible as the sample size becomes large.

From a practical standpoint, the primary difference between the ML approach and the Bayesian approach is that the former is concerned with calculating a point estimate of $\theta$, while the latter is concerned with mapping out the probability distribution of $\theta$ in the parameter space. The ML approach uses an estimate of the sampling distribution of $\theta$ to place constraints on the true value of $\theta$. In contrast, the Bayesian approach directly calculates the probability distribution of $\theta$, given the observed data, to place constraints on the true value of $\theta$.

In terms of LF evaluation, the LF estimate is related to the probability density of $(L, z)$

$$p(L, z) = \frac{1}{N} \phi(L, z) \frac{dV}{dz}, \qquad (9)$$

where N is the total number of sources in the observable Universe and is given by the integral of $\phi$ over L and $V(z)$. The quantity $p(L, z)dLdz$ is the probability of finding a source in the range L, L+dL and z, z+dz. Equation (9) separates the LF into its shape, given by $p(L, z)$, and its normalization, given by N. Once we have an estimate of $p(L, z)$, we can easily convert this to an estimate of $\phi(L, z)$ using equation (9).

In general, it is easier to work with the probability distribution of L and z instead of directly with the LF, because $p(L, z)$ is more directly related to the likelihood function. The function $\phi(L, z)$ can be described, as we have seen, by a parametric form with parameter $\theta$, so that we can derive the likelihood function for the observed data. The presence of flux limits and various other selection effects can make this difficult, since the observed data likelihood function is not simply given by equation (9). In this case, the set of luminosities and redshifts observed by a survey gives a biased estimate of the true underlying distribution, since only those sources with L above the flux limit at a given z are detected. In order to derive the observed data likelihood function, it is necessary to take the survey's selection method into account. This is done by first deriving the joint likelihood function of both the observed and unobserved data, and then integrating out the unobserved data. The probability $p(L, z)$ (as reported in Patel et al. 2013) then becomes

$$p(L, z|\theta) = \frac{\phi(L, z|\theta)p(\text{selected}|L, z) \, dV}{\lambda \, dz}, \qquad (10)$$

where $p(\text{selected}|L, z)$ stands for the probability connected with the selection factors of the survey and $\lambda$ is the expected number of sources, determined by

$$\lambda = \iint \phi(L, z|\theta) p(\text{selected}|L, z) d\log L \frac{dV}{dz} dz, \qquad (11)$$

where the integrals are taken over all possible values of redshift and luminosity.

This last equation gives the expected number of objects in a sample composed of sources of the same morphological type and collected in a single-field survey. For our purposes, we have to change the equation to the following:

$$\lambda = \sum_{\text{SED fields}} \iint \phi(L, z|\theta) p(\text{selected}|L, z) d\log L \frac{dV}{dz} dz, \qquad (12)$$

where we sum together the expected number of sources for each SED type, used for the SED fitting procedure, and survey areas that compose our HerMES Wide Fields sample.

Since the data points are independent, the likelihood function for all N sources in the Universe would be

$$p(L, z|\theta) = \prod_{i=1}^{N} p(L_i, z_i|\theta). \qquad (13)$$

Indeed, we do not know the luminosities and redshifts for all N sources, nor do we know the value of N, since our survey only covers a fraction of the sky and is subject to various selecting criteria. As a result, our survey only contains n sources. For this reason, the selection process must also be included in the probability model, and the total number of sources, N, is an additional parameter that needs to be estimated. Then the likelihood becomes

$$p(n|\theta) = p(N, \{L_i, z_i\}|\theta) = p(N|\theta)p(\{L_i, z_i\}|\theta), \qquad (14)$$

where $p(N|\theta)$ is the probability of observing N objects and $p(\{L_i, z_i\}|\theta)$ is the likelihood of observing a set of $L_i$ and $z_i$, both given the LF model. Is it possible to assume that the number of sources detected follows a Poisson distribution (Patel et al. 2013), where the mean number of detectable sources is given by $\lambda$? Then, the term $p(N, \{L_i, z_i\}|\theta)$ could be written as the product of individual source likelihood function, since each data point is independent:

$$p(N|\theta)p(\{L_i, z_i\}|\theta) = \frac{\lambda^N e^{-\lambda}}{N!} \prod_{i=1}^{N} \frac{\Phi(L, z|\{\theta\}) p(\text{selected}|L, z)}{\lambda} \frac{dV}{dz}. \qquad (15)$$

Then we can use the likelihood function for the LF to perform Bayesian inference by combining it with a prior probability distribution, $p(\theta)$, to compute the posterior probability distribution, $p(\theta|d_i)$, given by Bayes' theorem:

$$p(\theta|d_i) = \frac{p(\{d_i\}|\{\theta\})p(\{\theta\})}{\int p(\{d_i\}|\{\theta\})p(\{\theta\})d\theta}. \qquad (16)$$

The denominator of this equation represents the Bayesian evidence which is determined by integrating the likelihood over the prior parameter space. This last step is needed to normalize the posterior distribution.

Calculating the Bayesian evidence is computationally expensive, since it involves integration over m-dimensions for an m parameter LF model. Therefore, Monte Carlo Markov chain (MCMC) methods, used to examine the posterior probability, perform a random walk through the parameter space to obtain random samples from the posterior distribution. MCMC gives as a result the maximum of the likelihood, but an algorithm is needed to investigate in practice the region around the maximum. Kelly et al. (2008) suggested to use the Metropolis-Hastings algorithm (Metropolis et al. 1953; Hastings 1970) in which a proposed distribution is used to guide the variation of the parameters. The algorithm uses a proposal distribution which depends on the current state to generate a new proposal sample. The algorithm needs to be tuned according to the results and the number of iterations, as well as the parameter step size change. Once we obtain the posterior distribution, we have the best solution for each of the parameters describing the LF model that we have chosen at the beginning; we have the mean value and the standard deviation ($\sigma$) for each of the parameters that we can combine together to find the $\sigma$ of the parametric function chosen as the shape of our LF (see Section 4 for further details on our calculation).
## 3.3 A semi-parametric estimator

Schafer (2007) introduced the semi-parametric method in order to estimate LFs given redshift and luminosity measurements from an inhomogeneously selected sample of objects (e.g. a flux-limited sample). In such a limited sample, like ours, only objects with flux within some range are observable. When this bound on fluxes is transformed into a bound in luminosity, the truncation limits take an irregular shape as a function of redshift; additionally, the k-correction can further complicate this boundary.

We refer the reader to the original paper, Schafer (2007), for a complete description of the method; here we report only the main characteristics of it. This method shows various advantages in comparison with the other techniques previously described: it does not assume a strict parametric form for the LF (differently than the parametric MLE); it does not assume independence between redshifts and luminosities; it does not require the data to be split into arbitrary bins (unlike for the non-parametric MLE), and it naturally incorporates a varying selection function. This is obtained by writing the LF $\phi(z, L)$ as

$$ \log\phi(z, L) = f(z) + g(L) + h(z, L, \theta), \quad (17) $$

where $h(z, L, \theta)$ assumes a parametric form and is introduced to model the dependence between the redshift $z$, the luminosity $L$ and the real valued parameter $\theta$. The functions $f$ and $g$ are estimated in a completely free-form way.

Nevertheless, it is important to notice that this method assumes a complete data set in the untruncated region that requires some care when applying it to samples that may suffer some incompleteness. Discussion on how this issue may influence our results is reported in the later sections.

## 3.4 Parametrizing the LF

Using the classical ML technique (STY), as well the one based on Bayesian statistics, implies the assumption of a parametric form able to describe the LF. This choice is not straightforward and over the years the selected LF models varied. In this work, we decide to use the log-Gaussian function introduced by Saunders et al. (1990) to fit the IRAS IR LF and widely used for IR LF estimates (e.g. Gruppioni et al. 2010, 2013; Patel et al. 2013). Usually, this function is called the modified Schechter function since its formalism is very similar to the one introduced by Schechter (1976). This parametric function is defined as

$$ \Phi(L) = \Phi^* \left( \frac{L}{L^*} \right)^{1-\alpha} \exp \left[ -\frac{1}{2\sigma^2} \log^2 \left( 1 + \frac{L}{L^*} \right) \right], \quad (18) $$

where $\Phi^*$ is a normalization factor defining the overall density of galaxies, usually quoted in units of $h^3 Mpc^{-3}$, and $L^*$ is the characteristic luminosity. The parameter $\alpha$ defines the faint-end slope of the LF and is typically negative, implying relatively large numbers of galaxies with faint luminosities. We also checked whether another functional form was more suitable to describe our LFs, but we did not find any evidence of improvement or substantial differences by using e.g. a double power-law function (used by Rush & Malkan 1993 or Franceschini et al. 2001). We therefore decide to report and discuss the estimates obtained by using only the log-Gaussian function in order to be able to compare our results with other more recent results that use the same parametrization. This approach is well suited to describe the total galaxy population, but may be inadequate if we divide the population into subgroups according, for example, to their optical properties (see Section 5 for more details)

as done by other authors while studying the behaviour of the local mass functions of galaxies (e.g. Baldry 2012).

## 4 RESULTS

We estimate the LFs at SPIRE 250 µm as well as at SPIRE 350 and 500 µm by using the SPIRE 250 µm selected sample and extrapolating the luminosities from the SED fitting results. The higher sensitivity of the SPIRE 250 µm channel with respect to the 350 and 500 µm channels largely ensures that we do not miss sources detected only at these longer wavelengths. Additionally, we estimate the IR bolometric LFs using the integrated luminosity between 8 and 1000 µm and at 24, 70, 100 and 160 µm; these last monochromatic estimates are also used to check our procedure against other published LFs.

In Table 3, we report the values of the best parameter solutions of the parametric Bayesian ML procedure (explained in Section 3.2) using the log-Gaussian functional form (equation 18). In Fig. 6, we report the histograms representing the probability distribution of the best-fitting parameters produced by the MCMC procedures. To obtain these estimates, we run an MCMC procedure with $5 \times 10^6$ iterations. This procedure is a highly time-consuming process; thus, we focused our attention in the most local bin $0.02 < z < 0.1$ of our analysis where we want to obtain a precise estimate of the shape of the local LF observed by Herschel at 250 µm, which is our selection band. Such an estimate represents a fundamental benchmark to study the evolution of the LF (e.g. Vaccari et al., in preparation) as discussed later in Section 5.

As a summary, in Table A1 we report our $1/V_{max}$ LF values for each SPIRE band and the IR bolometric rest-frame luminosity per redshift bin. We exclude from the calculation the sources with $z < 0.02$, as explained in Section 2.2. The error associated with each value of $\Phi$ is estimated following Poissonian statistics, as shown in equation (3).

Since we use photometric redshifts in our sample, we quantify the redshift uncertainties that may affect our results by performing Monte Carlo simulations. We created 10 mock catalogues based on our actual sample, allowing the photometric redshift of each source to vary by assigning a randomly selected value according to the Gaussian SDSS photometric error. For each source in the mock catalogues, we performed the SED fitting and recomputed both the monochromatic and total IR rest-frame luminosities and the $V_{max}$-based LFs, using the randomly varied redshifts. The comparison between our real IR LF solution and the mean derived from the Monte Carlo simulations shows that the uncertainties derived from the use of the photometric redshifts do not significantly change the error bar estimated using the Poissonian approach and mainly alter the lower luminosity bins at the lower redshifts ($z < 0.1$). As an extra test we also check what happens if we estimate the LFs in each

<table><caption>**Table 3.** Best-fitting parameter solution and uncertainties for the local SPIRE 250 µm LF determined using the parametric Bayesian ML procedure. The redshift range for this solution is $0.02 < z < 0.1$.</caption><thead><tr><th>Parameter</th><th>⟨σ⟩</th></tr></thead><tbody><tr><td>log(L<sup>*</sup>) (L<sub>☉</sub>)</td><td>9.03<sup>+0.14</sup><sub>-0.13</sub></td></tr><tr><td>α</td><td>0.96<sup>+0.09</sup><sub>-0.07</sub></td></tr><tr><td>σ</td><td>0.39<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>log(Φ<sup>*</sup>) (Mpc<sup>-3</sup> dex<sup>-1</sup>)</td><td>-1.99<sup>+0.04</sup><sub>-0.02</sub></td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Probability histogram of the best-fitting parameters ($L^*, \alpha, \sigma$ and $\Phi^*$) for the SPIRE 250 µm LLF within $0.02 < z < 0.1$, determined using the MCMC parametric Bayesian procedure performing $5 \times 10^6$ iterations. The highlighted area is the $\pm 1\sigma$ confidence area for each parameter, as reported in Table 3.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** SPIRE 250 µm rest-frame LLF estimates. The black open circles are our $1/V_{\text{max}}$ estimates; the red dashed line is from the Fontanot, Cristiani & Vanzella (2012) model; the beige triple dot-dashed line is from the Negrello et al. (2007) model and the black dot-dashed and dashed lines are LLF prediction at 250 µm from Serjeant & Harrison (2005). The magenta shaded region is the $\pm 1\sigma$ best MCMC solution using the log-Gaussian functional form reported in the text. The magenta line in the right-hand panel is the mean from the MCMC solution plotted with the LF estimates in each field (colour-coded as reported in the legend; the colour-coded number reported in the plot below each field's name is the number of sources in each field in the considered redshift bin).</figcaption></figure>

field using only spectroscopic redshifts and correct the solutions for the incompleteness effect due to this selection. The resulting LFs are effectively undistinguishable and thus confirm that the uncertainties introduced by the use of photometric redshifts are of the order of the Poissonian ones.

Even though the differences are really small, the errors that we quote in Table A1 are the total errors, taking into account both Poissonian and redshift uncertainties associated with $\Phi$.

A summary of the results is reported in the following figures. In Figs 7 and 8, we report the SPIRE 250 µm rest-frame LF estimated
<figure><img src="image_10.png" /><figcaption>**Figure 8.** SPIRE 250 µm rest-frame LLF estimates from field to field. The colour-coded open circles are our $1/V_{\text{max}}$ results for each field (the black is the solution for all the five fields considered together); the red dashed line is the Fontanot et al. (2012) model; the beige triple dot-dashed line is the Negrello et al. (2007) model; the black dot-dashed and dashed lines are LLF predictions at 250 µm from Serjeant & Harrison (2005). Negrello et al. (2007) and Serjeant & Harrison (2005) estimates are reported at the same local ($z=0$) redshift in all panels.</figcaption></figure>
by using the $1/V_{max}$ and the parametric Bayesian ML, reporting both the solutions for the five fields together (see Table 1) and for each field separately. In Fig. 9, we compare our SPIRE 250 µm $1/V_{max}$ LF solution to the H-ATLAS results of Dye et al. (2010). The SPIRE LLFs in different fields do not show any field-to-field variations beyond what is expected from cosmic variance, i.e. about 15 per cent as predicted by theoretical models (Moster et al. 2011). To report the confidence area of our Bayesian ML solution, we estimate the standard deviation of the best-fitting model using the following equation:

$$ \sigma_{\Phi(x_1, x_2, \dots, x_n)}^2 = \sum_{j=1}^{n} \left( \frac{\partial \Phi}{\partial x_j} \sigma_{x_j} \right)^2 + 2 \sum_{j=1}^{n} \sum_{k=j+1}^{n} r_{x_j x_k} \left( \frac{\partial \Phi}{\partial x_j} \sigma_{x_j} \right) \left( \frac{\partial \Phi}{\partial x_k} \sigma_{x_k} \right). \quad (19) $$

This equation represents the general formula for the parametric standard deviation in the case of non-independent variables. The functional form of $\Phi$ is, as already stated, the log-Gaussian function described in equation (18), in which the parameters $L^*$, $\alpha$, $\sigma$ and $\Phi^*$ are in fact not independent from each other. Thus, $\Phi(x_1, x_2, \dots, x_n)$ reported in equation (19) can be translated, into our specific case, as $\Phi(L^*, \alpha, \sigma, \Phi^*)$, while $\sigma_{x_j}$ expresses the error associated with the $j$th parameter in the sum (and the same with $\sigma_{x_k}$ for the $k$th parameter).

In Figs 10 and 11 we report the SPIRE 250 µm luminosity distribution and the SPIRE 250 µm rest-frame LF estimated by using the semi-parametric method described in Section 3.3 and the modified $1/V_{max}$ estimates from Page & Carrera (2000) described in Section 3.1. In Figs 12, 13 and 14, we report the SPIRE 350/500 µm and IR bolometric rest-frame LFs, respectively. Finally, as a check on the robustness of our SPIRE 250 µm selected sample, we estimate the LFs also at other wavelengths, namely MIPS 24/70/160 µm and PACS 70/100/160 µm, and compare our results to others already published. In Fig. 15, we report the 24/70/90/160 µm rest-frame LFs compared with local predictions at these wavelengths given by different authors.

## 4.1 The IR LLD and the IR local SED

Once we obtain our LF solutions in each redshift bin and for each band, we can integrate them to find the luminosity density per redshift bin which is connected to the amount of energy emitted by the galaxies at each wavelength and at each instant. To obtain this information, we perform a $\chi^2$ fit to our $1/V_{max}$ estimates, using the modified Schechter function described in equation (18). Since we are limited to a local sample, at $z > 0.2$ we do not populate the low-luminosity bins of our LFs and for this reason we cannot really constrain the integration at higher redshift. We thus report in Figs 16, 17 and in Table 4 our luminosity density estimates for the SPIRE 250/350/500 µm and the IR bolometric luminosity within $z < 0.2$, reporting the results for three redshift bins whose mean redshifts are 0.05, 0.1 and 0.15.

In Fig. 18, we report the conversion of our luminosity density estimates at SPIRE 250/350/500 µm, as well at MIPS 24/70/160 µm wavelengths to the energy output and we compare our result to those reported by Driver et al. (2012). Our plotted estimates, together with others extrapolated at 90 and 170 µm, are reported in Table 5.

We find that, even though our sample is selected at 250 µm, we can reproduce the energy density at all the other considered FIR bands in the very local Universe. This confirms the shape of the

energy density published by Driver et al. (2012) estimated using the Galaxy And Mass Assembly survey I (GAMA I) data set combined with GALEX, SDSS and UKIRT.

## 4.2 The local SFR

The estimate of the LLF in the SPIRE bands is of fundamental importance for studying the evolution of the SPIRE LFs at higher redshift. In practice, LLF estimates guide the priors on the parameters that define the LF shape that is adopted when fitting the LF also at higher redshifts (Vaccari et al., in preparation). Additionally, thanks to the large volume sampled by shallow and wide-area surveys, these estimates allow us to calculate the star formation rate density (SFRD) in the local Universe with small uncertainties. By integrating the LF in different redshift bins, whenever the observed bands are related to the emission of the young stellar populations, like in this case, we can estimate the SFR at those redshifts. In this context, we can easily use the IR bolometric luminosity as a tracer of SFR and thus the IR bolometric luminosity density as a tracer of the SFRD.

We thus fit our $1/V_{max}$ LLF estimates with a modified Schechter function described in equation (18), obtaining the estimate of the LLD reported in Table 4. The lower and upper limits that we used in the LF integration to estimate the LLDs are $L = 10^8$ and $10^{14} L_\odot$, respectively. These limits guarantee that we account for the bulk of the IR luminosity emitted by our sources. We then convert the estimate of the luminosity density into SFRD using the Kennicutt (1998) relation [assuming a Salpeter initial mass function (IMF)]: $\psi(t) = \text{SFR} = k(\lambda)L(\lambda)$, where $k(\text{IR}) = 4.5 \times 10^{-44} [\text{M}_\odot \text{yr}^{-1} \text{W Hz}]$.

We used our SED fitting analysis and the IRAC colour-colour criteria by Lacy et al. (2004) and Donley et al. (2012) to quantify the possible AGN contamination in our sample, as discussed in Section 2.4. We find that in our sample the fraction of objects showing AGN-like IRAC colours and AGN-like SEDs is very small and even if we discard from our results the total luminosity contribution of these sources, our LFs and thus SFR estimates do not significantly deviate from the results obtained using our total sample. Even for these AGN-like sources (mainly located above $z \sim 0.25$), the vast majority of the IR luminosity is still contributed by dust emission associated with ongoing star formation. This is also confirmed by Hatziminaoglou, Fritz & Jarrett (2009), Hatziminaoglou et al. (2010) and Bothwell et al. (2011) that show how AGN contribution to the FIR emission of the general extragalactic population is rather small. For these reasons, we conclude that the AGN contribution does not significantly affect our LF and SFRD estimates.

The SFRD estimates we obtain from the IR bolometric luminosity density (estimated at $0.02 < z < 0.1$, $0.05 < z < 0.15$ and $0.1 < z < 0.2$) are reported in Table 6, together with other SFRD estimates obtained by various authors using different SFR tracers (all the results are converted to the same IMF and cosmology). These same data are also shown in Fig. 11. The uncertainties reported in Table 6 are percentage errors.

# 5 DISCUSSION

Using some of the widest area surveys performed by Spitzer and Herschel, in this paper we have studied in details the LLFs of SPIRE sources. Our LLFs at 250/350/500 µm (SPIRE) strongly constrain the LLD of the Universe throughout the FIR/submm wavelength range.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** SPIRE 250 µm rest-frame LLFs compared to the H-ATLAS estimate from Dye et al. (2010). The black open circles are our $1/V_{max}$ results; the blue open diamonds are the SDP H-ATLAS SPIRE 250 µm rest-frame LLF from Dye et al. (2010); the red open triangles are the SDP HerMES SPIRE 250 µm rest-frame LLF of Vaccari et al. (2010); the black open triangles are the SDP HerMES SPIRE 250 rest-frame LF of Eales et al. (2010); the red dashed line is the SPIRE 250 µm LF predicted by Fontanot et al. (2012); the beige triple dot-dashed line is the Negrello et al. (2007) model; the black dot-dashed and dashed lines are the LLF prediction at 250 µm from Serjeant & Harrison (2005). Negrello et al. (2007) and Serjeant & Harrison (2005) estimates are reported at the same local ($z=0$) redshift in all panels.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** SPIRE 250 µm luminosity distribution versus redshift plane, as reconstructed using the Schafer (2007) estimator. The red points are the data; the red dashed lines mark the flux limitations adopted in the application of the semi-parametric LF estimator by Schafer (2007) and the solid black lines are isodensity contours corresponding to the semi-parametric reconstructions of the source volume density as a function of luminosity and redshift.</figcaption></figure>

Our estimates mostly confirm and improve upon the HerMES SDP results published in Vaccari et al. (2010), thanks to our increased statistics; this is particularly visible in the 500 µm LF solution, which shows strongly reduced uncertainties. Dye et al. (2010) used Herschel SDP data to compute the H-ATLAS (Eales et al. 2010) SPIRE LLF. This analysis was carried out very early during the Herschel mission and relied on shallower SPIRE observations, much fewer ancillary data and smaller coverage area than our own. We thus judge that the H-ATLAS analysis is likely to have suffered from detection and cross-identification incompleteness and that the discrepancy we find between our results and theirs is therefore likely to be due to the H-ATLAS analysis. Our results are in fact broadly in agreement with their estimates in the highest luminosity bins where the uncertainties on the H-ATLAS SPIRE flux estimates were possibly smaller and their sample more complete, but at the lowest luminosities and redshifts they found LF values lower by up to 50 per cent, as shown in Fig. 9.

The semi-parametric method for the LF estimate of Schafer (2007, see Fig. 11) is in perfect agreement with other classical estimators at low redshifts, $z < 0.2–0.3$. At higher redshifts, the agreement becomes poorer, being acceptable at high luminosities but degrading at lower luminosity values, where the semi-parametric estimate is always in excess of the $1/V_{max}$ values. The precise origin of this problem is not fully understood, but it is clear that it happens in regions of the data space that are poorly sampled by the observations or where the data are scattered e.g. by the effects of the K-correction.

From the IR bolometric LF, we can estimate the SFRD of the local Universe in various redshift bins. In Fig. 17, we report our SFRD solutions and we compare them to others already published in the same redshift range. We see a large scatter in the local SFRD estimates using different SFR diagnostics. In particular, the H$\alpha$ measurements present the largest scatter between different published results. Our new data are entirely consistent with Vaccari et al. (2010) and show good agreement also with OII-based esti-

<figure><img src="image_13.png" /><figcaption>**Figure 11.** SPIRE 250 µm rest-frame LLF estimated using the semi-parametric method of Schafer (2007) and the modified $1/V_{max}$ approach of Page & Carrera (2000). Our classic $1/V_{max}$ estimate is shown in grey; in red is the estimate using the Page & Carrera (2000) method and in black the estimate using the Schafer (2007) approach.</figcaption></figure>

mates (except perhaps at $z = 0.2$ by Hogg et al. 1998). Instead, our SFRD based on the FIR/submm bolometric flux is systematically lower than the radio 1.4 GHz estimates and those combining IR and UV data by Martin et al. (2005) and Bothwell et al. (2011). In principle, the radio flux should be unaffected by dust extinction and thus a more faithful representation of the total SFR than either the IR or IR+UV values. Nevertheless, the radio flux can be more affected by the AGN activity than the IR/submm ones. If we include the UV-uncorrected portion of the SFRD mapped by short-wavelength UV spectral data to our FIR estimate, we find that our total SFRD UV+IR is comparable, within the errors, to the radio estimates, thus confirming that the UV+IR SFRD estimate is a good proxy for the total SFRD in the local Universe and the contamination from AGN in the radio derivation is negligible.

The analysis reported in this paper represents a fundamental local benchmark to study the evolution of the LF and, consequently, of the derived SFR with cosmic time. Studying the evolution of the LF requires very deep data that are then limited to very small areas of the sky and thus it is difficult to constrain the local shape of the LF where a large statistical sample of local galaxies (like ours) is required. This can be seen in Fig. 14 where we compare our local analysis with the one done using the deep COSMOS data (area 1.7 deg$^2$ and flux limited $S_{250} > 10$ mJy) that will be reported in Vaccari et al. (in preparation). Only the large area surveyed by our sample enables us to really study the local shape of the LF, while the deep sample allows us to populate only a few luminosity bins. On the other hand, deep data become more and more important with
<figure><img src="image_14.png" /><figcaption>**Figure 12.** SPIRE 350 µm rest-frame LLF estimates. The black open circles are our $1/V_{\text{max}}$; the red open triangles are the SDP HerMES SPIRE 350 µm rest-frame LLF from Vaccari et al. (2010); the green open triangles are the *Planck* 857 GHz or 350 µm LLF estimate from Negrello et al. (2013); the black dot-dashed and dashed lines are LLF prediction at 350 µm from Serjeant & Harrison (2005).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** SPIRE 500 µm rest-frame LLF estimates. The black open circles are our $1/V_{\text{max}}$; the red open triangles are the SDP HerMES SPIRE 500 µm rest-frame LLF from Vaccari et al. (2010); the green open triangles are the *Planck* 545 GHz or 550 µm LLF estimate from Negrello et al. (2013) converted to our wavelength by using a spectral index of $\alpha = 2.7$; the black dot-dashed and dashed lines are LLF prediction at 500 µm from Serjeant & Harrison (2005).</figcaption></figure>

increasing redshift where our sample soon starts being limited to the higher luminosity bins.

Our LF estimates show significant and rapid luminosity evolution already at low redshifts. In Fig. 19, we report our results about the redshift evolution of the parameters expressing the spatial density dependence of the LFs ($\Phi^*$) and the luminosity dependence ($L^*$) estimated for the IR bolometric at the 250 µm LFs. We found positive evolution in luminosity and negative evolution in density with $L_{\text{IR}}^* \propto (1+z)^{0.6 \pm 0.4}$, $\Phi_{\text{IR}}^* \propto (1+z)^{-2.1 \pm 0.4}$ for the IR bolometric LF and $L_{250}^* \propto (1+z)^{5.3 \pm 0.2}$, $\Phi_{250}^* \propto (1+z)^{-0.6 \pm 0.4}$ for the 250 µm LF. The high values of the evolution rates that we find (both positive and negative) for the luminosity and density parameters are however consistent with previous results based on previous and more limited data sets from *Spitzer* (Patel et al. 2013) and from *IRAS* (Hacking, Condon & Houck 1987; Lonsdale et al. 1990). Similar,

although slightly lower, trends for positive luminosity and negative density evolution are found by Gruppioni et al. (2013). Gruppioni et al. (2013) used a sample deeper and over a much smaller area than ours. Their sample includes sources as faint as ours but they are very few in the local Universe since they suffer from a small sample variance due to the little areas targeted. For this reason, we are able to get a more accurate estimate of the LFs down to similar luminosities in the local Universe.

Interesting for our analysis is the comparison with Negrello et al. (2013) reported in Figs 12 and 13 which show a steep LF in the lowest luminosity bins while our estimate remains flat down to $L_{350} \sim 10^8$ and $L_{500} \sim 10^7$ L⊙, respectively. In general, our low-$z$ LFs are computed at $z > 0.02$, while the *Planck* sources used by Negrello et al. (2013) are located at a mean redshift value of $z \sim 0.01$. This means that our analysis is based on a deeper sample somehow
<figure><img src="image_16.png" /><figcaption>**Figure 14.** The IR bolometric rest-frame LLFs. The black open circles are our $1/V_{\text{max}}$ results; the green open circles are the $1/V_{\text{max}}$ results using COSMOS data (area 1.7 deg$^2$ and flux limited $S_{250} > 10$ mJy; Vaccari et al., in preparation); the blue open squares are the SWIRE IR bolometric rest-frame LF of Patel et al. (2013) using an MIPS 70 and 160 µm selected sample in LH and XMM-LSS; the red open triangles are the IR bolometric rest-frame LF estimate of Vaccari et al. (2010); the red dashed line is the IR bolometric LF predicted of Fontanot et al. (2012); the pink open diamonds are the IRAS IR bolometric rest-frame LF of Sanders et al. (2003); the beige triple dot-dashed line is Negrello et al. (2007) model; the black dotted line is the Valiante et al. (2009) model. Sanders et al. (2003), Negrello et al. (2007) and Valiante et al. (2009) estimates are reported at the same local ($z=0$) redshift in all panels.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 15.** MIPS 24/70/60 µm and IRAS 60 µm LFs as derived by our SPIRE 250 µm sample. The black open circles are our $1/V_{max}$ in all the panels. Top left: the MIPS 24 µm LF estimate. The open red squares are the IRAS 25 µm LF from Shupe (1998); the open green triangles are the MIPS 24 µm LF from Marleau et al. (2007); the open pink hexagons are the MIPS 24 µm LF of Babbedge (2006); the blue asterisks are the 25 µm from IIFSCz by Wang et al. (2009) converted to MIPS 24 µm; the open light blue pentagons are the MIPS 24 µm LF from Rodighiero et al. (2010). Top right: the IRAS 60 µm LF estimate. The open green squares are the IRAS 60 µm LF from Saunders et al. (1990). Bottom left: the MIPS 70 µm LF estimate. The open blue squares are the MIPS 70 µm LF of Patel et al. (2013); the dot-dashed and dashed line are the LF estimates from Serjeant & Harrison (2005). Bottom right: the MIPS 160 µm LF estimate. The open blue squares are the MIPS 160 µm LF of Patel et al. (2013); the open black triangles are the ISO 170 µm LF from Takeuchi (2006) converted to MIPS 160 µm; the dot-dashed and dashed line are the LF estimates from Serjeant & Harrison (2005).</figcaption></figure>

complementary to the Planck's one. Our sample therefore does not suffer from contamination from either the Local Supercluster or the Virgo cluster [like Planck and thus potentially the Negrello et al. (2013) estimates] while representing the LF of typical galaxies in the not-so-nearby Universe (unlike Planck). Moreover, it can be argued that our measurement averages over any local inhomogeneity by sampling a larger cosmic volume than Planck (~10 times larger at $z \sim 0.2$ over 39 deg$^2$ than Planck at $z \sim 0.01$ over 30 000 deg$^2$).

Indeed, over a much smaller area, but with a much deeper sample, the flatness of the slope is also confirmed by Gruppioni et al. (2013) when measuring the $0 < z < 0.3$ IR LF. In any case, at values of $L_{350}$ brighter than $\sim 10^8 L_{\odot}$ and $L_{500}$ brighter than $\sim 10^7 L_{\odot}$, where we are ~100 per cent complete and where the Planck sample is less affected by the presence of local structures and inhomogeneities, we find that our results are in overall agreement with Negrello et al. (2013) at both 350 and 500 µm. Similar considerations can
<figure><img src="image_18.png" /><figcaption>**Figure 16.** SPIRE 250/350/500 µm rest-frame LFs evolution within 0.02 < z < 0.5 along with the luminosity density estimates. Left: the SPIRE 250/350/500 µm LF evolution within 0.02 < z < 0.5. The colour-coded filled points are our 1/Vmax solution in each redshift bin while the solid curves represent the best-fitting solution to the first three redshift bins reported in the legend by using a modified Schechter function whose best-fitting parameters are reported in the panel. Right: SPIRE 250/350/500 µm luminosity density resulting from the fit of the LFs in the first three redshift bins reported on the left.</figcaption></figure>

be made when we compare our IR bolometric LF with the previous estimate obtained by Sanders et al. (2003), which appears to be slightly steeper than ours in the lowest luminosity bins (see Fig. 14). The mean and median redshifts of the entire IRAS sample used by Sanders et al. (2003) are in fact $z = 0.0126$ and 0.0082, respectively, and their LF estimate can therefore be affected by the Virgo cluster in the same way as Planck's estimates discussed above.

Our ability to map the local LF with good precision has revealed a wiggle in the shapes of the functions, with a local maximum at $logL_{250} \sim 9.5$ and $logL_{IR} \sim 10.5$, respectively. This feature, which appears relatively stable with wavelength, is reminiscent of similar behaviour found in the local mass functions of galaxies (Baldry 2012; Ilbert et al. 2013; Moustakas 2013), and interpreted as due to the summed contributions of red and blue galaxies, having Schechter functions with different slopes and cutoff masses. Given the known relationship between stellar mass and IR luminosity, it may not come as a surprise that a similar feature appears in our IR LFs. To test this possibility, we have divided our sample into red and blue subpopulations, following the recipe of Baldry (2012),

<table><caption>**Table 4.** LLD estimates in the SPIRE 250/350/500 µm bands and for the IR bolometric luminosity using the local SPIRE sample within 0.02 < z < 0.1. The values are reported as log(LLD) and log/errors), expressed in L⊙ Mpc⁻³.</caption><thead><tr><th rowspan="2">⟨z⟩</th><th colspan="4">Local luminosity density</th></tr><tr><th>log(ρ<sub>L</sub>,σ)<sub>250</sub></th><th>log(ρ<sub>L</sub>,σ)<sub>350</sub></th><th>log(ρ<sub>L</sub>,σ)<sub>500</sub></th><th>log(ρ<sub>L</sub>,σ)<sub>IR</sub></th></tr></thead><tbody><tr><td>0.05</td><td>7.11, 0.02</td><td>6.64, 0.02</td><td>6.09, 0.02</td><td>7.92, 0.02</td></tr><tr><td>0.10</td><td>7.23, 0.02</td><td>6.75, 0.01</td><td>6.20, 0.01</td><td>8.02, 0.02</td></tr><tr><td>0.15</td><td>7.31, 0.02</td><td>6.82, 0.02</td><td>6.27, 0.02</td><td>8.07, 0.02</td></tr></tbody></table>

and separately calculated the LFs for the two classes. The results, reported in Fig. 20, confirm that red galaxies have an IR LF peaking at $log(L_{250}) \sim 9.5$ and $log(L_{IR}) \sim 10.5$ and decreasing at higher and lower L, while blue galaxies have steep Schechter slopes and lower characteristic luminosities. These are purely observational results and further analysis would be required to better constrain this feature, but this goes beyond the scope of this paper. At any
<figure><img src="image_19.png" /><figcaption>**Figure 17.** The IR bolometric rest-frame LF evolution within $0.02 < z < 0.5$ along with an illustration of the SFRD in the local Universe showing published results and from this work. Left: the IR bolometric LF within $0.02 < z < 0.5$, integrated in the first three redshift bins reported in the legend by using a modified Schechter function. Right: the derived SFRD in the local Universe. Black open circles are our data as a result of the integrations of the LFs on the left converted to SFRD by using the Kennicutt (1998) relation (assuming a Salpeter IMF) and black asterisks are our results plus the contribution of the UV SFRD as estimated by Wyder et al. (2005) at $(z) = 0.05$ and Budavári et al. (2005) at $(z) = 0.1, 0.15$. This sum should represent the total SFRD in the local Universe. The red open diamonds are $O\pi$ estimates by Gallego et al. (2002), Sullivan et al. (2000) and Hogg et al. (1998); the blue open triangles are $H\alpha$ estimates by Gallego et al. (1995), Tresse & Maddox (1998), Sullivan et al. (2000), Pérez-González (2003), Ly et al. (2007), Hanish et al. (2006), Brinchmann et al. (2004), Dale et al. (2010) and Westra et al. (2010); the green open square is radio 1.4 GHz estimates by Serjeant et al. (2002) and Condon 1989; the magenta crosses are FUV+IR estimates by Martin et al. (2005) and Bothwell et al. (2011); the cyan crosses are FUV estimates by Sullivan et al. (2000); the pink open squares are IR estimate from Vaccari et al. (2010); the black dashed line is from Hopkins & Beacom (2006).</figcaption></figure>

<table><caption>**Table 5.** Local energy output of the Universe at different wavelengths.</caption><thead><tr><th rowspan="2">λ<br>(μm)</th><th colspan="2">Local energy output</th></tr><tr><th>ρ<sub>L</sub>(λ) λ<br>(10<sup>33</sup> h W Mpc<sup>-3</sup>)</th><th>ρ<sub>L</sub>(λ) λ<br>(10<sup>33</sup> h W Mpc<sup>-3</sup>)</th></tr></thead><tbody><tr><td>24</td><td>3.91 ± 0.69</td><td></td></tr><tr><td>60</td><td>16.87 ± 3.47</td><td></td></tr><tr><td>70</td><td>22.18 ± 4.77</td><td></td></tr><tr><td>90</td><td>25.93 ± 5.59</td><td></td></tr><tr><td>100</td><td>27.10 ± 5.79</td><td></td></tr><tr><td>160</td><td>19.95 ± 4.27</td><td></td></tr><tr><td>170</td><td>18.54 ± 4.00</td><td></td></tr><tr><td>250</td><td>6.98 ± 1.45</td><td></td></tr><tr><td>350</td><td>2.32 ± 0.46</td><td></td></tr><tr><td>500</td><td>0.58 ± 0.14</td><td></td></tr></tbody></table>

rate, our findings seem to indicate that massive early-type spirals dominate the high-IR-luminosity end of the LF, while bluer lower mass late-type spirals and irregulars dominate its low-luminosity end.

We also performed a preliminary comparison with semi-analytical models of galaxy formation available in the literature, focusing our attention on the redshift range between $z = 0.02$ and 0.2. From these preliminary comparisons, we notice that the Fontanot et al. (2012) predictions (using the MORGANA code by Monaco, Fontanot & Taffoni 2007) seem to broadly reproduce the shape of the LF within the uncertainties, but they underestimate the LF at lower luminosities when compared to our IR bolometric LF estimates. Other predictions done by e.g. Negrello et al. (2007), Serjeant & Harrison (2005) and Valiante et al. (2009) at different wavelengths also show good agreement with our results at higher luminosities, but most of them seem to underestimate the LF when compared to what we obtain at lower luminosities. A more careful and systematic analysis of existing and improved models is required to properly address this issue (e.g. Gruppioni et al. 2015; Franceschini et al., in preparation).

<table><caption>**Table 6.** Star formation rate density in the local Universe: literature results and from this work. This table is an updated version of the one reported in Bothwell et al. (2011). The FUV unobscured SFRD added to our IR results and quoted in this table are from Wyder et al. (2005) at $z = 0.05$ and Budavári et al. (2005) at $z = 0.1$ and 0.15.</caption><thead><tr><th>Reference</th><th>SFR tracer</th><th>⟨z⟩</th><th>SFRD<br>(10<sup>-3</sup> M<sub>☉</sub> yr<sup>-1</sup> Mpc<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Gallego et al. (2002)</td><td>[O II]</td><td>0.025</td><td>9.3 ± 3</td></tr><tr><td>Sullivan et al. (2000)</td><td>[O II]</td><td>0.15</td><td>23 ± 3</td></tr><tr><td>Hogg et al. (1998)</td><td>[O II]</td><td>0.20</td><td>11 ± 4</td></tr><tr><td>Gallego et al. (1995)</td><td>H α</td><td>0.022</td><td>12 ± 5</td></tr><tr><td>Tresse & Maddox (1998)</td><td>H α</td><td>0.2</td><td>25 ± 4</td></tr><tr><td>Sullivan et al. (2000)</td><td>H α</td><td>0.15</td><td>14 ± 3</td></tr><tr><td>Pérez-González (2003)</td><td>H α</td><td>0.025</td><td>25 ± 4</td></tr><tr><td>Ly et al. (2007)</td><td>H α</td><td>0.08</td><td>13 ± 4</td></tr><tr><td>Hanish et al. (2006)</td><td>H α</td><td>0.01</td><td>16<sup>+2</sup><sub>-4</sub></td></tr><tr><td>Brinchmann et al. (2004)</td><td>H α</td><td>0.15</td><td>29 ± 5</td></tr><tr><td>Dale et al. (2010)</td><td>H α</td><td>0.16</td><td>10<sup>+6</sup><sub>-4</sub></td></tr><tr><td>Westra et al. (2010)</td><td>H α</td><td>0.05</td><td>6 ± 2</td></tr><tr><td>Westra et al. (2010)</td><td>H α</td><td>0.15</td><td>12 ± 3</td></tr><tr><td>Serjeant et al. (2002)</td><td>1.4 GHz</td><td>0.005</td><td>21 ± 5</td></tr><tr><td>Condon (1989)</td><td>1.4 GHz</td><td>0.005</td><td>21 ± 0.5</td></tr><tr><td>Sullivan et al. (2000)</td><td>FUV</td><td>0.150</td><td>39 ± 5</td></tr><tr><td>Martin et al. (2005)</td><td>FUV+IR</td><td>0.02</td><td>21 ± 2</td></tr><tr><td>Bothwell et al. (2011)</td><td>FUV+IR</td><td>0.05</td><td>25 ± 1.6</td></tr><tr><td>Vaccari et al. (2010)</td><td>IR</td><td>0.1</td><td>22.3 ± 8.2</td></tr><tr><td>This work</td><td>IR</td><td>0.05</td><td>14.11 ± 2.4</td></tr><tr><td>This work</td><td>IR</td><td>0.10</td><td>18.00 ± 2.9</td></tr><tr><td>This work</td><td>IR</td><td>0.15</td><td>20.10 ± 2.2</td></tr><tr><td>This work</td><td>FUV+IR</td><td>0.05</td><td>19.07 ± 2.4</td></tr><tr><td>This work</td><td>FUV+IR</td><td>0.10</td><td>22.53 ± 2.9</td></tr><tr><td>This work</td><td>FUV+IR</td><td>0.15</td><td>25.42 ± 2.2</td></tr></tbody></table>

## 6 CONCLUSIONS

The determination of the galaxy LF is often hampered by the difficulties of covering a wide area down to faint fluxes on the one hand, and determining counterparts and redshifts for detected sources in a complete and reliable manner on the other hand. In this work,
<figure><img src="image_20.png" /><figcaption>**Figure 18.** The multiwavelength energy output in the local Universe. The LLD at different wavelengths was computed by integrating the relevant monochromatic LLFs over the $0.02 < z < 0.1$ bin. Plotted values of this work, solid black circles, are reported in Table 5.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Evolution of $L^*$ and $\Phi^*$ as a function of $z$ [$L_* \propto (1+z)^{\alpha_L}$ and $\Phi^* \propto (1+z)^{\alpha_D}$] estimated for the LF at 250 µm (left-hand panels) and for the IR bolometric rest-frame LLF (right-hand panels) within $0.02 < z < 0.15$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20.** 250 µm and IR bolometric LLFs for the blue and red populations (reported with blue and red open circles, respectively) compared with the total of the two populations (black open circles) in the redshift range $0.02 < z < 0.2$.</figcaption></figure>
we have thus assembled and exploited the widest area *Spitzer* and *Herschel* extragalactic surveys to select IR galaxy samples in a complete and reliable manner, and the best UV/optical/NIR ancillary data to identify them. Thanks to *Spitzer* and *Herschel* observations, we are now able to reliably sample the IR bolo- metric luminosity of local sources and thus provide important in- sights into dust obscured star formation activity across cosmic time. Even with the best data sets, however, accurately constructing the LF remains a tricky pursuit, since the presence of observational selection effects due to e.g. detection thresholds in apparent mag- nitude, colour, surface brightness or some combination thereof can make any given galaxy survey incomplete and thus introduce bi- ases in the LF estimates. Only a comparison of results coming from different LF estimators applied to the same samples can ensure that we can assess the impact of these biases in a robust manner.

Armed with the *Spitzer* Data Fusion, we were able to describe the $0.02 < z < 0.5$ LLF of sources selected in wide fields by *Herschel* SPIRE imaging. We fully exploited the multiwavelength information collected within the *Spitzer* Data Fusion to perform an SED fitting analysis of SPIRE sources and thus estimate the monochromatic rest-frame luminosities at 250, 350 and 500 µm as well as the IR luminosity between 8 and 1000 µm. We then implemented a number of different statistical estimators to evaluate the LLFs of flux-limited samples in these bands: the classical $1/V_{max}$ estimator of Schmidt (1968) and the modified $1/V_{est}$ version of Page & Carrera (2000); a parametric ML technique based on a Bayesian approach as described in Kelly et al. (2008); and finally a semi- parametric approach introduced by Schafer (2007).

Our high-quality determinations of the IR LFs have revealed for the first time some previously unidentified features in their shape, which we interpret as due to the contributions of red (possibly early- type) and blue (possibly late-type) galaxy populations, with their different Schechter forms. By means of this analysis, we find that the LFs show significant and rapid luminosity evolution already at low redshifts, $0.02 < z < 0.2$. Converting our IR LD estimate into an SFRD, we can determine the SFRD of the local Universe up to redshift 0.2, where the integration of the LF solution is more reliable given that our data set fails to populate the low-luminosity bins of the LF at higher $z$. Summing over our IR SFRD estimate of the un- obscured contribution based on the UV dust-uncorrected emission from local galaxies, we estimate that $SFRD \simeq SFRD_0 + 0.08z$, where $SFRD_0 \simeq (1.9 \pm 0.03) \times 10^{-2} [M_{\odot} Mpc^{-3}]$ is our total SFRD estimate at $z \simeq 0.02$. This analysis represents a local bench- mark for studying the evolution of the IR LF and SFR function with cosmic time.

## ACKNOWLEDGEMENTS

Lucia Marchetti (LM) acknowledges support from the Science and Technology Facilities Council (STFC) under grant ST/J001597/1. Lucia Marchetti, Mattia Vaccari and Alberto Franceschini acknowl- edge support from ASI 'Herschel Science' Contracts I/005/07/1 and I/005/11/0. Mattia Negrello produced additional predictions based on his models. Julie Wardlow acknowledges the Dark Cosmol- ogy Centre funded by the Danish National Research Foundation. Mattia Vaccari acknowledges support from the Square Kilometre Array South Africa project, the South African National Research Foundation and Department of Science and Technology (DST/CON 0134/2014), the European Commission Research Executive Agency (FP7-SPACE-2013-1 GA 607254) and the Italian Ministry for For- eign Affairs and International Cooperation (PGR GA ZA14GR02).

Nicholas Seymour is the recipient of an ARC Future Fellowship. Anna Feltre acknowledges support from the ERC via an Advanced Grant under grant agreement no. 321323-NEOGAL. This work makes use of STILTS http://www.starlink.ac.uk/stilts/ and TOPCAT (Taylor 2005). SPIRE has been developed by a consortium of insti- tutes led by Cardiff University (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); Stockholm Observatory (Sweden); STFC (UK); and NASA (USA). The authors would like to thank the anonymous referee for helpful comments.

## REFERENCES

* Abazajian K. N. et al., 2009, ApJS, 182, 543

* Arnouts S., Cristiani S., Moscardini L., Matarrese S., Lucchin F., Fontana A., Giallongo E., 1999, MNRAS, 310, 540

* Babbedge T. S. R. et al., 2006, MNRAS, 370, 1159

* Baldry I. K. et al., 2012, MNRAS, 421, 621

* Bertin E., Arnouts S., 1996, A&A, 117, 393

* Bolton A. S. et al., 2012, AJ, 144, 144

* Bothwell M. S. et al., 2011, MNRAS, 415, 1815

* Brinchmann J., Charlot S., White S. D. M., Tremonti C., Kauffmann G., Heckman T., Brinkmann J., 2004, MNRAS, 351, 1151

* Budavári T. et al., 2005, ApJ, 619, L31

* Carliles S., Budavári T., Heinis S., Priebe C., Szalay A. S., 2010, ApJS, 712, 511

* Condon J. J., 1989, ApJ, 338, 13

* Csabai I., Dobos L., Trencséni M., Herczegh G., Józsa P., Purger N., Budavári T., Szalay A. S., 2007, Astron. Nachr., 328, 852

* Dale D. A. et al., 2010, ApJ, 712, L189

* Donley J. L. et al., 2012, ApJ, 748, 142

* Driver S. P. et al., 2012, MNRAS, 427, 3244

* Dye S. et al., 2010, A&A, 518, L10

* Eales S. et al., 2009, ApJ, 707, 1779

* Eales S. et al., 2010, PASP, 122, 499

* Eales S. A. et al., 2010, A&A, 518, L23

* Fernandez-Conde N., Lagache G., Puget J.-L., Dole H., 2008, A&A, 481, 885

* Fontanot F., Cristiani S., Vanzella E., 2012, MNRAS, 425, 1413

* Franceschini A., Aussel H., Cesarsky C. J., Elbaz D., Fadda D., 2001, A&A, 378, 1

* Franceschini A., Rodighiero G., Vaccari M., Berta S., Marchetti L., Mainetti G., 2010, A&A, 517, A74

* Frayer D. T. et al., 2006, AJ, 131, 250

* Gallego J., Zamorano J., Aragon-Salamanca A., Rego M., 1995, ApJ, 455, L1

* Gallego J., García-Dabó C. E., Zamorano J., Aragón-Salamanca A., Rego M., 2002, ApJ, 570, L1

* Griffin M. J. et al., 2010, A&A, 518, L3

* Gruppioni C. et al., 2010, A&A, 518, L27

* Gruppioni C. et al., 2013, MNRAS, 432, 23

* Gruppioni C. et al., 2015, MNRAS, 451, 3419

* Hacking P., Condon J. J., Houck J. R., 1987, ApJ, 316, L15

* Hanish D. J. et al., 2006, ApJ, 649, 150

* Hastings W., 1970, Biometrika, 57, 97

* Hatziminaoglou E., Fritz J., Jarrett T. H., 2009, MNRAS, 399, 1206

* Hatziminaoglou E. et al., 2010, A&A, 518, L33

* Hogg D. W., Cohen J. G., Blandford R., Pahre M. A., 1998, ApJ, 504, 622

* Hopkins A. M., Beacom J. F., 2006, ApJ, 651, 142

* Huynh M. T., Frayer D. T., Mobasher B., Dickinson M., Chary R.-R., Morrison G., 2007, ApJ, 667, L9
Ilbert O. et al., 2006, A&A, 457, 841

Ilbert O. et al., 2009, ApJ, 690, 1236

Ilbert O. et al., 2013, A&A, 556, A55

Kelly B. C., Fan X., Vestergaard M., 2008, ApJ, 682, 874

Kennicutt R. C., Jr, 1998, ApJ, 498, 541

Lacy M. et al., 2004, ApJS, 154, 166

Le Borgne D., Elbaz D., Ocvirk P., Pichon C., 2009, A&A, 504, 727

Le Floc'h E. et al., 2005, ApJ, 632, 169

Levenson L. et al., 2010, MNRAS, 409, 83

Lilly S. J., Tresse L., Hammer F., Crampton D., Le Fevre O., 1995, ApJ, 455, 108

Lonsdale C. J., Hacking P. B., Conrow T. P., Robinson M., 1990, ApJ, 358, 60

Lonsdale C. J. et al., 2003, PASP, 115, 897

Ly C. et al., 2007, ApJ, 657, 738

Magnelli B., Elbaz D., Chary R. R., Dickinson M., Le Borgne D., Frayer D. T., Willmer C. N. A., 2009, A&A, 496, 57

Makovoz D., Marleau F. R., 2005, PASP, 117, 1113

Marleau F. R., Fadda D., Appleton P. N., Noriega-Crespo A., Im M., Clancy D., 2007, ApJ, 663, 218

Martin D. C. et al., 2005, ApJ, 619, L59

Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E., 1953, J. Chem. Phys., 21, 1087

Monaco P., Fontanot F., Taffoni G., 2007, MNRAS, 375, 1189

Montero-Dorta A. D., Prada F., 2009, MNRAS, 399, 1106

Moster B. P., Somerville R. S., Newman J. A., Rix H.-W., 2011, ApJ, 731, 113

Moustakas J. et al., 2013, ApJ, 767, 50

Negrello M., Perrotta F., González-Nuevo J., Silva L., de Zotti G., Granato G. L., Baccigalupi C., Danese L., 2007, MNRAS, 377, 1557

Negrello M. et al., 2013, MNRAS, 429, 1309

Oliver S. J. et al., 2012, MNRAS, 424, 1614

Page M. J., Carrera F. J., 2000, MNRAS, 311, 433

Patel H., Clements D. L., Vaccari M., Mortlock D. J., Rowan-Robinson M., Pérez-Fournon I., Afonso-Luis A., 2013, MNRAS, 428, 291

Pérez-González P. G., 2003, PASP, 115, 1353

Pilbratt G. L. et al., 2010, A&A, 518, L1

Planck Collaboration VII, 2011, A&A, 536, A7

Poglitsch A. et al., 2010, A&A, 518, L2

Polletta M. et al., 2007, ApJ, 663, 81

Pozzi F. et al., 2004, ApJ, 609, 122

Rodighiero G. et al., 2010, A&A, 515, A8

Roseboom I. G. et al., 2010, MNRAS, 409, 48

Roseboom I. G. et al., 2012, MNRAS, 419, 2758

Rowan-Robinson M., Gonzalez-Solares E., Vaccari M., Marchetti L., 2013, MNRAS, 428, 1958

Rush B., Malkan M. A., 1993, BAAS, 25, 1362

Sandage A., Tammann G. A., Yahil A., 1979, ApJ, 232, 352

Sanders D. B., Mazzarella J. M., Kim D.-C., Surace J. A., Soifer B. T., 2003, AJ, 126, 1607

Saunders W., Rowan-Robinson M., Lawrence A., Efstathiou G., Kaiser N., Ellis R. S., Frenk C. S., 1990, MNRAS, 242, 318

Schafer C. M., 2007, ApJ, 661, 703

Schechter P., 1976, ApJ, 203, 297

Schmidt M., 1968, ApJ, 151, 393

Serjeant S., Harrison D., 2005, MNRAS, 356, 192

Serjeant S., Gruppioni C., Oliver S., 2002, MNRAS, 330, 621

Shupe D. L., Fang F., Hacking P. B., Huchra J. P., 1998, ApJ, 501, 597

Smith A. J. et al., 2012, MNRAS, 419, 377

Sullivan M., Treyer M. A., Ellis R. S., Bridges T. J., Milliard B., Donas J., 2000, MNRAS, 312, 442

Symeonidis M. et al., 2013, MNRAS, 431, 2317

Takeuchi T. T., Ishii T. T., Dole H., Dennefeld M., Lagache G., Puget J.-L., 2006, A&A, 448, 525

Taylor M. B., 2005, in Shopbell P., Britton M., Ebert R., eds, ASP Conf. Ser. Vol. 347, Astronomical Data Analysis Software and Systems XIV. Astron. Soc. Pac., San Francisco, p. 29

Tresse L., Maddox S. J., 1998, ApJ, 495, 691

Vaccari M. et al., 2010, A&A, 518, L20

Valiante E., Lutz D., Sturm E., Genzel R., Chapin E. L., 2009, ApJ, 701, 1814

Wang L., Rowan-Robinson M., 2009, MNRAS, 398, 109

Wang L. et al., 2014, MNRAS, 444, 2870

Werner M. W. et al., 2004, ApJ, 154, 1

Westra E., Geller M. J., Kurtz M. J., Fabricant D. G., Dell'Antonio I., 2010, ApJ, 708, 534

Wyder T. K. et al., 2005, ApJ, 619, L15

APPENDIX A

<table><caption>**Table A1.** SPIRE 250, 350, 500 µm and IR bolometric rest-frame $1/V_{max}$ LF estimates in the redshift ranges between 0.02 and 0.5, using the HerMES Wide Fields sample. L indicates $v L_v$ for the monochromatic LFs and $L_{IR}$ indicates the integrated luminosity between 8 and 1000 µm for the IR bolometric rest-frame LF. These L are expressed in units of $L_{\odot}$ and LLF estimates and their errors are in [Mpc$^{-3}$ dex$^{-1}$]. The quantity $\sigma$ is the total error (Poissonian error + redshift uncertainties, estimated as explained in the text) associated with $\Phi$ in each band and luminosity/redshift bin.</caption><thead><tr><th>log L</th><th>log (Φ, σ)<sub>250</sub></th><th>log (Φ, σ)<sub>350</sub></th><th>log (Φ, σ)<sub>500</sub></th><th>log (Φ, σ)<sub>IR</sub></th></tr><tr><td></td><td colspan="4" style="text-align: center;">0.02 &lt; z &lt; 0.1 LFs</td></tr></thead><tbody><tr><td>7.16</td><td>-</td><td>-</td><td>-2.19, -2.07</td><td>-</td></tr><tr><td>7.33</td><td>-</td><td>-</td><td>-1.90, -2.54</td><td>-</td></tr><tr><td>7.49</td><td>-</td><td>-</td><td>-2.06, -2.72</td><td>-</td></tr><tr><td>7.66</td><td>-</td><td>-</td><td>-2.06, -2.90</td><td>-</td></tr><tr><td>7.83</td><td>-</td><td>-</td><td>-2.17, -3.11</td><td>-</td></tr><tr><td>8.00</td><td>-</td><td>-</td><td>-2.12, -3.21</td><td>-</td></tr><tr><td>8.16</td><td>-2.10, -2.06</td><td>-2.09, -2.89</td><td>-2.15, -3.34</td><td>-</td></tr><tr><td>8.33</td><td>-1.96, -2.57</td><td>-2.15, -3.05</td><td>-2.30, -3.47</td><td>-</td></tr><tr><td>8.49</td><td>-2.08, -2.72</td><td>-2.11, -3.17</td><td>-2.52, -3.58</td><td>-</td></tr><tr><td>8.66</td><td>-2.06, -2.89</td><td>-2.14, -3.31</td><td>-2.75, -3.69</td><td>-</td></tr><tr><td>8.83</td><td>-2.17, -3.08</td><td>-2.23, -3.43</td><td>-3.34, -3.99</td><td>-</td></tr><tr><td>9.00</td><td>-2.15, -3.20</td><td>-2.48, -3.56</td><td>-3.44, -4.04</td><td>-2.10, -2.60</td></tr><tr><td>9.16</td><td>-2.11, -3.31</td><td>-2.66, -3.65</td><td>-4.34, -4.49</td><td>-2.11, -2.69</td></tr><tr><td>9.33</td><td>-2.26, -3.45</td><td>-3.05, -3.84</td><td>-</td><td>-1.96, -2.69</td></tr><tr><td>9.49</td><td>-2.49, -3.57</td><td>-3.49, -4.07</td><td>-</td><td>-2.17, -3.00</td></tr></tbody></table>
<table><caption>Table A1 – continued.</caption><thead><tr><th>log L</th><th>log (Φ, σ)<sub>250</sub></th><th>log (Φ, σ)<sub>350</sub></th><th>log (Φ, σ)<sub>500</sub></th><th>log (Φ, σ)<sub>IR</sub></th></tr></thead><tbody><tr><td>9.66</td><td>-2.69, -3.67</td><td>-3.86, -4.25</td><td>-</td><td>-2.12, -2.99</td></tr><tr><td>9.83</td><td>-3.12, -3.88</td><td>-</td><td>-</td><td>-2.13, -3.24</td></tr><tr><td>10.00</td><td>-3.53, -4.08</td><td>-</td><td>-</td><td>-2.25, -3.35</td></tr><tr><td>10.16</td><td>-4.04, -4.34</td><td>-</td><td>-</td><td>-2.35, -3.47</td></tr><tr><td>10.33</td><td>-</td><td>-</td><td>-</td><td>-2.50, -3.56</td></tr><tr><td>10.49</td><td>-</td><td>-</td><td>-</td><td>-2.91, -3.77</td></tr><tr><td>10.66</td><td>-</td><td>-</td><td>-</td><td>-3.06, -3.85</td></tr><tr><td>10.83</td><td>-</td><td>-</td><td>-</td><td>-3.28, -3.96</td></tr><tr><td>11.00</td><td>-</td><td>-</td><td>-</td><td>-3.94, -4.29</td></tr><tr><td>11.16</td><td>-</td><td>-</td><td>-</td><td>-4.16, -4.40</td></tr><tr><td colspan="5" style="text-align: center;">0.1 &lt; z &lt; 0.2 LFs</td></tr><tr><td>8.33</td><td>-</td><td>-</td><td>-2.31, -3.23</td><td>-</td></tr><tr><td>8.49</td><td>-</td><td>-</td><td>-2.31, -3.52</td><td>-</td></tr><tr><td>8.66</td><td>-</td><td>-</td><td>-2.44, -3.76</td><td>-</td></tr><tr><td>8.83</td><td>-</td><td>-</td><td>-2.70, -4.05</td><td>-</td></tr><tr><td>9.00</td><td>-</td><td>-</td><td>-3.08, -4.27</td><td>-</td></tr><tr><td>9.16</td><td>-</td><td>-</td><td>-3.57, -4.51</td><td>-</td></tr><tr><td>9.33</td><td>-2.28, -3.14</td><td>-2.63, -4.01</td><td>-4.21, -4.82</td><td>-</td></tr><tr><td>9.49</td><td>-2.28, -3.49</td><td>-2.91, -4.18</td><td>-4.76, -5.10</td><td>-</td></tr><tr><td>9.66</td><td>-2.41, -3.76</td><td>-3.39, -4.42</td><td>-</td><td>-</td></tr><tr><td>9.83</td><td>-2.64, -3.95</td><td>-4.02, -4.73</td><td>-</td><td>-</td></tr><tr><td>10.00</td><td>-2.98, -4.22</td><td>-4.55, -5.00</td><td>-</td><td>-</td></tr><tr><td>10.16</td><td>-3.45, -4.45</td><td>-5.45, -5.44</td><td>-</td><td>-2.36, -3.37</td></tr><tr><td>10.33</td><td>-4.05, -4.74</td><td>-</td><td>-</td><td>-2.41, -3.66</td></tr><tr><td>10.49</td><td>-4.61, -5.03</td><td>-</td><td>-</td><td>-2.48, -3.64</td></tr><tr><td>10.66</td><td>-</td><td>-</td><td>-</td><td>-2.69, -3.92</td></tr><tr><td>10.83</td><td>-</td><td>-</td><td>-</td><td>-3.01, -4.18</td></tr><tr><td>11.00</td><td>-</td><td>-</td><td>-</td><td>-3.40, -4.41</td></tr><tr><td>11.16</td><td>-</td><td>-</td><td>-</td><td>-3.64, -4.16</td></tr><tr><td>11.33</td><td>-</td><td>-</td><td>-</td><td>-4.15, -4.80</td></tr><tr><td>11.49</td><td>-</td><td>-</td><td>-</td><td>-4.76, -5.10</td></tr><tr><td colspan="5" style="text-align: center;">0.2 &lt; z &lt; 0.3 LFs</td></tr><tr><td>8.83</td><td>-</td><td>-</td><td>-2.72, -3.57</td><td>-</td></tr><tr><td>9.00</td><td>-</td><td>-</td><td>-2.88, -4.14</td><td>-</td></tr><tr><td>9.16</td><td>-</td><td>-</td><td>-3.21, -4.13</td><td>-</td></tr><tr><td>9.33</td><td>-</td><td>-</td><td>-3.69, -4.75</td><td>-</td></tr><tr><td>9.49</td><td>-</td><td>-2.78, -4.02</td><td>-4.19, -4.93</td><td>-</td></tr><tr><td>9.66</td><td>-</td><td>-3.09, -4.11</td><td>-5.25, -5.54</td><td>-</td></tr><tr><td>9.83</td><td>-</td><td>-3.52, -4.66</td><td>-5.85, -5.82</td><td>-</td></tr><tr><td>10.00</td><td>-2.81, -4.05</td><td>-4.02, -4.88</td><td>-</td><td>-</td></tr><tr><td>10.16</td><td>-3.13, -4.21</td><td>-4.85, -5.33</td><td>-</td><td>-</td></tr><tr><td>10.33</td><td>-3.56, -4.65</td><td>-5.55, -5.69</td><td>-</td><td>-</td></tr><tr><td>10.49</td><td>-4.13, -4.92</td><td>-</td><td>-</td><td>-</td></tr><tr><td>10.66</td><td>-4.89, -5.37</td><td>-</td><td>-</td><td>-2.81, -3.62</td></tr><tr><td>10.83</td><td>-5.55, -5.69</td><td>-</td><td>-</td><td>-2.97, -4.06</td></tr><tr><td>11.00</td><td>-</td><td>-</td><td>-</td><td>-3.208, -4.19</td></tr><tr><td>11.16</td><td>-</td><td>-</td><td>-</td><td>-3.48, -4.52</td></tr><tr><td>11.33</td><td>-</td><td>-</td><td>-</td><td>-3.77, -4.77</td></tr><tr><td>11.49</td><td>-</td><td>-</td><td>-</td><td>-4.27, -4.85</td></tr><tr><td>11.66</td><td>-</td><td>-</td><td>-</td><td>-4.85, -5.25</td></tr><tr><td colspan="5" style="text-align: center;">0.3 &lt; z &lt; 0.4 LFs</td></tr><tr><td>9.33</td><td>-</td><td>-</td><td>-3.10, -3.92</td><td>-</td></tr><tr><td>9.49</td><td>-</td><td>-</td><td>-3.58, -4.44</td><td>-</td></tr><tr><td>9.66</td><td>-</td><td>-</td><td>-4.31, -4.75</td><td>-</td></tr><tr><td>9.83</td><td>-</td><td>-3.03, -3.92</td><td>-4.88, -5.33</td><td>-</td></tr><tr><td>10.00</td><td>-</td><td>-3.43, -4.40</td><td>-6.09, -6.05</td><td>-</td></tr><tr><td>10.16</td><td>-</td><td>-4.04, -4.74</td><td>-5.79, -5.94</td><td>-</td></tr><tr><td>10.33</td><td>-3.02, -3.87</td><td>-4.72, -5.29</td><td>-</td><td>-</td></tr><tr><td>10.49</td><td>-3.47, -4.40</td><td>-5.62, -5.84</td><td>-</td><td>-</td></tr><tr><td>10.66</td><td>-4.13, -4.74</td><td>-5.79, -5.94</td><td>-</td><td>-</td></tr><tr><td>10.83</td><td>-4.80, -5.31</td><td>-</td><td>-</td><td>-</td></tr><tr><td>11.00</td><td>-5.79, -5.92</td><td>-</td><td>-</td><td>-3.06, -3.96</td></tr></tbody></table>
<table><caption>Table A1 – continued.</caption><thead><tr><th>log L</th><th>log (Φ, σ)<sub>250</sub></th><th>log (Φ, σ)<sub>350</sub></th><th>log (Φ, σ)<sub>500</sub></th><th>log (Φ, σ)<sub>IR</sub></th></tr></thead><tbody><tr><td>11.16</td><td>–</td><td>–</td><td>–</td><td>-3.25, -4.05</td></tr><tr><td>11.33</td><td>–</td><td>–</td><td>–</td><td>-3.42, -4.49</td></tr><tr><td>11.49</td><td>–</td><td>–</td><td>–</td><td>-3.73, -4.81</td></tr><tr><td>11.66</td><td>–</td><td>–</td><td>–</td><td>-4.16, -4.10</td></tr><tr><td>11.83</td><td>–</td><td>–</td><td>–</td><td>-4.63, -5.36</td></tr><tr><td colspan="5">0.4 &lt; z &lt; 0.5 LFs</td></tr><tr><td>9.66</td><td>–</td><td>–</td><td>-3.75, -4.43</td><td>–</td></tr><tr><td>9.83</td><td>–</td><td>-3.03, -3.92</td><td>-4.63, -4.71</td><td>–</td></tr><tr><td>10.00</td><td>–</td><td>-3.43, -4.40</td><td>-4.99, -5.38</td><td>–</td></tr><tr><td>10.16</td><td>–</td><td>-4.04, -4.74</td><td>-5.68, -5.89</td><td>–</td></tr><tr><td>10.33</td><td>–</td><td>-4.35, -4.70</td><td>–</td><td>–</td></tr><tr><td>10.49</td><td>-3.23, -4.11</td><td>-4.92, -5.37</td><td>–</td><td>–</td></tr><tr><td>10.66</td><td>-3.69, -4.41</td><td>-5.35, -5.76</td><td>–</td><td>–</td></tr><tr><td>10.83</td><td>-4.41, -4.70</td><td>–</td><td>–</td><td>–</td></tr><tr><td>11.00</td><td>-4.97, -5.38</td><td>–</td><td>–</td><td>–</td></tr><tr><td>11.16</td><td>-5.58, -5.85</td><td>–</td><td>–</td><td>–</td></tr><tr><td>11.33</td><td>–</td><td>–</td><td>–</td><td>-3.45, -4.42</td></tr><tr><td>11.49</td><td>–</td><td>–</td><td>–</td><td>-3.61, -4.42</td></tr><tr><td>11.66</td><td>–</td><td>–</td><td>–</td><td>-3.87, -4.65</td></tr><tr><td>11.83</td><td>–</td><td>–</td><td>–</td><td>-4.23, -5.10</td></tr><tr><td>12.00</td><td>–</td><td>–</td><td>–</td><td>-4.99, -5.60</td></tr></tbody></table>

This paper has been typeset from a TeX/LATEX file prepared by the author.

<footnote n="1">¹Department of Physical Sciences, The Open University, Milton Keynes MK7 6AA, UK</footnote>

<footnote n="2">²Dipartimento di Fisica e Astronomia, Università di Padova, vicolo Osservatorio 3, I-35122 Padova, Italy</footnote>

<footnote n="3">³Department of Physics and Astronomy, University of the Western Cape, Robert Sobukwe Road, 7535 Bellville, Cape Town, South Africa</footnote>

<footnote n="4">⁴INAF – Istituto di Radioastronomia, via Gobetti 101, I-40129 Bologna, Italy</footnote>

<footnote n="5">⁵Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="6">⁶Laboratoire AIM-Paris-Saclay, CEA/DSM/Irfu – CNRS – Université Paris Diderot, CE-Saclay, pt courrier 131, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="7">⁷Institut d'Astrophysique Spatiale (IAS), bâtiment 121, Université Paris-Sud 11 and CNRS (UMR 8617), F-91405 Orsay, France</footnote>

<footnote n="8">⁸ESO, Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany</footnote>

<footnote n="9">⁹California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA</footnote>

<footnote n="10">¹⁰Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="11">¹¹Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388 Marseille, France</footnote>

<footnote n="12">¹²Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK</footnote>

<footnote n="13">¹³Center for Astrophysics and Space Astronomy 389-UCB, University of Colorado, Boulder, CO 80309, USA</footnote>

<footnote n="14">¹⁴Herschel Science Centre, European Space Astronomy Centre, Villanueva de la Cañada, E-28691 Madrid, Spain</footnote>

<footnote n="15">¹⁵Department of Physics & Astronomy, University of California, Irvine, CA 92697, USA</footnote>

<footnote n="16">¹⁶Department of Physics, Virginia Tech, Blacksburg, VA 24061, USA</footnote>

<footnote n="17">¹⁷Sorbonne Universités, UPMC-CNRS, UMR7095, Institut d'Astrophysique de Paris, F-75014 Paris, France</footnote>

<footnote n="18">¹⁸Department of Astrophysical and Planetary Sciences, CASA 389-UCB, University of Colorado, Boulder, CO 80309, USA</footnote>

<footnote n="19">¹⁹School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff CF24 3AA, UK</footnote>

<footnote n="20">²⁰Instituto de Física y Astronomía, Universidad de Valparaíso, Avda. Gran Bretaña 1111, Valparaíso, Chile</footnote>

<footnote n="21">²¹UK Astronomy Technology Centre, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="22">²²Astronomy Centre, Department of Physics & Astronomy, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote n="23">²³Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey RH5 6NT, UK</footnote>

<footnote n="24">²⁴RAL Space, Rutherford Appleton Laboratory, Chilton, Didcot, Oxfordshire OX11 0QX, UK</footnote>

<footnote n="25">²⁵Instituto de Astrofísica de Canarias (IAC), E-38200 La Laguna, Tenerife, Spain</footnote>

<footnote n="26">²⁶Departamento de Astrofísica, Universidad de La Laguna (ULL), E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="27">²⁷Department of Physics, Denys Wilkinson Building, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="28">²⁸Infrared Processing and Analysis Center, MS 100-22, California Institute of Technology, JPL, Pasadena, CA 91125, USA</footnote>

<footnote n="29">²⁹Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada</footnote>

<footnote n="30">³⁰International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Perth, Australia</footnote>

<footnote n="31">³¹Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="32">³²SRON Netherlands Institute for Space Research, Landleven 12, NL-9747 AD Groningen, the Netherlands</footnote>

<footnote n="33">³³Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen, Denmark</footnote>