HerMES: COSMIC INFRARED BACKGROUND ANISOTROPIES AND THE CLUSTERING OF DUSTY STAR-FORMING GALAXIES^*

HerMES fields are organized, according to area and depth, into levels 1 through 7, with Level 1 maps being the smallest and deepest (∼310 arcmin²), and level 7 maps the widest and shallowest (∼270 deg²; Oliver et al. 2012).

This study focuses on a subset of the level 5 and 6 fields, totaling ∼70 deg², chosen for their large area and uniformity, and because they have a manageable level of Galactic cirrus contamination. The fields used for this study, and a summary of their properties, are given in Table 1. Combined, they represent an increase of more than four times the area of the initial HerMES study (Amblard et al. 2011). The largest of the HerMES fields, the HerMES Large-Mode Survey (HeLMS)—which was designed specifically to measure the power spectrum on large angular scales—is still in preparation, and will be the subject of a future study. Maps will be made available to the public through HeDaM³⁵ (Roehlly et al. 2011) as a part of data release 2.

The data obtained from the Herschel Science Archive were processed with a combination of standard ESA software and a customized software package SMAP. The maps themselves were then made using an updated version of SMAP/SHIM (Levenson et al. 2010), an iterative map-maker designed to optimally separate large-scale noise from signal. SMAP differs from HIPE (Ott 2010) in three fundamental ways which are relevant for power spectrum studies. First, the standard scan-by-scan temperature drift correction module within HIPE is overridden in favor of a custom correction algorithm which stitches together all of the time-ordered data (or timestreams), allowing us to fit to and remove a much longer noise mode. Further, the standard processing is modified such that a "sigma-kappa" deglitcher is used instead of a wavelet deglitcher, to improve performance in large blank fields. Lastly, imperfections from thermistor jumps, the "cooler burp" effect, and residual glitches are removed manually before map construction. Detailed descriptions of the updates to the SMAP pipeline of Levenson et al. (2010) are presented in Appendix A.

Following Amblard et al. (2011), we make maps with 10 iterations, fewer than SMAP's default of 20, in order to minimize the time needed to measure the transfer functions (Section 3.1.2) and uncertainties (Section 3.2) with Monte Carlo simulations.

Additionally, timestream data are divided into two halves and unique "jack-knife" map-pairs are made. These map-pairs are those ultimately used for estimating power spectra (Section 3.1). SMAP maps are natively made with pixel sizes of 6'', 833, and 12'', which is motivated by the beam size (sampling them by ∼1/3 FWHM). But because the cross-frequency power spectra calculations need maps of equivalent pixel sizes, three additional sets of map-pairs are made, with the extra sets having custom pixel sizes so that the cross-frequency power spectra can be performed at the pixel resolution native to the maps with the larger instrumental beams. In other words, the maps used to calculate the 250 × 350 μm spectra have identical 833 pixels, and those used for 250 × 500 and 350 × 500 μm have 12'' pixels.

At 100 μm, we use the Improved Reprocessing of the IRAS Survey³⁶ (IRIS; Miville-Deschênes & Lagache 2005), a data set which corrects the original plates for calibration, zero level, and striping problems. The resulting FWHM resolution and noise level are 43 ± 02 and 0.06 ± 0.02 MJy sr⁻¹, respectively, and the gain uncertainty is 13.5%. Data are available for up to three independent observations, or HCONs, although two of our fields were only observed twice. For fields in which tiles intersect, maps can be stitched together with custom software provided on their site.³⁷

The intensity in a given SPIRE map, I_map, can be approximated as

$$\begin{equation} I_{\rm map}=(T\otimes[ I_{\rm sky}\otimes B + N])W, \end{equation} \tag{ 3 }$$

where I_sky is the sky signal we wish to recover, T is the transfer function of the map-maker, B is the instrumental beam, N is the noise, and W is the window function, which includes the masking of map edges and of bright sources. We use ⊗ to represent a convolution in real-space. The instrumental noise, N, is made up of white noise, which dominates on angular scales k_θ ≳ 0.2 arcmin⁻¹, and 1/f noise. We note that as T has structure in 2D, and as the true beam may vary slightly across the map (particularly for bigger maps). These corrections are small enough that Equation (3) remains a reasonable approximation.

In the auto-correlation of a map (i.e., the auto-power spectrum of a single map), all of the power present—which includes both signal and noise—is correlated, while in the cross-correlation of jack-knife map-pairs, in principle the only signal correlated between them is the sky signal. Since we are interested in recovering the sky signal, we estimate the power spectrum from the cross-correlation of jack-knife map-pairs (discussed in Section 2.1). In practice, some correlated noise could exist between maps, particularly on large scales, as 1/f. To minimize this, we use map-pairs constructed by dividing the timestreams in half by time, which ensures that the maps are made from data taken at time intervals corresponding to very large scales. Note, this would not be true if, say, the data were split into those from even and odd detectors, since the same large-scale noise would be present in both maps. The remaining 1/f results from serendipitous alignment of large-scale noise (see Appendix B of Hajian et al. 2012), i.e., independent noise clumps in either map that happen to line up. Also note that excessive high-pass filtering of the time-ordered data before constructing maps would attenuate extragalactic large-scale signal along with unwanted signal, and is thus not employed.

The one-dimensional power spectrum is the azimuthal average of the (nearly isotropic) two-dimensional power spectrum of map-pairs in k-space. In order to recover the true power spectrum of the sky, the cross-spectrum of the map-pairs must be corrected for the transfer function, masking, and instrumental beam. We now describe these corrections in detail.

3.1.1. Masking and the Mode-coupling Matrix

Before calculating the power spectra, the maps are multiplied with windows whose values equal unity in the clean parts of the maps and taper to zero with a Gaussian profile (90'' FWHM) at the edges. Note that jack-knife maps typically do not cover identical patches of sky because the orientation of the telescope between observations changes, so that the windows of the two jack-knife map-pairs are different.

Next, because we are interested in the behavior of the power spectrum with the masking level of resolved point sources, we create a series of masks for each field and band to mask sources whose flux densities are greater than 50, 100, 200, and 300 mJy, in addition to extended sources. Extended sources are exceptional objects, e.g., local IRAS galaxies, which exceed 400 mJy but can be as bright as 1500 mJy at 250 μm. Mask positions are determined from catalogs of resolved sources which we construct. The source identification code first high-pass filters the maps in Fourier space to remove Galactic cirrus and other large-scale power, then convolves the maps by the instrumental beam, and finally measures all peaks with signal-to-noise ratio greater than 3σ. The total number of masked sources in each field and band is tabulated in Table 2. Note that this method could potentially suffer from Eddington bias (e.g., Chapin et al. 2009), a phenomenon where instrumental white noise systematically boosts faint sources above the detection limit. To ensure that this is not a significant problem, we compare the total number of sources above the cut to cumulative number counts from Glenn et al. (2010) and find that they are consistent.

Table 2. Number of Masked Sources

ObsID	300 mJy	200 mJy	100 mJy	50 mJy
	18	38	180	1567
bootes	4	9	46	723
	1	2	13	248
	9	23	165	1569
cdfs-swire	3	6	30	725
	1	1	6	141
	9	23	104	1073
elais-s1	1	3	22	437
	1	1	3	146
	16	42	221	2043
lockman-swire	4	9	57	930
	1	3	10	220
	23	66	368	2973
xmm-lss	5	13	85	1170
	1	2	16	437

Notes. Number of sources masked at 250, 350, and 500 μm is given from top to bottom in each panel. Not accounted for in the table are the extended sources masked in each field, of which there are 3, 2, 2, 4, and 4, respectively. A 50 mJy cut amounts to masking approximately 1.4%, 1.2%, and 0.5% of the pixels at 250, 350, and 500 μm, respectively.

Download table as:  ASCIITypeset image

Each source is masked by circles of 1.1 × FWHM in diameter, or 191, 277, and 403 at 250, 350, and 500 μm, respectively—chosen to cover the full first lobe of the beam, though we check that the exact size of the mask has a negligible effect on the spectra. For the data that concern this paper, unique masks are made for each field and each band (i.e., one at 250, one at 350, and one at 500 μm), so that only sources above the given cut in that band are masked. This means that when calculating the cross-frequency power spectra, not all of the sources masked at, e.g., 250 μm will also be masked at 350 μm, and vice versa. However, we additionally calculate an alternative set of spectra where we mask in all bands the sources identified at 250 μm (i.e., the same mask at 250, 350, and 500 μm), and the power spectrum pipeline is rerun. Plots and tables for this alternate masking scheme are presented in Appendix D, where we also show that the spectra at longer wavelengths are less sensitive to the level of source masking than at shorter wavelengths. Finally, we note that this method differs from that of Amblard et al. (2011), who instead masked all pixels above 50 mJy, as well as all neighboring pixels, the motivation being that a catalog-based masking scheme is better able to distinguish sources from spurious noise. The total number of masked sources in each field and band is tabulated in Table 2.

Masking in map-space can result in mode-coupling in Fourier space, which can bias the power spectrum. The coupling kernel, or mode-coupling matrix (Hivon et al. 2002), in the flat sky approximation is

$$\begin{equation} M_{kk^{\prime } } = \displaystyle \sum \limits _{\theta _k}\displaystyle \sum \limits _{\theta _{k^{\prime }}}\left| w_{kk^{\prime } }\right|^2/N(\theta _k), \end{equation} \tag{ 4 }$$

where $\left| w_{kk^{\prime } }\right|^2$ is the auto-power spectrum of the mask, and N(θ_k) is the number of modes in annulus of radius k. Since in our case the masks of the map-pairs are not identical, Equation (4) is generalized for different masks by replacing $\left| w_{kk^{\prime } }\right|^2$ with the cross-power spectrum of the masks, $\langle w^A_{kk^{\prime } } w^{^\ast B}_{kk^{\prime } } \rangle$ (see Tristram et al. 2005).

The mode-coupling matrix must be inverted in the final step in order to recover the de-coupled power spectrum. A unique mode-coupling matrix is calculated for each auto- and cross-frequency power spectrum, and for each flux cut, per field. Additionally, each $M_{kk^{\prime }}$ is tested on 1000 simulated maps with steep input spectra and found to be unbiased.

3.1.2. Transfer Function

The large-scale correlated noise of SPIRE is extremely low (e.g., Pascale et al. 2011). As a result, only minimal high-pass filtering is required to make well-behaved maps, and the resulting power spectrum can be measured out to relatively large scales (k_θ ≳ 0.01  arcmin⁻¹). This filtering is quantified by the transfer function of the map-maker, T, which must be accounted for in the final spectra.

We measure T with a Monte Carlo simulation whose steps are: (1) running the SMAP map-making pipeline on simulated observations of input maps with known power spectra resembling those of clustered DSFGs; (2) calculating the power spectrum of the output map with a pipeline identical to that used for the real data, including all masking, Fourier space filtering, and mode-coupling corrections; and (3) computing the average of the ratio of the output power spectrum to the known input spectrum. Transfer functions are calculated from 100 simulations of each field and wavelength, and are shown in Figure 2. We check that the spectra in step (1) are not sensitive to the steepness of the input spectra, and that the transfer functions have converged, with a mean error of ∼1%.

Zoom In Zoom Out Reset image size

As anticipated, all three bands converge on large scales for each field, meaning that the same filtering is performed in each band. And on small scales (k_θ > 0.3 arcmin⁻¹) all fields converge to the same three curves, which are the pixel window functions of the three bands. The angular scales on which the high-pass filtering occurs in each field are related to the average speeds and lengths of the maps scans: the former determines the scale in which 1/f noise is projected onto the timestreams, while the latter dictates the order of the polynomial removed from the timestreams. As described later in Section 4.2, the final spectra are weighted combinations of those in each field which are attenuated by less than 50%, corresponding to k_θ ⩾ 0.021, 0.009, 0.023, 0.015, and 0.009 arcmin⁻¹, for bootes, cdfs-swire, elais-s1, lockman-swire, and xmm-lss, respectively. Maps scanned with longer and faster scans have less attenuation on large scales, which is why the maps observed in fast-scan mode (cdfs-swire and lockman-swire) are also those which best measure the largest scales. Note, the excess power introduced on large scales by earlier versions of SMAP (Levenson et al. 2010; Amblard et al. 2011) is no longer present.

3.1.3. Instrument Beam

Present in each map-pair is correlated signal from the sky, and both correlated and uncorrelated noise. Three terms contribute to the uncertainties in the power spectrum: a non-Gaussian term due to the Poisson distributed compact sources, sample variance in the signal due to limited sky coverage, and the noise. In this order, the variance of the cross-spectra of maps A × B can be written as

$$\begin{eqnarray} \sigma ^2\left(\hat{P}_{b}^{A\times B}\right) & = & \frac{\sigma ^2_{\rm P}}{f_{\rm sky}} + \frac{2}{n_b}\left(\hat{P}_{b}^{A\times B}\right)^2 \nonumber \\ &&+ \frac{\hat{P}_{b}\left(\hat{N}_{b}^{A} + \hat{N}_{b}^{B}\right)+\hat{N}_{b}^{A}\hat{N}_{b}^{B}}{n_b}, \end{eqnarray} \tag{ 5 }$$

where $\hat{P}_{b}$ is the mean cross-spectrum of map-pairs, $\hat{N}_{b}$ is the average noise power spectrum of the map, n_b is the number of Fourier modes measured in bin b, and f_sky is the observed area divided by the solid angle of the full sky. The first term, $\sigma ^2_{\rm P}$, is given by the non-Gaussian part of the four-point function (as described in, e.g., Acquaviva et al. 2008; Hajian et al. 2012), and is particularly sensitive to the flux cut of the masked sources.

Shown in Figure 3 are the noise levels calculated from the power spectrum of the difference map of jack-knife map-pairs, which are consistent with the difference between the auto- and cross-power spectra of the maps (not shown). The noise behavior demonstrates the impressive performance and stability of the SPIRE instrument, with white noise in most cases nearly two orders of magnitude below the power from the sky signal. The noise spectra turn over on large scales due to the filtering performed by the map-maker, which is related to the length and speed of the scans.

Zoom In Zoom Out Reset image size

Uncertainties in the estimates of the power spectra are derived from Monte Carlo simulations of the pipeline on realistically simulated sky-maps. The maps include sources correlated between bands (i.e., the same sources appear in all three maps, but with different flux densities), which are necessary for estimating uncertainties in the cross-frequency power spectra, as well as both 1/f and white noise, and Galactic cirrus. Also included in the Monte Carlo simulations are systematic uncertainties arising from the beam and transfer function corrections, such that for each iteration, the beam and transfer function corrections are perturbed by the appropriate amount.

The ensemble of estimated output power spectra is used to measure, V, the covariance matrix

$$\begin{equation} {\bf V}_{bb^{\prime }}=\langle (P_b - \tilde{P}_b ) (P_{b^{\prime }} - \tilde{P}_{b^{\prime }})\rangle _{{\rm MC}}, \end{equation} \tag{ 6 }$$

where the tilde denotes the mean over every iteration in bin b. The resulting errors are

$$\begin{equation} \sigma _{P_{b}^{\rm map}}=\sqrt{{\bf V}_{bb}}. \end{equation} \tag{ 7 }$$

The non-Gaussian term emerges from the simulations as an offset in V. We check that this level is realistic by comparing it to the level of the four-point function estimated directly from data, with appropriate masking, following Fowler et al. (2010), which we found to be between 5% and 10% of the total error in the Poisson dominated regime, depending on flux cut of masked sources: more aggressive source masking results in a smaller non-Gaussian term.

In addition, there are ∼8% systematic errors due to absolute calibration uncertainty, of which ≲ 1% is due to beam area uncertainty, as described in Appendix B. Though they are accounted for when model fitting, they are not included in the reported error bars.

The most significant foreground for the extragalactic power spectrum is that from Galactic cirrus, which can dominate the signal on scales greater than ∼30'. Gautier et al. (1992) showed that the power spectrum of Galactic cirrus can be well approximated by a power law

$$\begin{equation} P_{k_{\theta }}^{\rm cirrus}=P_0\left(\frac{k}{k_0} \right)^{\alpha _{\rm c}}, \end{equation} \tag{ 8 }$$

whose amplitude, P₀, normalized at k₀ = 0.01 arcmin⁻¹, may vary from field to field, but whose index α_c ≈ −3.0. More recently, studies of various fields using BLAST and SPIRE data have found a much wider range of the index, e.g., α_c ∼ −2.4 (Bracco et al. 2011), −2.6 (Roy et al. 2010), −2.8 (Martin et al. 2010; Miville-Deschênes et al. 2010), −2.9 (Lagache et al. 2007), and even ranging from −2.5 to −3.6, depending on the field (Miville-Deschênes et al. 2007). Thus, rather than assume a power law with α = −3 (e.g., Viero et al. 2009), we proceed by treating each field independently.

In diffuse cirrus regions (i.e., column density, $N_{\rm H\,\scriptsize{I}} < 2 \times 10^{20}\, \rm cm^{-2}$ and brightness temperature, T_b < 12 K; e.g., Lockman & Condon 2005; Gillmon & Shull 2006), H i is a good tracer of dust, and the dust-to-gas ratio can be measured from the slope of the pixel-pixel scatter plot. Cirrus contamination can then be "cleaned" by scaling the H i maps by the dust-to-gas ratio and subtracting them directly from the dust maps in question. This approach has been used quite successfully by, e.g., Planck Collaboration et al. (2011b) or Pénin et al. (2012b), on maps which have high fidelity on large angular scales. Unfortunately, this technique is made complicated for SMAP-made SPIRE maps because the filtering of scales larger than ≳ 20 arcmin attenuates the very structure that the differencing with H i is meant to remove. Though H i maps can be filtered with the SMAP simulator to attenuate large scales, the remaining structures are faint with respect to the noise, and thus difficult to regress with SPIRE maps.

We instead adopt an approach similar to that used by Lagache et al. (2007), Viero et al. (2009), and Amblard et al. (2011), with some additional modifications. Diffuse Galactic cirrus emits as a modified blackbody proportional to ν^βB(ν), where B(ν) is the Planck function and β is the emissivity index. Typically, it has a temperature of ∼18 K and β ∼ 1.8 (e.g., Bracco et al. 2011), resulting in an SED which peaks at ∼170 μm (e.g., Martin et al. 2012). At 100 μm, cirrus emission has roughly the same amplitude as at 250 μm, but unlike in SMAP/SPIRE maps, the favorable large-scale properties of the IRIS maps make it possible to accurately measure the power spectra out to scales of ∼4°. Thus, assuming that the Galactic cirrus power spectrum is well described by a power law (Roy et al. 2010), we use the 100 μm power spectra, calculated from IRIS maps with sizes identical to their SPIRE counterparts, and with sources above 500 mJy masked, to estimate the best fit to the cirrus spectra in each field. Note that although larger size regions would better constrain the large-scale spectra, we intentionally use the exact same regions because these fields were chosen specifically because they were special places in the sky with low Galactic cirrus, and thus the spectra inside and the spectra surrounding the field are unlikely to be the same. Uncertainties in the power spectra are estimated analytically following Fowler et al. (2010). To distinguish between the power originating from cirrus and that from clustered galaxies, we include an estimate of the linear power (i.e., 2-halo) term constrained by the measured galaxy spectra of Pénin et al. (2012b). Also, we adopt the Béthermin et al. (2011) model to fix the Poisson level, which is unconstrained by data because of the 4' IRAS beam, although we note that on these angular scales the exact choice for the Poisson level has a negligible effect on the fit.

Next, assuming that the linear power spectrum from clustered DSFGs is independent of field, and after masking all sources above 300 mJy in the SPIRE bands, we estimate the contribution to the SPIRE spectra from cirrus by fitting the 100 μm and SPIRE auto-power spectra of all five fields simultaneously with: a Poisson term, 1- and 2-halo clustered galaxy terms, and a temperature (with fixed β = 1.8) which sets the band-to-band amplitudes, P₀.

Lastly, uncertainties are estimated with a Monte Carlo simulation where the slope and amplitudes of the best-fit power law at 100 μm are perturbed by an amount dictated by their errors, and the cirrus estimate pipeline described above is rerun with those values fixed. Results are given in Table 3. Note, gain uncertainties in the IRAS maps are not accounted for in the fit, as they would only act to increase the error of the best-fit temperature, but not the uncertainty in the best-fit power law. Also note that many of the indices are steeper than those of most previous analyses of cirrus power spectra, which may be due to the fact that these regions are specifically chosen as windows through the cirrus, and not representative of the mean. And although they are steep, they remain consistent with results at shorter wavelengths from Bazell & Desert (1988), or the extreme end of spectra found by Miville-Deschênes et al. (2007). Nevertheless, we check that fixing the slope of spectra in each field to α = −3 does not significantly alter the correction, and indeed find that the resulting cirrus-corrected data fall within the uncertainties.

We find that the method provides good constraints for the lockman-swire, cdfs-swire, and xmm-lss fields, while in the bootes and elais-s1 fields, because of the aggressive filtering in the SPIRE maps, the measured spectra only probe scales in which the cirrus contribution is negligible. Consequently, when later combining the spectra, the largest scale bins are constrained using a subset of the maps.

Individual spectra in each field, for each band and flux cut, are first rebinned following Amblard et al. (2011) into logarithmic intervals with width equal to Δk_θ/k_θ = 0.25 for k_θ ⩾ 0.033 (Δℓ/ℓ = 720), and linearly for larger scales, with bin widths of Δk_θ = 7.41 × 10⁻³ arcmin⁻¹ (Δℓ = 160). These bin-width values are chosen to ensure that, with the exception of correlations introduced on small scales from Poisson errors, the off-diagonals in the covariance matrix are always less than ∼10%. The rebinned uncertainties are given by

$$\begin{equation} \sigma _{P_b^{\rm sky}} ^2=\frac{ 1}{ \displaystyle \sum \limits _{i,j} \left({\bf V}_{i,j}^{-1}\right)}, \end{equation} \tag{ 9 }$$

where i, j span the entries of bin b and $({\bf V}_{i,j}^{-1})$ is the inverse of the subset of the covariance matrix V calculated from simulations (Section 3.2).

Next, the best estimates of the cirrus power spectra for each field (estimated in Section 3.3 and reported in Table 3) are subtracted in order to recover the power spectra of extragalactic sources, $P_{k_{\theta }}^{\rm exgal}$. Uncertainties in the cirrus estimate are propagated into the final uncertainties assuming that the errors are uncorrelated so that

$$\begin{equation} \sigma _{P_b^{\rm exgal}} ^2=\sigma _{P_b^{\rm sky}} ^2 + \sigma _{P_b^{\rm cirrus}} ^2. \end{equation} \tag{ 10 }$$

Finally, data in each band and flux cut are combined for the five fields, f, following Planck Collaboration et al. (2011b), where

$$\begin{equation} P_b^{\rm combined}= \displaystyle \sum \limits ^5_{f=1} W^f_b \times P_b^{f,\rm exgal}, \end{equation} \tag{ 11 }$$

and $W^f_b$ is the weight of each field and bin,

$$\begin{equation} W^f_b= \frac{\sigma ^{-2}_{P^{f,\rm exgal}_b}}{\displaystyle \sum \limits ^5_{f=1} \sigma _{P^{f,\rm exgal}_b}^{-2}}, \end{equation} \tag{ 12 }$$

which assumes that fields are far enough apart to be uncorrelated. Note that for each field, spectra at angular scales where the transfer function falls below 0.5 are omitted in the combined fit.

We measure signals in excess of Poisson noise in all auto- and cross-frequency power spectra, in each field. This excess signal originates from the clustering of DSFGs, and to varying degrees from Galactic cirrus on large scales. Total power spectra of the five fields (which includes power from Galactic cirrus and Poisson noise) with sources ⩾100 mJy masked are shown in Figure 4. Spectra with different levels of source masking behave similarly, and are thus here omitted for clarity. Also shown are dashed lines representing the best estimate of the Galactic cirrus, approximated as power laws. The spectra from each field agree within errors on angular scales where the power from Galactic cirrus is subdominant (≳ 0.06 arcmin⁻¹).

Zoom In Zoom Out Reset image size

The cirrus-subtracted power spectra in each field, which are combined using Equation (11), are presented in Figure 5 and Table 10. As expected, the spectra in each panel converge for increasing angular scales as the contribution from Poisson noise becomes subdominant. Consistent with expectations from Figure 1, the flux cut has a significant effect on the Poisson level at 250 μm, and a nearly negligible effect at 500 μm.

Zoom In Zoom Out Reset image size

Poisson levels are determined through a simultaneous fit to the combined spectra of the Poisson and clustered galaxy terms with templates adopted from the Viero et al. (2009) halo model, and are shown as a function of flux cut of masked sources in Figure 6 and tabulated in Table 8. We note that this estimate is subject to systematic uncertainties due to the mild degeneracy of the Poisson and 1-halo terms, more so for spectra with fewer masked sources, and that we account for those uncertainties in the estimate. As anticipated, shorter wavelengths are significantly more affected by removal of the brightest sources (see Figure 1). Previous measurements from Amblard et al. (2011) with 50 mJy sources masked are in relatively good agreement, with 250 μm higher by ∼5%, and 350 and 500 μm lower by ∼7% and 17%, respectively. The BLAST measurements, which masked sources with flux densities greater than 500 mJy (Viero et al. 2009; Hajian et al. 2012), appear to be higher than our values by ∼26%, 12%, and 11% at 250, 350, and 500 μm, respectively.

Zoom In Zoom Out Reset image size

Also plotted are estimates of the Poisson level derived from the Glenn et al. (2010) P(D) number counts using Equation (2). We find that the number count predictions overestimate our values by ∼16% at 250 μm, and underestimate our measured values by ∼12% and 16% at 350 and 500 μm, respectively. It should be noted that differences between different cuts in a single band are correlated, and that the shot noise levels fall within the calibration uncertainties.

Finally, we compare to the model predictions of Béthermin et al. (2011), a phenomenological model which to date is the best at reproducing the observed number counts from 15 μm to 1.1 mm. We find that the model is in very good agreement with the data at all flux-cut levels and at all but 350 μm, which underpredicts the data by approximately 1σ.

In Figure 7 we plot our combined auto-frequency power spectra along with a selection of recently published CIBA measurements. We show the two masking extremes of our data: those in which all sources greater than 50 mJy were masked (dark blue open circles), and those where only extended sources were masked (light blue open circles). We do this in order to adequately compare with the wide range of masking in the literature.

Zoom In Zoom Out Reset image size

Shown as black plus signs are the SPIRE auto-frequency power spectra of ∼15 deg² from Amblard et al. (2011) in which pixels greater than 50 mJy were masked, so that they should be compared to our dark blue circles. On larger scales (e.g., k_θ ≲ 0.08 arcmin⁻¹), particularly at shorter wavelengths, their spectra suffer from overcorrection of Galactic cirrus contamination (see discussion in Planck Collaboration et al. 2011b). On smaller scales we find that our spectra differ by factors of ∼0.91 ± 0.01, 1.09 ± 0.01, and 1.09 ± 0.01. Note that the calibration of the science demonstration phase maps used in Amblard et al. (2011) differed by 1.02, 1.05, and 0.94, at 250, 350, and 500 μm, respectively, but that those corrections were not applied here. These calibration differences, combined with the offsets resulting from the new estimate of the beam (Appendix B), may partially account for this difference. Also note that although the error bars on their data are comparable to ours at small angular scales, correlations due to the non-Gaussian term (first term in Equation (5)) were not included when they rebinned into log bins, thus artificially deflating their errors.

Shown as red squares at 350 and 500 μm are results from the Planck Collaboration et al. (2011b), which should be compared to our light blue circles. Note, comparisons of the two spectra must be made with caution, bearing in mind that the flux density of masked sources in Planck is much higher (710 and 540 mJy at 350 and 550 μm). In addition, because the passbands of the two instruments are not the same, Planck data at 857 and 545 GHz (350 and 550 μm) are color corrected by multiplying them by factors of 0.99 and 1.30, respectively.

We find possible inconsistencies between the two sets of measurements. On scales k_θ ≲ 0.04 arcmin⁻¹, the Planck spectra appear to be offset high by factors of 1.26 ± 0.06 and 1.40 ± 0.06 at 350 and 500 μm, respectively, while on smaller angular scales there is an apparent excess of power in the Planck data, though it should be noted that, point by point, the 350 μm values do agree within errors. Potential explanations for this discrepancy include calibration, excess Poisson, or reconstruction systematic errors in the Planck beam. We discuss these scenarios in more detail in Section 6.2.

In Figure 8 we plot our combined auto-frequency power spectra next to a selection of published halo models. As in Figure 7, we show the two masking extremes of our data in order to adequately compare with the wide range of masking in the literature.

Zoom In Zoom Out Reset image size

The Viero et al. (2009) models (brown dotted lines), which were fit to BLAST data with sources greater than 500 mJy masked and appeared to be a good match to the Planck Collaboration et al. (2011b) data, here do a poor job of describing the SPIRE measurements, overestimating the power on scales greater than ∼40''.

The Amblard et al. (2011) models, which assumed a masking level of 50 mJy and were fit to each band individually, are shown as gray dot-dashed lines. Similar to the differences in the data, the large-scale power is underestimated due to the overcorrection for the contribution from Galactic cirrus, while on small scales, the small differences may be due to difference in the calibration.

The halo models from the Planck Collaboration et al. (2011b) are shown as red dashed lines at 350 and 500 μm. They assumed a masking level consistent with the masking level of Planck data, as do all following models. Their published data must be corrected for calibration by dividing them by factors of 1.14 and 1.30 at 350 and 500 μm, respectively (Planck Collaboration et al. 2013), after which they agree very well with our data.

The Xia et al. (2012) halo model, shown as orange triple-dot-dashed lines in Figure 8, was fit to Planck and corrected SPIRE data from Planck Collaboration et al. (2011b, Section 5.3). It adopts a description for the source population from Lapi et al. (2011) (an update of Granato et al. 2004). It appears to be consistent with the overall amplitude of the data, but has a bump of excess emission at around k_θ ∼ 0.1–0.03 that the data do not show. This evolutionary model assumes that the steep part of the source counts at submillimeter wavelengths is dominated by massive, proto-spheroidal galaxies in the process of forming most of their stars. The model also includes small contributions from late-type and starburst galaxies. Notable in this model is that the redshift distribution of the emission peaks at slightly higher redshifts, broadly around z ∼ 1.7–2.2, increasing with increasing wavelength; this is in distinction to other models which are strongly peaked at z ∼ 1. It was fit to data from Herschel/SPIRE (Amblard et al. 2011), Planck/HFI (Planck Collaboration et al. 2011b), SPT (Shirokoff et al. 2011), and ACT (Dunkley et al. 2011), and as it is physically based, it is much more constrained than the phenomenological models used by, e.g., Viero et al. (2009) or Shang et al. (2012). That two populations are represented is evidenced by the clear feature at ∼0.2 arcmin⁻¹, a feature which is not apparent in the data, suggesting an overestimate of the contribution of late-type galaxies to the total spectrum.

Lastly, shown as green dashed lines is the Shang et al. (2012) model, whose main feature is to implement a luminosity–mass (L–M) relation, such that more massive halos host more luminous sources. Though the model was fit primarily to Planck data, it appears to fit our spectra at 500 μm quite well. The fit is less good at 250 and 350 μm; though the shape is in good agreement, the curves are high by ∼30%. Despite this success, the model has some points of concern. In particular, it underpredicts the contribution of lower redshift sources to the CIB (e.g., they are significantly below the lower limits measured from stacking 24 μm selected sources; Jauzac et al. 2011), and we also find an uncharacteristically large β (discussed further in Section 6.3).

Finally, we highlight the feature that appears in the data at k_θ ∼ 0.03–0.04 arcmin⁻¹, particularly in all auto- and cross-frequency spectra, but predominantly at 500 μm. A similar feature was visible in the SPIRE data of Amblard et al. (2011) as well. It is not clear if this is a real feature in the sky, which would be unexpected and is unlikely, or due to noise present in the data. Similar features do not appear in the transfer function, or in the simulations of the pipeline which tested for potential biases.

4.4.1. Clustered Galaxy Power Spectra

Ultimately, we are interested in the power spectra of clustered DSFGs, which we estimate by removing the Poisson noise from the cirrus-subtracted, combined power spectra of Figure 5. Results are shown in Figure 9. The fraction anisotropy is calculated as

$$\begin{equation} \delta I/I = \sqrt{2\pi k_{\theta }^2 P(k_{\theta })}/I_{\rm CIB}, \end{equation} \tag{ 13 }$$

where I_CIB is the overall amplitude of the CIB, measured to be 0.71 ± 0.17, 0.59 ± 0.14, and 0.38 ± 0.10 MJy sr⁻¹ at 250, 350, and 500 μm, respectively (e.g., Lagache et al. 2000; Marsden et al. 2009), and k_θ is converted to sr. We find δI/I = 14% ± 4%, consistent with findings from Viero et al. (2009).

Zoom In Zoom Out Reset image size

We fit a simple power law to the clustering spectra over the range k_θ = 0.01–1.4 arcmin⁻¹, finding a mild change in the slope with changes in the masking level. Specifically, from most (S_cut > 50 mJy) to least aggressively masked (only extended sources), we find slopes of −1.60 ± 0.05 to −1.50 ± 0.07 at 250 μm, −1.52 ± 0.05 to −1.36 ± 0.06 at 350 μm, and −1.52 ± 0.06 to −1.47 ± 0.06 at 500 μm. The χ² of these fits for 15 degrees of freedom (dof) are 7, 10, and 3 (reduced χ² ∼ 0.5, 0.7, and 0.3) at 250, 350, and 500 μm, respectively.

The next notable feature is the reduction of 1-halo power with flux cut of masked sources, particularly at 250 μm, which is shown as dashed lines in Figure 9, whereas the 2-halo power, shown as dotted lines, remains relatively unchanged. We demonstrate in Section 6.1 how this result can be interpreted as more luminous sources residing in more massive halos, motivating the use of a model later in the paper in which an L–M relationship is invoked (e.g., Sheth 2005; Skibba et al. 2006; Shang et al. 2012). We show that attempting to account for the reduction of power entirely with the Poisson term leads to significant tension in the fit. Though the reduction of power is much less significant at 500 μm, we remind the reader that there are far fewer sources at each given flux-cut level at 500 μm than at 250 or 350 μm (see Table 2). The capability to mask fainter sources reliably at 500 μm would require either maps with higher angular resolution (i.e., less confusion noise) or a way to probe deeper into the confusion using ancillary data (e.g., XID; Roseboom et al. 2010).

The cross-correlation power spectrum is defined as

$$\begin{equation} C_{\rm A\times B} = \frac{P_{k_{\theta }}^{\rm A\times B}}{\sqrt{P_{k_{\theta }}^{\rm A}\,\cdot \,P_{k_{\theta }}^{\rm B}}}, \end{equation} \tag{ 14 }$$

i.e., the ratio of the cross-frequency power spectra to the geometric mean of the two auto-frequency power spectra. Identical maps would thus have a cross-correlation of unit amplitude, as would maps containing sources located at identical redshifts and with identical colors. Departures from unity would be an indication that sources are not all at the same redshift, or that their colors (or average temperatures) are variable, and the strength and shape of the cross-correlation signal would depend on the level of correlation between maps. Consequently, the cross-correlation provides strong constraints for source population models.

We show the cross-correlations as functions of the flux cut of masked sources in Figure 10. The measurement becomes very uncertain at angular scales k_θ ⩽ 0.1 arcmin⁻¹. At larger k_θ, with the exception of the 50 mJy cut, we find similar levels of correlation for all levels of masking, which can be approximated as horizontal lines at 0.95 ± 0.04, 0.86 ± 0.04, and 0.95 ± 0.03 for A × B = 250 × 350, 250 × 500, and 350 × 500, respectively. Cross-correlations of maps with sources masked at 50 mJy appear to be less correlated, with hints of a reduction in the correlation with increasing k_θ, which, as first predicted by Knox et al. (2001), would be an indication that longer wavelengths are more sensitive to higher z.

Zoom In Zoom Out Reset image size

These results compare favorably with the Planck Collaboration et al. (2011b) measurements of the cross-correlation, who found 0.89 and 0.91 for two different fields at 350 × 550 μm, and are also consistent within errors with the cross-correlations measured by Hajian et al. (2012).

The power spectrum of CIBA in the halo model formalism is written as the sum of three terms: the linear (or 2-halo) term, which accounts for pairs of galaxies in separate halos and dominates the spectrum on large scales; the nonlinear (or 1-halo) term, which describes pairs of galaxies residing in the same halo and is the dominant term on small scales; and the Poisson (or shot) noise term:

$$\begin{equation} P_{\nu \nu ^{\prime }} (k, z) = P_{\nu \nu ^{\prime }}^{\rm 1h} (k, z) + P_{\nu \nu ^{\prime }}^{\rm 2h} (k, z) + P_{\nu \nu ^{\prime }}^{\rm shot} (k, z). \end{equation} \tag{ 18 }$$

Common to most halo models is that a distinction is made between central and satellite galaxies, with N^gal = N^cen + N^sat. All halos above a minimum mass M_min host a galaxy at their center,

$$\begin{equation} N^{\rm cen}(M) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}0 & M < M_{\rm min}, \\ 1 & M \ge M_{\rm min}, \end{array}\right. \end{equation} \tag{ 19 }$$

while any additional galaxies in the same halo would be designated as satellites which trace the dark matter density profile (e.g., Zheng et al. 2005). Halos host satellites when their mass exceeds the pivot mass M₁ (also known as M_sat in the literature), and the number of satellites is an exponential function of halo mass:

$$\begin{equation} N^{\rm sat}(M)=\left(\frac{M}{M_1} \right)^{\alpha }. \end{equation} \tag{ 20 }$$

In earlier halo models, galaxies were assumed to contribute equally to the emissivity density (e.g., Viero et al. 2009; Amblard et al. 2011; Planck Collaboration et al. 2011b). Therefore, assuming that the CIB originates from galaxies, spatial variations in the specific emission coefficient j_ν directly trace fluctuations in the galaxy number density

$$\begin{equation} \delta j_{\nu }/\bar{j}_{\nu } = \delta n^{\rm gal} / \bar{n}^{\rm gal}. \end{equation} \tag{ 21 }$$

The linear, 2-halo term dominates on large scales and is given by the clustering of galaxies in separate dark matter halos:

$$\begin{equation} P_{\nu \nu ^{\prime }}^{\rm 2h} (k, z) = \frac{1}{\bar{j}_{\nu }\bar{j}_{\nu }^{\prime }} P_{\rm lin}(k, z) D_{\nu }(k, z) D_{\nu }^{\prime }(k, z), \end{equation} \tag{ 22 }$$

with

$$\begin{eqnarray} D_{\nu } (k, z) & = & \int dM \frac{dN}{dM} (z) b(M, z) u_{\rm gal}(k, z, M) \nonumber \\ &&\times [ N_{\nu }^{\rm cen}(M, z) + N_{\nu }^{\rm sat}(M, z)], \end{eqnarray} \tag{ 23 }$$

where P_lin(k, z) is the linear dark matter power spectrum, b(M, z) is the linear large-scale bias, and u_gal(k, z, M) is the normalized Fourier transform of the galaxy density distribution within a halo, which is assumed to equal the dark matter density profile, i.e., u_gal(k, z, M) = u_DM(k, z, M). The nonlinear, 1-halo term dominates on small scales, and is written as

$$\begin{eqnarray} P_{\nu \nu ^{\prime }}^{\rm 1h} (k, z) & = &\frac{1}{ \bar{j}_{\nu } {\bar{j}_{\nu }^{\prime }}} \int dM \frac{dN}{dM} (z)\nonumber \\ && \times \left[N_{\nu }^{\rm cen}(M,z) N_{\nu ^{\prime }}^{\rm sat}(M,z) u_{\rm gal}(k, z, M)\right. \nonumber \\ &&+ N_{\nu ^{\prime }}^{\rm cen}(M,z) N_{\nu }^{\rm sat}(M,z) u_{\rm gal}(k, z, M) \nonumber \\ &&+\left. N_{\nu ^{\prime }}^{\rm sat}(M, z) N_{\nu }^{\rm sat}(M, z) u^2_{\rm gal}(k, z, M) \right ] \end{eqnarray} \tag{ 24 }$$

(Cooray & Sheth 2002), where dN/dM is the halo mass function.

However, conceptually it is wrong to assume that galaxies of different luminosities have equal weight in contributing to the power spectrum of the intensity fluctuations. A consequence of this assumption is that the excess signal on small angular scales from galaxies in massive halos can only be reproduced by having more satellite galaxies, leading to previous estimates of α which exceed predictions for subhalo indices from semi-analytic models (α ⩽ 1; e.g., Gao et al. 2004; Hansen et al. 2009) and a significant overabundance of satellites. For example, previous halo model fits to SPIRE data from Amblard et al. (2011) found α ∼ 1.7 at 250 μm and ∼1.8 at 350 and 500 μm, albeit with large errors.

A way to overcome this excess satellite problem, while still producing enough 1-halo power, is to have a model with fewer but more luminous satellites and weight galaxies by their luminosities. One such model is that of Shang et al. (2012), which invokes an L–M relation to tie the emissivity from galaxies to their host halo masses (also see, e.g., Yang et al. 2003; Vale & Ostriker 2004). The advantage of this model is that in principle one can predict the abundance as well as the clustering of galaxies observed at different frequency bands simultaneously, while in previous halo models it was impossible to predict the power spectrum across different frequency bands at the same time. As we will now show, this new formalism is similar to previous ones up until the numbers of central and satellite galaxies are substituted for luminosity-weighted quantities.

Hereafter, we follow the formalism of Shang et al. (2012), with some modifications. Novel to our implementation is the simultaneous fitting to the power spectra in all three bands, and to the number counts of sources above about 0.1 mJy (Glenn et al. 2010). In this implementation of the halo model, the mean comoving specific emission coefficient is

$$\begin{equation} \bar{j}_{\nu } (z) = \int dL \frac{dn}{dL}(L,z) \frac{L_{\nu }[(1+z)\nu ]}{4\pi }, \end{equation} \tag{ 25 }$$

where L is the luminosity and dn/dL is the luminosity function of DSFGs. What this means is that, unlike the earlier models, which assume that all galaxies contribute equally to the emissivity density (i.e., have the same luminosity), here the emissivities of galaxies in a given halo, by way of their luminosities, depend on the redshift, halo mass, and frequency:

$$\begin{equation} L_{\nu }[(1+z)\nu ] = L_0 (1+z)^{\eta } \Sigma (M)\Theta \left[ \left(1+z \right) \nu \right], \end{equation} \tag{ 26 }$$

where L₀ is the overall normalization factor, η describes the redshift evolution, Σ(M) describes the relation between infrared luminosity and halo mass (the L–M relation), and Θ(v) describes the shape of the infrared SED. Note that here M represents both the mass of the main halo and the infall mass of the subhalo.

Furthermore, the effective, luminosity-weighted number of central and satellite galaxies is

$$\begin{equation} f_{\nu }^{\rm cen} = N^{\rm cen} \frac{L_{(1+z)\nu }^{\rm cen}(M, z)}{4\pi }, \end{equation} \tag{ 27 }$$

$$\begin{equation} f_{\nu }^{\rm sat} = \int dm \frac{dn}{dm} (M, z) \frac{L_{(1+z)\nu }^{\rm sat}(m, z)}{4\pi }, \end{equation} \tag{ 28 }$$

where dn/dm(M, z) is the subhalo mass function of the main halo whose mass is M.

The terms in Equations (23) and (24) remain the same, except the numbers of central and satellite galaxies are substituted with their luminosity-weighted counterparts, such that D_ν in the 2-halo term (Equation (23)) becomes

$$\begin{eqnarray} D_{\nu } (k, z) & = & \int {dM} \frac{dN}{dM} (z) b(M, z) u_{\rm gal}(k, z, M) \nonumber \\ && \times [f_{\nu }^{\rm cen}(M, z) + f_{\nu }^{\rm sat}(M, z) ], \end{eqnarray} \tag{ 29 }$$

and the 1-halo term (Equation (24)) becomes

$$\begin{eqnarray} P_{\nu \nu ^{\prime }}^{\rm 1h} (k, z) & = &\frac{1}{\bar{j}_{\nu } \bar{j}_{\nu ^{\prime }}}\int dM \frac{dN}{dM} (z) \nonumber \\ && \times \left [ f_{\nu ^{\prime }}^{\rm cen}(M, z) f_{\nu }^{\rm sat}(M, z) u_{\rm gal}(k, z, M)\right. \nonumber \\ && +f_{\nu }^{\rm cen}(M, z) f_{\nu ^{\prime }}^{\rm sat}(M, z) u_{\rm gal}(k, z, M) \nonumber \\ &&\left.+f_{\nu ^{\prime }}^{\rm sat}(M, z) f_{\nu }^{\rm sat}(M, z)u^2_{\rm gal}(k, z, M) \right ]. \end{eqnarray} \tag{ 30 }$$

We define halos here as overdense regions whose mean density is 200 times the mean background density of the universe according to the spherical collapse model, and we adopt the density profile of Navarro et al. (1997) with the concentration parameter of Bullock et al. (2001), and the fitting function of Tinker et al. (2008) for the halo mass function and its associated prescription for the halo bias (Tinker et al. 2010). For the subhalo mass function, we use the fitting function of Tinker & Wetzel (2010). We will now describe each of the terms in Equation (26) in more detail. Note that using instead the concentration parameter of Duffy et al. (2008) leads to very little change in the final best-fit parameters.

5.2.1. (1 + z)^η: The Luminosity Evolution

5.2.2. Σ(M): The L–M Relation

Observationally it is clear that some halos are more efficient than others at hosting star formation (e.g., Béthermin et al. 2012b; Wang et al. 2012), and that the halo mass of most efficient star formation evolves with redshift (i.e., downsizing; e.g., Cowie et al. 1996; Bundy et al. 2006). It is also clear that star formation in halos is suppressed by several plausible mechanisms at the high-mass (e.g., accreting black holes; Birnboim & Dekel 2003; Kereš et al. 2005) and low-mass (e.g., feedback from supernovae, photoionization heating; Dekel & Silk 1986; Thoul & Weinberg 1996) extremes. Thus, following Shang et al. (2012), we assume that the L–M relation, Σ(m), can be parameterized by a simple log-normal distribution

$$\begin{equation} \Sigma (m) = m \frac{1}{\sqrt{2\pi \sigma ^2_{L/m}}} \exp \left[-\frac{(\log m - \log M_{\rm peak})^2}{2\sigma ^2_{L/m}}\right], \end{equation} \tag{ 31 }$$

where M_peak describes the peak of the specific IR emissivity per unit mass, and $\sigma ^2_{L/m}$ describes the range of halo masses in which galaxies producing IR emission reside. The minimum halo mass to host a galaxy, M_min, is left as a free parameter, but we place a lower limit on it such that L = 0 at M < M_min.

Note that we have implicitly assumed that the shape of the relation between halo mass and infrared luminosity is redshift-independent and identical for both central and satellite galaxies. Equation (25) can then be recast as

$$\begin{eqnarray} \bar{j}_{\nu } (z) &=& \int {dM} \frac{dN}{dM}(z) \frac{1}{4\pi } \left [ N^{\rm cen}L_{(1+z)\nu }^{\rm cen} \right. \nonumber \\ & &\left. + \int dm \frac{dn}{dm} (M, z) L_{(1+z)\nu }^{\rm sat} \right ], \end{eqnarray} \tag{ 32 }$$

where m is the subhalo mass at the time of accretion (e.g., Wetzel & White 2010; Shang et al. 2012) and dn/dm is the subhalo mass function in a host halo of mass M at a given redshift.

5.2.3. Θ(ν): The Model SED

Following Hall et al. (2010) and Shang et al. (2012), we begin with the simplest model SED, a single modified blackbody:

$$\begin{equation} \Theta (\nu) \propto \nu ^{\beta } B(\nu ,T_{\rm d}), \end{equation} \tag{ 33 }$$

where B(ν, T_d) is the blackbody spectrum (or Planck function), with effective dust temperature T_d, and β is the emissivity index. Both T_d and β are free variables with no redshift evolution. We refer to this as Model 1.

Next, in attempting to address the growing observational evidence for evolving temperature with redshift (e.g., Pascale et al. 2009; Amblard et al. 2010; Viero et al. 2013), we introduce an additional parameter, T_z, such that $T_{\rm dust} \propto (1+z)^{T_z}$. Note that the β parameter remains a free variable without redshift evolution in this model, which we refer to as Model 2.

Lastly, motivated by the findings of, e.g., Dunne & Eales (2001) or Elbaz et al. (2011)—who found that a typical DSFG spectrum is better fit by a linear combination of two SEDs: (1) those from hotter (T_{d, warm} ∼ 50 K), star-forming regions, and (2) those from colder (T_{d, cold} ∼ 20 K) regions of diffuse interstellar medium—we adopt a two-component SED. The ratio of the masses in the two components is defined as ξ =log(N_c/N_w), and is independent of redshift. Here β is redshift independent and fixed to equal 2. We refer to this model as Model 3.

5.2.4. Markov Chain Monte Carlo

The best-fit parameters for Models 1, 2, and 3 are tabulated in Tables 4, 5, and 6, respectively, and shown in Figure 11. The respective best fits to the counts are shown in Figure 12. The corresponding χ² (dof) are 368 (225), 357 (224), and 371 (223), or $\chi _{\rm reduced}^2=1.6$, 1.6, and 1.7, for Models 1, 2, and 3, respectively. That the addition of a parameter does not significantly improve the fits is discussed in Section 6.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Model 1: Best-fit Parameters and Corresponding Correlation Matrix

Parameter	log(Mmin)	log(Mpeak)	T	β	$\sigma ^2_{L/m}$	log(L0)	η
log(Mmin)	9.8 ± 0.5	0.16	−0.09	0.06	−0.15	−0.23	0.23
log(Mpeak)	...	12.2 ± 0.5	−0.11	0.16	−1.00	−0.42	0.49
T	...	...	23.1 ± 1.3	−0.97	0.10	0.74	−0.29
β	...	...	...	1.4 ± 0.1	−0.16	−0.62	0.25
$\sigma ^2_{L/m}$	...	...	...	...	0.4 ± 0.0	0.41	−0.50
log(L0)	...	...	...	...	...	−1.7 ± 0.1	−0.74
η	...	...	...	...	...	...	2.0 ± 0.1

Notes. Best-fit parameters and correlations between them in Model 1. Off-diagonal values of +1 or −1 mean that the parameters are highly correlated or anti-correlated, respectively, while values near 0 mean that they are independent of one another.

Download table as:  ASCIITypeset image

Table 5. Model 2: Best-fit Parameters and Corresponding Correlation Matrix

Parameter	log(Mmin)	log(Mpeak)	T	Tz	β	$\sigma ^2_{L/m}$	log(L0)	η
log(Mmin)	10.1 ± 0.5	−0.02	0.20	−0.27	−0.10	0.02	0.26	−0.25
log(Mpeak)	...	12.3 ± 0.5	0.21	−0.01	−0.23	−1.00	−0.18	0.20
T	...	...	20.7 ± 1.2	−0.66	−0.92	−0.21	0.76	−0.56
Tz	...	...	...	0.2 ± 0.0	0.38	0.01	−0.81	0.89
β	...	...	...	...	1.6 ± 0.1	0.23	−0.53	0.31
$\sigma ^2_{L/m}$	...	...	...	...	...	0.3 ± 0.0	0.18	−0.20
log(L0)	...	...	...	...	...	...	−1.8 ± 0.1	−0.90
η	...	...	...	...	...	...	...	2.4 ± 0.1

Notes. Best-fit parameters and correlations between them in Model 2. Off-diagonal values of +1 or −1 mean that the parameters are highly correlated or anti-correlated, respectively, while values near 0 mean that they are independent of one another.

Download table as:  ASCIITypeset image

Table 6. Model 3: Best-fit Parameters and Corresponding Correlation Matrix

Parameter	log(Mmin)	log(Mpeak)	Twarm	Tcold	ξ	Tz	$\sigma ^2_{L/m}$	log(L0)	η
log(Mmin)	10.1 ± 0.6	−0.39	0.40	0.40	−0.23	−0.43	0.39	0.42	−0.41
log(Mpeak)	...	12.1 ± 0.5	−0.75	−0.91	0.02	0.89	−1.00	−0.79	0.90
Twarm	...	...	26.6 ± 2.8	0.80	−0.05	−0.90	0.75	0.94	−0.90
Tcold	...	...	...	14.2 ± 1.0	0.05	−0.93	0.90	0.80	−0.92
ξ	...	...	...	...	1.8 ± 0.1	0.18	−0.02	−0.27	0.19
Tz	...	...	...	...	...	0.4 ± 0.1	−0.88	−0.95	0.99
$\sigma ^2_{L/m}$	...	...	...	...	...	...	0.4 ± 0.0	0.79	−0.89
log(L0)	...	...	...	...	...	...	...	−1.9 ± 0.1	−0.96
η	...	...	...	...	...	...	...	...	2.7 ± 0.2

Notes. Best-fit parameters and correlations between them in Model 3. Here ξ is the ratio of the masses of the cold and warm components. Off-diagonal values of +1 or −1 mean that the parameters are highly correlated or anti-correlated, respectively, while values near 0 mean that they are independent of one another.

Download table as:  ASCIITypeset image

Correlations between parameters are presented in the off-diagonal entries of Tables 4, 5, and 6 for Models 1, 2, and 3, respectively. As pointed out by, e.g., Shang et al. (2012), certain pairs of parameters exhibit very high levels of correlation: in particular, the SED parameters (β and T in Models 1 and 2, T_cold and T_warm in Model 3), the luminosity parameters η and L₀, and the mass parameters $\sigma ^2_{L/m}$.

We find that star formation is most efficient in halos ranging from log(M_☉) = 11.7 to 12.5, peaking at log(M_☉) ∼ 12.1, which is consistent with several recent results from observations and simulations (e.g., Moster et al. 2010; Wang et al. 2012; Behroozi et al. 2013). In Cooray et al. (2010), where a halo model is developed to fit the angular correlation functions of galaxies brighter than 30 mJy, the minimum halo mass scale is log (M_min/M_☉) = 12.6, 12.9, and 13.5 at 250, 350, and 500 μm, respectively. These values are much higher than our best-fit log (M_min/M_☉) = 10.1 ± 0.6, which is due to the fact that faint galaxies (around 5 mJy) dominate the power spectrum of the intensity fluctuations (see Figure 1). In Amblard et al. (2011), the minimum halo mass scale is log (M_min/M_☉) = 11.1, 11.5, and 11.8 at 250, 350, and 500 μm, respectively, which are higher than our best-fit value. Furthermore, the evolution of the dust temperature, characterized by T_z, is in very good agreement with stacking measurements found by, e.g., Pascale et al. (2009) and Viero et al. (2013).

The redshift distribution of the emissivity—which in previous halo models has been either parameterized or adopted from galaxy population models—is here an output of the L–M relation, and shown compared to a selection of models (Valiante et al. 2009; Béthermin et al. 2011, 2012a) and previous estimates (Amblard et al. 2011) in Figure 13.

Zoom In Zoom Out Reset image size

Lastly, we plot the absolute CIB in each band output by our model, along with several measurements from the literature, in Figure 14. We find that our models are consistent with the fiducial FIRAS values (Fixsen et al. 1998; Lagache et al. 2000), though we note that the uncertainties in the fiducial measurements are of order 30%.

Zoom In Zoom Out Reset image size

Since the 1-halo term originates from multiple galaxies occupying the same halo, the reduction of the 1-halo term with more aggressive source masking suggests that some fraction of the more luminous resolved SPIRE sources are satellites in massive halos or cluster members, though it should be noted that some fraction of bright galaxies could be close pairs, which would make the same imprint on the spectra. This is consistent with the interpretation from clustering measurements of resolved SPIRE sources which claims that ∼14% of sources with S₂₅₀ > 30 mJy appear as satellites (Cooray et al. 2010).

One may wonder if the reduction in power can be solely attributed to a reduction in the Poisson level, but this is unlikely. We show this by fitting the unmasked 250 μm spectra with a 2-halo, 1-halo, and Poisson term, and then fitting the same terms to spectra with sources above 50 mJy masked. The 2-halo term is fixed in both for both levels of masking, while in the latter fit the 1-halo term is first fixed and then allowed to float, in order that the two fits can be compared. For a floating 1-halo term, χ² = 10.0 for 15 dof ($\chi _{\rm reduced}^2=0.7$), while for a fixed 1-halo term χ² = 98.3 for 16 dof ($\chi _{\rm reduced}^2=6.1$), thus ruling out the possibility that the 1-halo term has not been affected by masking.

That some fraction of more luminous sources are found in satellite halos is also expected from semi-analytic models of DSFGs. For example, González et al. (2011) predict that 38% of all DSFGs (defined there as S ≳ 1 mJy at 850 μm) and 24% of the most luminous sources (S ≳ 5 mJy at 850 μm) are satellites. Furthermore, observations of individual groups and clusters universally find that when star-forming galaxies are present they are located on the outskirts of the massive halos (typically defined as the volume between R₅₀₀ and R₂₀₀; e.g., Saintonge et al. 2008; Tran et al. 2009; Braglia et al. 2011). Even more support for this scenario comes from stacking in the submillimeter at positions of brightest cluster galaxies, which shows a bump in emission at 0.8 Mpc from the central galaxy (Coppin et al. 2011). In light of these observations, the reduction of 1-halo power with masking is unsurprising.

Released to the archive at about the same time, initial CIBA power spectra from Herschel/SPIRE (Amblard et al. 2011) and Planck/HFI (Planck Collaboration et al. 2011b) were found to be discrepant by more than 15%. The updated, published Planck Collaboration et al. (2011b) paper explored this difference in detail by comparing to power spectra of SPIRE maps that have had no masking applied (rather than the published spectra which masked all pixels greater than 50 mJy). They found that a discrepancy still remained, with SPIRE measurements lower by factors of ∼1.7 and ∼1.2 at 857 and 545 GHz (350 and 550 μm) over the angular range 0.02 < k_θ/arcmin < 0.07 (400 < ℓ < 1500). If we compare our spectra, having only masked extended sources, to the Planck/HFI spectra over the same angular range, we find our spectra are also low, but by factors of ∼1.3 ± 0.06 and 1.4 ± 0.06 at 350 and 500 μm, respectively.

They next proposed a way to resolve the discrepancy. They showed that the spectra of the two groups could be brought into agreement if: (1) no Galactic cirrus is removed from the SPIRE data; and (2) the beam surface area—which they claim to be overestimated by ∼4% and 9% at 350 and 500 μm, respectively—is corrected. Indeed, we also find that the previous cirrus values were overestimated; and from a careful estimate of the beam area (described in detail in Appendix B), we also find a correction, though equaling ∼2% and 8%. These corrections have brought the overall offset, particularly on scales k_θ ≲ 0.04 arcmin⁻¹ (ℓ ≲ 900), into agreement; however, on small scales they have been unable to close the gap entirely. The remaining offset on the largest scales can be attributed to systematic calibration uncertainties of the two instruments: 7% for SPIRE (Appendix B), and 7% for Planck/HFI (Planck Collaboration et al. 2011a), which may be a product of the very different calibration strategies of the two instruments, as well as potential systematic uncertainties from cirrus removal. But it is on scales k_θ ≳ 0.04 arcmin⁻¹ (ℓ ≳ 900)—scales on which the contribution from Galactic cirrus is negligible—that excess powers in the Planck curves are still either in tension (350 μm) or do not agree (500 μm).

Ultimately, since this paper first appeared it was determined that the Planck calibration of the 857 and 545 GHz (350 and 550 μm, respectively) channels was off by 7% and 15%, meaning that to properly compare the data their spectra must be divided by 1.14 and 1.30 (see Planck Collaboration et al. 2013). This correction—which has been applied to the data in Figure 7—indeed brings the two curves into generally good agreement.

Our halo models represent a step forward by being the first to fit the auto- and cross-frequency power spectra, number counts, and absolute CIB levels simultaneously. The added complexity introduced by these models is justified as simpler halo models have been unable to simultaneously fit published power spectra (e.g., Planck Collaboration et al. 2011b). Yet, considering the wide range of data fit, our models remain relatively simple; e.g., they have a comparable number of parameters to the model presented in Amblard et al. (2011, which varied the Poisson level, M_min, M₁, α, and 4 dS/dz nodes, but fit each spectra independently and had no L–M relation). However, with reduced χ² of ∼1.6, our models cannot formally be claimed to be good fits, and in fact there are some obvious problems with these models.

Firstly, the addition of parameters only marginally improves the fits, if at all. This is because the χ² is dominated by the poor fit to the counts, particularly at the bright end (i.e., S ≳ 20 mJy). Considering that the clustering power is dominated by the sources below that flux density level (Figure 1), this tension between the fits to the counts and the fits to the spectra is not terribly surprising. It may be that this tension is attributable to there being two types of submillimeter-emitting galaxies (i.e., hotter "starburst" type galaxies, and more "normal main sequence" galaxies; e.g., Elbaz et al. 2011). Though we explore the possibility of two populations in Model 3, we do not decouple their contributions in redshift. If local starbursts are indeed responsible for a separate, lower redshift nonlinear term, then our model would struggle to satisfy the linear and nonlinear components at all redshifts, as it appears to do.

Second, there are significant degeneracies, as well as possibly questionable assumptions, built into this model. The strong coupling of the SED parameters T, β, and temperature evolution T_z renders the resulting best-fit parameters difficult to meaningfully interpret. Also, the very strong anti-correlation of luminosity evolution η and normalization L₀ likewise makes interpreting the evolution of the galaxy luminosity or galaxy bias difficult. Furthermore, the assumption that the luminosity evolution increases to z = 2 and then abruptly flattens, though seen in several observations (e.g., Stark et al. 2009; González et al. 2010), is likely to be extreme and prone to galaxy selection effects (e.g., Weinmann et al. 2011). Getting this wrong would result in compensating for the discrepant high-z power by the other parameters in unpredictable ways. Lastly, the observed quenching of star formation in the cores of the most massive halos (e.g., Cattaneo et al. 2006) is not treated by the model, but may have a significant impact on the nonlinear component of the power spectrum.

Thus, our model is far from the final say in the interpretation of CIBA, as it appears that the quality of the data demands a model with additional levels of sophistication. Future models can address many of these limitations by carefully implementing observations. For example, strong constraints on the stellar-mass to halo-mass relationship (e.g., Behroozi et al. 2013), star-formation to halo-mass relationship (e.g., Wang et al. 2012), and infrared luminosity to stellar-mass relationship (Viero et al. 2013) now exist. That, combined with knowledge of the quiescent fraction of galaxies with redshift (e.g., Quadri et al. 2012), would result in more constrained models with fewer parameters. Furthermore, future models could be extended to fit not only the three SPIRE bands, but also longer wavelength data from, e.g., ACT, Planck, and SPT, as well as at shorter wavelengths from, e.g., Spitzer, IRAS, and WISE. By combining long-wavelength multi-band studies of map-based power spectra with discrete object correlations at shorter wavelengths, we should be able to build a much more complete picture of the relationship between stars, star formation, and dark matter halos.