Submillimeter number counts at 250 μm, 350 μm and 500 μm in BLAST data

M. Béthermin; H. Dole; M. Cousin; N. Bavouzet

doi:10.1051/0004-6361/200913910

Home

All issues

Volume 516 (June-July 2010)

A&A, 516 (2010) A43

Full HTML

Free Access

Issue		A&A Volume 516, June-July 2010


Article Number		A43
Number of page(s)		12
Section		Cosmology (including clusters of galaxies)
DOI		https://doi.org/10.1051/0004-6361/200913910
Published online		24 June 2010

A&A 516, A43 (2010)

Submillimeter number counts at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m in BLAST data

M. Béthermin - H. Dole - M. Cousin - N. Bavouzet

Institut d'Astrophysique Spatiale (IAS), Université Paris-Sud 11 and CNRS (UMR8617), Bât. 121, 91405 Orsay, France

Received 18 December 2009 / Accepted 6 March 2010

Abstract
Context. The instrument BLAST (Balloon-borne Large-Aperture Submillimeter Telescope) performed the first deep and wide extragalactic survey at 250, 350 and 500 $\mu$ m. The extragalactic number counts at these wavelengths are important constraints for modeling the evolution of infrared galaxies.
Aims. We estimate the extragalactic number counts in the BLAST data, which allow a comparison with the results of the P(D) analysis of Patanchon et al. (2009).
Methods. We use three methods to identify the submillimeter sources. 1) Blind extraction using an algorithm when the observed field is confusion-limited and another one when the observed field is instrumental-noise-limited. The photometry is computed with a new simple and quick point spread function (PSF) fitting routine (FASTPHOT). We use Monte-Carlo simulations (addition of artificial sources) to characterize the efficiency of this extraction, and correct the flux boosting and the Eddington bias. 2) Extraction using a prior. We use the Spitzer 24 $\mu$ m galaxies as a prior to probe slightly fainter submillimeter flux densities. 3) A stacking analysis of the Spitzer 24 $\mu$ m galaxies in the BLAST data to probe the peak of the differential submillimeter counts.
Results. With the blind extraction, we reach 97 mJy, 83 mJy and 76 mJy at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m respectively with a 95% completeness. With the prior extraction, we reach 76 mJy, 63 mJy, 49 mJy at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m respectively. With the stacking analysis, we reach 6.2 mJy, 5.2 mJy and 3.5 mJy at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m respectively. The differential submillimeter number counts are derived, and start showing a turnover at flux densities decreasing with increasing wavelength.
Conclusions. There is a very good agreement with the P(D) analysis of Patanchon et al. (2009). At bright fluxes (>100 mJy), the Lagache et al. (2004) and Le Borgne et al. (2009) models slightly overestimate the observed counts, but the data agree very well near the peak of the differential number counts. Models predict that the galaxy populations probed at the peak are likely $z\sim 1.8$ ultra-luminous infrared galaxies.

Key words: cosmology: observations - galaxies: statistics - galaxies: evolution - galaxies: photometry - infrared: galaxies

1 Introduction

Galaxy number counts, a measurement of the source surface density as a function of flux density, are used to evaluate the global evolutionary photometric properties of a population observed at a given wavelength. These photometric properties mainly depend on the source redshift distribution, spectral energy distribution (SED), and luminosity distribution in a degenerate way for a given wavelength. Even though this is a rather simple tool, measurements of number counts at different observed wavelengths greatly help in constraining those degeneracies. Backward evolution models, among these Pearson & Khan (2009); Le Borgne et al. (2009); Lagache et al. (2004); Gruppioni et al. (2005); Rowan-Robinson (2009); Franceschini et al. (2009); Valiante et al. (2009); Chary & Elbaz (2001) are able to broadly reproduce (with different degrees of accuracy) the observed number counts from the near-infrared to the millimeter spectral ranges, in addition to other current constraints, like such as measured luminosity functions and the spectral energy distribution of the Cosmic Infrared Background (CIB) (Kashlinsky 2005; Hauser et al. 1998; Hauser & Dwek 2001; Fixsen et al. 1998; Puget et al. 1996; Dole et al. 2006; Lagache et al. 1999,2005; Gispert et al. 2000). In the details, however, the models disagree in some aspects like the relative evolution of luminous and ultra-luminous infrared galaxies (LIRG and ULIRG) and their redshift distributions, or the mean temperature or colors of galaxies, as is shown for instance in LeFloc'h et al. (2009) from Spitzer 24 $\mu$ m deep observations.

One key spectral range lacks valuable data to get accurate constraints as yet: the sub-millimeter range, between 160 $\mu$ m and 850 $\mu$ m, where some surveys were conducted on small areas. Fortunately this spectral domain is intensively studied with the BLAST balloon experiment (Devlin et al. 2009) and the Herschel and Planck space telescopes. This range, although it is beyond the maximum of the CIB's SED in wavelength, allows us to constrain the poorly-known cold component of galaxy SED at a redshift greater than a few tenths. Pioneering works have measured the local luminosity function (Dunne et al. 2000) and shown that most milli-Jansky sources lie at redshifts z>2 (Pope et al. 2006; Ivison et al. 2005,2002; Chapman et al. 2005; Pope et al. 2005; Chapman et al. 2003a). Other works showed that the galaxies SED selected in the submillimeter range (Kovács et al. 2006; Benford et al. 1999; Chapman et al. 2003b; Sajina et al. 2003; Michaowski et al. 2010; Lewis et al. 2005; Sajina et al. 2006; Beelen et al. 2006) can have typically warmer temperatures and higher luminosities than galaxies selected at other infrared wavelengths.

Data in the submillimeter wavelength with increased sensitivity are thus needed to match the depth of infrared surveys, conducted by Spitzer in the mid- and far-infrared with the MIPS instrument (Rieke et al. 2004) at 24 $\mu$ m, 70 $\mu$ m and 160 $\mu$ m (Frayer et al. 2006a; Marleau et al. 2004; Frayer et al. 2009,2006b; Rodighiero et al. 2006; Chary et al. 2004; Béthermin et al. 2010; Shupe et al. 2008; LeFloc'h et al. 2009; Papovich et al. 2004; Dole et al. 2004) as well as the near-infrared range with the IRAC instrument (Fazio et al. 2004b) between 3.6 $\mu$ m and 8.0 $\mu$ m (Sullivan et al. 2007; Magdis et al. 2008; Ashby et al. 2009; Franceschini et al. 2006; Fazio et al. 2004a; Barmby et al. 2008). Infrared surveys have allowed the resolution of the CIB by identifying the contributing sources - directly at 24 $\mu$ m and 70 $\mu$ m, or indirectly trough stacking at 160 $\mu$ m (Béthermin et al. 2010; Dole et al. 2006).

Although large surveys cannot solve by themselves all the unknowns about the submillimeter SED of galaxies, the constraints given by the number counts can greatly help in unveiling the statistical SED shape of submillimeter galaxies as well as the origin of the submillimeter background.

The instrument BLAST (Balloon-borne Large-Aperture Submillimeter Telescope, Pascale et al. 2008) performed the first wide and deep survey in the 250-500 $\mu$ m range (Devlin et al. 2009) before the forthcoming Herschel results. Marsden et al. (2009) show that sources detected by Spitzer at 24 $\mu$ m emit the main part of the submillimeter background. Khan et al. (2009) claimed that only 20% of the CIB is resolved by the sources brighter than 17 mJy at 350 $\mu$ m. Patanchon et al. (2009) has performed a P(D) fluctuation analysis to determine the counts at BLAST wavelength (250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m). In this paper we propose another method to estimate the number counts at these wavelengths and compare the results with those of Patanchon et al. (2009).

2 Data

2.1 BLAST sub-millimeter public data in the Chandra Deep Field South (CDFS)

The BLAST holds a bolometer array, which is the precursor of the spectral and photometric imaging receiver (SPIRE) instrument on Herschel, at the focus of a 1.8 m diameter telescope. It observes at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, with a 36'', 42'' and 60'' beam, respectively (Truch et al. 2009).

An observation of the Chandra Deep Field South (CDFS) was performed during a long duration flight in Antarctica in 2006, and the data of the two surveys are now public: a 8.7 deg² shallow field and a 0.7 deg² confusion-limited (Dole et al. 2004) field in the center part of the first one. We use the non-beam-smoothed maps and associated point spread function (PSF) distributed on the BLAST website. The signal and noise maps were generated by the SANEPIC algorithm (Patanchon et al. 2008).

2.1 Spitzer 24 $\mu$ m data in the CDFS

Several infrared observations were performed in the CDFS. The Spitzer Wide-Field InfraRed Extragalactic (SWIRE) survey overlaps the CDFS BLAST field at wavelengths between 3.6 $\mu$ m and 160 $\mu$ m. We used only the 24 $\mu$ m band, which is 80% complete at 250 $\mu$ Jy. The completeness is defined as the probability to find a source of a given flux in a catalog. The Far-Infrared Deep Extragalactic Legacy (FIDEL) survey is deeper but narrower (about 0.25 deg²) than SWIRE and 80% complete at 57 $\mu$ Jy at 24 $\mu$ m. We used the Béthermin et al. (2010) catalogs constructed from these two surveys. These catalogs were extracted with SExtractor (Bertin & Arnouts 1996) and the photometry was performed with the allstar routine of the DAOPHOT package (Stetson 1987). The completeness of this catalog was characterized with Monte-Carlo simulations (artificial sources added on the initial map and extracted).

3 Blind source extraction and number counts

We started with a blind source extraction in the BLAST bands. Each wavelength was treated separately. For each wavelength we defined two masks: a shallow zone (about 8.2 deg²) covering the whole field except the noisier edge; and a deep zone (about 0.45 deg²) in the center of the confusion-limited area. We used different extraction methods in the shallow zone and the deep one, but the photometry and the corrections of the extraction bias were the same.

3.1 Detector noise-limited extraction (shallow zone)

In the shallow zone we used the non-smoothed map and the corresponding map of the standard deviation of the noise. The map was then cross-correlated by the PSF. The result of this cross-correlation is

$\begin{displaymath}m_{\rm conv}(i_{0},j_{0}) = \sum_{i=-N}^{+N} \sum_{j=-N}^{+N} m(i_0+i,j_0+j) \times {{\it PSF}}(i,j), \end{displaymath}$

(1)

where $m_{\rm conv}(i_{0},j_{0})$ is the flux density in the pixel (i₀,j₀) of the cross-correlated map, m(i,j) the flux density in the pixel (i,j) of the map, and PSF(i,j) the value of the normalized PSF in the pixel (i,j) (the center of the PSF is in the center of the pixel (0,0)). The PSF size is (2N+1) $\times$ (2N+1) pixels. The standard deviation of the noise in the cross-correlated map is thus

$\begin{displaymath}n_{\rm conv}(i_{0},j_{0}) = \sqrt{\sum_{i=-N}^{+N} \sum_{j=-N}^{+N} n^2(i_0+i,j_0+j) \times {{\it PSF}}^2(i,j)}, \end{displaymath}$

(2)

where n ( $n_{\rm conv}$ ) is the initial (cross-correlated) map of the standard deviation of the noise.

We found the pixels where $m_{\rm conv}/n_{\rm conv}>3$ and kept the local maxima. The precise center of the detected sources was computed by a centroid algorithm. This low threshold caused lots of spurious detections, but helped to deblend the fluxes of 3 to 4-sigma sources and avoided to overestimate their fluxes. We could thus limit the flux boosting effect. A final cut in flux after the PSF fitting photometry eliminated the main part of these sources. We performed the extraction algorithm on the flipped map (initial map multiplied by a factor of -1) to check it. We found few spurious sources brighter than the final cut in flux determined in the Sect. 3.4. We found a spurious rate of 12%, 11% and 25% at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively.

3.2 Confusion-limited extraction (deep zone)

In the confusion-limited zone we also used a non-smoothed map. In this region the noise is dominated by the confusion and not by the instrumental noise. Consequently, the method based on instrumental noise presented in the Sect. 3.1 is not relevant. We used an atrou wavelet filtering (Starck et al. 1999; Dole et al. 2001) to remove fluctuations at scales larger than 150''. Then we divided the resulting map by $\sigma_{\rm filtered~map}$ , which is the standard deviation of the pixel values on the filtered map in the working area. We finally kept local maxima with a signal greater than 3. The center of the sources was also determined by a centroid algorithm. The initial map and the cleaned map are shown in Fig. 1. When we flip the map, we find no spurious source brighter than the final cut in flux determined in Sect. 3.4.

$\begin{figure} \par\includegraphics[width=8cm,clip]{13910f1.eps} \end{figure}$	Figure 1: Position of sources brighter than the 95% completeness flux at 250 $\mu$ m in deep zone. Top: initial map. Bottom: residual map. The area out of the mask are represented darker. These 1 $^\circ \times$ 1 $^\circ$ map are centered on the coordinates (RA, Dec) = (3h32min30s, -27 $^\circ$ 50'). The horizontal axis is aligned with the right ascension.
Open with DEXTER

3.3 A simple and quick PSF fitting routine: FASTPHOT

For both noise- and confusion-limited extraction, we apply the same quick and simple PSF fitting routine on the non-beam-smoothed map. This routine fits all the detected sources at the same time and is consequently efficient for deblending (although no source was blended in this case; but source-blending will be an issue for an extraction using a prior, detailed in Sect. 4. We suppose that the noise is Gaussian and the position of sources is known. We then maximize the likelihood

$\begin{displaymath}L(m\vert S) = \prod_{\rm pixels} C(n) \times \exp \left[ - \f... ...} {{\it PSF}}_{x_i, y_i} \times S_i \right)^2}{2 n^2} \right], \end{displaymath}$

(3)

where m and n are the map and the noise map. PSF _{x_i,
y_i} is a unit-flux PSF centered at the position (x_i, y_i), which are the coordinates of the ith source. These coordinates are not necessarily integers. C(n) is a normalization constant and depends only of the value of the noise map. ${\vec S}$ is a vector containing the flux of the sources.

The value of S, which maximizes the likelihood, satisfies the following linear equation stating that the derivative of the likelihood logarithm equals zero

$\begin{displaymath}\forall i, 0 = \frac{\partial{\log \bigl(L(m\vert S)\bigl)}}{\partial S_i} = A.S + B, \end{displaymath}$

(4)

where A is a matrix and B a vector defined by

$\begin{displaymath}A = (a_{ij}) = - \sum_{\rm pixels} \frac{{{\it PSF}}_{x_i,y_i} \times {{\it PSF}}_{x_j,y_j}}{n^2} \end{displaymath}$

(5)

$\begin{displaymath}B = (b_i) = \sum_{\rm pixels} \frac{{{\it PSF}}_{x_i,y_i} \times map}{n^2}\cdot \end{displaymath}$

(6)

To perform this operation fast, we used a $70'' \times 70''$ (respectively $90'' \times 90''$ and $110'' \times 110''$ ) PSF at 250 $\mu$ m (respectively 350 $\mu$ m and 500 $\mu$ m). This PSF, provided by the BLAST team, is the response for a unit-flux source and takes into account all the filtering effects. We used the conjugate gradient method to solve the Eq. (4) quickly.

This routine was tested with $200 \times 200$ pixels simulated maps containing 400 sources at a known positions with a beam of 10 pixels FWHM. The flux of all sources was perfectly recovered in the case where no noise was added. This routine (FASTPHOT) performs simultaneous PSF fitting photometry of 1000 sources in less than 1 s. It is publicly available.

Table 1: 95% completeness flux density and photometric noise for different depths at different wavelengths.

3.4 Completeness and photometric accuracy

The completeness is the probability to detect a source of a given flux density. We measured it with a Monte-Carlo simulation. We added artificial point sources (based on PSF) on the initial map at random positions and performed the same source extraction and photometry algorithm as for the real data. A source was considered to be detected if there was a detection in a 20'' radius around the center of the source. Table 1 gives the 95% completeness flux density (for which 95% of sources at this flux are detected) for different wavelengths and depths.

The photometric noise was estimated with the scatter of the recovered fluxes of artificial sources. We computed the standard deviation of the difference between input and output flux. This measurement includes instrumental and confusion noise ( $\sigma_{\rm tot} = \sqrt{\sigma_{\rm instr}^2 + \sigma_{\rm conf}^2}$ ). The results are given in Table 1. In the deeparea, the photometric uncertainties are thus dominated by the confusion noise. The estimations of the confusion noise between the deep and shallow areas are consistent. It shows the accuracy and the consistency of our method.

Note that the uncertainties on flux densities in the Dye et al. (2009) catalog (based only on instrumental noise) are consequently largely underestimated in the confusion-limited area. Indeed, their 5 $\sigma$ detection threshold (based only on instrumental noise) at 500 $\mu$ m in the deep zone corresponds to 1.76 $\sigma$ if we also include the confusion noise.

$\begin{figure} \par\includegraphics[width=8.6cm,clip]{13910f2.eps} \end{figure}$

Figure 2:

Extragalactic number counts at 250 $\mu$ m in the BLAST data. Diamond: counts deduced from the source catalog on the whole shallow field; square: counts deduced from the catalog of the deep part of the field; triangle: counts deduced from catalog of the deep part of the field with a 24 $\mu$ m prior (this measurement gives only a lower limits to the counts); cross: counts computed with a stacking analysis; grey asterisk: counts computed with a P(D) analysis by Patanchon et al. (2009); grey short dashed line: Lagache et al. (2004) model prediction; grey long dashed line and grey area: Le Borgne et al. (2009) model prediction and 1- $\sigma$ confidence area.

Open with DEXTER

The faint flux densities are overestimated due to the classical flux boosting effect. This bias was measured for all bands for 60 flux densities between 10 mJy and 3 Jy with the results of the Monte-Carlo simulations. The measured fluxes were deboosted with this relation. We cut the catalogs at the 95% completeness flux, where the boosting factor is at the order of 10%. Below this cut, the boosting effect increases too quickly to be safely corrected. We also observed a little underestimation at high flux of 1%, 0.5% and 0.5% at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m. It is due to FASTPHOT, which assumes that the position is perfectly known, which is not true, especially for a blind extraction.

3.5 Number counts

We computed number counts with catalogs corrected for boosting. For each flux density bin we subtracted the number of spurious detections estimated in the Sects. 3.1 and 3.2 to the number of detected sources and divided the number of sources by the size of the bin, the size of the field and the completeness.

We also applied a corrective factor for the Eddington bias. We assumed a distribution of flux densities in ${\rm d}N/{\rm d}S \propto~S^{-r}$ with $r=3\pm0.5$ . This range of possible values for r was estimated considering the Patanchon et al. (2009) counts and the Lagache et al. (2004) and Le Borgne et al. (2009) model predictions. We then randomly kept sources with a probability given by the completeness and added a random Gaussian noise to simulate photometric noise. Finally we computed the ratio between the input and output number of sources in each bin. We applied a correction computed for r = 3 to each point. We estimated the uncertainty on this correction with the difference between corrections computed for r = 2.5 and r = 3.5. This uncertainty was quadratically combined with a Poissonian uncertainty (clustering effects are negligible due to the little number of sources in the map, see Appendix A).

The calibration uncertainty of BLAST is 10%, 12% and 13% at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m respectively (Truch et al. 2009). This uncertainty is combined with other uncertainties on the counts. The results are plotted in Fig. 2 and given in Table 2 and interpreted in Sect. 6.

3.6 Validation

We used simulations to validate our method. We generated 50 mock catalogs based on the Patanchon et al. (2009) counts, and which covered 1 deg² each. These sources are spatially homogeneously distributed. We then generated the associated maps at 250 $\mu$ m. We used the instrumental PSF, and added a Gaussian noise with the same standard deviation as in the deepest part of real map.

We performed an extraction of sources and computed the number counts with the method used in the confusion limited part of the field (Sect. 3.2). We then compared the output counts with the initial counts (Fig. 3). We used two flux density bins: 100-141 mJy and 141-200 mJy. We found no significant bias. The correlation between the two bins is 0.46. The neighbor points are thus not anti-correlated as in the Patanchon et al. (2009) P(D) analysis.

The same verification was done on 20 Fernandez-Conde et al. (2008) simulations (based on the Lagache et al. 2004 model). These simulations include clustering. This model overestimates the number of the bright sources, and the confusion noise is thus stronger. The 95% completeness is then reach at 200 mJy. But there is also a very good agreement between input and output counts in bins brighter than 200 mJy. We found a correlation between two first bins of 0.27.

4 Source extraction using Spitzer 24 $\mu$ m catalog as a prior

In addition to blind source extraction in the BLAST data (Sect. 3) we also performed a source extraction using a prior.

4.1 PSF fitting photometry at the position of the Spitzer 24 $\mu$ m

The catalogs of infrared galaxies detected by Spitzer contain more sources than the BLAST catalog. The 24 $\mu$ m Spitzer PSF has a Full Width at Half Maximum (FWHM) of 6.6''. It is smaller than the BLAST PSF (36'' at 250 $\mu$ m). Consequently, the position of the Spitzer sources is known with sufficient accuracy when correlating with the BLAST data.

We applied the FASTPHOT routine (Sect. 3.3) at the positions of 24 $\mu$ m sources. We used the Béthermin et al. (2010) SWIRE catalog cut at $S_{24}=250~\mu$ Jy (80% completeness). In order to avoid software instabilities, we kept in our analysis only the brightest Spitzer source in a 20'' radius area (corresponding to 2 BLAST pixels). The corresponding surface density is 0.38, 0.49 and 0.89 Spitzer source per beam at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively.

This method works only if there is no astrometrical offset between the input 24 $\mu$ m catalog and the BLAST map. We stacked the BLAST sub-map centered on the brightest sources of the SWIRE catalog and measured the centroid of the resulting artificial source. We found an offset of less than 1''. It is negligible compared to the PSF FWHM (36'' at 250 $\mu$ m).

We worked only in the central region of the deep confusion-limited field (same mask as for blind extraction), where the photometric noise is low.

4.2 Relevance of using Spitzer 24 $\mu$ m catalog as a prior

The S₂₅₀/S₂₄ ( S₃₅₀/S₂₄ or S₅₀₀/S₂₄) color is not constant, and some sources with a high color ratio could have been missed in the prior catalog (especially high-redshift starbursts). We used the Lagache et al. (2004) and Le Borgne et al. (2009) models to estimate the fraction of sources missed. We selected the sources in the sub-mm flux density bin and computed the 24 $\mu$ m flux density distribution (see Fig. 4). According to the Lagache et al. (2004) model, 99.6%, 96.4% and 96.9% of the sub-mm selected sources⁴ are brighter than $S_{24}=250~\mu$ Jy for a selection at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively The Le Borgne et al. (2009) model gives 99.8%, 98.3% and 95.0%, respectively.

Table 2: Number counts deduced from source extraction. The not normalized counts can be obtained dividing the S^2.5.dN/dS column by the $S_{\rm mean}$ column.

4.3 Photometric accuracy

The photometric accuracy was estimated with Monte-Carlo artificial sources. We added five sources of a given flux at random positions on the original map and add them to the 24 $\mu$ m catalog. We then performed a PSF fitting and compared the input and output flux. We did this 100 times per tested flux for 10 flux densities (between 10 and 100 mJy). In this simulation we assumed that the position of the sources is exactly known. It is a reasonable hypothesis due to the 24 $\mu$ m PSF FWHM (6.6'') compared to the BLAST one (36'' at 250 $\mu$ m).

We did not detect any boosting effect for faint flux densities as expected in this case of detection using a prior. For a blind extraction there is a bias of selection toward sources located on peaks of (instrumental or confusion) noise. This is not the case for an extraction using a prior, for which the selection is performed at another wavelength.

The scatter of output flux densities is the same for all the input flux densities. We found a photometric noise $\sigma_S$ of 21.5 mJy, 18.3 mJy and 16.6 mJy at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively. It is slightly lower than for the blind extraction, for which the position of source is not initially known.

4.4 Estimation of the number counts

From the catalog described in Sect. 4.1 we give an estimation of the submillimeter number counts at flux densities fainter than reached by the blind-extracted catalog. We cut the prior catalog at 3 $\sigma_S$ , corresponding to 64 mJy, 54 mJy and 49 mJy at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively. We worked in a single flux density bin, which is defined to be between this value and the cut of the blind-extracted catalog. There is no flux boosting effect, but we needed to correct the Eddington bias. The completeness could not be defined in the same way as for the blind extraction, because the selection was performed at another wavelength. We thus cannot suppose power-law counts, because the selection function is unknown and the distribution of the extracted sources cannot be computed.

$\begin{figure} \par\includegraphics[width=8cm,clip]{13910f3.eps} \end{figure}$	Figure 3: Number counts at 250 $\mu$ m deduced from a blind extraction for 50 realizations of a simulation based on the Patanchon et al. (2009) counts. The horizontal solid line represents the input count value. The lower panel is the result of the 100-141 mJy bin, the upper panel is the 141-200 mJy bin.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8cm,clip]{13910f4.eps} \end{figure}$	Figure 4: Flux density distribution at 24 $\mu$ m of the 54 mJy <S₃₅₀< 83 mJy sources. The Lagache et al. (2004) model is plotted in black and the Le Borgne et al. (2009) model is plotted in grey. The dashed line represents the cut of our catalog at 24 $\mu$ m.
Open with DEXTER

The Eddington bias was estimated with another method. We took the sub-mm flux of each of the sources selected at 24 $\mu$ m and computed how many sources lie in our count bin. We added a Gaussian noise $\sigma_S$ to the flux of each source to simulate the photometric errors. We computed the number of sources in the counts bin for the new fluxes. We then compute the mean of the ratio between the input and output number of sources in the selected bin for 1000 realizations. The estimated ratios are 0.42, 0.35 and 0.21 at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively. These low values indicate that on average the photometric noise introduces an excess of faint sources in our flux bin. This effect is strong because of the steep slope of the number counts, implying more fainter sources than brighter sources. The results are interpreted in the Sect. 6.

4.5 Sub-mm/24 color

In this part we work only on $S > 5 \sigma_S$ sources of the catalog described in Sect. 4.1 to avoid bias due to the Eddington bias in our selection. At 250 $\mu$ m, we have two sources verifying this criterion with a S₂₅₀/S₂₄ color of 16 and 60. No sources are brighter than 5 $\sigma_S$ at larger wavelengths. For this cut in flux ( $S_{250}> 5\sigma_S$ ), the Lagache et al. (2004) and Le Borgne et al. (2009) models predict a mean S₂₅₀/S₂₄ color of 39 and 41, respectively. The two models predict a mean redshift of 0.8 for this selection, and the K-correction effect explains these high colors.

5 Non-resolved source counts by stacking analysis

5.1 Method

In order to probe the non-resolved source counts, we used same method as Béthermin et al. (2010), i.e. the stacking analysis applied to number counts (hereafter ``stacking counts''). We first measured the mean flux at 250 $\mu$ m, 350 $\mu$ m or 500 $\mu$ m as a function of the 24 $\mu$ m flux ( $\overline{S_{250,~350~{\rm or}~500}} = f(S_{24})$ ). This measurement was performed by stacking in several S₂₄ bins. We used the Béthermin et al. (2010) catalog at 24 $\mu$ m of the FIDEL survey. It is deeper than the SWIRE one used in Sect. 4, but covers a smaller area (0.25 deg²). The photometry of stacked images was performed with the PSF fitting method (Sect. 3.3), and the uncertainties on the mean flux are computed with a bootstrap method (Bavouzet 2008). We then computed the counts in the sub-mm domain with the following formula:

$\begin{displaymath}\frac{{\rm d}N}{{\rm d}S_{\rm submm}}\biggl \vert _{S_{\rm su... ...\rm d}S_{\rm submm}}{{\rm d}S_{24}}\biggl \vert _{S_{24}}\cdot \end{displaymath}$

(7)

We show in Appendix B that the clustering effect can be neglected. The results are given in Table 3 and are plotted in Fig. 2.

Table 3: Number counts deduced from stacking.

5.2 Validity of the stacking analysis in the sub-mm range

There are 1.8, 2.4 and 4.5 S₂₄>70 $\mu$ Jy sources per BLAST beam at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively. We thus stacked several sources per beam. Béthermin et al. (2010) showed that the stacking analysis is valid at 160 $\mu$ m in the Spitzer data, where the size of the beam is similar to the BLAST one.

To test the validity of the stacking analysis in the BLAST data from a Spitzer 24 $\mu$ m catalog, we generated a simulation of a 0.25 deg² with a Gaussian noise at the same level as for the real map and with source clustering, following Fernandez-Conde et al. (2008). We stacked the 24 $\mu$ m simulated sources per flux bin in the BLAST simulated maps. We measured the mean BLAST flux for each 24 $\mu$ m bin with the same method as applied on the real data. At the same time we computed the mean sub-mm flux for the same selection from the mock catalog associated to the simulation. We finally compared the mean BLAST fluxes measured by stacking with those directly derived from the mock catalog to estimate the possible biases (see Fig. 5). The stacking measurements and expected values agree within the error bars. We notice a weak trend of overestimation of the stacked fluxes at low flux density (S₂₄< 200 $\mu$ Jy) however, but it is still within the error bars. We can thus stack 24 $\mu$ m Spitzer sources in the BLAST map.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13910f5.eps} \end{figure}$	Figure 5: Ratio between the mean flux density at 250 $\mu$ m found by the stacking analysis and the expected flux for different S₂₄ bins. It is based on a Fernandez-Conde et al. (2008) simulation of a 0.25 deg² field with a noise and PSF similar to the BLAST deep region.
Open with DEXTER

5.3 Mean 24 $\mu$ m to sub-mm color deduced by stacking analysis

The stacking analysis allowed to measure the mean 24 $\mu$ m to sub-mm colors of undetected sub-mm galaxies. These colors depends on the SED of galaxies (or K-correction) and the redshift distribution in a degenerate way. The $S_{\rm submm}/S_{24}$ color and ${\rm d}S_{\rm submm}/{\rm d}S_{24}$ as a function of S₂₄ are plotted in Fig. 6.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13910f6.eps} \end{figure}$	Figure 6: Black solid line: ${\rm d}S_{\rm submm}/{\rm d}S_{24}$ as a function of S₂₄. Grey dashed line: $S_{\rm submm}/S_{24}$ color as a function of S₂₄.
Open with DEXTER

The colors are higher for the fainter 24 $\mu$ m flux ( $S_{24} < 100~\mu$ Jy). This behavior agrees with the model expectations: the faint sources at 24 $\mu$ m lie at a higher mean redshift than the brighter ones. Due to the K-correction, the high-redshift sources have a brighter sub-mm/24 color than local ones.

The colors found by the stacking analysis are lower than those obtained by an extraction at 250 $\mu$ m (Sect. 4.5). It is an effect of selection. The mid-infrared is less affected by the K-correction than the sub-mm, and a selection at this wavelength selects lower redshift objects. We thus see lower colors because of the position of the SED peak (around 100 $\mu$ m rest-frame).

We also investigated the evolution of the derivative d $S_{\rm submm}/{\rm d}S_{24}$ as a function of S₂₄, which explicits how the observed sub-mm flux increases with the 24 $\mu$ m flux densities. At high 24 $\mu$ m flux densities (S₂₄> 400 $\mu$ Jy) the derivative is almost constant and small (<20 and compatible with zero), meaning that the observed sub-mm flux density does not vary much with S₂₄. For these flux bins we select only local sources and do not expect a strong evolution of the color. At fainter 24 $\mu$ m flux densities the observed decrease can be explained by redshift and K-correction effects, as above.

The color in the faintest 24 $\mu$ m flux density bin (70 to 102 $\mu$ Jy) is slightly fainter than in the neighboring points. It can be due to the slight incompleteness of the 24 $\mu$ m catalog (about 15%), which varies spatially across the field: the sources close to the brightest sources at 24 $\mu$ m are hardly extracted. The consequence is a bias to the lower surface density regions, leading to a slight underestimation of the stacked flux measurement.

5.4 Accuracy of the stacking counts method on BLAST with a Spitzer 24 $\mu$ m prior

Béthermin et al. (2010) showed that the stacking counts could be biased: the color of sources can vary a lot as a function of the redshift. The assumption of a single color for a given S₂₄ is not totally realistic and explains some biases. We used two simulated catalogs (containing for each source S₂₄, S₂₅₀, S₃₅₀ and S₅₀₀) to estimate this effect: a first one based on the Lagache et al. (2004) model that covered $20 \times 2.9$ deg² and a second one based on the Le Borgne et al. (2009) model and that covered 10 deg². The large size of these simulations allows us to neglect cosmic variance.

In order to compute the stacking counts, we first computed the counts at 24 $\mu$ m from the mock catalog. Then we computed the mean S_{250, 350 or 500} flux density (directly in the catalog) in several S₂₄ bins to simulate a stacking. We finally applied the Eq. (7) to compute stacking counts at the BLAST wavelengths.

The ratio between the stacking counts and the initial counts is plotted in Fig. 8 for the two mock catalogs. Between 1 mJy and 10 mJy we observe an oscillating bias. This bias is less than 30% at 250 $\mu$ m and 50% at other wavelengths. When the flux becomes brighter than 25 mJy at 250 $\mu$ m (18 mJy at 350 $\mu$ m and 7.5 mJy at 350 $\mu$ m), we begin to strongly underestimate the counts. The analysis of real data also shows a very strong decrease in the counts around the same fluxes (see Fig. 7). Consequently, we cut our stacking analysis at these fluxes and we applied an additional uncertainty to the stacking counts of 30% at 250 $\mu$ m (50% at 350 $\mu$ m and 500 $\mu$ m).

Using the 24 $\mu$ m observations as a prior to stack in the BLAST bands seems to give less accurate results than in the Spitzer MIPS bands. For a given S₂₄ flux, the sub-mm emission can vary a lot as a function of the redshift. But the simulations shows that this method works for faint flux densities. It is due to the redshift selection which is similar for faint flux densities (see Fig. 9) and very different at higher flux densities (see Fig. 10). For example, $S_{24} \sim 100~\mu$ Jy sources are distributed around z=1.5 with a broad dispersion in redshift. $S_{350} \sim 4$ mJy (based on averaged colors, 4 mJy at 350 $\mu$ m corresponds to $S_{24} \sim 100~\mu$ Jy) sources have quite a similar redshift distribution except an excess for z > 2.6. At higher flux densities (around 2 mJy at 24 $\mu$ m) the distribution is very different. The majority of the 24 $\mu$ m-selected sources lies at z < 1 and the distribution of 350 $\mu$ m-selected sources peaks at $z \sim 1.5$ . Another possible explanation is that fainter sources lies near z=1 and are thus selected at the 12 $\mu$ m rest-frame, which is a very good estimator of the infrared bolometric luminosity according to Spinoglio et al. (1995).

In order to limit the scatter of the sub-mm/24 color, we tried to cut our sample into two redshift boxes following the Devlin et al. (2009) IRAC color criterion ([3.6]-[4.5] = 0.068([5.8]-[8.0])-0.075). But we had not enough signal in the stacked images to perform the analysis.

$\begin{figure} \par\includegraphics[width=14cm,clip]{13910f7.eps} \end{figure}$

Figure 7:

Number counts at BLAST wavelengths coming from the data (points) and the models (lines). Short dashed line: initial (black) and stacking (grey) counts from the Lagache et al. (2004) mock catalog; long dashed line: initial (black) and stacking (grey) counts from Le Borgne et al. (2009); diamond: stacking counts built with the FIDEL catalog; square: stacking counts built with the SWIRE catalog; grey vertical dot-dash line: flux cut for stacking counts.

Open with DEXTER

$\begin{figure} \par\includegraphics[width=14cm,clip]{13910f8.eps} \end{figure}$	Figure 8: Ratio between stacking counts and initial counts for two mock catalogs; short dashed line: Lagache et al. (2004) catalog; long dashed line: Le Borgne et al. (2009); grey vertical dot-dash line: flux cut for stacking counts; horizontal dot line: estimation of the uncertainty intrinsic to the stacking method.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8cm,clip]{13910f9.eps} \end{figure}$	Figure 9: Solid line: distribution in redshift of the sources with 102 $\mu$ Jy <S₂₄< 150 $\mu$ Jy for the mock catalogs generated with the Lagache et al. (2004) (black) and the Le Borgne et al. (2009) (grey) models; dashed line: distribution in redshift of the sources with 4 mJy <S₂₅₀< 6 mJy (determined using the mean 250/24 color).
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8cm,clip]{13910f10.eps} \end{figure}$	Figure 10: Solid line: distribution in redshift of the sources with 2213 $\mu$ Jy <S₂₄< 3249 $\mu$ Jy for the mock catalogs generated with the Lagache et al. (2004) (black) and Le Borgne et al. (2009) (grey) models; dashed line: distribution in redshift of sources with 47 mJy <S₂₅₀< 62 mJy (determined using the mean 250/24 color).
Open with DEXTER

6 Interpretation

6.1 Contribution to the CIB

We integrated our counts assuming power-law behavior between our points. Our points are not independent (especially the stacking counts), and we thus combined errors linearly. The contribution of the individually detected sources ( S₂₅₀ > 64 mJy, S₃₅₀ > 54 mJy, S₅₀₀ > 49 mJy) is then 0.24 ^+0.18_-0.13 nW.m².sr^-1, 0.06 ^+0.05_-0.04 nW.m².sr^-1 and 0.01 ^+0.01_-0.01 nW.m².sr^-1 at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively. Considering the total CIB level of Fixsen et al. (1998) (FIRAS absolute measurement), we resolved directly only 2.3%, 1.1% and 0.4% at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively.

The populations probed by the stacking counts ( S₂₅₀ > 6.2 mJy, S₃₅₀ > 5.2 mJy, S ₅₀₀ > 3.5 mJy) emit 5.0 ^+2.5_-2.6 nW m² sr^-1, 2.8 ^+1.8_-2.0 nW m² sr^-1 and 1.4 ^+2.1_-1.3 nW m² sr^-1 at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively. This corresponds to about 50% of the CIB at these three wavelengths.

6.2 Comparison with Patanchon et al. (2009)

The agreement between our resolved counts built from the catalogs and the P(D) analysis of Patanchon et al. (2009) is excellent (Fig. 2). We confirm the efficiency of the P(D) analysis to recover number counts without extracting sources. The stacking counts probe the flux densities between 6 mJy and 25 mJy at 250 $\mu$ m (between 5 mJy and 13 mJy at 350 $\mu$ m and 3 mJy and 7 mJy at 500 $\mu$ m). In this range there is only one P(D) point. At the three BLAST wavelengths the P(D) points agree with our stacking counts (Fig. 2). Our results thus confirm the measurement of Patanchon et al. (2009) and give a better sampling in flux.

6.3 Comparison with ground-based observations

We compared our results with sub-mm ground-based observations of SHARC. Khan et al. (2007) estimated a density of S₃₅₀ > 13 mJy sources of 0.84_-0.61^+1.39 arcmin^-2. For the same cut, we found $0.26 \pm 0.13$ arcmin^-1, which agrees with their work. Our measurement ( $175 \pm 75$ sources deg^-2 brighter than 25 mJy) agrees also with that of Coppin et al. (2008) ones at the same wavelength (200-500 sources deg^-2 brighter than 25 mJy).

We also compared our results at 500 $\mu$ m with the SCUBA ones at 450 $\mu$ m. Borys et al. (2003) find 140 _-90⁺¹⁴⁰ gal deg^-2for S₄₅₀ > 100 mJy. We found $1.2 \pm 1.0$ gal.deg^-2. We significantly disagree with them. Borys et al. (2003) claim 5 4- $\sigma$ detections in a 0.046 deg² field in the Hubble deep field north (HDFN). These five sources are brighter than 100 mJy. We find no source brighter than 100 mJy in a 0.45 deg² field at 350 $\mu$ m nor at 500 $\mu$ m. The cosmic variance alone thus cannot explain this difference. A possible explanation is that they underestimated the noise level and their detections are dominated by spurious sources. It could also be due to a calibration shift (by more than a factor 2). The observation of the HDFN by Herschel will allow us to determine whether that these bright sources might be spurious detections.

We also compared our results with the estimations based on lensed sources at 450 $\mu$ m with SCUBA (Knudsen et al. 2006; Smail et al. 2002). For example, Knudsen et al. (2006) find 2000-50 000 sources deg² brighter than 6 mJy. It agrees with our 3500 ⁺⁷⁷⁰⁰_-3400 sources deg².

6.4 Comparison with the Lagache et al. (2004) and Le Borgne et al. (2009) models

At 250 $\mu$ m and 350 $\mu$ m the measured resolved source counts are significantly lower (by about a factor of 2) than the Lagache et al. (2004) and Le Borgne et al. (2009) models. Nevertheless, our counts are within the confidence area of Le Borgne et al. (2009). The same effect (models overestimating the counts) was observed at 160 $\mu$ m (Béthermin et al. 2010; Frayer et al. 2009). It indicates that the galaxies' SED or the luminosity functions used in both models might have to be revisited. At 500 $\mu$ m our counts and both models agree very well, but our uncertainties are large, which renders any discrimination difficult.

Concerning the stacking counts, they agree very well with the two models. Nevertheless, our uncertainties are larger than 30%. We thus cannot check if the disagreement observed between the Lagache et al. (2004) model and the stacking counts at 160 $\mu$ m (Béthermin et al. 2010) of 30% at S₁₆₀ =20 mJy still holds at 250 $\mu$ m.

6.5 Implications for the probed populations and the models

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13910f11.eps} \end{figure}$	Figure 11: Mean redshift (solid line) of sources for different fluxes at 350 $\mu$ m for the Lagache et al. (2004) (black) and Le Borgne et al. (2009) (grey) models and corresponding infrared luminosity defined in Sect. 6.5 (dashed line).
Open with DEXTER

We showed that the two models nicely reproduce the sub-mm counts, especially below 100 mJy. We can thus use them to estimate which populations are constrained by our counts. For each flux density bin we computed the mean redshift of the selected galaxies in both models. We then used the SEDs given by the models at that mean redshift and at that flux bin and derived the infrared bolometric luminosity. The luminosities are shown in dashed lines in Fig. 11 for 350 $\mu$ m as an example, and the redshift is given in solid lines.

The stacking counts reach 6.2 mJy, 5.3 mJy and 3.5 mJy at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively. This corresponds to faint ULIRGs ( $L_{\rm IR} \approx 1.5 \times 10^{12}~L_{\hbox{$\odot$ }}$ ) around z = 1.5, 1.8 and 2.1 at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively. Our measurements show that the predicted cold-dust emissions (between 100 $\mu$ m and 200 $\mu$ m rest frame) of this population in the models are believable.

At 250 $\mu$ m and 350 $\mu$ m the resolved sources ( S₃₅₀>85 mJy) are essentially $z \sim 1$ ULIRGs ( $L_{\rm IR} > 10^{12} L_{\hbox{$\odot$ }}$ ) and HyLIRGs ( $L_{\rm IR} > 10^{13}~L_{\hbox{$\odot$ }}$ ) according to the models. In Lagache et al. (2004) the local cold-dust sources contribute at very bright flux (>200 mJy). This population is not present in the Le Borgne et al. (2009) model. It explains the difference between the two models for fluxes brighter than 100 mJy at 350 $\mu$ m (see Fig. 11). At 500 $\mu$ m, Lagache et al. (2004) predict that bright counts are dominated by local cold-dust populations and Le Borgne et al. (2009) that they are dominated by medium redshift HyLIRGs. Nevertheless, there is a disagreement with the observations for this flux density range, suggesting that there could be less HyLIRGs than predicted. But these models do not currently include any AGN contribution, which is small except at luminosities higher than 10 $^{12}~L_\odot$ (Lacy et al. 2004; Valiante et al. 2009; Daddi et al. 2007).

7 Conclusion

Our analysis provides new stacking counts, which can be compared with the Patanchon et al. (2009) P(D) analysis. We have a good agreement between the different methods. Nevertheless, some methods are more efficient in a given flux range.

The blind extraction and the extraction using a prior give a better sampling in flux and slightly smaller error bars. The P(D) analysis uses only the pixel histogram and thus looses the information on the shape of the sources. The blind extraction is a very efficient method for extracting the sources, but lots of corrections must be applied carefully. When the confusion noise totally dominates the instrumental noise, the former must be determined accurately, and the catalog flux limit must take this noise (Dole et al. 2003) into account.

Estimating the counts from a catalog built using a prior is a good way to deal with the flux boosting effect. This method is based on assumptions however. We assume that all sources brighter than the flux cut at the studied wavelength are present in the catalog extracted using a prior. We also assume a flux distribution at the studied wavelength for a selection at the prior wavelength to correct for the Eddington bias. Consequently an extraction using a prior must be used in a flux range where the blind extraction is too affected by the flux boosting to be accurately corrected.

P(D) analysis and stacking counts estimate the counts at flux densities below the detection limit. These methods have different advantages. The P(D) analysis fits all the fluxes at the same time, where the stacking analysis flux depth depends on the prior catalog's depth (24 $\mu$ m Spitzer for example). But the P(D) analysis with a broken power-law model is dependent on the number and the positions of the flux nodes. The uncertainty due to the parameterization was not evaluated by Patanchon et al. (2009). The stacking counts on the other hand are affected by biases due to the color dispersion of the sources. The more the prior and stacked wavelength are correlated, the less biased are the counts. A way to overcome this bias would be to use a selection of sources (in redshift slices for example), which would reduce the color dispersion, and the induced bias; we did not use this approach here because of a low signal-to-noise ratio.

The stacking and P(D) analysis are both affected by the clustering in different ways. For the stacking analysis this effect depends on the size of the PSF. This effect is small for BLAST and will be smaller for SPIRE. The clustering broadens the pixel histogram. Patanchon et al. (2009) show that it is negligible for BLAST. Clustering will probably be an issue for SPIRE. The cirrus can also affect the P(D) analysis and broaden the peak. Patanchon et al. (2009) use a high-pass filtering that reduces the influence of these large scale structures.

The methods used in this paper will probably be useful to perform the analysis of the Herschel SPIRE data. The very high sensitivity and the large area covered will reduce the uncertainties and increase the depth of the resolved source counts. Nevertheless, according to the models (e.g. Le Borgne et al. (2009)), the data will also be quickly confusion-limited and it will be very hard to directly probe the break of the counts. The P(D) analysis of the deepest SPIRE fields will allow us to constrain a model with more flux nodes and to better sample the peak of the normalized differential number counts. The instrumental and confusion noise will be lower, and a stacking analysis per redshift slice will probably be possible. These analyses will give stringent constraints on the model of galaxies and finally on the evolution of the infrared galaxies.

Acknowledgements

We warmly acknowledge Guillaume Patanchon for his precious comments and discussions. We thank Damien Le Borgne and Guilaine Lagache for distributing their model and their comments. We also acknowledge Alexandre Beelen and Emeric Le Floc'h for their useful comments. We thank Maxime Follin, for his help during his Licence 3 training at the Université Paris Sud 11. We thank the BLAST team for the well-documented public release of their data. We warmly thank the referee Steve Willner, who helped a lot to improve the quality of this paper. This work is based in part on archival data obtained with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Support for this work was provided by an award issued by JPL/Caltech.

Appendix A: Effect of clustering on the uncertainties of number counts

Béthermin et al. (2010) showed how the clustering is linked with the uncertainties of the counts. We used the formalism of Béthermin et al. (2010) to estimate the effect of the clustering on our BLAST counts. There are few sources detected at 250 $\mu$ m and the BLAST coverage is inhomogeneous. It is consequently very hard to estimate the clustering of the resolved population. We thus used the clustering measured at 160 $\mu$ m by Béthermin et al. (2010) and assumed a 250/160 color equal to unity. We then used the same method to compute the uncertainties. We then compare the uncertainties with and without clustering. Neglecting the clustering implies an underestimation of the uncertainties on the counts of 35% in the 203-336 mJy bin at 250 $\mu$ m, and less than 20% in the other bins. We can thus suppose a Poissonian behavior, knowing that the Poisson approximation underestimates the error bars for the 203-336 mJy bin at 250 $\mu$ m. Nevertheless, our model of clustering at 250 $\mu$ m has strong assumptions (single 250/160 color, same clustering at 250 $\mu$ m as measured at 160 $\mu$ m), and it would be more conservative to update it with Herschel clustering measurements.

Appendix B: Effect of clustering on stacking

B.1 A formalism to link clustering and stacking

The clustering can bias the results of a stacking. We present a formalism based on Bavouzet (2008) work.

The expected results for mean stacking of an N non-clustered populations is

$\begin{displaymath}M(\theta) = \overline{S_{\rm s}} \times {{\it PSF}}(\theta) + \int_{0}^{\infty} S \frac{{\rm d}N}{{\rm d}S} {\rm d}S, \end{displaymath}$

(B.1)

where M is the map resulting from stacking, $\theta$ the distance to the center of the cutout image, $\overline{S_{\rm s}}$ the mean flux of the stacked population. The integral is an approximation because the central source is treated in the first term. This approximation is totally justified in a strongly confused field where the number of sources is enormous. PSF is the instrumental response and is supposed to be invariant per rotation ( $\theta=0$ corresponds to the center of this PSF). $\frac{{\rm d}N}{{\rm d}S}$ is the number of the source per flux unit and per pixel. We assume an absolute calibration. The integral in the Eq. (B.1) is equal to the CIB brightness

$\begin{displaymath}I_{\rm CIB} = \int_{0}^{\infty} S \frac{{\rm d}N}{{\rm d}S} {\rm d}S. \end{displaymath}$

(B.2)

This term is constant for all pixels of the image and corresponds to a homogeneous background.

The stacked sources can actually be autocorrelated. The probability density to find a stacked source in a given pixel and another in a second pixel separated by an angle $\theta$ ( $p(\theta)$ ) is linked with the angular autocorrelation function ( $\omega(\theta)$ ) by

$\begin{displaymath}p(\theta) = \rho_{\rm s}^2 (1+\omega(\theta)), \end{displaymath}$

(B.3)

where $\rho_{\rm s}$ is the number density of the stacked source.

If we assume that there is no correlation with other populations, the results of the stacking of N autocorrelated sources is

$\begin{displaymath}M(\theta) = \overline{S_{\rm s}} \times {{\it PSF}}(\theta) +... ... (1+\omega(\theta)) \ast {{\it PSF}}(\theta) + I_{\rm CIB,ns}, \end{displaymath}$

(B.4)

where $I_{{\rm CIB,s}}$ and $I_{\rm CIB,ns}$ is the CIB contribution of stacked and non-stacked sources. If we subtract the constant background of the image, we find

$\begin{displaymath}M(\theta) = \overline{S_{\rm s}} \times {{\it PSF}}(\theta) + I_{\rm CIB,s} \times \omega(\theta) \ast {{\it PSF}}(\theta). \end{displaymath}$

(B.5)

The second term of this equation corresponds to an excess of flux due to clustering. This signal is stronger in the center of the stacked image. The central source appears thus brighter than expected, because of the contribution due to clustering.

The flux of the central stacked source computed by PSF-fitting photometry is

$\begin{displaymath}S_{\rm mes} = \frac{\int \int M \times {\it PSF} {\rm d}\Omeg... ...PSF}^2 {\rm d}\Omega} = \overline{S_{\rm s}} + S{_{\rm clus}}, \end{displaymath}$

(B.6)

where $S_{\rm clus}$ , the overestimation of flux due to clustering is given by

$\begin{displaymath}S_{\rm clus} = I_{{\rm CIB,s}} \times \frac{\int \int \big ( ... ...big ) {\rm d}\Omega}{\int \int {\it PSF}^2 {\rm d}\Omega}\cdot \end{displaymath}$

(B.7)

Basically, the stronger the clustering, the larger the bias. In addition, the wider the PSF, the larger the overestimation. The stacked signal can be dominated by the clustering, if the angular resolution of the instrument is low compared to the surface density of the sources (like Planck, cf. Fernandez-Conde et al. (2010)) or if strongly clustered populations are stacked.

B.2 Estimation of the bias due to clustering

The estimation of $S_{\rm clus}$ with Eq. (B.7) requires particular hypotheses. The stacked population is $S_{24}> 70~\mu$ Jy sources detected by Spitzer. Their contribution to the CIB is 5.8 nW m^-2 sr^-1, 3.4 nW m^-2 sr^-1 and 1.4 nW m^-2 sr^-1 at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively (estimated by direct stacking of all the sources). Following the clustering of 24 $\mu$ m sources estimated by Béthermin et al. (2010), we suppose the following autocorrelation function:

$\begin{displaymath}\omega(\theta) = 2.3\times10^{-4} \times \big(\frac{\theta}{\rm deg}\big)^{-0.8}\cdot \end{displaymath}$

(B.8)

The excess of flux due to clustering ( $S_{\rm clus}$ ) is then 0.44 mJy, 0.35 mJy and 0.16 mJy at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m, respectively. This is significantly lower than the bootstrap uncertainties on these fluxes. We can thus neglect the clustering.

B.3 Measurement of the angular correlation function by stacking

This new formalism provides a simple tool to measure the angular autocorrelation function (ACF) from a source catalog. This method uses a map called ``density map''. One pixel of this map contains the number of sources centered on it. It is equivalent of a map of unit flux sources with the ${\it PSF} = \delta$ (Dirac distribution). The result of the stacking is thus

$\begin{displaymath}M(\theta) = \rho_{\rm s} \times \delta(\theta) + \rho_{\rm s} (1+\omega(\theta)). \end{displaymath}$

(B.9)

The ACF can then be easily computed with

$\begin{displaymath}\forall \theta \neq 0, \omega(\theta) = \frac{M(\theta)}{\rho_{\rm s}}-1. \end{displaymath}$

(B.10)

References

Ashby, M. L. N., Stern, D., Brodwin, M., et al. 2009, ApJ, 701, 428 [NASA ADS] [CrossRef] [Google Scholar]
Barmby, P., Huang, J., Ashby, M. L. N., et al. 2008, ApJS, 177, 431 [NASA ADS] [CrossRef] [Google Scholar]
Bavouzet, N. 2008, Ph.D. Thesis, Université Paris-Sud 11 http://tel.archives-ouvertes.fr/tel-00363975/ [Google Scholar]
Beelen, A., Cox, P., Benford, D. J., et al. 2006, ApJ, 642, 694 [NASA ADS] [CrossRef] [Google Scholar]
Benford, D. J., Cox, P., Omont, A., Phillips, T. G., & McMahon, R. G. 1999, ApJ, 518, L65 [NASA ADS] [CrossRef] [Google Scholar]
Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Béthermin, M., Dole, H., Beelen, A., et al. 2010, A&A, 512, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Borys, C., Chapman, S., Halpern, M., et al. 2003, MNRAS, 344, 385 [NASA ADS] [CrossRef] [Google Scholar]
Chapman, S. C., Blain, A. W., Ivison, R. J., et al. 2003a, Nature, 422, 695 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Chapman, S. C., Helou, G., Lewis, G. F., et al. 2003b, ApJ, 588, 186 [NASA ADS] [CrossRef] [Google Scholar]
Chapman, S. C., Blain, A. W., Smail, I., et al. 2005, ApJ, 622, 772 [NASA ADS] [CrossRef] [Google Scholar]
Chary, R., & Elbaz, D. 2001, ApJ, 556, 562 [NASA ADS] [CrossRef] [Google Scholar]
Chary, R., Casertano, S., Dickinson, M. E., et al. 2004, ApJS, 154, 80 [NASA ADS] [CrossRef] [Google Scholar]
Coppin, K., Halpern, M., Scott, D., et al. 2008, MNRAS, 384, 1597 [NASA ADS] [CrossRef] [Google Scholar]
Daddi, E., Alexander, D. M., Dickinson, M., et al. 2007, ApJ, 670, 173 [NASA ADS] [CrossRef] [Google Scholar]
Devlin, M. J., Ade, P. A. R., Aretxaga, I., et al. 2009, Nature, 458, 737 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Dole, H., Gispert, R., Lagache, G., et al. 2001, A&A, 372, 364 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Dole, H., Lagache, G., & Puget, J.-L. 2003, ApJ, 585, 617 [NASA ADS] [CrossRef] [Google Scholar]
Dole, H., Le Floc'h, E., Pérez-González, P. G., et al. 2004, ApJS, 154, 87 [NASA ADS] [CrossRef] [Google Scholar]
Dole, H., Lagache, G., Puget, J.-L., et al. 2006, A&A, 451, 417 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
Dunne, L., Eales, S., Edmunds, M., et al. 2000, MNRAS, 315, 115 [NASA ADS] [CrossRef] [Google Scholar]
Dye, S., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 703, 285 [NASA ADS] [CrossRef] [Google Scholar]
Fazio, G. G., Ashby, M. L. N., Barmby, P., et al. 2004a, ApJS, 154, 39 [NASA ADS] [CrossRef] [Google Scholar]
Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004b, ApJS, 154, 10 [NASA ADS] [CrossRef] [Google Scholar]
Fernandez-Conde, N., Lagache, G., Puget, J.-L., et al. 2008, A&A, 481, 885 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Fernandez-Conde, N., Lagache, G., Puget, J.-L., et al. 2010, A&A, 515, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Fixsen, D. J., Dwek, E., Mather, J. C., Bennett, C. L., & Shafer, R. A. 1998, ApJ, 508, 123 [NASA ADS] [CrossRef] [Google Scholar]
Franceschini, A., Rodighiero, G., Cassata, P., et al. 2006, A&A, 453, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Franceschini, A., Rodighiero, G., Vaccari, M., Marchetti, L., & Mainetti, G. 2009, A&A, accepted [arXiv:0906.4264] [Google Scholar]
Frayer, D. T., Fadda, D., Yan, L., et al. 2006a, AJ, 131, 250 [NASA ADS] [CrossRef] [Google Scholar]
Frayer, D. T., Huynh, M. T., Chary, R., et al. 2006b, ApJ, 647, L9 [NASA ADS] [CrossRef] [Google Scholar]
Frayer, D. T., Sanders, D. B., Surace, J. A., et al. 2009, AJ, 138, 1261 [NASA ADS] [CrossRef] [Google Scholar]
Gispert, R., Lagache, G., & Puget, J. L. 2000, A&A, 360, 1 [NASA ADS] [Google Scholar]
Gruppioni, C., Pozzi, F., Lari, C., Oliver, S., & Rodighiero, G. 2005, ApJ, 618, L9 [NASA ADS] [CrossRef] [Google Scholar]
Hauser, M. G., & Dwek, E. 2001, ARA&A, 39, 249 [NASA ADS] [CrossRef] [Google Scholar]
Hauser, M. G., Arendt, R. G., Kelsall, T., et al. 1998, ApJ, 508, 25 [NASA ADS] [CrossRef] [Google Scholar]
Ivison, R. J., Greve, T. R., Smail, I., et al. 2002, MNRAS, 337, 1 [NASA ADS] [CrossRef] [Google Scholar]
Ivison, R. J., Smail, I., Dunlop, J. S., et al. 2005, MNRAS, 364, 1025 [NASA ADS] [CrossRef] [Google Scholar]
Kashlinsky, A. 2005, Phys. Rep., 409, 361 [NASA ADS] [CrossRef] [Google Scholar]
Khan, S. A., Chanial, P. F., Willner, S. P., et al. 2009, ApJ, 706, 319 [NASA ADS] [CrossRef] [Google Scholar]
Khan, S. A., Shafer, R. A., Serjeant, S., et al. 2007, ApJ, 665, 973 [NASA ADS] [CrossRef] [Google Scholar]
Knudsen, K. K., Barnard, V. E., van der Werf, P. P., et al. 2006, MNRAS, 368, 487 [NASA ADS] [CrossRef] [Google Scholar]
Kovács, A., Chapman, S. C., Dowell, C. D., et al. 2006, ApJ, 650, 592 [NASA ADS] [CrossRef] [Google Scholar]
Lacy, M., Storrie-Lombardi, L. J., Sajina, A., et al. 2004, ApJS, 154, 166 [NASA ADS] [CrossRef] [Google Scholar]
Lagache, G., Abergel, A., Boulanger, F., Désert, F. X., & Puget, J.-L. 1999, A&A, 344, 322 [NASA ADS] [Google Scholar]
Lagache, G., Dole, H., Puget, J.-L., et al. 2004, ApJS, 154, 112 [NASA ADS] [CrossRef] [Google Scholar]
Lagache, G., Puget, J.-L., & Dole, H. 2005, ARA&A, 43, 727 [NASA ADS] [CrossRef] [Google Scholar]
Le Borgne, D., Elbaz, D., Ocvirk, P., et al. 2009, A&A, 504, 727 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
LeFloc'h, E., Aussel, H., Ilbert, O., et al. 2009, ApJ, 703, 222 [NASA ADS] [CrossRef] [Google Scholar]
Lewis, G. F., Chapman, S. C., & Helou, G. 2005, ApJ, 621, 32 [NASA ADS] [CrossRef] [Google Scholar]
Magdis, G. E., Rigopoulou, D., Huang, J., et al. 2008, MNRAS, 386, 11 [NASA ADS] [CrossRef] [Google Scholar]
Marleau, F. R., Fadda, D., Storrie-Lombardi, L. J., et al. 2004, ApJS, 154, 66 [NASA ADS] [CrossRef] [Google Scholar]
Marsden, G., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 707, 1729 [NASA ADS] [CrossRef] [Google Scholar]
Michałowski, M. J., Hjorth, J., & Watson, D. 2010, A&A, 514, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Papovich, C., Dole, H., Egami, E., et al. 2004, ApJS, 154, 70 [NASA ADS] [CrossRef] [Google Scholar]
Pascale, E., Ade, P. A. R., Bock, J. J., et al. 2008, ApJ, 681, 400 [NASA ADS] [CrossRef] [Google Scholar]
Patanchon, G., Ade, P. A. R., Bock, J. J., et al. 2008, ApJ, 681, 708 [NASA ADS] [CrossRef] [Google Scholar]
Patanchon, G., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 707, 1750 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pearson, C., & Khan, S. A. 2009, MNRAS, 399, L11 [NASA ADS] [CrossRef] [Google Scholar]
Pope, A., Borys, C., Scott, D., et al. 2005, MNRAS, 358, 149 [NASA ADS] [CrossRef] [Google Scholar]
Pope, A., Scott, D., Dickinson, M., et al. 2006, MNRAS, 370, 1185 [NASA ADS] [CrossRef] [Google Scholar]
Puget, J.-L., Abergel, A., Bernard, J.-P., et al. 1996, A&A, 308, L5 [NASA ADS] [Google Scholar]
Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25 [NASA ADS] [CrossRef] [Google Scholar]
Rodighiero, G., Lari, C., Pozzi, F., et al. 2006, MNRAS, 371, 1891 [NASA ADS] [CrossRef] [Google Scholar]
Rowan-Robinson, M. 2009, MNRAS, 394, 117 [NASA ADS] [CrossRef] [Google Scholar]
Sajina, A., Borys, C., Chapman, S., et al. 2003, MNRAS, 343, 1365 [NASA ADS] [CrossRef] [Google Scholar]
Sajina, A., Scott, D., Dennefeld, M., et al. 2006, MNRAS, 369, 939 [NASA ADS] [CrossRef] [Google Scholar]
Shupe, D. L., Rowan-Robinson, M., Lonsdale, C. J., et al. 2008, AJ, 135, 1050 [NASA ADS] [CrossRef] [Google Scholar]
Smail, I., Ivison, R. J., Blain, A. W., et al. 2002, MNRAS, 331, 495 [NASA ADS] [CrossRef] [Google Scholar]
Spinoglio, L., Malkan, M. A., Rush, B., Carrasco, L., & Recillas-Cruz, E. 1995, ApJ, 453, 616 [NASA ADS] [CrossRef] [Google Scholar]
Starck, J. L., Aussel, H., Elbaz, D., Fadda, D., & Cesarsky, C. 1999, A&AS, 138, 365 [Google Scholar]
Stetson, P. B. 1987, PASP, 99, 191 [NASA ADS] [CrossRef] [Google Scholar]
Sullivan, I., Cooray, A., Chary, R., et al. 2007, ApJ, 657, 37 [NASA ADS] [CrossRef] [Google Scholar]
Truch, M. D. P., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 707, 1723 [NASA ADS] [CrossRef] [Google Scholar]
Valiante, E., Lutz, D., Sturm, E., Genzel, R., & Chapin, E. 2009, ApJ, 701, 1814 [NASA ADS] [CrossRef] [Google Scholar]

Footnotes

... website: http://www.blastexperiment.info
... available: On the IAS website http://www.ias.u-psud.fr/irgalaxies/
... beam: The beam solid angles are taken as 0.39 arcmin², 0.50 arcmin² and 0.92 arcmin² at 250 $\mu$ m, 350 $\mu$ m and 500 $\mu$ m respectively.
... catalog: The bins are defined as 64 to 97 at 250 $\mu$ m, 54 to 83 at 350 $\mu$ m and 49 to 76 at 500 $\mu$ m.

All Tables