© ESO, 2015

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 comoving star-formation rate density (SFRD), because the optical/ultraviolet emission originating from stars must be corrected for IR extinction. Moreover, if the bulk of far-IR (FIR) luminosity 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 λ> 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; Boquien 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 ~ 1, luminous IR galaxies contribute about 50% to the SFRD (Murphy et al. 2011), while the SFRD increases at higher redshift, z ~ 2 because we have an important contribution from more active star-forming galaxies, or ultraluminous IR galaxies (L_IR ≥ 10¹²L_⊙, 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 ≥ 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.

Fig. 1 HeViCS image at 250 μm with detected sources (red points). The blue squares indicate regions with the full coverage.

2. Methodology

The determination of the position and the flux density of the sources have been performed in different steps. From the starting map, presented in Sect. 2.1, an estimate of the background has been subtracted (Sect. 2.2). The source positions have been determined from this background-subtracted map using 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).

2.1. Data

The four HeViCS fields have been observed in two orthogonal scans at a sampling rate of 60′′ s^-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 (naiveScanmapper 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)¹ 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². 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 σ 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 × M_med − 1.5 × M_avg, where M_med and M_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.

Fig. 2 V2 field region before (left panel) and after (right panel) background subtraction.

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.

To test the reliability of this background estimation, we selected 45 boxes of 0.5 × 0.5 deg² 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 ≤ N_bef ≤ 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 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σ 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 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³) 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 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 sourceExtractorSussextractor, 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 sourceExtractorSussextractor 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σ 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{\mathrm{70}}\mathrm{\times }\mathrm{1.2}\mathrm{\approx }\mathrm{10}\mathit{\sigma }\mathrm{\approx }\mathrm{60}$ mJy and larger than 0.7 arcmin²; we reject all sources detected by sourceExtractorSussextractor 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 sourceExtractorSussextractor.

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: ${\mathit{M}}_{\mathrm{diff}}\mathrm{=}\mathrm{\pm }\frac{{\mathit{M}}_{\mathrm{1}}\mathrm{-}{\mathit{M}}_{\mathrm{2}}}{\mathrm{2}}\mathit{,}$(1)where M₁ and M₂ 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 sourceExtractorSussextractor, 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 × 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 sourceExtractorSussextractor 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σ 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σ 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 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.

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.

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.

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⁴. 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 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 (Baloon-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 ≲ 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 ≲ F ≲ 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 beamsizes. 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 ≥ 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 λ ≥ 350 μm.

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

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 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 ≳500 μm (e.g. Negrello et al. 2010).

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.

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

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 an AGN contribution in shaping the galaxy SEDs, although the source fluxes ≳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 ~ 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 ≳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 (ℛ), 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 ℛ, 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² 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 ≈ 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 minimum error radii for the HeViCS 250, 350, and 500 μm sources, respectively. For each HeViCS source, the r-band identification with the ℛ > 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 ~ 0.3 ± 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.

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 availables, photometric redshift are considered otherwise.

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 ~ 0.3.

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.

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

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 Δ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₂₅₀/F₃₅₀> 1.7. The sources at higher redshift mostly populate the region with F₂₅₀/F₃₅₀< 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₂₅₀/F₃₅₀ 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₂₅₀/F₃₅₀ ratio appears high and the F₅₀₀/F₃₅₀ 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₂₅₀/F₃₅₀ 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₅₀₀/F₃₅₀ ratio increases, as the 500 μm band approach the peak of the SED.

Most of the sources occupy a region with F₂₅₀/F₃₅₀> 0.6 and F₅₀₀/F₃₅₀< 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 β = 2, 10 ≤ T ≤ 30 and a redshift equal to the average of the bin. Figure 12 show 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₂₅₀/F₃₅₀< 0.8; a region with 0.8 <F₂₅₀/F₃₅₀< 1.7; and a region with F₂₅₀/F₃₅₀> 1.7. For the range of temperature and β parameters chosen, the region with lower F₂₅₀/F₃₅₀ has on average redshift z ≥ 2 with a higher F₅₀₀/F₃₅₀ 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₂₅₀/F₃₅₀< 1.7 is the most populated, because there is a large set of temperature, β, and redshifts that fall in this range. The region with F₂₅₀/F₃₅₀> 1.7 on average is populated by sources at z< 1. The FIR sources exhibit a variation with redshift that follows the trend from bottom right to top left along a curve with increasing redshift.

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 β = 2, 10 ≤ T ≤ 30, and the redshift equals to the average of the bin. Rhombus, triangles, and stars indicate a modified black body with β = 2 and temperature of 10, 20, and 30 K, respectively.

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 ≲ 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 ≥ 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 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 ~ 0.3 ± 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₂₅₀/F₃₅₀> 0.6 and F₅₀₀/F₃₅₀< 1. We show that most of our sources are consistent with the chosen ranges of temperature, β, 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₂₅₀/F₃₅₀> 1.7), and sources with high redshift occupy the top left region of the diagram (F₂₅₀/F₃₅₀< 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.

Acknowledgments

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 Astrofisica 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

All Tables

All Figures

	Fig. 1 HeViCS image at 250 μm with detected sources (red points). The blue squares indicate regions with the full coverage.
In the text

	Fig. 2 V2 field region before (left panel) and after (right panel) background subtraction.
In the text

	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.
In the text

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.

In the text

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.

In the text

	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).
In the text

	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.
In the text

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

In the text

	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 availables, photometric redshift are considered otherwise.
In the text

	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.
In the text

	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).
In the text

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 β = 2, 10 ≤ T ≤ 30, and the redshift equals to the average of the bin. Rhombus, triangles, and stars indicate a modified black body with β = 2 and temperature of 10, 20, and 30 K, respectively.

In the text