THE ATACAMA COSMOLOGY TELESCOPE: A MEASUREMENT OF THE COSMIC MICROWAVE BACKGROUND POWER SPECTRUM AT 148 AND 218 GHz FROM THE 2008 SOUTHERN SURVEY

Accurate measurements of the arcminute scale temperature anisotropies in the millimeter-wave sky are uncovering a complex, yet revealing picture (Lueker et al. 2010; Fowler et al. 2010, hereafter F10). On intermediate scales (500 ≲ ℓ ≲ 3000), the primordial acoustic features imprinted on the cosmic microwave background (CMB) at last scattering (z ≃ 1100) dominate, with subtle distortions expected from gravitational lensing by intervening large-scale structure. This intermediate range of multipoles is often called the damping tail of the CMB, as the acoustic oscillations are exponentially damped due to photon diffusion (Silk 1968; Bond & Efstathiou 1984). On the smallest scales (ℓ ≳ 3000), emission from radio- and dust-enshrouded star-forming galaxies (SFGs), together with the Sunyaev–Zel'dovich (SZ) effect (Zeldovich & Sunyaev 1969), dominates over primary CMB fluctuations.

The damping tail measurements are an additional test of the predictions of the ΛCDM cosmological model, a model that is an excellent fit to current CMB data (e.g., Larson et al. 2011; Reichardt et al. 2009b; Brown et al. 2009), large-scale structure measurements (e.g., Reid et al. 2010; Percival et al. 2010), supernova observations (e.g., Kessler et al. 2009; Amanullah et al. 2010), and a host of other astronomical observations (see, e.g., Spergel et al. 2007; Komatsu et al. 2011 for reviews). The amplitude of the fluctuations in the damping tail is a sensitive probe of matter fluctuations at k ≃ 0.1–0.25 Mpc⁻¹—thus, precision measurements constrain the spectral index of the primordial curvature perturbations, n_s, and its variation with scale. Because the positions of the high-order acoustic peaks are sensitive to the evolution of the sound speed of the universe and its composition, these measurements constrain the primordial helium fraction and the number of relativistic species including neutrinos (see, e.g., White 2001; Bashinsky & Seljak 2004; Trotta & Hansen 2004; Ichikawa et al. 2008a, 2008b; Komatsu et al. 2011). While the damping tail gives us leverage on early universe cosmology, the composite tail (ℓ ≳ 3000) is sensitive to a variety of astrophysical phenomena, including the properties of the intracluster medium, the redshift distribution and clustering of dusty submillimeter galaxies (Hall et al. 2010), and the physics of reionization (Huffenberger & Seljak 2005; Sehgal et al. 2007). Sehgal et al. (2010b) describe theoretical expectations for these small-scale measurements.

In this paper, we present a measurement of the CMB power spectrum over the range of multipoles 500 < ℓ< 10,000 from the Atacama Cosmology Telescope (ACT) using the 148 GHz and 218 GHz channel data from the 2008 observing season. ACT is a millimeter-wave, arcminute-resolution telescope (Fowler et al. 2007; Swetz et al. 2010) custom built to make precise observations of the microwave sky over the damping and composite tail regimes. In recent years, several groups have reported rapidly improving measurements of fluctuations over various portions of this multipole range: Bolocam (Sayers et al. 2009), QUaD (Brown et al. 2009; Friedman et al. 2009), APEX-SZ (Reichardt et al. 2009a), ACBAR (Reichardt et al. 2009b), Sunyaev Zel'dovich Array (SZA; Sharp et al. 2010), BIMA (Dawson et al. 2006), Cosmic Background Imager (CBI; Sievers et al. 2009), South Pole Telescope (SPT; Lueker et al. 2010; Hall et al. 2010), and ACT (F10). This paper enhances the results of F10 in a few ways. First, we augment our 148 GHz data with results from our 218 GHz channel. Second, the power spectrum estimation methods have been revised to provide enough angular frequency resolution to detect the acoustic features on the damping tail. Third, the area used for the power spectrum analysis has been increased from ≃220 deg² to ≃300 deg².

This paper is one in a set of papers describing ACT and its 2008 Southern survey. The ACT instrument is described in Swetz et al. (2010), a 148 GHz point-source catalog is presented in Marriage et al. (2010b, hereafter M10), and a 148 GHz SZ cluster catalog is presented in Marriage et al. (2010a). Menanteau et al. (2010) discuss the multi-wavelength follow-up of ACT clusters, while Sehgal et al. (2010a) present σ₈ constraints from SZ cluster detections. On the power spectrum side, Hajian et al. (2010) report on the calibration of ACT maps using cross-correlations with Wilkinson Microwave Anisotropy Probe (WMAP) seven-year maps, and Dunkley et al. (2010) present the constraints on cosmological parameters derived from the power spectrum presented here.

The paper is organized as follows. In Section 2, we briefly review the ACT instrument and observations made so far, touching upon mapmaking and beam determination techniques, and the calibration of our data. In Section 3, we describe the method used in power spectrum estimation. Simulations used to test our pipeline are discussed in Section 4. We present our results in Section 5 and discuss tests performed to validate our results. The main sources of foreground contamination and the methods used to treat them are discussed in Section 6. We discuss the lensing contribution to the power spectrum in Section 7, before concluding in Section 8.

ACT is a 6 m off-axis Gregorian telescope (Fowler et al. 2007) situated at an elevation of 5190 m on Cerro Toco in the Atacama desert in northern Chile. ACT has three frequency bands centered at 148 GHz (2.0 mm), 218 GHz (1.4 mm), and 277 GHz (1.1 mm) with angular resolutions of roughly 14, 10, and 09, respectively. The high-altitude site in the arid desert is excellent for microwave observations due to low precipitable water vapor and stability of the atmosphere. The tropical location of ACT permits observations on both the northern and southern celestial hemispheres. Further details on the instrument are presented in Swetz et al. (2010), F10, and references therein.³² As of this writing, ACT has completed three observing seasons (2007, 2008, and 2009) surveying two stripes on the sky: a 9°-wide stripe centered on declination −535 and a 5°-wide stripe centered on the celestial equator. The power spectrum presented here is derived from the southern stripe data from the two lower frequency channels in 2008. An important aspect of ACT's scanning strategy is the cross-linking of observations. Every point in the survey area is scanned along two different directions during the night. We use constant elevation scans centered on two different azimuths—once when the survey area is rising on the sky and again when it is setting toward the end of the observing night. In principle, this cross-linking allows an unbiased reconstruction of all modes in the map and permits recovery of the power spectrum for multipoles as low as a few hundred with errors completely dominated by cosmic variance. A more detailed account of the observing strategy of ACT can be found in F10 and M10.

2.1. Mapmaking

For the 148 GHz data, we use the maps from F10, in which details of the data reduction and mapping can be found. We model the data as d = Pm + n, where P projects the map m into the time stream and n is the noise, with covariance matrix N. There are no constraints on the contents of P and m, which can contain multiple components. For ACT, in addition to solving for the map of the sky, we solve simultaneously for noise correlated amongst the detectors. Under this model, the data can be described by

$$\begin{equation} \bm d = {\rm P}{\bm m}+ {\rm A}{\bm c} +\bm n, \end{equation} \tag{ 1 }$$

where A are (assumed constant) patterns of correlation across the array (such as a common mode) and c are the time streams associated with each pattern in A. This is mathematically equivalent to having separate blocks in a generalized projection matrix and a generalized map solution. The maximum-likelihood solution, ${\tilde{\bm m}}$, is then given by solving the standard mapping equation:

$$\begin{equation} { {\rm P}^{\rm T} {\rm N}^{-1}{\rm P} {\tilde{\bm {m}}}={\rm P}^{\rm T}{\rm N}^{-1} \bm d}. \end{equation} \tag{ 2 }$$

We solve for ${\tilde{\bm m}}$ iteratively using a preconditioned conjugate gradient (PCG) method (Press et al. 2007; Hinshaw et al. 2007).

The higher atmospheric noise in the 218 GHz data required some changes to the mapping pipeline. For the 148 GHz data, we found that taking the array patterns A to be the eigenvectors corresponding to the 10 largest eigenvalues of the data covariance matrix for each 15 minute chunk of time-ordered data (TOD) worked well for correlated noise rejection. For the 218 GHz data, substantial atmospheric power remains with this prescription, and so we altered it as follows: first, for each TOD we take the band-limited data between 0.25 and 4 Hz, find the eigenvalues and eigenvectors of the corresponding covariance matrix, and keep all modes with eigenvalues larger than 3.5² times the median eigenvalue (equivalently modes with time stream amplitudes larger than 3.5 times the median). We typically find 30–50 modes. Then, we create the covariance matrix from the data high-pass filtered above 4 Hz, project out the modes already found in the 0.25–4 Hz band, and keep all remaining modes with eigenvalue larger than 2.5² times the median. We typically find 1–2 additional modes in this step. Of the several different mode removal schemes we tried, we found this fairly aggressive one gave the best signal-to-noise ratio (S/N) on intermediate and small angular scales, where the 218 GHz data are most valuable, at the price of worse signal to noise and slower convergence of the mapper on large scales. Since the method estimates both the correlated modes and the map of the sky simultaneously, mode removal does not bias the maps, although it makes some sky map modes noisier.

The second change in the mapping is the use of a preconditioner with the conjugate gradient technique. In general, a preconditioner is an approximation to $(\rm {P^T N^{-1} P})^{-1}$ that is used to speed the convergence of the conjugate gradient solution. In general, we would use a two-piece preconditioner, with separate parts for the sky solution and the correlated noise component. However, if the columns of A are diagonal (which they are by construction), the natural preconditioner for the correlated part of the noise is simply the identity matrix and so we neglect it. The sky preconditioner is a diagonal matrix with elements equal to the inverse of the "n_obs" map—a map of the number of observations in each pixel. We find the use of this preconditioner greatly aids in the convergence of the conjugate gradient solution for 218 GHz data, unlike for the 148 GHz data where it makes little difference.

2.2. Beam Measurements

2.3. Calibration

Das et al. (2009) describe the basis of our power spectrum method. This paper improves on the analyses of Fowler et al. (2010) by optimizing the angular frequency resolution of the spectrum over the damping tail to resolve the acoustic features, as well as by extending the multipole range of the power spectrum to ℓ = 10, 000.

3.1. Fields Used for Spectra

The power spectrum is computed from independent analysis of four contiguous patches in the ACT southern stripe (see Figure 1) that together cover a rectangular area of 296deg² from α = 00^h22^m to 06^h52^m (55 to 103°) in right ascension and from δ = −55° to −50° in declination. Each patch is 5° × 148 in size. The full data set is also split into four subsets by distributing data from each of four successive nights into different subsets, thereby producing four independent maps of identical coverage and roughly equal depth. We refer to each subset as a split. Therefore, for each of the two frequency bands, we have four patches with four splits in each patch. Each split is properly cross-linked, i.e., it represents the maximum-likelihood solution of a data set containing both rising and setting observations. A typical white noise level in the 148 GHz season map (all four splits taken together) range from 25 to 50 μK–arcmin, while for the 218 GHz map it varies between 40 and 70 μK–arcmin, with the noise being higher toward the edges of the maps (Figure 1). The noise power spectra of the 148 GHz and 218 GHz maps are shown in Figure 2. At small scales, the noise is typically white, while at larger scales there is a plateau of noise with a red tail toward low multipoles, which is primarily attributable to atmospheric contamination.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Preprocessing of Maps

3.3. The Data Window

After the maps are prewhitened, four patches are cut from each split. Each split patch is then multiplied with a window. The window is a product of three components—a point-source mask, the n_obs map, and a tapering function. The details of the point-source mask are given in Section 6. To simplify the application, we generate a single n_obs map per patch for each frequency by summing over the four splits for that patch. For the cross-frequency spectra, we again generate a single n_obs map by summing the n_obs maps from the two frequencies. Multiplication by the n_obs map essentially downweights the relatively scantily observed and poorly cross-linked pixels toward the top and bottom edges of the map. Finally, to avoid ringing from the sharp edges of a patch, each patch is multiplied by a taper that gently rolls to zero at the edge of the map. We use a simple taper generated by taking a tophat function which is unity in the center and zero over 10 pixels at the edges, and convolving this with a Gaussian kernel of FWHM ≃25.

We will denote the two frequency channels as A and B and use i, j, k, l to denote the sub-season splits. The patches will be denoted by Greek indices. We will distinguish the data windows only by their patch indices and denote them by $W_{\alpha }(\hat{\mathbf {n}})$, suppressing the channel labels. Thus, split i of the patch α of frequency A, after windowing, becomes:³⁴

$$\begin{equation} \widetilde{T}^{i\!A}_{\alpha }(\hat{\mathbf {n}}) \equiv W_{\alpha }(\hat{\mathbf {n}}) T^{i\!A}_{F;\alpha }(\hat{\mathbf {n}}), \end{equation} \tag{ 3 }$$

where $T^{i\!A}_{F;\alpha }(\hat{\mathbf {n}})$ denotes the map after filtering and prewhitening and is related to the original map $T^{i\!A}_{\alpha }(\hat{\mathbf {n}})$ through the Fourier-space relation:

$$\begin{equation} T^{i\!A}_{F;\alpha }({\bm \ell }) = F_\ell T^{i\!A}_{\alpha }({\bm \ell }), \end{equation} \tag{ 4 }$$

where F_ℓ is a product of the prewhitening and beam transfer functions, the pixel window function and the high-pass sine-squared filter. Here, and in what follows, the boldface ℓ is meant to represent a two-dimensional (2D) wavevector in Fourier space, while we use regular ℓ to represent its absolute value. Symbols subscripted by ℓrepresent quantities that are isotropic, while ones shown with an argument (ℓ) or subscripted by ℓ stand for explicitly 2D quantities in Fourier space.

3.4. Azimuthal Fourier-mode Weights and Binning

After windowing the data, a 2D pseudo-spectrum is computed as

$$\begin{equation} \widetilde{C}_{\bm \ell }^{i\!A\times j\!B} = \mathrm{Re}\big[\widetilde{T}^{*i\!A}_{{\bm \ell }} \widetilde{T}^{j\!B}_{{\bm \ell }} \big], \end{equation} \tag{ 5 }$$

where the patch index has been suppressed to simplify notation. At this stage, we crop out and retain a rectangular area of the 2D spectrum defined by −10, 000 < ℓ_x < 10, 000 and −10, 000 < ℓ_y < 10, 000, throwing out all information at ℓ_x, ℓ_y > 10, 000. This trimming reduces the number n_F of Fourier-space pixels by a factor of ∼4 and is crucial for several subsequent operations, especially the direct computation of the mode-coupling matrix (see Section 3.5), which involves a computational step that scales as n²_F. The one-dimensional (1D) binned spectrum $\widetilde{C}_b$ results from averaging the 2D spectrum in annular bins:

$$\begin{equation} \widetilde{C}_b^{i\!A\times j\!B} = \sum _{{\bm \ell }} P_{b{\bm \ell }} \widetilde{C}_{\bm \ell }^{i\!A\times j\!B}, \end{equation} \tag{ 6 }$$

where P_bℓ is the binning matrix, whose value is zero when ℓ is outside the annulus defined by bin index b. Because the noise properties are not isotropic in Fourier space, we employ a 2D-weighted average, with the weights determined by the noise power spectrum. We briefly describe here how the weights are computed. For a given frequency, say, A, we compute the mean 2D auto-spectrum from the four split patches and subtract the mean 2D cross-power spectrum from it. This gives an estimate of the 2D noise power spectrum for the split patch N^AA_ℓ. We assume that each split has roughly equal noise and approximate the season noise power spectrum as N^AA_ℓ/4. We estimate the 2D variance of the spectral estimator C^iA×jB_ℓ as

$$\begin{equation} \sigma ^2\big(C_{\bm \ell }^{i\!A\times j\!B}\big) \simeq \big(C_{\bm \ell }^{\rm th} {b^A_\ell }^2 +N_{\bm \ell }^{A\!A}/4\big) \big(C_{\bm \ell }^{\rm th} {b^B_\ell }^2 +N_{\bm \ell }^{B\!B}/4\big), \end{equation} \tag{ 7 }$$

where C^th_ℓ is an azimuthally symmetric theoretical CMB spectrum and b_ℓ represents the beam transfer function. The first estimate of the azimuthal weights is chosen to be the inverse of the variance (Equation (7)):

$$\begin{equation} w^{i\!A\times j\!B}_{\bm \ell }= \frac{1}{ \sigma ^2\big(C_{\bm \ell }^{i\!A\times j\!B}\big) }. \end{equation} \tag{ 8 }$$

The resulting 2D Fourier-space map³⁵ of weights is typically noisy, and before using it to bin the power spectrum, we treat it in the following manner. In every annulus, defined by the bins b, we identify outlier pixels whose values are greater than a threshold α_b (typically chosen to be 10 times the median in that annulus) and set their values to α_b. We then divide the pixels in each annulus by their mean and smooth the resulting 2D Fourier-space map with a Gaussian with FWHM of three pixels. These operations are designed to produce well-behaved Fourier-space weight maps, with emphasis on bringing out the overall azimuthal variation in the weights by removing outliers and radial dependence. The Fourier-space weight map used in the 148 GHz × 148 GHz spectrum is shown in Figure 3. A comparison of this Fourier-space weight map with the signal and noise levels in Figure 2 shows the expected behavior of the weights being isotropic in the signal-dominated regime at low ℓ, with the anisotropy becoming more pronounced as the noise begins to dominate. The main sources of anisotropy in the weights are the rays of excess noise along and perpendicular to the scan directions, leading to two slightly rotated (because the rising and setting scans are not exactly orthogonal) "X-"shaped patterns of low weight regions. Another significant noise term is attributable to scan synchronous noise, likely caused by instabilities induced by acceleration at scan turnarounds. The noise leads to horizontal striping in the maps and manifests itself as an excess of power in a vertical strip −90 < ℓ < 90 in the Fourier space. We set the pixels inside this strip in the 2D Fourier-space weight map to zero. In terms of the 2D weights, the binning matrix can be expressed as

$$\begin{equation} P_{b{\bm \ell }} = \left. \frac{w_{\bm \ell }}{\sum _{\bm \ell }{w_{\bm \ell }}} \right|_{{\bm \ell }\in b}. \end{equation} \tag{ 9 }$$

It is noteworthy that both the n_obs weighting in real space and the azimuthal weighting in Fourier space lead to a 5%–10% lowering of the uncertainties in the final spectral estimates.

Zoom In Zoom Out Reset image size

3.5. Mode Coupling

From Equations (3) to (6), we can express the binned pseudo-spectrum $\widetilde{C}_b$ in terms of the underlying spectrum C_ℓ:

$$\begin{equation} \widetilde{C}_b^{i\!A\times j\!B} = \sum _{{\bm \ell },{\bm \ell }^{\prime } } P_{b{\bm \ell }} |W({\bm \ell }-{\bm \ell }^{\prime })|^2 F_{\ell ^{\prime }}^2 C^{i\!A\times j\!B}_{{\bm \ell }^{\prime }}. \end{equation} \tag{ 10 }$$

This quantity can be related to a binned version of the true spectrum C_b through an inverse binning operator Q_ℓb, which is unity when ℓ ∈ b and zero otherwise:

$$\begin{eqnarray} \nonumber \widetilde{C}_b^{i\!A\times j\!B} & =& \sum _{{\bm \ell },{\bm \ell }^{\prime },b^{\prime } } P_{b{\bm \ell }} |W({\bm \ell }-{\bm \ell }^{\prime })|^2 F_{\ell ^{\prime }}^2 Q_{{\bm \ell }^{\prime } b^{\prime }}C^{i\!A\times j\!B}_{b^{\prime }}\\ & \equiv & \sum _{b^{\prime }} M_{bb^{\prime }} C^{i\!A\times j\!B}_{b^{\prime }}, \end{eqnarray} \tag{ 11 }$$

where $M_{bb^{\prime }}$ is the mode-coupling matrix. We compute the $M_{bb^{\prime }}$ exactly without resorting to any of the 1D approximations often used in the flat-sky case (Hivon et al. 2002; Das et al. 2009; Lueker et al. 2010). The mode-coupling matrices are well behaved and stable to inversion. We define the unbiased, decorrelated estimator of the power spectrum (indicated with the circumflex):

$$\begin{equation} \widehat{C}_b^{i\!A\times j\!B} =\sum _{b^{\prime }} M^{-1}_{bb^{\prime }} \widetilde{C}_{b^{\prime }}^{i\!A\times j\!B}. \end{equation} \tag{ 12 }$$

For a single-frequency spectrum, we compute the cross-power spectrum in each patch as the mean of the six cross-power spectra (from four splits). We also compute the auto-spectrum as the mean of four auto-spectra from the splits. The variance for a given patch is determined by combining the noise estimate (obtained as one-fourth the difference of the mean auto and the mean cross-spectrum) in that patch with a theoretical signal power spectrum. We combine the four mean cross-power spectra from the patches with inverse variance weighting to obtain the final spectrum.

For the cross-frequency power spectra, we compute 12 cross-power spectra in each patch by crossing each split for a frequency with the three other splits from the other frequency. We do not cross the same splits (consisting of data from the same nights) between the two frequencies because of the possible contamination from covariant atmospheric noise among the two channels. We proceed by averaging those 12 per patch and then combining the cross-spectra from the patches with inverse weighting as above. In this case, we estimate the variance in a patch as the square root of the product of the single-frequency variances in that patch.

For parameter estimation (see Dunkley et al. 2010), we wish to compare our spectral estimates $\widehat{C}_b$ to theoretical power spectra C^th_ℓ. To do so, we find the bandpower window function, B_bℓ which converts C^th_ℓ to binned theoretical spectra C^th_b. (We suppress the superscripts on the Cs for clarity.) We begin with Equation (10) which relates the binned spectrum to the underlying 2D spectrum and introduce a matrix I_ℓℓ that describes the mapping of a theoretical spectrum defined at integer multipoles onto our discrete Fourier pixels. The binned theoretical spectrum follows from Equation (12):

$$\begin{eqnarray} \nonumber C_b^{\rm th} &=& \sum _{b^{\prime },{\bm \ell },^{\prime \prime }{\bm \ell }^{\prime }} M^{-1}_{bb^{\prime }} P_{b^{\prime }{\bm \ell }^{\prime \prime }} |W({\bm \ell }^{\prime \prime } -{\bm \ell }^{\prime })|^2 F_{\ell ^{\prime }}^2 I_{{\bm \ell }^{\prime } \ell } C_\ell ^{\rm th} \\ & \equiv & B_{b\ell } C_\ell ^{\rm th}. \end{eqnarray} \tag{ 13 }$$

Since each rectangular Fourier-space pixel spans a range of integral multipoles, it gets contributions from several multipoles in the theoretical spectrum. We choose I_ℓℓ as a Gaussian with its peak at the pixel center (ℓ_x, ℓ_y), and width corresponding to the pixel size, so as to give more weight toward the center of the pixel. The bandpower window functions are insensitive to the exact form of this response function. As with the mode-coupling matrix, we compute B_bℓ exactly.

3.6. Summary of Power Spectrum Method

Simulations play an important role both in verifying that our mapper is working correctly and in understanding our spectrum. Specifically, we use simulations to test for map convergence, estimate bandpower covariances, and confirm analytic estimates of the power spectrum uncertainties.

To test for convergence, we perform a set of end-to-end simulations of the 148 GHz and 218 GHz maps. We simulate signal-only time streams by mock observations of input celestial maps (from Sehgal et al. 2010b) and generate noise time streams which capture the main features of the noise in the data. To test for convergence, we simulate a signal-plus-noise map, subtract from it a noise-only map at the same iteration, and compare the power spectrum of the resulting difference map with that from a signal-only map. For 148 GHz, the maps are well converged to multipoles as low as ℓ ≳ 200 by iteration 1000, and we keep ℓ ⩾ 500 in our final spectrum. For the 218 GHz maps, where higher atmospheric noise warranted a slightly different mapmaking approach (see Section 2.1), convergence at large scales is slow. After 1000 iterations, the 218 GHz maps are converged down to ℓ = 1000, but not at smaller ℓ. Moreover, the 218 GHz spectrum is highly noise dominated at low multipoles. Thus, we limit spectra involving 218 GHz to ℓ ⩾ 1500 in our analyses.

In order to investigate the bandpower covariance and to validate the analytic prescription for uncertainties in the spectrum, we ran a large set of Monte Carlo simulations that contain key aspects of the noise properties manifested in the data. Due to the iterative nature of our map-making pipeline, it is prohibitively expensive to run a large number of end-to-end simulations starting from simulated time streams. Therefore, we adopt an intermediate approach, in which we find a prescription from the data for generating realistic realizations of noise in map space, to which we add signal realizations. We start with two quarter-season split maps, $T_0(\hat{\mathbf {n}})$ and $T_1(\hat{\mathbf {n}})$, with corresponding n_obs maps labeled $n_0(\hat{\mathbf {n}})$ and $n_1(\hat{\mathbf {n}})$. We use the difference of these maps to estimate a one-hit noise map $D_{01}(\hat{\mathbf {n}})$ (i.e., a noise map representative of the variance in each pixel had it been observed only once),

$$\begin{eqnarray} &&D_{01}(\hat{\mathbf {n}}) = \frac{T_0(\hat{\mathbf {n}})-T_1(\hat{\mathbf {n}})}{\sqrt{{n_{0}^{-1}(\hat{\mathbf {n}})}+{n_{1}^{-1}(\hat{\mathbf {n}})}}}. \end{eqnarray} \tag{ 14 }$$

With the other two splits, we compute another one-hit noise map $D_{23}(\hat{\mathbf {n}})$. We compute the average of the 2D power spectra of D₀₁ and D₂₃, and use its square root as the amplitude for a Gaussian random field, and generate random one-hit noise maps by randomizing the phases and converting back to real space. Finally, we divide one such map by $\ssty \sqrt{n_i(\hat{\mathbf {n}})}$ to generate a realization of the noise for split i. This procedure captures the non-white and anisotropic aspects of map noise, as shown in Figure 4.

Zoom In Zoom Out Reset image size

We also simulate a random residual point-source component in our Monte Carlo realizations. The point-source model for 148 GHz consists of two components—a radio source population and a dusty submillimeter population, both based on Toffolatti et al. (1998). Source maps are generated assuming a Poisson distribution. For each flux bin S, spanning the range 0.01–20 mJy, we generate the number of point sources in that bin as a Poisson realization from the model flux distribution dN/dS (the upper limit is chosen to be the approximate 5σ detection threshold at 148 GHz). Once the number of sources in each bin is generated, their individual fluxes are dithered to redistribute the fluxes within each bin, and the sources are laid down at random positions in the map. For 218 GHz, we scale the radio and IR components by their appropriate spectral indices, based on the findings of Dunkley et al. (2010). With these settings, we closely reproduce the level of the Poisson point-source signal seen in the high multipole regime of the data power spectra. We estimate the Poisson contribution to the error bars on bandpowers from the scatter of 500 simulations of point-source-only maps and add that as a correction to the (signal+noise)-only error budget obtained from the signal. In these point-source simulations, we have neglected the contribution from the correlated IR component since that term is important only at ℓ < 2500 where the primary CMB and noise are several orders of magnitude larger. At the current level of noise in our maps, any point-source correction, correlated or Poissonian, is largely negligible, and we ignore their contribution in the analytic error bar estimation and bandpower covariance.

We estimate 148 GHz and 218 GHz spectra for 960 signal+noise realizations. Each realization involves simulating four patches and four splits in each patch for each frequency, and estimating the spectra with exactly the same methods as used on the data. The resulting spectra are used for quantifying the covariance between bandpowers. Figure 5 displays the normalized covariance matrix, ${\rm Cov}(i,j)/\sqrt{{\rm Cov}(i,i) {\rm Cov}(j,j)}$, where i, j denote bin indices, showing that correlations between adjacent bins are insignificant at the 10% level.

Zoom In Zoom Out Reset image size

The uncertainties in the bandpowers are evaluated using the analytic prescription described in the Appendix. First, a large set of isotropic white noise simulations are generated to validate the analytic formula for the single-frequency and the cross-frequency error bars. These simulated error bars are found to be in agreement with analytic prescription to less than a percent. The next set of simulations are run with anisotropic noise, Poisson point-source realizations, and with real-space and Fourier-space weighting as described in Section 3, and small corrections to analytic prescription due to the anisotropic nature of the noise and weighting are evaluated against the Monte Carlo simulations.

Applying the methods described in Section 3, we compute the decorrelated bandpowers for the 148, 218, and 148 × 218 GHz spectra from the unbiased map solutions. The bandpowers are presented in Table 1 and are displayed in Figures 6 and 7. In Figure 6, we have plotted the lower multipole portion of the power spectrum 500 < ℓ < 3000 with an ℓ⁴ scaling to emphasize the higher-order acoustic peaks in the primary CMB spectrum. In the 148 GHz spectrum, we clearly resolve the second through the seventh peaks of the CMB. As discussed in Dunkley et al. (2010), the higher-order peaks provide leverage on cosmological parameters such as the spectral index and its running, the primordial helium fraction and the number of relativistic species. Figure 7, on the other hand, emphasizes the high multipole tails of the spectra, where the signal is dominated by emission from dusty SFGs and unresolved radio sources. The intermediate range 2000 < ℓ < 4000 gets a significant contribution from the thermal SZ effect and lets us constrain the amplitude of the SZ spectrum.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Single Frequency Bandpowers, ${\cal B}_b = \ell _b (\ell _b+1) C_b/2\pi$ (μK²)

ℓ Range	Central ℓb	148 GHz		148 GHz × 218 GHz		218 GHz
		${\cal B}_b$	$\sigma ({\cal B}_b)$	${\cal B}_b$	$\sigma ({\cal B}_b)$	${\cal B}_b$	$\sigma ({\cal B}_b)$
491–540	515	2423.0	270.8	...	...	...	...
541–590	565	2701.4	274.8	...	...	...	...
591–640	615	1952.1	205.4	...	...	...	...
641–690	665	1701.9	167.7	...	...	...	...
691–740	715	1895.7	179.5	...	...	...	...
741–790	765	2386.5	202.0	...	...	...	...
791–840	815	2415.5	200.5	...	...	...	...
841–890	865	2168.0	178.2	...	...	...	...
891–940	915	1567.5	138.3	...	...	...	...
941–990	965	1353.6	109.9	...	...	...	...
991–1040	1015	976.0	89.2	...	...	...	...
1041–1090	1065	1089.5	92.9	...	...	...	...
1091–1140	1115	1229.6	100.2	...	...	...	...
1141–1190	1165	1169.7	88.2	...	...	...	...
1191–1240	1215	861.9	73.4	...	...	...	...
1241–1290	1265	819.3	68.7	...	...	...	...
1291–1340	1315	740.8	60.6	...	...	...	...
1341–1390	1365	746.8	60.5	...	...	...	...
1391–1440	1415	782.4	64.4	...	...	...	...
1441–1490	1465	718.8	55.8	...	...	...	...
1491–1540	1515	619.3	50.1	559.6	43.6	585.9	87.6
1541–1590	1565	495.6	44.5	489.0	40.2	508.0	83.0
1591–1640	1615	467.4	38.8	428.3	35.3	480.1	77.9
1641–1690	1665	380.4	34.0	382.0	32.7	446.5	74.0
1691–1740	1715	351.6	33.6	347.8	32.1	436.8	75.0
1741–1790	1765	359.6	32.3	336.7	30.7	470.6	74.1
1791–1840	1815	334.4	29.4	336.9	28.3	396.1	65.3
1841–1890	1865	274.8	25.3	276.6	25.8	371.5	65.0
1891–1940	1915	226.3	23.3	198.8	23.7	199.7	62.1
1941–1990	1965	259.7	23.1	254.6	23.9	310.2	61.9
1991–2040	2015	235.8	21.6	226.0	22.5	340.2	60.6
2041–2090	2065	242.8	22.8	227.7	23.9	283.9	63.5
2091–2140	2115	168.8	18.6	206.9	21.2	268.3	57.8
2141–2190	2165	164.7	17.8	168.4	19.7	203.4	54.9
2191–2240	2215	152.2	17.2	151.0	19.7	269.7	60.5
2241–2290	2265	132.2	15.8	114.6	18.1	206.7	56.3
2291–2340	2315	114.9	14.7	136.5	18.1	135.8	54.2
2341–2390	2365	106.3	15.4	135.3	19.1	222.3	59.2
2391–2700	2550	82.9	5.4	99.7	6.7	157.9	21.2
2701–3000	2850	57.5	5.3	72.6	6.5	189.1	20.8
3001–3400	3200	37.2	4.8	54.0	5.5	158.3	16.9
3401–3800	3600	39.8	5.6	50.9	5.9	145.1	16.6
3801–4200	4000	48.1	6.8	62.9	6.7	204.5	17.9
4201–4600	4400	34.9	8.3	74.6	7.8	198.7	19.7
4601–5000	4800	36.4	10.2	98.0	9.3	248.7	22.2
5001–5900	5450	55.8	9.6	106.3	8.3	281.7	18.9
5901–6800	6350	56.6	15.4	131.7	12.5	378.3	26.8
6801–7700	7250	127.8	25.6	127.7	19.1	463.7	38.1
7701–8600	8150	115.3	42.4	228.0	30.0	620.6	55.6
8601–9500	9050	208.3	73.1	175.9	48.6	812.8	85.0
9501–9900	9750	316.9	172.4	152.4	108.1	695.0	176.1

Download table as:  ASCIITypeset image

The data are checked for consistency by performing various null tests, in which difference maps are constructed to cancel true sky signals, and their power spectra examined. We describe the suite of tests here. As described in Section 3.1, the TODs are split into four subsets, generating four maps with independent detector noise. We verify that these maps are consistent with each other by generating difference maps from each pair and computing the two-way cross-spectra from the three independent pairs of difference maps:

$$\begin{eqnarray} &&C_{\bm \ell }^{(i-j)\times (k-l)} = \langle \Delta T^{(i-j)*} ({\bm \ell }) \Delta T ^{(k-l)} ({\bm \ell })\rangle, \end{eqnarray} \tag{ 15 }$$

where ΔT^(i−j) ≡ Tⁱ − T^j and the indices i, j, k, l range from 0 to 3. The same point-source mask used to mask the full data set is applied to the difference maps before calculating the spectra; all other settings also remain the same. The difference maps, ΔT^(i−j), are downweighted by the same n_obs maps summed over the four splits used in the windowing of the CMB data (see Section 3.3). Similarly, when binning the power spectra, we use the same azimuthal weights in each patch described in Section 3.3. The three spectra are shown in the top two panels of Figure 8 for the 148 (top left) and 218 (top right) GHz maps. We compute the bandpowers in the range 500 < ℓ < 10, 000 (1500 < ℓ < 10, 000) for 148 (218) GHz. The mean of the three spectra is consistent with null with χ² = 43.5(19.5) for 148 (218) GHz, with 51(31) degrees of freedom.

Zoom In Zoom Out Reset image size

Another null test probes the consistency of data taken with the telescope accelerating as it reverses direction at the ends of the scan (turnarounds). Note that for the standard maps, no turnaround cut is applied. Four new split maps are made cutting data near the turnarounds, amounting to losing ≃24% of the total data. Two difference maps are made by pairing one split of the standard map with a different split of the new maps (we avoid differencing the same splits as they have very similar noise structure), and a two-way cross-power spectrum is produced. Any artifact due to the turnaround would be left in these difference maps and might produce excess power. However, we find the resulting spectra, shown in the second row of Figure 8, to be consistent with null with a χ² = 62 and 18, respectively, for the 148, 218 maps (again for 51 and 31 degrees of freedom).

The last two null tests probe possible systematics associated with the focal plane and from diurnal effects. The first compares maps made from only the top half of the detector array with those from the bottom half. The second compares maps processed with data only from the middle half of the night to the standard maps. The resulting spectra are shown in the final two panels of Figure 8 and are also consistent with null.

The main foregrounds in the 148 GHz and 218 GHz bands at the angular scales considered here are point sources. We mask the detected ones as described below. Another foreground is the diffuse Galactic dust, which is discussed below in some detail. In the companion paper, Dunkley et al. (2010), we also consider the contributions to the power spectra from the thermal and kinetic SZ effects and the clustering of infrared point sources.

The point sources include radio sources and dusty SFGs. The radio sources are typically active galactic nuclei exhibiting synchrotron-dominated spectra, with emission extending down to lower radio frequencies. The dusty SFGs are characterized by absorption of ultraviolet photons from star-forming regions by dust. The dust reemits into graybody radiation that rises with increasing frequency into submillimeter bands.

For the purpose of extracting the underlying CMB signal, it is important to identify and mask bright point sources in the maps before the power spectrum is computed. A detailed study of the point-source population detected in the ACT 148 GHz 2008 survey is presented in M10. Point sources are detected with the matched filter method (Tegmark & de Oliveira-Costa 1998; Wright et al. 2009; Vieira et al. 2010). To construct the point-source masks, we consider detections with S/N >5 in either frequency. Over the 296deg² area used for power spectrum estimation, we find a total of 164 point sources, taking the union of the detections at the two frequencies. Before computing the spectrum, we mask out a 10' diameter region around each of the 164 point sources. This amounts to masking about 1% of the total area. As discussed in Section 3.2, the prewhitening operation effectively guards against the leakage of power due to the application of this mask.

We estimate the Galactic dust contribution to our power spectra by cross-correlating ACT maps with the predictions for infrared cirrus emission at CMB frequencies based on the multicomponent dust model "8" of Finkbeiner et al. (1999, hereafter FDS). The FDS-based dust maps at 148 GHz and 219 GHz are available as a part of the Sehgal et al. (2010b) simulations.³⁶ We resample these maps into ACT pixelization and cut out the ACT patch regions. The 219 GHz FDS map is shown in Figure 9. Patches 0 through 3 labeled in the figure are situated off the Galactic plane with central Galactic latitudes of −64°, −56°, −43°, and −28°, respectively. The power spectra of the predicted dust signal in each patch are shown in Figure 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We proceed with the cross-correlation as follows. We express the ACT map as a sum of the cosmic and galactic dust components: T_ACT = T_CMB + ξT_FDS, where ξ is the dust amplitude, predicted to be unity by FDS. We take two approaches to estimate ξ: one based on real-space operations and the other based on power spectra.

In the real-space approach, we prewhiten both the ACT and FDS maps, convolve the ACT maps with a 61 FWHM Gaussian beam appropriate to the FDS resolution, apply the ACT point-source mask described above, filter out modes below ℓ = 500 (1500) for 150 GHz (220 GHz), and produce the maps ${\cal T}^{i}_{\rm {ACT}}(\hat{\mathbf {n}})$ and ${\cal T}^{i}_{\rm {FDS}}(\hat{\mathbf {n}})$ in patch i. Then we estimate ξ in patch i as

$$\begin{equation} \xi ^{i}_{\rm real} =\frac{ \int d^2 \hat{\mathbf {n}}{\cal T}^{i}_{\rm {ACT}}(\hat{\mathbf {n}}) {\cal T}^{i}_{\rm {FDS}} (\hat{\mathbf {n}}) }{\int d^2 \hat{\mathbf {n}}{\cal T}^{i}_{\rm {FDS}}(\hat{\mathbf {n}}) {\cal T}^{i}_{\rm {FDS}}(\hat{\mathbf {n}})}. \end{equation} \tag{ 16 }$$

The error on the above estimate is obtained by cross-correlating 960 random CMB plus noise simulations (see Section 4) with the dust maps. In the spectrum-based approach, we estimate ξ in each multipole bin b as the ratio of the ACT–FDS cross-spectrum to the auto-spectrum of the FDS map for the same patch:

$$\begin{equation} \xi _b^i = \frac{C_b^{{\rm ACT}\times {\rm FDS};i}}{C_b^{{\rm FDS} \times {\rm FDS};i}}. \end{equation} \tag{ 17 }$$

The ACT-FDS cross-spectrum is calculated using both pixel-space weighting and Fourier-space azimuthal weighting, as with the ACT power spectrum, with the FDS and ACT beams deconvolved. The final estimate, ξⁱ_spec, is obtained as the mean across the bins over the range 500 < b < 3000 (1500 < b < 3000) for 150 GHz (220 GHz) where the upper limit is dictated by the FDS resolution. The uncertainty on ξⁱ_spec is estimated as the scatter in the same quantity computed from the cross-spectra between the FDS maps and 960 random CMB plus noise realizations. As a check for systematics, we circularly shift the FDS patches by one patch, such that patch_0 of ACT is crossed with patch_3 of FDS, patch_1 of ACT with patch_0 of FDS, etc., and recompute the real- and Fourier-space estimates of ξ. The results are shown in Figure 11. The errors shown arise from treating the FDS template as deterministic.

Zoom In Zoom Out Reset image size

This figure shows that a dust-correlated signal is observed. More work will be needed to understand how well the FDS template represents faint dust well off the Galactic plane (see, e.g., Veneziani et al. 2010). The fact that such a small signal (ℓ²C_ℓ/2π ∼ 0.1–1.0 μK² at ℓ = 3000 and 148 GHz) can be recovered with our map-based technique demonstrates the power of our maximum-likelihood map estimation. Recently, Hall et al. (2010) have also reported a measurement of infrared cirrus in the SPT maps.

Subtracting the FDS template from our maps has negligible effect on the power spectra, consistent with the minuteness of the dust signal.

The CMB photons are deflected by large-scale structure potentials along their path from the last scattering surface at z ≃ 1100 to us (see Lewis & Challinor 2006 for a review). The typical (rms) deflection in the ΛCDM model is about 27 and the deflection pattern is coherent over degree scales, comparable to the size of the acoustic features on the primary CMB. These coherent deflections produce distortions of the hot and cold spots on the CMB, leading to a broadening of their size distribution. In the power spectrum, this effect manifests itself as the smoothing of the acoustic peaks, which can be used as a signal to look for lensing. The first attempt at a detection of lensing in the power spectrum was made in Reichardt et al. (2009b) using the data from the ACBAR experiment. They quantified the effect of including lensing in their analysis through the log ratio of the lensed to unlensed Bayesian evidence ($\Delta {\ln \cal E}$), and by combining WMAP5 and ACBAR data sets found $\Delta {\ln \cal E}=2.63$. Calabrese et al. (2008) analyzed the ACBAR data with a different approach, where they introduced a scaling of the power spectrum of the lensing potential (C^Ψ_ℓ → A_LC^Ψ_ℓ, with A_L = 0 corresponding to no lensing and A_L = 1 to the standard ΛCDM expectation). With WMAP5+ACBAR, they found A_L = 3.0^+0.9_−0.9 (68% confidence level). According to footnote 17 of Reichardt et al. (2009b), a reanalysis of the ACBAR data using this parameterization yields A_L = 1.60^+0.55_−0.26.

Our 148 GHz power spectrum is shown against lensed and unlensed models in Figure 12. We use the parameterization of Calabrese et al. (2008), and using the parameter estimation methodology described in the companion paper Dunkley et al. (2010), we constrain the lensing parameter A_L based on WMAP7 and ACT power spectra. Figure 13 shows the marginalized 1D likelihood for A_L using WMAP7+ACT. We find A_L = 1.3^+0.5(+1.2)_{−0.5(−1.0)} at 68%(95%) confidence, with the best-fit lensed CMB spectrum with A_L = 1 having an improved goodness-of-fit to the WMAP7+ACT data of Δχ² = 8 less than the unlensed model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We check for systematics that might have given rise to a spurious lensing signal. The projection scheme (cylindrical-equal-area) used for the ACT maps is not particularly optimized for lensing studies—so we test whether this projection could introduce a lensing-like signal. We simulate a low-noise unlensed CMB signal and run it through the mapmaking pipeline and try to reconstruct a "lensing convergence" in the resulting map using standard quadratic estimator techniques (Hu & Okamoto 2002). We find the reconstructed convergence power spectrum to be consistent with null (there is a small known bias at high multipoles that we have entirely traced to mode coupling due to the finiteness of the patch), showing that the projection does not introduce any significant lensing-like signal in our patches. To further test if any other step in our pipeline could produce spurious peak smearing in the power spectrum, we generate an end-to-end simulation of a noisy map exactly as described in Section 4, only this time replacing the lensed CMB signal time stream with its unlensed version. The resulting maps are then processed through the power spectrum pipeline, and the power spectrum is then analyzed with the parameter estimation method described in Dunkley et al. (2010). We find the lensing amplitude parameter, A_L, described above to be consistent with zero (see Figure 13).

Here we have used the smearing of the acoustic features to look for the lensing signal. More promising ways of extracting this signal involve using information beyond the two-point correlation (optimal quadratic estimators) and cross-correlation with tracers of large-scale structure. Early efforts cross-correlating the lensing reconstruction in WMAP three-year maps with luminous red galaxies, radio sources, and quasars (Smith et al. 2007; Hirata et al. 2008) yielded ∼3σ evidence for lensing. Data from WMAP are not the best suited for this kind of study because of its limited angular resolution. With arcminute-resolution CMB data, such as from ACT, there is a much higher potential for a detection. Such efforts are underway with the ACT maps.

We have presented a measurement of the power spectrum of the CMB sky observed with the ACT at 148 GHz and 218 GHz. The 148 GHz spectrum spans a large dynamic range from ℓ = 500 to ℓ = 10, 000, covering the damping tail, where the primary anisotropies with the higher-order acoustic peaks dominate, to the composite high multipole tail of the CMB where emission from dusty galaxies and radio sources and the SZ effect contribute. The second through the seventh acoustic peaks of the CMB are clearly visible in this spectrum. For the 148 × 218 GHz and 218 × 218 GHz spectra, we present measurements from 1500 < ℓ < 10, 000. The high multipole information from these spectra helps us constrain the nature of the point-source populations that contribute to that range, and in turn, constrain the SZ contribution to the 148 GHz spectrum. We measure the presence of a very faint dust signal at high Galactic latitudes. Recovering such a faint signal gives us confidence in the fidelity of the maps. We find evidence of gravitational lensing of the CMB in the power spectrum at the 2.8σ level. Constraints on cosmological parameters and the interpretation of high multipole spectra in terms of point-source populations and the SZ effect are presented in the companion paper by Dunkley et al. (2010).

ACT is on the Chajnantor Science preserve, which was made possible by the Chilean Comisión Nacional de Investigación Científica y Tecnológica. We are grateful for the assistance we received at various times from the ALMA, APEX, ASTE, CBI/QUIET, and NANTEN2 groups. The PWV data come from the public APEX weather Web site. Field operations were based at the Don Esteban facility run by Astro-Norte. Reed Plimpton and David Jacobson worked at the telescope during the 2008 season. We thank Norm Jarosik for support throughout the project. This work was supported by the US National Science Foundation through awards AST-0408698 for the ACT project, and PHY-0355328, AST-0707731, and PIRE-0507768. Funding was also provided by Princeton University and the University of Pennsylvania. The PIRE program made possible exchanges between Chile, South Africa, Spain, and the US that enabled this research program. Computations were performed on the GPC supercomputer at the SciNet HPC Consortium. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada, the Government of Ontario, Ontario Research Fund-Research Excellence, and the University of Toronto. We thank the referee for useful comments. S.D. acknowledges support from the Berkeley Center for Cosmological Physics Fellowship. S.D. thanks Christian Reichardt and Oliver Zahn for useful discussions. We thank Bruce Bassett for suggestions on testing lensing in the power spectrum. S.D., A.H., and T.M. were supported through NASA grant NNX08AH30G. J.D. acknowledges support from an RCUK Fellowship. R.H. received funding from the Rhodes Trust. A.D.H. received additional support from a Natural Science and Engineering Research Council of Canada (NSERC) PGS-D scholarship. A.K. and B.P. were partially supported through NSF AST-0546035 and AST-0606975, respectively, for work on ACT. L.I. acknowledges partial support from FONDAP Centro de Astrofísica. R.D. was supported by CONICYT, MECESUP, and Fundación Andes. E.S. acknowledges support by NSF Physics Frontier Center grant PHY-0114422 to the Kavli Institute of Cosmological Physics. K.M., M. Hilton, and R.W. received financial support from the South African National Research Foundation (NRF), the Meraka Institute via funding for the South African Centre for High Performance Computing (CHPC), and the South African Square Kilometer Array (SKA) Project. Y.-T.L. acknowledges support from the World Premier International Research Center Initiative, MEXT, Japan. N.S. is supported by the US Department of Energy contract to SLAC no. DE-AC3-76SF00515. We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA). Support for LAMBDA is provided by the NASA Office of Space Science. The data will be made public through LAMBDA (http://lambda.gsfc.nasa.gov/) and the ACT Web site (http://www.physics.princeton.edu/act/).

Here we derive an analytic expression for the expected error bars on the cross-frequency cross-power spectrum, assuming that the signal and the noise can be approximated as Gaussian random fields. We do not treat the Poisson component from point sources here, as the Poisson noise correction is negligible at our current noise level. As before, we denote the two frequency channels as A and B and use lowercase romans i, j, k, l, etc., to denote the sub-season data splits.

The estimator for the cross-frequency cross-power spectrum in bin b is given by

$$\begin{equation} \widehat{C}^{A\times B}_b = \frac{1}{n_d(n_d-1)} \sum _{i,j; i \ne j} \widehat{C}_b ^{(i\!A \times j\!B)}, \end{equation} \tag{ A1 }$$

where n_d is the number of sub-season data splits (four in our case), and

$$\begin{equation} \widehat{C}_b^{(i\!A \times j\!B)} = \frac{1}{\nu _b} \sum _{{\bm \ell }\in b} T^{*i\!A}_{{\bm \ell }} T^{j\!B}_{{\bm \ell }}, \end{equation} \tag{ A2 }$$

where the sum is over pixels in the annulus in Fourier space contained in the bin b and ν_b is the effective number of pixels in the annulus over which the sum is taken. Evidently,

$$\begin{eqnarray} \big\langle C_b ^{(i\!A \times j\!B)} \big\rangle &=& C_b, \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} \big\langle \widehat{C}^{A\times B}_b \big\rangle &=& C_b. \end{eqnarray} \tag{ A4 }$$

First, we will compute the expected covariance between a pair of estimators:

$$\begin{eqnarray*} \Theta ^{(i\!A\times j B);(k\!A\times l B)}_{bb} &=& \big\langle \big(\widehat{C}_b ^{(i\!A \times j\!B)}-\big\langle \widehat{C}_b ^{(i\!A \times j\!B)}\big\rangle \big) \nonumber \\ && \times \big(\widehat{C}_b ^{(k\!A \times l\!B)} - \big\langle \widehat{C}_b ^{(k\!A \times l\!B)}\big\rangle \big)\big\rangle, \end{eqnarray*}$$

where i ≠ j and k ≠ l. From the definition (A2), we can write

$$\begin{eqnarray} \Theta ^{(i\!A\times j B);(k\!A\times l B)}_{bb} &=& \frac{1}{\nu _b^2} \sum _{{\bm \ell }\in b} \sum _{{\bm \ell }^{\prime } \in b} \big[ \big\langle T^ {*i\!A}_{{\bm \ell }} T^{j\!B}_{{\bm \ell }} T^{*k\!A}_{{\bm \ell }^{\prime }} T^{l\!B}_{{\bm \ell }^{\prime }} \big\rangle \big] \nonumber \\ && - \big\langle C_b ^{(i\!A \times j\!B)} \big\rangle \big\langle C_b ^{(k\!A \times l\!B)} \big\rangle. \end{eqnarray} \tag{ A5 }$$

Applying Wick's Theorem to the above formula, we get

$$\begin{eqnarray} \Theta ^{(i\!A\times j B);(k\!A\times l B)}_{bb} &=& \frac{1}{\nu _b^2} \big[ \nu _b^2 \big\langle C_b ^{(i\!A \times j\!B)} \big\rangle \big\langle C_b ^{(k\!A \times l\!B)} \big\rangle \nonumber \\ && + \nu _b \big\langle C_b ^{(i\!A \times k\!A)} \big\rangle \big\langle C_b ^{(j\!B \times l\!B)} \big\rangle \nonumber \\ && + \nu _b \big\langle C_b ^{(i\!A \times l\!B)} \big\rangle \big\langle C_b ^{(k\!A \times j\!B)} \big\rangle \big] \nonumber \\ && - \big\langle C_b ^{(i\!A \times j\!B)} \big\rangle \big\langle C_b ^{(k\!A \times l\!B)} \big\rangle \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} \hspace{4.5pc}& = & \frac{1}{\nu _b} \big[ \big\langle C_b ^{(i\!A \times k\!A)} \big\rangle \big\langle C_b ^{(j\!B \times l\!B)} \big\rangle \nonumber \\ \hspace{4.5pc}&& + \big\langle C_b ^{(i\!A \times l\!B)} \big\rangle \big\langle C_b ^{(k\!A \times j\!B)} \big\rangle \big]. \end{eqnarray} \tag{ A7 }$$

Now we will examine various cases of the above variance.

Type 1. All four splits different, i ≠ j ≠ k ≠ l. In this case, all four ensemble averages are cross-power spectra and evaluate to the underlying spectrum C_b, giving

$$\begin{equation} \Theta ^{(i\!A\times j B);(k\!A\times l B); {\rm Type 1}}_{bb} = \frac{2}{\nu _b} C_b^2. \end{equation} \tag{ A8 }$$

Type 2a. One split in common, of channel A (i = k ≠ l ≠ j). Here, the first ensemble average in Equation (A7) is an auto-frequency auto-power spectrum which evaluates out to (C_b + N^AA_b), where N^AA_b is the noise in the channel A data split. This gives

$$\begin{equation} \Theta ^{(i\!A\times j B);(k\!A\times l B); {\rm Type 2a}}_{bb} = \frac{1}{\nu _b} \big(2 C_b^2+ C_b N_b^{A\!A}\big). \end{equation} \tag{ A9 }$$

Type 2b. One split in common, of channel B, (l = j ≠ k ≠ l):

$$\begin{equation} \Theta ^{(i\!A\times j B);(k\!A\times l B) ; {\rm Type 2b}}_{bb} = \frac{1}{\nu _b} \big(2 C_b^2+ C_b N_b^{B\!B}\big). \end{equation} \tag{ A10 }$$

Type 2c. One split in common, one of channel A, the other of channel B (i = l or k = j):

$$\begin{equation} \Theta ^{(i\!A\times j B);(k\!A\times l B); {\rm Type 2c}}_{bb} = \frac{1}{\nu _b} \big(2 C_b^2+ C_b N_b^{A\!B}\big), \end{equation} \tag{ A11 }$$

where N^AB_b represents any noise that is correlated across the two frequency channels.

Type 3a. Two splits in common (i = k, j = l). In this case, we have

$$\begin{eqnarray} \Theta ^{(i\!A\times j B);(k\!A\times l B); {\rm Type 3a}}_{bb} &=& \frac{1}{\nu _b} \big[2 C_b^2+ C_b \big(N_b^{A\!A} +N_b^{B\!B}\big) \nonumber \\ && + N_b^{A\!A} N_b^{B\!B}\big]. \end{eqnarray} \tag{ A12 }$$

Type 3b. Two splits in common (i = l, j = k). This evaluates to

$$\begin{equation} \Theta ^{(i\!A\times j B);(k\!A\times l B); {\rm Type 3b}}_{bb} = \frac{1}{\nu _b} \left[2 C_b^2+ 2 C_b N_b^{A\!B} + N_b^{A\!B} N_b^{A\!B}\right]. \end{equation} \tag{ A13 }$$

Next, we turn our attention to the variance of the mean spectrum estimator, $\widehat{C}_b^{A\times B}$, defined in Equation (A1):

$$\begin{equation} \Theta ^{{\rm mean;mean}} = \big\langle \big(\widehat{C}^{A\times B}_b -\big\langle \widehat{C}_b^{A\times B}\big\rangle \big) \big(\widehat{C}_b^{A\times B} -\big\langle \widehat{C}_b^{A\times B} \big\rangle \big)\big\rangle. \end{equation} \tag{ A14 }$$

Expanding one of the terms out, we have

$$\begin{eqnarray} \Theta ^{{\rm mean;mean}} &=& \frac{1}{n_d(n_d-1)} \sum _{i,j; i\ne j }\big\langle \big(\widehat{C}_b^{(i\!A\times j B)} -\big\langle \widehat{C}_b^{(i\!A\times j B)}\big\rangle \big) \nonumber \\ && \times \big(\widehat{C}_b ^{A\times B} -\big\langle \widehat{C}_b^{A\times B} \big\rangle \big)\big\rangle. \end{eqnarray} \tag{ A15 }$$

Each of the n_d(n_d − 1) terms in the sum is statistically equivalent to the other and should evaluate to be the same. Hence,

$$\begin{equation} \Theta ^{{\rm mean;mean}} = \Theta ^{(i\!A\times j B);{\rm mean}}, \end{equation} \tag{ A16 }$$

which can be expanded out as

$$\begin{equation} \Theta ^{{\rm mean;mean}} = \frac{1}{n_d(n_d-1)} \sum _{k,l; k\ne l }\Theta ^{(i\!A \times j B);(k\!A \times l B)}. \end{equation} \tag{ A17 }$$

The last piece of the calculation is to figure out how many times each of the types defined above appear in this sum. Since i and j are fixed, it is apparent that for Type 1, k can be chosen in (n_d − 2) ways and l in (n_d − 3)ways. So, Type 1 appears (n_d − 2)(n_d − 3) times. For Type 2a, k is fixed by equality with i, so l can be chosen only in (n_d − 2) ways. Therefore, Type 2a appears (n_d − 2) times. From symmetry, Type 2b also appears (n_d − 2) times. For Type 2c, either i = l or k = j which gives 2(n_d − 2) terms. Finally, for each of Types 3a and 3b, all four indices are completely fixed, so they appear only once. Putting these together, we have

$$\begin{eqnarray*} \nonumber \Theta ^{{\rm mean;mean}} & = & \frac{1}{n_d(n_d-1)} \left[ (n_d-2)(n_d-3) \frac{2}{\nu _b} C_b^2 \right.\\ \nonumber & & + (n_d-2) \frac{1}{\nu _b} \big(2 C_b^2+ C_b N_b^{A\!A}\big)\\ \nonumber & & + (n_d-2) \frac{1}{\nu _b} \big(2 C_b^2+ C_b N_b^{B\!B}\big)\\ \nonumber & & + 2(n_d-2) \frac{1}{\nu _b} \big(2 C_b^2+ C_b N_b^{A\!B}\big)\\ \nonumber & & + \frac{1}{\nu _b} \left\lbrace 2 C_b^2+ C_b \big(N_b^{A\!A}+N_b^{B\!B}\big) + N_b^{A\!A} N_b^{B\!B}\right\rbrace \\ \nonumber & & + \left.\frac{1}{\nu _b} \left\lbrace 2 C_b^2+ 2 C_b N_b^{A\!B} + N_b^{A\!B} N_b^{A\!B}\right\rbrace \right] \\ \nonumber & = & \frac{1}{\nu _b} \frac{1}{n_d(n_d-1)} \big[ n_d(n_d-1) 2 C_b^2 \\ \nonumber && + (n_d-1) C_b \big(N_b^{A\!A}+N_b^{B\!B}\big) + 2 (n_d-1) C_b N_b^{A\!B} \\ \nonumber && +\big(N_b^{A\!A} N_b^{B\!B} + N_b^{A\!B} N_b^{A\!B}\big)\big], \end{eqnarray*}$$

which simplifies to

$$\begin{eqnarray} \Theta ^{{\rm mean;mean}} & = & \frac{1}{\nu _b} \bigg[ 2 C_b^2 + \frac{C_b}{n_d} \big(N_b^{A\!A}+N_b^{B\!B}\big) + \frac{2}{n_d} C_b N_b^{A\!B} \nonumber \\ && + \frac{\big(N_b^{A\!A}N_b^{B\!B}+N_b^{A\!B}N_b^{A\!B}\big)}{n_d (n_d-1)} \bigg]. \end{eqnarray} \tag{ A18 }$$

Note that if there is only one channel A ≡ B, this reduces to the familiar form (see Equation (9) in Fowler et al. 2010)

$$\begin{eqnarray} \Theta ^{{\rm mean;mean}} & = & \frac{1}{\nu _b} \left[ 2 C_b^2 + 4 \frac{C_b}{n_d} N_b^{A\!A} + 2 \frac{\big(N_b^{A\!A}\big)^2}{n_d (n_d-1)} \right].\qquad \end{eqnarray} \tag{ A19 }$$

Finally, if we have n_p patches with equal noise, the above variance should be divided by n_p. In practice, a weighted combination of the individual patch variances is used. Under the assumption of isotropic noise and filtering, the effect of the data window, W, can be taken into account by correcting the number of modes as ν_b → ν_bw²₂/w₄, where w_n represents the nth moment of the window. There is a small additional correction due to the anisotropic nature of the noise and filtering, which we calibrate against Monte Carlo simulations.

For estimating parameters from the three spectra, C^148×148_b, C^148×218_b, and C^218×218_b a joint likelihood function is written as

$$\begin{eqnarray} - 2 \ln {\it L} &=& \left(\begin{array}{@{}c@{}} C_b^{148\times 148} - C_b^{148\times 148, \rm {theory}}\\[3pt] C_b^{148\times 218} - C_b^{148\times 218, \rm {theory}}\\[3pt] C_b^{218\times 218} - C_b^{218\times 218, \rm {theory}} \end{array} \right)^T \nonumber \\ &&\times {\bf \Sigma }^{-1} \left(\begin{array}{@{}c@{}}C_b^{148\times 148} - C_b^{148\times 148, \rm {theory}}\\[3pt] C_b^{148\times 218} - C_b^{148\times 218, \rm {theory}}\\[3pt] C_b^{218\times 218} - C_b^{218\times 218, \rm {theory}} \end{array} \right),\qquad \end{eqnarray} \tag{ A20 }$$

where

$$\begin{equation} {\bf \Sigma } = \left(\begin{array}{@{}ccc@{}}\Sigma _b^{148-148,148-148} & \Sigma _b^{148-148,148-218} & \Sigma _b^{148-148,218-218} \\[3pt] \Sigma _b^{148-218,148-148} & \Sigma _b^{148-218,148-218} & \Sigma _b^{148-218,218-218} \\[3pt] \Sigma _b^{218-218,148-148} & \Sigma _b^{218-218,148-218} & \Sigma _b^{218-218,218-218} \end{array} \right) \end{equation} \tag{ A21 }$$

with

$$\begin{equation} \Sigma _b^{A-B,C-D} = \big\langle \big(\widehat{C}^{A\times B}_b -\big\langle \widehat{C}_b^{A \times B}\big\rangle \big)\big(\widehat{C}_b^{C\times D} -\big\langle \widehat{C}_b^{C\times D}\big)\big\rangle \big\rangle, \end{equation} \tag{ A22 }$$

and can be computed in a similar fashion as above.