# The Herschel Virgo Cluster Survey XVII. SPIRE point-source catalogs and number counts

Ciro Pappalardo, George J. Bendo, Simone Bianchi, Leslie Hunt, Stefano Zibetti, Edvige Corbelli, Sperello Di Serego Alighieri, Marco Grossi, Jonathan Davies, Maarten Baes, et al.

► To cite this version:

Ciro Pappalardo, George J. Bendo, Simone Bianchi, Leslie Hunt, Stefano Zibetti, et al. The Her-
schel Virgo Cluster Survey XVII. SPIRE point-source catalogs and number counts. Astronomy and
Astrophysics - A&A, 2015, 573, 10.1051/0004-6361/201424549. hal-01439923

HAL Id: hal-01439923

https://hal.science/hal-01439923

Submitted on 15 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
# The Herschel* Virgo Cluster Survey
## XVII. SPIRE point-source catalogs and number counts**

Ciro Pappalardo<sup>1,2</sup>, George J. Bendo<sup>3</sup>, Simone Bianchi<sup>4</sup>, Leslie Hunt<sup>4</sup>, Stefano Zibetti<sup>4</sup>, Edvige Corbelli<sup>4</sup>, Sperello di Serego Alighieri<sup>4</sup>, Marco Grossi<sup>1,2</sup>, Jonathan Davies<sup>5</sup>, Maarten Baes<sup>6</sup>, Ilse De Looze<sup>6</sup>, Jacopo Fritz<sup>6</sup>, Michael Pohlen<sup>5</sup>, Matthew W. L. Smith<sup>5</sup>, Joris Verstappen<sup>6</sup>, Médéric Boquien<sup>7</sup>, Alessandro Boselli<sup>7</sup>, Luca Cortese<sup>8,9</sup>, Thomas Hughes<sup>6</sup>, Sebastien Viaene<sup>6</sup>, Luca Bizzocchi<sup>1,2</sup>, and Marcel Clemens<sup>10</sup>

Received 7 July 2014 / Accepted 8 September 2014

## ABSTRACT

**Aims.** We present three independent catalogs of point-sources extracted from SPIRE images at 250, 350, and 500 µm, acquired with the Herschel Space Observatory as a part of the Herschel Virgo Cluster Survey (HeViCS). The catalogs have been cross-correlated to consistently extract the photometry at SPIRE wavelengths for each object.

**Methods.** Sources have been detected using an iterative loop. The source positions are determined by estimating the likelihood to be a real source for each peak on the maps, according to the criterion defined in the *sourceExtractorSussextractor* task. The flux densities are estimated using the *sourceExtractorTimeline*, a timeline-based point source fitter that also determines the fitting procedure with the width of the Gaussian that best reproduces the source considered. Afterwards, each source is subtracted from the maps, removing a Gaussian function in every position with the full width half maximum equal to that estimated in *sourceExtractorTimeline*. This procedure improves the robustness of our algorithm in terms of source identification. We calculate the completeness and the flux accuracy by injecting artificial sources in the timeline and estimate the reliability of the catalog using a permutation method.

**Results.** The HeViCS catalogs contain about 52 000, 42 200, and 18 700 sources selected at 250, 350, and 500 µm above $3\sigma$ and are ~75%, 62%, and 50% complete at flux densities of 20 mJy at 250, 350, 500 µm, respectively. We then measured source number counts at 250, 350, and 500 µm and compare them with previous data and semi-analytical models. We also cross-correlated the catalogs with the Sloan Digital Sky Survey to investigate the redshift distribution of the nearby sources. From this cross-correlation, we select ~2000 sources with reliable fluxes and a high signal-to-noise ratio, finding an average redshift $z \sim 0.3 \pm 0.22$ and 0.25 (16–84 percentile).

**Conclusions.** The number counts at 250, 350, and 500 µm show an increase in the slope below 200 mJy, indicating a strong evolution in number of density for galaxies at these fluxes. In general, models tend to overpredict the counts at brighter flux densities, underlying the importance of studying the Rayleigh-Jeans part of the spectral energy distribution to refine the theoretical recipes of the models. Our iterative method for source identification allowed the detection of a family of 500 µm sources that are not foreground objects belonging to Virgo and not found in other catalogs.

**Key words.** galaxies: evolution – galaxies: statistics – infrared: galaxies – submillimeter: galaxies

## 1. Introduction

At least half of the energy emitted by the stars and active galactic nuclei (AGN) in a Hubble time has been absorbed by dust in the interstellar medium and re-emitted as IR light (Puget et al. 1996; Fixsen et al. 1998; Lagache et al. 1999). The abundance and properties of dust are then crucial to determine the cosmic

<footnote n="1">¹ Centro de Astronomia e Astrofísica da Universidade de Lisboa, Observatório Astronómico de Lisboa, Tapada da Ajuda, 1349-018 Lisbon, Portugal
e-mail: ciro@oal.ul.pt</footnote>

<footnote n="2">² Department of Physics, Faculty of Sciences, University of Lisbon, Campo Grande, 1749-016 Lisbon, Portugal</footnote>

<footnote n="3">³ Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="4">⁴ Osservatorio Astrofisico di Arcetri – INAF, Largo E. Fermi 5, 50125 Firenze, Italy
e-mail: cpappala@arcetri.astro.it</footnote>

<footnote n="5">⁵ Department of Physics and Astronomy, Cardiff University, The Parade, Cardiff, CF24 3AA, UK</footnote>

<footnote n="6">⁶ Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, 9000 Gent, Belgium</footnote>

<footnote n="7">⁷ Laboratoire d'Astrophysique de Marseille - LAM, Université d'Aix-Marseille & CNRS, UMR 7326, 38 rue F. Joliot-Curie, 13388 Marseille Cedex 13, France</footnote>

<footnote n="8">⁸ Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Mail H30 - PO Box 218, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="9">⁹ European Southern Observatory, Karl-Schwarzschild Str. 2, 85748 Garching bei Muenchen, Germany</footnote>

<footnote n="10">¹⁰ Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy</footnote>

<footnote>* Herschel is an ESA space observatory with science instruments provided by a European-led principal investigator consortia and with an important participation from NASA.
** The 250, 350, 500 µm, and the total catalogs are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/573/A129</footnote>
comoving star-formation rate density (SFRD), because the opti-
cal/ultraviolet emission originating from stars must be corrected
for IR extinction. Moreover, if the bulk of far-IR (FIR) lumi-
nosity is due to the dust heated by massive stars in star-forming
regions, the IR luminosity could be used as direct star formation
tracer. However, at $\lambda > 250$ µm, much of the dust emission in
nearby galaxies may be dominated by dust heated by the total
stellar population, mainly formed by evolved stars (e.g., Bendo
et al. 2010, 2012; Boqujen et al. 2011).

Previous observations with IRAS (Saunders et al. 1990; Oliver et al. 1992), ISO (Oliver et al. 2002), and Spitzer (Shupe et al. 2008; Frayer et al. 2009) have shown a strong evolution of galaxy luminosities in mid-IR and FIR wavelengths because of the decline of the average SFRD with time. At redshift $z \sim 1$, luminous IR galaxies contribute about 50% to the SFRD (Murphy et al. 2011), while the SFRD increases at higher redshift, $z \sim 2$ because we have an important contribution from more active star-forming galaxies, or ultraluminous IR galaxies ($L_{IR} \ge 10^{12} L_{\odot}$, Chapman et al. 2005; Daddi et al. 2007; Papovich et al. 2007; Magdis et al. 2012; Magnelli et al. 2009), which are considered to be the ancestors of massive quiescent galaxies in the present epoch (e.g. Le Floc'h et al. 2005).

To determine the contribution of IR light from less luminous galaxies to the SFRD, deep IR observations are needed. Previous works have shown that confusion is the most important component of noise at FIR wavelengths (Blain et al. 1998). Faint unresolved sources are collected inside the observing beam, and the sky background in the resulting map is produced mostly by these objects. The bulk of this unresolved emission consists of dusty star-forming galaxies at a redshift $z \ge 1.2$ (Devlin et al. 2009) and forms the FIR background (FIRB, Puget et al. 1996; Fixsen et al. 1998; Lagache et al. 1999).

The Herschel Space Observatory (Pilbratt et al. 2010) with its two spectro-photometers, SPIRE (Griffin et al. 2010) and PACS (Poglitsch et al. 2010), has the necessary sensitivity to investigate the FIRB in greater detail, directly measuring the IR luminosity function at higher redshift. The starting point for such an investigation is the construction of a list of unresolved sources down to the confusion noise limit. These data can be used to infer statistical properties of galaxy populations to determine their evolution. A tool widely used to investigate galaxy evolution over Hubble time is the number density of galaxies as a function of flux density, or the number counts. They are a direct tracer of the past galaxy evolution that can be compared to theoretical models.

In this paper, we present three catalogs of point-sources extracted from the Herschel Virgo Cluster Survey (HeViCS, Davies et al. 2010), a survey that observed a region that is about 84 square degrees centered in the Virgo cluster, of which 55 square degrees are observed at uniform depth. In this survey, Herschel observed with SPIRE and PACS four overlapping fields in fast parallel-mode, labeled as V1, V2, V3, and V4 in Fig. 1, which are 16 square degrees each (see Auld et al. 2013, for details about the survey and the data reduction). In this paper we considered only SPIRE data in the analysis, leaving the PACS catalog and data analysis to a future companion paper (Fritz et al., in prep.).

Although the primary goal of HeViCS was to study the galaxies within the Virgo Cluster itself, the survey serendipitously detected a substantial number of background galaxies, making it an excellent survey for studying the FIRB.

The paper is organized as follows: in Sect. 2, we describe the method for the determination of the source flux and position and show how we estimate the background. We also investigate the

goodness of our catalog through the reliability, the completeness,
and the accuracy in flux estimation. In Sect. 3, we discuss the
number counts of the HeViCS catalogs with an analysis of the
redshift distribution of the sources, which are obtained by cross-
correlating the catalog with Sloan Digital Sky Survey catalog
(SDSS - DR10, Ahn et al. 2014). Our conclusions are given in
Sect. 4.

**2. Methodology**

The determination of the position and the flux density of the
sources have been performed in different steps. From the start-
ing map, presented in Sect. 2.1, an estimate of the back-
ground has been subtracted (Sect. 2.2). The source positions
have been determined from this background-subtracted map us-
ing the *sourceExtractorSussextractor* software. Finally,
the flux density in each position was estimated using the
*sourceExtractorTimeline* software, and each source has
been subtracted from the timeline data once determined its flux
density (Sect. 2.3).

<figure><img src="image_2.png" /><figcaption>Fig. 1. HeViCS image at 250 µm with detected sources (red points). The blue squares indicate regions with the full coverage.</figcaption></figure>
## 2.1. Data

The four HeViCS fields have been observed in two orthogonal scans at a sampling rate of 60'' s⁻<fnref n="1" /> with four repetitions for a total of eight scans. The final map is obtained by averaging the bolometer timeline measurements that sampled the region in each map pixel (naiveScan mapping task in HIPE 9.0.0, Ott 2010), as described in detail in Auld et al. (2013). The main difference between our pipeline and the standard one is that we did not run the default temperatureDrift-Correction task and the residual, median baseline subtraction. Instead we use the BRIght Galaxy ADaptive Element (BriGAdE, Smith et al., in prep.) which uses an alternative technique for correcting temperature drifts. The SPIRE instrument (Griffin et al. 2010)<fnref n="1" /> has three bolometer arrays that observed the sky at 250, 350, and 500 µm. The SPIRE calibration tree v8 was used and the full width at half-maximum (FWHM) measurements of the beams in the final maps for point sources are 17.''5, 23.''9 and 35.''1 at 250, 350, and 500 µm, respectively<fnref n="2" />. The maps we used have pixel sizes of 6'', 8'', and 12'' at 250, 350, and 500 µm, respectively. The best spatial resolution and the deepest data are at 250 µm, and for this reason, we chose the source positions at this wavelength as a reference for HeViCS catalogs.

## 2.2. Background subtraction

The SPIRE data are contaminated by Galactic cirrus emission, particularly in the V2 field. The cirrus emission does not significantly affect the flux densities measured using the timeline-based source fitter, but it does cause `sourceExtractorSussextractor` to identify more sources in regions with strong cirrus emission. This is because the diffuse component in the convolved map is considered by `sourceExtractorSussextractor` as a sequence of peaks with each one a high likelihood to be a source (see Sect. 2.3 for details). To mitigate this effect, we estimate the background using a method developed in `Sextractor` (Bertin & Arnouts 1996). In this method, the map is regridded into cells larger than the pixel size, and then we iteratively estimate the mean and the standard deviation of the distribution of pixel values in each cell. At each step, the outlier pixels are removed until the variation of $\sigma$ in an iteration is small (see Bertin & Arnouts 1996, for details). The value of the background in each block is then estimated by using $2.5 \times M_{\text{med}} - 1.5 \times M_{\text{avg}}$, where $M_{\text{med}}$ and $M_{\text{avg}}$ are the median and the mean inside each block. Finally, the background map is obtained by linearly interpolating between the cells. The size of the blocks fixes the size of the structure that we are able to remove. If the cell size is too large, it results in a map in which the diffuse emission is not well constrained, while a cell size that is too small results in a map that reproduces the diffuse emission well but does not remove the extended emission due to the resolved extragalactic sources larger than the beam. After performing a series of tests, we chose to fix the block size to 8 pixels, slightly smaller than the beam size (9 pixels). Figure 2 shows a region of the V2 field before and after cirrus subtraction. This background is subtracted from the map data that is used as input by `sourceExtractorSussextractor` in the source identification step described in Sect. 2.3.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** V2 field region before (left panel) and after (right panel) background subtraction.</figcaption></figure>

To test the reliability of this background estimation, we selected 45 boxes of 0.5 × 0.5 deg<fnref n="2" /> that cover different regions in the V2 field with strong and weak cirrus emission (top panel of Fig. 3). In each box, we estimate the number of sources detected before and after the background subtraction. We fix a low threshold for the source-detection (~5 mJy) because we want to also consider the sources with flux densities well below the limit chosen for the catalog for this test. Selecting a higher threshold would introduce a bias due to sources with flux densities around the detection limit. If the background subtraction has worked properly, we should measure the same number of sources in each box with relatively small variations in the total number of sources detected in each box. Results are shown in bottom panel of Fig. 3. Before the subtraction, the number of sources detected with `sourceExtractorSussextractor` spans a high range of values, $188 \le N_{\text{bef}} \le 1199$ with an average value of 600 ± 242.7 (bottom panel of Fig. 3). After the background subtraction the number of sources is quite homogeneous for regions with either strong or weak cirrus contamination with an average value of 697 sources and a significantly lower standard deviation of 37.9 (blue dashed line in bottom panel of Fig. 3). The number of detections above 20 mJy are consistent with the results obtained with the previous test. They reproduce a similar trend but are rescaled, because of the smaller number of sources due to the higher threshold. The background-subtraction process reduces the standard deviation of the detections above 20 mJy from 50% (before subtraction) to 10% (after the subtraction). This reduction in the standard deviation and the homogeneity of the number of detections after the background subtraction confirm the effectiveness of our cirrus removal methods.

## 2.3. Source identification and flux density measurement

Source detection is performed in an iterative loop. A list of initial positions of point sources is produced using the implementation of `sourceExtractorSussextractor` (Savage & Oliver 2007) in HIPE 9.0.0 (Ott 2010) that is applied to a map where the foreground cirrus structure has been subtracted (see Sect. 2.2). The `sourceExtractorSussextractor` software convolves the reduced map with a bidimensional Gaussian kernel having a FWHM of 17.''5, 23.''9, and 35.''1 at 250, 350, and 500 µm, respectively. The value of each pixel in the convolved map corresponds to the maximum likelihood estimate of the flux density that a source should have if centered on that pixel. The extractor routine finally searches for local maxima in pixel regions with a size comparable to the beam, which are the neighboring eight pixels. For this analysis, we only used sources from locations for which we have a redundancy of ≥35 bolometer samplings

<footnote n="1">¹ http://herschel.esac.esa.int/Docs/SPIRE/pdf/spire.om.pdf</footnote>

<footnote n="2">² http://herschel.esac.esa.int/twiki/pub/Public/SpireCalibrationWeb/beam_release note v1-1.pdf and http://herschel.esac.esa.int/twiki/bin/view/Public/SpirePhotometerBeamProfileAnalysis</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Top panel: boxes used in V2 field to test cirrus subtraction. Bottom panel: distribution of the number of sources detected in the V2 boxes shown in the top panel before (black solid line) and after (blue dashed line) background subtraction.</figcaption></figure>

per pixel, or eight coverages (see Table 1 and Fig. 1). Since the main source of noise in our maps is the confusion given by the FIRB (Puget et al. 1996), we only use sources with initial flux density measurements of >20 mJy. This threshold corresponds approximately to $3\sigma$ above the confusion noise (estimated in ~6–7 mJy beam$^{-1}$ at 250, 350, and 500 µm, Nguyen et al. 2010).

Each source is then fitted with a Gaussian function using `sourceExtractorTimeline`, a timeline-based point source fitter (see Bendo et al. 2013) implemented in HIPE 9.0.0 (Ott 2010). The method fits timeline data from all bolometers within an individual array with a two-dimensional Gaussian

<table><caption>**Table 1.** Regions fully sampled, which are selected in each HeViCS field for the source extraction with flux densities >20 mJy.</caption><thead><tr><th>Field<br>(1)</th><th>RA<br>deg<br>(2)</th><th>Dec<br>deg<br>(3)</th><th>Pos. angle<br>deg<br>(4)</th><th>Area<br>deg<sup>2</sup><br>(5)</th><th>250 µm<br>sources<br>(6)</th></tr></thead><tbody><tr><td>V1</td><td>184.3</td><td>14.15</td><td>20</td><td>4 × 4</td><td>14 061</td></tr><tr><td>V2</td><td>187.2</td><td>12.07</td><td>20</td><td>4 × 4</td><td>12 956</td></tr><tr><td>V3</td><td>187.1</td><td>8.6</td><td>25</td><td>4 × 4</td><td>13 828</td></tr><tr><td>V4</td><td>185.0</td><td>5.25</td><td>20</td><td>4 × 4</td><td>14 494</td></tr><tr><td>Combined</td><td>—</td><td>—</td><td>—</td><td>55</td><td>52 020</td></tr></tbody></table>

**Notes.** Column 1: field; Columns 2, 3: center coordinates; Column 4: position angle (counterclockwise); Column 5: area covered; Column 6: number of sources detected.

function. As an input, the fitter uses the source position given by `sourceExtractorSussextractor`, an aperture containing the target data (using apertures of 22", 30", and 42" at 250, 350, and 500 µm, respectively), and an aperture identifying the background annulus (set to a radius of 700–720"). The output consists of parameters describing the amplitude, position, the FWHM, and the parameters describing the background flux density.

Two fits are performed to the timeline data. We first fit a circular Gaussian function in which the FWHM is allowed to vary. We then use the FWHM to determine whether the source is resolved or unresolved and reject sources that are either too narrow (which may be unremoved glitches) or too broad (which are probably extended sources). To test at which distance two sources are distinguishable, we injected in the timeline data point-like objects with the same flux density and an increasing distance from each other. We found that two sources are resolved only if they are at a distance above ~22", 30" and 46" at 250, 350, and 500 µm, respectively. Fixing the distance between the sources to these values, we then injected in the timeline data point-like objects at different flux densities to investigate the FWHM recovered by `sourceExtractorTimeline`. Based on these tests, we have found that artificial sources with FWHM corresponding to the telescope beam (17.5, 23.9, and 35.1 at 250, 350, and 500 µm, respectively<fnref n="3" />) added to timeline data may have FWHM between 10" and 30" at 250 µm, 13.3 and 40" at 350 µm, and 20" and 60" at 500 µm. Sources, whose FWHM are determined via `sourceExtractorTimeline`, and were outside these values have been excluded from the analysis.

The flux densities of resolved sources in this range of FWHM are calculated fitting a circular Gaussian function with the FWHM set to the measured SPIRE beam; we explain how this produces the most reliable flux densities for our targets in Sect. 2.5.2.

Although `sourceExtractorSussextractor` produces flux density measurements for the sources it detects and although other software packages have been developed for source extraction within confused extragalactic fields (including software developed specifically for Herschel), we prefer to use the timeline-based source extraction for several reasons. First of all, the SPIRE data are flux calibrated at the timeline level using timeline-based PSF-fitting techniques (Bendo et al. 2013), so the methods that we have used are consistent with the flux calibration measurements themselves and are therefore expected to be more accurate. Second, when the maps are

<footnote n="3">³ http://herschel.esac.esa.int/twiki/pub/Public/SpireCalibrationWeb/beam_release note v1-1.pdf</footnote>
produced from the timeline data, multiple measurements from the timelines are averaged together to produce the surface brightness measured in each pixel. The mapmaking process effectively leads to the loss of information about the emission from the source. Because of this, it is preferable to use the timeline data for performing photometry measurements through point-spread function (PSF) fitting down to a level of 30 mJy, as the timeline data preserves the surface brightness and positional information of each individual measurement (see also, Pearson et al. 2013, and their tests with different methods). We also note that sourceExtractorSussextractor itself is subject to some systematic bias issues related to how it uses the error map for weighting (Smith et al. 2012).

Once we calculated the flux densities of each resolved source, we removed them from the timeline data using the flux density and the width of the Gaussian obtained in the fit. This subtraction step removes bright sources from the data that may either cause problems when fitting fainter sources or possibly hide other objects. We repeat the source identification with sourceExtractorSussexextractor, the source fitting with sourceExtractorTimeline, and the source subtraction for a total of three times on the data. In the second and third iterations, we work with the images where we had subtracted all the sources identified in previous iterations. This procedure increases the number of sources detected by ≈10%, as it allows us to detect sources that may be missed in just one iteration.

## 2.4. Extended sources

Another problem with our method is that sourceExtractorSussexextractor tends to break up extended objects (generally Virgo cluster galaxies) into “families” of point sources, which could potentially contaminate our catalog. To address this problem, we performed a series of tests to determine how to remove bright extended sources a priori. The method we developed defines a mask using the approach of Sextractor (Bertin & Arnouts 1996), which detects a source when a fixed number of contiguous pixels is above a σ-threshold estimated from the background map. We estimate that 70 contiguous pixels above $1.2\sigma$ is a reasonable threshold to remove most of extended sources. We therefore reject sources with a total of 250 µm flux density higher than $\sqrt{70} \times 1.2 \approx 10\sigma \approx 60$ mJy and larger than 0.7 arcmin²; we reject all sources detected by sourceExtractorSussexextractor that conform to these criteria. In this sense, these maps are not really masks in the usual sense of the words but a way to define regions where the sources are rejected.

## 2.5. Quality assessment

In this section, we determine the reliability and the completeness of the catalogs with different tests. Then, we also verify the flux accuracy in the sourceExtractorTimeline fitter, and the precision in the determination of the source position in sourceExtractorSussexextractor.

### 2.5.1. Reliability

The reliability of a catalog can be quantified by the number of false detections found at a given flux density. When the dominant sources of noise are the background fluctuations, a false detection is simply a source that has a nominal signal-to-noise ratio below a fixed threshold. However, when the noise is dominated

by the confusion, or the number of faint unresolved sources inside a beam is high, such a criterion is no longer valid. In a crowded field, the estimated flux density is uncertain, since each detection can be contaminated by neighboring sources.

To quantify the number of false detections due to the instrumental noise, we use a permutation method similar to that used by Smith et al. (2012). We build two maps using the two halves of the eight HeViCS-scans and calculate two difference maps:

$$M_{\text{diff}} = \pm \frac{M_1 - M_2}{2}, \qquad (1)$$

where $M_1$ and $M_2$ are the maps built with the first and the last four scans. This procedure removes the confusion noise and considers only the uncertainties due to the instrument. The subtraction is done in two different permutations because a cosmic ray can give a positive or negative spike depending on the map considered. We search for sources using sourceExtractorSussexextractor, which found 13, 21, 3, and 8 sources in V1, V2, V3, and V4, respectively. However the position of the sources are confined to the external regions of the map in all cases, where the number of bolometer samplings is low. In the catalog, we considered only regions of the field with a sampling redundancy ≥35, therefore ignoring all the sources on the borders. This implies that we can reasonably assume that there are no false detections due to glitches or different noise properties in different scans in the regions considered for the source extraction.

### 2.5.2. Completeness and flux accuracy

To determine the completeness and the flux accuracy of the catalog, we artificially added a Gaussian source at a random position to the timeline data with a FWHM equal to the beam for the appropriate photometric filter. In each HeViCS field, we define a box of $0.5 \times 0.5$ deg² and add a source with flux densities of 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, and 100 mJy. For each flux value, the test is repeated 500 times.

We then extract a list of point sources in the selected box using sourceExtractorSussexextractor and search for the nearest source at a distance equal to at most half the beam. We avoid a larger distances because the fields are crowded, and we can easily find an unrelated source close to the input coordinate.

As a preliminary test, we verified that the completeness and flux accuracies are consistent all over for the four HeViCS fields. Table 2 and Fig. 4 show the results of this test at 250 µm: Cols. 1 and 2 indicate the field selected and the flux density of the source added in the timeline; Cols. 3 and 4 give the percentage of sources detected and not detected; and Col. 5 reports the statistical outliers, defined as the sources with flux densities that differ by more than $3\sigma$ from the measured average. In our case, this value ranges from 0.2% to 2.2%, indicating that the flux density calculated with the sourceExtractorTimeline is robust when the source is detected.

The fraction of high-flux sources detected at 250 µm has more than 95% of the sources that are down to 40 mJy, independent of the chosen fields (left panel of Fig. 4). Below this value, the completeness level decreases rapidly to ~88% at 30 mJy and to ~74% at 20 mJy, as expected by considering the confusion noise level of ~20 mJy at $3\sigma$ found by Nguyen et al. (2010). The flux estimation of sourceExtractorTimeline is quite stable until 20 mJy, despite the rapid decrease of the completeness to 74%, where the recovered flux is on average 21 ± 5 mJy, and agrees with the input value (right panels of Fig. 4). However, below this threshold, the fluxes of sources between 5 and 10 mJy
<table><caption>**Table 2.** Test of completeness and flux accuracy at 250 µm.</caption><thead><tr><th>Field<br>(1)</th><th>F<sub>inp</sub><br>mJy<br>(2)</th><th>N<sub>det</sub><br>%<br>(3)</th><th>N<sub>ndet</sub><br>%<br>(4)</th><th>N<sub>out</sub><br>%<br>(5)</th><th>F<sub>avg</sub><br>mJy<br>(6)</th><th>F<sub>med</sub><br>mJy<br>(7)</th><th>F<sub>err</sub><br>mJy<br>(8)</th></tr></thead><tbody><tr><td rowspan="10">V1</td><td>100</td><td>99.0</td><td>1.0</td><td>0.8</td><td>100</td><td>99</td><td>6</td></tr><tr><td>90</td><td>99.6</td><td>0.4</td><td>2.2</td><td>90</td><td>89</td><td>5</td></tr><tr><td>80</td><td>99.2</td><td>0.8</td><td>2.2</td><td>80</td><td>79</td><td>5</td></tr><tr><td>70</td><td>99.4</td><td>0.6</td><td>2.2</td><td>70</td><td>69</td><td>5</td></tr><tr><td>60</td><td>98.8</td><td>1.2</td><td>0.6</td><td>61</td><td>59</td><td>7</td></tr><tr><td>50</td><td>97.8</td><td>2.2</td><td>2.2</td><td>50</td><td>50</td><td>5</td></tr><tr><td>40</td><td>94.8</td><td>5.2</td><td>2.4</td><td>40</td><td>40</td><td>6</td></tr><tr><td>30</td><td>87.2</td><td>12.8</td><td>2.2</td><td>30</td><td>29</td><td>5</td></tr><tr><td>20</td><td>75.0</td><td>25.0</td><td>1.6</td><td>21</td><td>20</td><td>5</td></tr><tr><td>10</td><td>32.0</td><td>68.0</td><td>1.0</td><td>17</td><td>16</td><td>8</td></tr><tr><td>5</td><td>13.2</td><td>86.8</td><td>0.2</td><td>23</td><td>16</td><td>19</td></tr><tr><td rowspan="10">V2</td><td>100</td><td>100.0</td><td>0.0</td><td>1.6</td><td>100</td><td>100</td><td>6</td></tr><tr><td>90</td><td>100.0</td><td>0.0</td><td>1.8</td><td>90</td><td>90</td><td>5</td></tr><tr><td>80</td><td>99.4</td><td>0.6</td><td>1.0</td><td>81</td><td>80</td><td>6</td></tr><tr><td>70</td><td>99.8</td><td>0.2</td><td>1.4</td><td>70</td><td>69</td><td>6</td></tr><tr><td>60</td><td>98.8</td><td>1.2</td><td>2.2</td><td>60</td><td>60</td><td>6</td></tr><tr><td>50</td><td>97.2</td><td>2.8</td><td>1.8</td><td>50</td><td>49</td><td>5</td></tr><tr><td>40</td><td>95.2</td><td>4.8</td><td>1.6</td><td>41</td><td>40</td><td>6</td></tr><tr><td>30</td><td>91.8</td><td>8.2</td><td>1.0</td><td>31</td><td>30</td><td>6</td></tr><tr><td>20</td><td>75.8</td><td>24.2</td><td>1.6</td><td>21</td><td>21</td><td>5</td></tr><tr><td>10</td><td>40.4</td><td>59.6</td><td>1.4</td><td>15</td><td>14</td><td>6</td></tr><tr><td>5</td><td>16.0</td><td>84.0</td><td>0.2</td><td>19</td><td>14</td><td>16</td></tr><tr><td rowspan="10">V3</td><td>100</td><td>100.0</td><td>0.0</td><td>1.2</td><td>100</td><td>99</td><td>5</td></tr><tr><td>90</td><td>99.8</td><td>0.2</td><td>2.0</td><td>90</td><td>90</td><td>6</td></tr><tr><td>80</td><td>99.8</td><td>0.2</td><td>2.6</td><td>80</td><td>79</td><td>6</td></tr><tr><td>70</td><td>100.0</td><td>0.0</td><td>1.4</td><td>70</td><td>70</td><td>6</td></tr><tr><td>60</td><td>98.4</td><td>1.6</td><td>2.0</td><td>60</td><td>59</td><td>6</td></tr><tr><td>50</td><td>98.0</td><td>2.0</td><td>2.4</td><td>50</td><td>49</td><td>5</td></tr><tr><td>40</td><td>95.2</td><td>4.8</td><td>1.6</td><td>40</td><td>39</td><td>6</td></tr><tr><td>30</td><td>84.2</td><td>15.8</td><td>1.4</td><td>31</td><td>30</td><td>6</td></tr><tr><td>20</td><td>68.8</td><td>31.2</td><td>1.6</td><td>22</td><td>21</td><td>6</td></tr><tr><td>10</td><td>35.6</td><td>64.4</td><td>0.6</td><td>19</td><td>16</td><td>11</td></tr><tr><td>5</td><td>13.4</td><td>86.6</td><td>0.2</td><td>23</td><td>16</td><td>28</td></tr><tr><td rowspan="10">V4</td><td>100</td><td>100.0</td><td>0.0</td><td>2.0</td><td>100</td><td>99</td><td>6</td></tr><tr><td>90</td><td>100.0</td><td>0.0</td><td>2.2</td><td>89</td><td>89</td><td>5</td></tr><tr><td>80</td><td>99.4</td><td>0.6</td><td>1.8</td><td>80</td><td>80</td><td>5</td></tr><tr><td>70</td><td>99.8</td><td>0.2</td><td>1.2</td><td>70</td><td>69</td><td>5</td></tr><tr><td>60</td><td>99.6</td><td>0.4</td><td>2.0</td><td>60</td><td>59</td><td>6</td></tr><tr><td>50</td><td>98.8</td><td>1.2</td><td>1.8</td><td>50</td><td>49</td><td>6</td></tr><tr><td>40</td><td>96.4</td><td>3.6</td><td>1.4</td><td>40</td><td>39</td><td>5</td></tr><tr><td>30</td><td>90.8</td><td>9.2</td><td>1.6</td><td>30</td><td>30</td><td>5</td></tr><tr><td>20</td><td>77.4</td><td>22.6</td><td>1.0</td><td>22</td><td>21</td><td>6</td></tr><tr><td>10</td><td>37.8</td><td>62.2</td><td>0.2</td><td>19</td><td>16</td><td>14</td></tr><tr><td>5</td><td>17.6</td><td>82.4</td><td>1.2</td><td>19</td><td>15</td><td>14</td></tr></tbody></table>

**Notes.** Column 1: field; Col. 2: input flux in mJy; Cols. 3, and 4: percentage of sources detected and not detected; Col. 5: percentage of statistical outliers; Cols. 6–8: average, median, and standard deviation of the flux densities recovered, respectively, excluding the statistical outliers.

are systematically overestimated because the effect called “flux boosting” becomes relevant at these values. This last result also confirms that we find only a “sea” of faint sources, which constitutes the unresolved FIR background, below a flux density of ~20 mJy.

Once verified, the consistency of our results with respect to the chosen field we calculated the completeness and the flux accuracy of our method at all SPIRE wavelengths and for all

the fields, as shown in Fig. 5. The statistics for the tests done at 250, 350, and 500 µm are also given in Table 3. The completeness at all SPIRE bands is around 95% for flux densities of 40 mJy, but the completeness decreases at longer wavelengths at lower fluxes. At 30 mJy, we have a completeness of 88%, 84%, and 82% at 250, 350, and 500 µm, while the completeness at 20 mJy is reduced to 72%, 61%, and 57% at 250, 350, and 500 µm. The errors in position are also quite similar, and the differences in absolute value are due to the different pixel sizes, 6”, 8”, and 12” at 250, 350, and 500 µm, respectively (middle panel of Fig. 5).

## 2.6. Final catalog

The final catalogs consist of 52 020, 42 278, and 18 691 unique sources with fluxes above 20 mJy selected at 250, 350, and 500 µm over an area of ~55 deg². The machine-readable version of the catalog is available online<fnref n="4" />. In the overlapping regions of the different fields, the fluxes have been averaged, and the uncertainties have been added in quadrature.

We considered counterparts to the 250 µm sources at 350 µm and 500 µm, and their fluxes and positions were determined by positional coincidence with the 250 µm sources; positional offsets were required to be within a radius equal to half a beam, which are at 12” and 18” at 350 µm and 500 µm, respectively. As shown in Fig. 6, at 350 and 500 µm about 90% and 75% of the sources have an offset from the 250 µm counterpart below 6” and 9”, respectively. These offsets are smaller than the size of a pixel (8” at 350, and 12” at 500 µm) indicating that the number of spurious objects can be considered negligible. The on line version of the catalogs at 350 and 500 µm has also a flag that indicates if there is an object with flux above 20 mJy within a radius of 12” and 18” at 350 µm and 500 µm at 250 µm, respectively. About 70% at 350 µm and 65% at 500 µm of the sources have a counterpart at 250 µm, and within these, 68% and 46% lie at a radial distance of 9” (about half beam at 250 µm).

# 3. Analysis

## 3.1. Number counts

As a preliminary test, we analyze the number counts for each different field separately to investigate the reliability of the procedures used to extract sources and to perform the photometry. Since the effect of clustering due to the background sources can be considered negligible in our maps because of the large map size, we expect not to find significant variations among different fields, at least, if a significant fraction of the point sources are not within the cluster. This is confirmed in Fig. 7, which shows the number counts of the four HeViCS fields, separately. The units in the ordinate are used to have a direct representation of the evolution in luminosity of the sources, since the number counts in an isotropic Euclidean Universe with no evolution would appear as a straight horizontal line (Longair 1966). There is no relevant differences between the fields with a high number of extended sources (in particular the V2 field) and the one in which this component is negligible. This consistency is verified at each flux value, in particular for faint fluxes, indicating that Virgo cluster galaxies give a negligible contribution to the total number counts. This also suggests that the Virgo cluster does not contain an infrared-bright, optically-faint population of galaxies. There is some difference between different fields for sources at fluxes

<footnote n="4">⁴ http://www.hevics.org/</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Completeness and flux-position accuracy in the catalog at 250 µm. In each panel, the field V1, V2, V3, and V4 are shown as solid, dashed, dotted, and dash-dotted line. **Left panel:** ratio between the numbers of objects detected and the number of sources injected. **Middle panel:** radial offset between the sources injected and the source recovered as a function of the input flux density. **Right panel:** logarithm of the ratio between the average recovered flux density and the input one, as a function of the injected fluxes.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Completeness and flux-position accuracy in the catalog at SPIRE wavelengths. In each panel, the black triangles with solid line, blue stars with dash line, and red circle with dot line represent the results at 250, 350, and 500 µm. **Left panel:** ratio between the numbers of objects detected and the number of sources injected. **Middle panel:** radial offset between the sources injected and the source recovered as a function of the input flux density. **Right panel:** logarithm of the ratio between the average recovered flux density and the input one, as a function of the injected fluxes.</figcaption></figure>

above 150 mJy, but these differences are mainly driven by the relatively small number of sources at those fluxes.

The HeViCS number counts for the combined fields corrected for the incompleteness are shown in Fig. 8 with the result of other analyses: Herschel Astrophysical TeraHertz Large Area Survey (H-ATLAS, Clements et al. 2010), Herschel Multi-Tiered Extragalactic Survey (HerMES) SDP (Oliver et al. 2010) and DR1 (Béthermin et al. 2012a), together with a statistical analysis from Glenn et al. (2010). Finally, we also include the results of (Balloon-borne Large-Aperture Submillimeter Telescope BLAST, Pascale et al. 2008), an experiment realized using a precursor of the SPIRE instrument.

A characteristic feature at all the wavelengths is an increase in the number counts at $F \le 200$ mJy at 250, 350, and 500 µm, indicating a strong evolution for the galaxy populations at faint fluxes. Comparison of the number counts at different wavelengths shows that the slope of the faint end where the statistics are good (40 $\le F \le 150$ mJy) steepens going from 250 µm to the longer wavelengths. Below this flux level, the number counts flatten and decrease because of the confusion.

There is quite good agreement between our data and the results of other Herschel surveys at faint fluxes below 200

and 180 mJy at 250 and 350 µm, respectively. However at 500 µm, the number of sources at low flux densities is higher with respect to other surveys, which is probably an effect of our source extraction method. At flux densities below 30 mJy, we found ~90% more sources with respect to the HerMES SDP survey; this number decreases to ~40% when we consider sources with a flux above 60 mJy. A similar tendency is seen comparing our results to the HerMES DR1 (Béthermin et al. 2012a) data and the deeper HerMES statistical analysis done by Glenn et al. (2010), which resolved 43% of the FIR background at 500 µm. The subtraction of the sources from the timeline allowed us to recover more objects with respect to other methods, and this effect is more relevant when the beam is larger, because the probability in this case that sources are blended together, or partially covered by a brighter source is higher with respect to smaller beam-sizes. The increased slope of the number counts at wavelengths above 350 µm might be due to the presence of proto-spheroidal galaxies, which are mostly evolving passively at lower redshift, and became more active at redshift $z \ge 1.5$. Because of the strongly negative k-correction at sub-mm wavelengths, the number of these galaxies is weakly dependent on the redshift, resulting in a steep rise of the number counts at $\lambda \ge 350$ µm.
<table><caption>**Table 3.** Test of completeness and flux accuracy in the HeViCS fields.</caption><thead><tr><th>Wlg<br>μm<br>(1)</th><th>F<sub>inp</sub><br>mJy<br>(2)</th><th>N<sub>det</sub><br>%<br>(3)</th><th>N<sub>ndet</sub><br>%<br>(4)</th><th>N<sub>out</sub><br>%<br>(5)</th><th>F<sub>avg</sub><br>mJy<br>(6)</th><th>F<sub>med</sub><br>mJy<br>(7)</th><th>F<sub>err</sub><br>mJy<br>(8)</th></tr></thead><tbody><tr><td rowspan="10">250</td><td>100</td><td>99.8</td><td>0.2</td><td>1.4</td><td>100</td><td>99</td><td>6</td></tr><tr><td>90</td><td>99.9</td><td>0.1</td><td>2.1</td><td>90</td><td>89</td><td>5</td></tr><tr><td>80</td><td>99.5</td><td>0.5</td><td>1.9</td><td>80</td><td>79</td><td>6</td></tr><tr><td>70</td><td>99.8</td><td>0.2</td><td>1.6</td><td>70</td><td>69</td><td>6</td></tr><tr><td>60</td><td>98.9</td><td>1.1</td><td>1.7</td><td>60</td><td>59</td><td>6</td></tr><tr><td>50</td><td>98.0</td><td>2.1</td><td>2.1</td><td>50</td><td>49</td><td>5</td></tr><tr><td>40</td><td>95.4</td><td>4.6</td><td>1.8</td><td>40</td><td>40</td><td>6</td></tr><tr><td>30</td><td>88.5</td><td>11.5</td><td>1.6</td><td>30</td><td>30</td><td>6</td></tr><tr><td>20</td><td>74.2</td><td>25.8</td><td>1.5</td><td>21</td><td>21</td><td>5</td></tr><tr><td>10</td><td>36.4</td><td>63.5</td><td>0.8</td><td>17</td><td>15</td><td>10</td></tr><tr><td>5</td><td>15.0</td><td>85.0</td><td>0.4</td><td>21</td><td>16</td><td>19</td></tr><tr><td rowspan="10">350</td><td>100</td><td>99.9</td><td>0.1</td><td>1.7</td><td>100</td><td>99</td><td>6</td></tr><tr><td>90</td><td>99.9</td><td>0.1</td><td>1.6</td><td>90</td><td>89</td><td>6</td></tr><tr><td>80</td><td>99.7</td><td>0.3</td><td>1.7</td><td>80</td><td>79</td><td>6</td></tr><tr><td>70</td><td>99.5</td><td>0.5</td><td>1.9</td><td>70</td><td>69</td><td>6</td></tr><tr><td>60</td><td>98.9</td><td>1.1</td><td>1.5</td><td>60</td><td>59</td><td>6</td></tr><tr><td>50</td><td>97.3</td><td>2.7</td><td>1.5</td><td>50</td><td>49</td><td>6</td></tr><tr><td>40</td><td>93.0</td><td>7.0</td><td>1.9</td><td>40</td><td>39</td><td>6</td></tr><tr><td>30</td><td>84.2</td><td>15.8</td><td>1.6</td><td>31</td><td>30</td><td>6</td></tr><tr><td>20</td><td>61.5</td><td>38.5</td><td>1.1</td><td>23</td><td>22</td><td>7</td></tr><tr><td>10</td><td>25.1</td><td>74.9</td><td>0.5</td><td>18</td><td>17</td><td>8</td></tr><tr><td>5</td><td>11.7</td><td>88.3</td><td>0.5</td><td>20</td><td>18</td><td>10</td></tr><tr><td rowspan="10">500</td><td>100</td><td>100.0</td><td>0.0</td><td>1.3</td><td>99</td><td>99</td><td>6</td></tr><tr><td>90</td><td>100.0</td><td>0.0</td><td>0.9</td><td>90</td><td>89</td><td>6</td></tr><tr><td>80</td><td>99.9</td><td>0.1</td><td>1.6</td><td>79</td><td>79</td><td>6</td></tr><tr><td>70</td><td>99.6</td><td>0.4</td><td>1.3</td><td>70</td><td>69</td><td>6</td></tr><tr><td>60</td><td>99.3</td><td>0.7</td><td>1.1</td><td>60</td><td>59</td><td>6</td></tr><tr><td>50</td><td>97.5</td><td>2.5</td><td>1.5</td><td>50</td><td>49</td><td>6</td></tr><tr><td>40</td><td>94.1</td><td>5.9</td><td>1.3</td><td>40</td><td>40</td><td>6</td></tr><tr><td>30</td><td>82.1</td><td>17.9</td><td>1.2</td><td>31</td><td>30</td><td>6</td></tr><tr><td>20</td><td>57.9</td><td>42.1</td><td>0.7</td><td>23</td><td>23</td><td>7</td></tr><tr><td>10</td><td>24.1</td><td>75.9</td><td>0.5</td><td>21</td><td>19</td><td>12</td></tr><tr><td>5</td><td>10.5</td><td>89.5</td><td>0.5</td><td>21</td><td>18</td><td>11</td></tr></tbody></table>

**Notes.** Column 1: wavelength; Col. 2: input flux in mJy; Cols. 3, and 4: percentage of sources detected and not detected; Col. 5: percentage of statistical outliers; Cols. 6–8: average, median, and standard deviation of the flux densities recovered excluding statistical outliers, respectively.

Together with the results of other analyses, we compare the HeViCS number of counts with other theoretical models. There are two main approaches to theoretically understanding the number counts of galaxies. A first method, “forward evolution”, is based on semi-analytic models of galaxy formation and evolution. In these models, the time evolution of the observational properties of galaxies are predicted using prescriptions for the physics of the baryonic components that are embedded in the large structure skeleton of the dark matter, following the standard ΛCDM cosmology. An example of this method is shown in Lacey et al. (2008), which has been applied with success to Spitzer data, and extended to Herschel FIR wavelengths (Lacey et al. 2010).

The second method, “backward evolution” (e.g. Béthermin et al. 2011), starts from the observed luminosity and density function at $z=0$ and then assumes a parametric form for the luminosity evolution. A set of spectral energy distribution (SED) templates representative of different galaxy populations are made to evolve in density and luminosity according

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Normalized histograms of the offset between sources at 250 µm and the associated counterparts at 350 (black solid line) and 500 µm (red dashed line).</figcaption></figure>

to the luminosity function parameterization. This approach can give different results according to the different choices of the galaxy SED templates, and the different physical processes that are taken into account in the modeling. We choose the model of Valiante et al. (2009) for our comparison with the data, which includes an evolution with the redshift for the AGN contribution to the total FIR luminosity, and the model of Negrello et al. (2007), which includes the contribution of gravitational lensing to the observed fluxes, a contribution not negligible at wavelengths $\ge 500$ µm (e.g. Negrello et al. 2010).

Analyzing GOODS-Herschel data, Elbaz et al. (2011) found that the evolution of IR galaxies follows two distinct star-formation mechanisms. Most galaxies follow a main sequence with a tight correlation between stellar mass and star formation, independent of luminosity and redshift. A smaller galaxy population, less than 20%, shows a starburst phase, with high specific star-formation rate (e.g., Magnelli et al. 2012). Supported by these results, Béthermin et al. (2012b) proposed a model that considers the AGN contribution, gravitational lensing, and assumed only two types of galaxy SEDs for the main sequence and the starburst population. They also introduce an evolution with redshift for the SED templates, following the recent results of Magdis et al. (2012), which showed that the mean radiation field of a galaxies is stronger at higher redshift. There is a substantial amount of observational evidence that suggests a redshift evolution of the total IR luminosity and dust temperature (e.g. Amblard et al. 2010; Elbaz et al. 2010; Hwang et al. 2010; Casey et al. 2012; Magnelli et al. 2014).

We compare our observed number counts with the predicted number counts at each SPIRE wavelength in Fig. 8, according to the models of Lacey et al. (2010); Negrello et al. (2007); Valiante et al. (2009) and Béthermin et al. (2012b). The Lacey et al. (2010) model predictions tend to exceed the observed ones at all fluxes. This could be because the model employs a top-heavy stellar initial mass function (similarly to Baugh et al. 2005) that produces more dust per unit star-formation rate. In this way, they match the bright end of galaxy counts at 850 µm, but the number counts are systematically overestimated at lower fluxes. The model of Valiante et al. (2009) gives the best agreement with our data (see also Glenn et al. 2010), indicating the importance of
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Number counts estimated in each field at 250 (**left**), 350 (**middle**), and 500 (**right**) µm. Green diamonds, black circles, blue triangles, and red stars indicate V1, V2, V3, and V4 field, respectively.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Number counts (black circles) at 250 (left), 350 (middle), and 500 (right) µm. Magenta triangles, red stars, blue diamonds, green triangles, and cyan hexagons show the observations from H-ATLAS (Clements et al. 2010), HerMES SDP (Oliver et al. 2010), BLAST (Patanchon et al. 2009), HerMES DR1 (Béthermin et al. 2012a), and the results of a *P*(D) (probability of deflection) analysis applied to the HerMES SDP data from Glenn et al. (2010), respectively. The green, black, blue, and red solid lines show the models of Lacey et al. (2010), Negrello et al. (2007), Valiante et al. (2009), and Béthermin et al. (2012b).</figcaption></figure>

an AGN contribution in shaping the galaxy SEDs, although the source fluxes $\ge 100$ mJy at 350 and 500 µm tend to be overestimated. The models by Negrello et al. (2007) reproduce quite well the slope and the evolution at lower fluxes but overestimate the number counts at 250 and 350 µm. At 500 µm, the comparison with our data is good, but the increase in the number of counts occurs at lower fluxes with respect to our data.

Generally, the models tend to overpredict the observed counts at bright flux and underpredict them at faint ones. These discrepancies between model predictions and observations can arise from several factors. For example, the evolution of the SEDs with the redshift is poorly known, and the SEDs library is well calibrated only in the local Universe, where the *k*-correction is small. Another effect that can play a role is the AGN contribution to the dust heating (e.g. Almaini et al. 1999; Granato et al. 1997; Farrah et al. 2007; Sajina et al. 2008), although recent results suggest that dust emission due to heating by stars dominates that from the heating by AGNs at wavelengths larger than ~30 µm (Hatziminaoglou et al. 2010). The Chandra X-ray images of luminous BzK galaxies at $z \sim 2$ (Daddi et al. 2007) exhibit soft-to-hard X-ray ratio typical of heavily obscured AGN, but deeper observations have shown that only 25% of luminous BzK galaxies host heavily obscured AGN (Alexander et al. 2011) and the rest are star-forming galaxies with unobscured AGN. In any case, the difficulties in reproducing the number counts at fluxes $\ge 100$ mJy underline the importance of the Rayleigh-Jeans side of the SED to characterize the

evolution of FIR galaxy population, a region for which Herschel is fundamental.

### 3.1.1. Contribution to the number counts from different redshift bins

We searched for the optical counterparts of the HeViCS FIR sources in the *r*-band of the SDSS DR10 survey (Ahn et al. 2014) using the likelihood ratio method described in Sutherland & Saunders (1992). This method associates a parameter, the reliability ($R$), to each optical source, that quantifies the probability that the SDSS source considered is the true optical counterpart of the FIR object of reference. To determine $R$, we account for the possibility that the counterpart is below the optical magnitude limit with the probability to find serendipitously a spurious source in a given position. To define the optical matching sample, we extracted a 55 deg$^2$ cutout from the *r*-band SDSS catalog, perfectly overlapping with the HeViCS surveyed area. This contains 1 714 999 sources with a completeness limit of 95% at $r \approx 22.2$ AB mag in optical imaging (Ahn et al. 2014). The Milky Way extinction corrected SDSS *r*-magnitudes were used to compute the surface density distribution of the optical objects.

The positional uncertainties in the sky coordinates were derived for each infrared source from the timeline fitting; however, we further considered that error radii less than a few arcsec are unrealistic given the possibility of systematic offsets (see Smith et al. 2012, for details). We thus assumed 2.′′5, 3.′′4, and 5′′ as
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Redshift histogram for the sources that have an optical counterpart in the SDSS catalog with a reliability above 80% at 250 (black dashed line), 350 (blue dotted line), and 500 (red dot-dashed line) µm. Spectroscopic redshifts have been considered where available, photometric redshift are considered otherwise.</figcaption></figure>

minimum error radii for the HeViCS 250, 350, and 500 µm sources, respectively. For each HeViCS source, the r-band identification with the $R > 80\%$ was taken as the real counterpart, corresponding to 40%, 25%, and 19% of the total sources at 250, 350, and 500 µm, respectively. The average redshift of the catalog having significant detections at all SPIRE wavelengths is $z \sim 0.3 \pm 0.22$ and 0.25 (16–84 percentile). The relatively large spread arises because ~80% of the sources are at $z < 0.5$, while the remaining 20% are at higher redshift with $0.5 < z < 3.5$, as shown in the histograms of Fig. 9. In this figure, we have considered the spectroscopic redshift where available. For the others, we use the photometric redshift obtained with a robust fit on spectroscopically observed objects with similar colors and inclination angle. In Fig. 10, we show the offset between the HeViCS source and the associated optical counterparts. Most of the sources are within a radius smaller than 6", which is the smallest pixel size of SPIRE instrument, confirming the reliability of this procedure.

Clearly, the redshift distribution considered for the sources with an optical counterpart in SDSS is biased towards low redshift galaxies. Because the SDSS catalog is 95% complete for $r = 22.2$ AB mag, the cross correlation with our FIR catalogs give us the possibility to estimate the contribution of the low redshift component to the total number counts. In Fig. 11, we show the number of counts of the HeViCS sources, the subsample with optical counterparts and redshift estimation, and the component without redshift. Above 100 mJy, the number counts are dominated by the component with the optical counterpart that forms the bulk of the bright end of the distribution, while the component without redshift becomes more important at flux densities below ~60 mJy.

The number counts of objects with no optical counterparts for any SPIRE wavelength are steeper than the sample with optical counterparts, indicating a more rapid density evolution for the former component. Moreover, the contribution from the component with no optical counterparts to the total number counts becomes more relevant for an increasing wavelength. This is because the peak of IR SED is around 200 µm, and by increasing the wavelength, we consider sources statistically at higher redshift, where the SDSS is less efficient in detecting

sources. At fluxes above 100 mJy, the number counts at all SPIRE wavelengths, as expected, is dominated by the component with optical counterpart, which are on average star-forming galaxies at $z \sim 0.3$.

## 3.2. Far-infrared colors

In this section, we investigate the possibility, as seen, for example, in H-ATLAS (Amblard et al. 2010), to use the color-color diagram at SPIRE wavelengths to statistically constrain the redshift distribution of the sources, as suggested in Hughes et al. (2002). Figure 12 shows a color-color diagram for the HeViCS sources that have an optical counterpart at 250 µm with a reliability above 80%, a flux density at 250 µm above 30 mJy, and a magnitude in r-band below 22.2. We also construct a subsample of ~2000 sources with a reliable flux estimation at all SPIRE wavelengths. We divided this subsample in redshift bins of $\Delta z = 0.1$ to have a comparable number of sources in each bin, except for the highest ones, where we are limited by small-number statistics. The sources at lower redshift span a large region of the diagram, showing an overdensity in the region with $F_{250}/F_{350} > 1.7$. The sources at higher redshift mostly populate the region with $F_{250}/F_{350} < 1.7$. However, because of the range in temperature and emissivity, this region of the diagram is also heavily populated by sources at all redshifts.

The shift of the distribution of sources to lower $F_{250}/F_{350}$ at increasing redshift can be explained by the following considerations: at low redshifts, the SED of an average galaxy peaks at wavelengths below 250 µm. The $F_{250}/F_{350}$ ratio appears high and the $F_{500}/F_{350}$ ratio appears low, as both ratios sample the Rayleigh-Jeans side of the SED. As redshift increases, this peak moves beyond the 250 µm band. Thus, the $F_{250}/F_{350}$ ratio decreases as the bands shift with redshift into the region near the peak of the SED and onto the Wien side of the SED, and the $F_{500}/F_{350}$ ratio increases, as the 500 µm band approach the peak of the SED.

Most of the sources occupy a region with $F_{250}/F_{350} > 0.6$ and $F_{500}/F_{350} < 1$. This is because sources at lower redshift are on average detected at lower SPIRE wavelengths, whereas higher redshift sources have statistically higher flux values at longer wavelengths (see for example Hughes et al. 2002, and their analysis). To better investigate this hypothesis in Fig. 12, we showed a reference curve for a modified black body with $\beta = 2$, $10 \le T \le 30$ and a redshift equal to the average of the bin. Figure 12 shows that at fixed temperature sources move towards the top left side of the diagram. However, increasing the redshift shows a shift of the curve towards the same direction, reinforcing the idea that increasing the redshift on average sources tend to occupy the region in the top left of the diagram (as shown also in Amblard et al. 2010; Hughes et al. 2002).

From the color-color diagram, we can infer statistical information about the redshift distribution of our sources. Three different regions can be defined: a region with $F_{250}/F_{350} < 0.8$; a region with $0.8 < F_{250}/F_{350} < 1.7$; and a region with $F_{250}/F_{350} > 1.7$. For the range of temperature and $\beta$ parameters chosen, the region with lower $F_{250}/F_{350}$ has on average redshift $z \ge 2$ with a higher $F_{500}/F_{350}$ because of the increased emission redshifted to 500 µm. In the subset of our sample with SDSS counterparts, the number of sources with these values of the ratios is negligible. The region with $0.8 < F_{250}/F_{350} < 1.7$ is the most populated, because there is a large set of temperature, $\beta$, and redshifts that fall in this range. The region with $F_{250}/F_{350} > 1.7$ on average is populated by sources at $z < 1$. The FIR sources exhibit a variation with redshift that follows the
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Offsets between the HeViCS coordinates and its associate optical counterpart. In the *left*, *middle*, and *right* panels, the offsets at 250, 350, and 500 µm are shown, respectively.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Number counts (black circles) at 250 (left), 350 (middle), and 500 (right) µm with the component given by the sources with SDSS counterparts (blue stars) and no-SDSS counterpart (red triangles).</figcaption></figure>

trend from bottom right to top left along a curve with increasing redshift.

**4. Conclusions**

In this paper, we present three point source catalogs selected at 250, 350, and 500 µm extracted from the HeViCS survey (Davies et al. 2010) over an area of ~55 deg². The FWHM and flux density of the sources have been estimated using the implementation in HIPE 9.0.0 (Ott 2010) of *sourceExtractorTimeline*, a timeline-based point source fitter. Then we subtracted the sources from the timeline data using the same Gaussian function that we obtained in the fit. This subtraction step allowed us to detect objects that are fainter or hidden by other close sources, improving the robustness of our algorithm in terms of identifying sources.

We performed various tests to verify the goodness, reliability, and flux accuracy of the catalogs. From these analyses, we found that the number of false detections due to background variations or glitches is negligible in the region selected for the source extraction. The SPIRE catalogs are about 95% complete at 40 mJy, but the completeness decreases at longer wavelengths at lower fluxes. We find a completeness of 88%, 84%, and 82% at 250, 350, and 500 µm at 30 mJy, while the completeness drops to 72%, 61%, and 57% at 250, 350, and 500 µm at 20 mJy, confirming that the confusion noise becomes relevant both for flux estimation and for source detection at this threshold.

We analyze the number counts in the four HeViCS fields, and the results are consistent at each flux value, in particular for faint fluxes, indicating that Virgo cluster galaxies have a negligible contribution to the total number counts. Another implication of these results is that the Virgo cluster does not contain a significant population infrared-bright, optically-faint of galaxies.

We find an increase of the slope in the number counts at $F \le 200$ mJy at all SPIRE wavelengths, indicating a strong evolution in luminosities for these galaxies, which agrees with previous studies. Comparison of the number counts at different wavelengths shows that the slope of the faint end steepens going from 250 µm to longer wavelengths. Below 40 mJy, the number counts flatten and decrease because of the confusion. This increase in slope might be due to the presence of actively star-forming galaxies at $z \ge 1.5$ that are passive at low redshift. No model is able to reproduce the bright end of the distribution at all wavelengths, which shows the importance of the Rayleigh-Jeans side of the SED in characterizing the evolution of FIR galaxies population and the importance of the Herschel SPIRE instrument that investigates this region. Nevertheless, there is a relatively small number of sources at the bright end of the number counts, thus precluding definitive conclusions. A point to be better addressed in future modeling efforts is the impact of gravitational lensing on the number counts, which is mostly relevant at 500 µm, and the possibility of an evolving SEDs for FIR galaxies with redshift, as found by Magnelli et al. (2011) and others.

We cross-correlated the SPIRE catalogs with the SDSS DR10 survey (Ahn et al. 2014) and found optical
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Position in a color–color diagram at different redshift bins of the HeViCS sources with photometric (black stars), spectroscopic (red stars) redshift, and reliable optical counterparts. The red solid line shows a reference curve for a modified black body with $\beta = 2$, $10 \le T \le 30$, and the redshift equals to the average of the bin. Rhombus, triangles, and stars indicate a modified black body with $\beta = 2$ and temperature of 10, 20, and 30 K, respectively.</figcaption></figure>

counterparts for a subsample of the HeViCS sources. For sources with both photometric and spectroscopic redshift determinations, most sources lie at $z < 0.5$, with an average redshift of $z \sim 0.3 \pm 0.22$ and 0.25 (16–84 percentile).

We divided the catalog into sources with and without optical counterparts and estimated the contribution at each SPIRE wavelength to the total number counts from the low redshift component. Above 100 mJy, the number of counts is dominated by the component with the optical counterpart (Fig. 11), which forms the bulk of the bright end of the distribution. Sources with

no optical counterpart underwent a more rapid evolution with respect to the sources with SDSS counterpart.

We then constructed a color–color diagram of the sources with a known redshift. Most of the sources occupy a region with $F_{250}/F_{350} > 0.6$ and $F_{500}/F_{350} < 1$. We show that most of our sources are consistent with the chosen ranges of temperature, $\beta$, and $z$. However, this variation must be considered in a statistical sense, because sources with a lower redshift populate larger regions with an overdensity towards the bottom right corner of the diagram ($F_{250}/F_{350} > 1.7$), and sources with high redshift
occupy the top left region of the diagram ($F_{250}/F_{350} < 0.8$), as shown in Fig. 12. For sources with known redshift, their position in the SPIRE color–color diagram tends to move towards the top-left corner as redshift increases.

**Acknowledgements.** We warmly thank the referee for his/her constructive comments and suggestions. We would like to thank D. Munro for freely distributing his Yorick programming language (available at http://www.maumae.net/yorick/doc/index.html). S.B., L.H., S.Z., S.diS.A. are supported through the ASI-INAF agreement I/016/07/0 and I/009/10/0. C.P. was also supported by PRIN-INAF 2009/11 grant (extended to 2012). C.P. acknowledges support from the Science and Technology Foundation (FCT, Portugal) through the Postdoctoral Fellowship SFRH/BPD/90559/2012, PEst-OE/FIS/UI2751/2014, and PTDC/FIS-AST/2194/2012. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the US Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofísica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

## References

Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543

Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2014, ApJS, 211, 17

Alexander, D. M., Bauer, F. E., Brandt, W. N., et al. 2011, ApJ, 738, 44

Almaini, O., Lawrence, A., & Boyle, B. J. 1999, MNRAS, 305, L59

Amblard, A., Cooray, A., Serra, P., et al. 2010, A&A, 518, L9

Aretxaga, I., Hughes, D. H., Coppin, K., et al. 2007, MNRAS, 379, 1571

Auld, R., Bianchi, S., Smith, M. W. L., et al. 2013, MNRAS, 428, 1880

Baugh, C. M., Lacey, C. G., Frenk, C. S., et al. 2005, MNRAS, 356, 1191

Bendo, G. J., Wilson, C. D., Pohlen, M., et al. 2010, A&A, 518, L65

Bendo, G. J., Boselli, A., Dariush, A., et al. 2012, MNRAS, 419, 1833

Bendo, G. J., Griffin, M. J., Bock, J. J., et al. 2013, MNRAS, 433, 3062

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

Béthermin, M., Dole, H., Lagache, G., Le Borgne, D., & Penin, A. 2011, A&A, 529, A4

Béthermin, M., Le Floc'h, E., Ilbert, O., et al. 2012a, A&A, 542, A58

Béthermin, M., Daddi, E., Magdis, G., et al. 2012b, ApJ, 757, L23

Binggeli, B., Sandage, A., & Tammann, G. A. 1985, AJ, 90, 1681

Blain, A. W., Ivison, R. J., & Smail, I. 1998, MNRAS, 296, L29

Boquien, M., Calzetti, D., Combes, F., et al. 2011, AJ, 142, 111

Casey, C. M., Berta, S., Béthermin, M., et al. 2012, ApJ, 761, 140

Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772

Ciesla, L., Boselli, A., Smith, M. W. L., et al. 2012, A&A, 543, A161

Clements, D. L., Rigby, E., Maddox, S., et al. 2010, A&A, 518, L8

Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156

Davies, J. I., Baes, M., Bendo, G. J., et al. 2010, A&A, 518, L48

Devlin, M. J., Ade, P. A. R., Aretxaga, I., et al. 2009, Nature, 458, 737

Eales, S., Dunne, L., Clements, D., et al. 2010, PASP, 122, 499

Elbaz, D., Hwang, H. S., Magnelli, B., et al. 2010, A&A, 518, L29

Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119

Farrah, D., Bernard-Salas, J., Spoon, H. W. W., et al. 2007, ApJ, 667, 149

Ferrarese, L., Côté, P., Cuillardre, J.-C., et al. 2012, ApJS, 200, 4

Fixsen, D. J., Dwek, E., Mather, J. C., Bennett, C. L., & Shafer, R. A. 1998, ApJ, 508, 123

Frayer, D. T., Sanders, D. B., Surace, J. A., et al. 2009, AJ, 138, 1261

Glenn, J., Conley, A., Béthermin, M., et al. 2010, MNRAS, 409, 109

Granato, G. L., Danese, L., & Franceschini, A. 1997, Mem. Soc. Astron. It., 68, 39

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

Hatziminaoglou, E., Omont, A., Stevens, J. A., et al. 2010, A&A, 518, L33

Herschel, W. 1800, Roy. Soc. Lond. Philosoph. Trans. Ser. I, 90, 255

Hughes, D. H., Aretxaga, I., Chapin, E. L., et al. 2002, MNRAS, 335, 871

Hwang, H. S., Elbaz, D., Magdis, G., et al. 2010, MNRAS, 409, 75

Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2008, MNRAS, 385, 1155

Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2010, MNRAS, 405, 2

Lagache, G., Puget, J.-L., & Gispert, R. 1999, Ap&SS, 269, 263

Lagache, G., Dole, H., Puget, J.-L., et al. 2004, ApJS, 154, 112

Le Floc'h, E., Papovich, C., Dole, H., et al. 2005, ApJ, 632, 169

Longair, M. S. 1966, MNRAS, 133, 421

Low, F. J., Rieke, G. H., & Gehrz, R. D. 2007, ARA&A, 45, 43

Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6

Magnelli, B., Elbaz, D., Chary, R. R., et al. 2009, A&A, 496, 57

Magnelli, B., Elbaz, D., Chary, R. R., et al. 2011, A&A, 528, A35

Magnelli, B., Lutz, D., Santini, P., et al. 2012, A&A, 539, A155

Magnelli, B., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, A86

Murphy, E. J., Chary, R.-R., Dickinson, M., et al. 2011, ApJ, 732, 126

Niemi, S.-M., Somerville, R. S., Ferguson, H. C., et al. 2012, MNRAS, 421, 1539

Negrello, M., Perrotta, F., González-Nuevo, J., et al. 2007, MNRAS, 377, 1557

Negrello, M., Hopwood, G., De Zotti, G., et al. 2009, Science, 330, 800

Nguyen, H. T., Schulz, B., Levenson, L., et al. 2010, A&A, 518, L5

Oliver, S. J., Rowan-Robinson, M., & Saunders, W. 1992, MNRAS, 256, 15P

Oliver, S., Mann, R. G., Carballo, R., et al. 2002, MNRAS, 332, 536

Oliver, S. J., Wang, L., Smith, A. J., et al. 2010, A&A, 518, L21

Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614

Ott, S. 2010, Astronomical Data Analysis Software and Systems XIX, 434, 139

Papovich, C., Rudnick, G., Le Floc'h, E., et al. 2007, ApJ, 668, 45

Pascale, E., Ade, P. A. R., Bock, J. J., et al. 2008, ApJ, 681, 400

Patanchon, G., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 707, 1750

Pearson, C., Lim, T., North, C., et al. 2013, Exp. Astron., 37

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Puget, J.-L., Abergel, A., Bernard, J.-P., et al. 1996, A&A, 308, L5

Sajina, A., Yan, L., Lutz, D., et al. 2008, ApJ, 683, 659

Saunders, W., Rowan-Robinson, M., Lawrence, A., et al. 1990, MNRAS, 242, 318

Savage, R. S., & Oliver, S. 2007, ApJ, 661, 1339

Shupe, D. L., Rowan-Robinson, M., Lonsdale, C. J., et al. 2008, AJ, 135, 1050

Smith, A. J., Wang, L., Oliver, S. J., et al. 2012a, MNRAS, 419, 377

Smith, D. J. B., Dunne, L., da Cunha, E., et al. 2012b, MNRAS, 427, 703

Sutherland, W., & Saunders, W. 1992, MNRAS, 259, 413

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