Online material

Appendix A: Estimation and correction on the bias caused by the galaxy clustering on the stacking results

As explained in Sect. 3.2, the standard stacking technique can be strongly affected by the bias caused by the clustering of the galaxies. We use two independent methods to estimate and correct it.

Appendix A.1: Estimation of the bias using a simulation based on the real catalog

We performed an estimate of the bias induced by the clustering using a realistic simulation of the COSMOS field based on the positions and stellar masses of the real sources. The flux of each source in this simulation is estimated using the ratio between the mean far-IR/(sub-)mm fluxes and the stellar mass found by a first stacking analysis. The galaxies classified as passive are not taken into account in this simulation. This technique assumes implicitly a flat sSFR-M_⋆ relation, since we use a constant SFR/M_⋆ ratio versus stellar mass at fixed redshift. However, we checked that using a more standard sSFR relation (e.g., Rodighiero et al. 2011) has a negligible impact on the results. We applied no scatter around this relation in our simulation for simplicity. As mean stacking is a linear operation, the presence or not of a scatter has no impact on the results (Béthermin et al. 2012b).

A simulated map is thus produced using all the star-forming galaxies of the Ilbert et al. (2013) catalog. In order to avoid edge effects (absence of sources and thus a lower background caused by the faint unresolved sources in the region covered by the optical/near-IR data), we fill the uncovered regions drawing with replacement sources from the UltraVISTA field and putting them at a random position. The number of drawn sources is chosen to have exactly the same number density inside and outside the UltraVISTA field.

Finally, we measured the mean fluxes of the M_⋆> 3 × 10¹⁰M_⊙ sources by stacking in the simulated maps, using exactly the same photometric method as for the real data. We finally computed the relative bias between the recovered flux and the input flux (S_out/S_in − 1). The results are shown Fig. A.1 (blue triangles). The uncertainties are computed a bootstrap method. As expected, the bias increases with the size of the beam. We can see a rise of the bias with redshift up to z ~ 2. This trend can be understood considering the rise of the clustering of the galaxy responsible for the cosmic infrared background (Planck Collaboration XXX 2014) and a rather stable number density of emitters especially below z = 1 (Béthermin et al. 2011; Magnelli et al. 2013; Gruppioni et al. 2013). At higher redshift, we found a slow decrease. This trend is probably driven by the decrease in the infrared luminosity density at high redshift (Planck Collaboration XXX 2014; Burgarella et al. 2013) combined with the decrease in the number density of infrared emitters (Gruppioni et al. 2013).

Fig. A.1 Relative bias induced by the clustering as a function of redshift at the various wavelengths we used in our analysis. The FWHM of the beam is provided in brackets. The blue triangles are the estimations from the simulation (Sect. A.1) and the red diamonds are provided by the fit of the clustering component in map space (Sect. A.2). These numbers are only valid for a complete sample of M⋆> 3 × 1010M⊙ galaxies.

Open with DEXTER

Appendix A.2: Estimation of the bias fitting the clustering contribution in the stacked images

The method presented in the previous section only takes into account the contamination of the stacks by known sources. However, faint galaxy populations could have a non negligible contribution, despite their total contribution to the infrared luminosity and their clustering are expected to be small. We thus used a second method to estimate the bias caused by the clustering which takes into account a potential contamination by these low-mass galaxies. This method is based on a simultaneous fit in the stacked images of three components: a point source at the center of the image, a clustering contamination, and a background. This technique was already successfully used by several previous works based on Herschel and Planck data (Béthermin et al. 2012b; Heinis et al. 2013, 2014; Welikala et al. 2014).

In presence of clustering, the outcome of a stacking is not only a PSF with the mean flux of the population and a constant background (corresponding to the surface brightness of all galaxy populations i.e., the cosmic infrared background). There is in addition a signal coming from the greater probability of finding another neighboring infrared galaxy compared to the field because of galaxy clustering. The signal in the stacked image can thus be modeled by Bavouzet (2008) and Béthermin et al. (2010b)(A.1)where m is the stacked image, PSF the point spread function, and w the auto-correlation function. The symbol ∗ represents the convolution. α, β, and γ are free parameters corresponding to the intensity of the mean flux of the population, the clustering signal, and the background, respectively. This method works only if the PSF is well-known, the extension of the sources is negligible compared to the PSF, and the effects of the filtering are small at the scale of the stacked image. Consequently, we applied this method only to the SPIRE data for which these hypotheses are the most solid. The uncertainties on the clustering bias (β/α for the photometry we chose to use for SPIRE data) are estimated fitting the model described previously on a set of stacked images produced from 1000 bootstrap samples. The results are shown in Fig. A.1 (red diamonds).

Appendix A.3: Corrections of the measurements

In Fig. A.1, we can see that the two methods provide globally consistent estimates. This confirms that the low-mass galaxies not included in the UltraVISTA catalog have a minor impact. We found few outliers for which the two methods disagree. In particular, in the 1.5 <z< 1.75 bin, the estimation from the simulation is higher than the trend of the redshift evolution at all wavelengths, and the results from the profile fitting are lower. This could be caused, as instance, by a structures in the field or a systematic effect in the photometric redshift. Because of these few catastrophic outliers, we chose to use a correction computed

from a fit of the redshift evolution of the bias instead of an individual estimate in each redshift slice.

The evolution of the bias with redshift is fitted independently at each wavelength. We chose to use a simple, second-order, polynomial model (az² + bz + c). We used only the results from the simulation to have a consistent treatment of the various wavelengths. The scatter of the residuals is larger than the residuals, probably because bootstrap does not take into account the variance coming from the large-scale structures. We thus used the scatter of the residuals to obtain a conservative estimate of the uncertainties on the bias. In Fig. A.1, the best fit is represented by a solid line and the 1σ confidence region by a dashed line.

In a few case, the bias at z> 3 can converge to unphysical negative values. We then apply no corrections, but combine the typical uncertainty on the bias to the error bars. A special treatment is also applied to the samples of strong starbursts. Their flux is typically 10 times brighter in infrared by construction (their sSFR is 10 times larger than the main sequence). In contrast, the clustering signal is not expected to be significantly stronger, because the clustering of massive starbursts and main-sequence galaxies is relatively similar (Béthermin et al. 2014). We thus divide the bias found for the full population of galaxy by a factor of 10 to estimate the one of the starbursts for simplicity.

Appendix A.4: Testing another method

We also tried to apply the simstack algorithm (Viero et al. 2013) to our data. This algorithm is adapted from Kurczynski & Gawiser (2010) and uses the position of the known sources to deblend their contamination. Contrary to Kurczynski & Gawiser (2010), simstack can consider a large set of distinct galaxy populations. The mean flux of the each population is used to estimate how sources contaminate their neighbors. All populations are treated simultaneously. This is the equivalent of PSF-fitting codes but applied to a full population instead of each source individually. Unfortunately, this method is not totally unbiased in our case. We found biases up to 15% running simstack on the simulation presented in Sect. A.1, probably because the catalog of mass-selected sources is not available around bright sources. At the edge of the optical/near-IR-covered region, the flux coming from the sources outside the covered area is not corrected, when the flux from all neighbors is taken into account at the middle of zone where the mass catalog is extracted. Indeed, the algorithm works correctly if we put on the simulation only sources present in the input catalog.

Appendix B: Fit residuals

Figures B.1 and B.2 shows the residuals of the fits of our mean SEDs derived by stacking. We did not find any systematic trend, except a 2σ underestimation of the millimeter data in main-sequence galaxies at z> 3.

	Fig. B.1 Residuals of our fit of mean SEDs of main-sequence galaxies by the Draine & Li (2007) model.
Open with DEXTER

	Fig. B.2 Residuals of our fit of mean SEDs of strong starbursts by the Draine & Li (2007) model.
Open with DEXTER

© ESO, 2015