BLAST: CORRELATIONS IN THE COSMIC FAR-INFRARED BACKGROUND AT 250, 350, AND 500 μm REVEAL CLUSTERING OF STAR-FORMING GALAXIES

With the discovery of the cosmic far-infrared background (CIB; Puget et al. 1996; Fixsen et al. 1998) and subsequent studies from the mid-infrared to the millimeter, it has been established that the peak epoch of star formation lies between 1 ≲ z ≲ 3 (Dickinson et al. 2003; Hopkins 2004).

Where star formation occurs with respect to the underlying dark matter distribution is less well understood. In the local universe, the sites of most active star formation occur far from the densest environments, a property which reverses at z ∼ 1 (Elbaz et al. 2007). This points to the environment contributing to the mechanisms that trigger or quench star formation.

Measurements of the clustering of star-forming galaxies on large and small scales can be used to directly relate star formation to the environment. In the regime of linear growth of structure, the clustering amplitude relative to that of the underlying dark matter is described by a bias parameter, b, which describes how strongly star formation traces the underlying dark matter. On smaller scales, the distribution of galaxies hosted by a dark matter halo is described by the halo occupation distribution (HOD; Peacock & Smith 2000). It contains information on the abundance of star-forming sources within individual halos. Targeting the submillimeter band is particularly efficient for observing star formation directly. The submillimeter background results from the thermal emission of interstellar dust in high-redshift star-forming galaxies, which is heated by optical and ultraviolet radiation from stars and to a lesser extent active galactic nuclei (see Blain et al. 2002). A substantial effort has been devoted to surveys of these galaxies (e.g., Smail et al. 1997; Hughes et al. 1998; Barger et al. 1998; Borys et al. 2003); thus, in principle, a clustering signal could be measured from these sources directly.

However, direct measurement of the clustering properties of resolved submillimeter galaxies has been elusive. Due to limited mapping speeds, the areal coverage of even the most ambitious submillimeter surveys has been relatively small (e.g., SHADES mapped approximately a quarter square degree at 850 μm; Coppin et al. 2006). In addition, due to steeply falling counts and modest resolutions of single-dish submillimeter telescopes, source confusion has made it difficult to resolve any sources other than those with very high signal-to-noise ratio. These sources span a relatively wide redshift range, roughly 1 ⩽ z ⩽ 4, peaking at z ∼ 2.4 (Chapman et al. 2005), which has the effect of washing out the angular clustering signal, further complicating measurements. This difficulty was confirmed by Scott et al. (2006), who re-analyzed all the SCUBA fields and found tentative evidence of strong angular clustering, but with errors too large to adequately constrain the spatial correlation length. Furthermore, relatively large beams have made it difficult to efficiently identify direct counterparts (Barger et al. 1999; Ivison et al. 2000) in order to obtain spectroscopic or photometric redshifts. Blain et al. (2004) attempted to measure the spatial clustering properties combining 73 sources with spectroscopic information, but again were only able to tentatively measure the clustering length. Other SCUBA galaxy clustering measurements, some tentatively detecting clustering, have been made by Webb et al. (2003) and Blake et al. (2006); meanwhile Almaini et al. (2003) claim to have found evidence for strong angular clustering between X-ray and submillimeter populations, although Borys et al. (2004) were unable to confirm the result using another, seemingly less biased, estimator. Using a nearest neighbor analysis, Greve et al. (2004) found that most significant MAMBO (1.2 mm) sources come in pairs, separated by ∼23''. Put together, these limitations have made it extremely difficult to measure the clustering signal robustly.

To study the clustering properties of submillimeter galaxies, it is more powerful to consider the statistics of the unresolved CIB, which contains the full intensity, rather than working from a limited catalog. In other words, instead of measuring a correlation among numbers of galaxies, use the background fluctuations of the total intensity to measure correlations among brightnesses of galaxies. Devlin et al. (2009) and Marsden et al. (2009) have demonstrated that the CIB is composed of emission by discrete sources. Since submillimeter galaxies are optically thin, this signal will be proportional to the total star formation rates of those sources. Correlations in the CIB will have a contribution in excess of white noise—which arises from Poisson sampling of a background made up of a finite number of sources—in the presence of clustering, with an amplitude that should be detectable with current surveys (Scott & White 1999; Haiman & Knox 2000; Knox et al. 2001; Magliocchetti et al. 2001; Perrotta et al. 2003; Amblard & Cooray 2007; Negrello et al. 2007). Initial attempts to detect correlations by Peacock et al. (2000), for the Hubble Deep Field observed by SCUBA at 850 μm, and by Lagache & Puget (2000) for a 0.25 deg² Infrared Space Observatory (ISO) field at 170 μm, were only able to measure a signal consistent with the Poisson contribution. More recently, Grossan & Smoot (2007) and Lagache et al. (2007) reported the weak detection of a clustering component in 160 μm data from ∼9 deg² Spitzer fields.

In this paper, we report the detection of correlations in the submillimeter part of the CIB due to the clustering of star-forming galaxies, in a 6 deg² field centered on the Great Observatories Origins Deep Survey South field (GOODS-South; Giavalisco et al. 2004). These data were collected by the Balloon-borne Large Aperture Submillimeter Telescope (BLAST; Devlin et al. 2009), which is designed to bracket the peak of redshifted thermal emission from dust by observing at 250, 350, and 500 μm. Operating above most of the atmosphere, BLAST is able to make observations in bands which are difficult or impossible to observe from the ground. A detailed description of the instrument and calibration can be found in Pascale et al. (2008) and Truch et al. (2009).

This paper is organized as follows. In Part I, we describe how we make the measurement—from map preparation to power spectrum calculation. We address each contribution to the total power spectrum, and how they are removed to uncover the clustering signal. We show that the observed spectra are consistent across BLAST bands and correspond to a modest bias. In Part II, we assess the plausibility of the detection by fitting models to the data: beginning with a simple linear bias model, followed by a more detailed halo model (Mo & White 1996; Cooray & Sheth 2002) and HOD (Peacock & Smith 2000). When required we adopt the concordance model, a flat ΛCDM cosmology with Ω_M = 0.274, Ω_Λ = 0.726, H₀ = 70.5 km s⁻¹ Mpc⁻¹, and σ₈ = 0.81 (Hinshaw et al. 2009).

Galaxy clustering can be expressed in a number of ways, the most common being the two-point correlation function, w(θ), which measures the number of pairs at a given distance in excess of what would be expected of a Poisson distribution. Alternatively, the clustering of galaxies can be expressed as a power spectrum in excess of Poisson noise, $\tilde{P}(k_{\theta })$, which can be expressed in dimensionless units as the fractional variance per logarithmic increment of wavenumber (see Peacock 1999), i.e.,

$$\begin{equation} \tilde{\Delta }_{k_{\theta }}^2 \equiv 2\pi k_{\theta }^2 \tilde{P}(k_{\theta }), \end{equation} \tag{ 1 }$$

where k_θ is the angular wavenumber, which is also known as σ in the literature, and is expressed in inverse angular scale as k_θ = 1/λ. It is related to the multipole index, ℓ, by ℓ = 2πk_θ. The tildes over $\Delta _{k_{\theta }}^2$ and P(k_θ) denote that these quantities refer to galaxy locations rather than intensities.

Naturally, the correlation function and the power spectrum of galaxy clustering are related; they form a Hankel transform pair. Explicitly, for small surveys,

$$\begin{equation} \begin{array}{c} \dsty \tilde{\Delta }^2_{k_{\theta }} = (2\pi k_{\theta })^2 \int ^{\infty } _0 w (\theta) J_0 (2\pi k_{\theta } \theta) \theta \ d\theta, \\ \ms \dsty w (\theta) = \int ^{\infty } _0 \tilde{\Delta }^2_{k_{\theta }} J_0 (2 \pi k_{\theta } \theta)\ dk_{\theta } / k_{\theta }.\end{array} \end{equation} \tag{ 2 }$$

For small separation correlations of local galaxies, angular clustering is often described as a power law, w(θ) = (θ/θ₀)⁻, where ≃ 0.8 is the canonical slope (e.g., Giavalisco et al. 1998). The power spectrum analog for small areas (see Peacock et al. 2000) is

$$\begin{equation} \tilde{\Delta }^2_{k _{\theta }} (k_{\theta }) = (2\pi k_{\theta } \theta _0)^{\epsilon } 2^{1-\epsilon } \frac{\Gamma (1- \epsilon /2)}{\Gamma (\epsilon /2)}, \end{equation} \tag{ 3 }$$

which is equal to 3.35(k_θθ₀)^0.8 for ≃ 0.8.

As previously stated, this holds for correlations of galaxy locations. On the other hand, the power spectrum of background fluctuations, which is what we calculate, measures correlations of galaxy intensities. The redshift distribution of the cumulative flux contributed by the background sources is represented as

$$\begin{equation} \frac{d{\cal S}}{dz} = \int _{0}^{\infty }S\frac{dN}{dS\ dz}(S,z)\ dS, \end{equation} \tag{ 4 }$$

where dN/(dSdz) is the number density of sources per unit flux density and redshift interval. The measured power spectrum of background fluctuations, P(k_θ), is the two-dimensional, flux weighted, projection of the three-dimensional galaxy clustering spectrum, P_3D(k_θ). For k_θ ≫ 1 and a flat cosmology, their relationship can be approximated as

$$\begin{equation} P(k_{\theta }) = \int _{z_{\rm min}}^{z_{\rm max}}P_{\rm 3D}(2\pi k_{\theta }/x(z),z)\left(\frac{d{\cal S}}{dz}(z)\right)^{2}\frac{1}{dV_{\mathrm{c}}(z)}\ dz, \normalsize \end{equation} \tag{ 5 }$$

(e.g., Tegmark et al. 2002), where x(z) is the comoving radial distance and dV_c(z) is the comoving volume element, i.e., dV_c(z) = x(z)²dx/dz (assuming a flat cosmology), where dN/(dSdz) is the number density of sources per unit flux density and redshift interval.

P(k_θ) can expressed in dimensionless units as

$$\begin{equation} \Delta _{k_{\theta }}^2 \equiv 2\pi k_{\theta }^2 P(k_{\theta })\big/I^2_{\nu }, \end{equation} \tag{ 6 }$$

where I_ν is the intensity of the background in Jy. Although $\Delta _{k_{\theta }}^2$ and $\tilde{\Delta }_{k_{\theta }}^2$ have the same units and similar meaning, they differ because Equation 6 deals with an intensity weighted power spectrum.

In addition to the clustering signal, the total power spectrum has contributions from instrumental noise, Poisson noise from individual background galaxies, and cirrus emission. We will address each contribution individually.

3.1. Map Preparation

BLAST observed a wide 8.7 deg² patch, centered on the GOODS-South field (3^h32^m35^s, − 28°15'; hereafter BGS-Wide), with mean 1σ sensitivities of 36, 31, and 20 mJy beam⁻¹ at 250, 350, and 500 μm, as well as a deep, nested field of 0.8 deg², centered on (3^h32^m30^s, − 27°48'; hereafter BGS-Deep) with mean 1σ sensitivities of 11, 9, and 6 mJy beam⁻¹, respectively.¹⁶ A 6 deg² region, centered on BGS-Wide, is selected from the map, because of its uniformity in observed depth. The BLAST bolometers are prone to drifts on timescales greater than ∼10 s. To retain as much large angular scale signal as possible, a fast scan rate is preferred because larger scales will survive the high-pass filtering, designed to remove noise below the $0.1\ \rm Hz$ 1/f knee. For this analysis, we use only the data from the wide region of the map, where the scan rate is 0.1 deg s⁻¹, and the rms of the maps is 2.4, 1.9, 1.0 MJy sr⁻¹, at 250, 350, and 500 μm, which is applicable to calculating uncertainties of point-source flux densities. The data for the nested deep region, whose scan rate is only 0.05 deg s⁻¹, are not included. Large-scale noise is removed by high-pass filtering the time-streams at 0.2 Hz. Correlated noise—a drift of multiple detectors in unison—is not removed because it would inevitably suppress large-scale signal as well. Instead, a cross-correlation of a subset of the maps is used to remove large-scale correlated noise (described below). Although fully optimized map-makers are available (e.g., SANEPIC; Patanchon et al. 2008), the long time to convergence (typically 24 hr on 10 processors for this particular map) makes it impractical to use them for our Monte Carlo simulations. Full analysis with SANEPIC maps, including the nested deep field, will be the subject of a future paper. Here, we instead use OPTBIN (Pascale et al. 2009) a fast, naive, map-maker whose transfer function is calculated using a Monte Carlo simulation (see Section 3.2 and Figure 1 for details), and we use SANEPIC maps for consistency checks. Parallel power spectrum analysis with both map-makers show agreement well within the simulated uncertainties.

Zoom In Zoom Out Reset image size

3.2. Power Spectrum Calculation

The map intensity, S_map, can be written as

$$\begin{equation} S_{\mathrm{map}} = (T\otimes [S_{\mathrm{sky}}\otimes B + N]) W, \end{equation} \tag{ 7 }$$

where we use ⊗ to represent a convolution, S_sky is the true sky surface brightness, Tis the transfer function of the map-maker, B is the measured instrumental beam, N is the instrumental noise, and W is the "aperture function," which is zero beyond the region of interest. The autocorrelation of a map will contain a contribution from detector noise. To suppress this instrumental noise, cross-power spectra are taken among a set of four maps which are made by dividing the time-stream into four roughly equal parts and then making four separate maps (hereafter referred to as sub-maps). The time-streams, which are made up of numerous chunks, are divided into every fourth chunk (e.g., 1, 5, 9,..., and 2, 6, 10,..., etc.), so that the sub-maps have as similar coverage as possible. The number of sub-maps chosen maximizes the number that can be made while maintaining uniformity in hits, retaining some cross-linking, and avoiding holes in the maps. The rms of the resulting sub-maps is 4.6, 3.6, and 2.0 MJy sr⁻¹, at 250, 350, and 500 μm. In the cross-spectrum, noise which is uncorrelated between sub-maps averages to zero. Consequently, the spectrum does not depend on modeling the potentially complicated or non-stationary noise.

We prepare the maps before calculating the power spectrum by removing their means, apodizing them with a Welch window (Press 2002, Chapter 13.4), and zero-padding them with a width on each side equal to half the map. The cross-correlation two-dimensional power spectrum of each pair of maps is calculated. The azimuthal average of the amplitudes (which in two-dimensional k-space appears to be isotropic) is taken to find the one-dimensional power spectrum, P(k). The resulting spectra are averaged and divided by the power spectrum of the beam and the transfer function. The transfer function is calculated with a Monte Carlo simulation from simulated maps made with the BLAST simulator,¹⁷ and is shown in Figure 1. At angular scales, which are large compared to the high-pass filter (<0.03 arcmin⁻¹) and approaching the cutoff of the beam (>0.9 arcmin⁻¹), the transfer function is unreliable. We are interested in the scales bracketed by these limits.

3.2.1. Jackknife Tests

We have performed jackknife tests in which a correlation is calculated between two distinct difference maps. Specifically, we take the cross-correlation of: (sub − map 1 − sub − map 2) and (sub − map 3 − sub − map 4). If the cross-correlation measures only signal, then taking a difference between sub-maps should cancel sky signal and results in a cross-correlation power spectrum consistent with zero. And indeed, our results are consistent with zero, in the range of interest, in all three bands.

As a further check, the cross-power spectrum is compared to the difference of the auto-power spectrum (of a map made from the entire time-stream) and the measured noise. The noise is estimated from odd–even pixel jackknife maps, and is approximately white over the angular scales of interest. The resulting two spectra are in excellent agreement.

3.2.2. Poisson Noise

Poisson (or shot) noise arises from the finite number of galaxies per unit area. It can be calculated analytically from the source counts as

$$\begin{equation} P_{\mathrm{shot}} = \int _0 ^{\infty } S^2 \frac{dN}{dS} dS. \end{equation} \tag{ 8 }$$

Alternatively, Poisson noise levels can be estimated from Monte Carlo simulations using mock-maps that are populated with uncorrelated sources whose fluxes are drawn from the measured BLAST counts (Patanchon et al. 2009). This has the added advantage that realistic uncertainties can be calculated as well. Since the counts are so steep, care must be taken to reproduce the bright end of the counts faithfully, so that extremely rare bright sources do not unrealistically appear in the realizations. We find that the two methods of estimating the level of shot noise agree within the 1σ uncertainties.

Due to the steep nature of the source counts, the Poisson noise is dominated by the contributions of the fainter population. Furthermore, Marsden et al. (2009) find only 15% of the total sky intensity is associated with a 3σ catalog. We find that removal of only the brightest sources results in an approximately 5% reduction in Poisson noise at 250 μm, and that more aggressive cuts lead to removal of correlations beyond just Poisson noise. We find behavior consistent with this in our simulations. Therefore, we subtract only five sources above 500 mJy at 250 μm, two sources above 400 mJy at 350 μm, and no sources at 500 μm. To subtract sources, first we make a source list by performing a noise-weighted convolution of the maps with the effective BLAST point-spread function (PSF) and identify local maxima in the smoothed map. We then subtract a scaled effective PSF with amplitude taken from the source list. We perform the same operation on our mock-maps, from which we calculate Poisson levels of 11.4 ± 1.0 × 10³, 6.3 ± 0.5 × 10³, and 2.7 ± 0.2 × 10³Jy²sr⁻¹ at 250, 350, and 500 μm. These are shown as dashed lines in Figure 2.

Zoom In Zoom Out Reset image size

3.2.3. Estimates of Galactic Foregrounds

Gautier et al. (1992) show that the power spectrum of Galactic cirrus can be approximated by a power law,

$$\begin{equation} P_{\mathrm{cirrus}}(k_{\theta }) = P_0\left(\frac{k_{\theta }}{k_0} \right) ^{-\alpha }, \end{equation} \tag{ 9 }$$

where k_θ is the angular wavenumber in inverse arcminutes, and P₀ is the power spectrum value at k₀ = 0.01 arcmin⁻¹. We measure the cirrus component at 100 μm from the cross-correlation of the three co-added IRIS (HCON¹⁸) maps (reprocessed IRAS: Miville-Deschênes & Lagache 2005) for a ∼15 deg²region surrounding the BGS-Wide field. The amplitude of the observed power spectrum has a contribution from cirrus emission, which is highly variable on the sky. The GOODS region was specifically chosen because it is a low cirrus region, with a mean intensity of 1.39 ± 0.18 MJy sr^-1. At scales 0.008 < k_θ < 0.03 arcmin⁻¹, which is larger than the scales probed by BLAST, we find that the cirrus is well approximated by a power law, P₀ = (0.47 ± 0.18) × 10⁶ Jy²sr⁻¹and α = 2.91 ± 0.11, values which are similar to those found in low cirrus emission regions measured by Lagache et al. (2007) and Miville-Deschênes et al. (2007). We assume that the power spectrum continues to smaller angular scales, and scale it to the BLAST bands using the average dust emission color (I_BLAST/I₁₀₀)², which is found by assuming cirrus emission behaves as a modified blackbody, ν^βB(ν), where B(ν) is the Planck function, and β is the emissivity index (Draine & Lee 1984). The scaled power spectrum and errors are calculated with a Monte Carlo simulation, varying temperature (17.5 ± 1.5 K) and β (1.9 ± 0.2) (Boulanger et al. 1996). The resulting power-law approximations and uncertainties are illustrated as dotted lines in Figure 2, which show that we are negligibly affected by Galactic cirrus on all recovered angular scales.

3.2.4. Uncertainties

To estimate uncertainties in our angular power spectra, clustered signal plus noise simulated maps are analyzed with the same pipeline as the astronomical data. We include clustering in the simulations in case the amplitude or shape of the transfer function depends on the input, but we found that not to be the case. We follow Almaini et al. (2005, see appendix for algorithm) to introduce correlations to the simulated maps, with an angular correlation length, θ₀ = 30. This angle was chosen from the measured upper limit of Peacock et al. (2000). Realizations with stronger clustering clearly do not look like BLAST maps.

Figure 2 shows an unambiguous signal in excess of Poisson noise on scales of 0.04–0.2 arcmin⁻¹, which cannot be explained by Galactic cirrus, and which we interpret as correlations from clustered star-forming galaxies. The vertical dashed line indicates the angle at which the distribution of 24 μm selected FIDEL galaxies (Magnelli et al. 2009) begin to show variance in excess of Poisson. Indeed, the similarity to Figure 3 of Marsden et al. (2009) is striking.

The CIB has an amplitude I^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 (Marsden et al. 2009). Figure 3 shows the clustering component of the power spectrum normalized by I^CIB_ν in the three BLAST bands, where the 250 μm and 500 μm data have been displaced slightly for visual clarity. The contributions of cirrus and Poisson noise have been subtracted, and the data rebinned in logarithmic intervals. These results are also listed in Table 1. We use the BLAST estimates for the CIB because they are the most precise estimate available of the CIB in these wavelength bands. Doing so has the additional benefit in this case that calibration uncertainties completely vanish in the ratio. The fit to a single power law, and also the agreement in amplitude across the BLAST bands, once the power spectra are normalized to the sky intensity, are excellent.

Zoom In Zoom Out Reset image size

Table 1. CIB Normalized Clustering Power Spectra, P(k_θ)/I²_ν

kθ	BLAST 250	BLAST 350	BLAST 500
$(\rm arcmin^{-1})$	$(\rm steradian^{-1})$	$(\rm steradian^{-1})$	$(\rm steradian^{-1})$
0.044	(1.39 ± 0.52) × 10−7	(1.33 ± 0.47) × 10−7	(2.47 ± 1.17) × 10−7
0.070	(5.94 ± 2.02) × 10−8	(4.98 ± 1.83) × 10−8	(5.71 ± 2.06) × 10−8
0.113	(2.99 ± 0.95) × 10−8	(2.62 ± 0.90) × 10−8	(2.98 ± 1.00) × 10−8
0.183	(9.81 ± 7.51) × 10−9	(9.94 ± 6.92) × 10−9	(6.91 ± 4.62) × 10−9

Notes. The errors do not include uncertainties in the CIB. CIB values are listed in Section 4.

Download table as:  ASCIITypeset image

The relative variance of the CIB, $\Delta ^2_{k_{\theta }}$, formed by dividing P(k_θ) by 2πk²_θ, is shown in the bottom panel of Figure 3, where the 250 μm and 500 μm data have been displaced slightly for visual clarity. (Amplitudes have been additionally converted from $\rm arcmin^{-2}$ to $\rm steradian^{-1}$ in the figure and in the discussion below.) The best fits of Equation (3) are shown in Table 2. The data are consistent (to within 1σ) with Δ²_k = Constant, with the same amplitude in all three bands, even taking into account that the error bars shown contain a large component of uncertainty which is common mode between channels, and thus overestimate the anticipated scatter. The best-fit amplitude for a power law with slope of −2 is shown as a dotted line in the bottom panel, corresponding to $\Delta _{k_{\theta }} = \delta I/I= 15.1\% \pm 1.7\%$. This is directly analogous to the square root of the "band power" measured in cosmic microwave background (CMB) anisotropy experiments, where typically δT/T ∼ 10⁻⁵. From this point of view the CIB is much clumpier than the CMB, as one would expect since the galaxies are observed after a much longer period of linear growth. The dotted line in the upper panel of Figure 3 is $2\pi k_{\theta }^2\Delta ^2_{k_{\theta }}$, corresponding exactly to the fit in the lower panel.

Table 2. Best-fit Values Obtained for θ₀ and

BAND	θ0 (arcsec)
250 μm	0.017 ± 0.020	0.27 ± 0.19
350 μm	0.008 ± 0.006	0.26 ± 0.19
500 μm	0.0	0.0

Download table as:  ASCIITypeset image

The values of the parameters and θ₀ (see Equation (3)) that best fit the data are shown in Table 2. The large uncertainty quoted for θ₀ is partly due to the awkwardness of this parameterization near zero slope. The bottom panel of Figure 3 shows the best fit for 250 μm data as a dashed line. Such a small θ₀, despite significant power in excess of Poisson noise, requires that the sources which make up the CIB are distributed over a wide range of redshifts, and has implications for future clustering measurement strategies (e.g., with Herschel or Planck), which we discuss in Section 11.

We have calculated the correlations between different bands (e.g., 250 × 350) using the same pipeline and set of sub-maps. The cross-spectra are normalized by the square root of the auto-spectra of the two bands, so that the final curve would be unity at all scales for identical maps, and zero at all scales for two completely different maps. Results show that the cross-correlations are 0.95 ± 0.06, 1.06 ± 0.09, and 0.92 ± 0.04, for 250 × 350, 350 × 500, and 250 × 500, respectively, over the range of angular scales 0.04 < k_θ < 0.5 arcmin⁻¹. Neighboring bands are more correlated with each other than are the 250 and 500 μm bands. While we find the same spectrum in all three bands, the phases are different, as one would expect if the three BLAST bands have different selection functions and sample the galaxies in the CIB at different redshifts. While the cross-band correlation provides a powerful tool for testing the redshift distributions and spectral energy densities of source population models, a more detailed analysis is required before making any strong conclusions. This will be the subject of a future study using the SANEPIC maps and including the BGS-Deep data.

We have measured the variance projected onto two dimensions, but galaxies are of course distributed in three dimensions. Knowledge of the redshift distribution of the galaxies in the CIB allows interpretation of these power spectra in terms of a bias factor quantifying the comparison to cold dark matter (CDM) spectra, which we find in the next section. As the angular scales we probe get smaller, non-linear growth eventually sets in, which we examine by fitting halo models to these spectra later in the paper.

To interpret our detection of correlations requires comparison to an underlying model whose parameters the data constrain. Such a model could contain details of the complete source population, including number counts, i.e., intensity distributions and redshift distributions, as well as a framework for describing the linear and non-linear clustering regimes. This latter part might involve, for example, a HOD that accounts for the galaxy distribution in a given dark matter halo as a function of luminosity (i.e., the so-called conditional luminosity function; Cooray 2006). For a background of primarily unresolved sources, the conditional luminosity function model is an improvement on the simple HOD; however, its increase in complexity comes with an increase of free parameters.

Given the very good fit to a constant $\Delta ^2_{k_{\theta }}$, we first explore the physical meaning of the BLAST power spectra in the context of a purely linear model. We assume that the galaxies which comprise the CIB have a power spectrum which is scaled from the power spectrum of dark matter at every redshift, P(k) = b²P_DM(k). Because P_DM is redshift dependent, estimating b² requires knowledge of the redshift distribution of the galaxies, which comprise the BLAST signals. To find this distribution, we adopt the model of Lagache et al. (2004), described in Section 6. It is worth noting that since we are really measuring an emission-weighted bias at rest-frame far-IR wavelengths, there will be some degeneracy between the value of b and changes in the redshift distribution or far-IR spectral shapes assumed for the sources.

We extend the fit into the non-linear regime to study whether the BLAST data can constrain parameters of the HOD. We can check whether our correlations are consistent with the model (e.g., by judging whether the same bias fits all three BLAST wave bands, and that the cross-band correlations of simulations made with the model agree with the data), but we will not explore how the model might be improved. It is important to understand that if the distribution of source redshifts or counts were different, then we would infer a different bias level. However, since the model we have adopted agrees with a large body of multi-wavelength observations, we are confident that these prescriptions are a reasonable first attempt.

Knowledge of the full redshift distribution of BLAST sources would allow us to estimate the redshift distribution of the measured clustering signal, i.e., the redshift range probed by the power spectrum of correlations due to clustering described by the P × dS/dz distribution. Devlin et al. (2009), and in more detail Patanchon et al. (2009), find the number counts of the BLAST sources, and Pascale et al. (2009) divides the redshift distribution roughly into four redshift bins. Precisely understanding the redshift distribution of those sources is a work in progress (e.g., E. L. Chapin et al. 2010, in preparation). In the meantime, we adopt the model of Lagache et al. (2004)—a model which approximates the counts and redshift distributions of the populations expected to make up the CIB—to describe the underlying source population.¹⁹ Specifically, it is a phenomenological model which extrapolates the local 60 μm luminosity function of IRAS sources (divided into "regular" and "star-forming" components) to longer wavelengths, assuming a set of spectral energy distribution templates. The extrapolation to higher redshift is parameterized as a mixture of luminosity and density evolution, and is constrained to be consistent with most of the available data on source counts, redshift distributions, and the far-infrared background intensity. It becomes less reliable with greater extrapolation to longer wavelengths, and therefore should be considered an approximation; however, it does reproduce the BLAST (Patanchon et al. 2009) counts at a level which is sufficient for our purposes. Furthermore, since the clustering signal is dominated by the faint source population, it is critical that the model is consistent with the level of the background at the BLAST wavelengths, and that the cross-band correlations of simulated maps agree with those measured, which we find to be the case.

The redshift distribution of the source population model at each wavelength is shown in Figure 4. Most of the clustering signal is coming from relatively high redshifts. The medians of the redshift distributions in Figure 5 are z = 1.61, 1.88, and 2.42, at 250, 350, and 500 μm, respectively, and the upper and lower quartiles of the distributions are z=(1.3, 2.2), (1.5, 2.7), and (1.7, 3.2), respectively. At the representative scale of 0.1 arcmin⁻¹ (red solid line in Figure 5), 95% of the background originates from sources 1.1 < z, 1.2 < z, and 1.4 < z in the three bands.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We first carry out a simple fit to a scaled version of the linear theory power spectrum of the dark matter P_DM. As can be seen in Figure 6, a simple biasing prescription provides a good fit to the BLAST data. The required bias levels are b = 3.8 ± 0.6, 3.9 ± 0.6, and 4.4 ± 0.7 at 250, 350, and 500 μm, respectively, with a reduced χ²_min ∼ 0.4 (with 10 degrees of freedom) in all three bands.

Zoom In Zoom Out Reset image size

More detailed modeling could be attempted. In principle, the cross-band measurements and wavelength dependence of the measured correlation amplitudes could be used to estimate the variation of bias with redshift, as discussed by Knox et al. (2001). However, we leave this to a future study.

As illustrated in Figure 6, a simple biasing prescription provides a good fit to our data. The main drawback of this simple fit is that it may not realistically account for the 1-halo, non-linear clustering component. Our data are expected to bracket the physical scales corresponding to the transition from linear to non-linear clustering regimes, and though these two contributions appear to have combined to look very much like a scaled linear bias, to conclude that there is no contribution from a 1-halo term may be unphysical. Therefore, in this section we use a particular implementation of the "halo model," which assigns galaxies to halos as a function of mass, and consequently probes into the territory of non-linear fluctuations, to explore how the 1-halo term could contribute power on small scales. This also allows us to discuss our results in the context of other measurements of galaxy clustering.

The halo model of large-scale structure has proven to be a powerful tool for describing the clustering properties of cosmic objects (for a review, see Cooray & Sheth 2002). Its main ingredient is the parameterization of the HOD (Peacock & Smith 2000), which describes how galaxies populate dark matter halos as a function of halo mass. The power spectrum of galaxies is written as the sum of two components: the 1-halo term, P_1h, which describes pairs of objects within the same dark matter halo, and the 2-halo term, P_2h, which accounts for pairs of objects in different halos, resulting in P(k, z) = P_1h(k, z) + P_2h(k, z). The number of pairs of galaxies within an individual halo is related to the variance of the HOD, σ²(M, z) = 〈N_gal(N_gal − 1)〉, while the number of pairs of galaxies in separate halos is simply the square of the mean halo occupation number, N(M, z) = 〈N_gal〉, (HON, hereafter). We model the HON using a central-satellite formalism (see, e.g., Zheng et al. 2005): this assumes that the first galaxy to be hosted by a halo lies at its center, while any remaining galaxies are classified as satellites and are distributed in proportion to the halo mass profile. Different HODs for central and satellite galaxies are then applied. For central galaxies, the mean HON, N_cen, is described by a step function such that halos above a minimum mass threshold M_min contain a single central galaxy and halos below this threshold contain no galaxies. For satellite galaxies, a power law in mass describes their mean HON (e.g., Zehavi et al. 2005)

$$\begin{equation} N_{\rm sat}(M) = \left(\frac{M}{M_1} \right) ^{\alpha }, \end{equation} \tag{ 10 }$$

where M₁ is the mass-scale at which a halo hosts exactly one satellite galaxy (in addition to the central galaxy). Both semi-analytic models (e.g., Berlind et al. 2000) and hydrodynamic simulations (e.g., Berlind et al. 2003) show that the distribution of galaxies within a halo is close to Poisson in the high-occupancy regime, i.e., when N_sat ≫ 1, and (strongly) sub-Poissonian in the low-occupancy regime. In order to agree with these results, satellite galaxies are assumed to be Poisson distributed at fixed halo mass. The distinction between central and satellite galaxies then automatically accounts for the sub-Poissonian behavior of the HOD in the low-occupancy regime (Zheng et al. 2005).

The 1- and 2-halo power spectra are

$$\begin{eqnarray} P_{\rm 1h}(k,z) & = & \int _{\cal {M}}n_{\rm halo}(M,z)\big[2N_{\rm cen}(M)N_{\rm sat}(M) u_{\rm DM}(k,z|M) \nonumber \\ & &+ N^2 _{\rm sat}(M) u^2 _{\rm DM}(k,z|M) \big] dM\big/n_{\rm gal}^{2}(z),\nonumber \\ P_{\rm 2h}(k,z) & = & P_{\rm DM}(k,z) \bigg [ \int _{\cal {M}}n_{\rm halo}(M,z)N_{\rm gal}(M,z) \nonumber \\ & &\times b(M,z)u_{\rm DM}(k,z|M)dM \bigg ]^{2}\bigg/n_{\rm gal}^2(z). \end{eqnarray} \tag{ 11 }$$

The meaning of the symbols here is as follows: P_DM is the linear power spectrum of dark matter, derived using the recipes of Eisenstein & Hu (1998) for the matter transfer function; n_halo is the halo-mass function (see Sheth et al. 2001); b is the linear bias parameter; u_DM is the normalized dark matter halo density profile in Fourier space; and n_gal is the mean number of galaxies per unit comoving volume at redshift z,

$$\begin{equation} n_{\rm gal}(z) = \int _{\cal {M}}n_{\rm halo}(M,z)\Bigg {[} 1 + \left(\frac{M}{M_1} \right) ^{\alpha }\Bigg {]}dM. \normalsize \end{equation} \tag{ 12 }$$

The expression for the 1-halo term implicitly assumes that the distribution of galaxies traces that of the dark matter, for which we have adopted the profile of Navarro et al. (1997,; hereafter NFW), with the same concentration parameter as Bullock et al. (2001). Since the NFW profile formally extends to infinity, it is necessary to artificially truncate the distribution at some radius, r_cut. Typically, this is chosen to be the virial radius of the halo; however, this may not necessarily be realistic. We address this by first adopting the assumption that r_cut = r_vir, and then exploring the consequences of relaxing that requirement, so that galaxies are allowed to lie further out. On large scales, where clustering is predominantly linear, u_DM ∼ 1, so that the 2-halo power spectrum simplifies to P_2h = b²_eff(z)P_DM(k, z), where b_eff(z) is the effective large-scale bias,

$$\begin{equation} b_{\rm eff}(z) = \int _{\cal {M}}n_{\rm halo}(M,z)N_{\rm gal}(M)b(M,z)\ dM/n_{\rm gal}(z). \end{equation} \tag{ 13 }$$

Our model has two free parameters, M_min and α, which we vary through 0 ⩽ α ⩽ 2 and 10 ⩽ log(M_min/M_☉) ⩽ 16, with steps of 0.05 in both log M and α. The best-fit values of the parameters are determined through a χ² minimization technique by fitting the observed power spectrum at each of the three BLAST bands simultaneously.

Throughout we assume that both M_min and α remain constant in time, although, in principle, they are functions of redshift. Whether these parameters evolve with redshift would be difficult to constrain from our data alone; nevertheless, our assumption is consistent with what is observed for other classes of high-redshift sources (e.g., quasars, see Porciani et al. 2004). For each M_min–α pair, the mass-scale M₁ is fixed by requiring that at every redshift, z, the number density of the background sources derived from the halo model formalism matches that predicted by the adopted source population model, i.e.,

$$\begin{equation} \int _{0}^{\infty }\frac{dN}{dS\ dz}(S,z)\ dS = n_{\rm gal}(z)\ dV_{\rm c}(z). \end{equation} \tag{ 14 }$$

Under the assumption that the dark matter halos are described by an NFW profile truncated at the virial radius, i.e., where r_cut = r_vir, and that galaxies within the halo trace the underlying dark matter distribution, the best fit to the observed angular power spectrum gives log(M_min/M_☉) = 11.5^+0.4_−0.1 and α ⩽ 1.0, with χ²_min = 16.3 with 10 degrees of freedom; i.e., the model is marginally consistent at the 2σ level.

While this model is not formally ruled out by the data, as shown by a dotted line in Figure 6, it poorly reproduces the shape of the observed power spectrum, which exhibits a steeper slope. If we were to weight this model to fit the large-scale power spectrum (e.g., $k \lesssim 0.08\; \rm arcmin^{-1}$), then it would over-predict the power on small scales. Arguing that perhaps the model describes the small scale, 1-halo term correctly, but is under-predicting the large scales, is not a good description because: (1) the 2-halo term is less sensitive to the underlying assumptions than the 1-halo term (for a discussion see Tinker et al. 2009); and (2) the observed small-scale power spectrum is still too steep to be accounted for by the shallower 1-halo term. Another possibility is that the redshift distribution of the background cumulative flux predicted by the adopted source count model is incorrect. In order to reproduce the observed shape of the angular power spectrum, the bulk of the background would have to originate from sources at z < 1, which is ruled out by Devlin et al. (2009); Marsden et al. (2009); and Pascale et al. (2009). However, a full investigation of the leeway in changing the redshift distribution (and degeneracies with other changes to the model) are beyond the scope of the present study.

In light of this, we explore the possibility that the discrepancy between the predicted power spectrum and the observed one is related to the modeling of the 1-halo term. There are two obvious ways to modify the shape and normalization of the 1-halo power spectrum. One is to allow the dark matter halos, although still following an NFW profile, to be truncated at a scale r_cut>r_vir. Thus, satellite galaxies are distributed over a larger volume. This idea is not new; Magliocchetti & Porciani (2003) find that in order to adequately fit the 1- and 2-halo term to the 2dF Galaxy Redshift Survey data set, it is necessary that the galaxies are allowed to reside out to twice the virial radius. Furthermore, from semi-analytic models Diaferio et al. (1999) show that blue (and hence star-forming) galaxies tend to reside in the outskirts of their host halos, while red galaxies are found closer to the halo center.

The second possibility is that the distribution of galaxies within the halos does not follow that of the underlying dark matter. For example, a power-law distribution ρ(r) ∝ r^−γ with γ < 2 would make the 1-halo angular power spectrum steeper than that predicted by an NFW profile. Similar arguments have been made by Watson et al. (2009), who found that by allowing the inner slope of the density profile to vary, they fit the small-scale clustering of luminous red galaxies quite well.

Here we only explore the first possibility, examining r_cut = 1, 2, 3, and 4 × r_vir. The results are shown in Figure 6, and summarized in Table 3. The best-fit model angular power spectrum for the case r_cut = 3 × r_vir is shown in Figure 7, while the corresponding HON and large-scale effective bias are shown in Figure 8. Note that while the value of the reduced χ²_min approaches unity for increasing values of r_cut, the best-fit values of M_min are negligibly affected by the changes, while α only marginally increases with increasing r_cut. The effective mass of the halo, M_eff, is the weighted mean over the halo-mass distribution:

$$\begin{equation} M_{\rm eff}(z) = \int _{\cal {M}}n_{\rm halo}(M,z)N(M,z)MdM/n_{\rm gal}(z). \normalsize \end{equation} \tag{ 15 }$$

The results for the large-scale effective bias b, and mass log(M/M_☉), at the respective medians of the redshift distributions of the sources contributing to the background in each band (see Figure 5) are 2.2 ± 0.2, 2.4 ± 0.2, and 2.6 ± 0.2; and 12.9 ± 0.3, 12.8 ± 0.2, and 12.7 ± 0.2, at 250, 350 and 500 μm, respectively. They are minimally affected by the change in r_cut, changing by less than 8% for the bias, and 3% for the mass, over the full range of r_cut explored.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows the contribution to the total clustering power spectrum from sources in different redshift slices, within the assumed source population model. As expected, the power is dominated by the contribution from sources in the range 0.7 < z < 1.5, with an increasing contribution from sources at z = 1.5–3.0 for increasing wavelengths, as is expected from Figure 5. This is consistent with the findings of Devlin et al. (2009); Marsden et al. (2009); and Pascale et al. (2009), who through stacking show that of the sources making up the CIB, the fraction at z>1.2 increases from 40% at 250 μm, to 50% and 60% at 350, and 500 μm.

Zoom In Zoom Out Reset image size

Finally, we can use the model to interpret the clustering in terms of a three-dimensional spatial correlation length, r₀. To do that, we Fourier transform the best-fit power spectrum to find the spatial correlation function, ξ(r), from which r₀ is then simply the linear comoving scale at which the correlation function equals 1 at each redshift. This is a model-dependent approach to estimating r₀, and as such should be considered an approximation. The more typical approach, which involves finding the angular correlation length, the redshift distribution, and deprojecting the signal by inverting the Limber equation, would result in very large uncertainties. The model-dependent r₀ is only mildly sensitive to the choice of r_cut, varying by 10% over the full range. The model-dependent values for r₀ are illustrated in Figure 10 as a solid line with a shaded area representing 3σ uncertainties. The three BLAST points are plotted at the flux-weighted median redshifts of the unique redshift distributions probed by the three bands (i.e., the distributions shown in Figure 4), and should not be interpreted as the locations of all the sources contributing power in those bands. At these effective redshifts, (i.e., z = 1.60, 1.86, and 2.15), we find r₀ = 4.9, 5.0, and $5.2 \pm 0.6\, h^{-1}\, \rm Mpc$ at 250, 350, and 500 μm, respectively.

Zoom In Zoom Out Reset image size

10.1. Comparison with other Observations

10.2. Clustering in Context

We report the detection of correlations in excess of Poisson noise in the CIB, over scales of approximately 5'–25', with BLAST at 250, 350, and 500 μm, at a level with respect to the CIB of ΔI/I = 15.1% ± 1.7%.

The CIB is made almost entirely out of individual sources distributed over a wide range of redshifts. We find that within the context of a reasonable model for the source population, the signal originates from galaxies in the redshift ranges 1.3 ⩽ z ⩽ 2.2, 1.5 ⩽ z ⩽ 2.7, and1.7 ⩽ z ⩽ 3.2, with median redshifts z = 1.61, 1.88, and 2.42, at 250, 350, and 500 μm, respectively. Fitting to the linear theory power spectrum, we find that the BLAST galaxies responsible for the CIB fluctuations have a bias parameter, b = 3.8 ± 0.6, 3.9 ± 0.6, and 4.4 ± 0.7.

We further interpret our results in terms of the halo model. We find that the simplest prescription does not fit very well. One way to improve the fit is to increase the radius at which we artificially truncate dark matter halos to well outside the virial radius. This may imply that the star-forming galaxies that we are seeing at z ∼ 1 are preferentially found in the outskirts of groups and clusters. This is consistent with related phenomena that have been observed at other wavelengths (Magliocchetti et al. 2004; Marcillac et al. 2007), as well as in simulations (Diaferio et al. 1999).

For a HOD with "satellite" galaxies occupying halos out as far as r_cut = 3r_vir, we find parameters log(M_min/M_☉ = 11.5^0.4_−0.1, and α = 1.1^+0.8_−0.1, resulting in effective biases b_eff = 2.2 ± 0.2, 2.4 ± 0.2, and 2.6 ± 0.2, and effective masses $\mathrm{log}(M_{\mathrm{\rm eff}}/M_{\odot }) = 12.9\pm 0.3$, 12.8 ± 0.2, and 12.7 ± 0.2 at 250, 350, and 500 μm, corresponding to spatial correlation lengths of r₀ = 4.9, 5.0, and $5.2 \pm 0.7 h^{-1} \rm Mpc$, respectively.

In the context of the model, we see that star formation is highly biased at z ≳ 1, unlike in the local universe, where analogous galaxy populations, such as IRAS galaxies, are found to be mildly anti-biased (e.g., M_eff ≲ 10¹¹ M_☉ and b = 0.86; Saunders et al. 1992).

We find relatively small values for θ₀, further confirming that the sources which make up the CIB are distributed over a wide range of redshifts, which has implications for planning future submillimeter clustering measurements. For example, as Knox et al. (2001) have argued, to match the precision of the bias measurement made from correlations in the background by using discrete sources will be very challenging. To achieve ∼5% accuracy would effectively require thousands of sources with exact redshifts, over tens of square degrees. When redshifts are only approximately known, this increases to hundreds of thousands of sources over hundreds of square degrees, which is only possible with instruments whose beams are smaller than 5''. Thus, while measuring the clustering from resolved sources has numerous advantages—for example, studying the clustering properties of a subset of galaxies—accurately measuring the large-scale bias is best achieved through correlation analysis of the background fluctuations. This will be a focus of future studies with BLAST, as well as with the Herschel and Planck satellites, and SCUBA-2.

We acknowledge the support of NASA through grant numbers NAG5-12785, NAG5-13301, and NNGO-6GI11G, the NSF Office of Polar Programs, the Canadian Space Agency, the Natural Sciences and Engineering Research Council (NSERC) of Canada, and the UK Science and Technology Facilities Council (STFC). C.B.N. acknowledges support from the Canadian Institute for Advanced Research. We thank the anonymous referee for her/his excellent comments, which improved this paper greatly, and Guilaine Lagache, Ravi Sheth, and Adam Muzzin for useful discussions and helpful advice. M.P.V. extends his warm thanks to Olivier Doré for his help and infinite patience, and Silvia Bonoli for kindly providing simulated clustering lengths of dark matter halos.