text:
# First Year Wilkinson Microwave Anisotropy Probe (WMAPWMAP is the result of a partnership between Princeton University and NASA’s Goddard Space Flight Center. Scientific guidance is provided by the WMAP Science Team. ) Observations: The Angular Power Spectrum

## 1 INTRODUCTION

The Wilkinson Microwave Anisotropy Probe (WMAP) mission was designed to measure the CMB anisotropy with unprecedented precision and accuracy on angular scales from the full sky to several arc minutes by producing maps at five frequencies from 23 to 94 GHz. The WMAP satellite mission (Bennett et al., 2003a) employs a matched pair of 1.4m telescopes (Page et al., 2003c) observing two areas on the sky separated by $``$141. A differential radiometer (Jarosik et al., 2003a) with a total of 10 feeds for each of the two sets of optics (Barnes et al., 2003; Page et al., 2003c) measures the difference in sky brightness between the two sky pixels. The satellite is deployed at the Earth-Sun Lagrange point, $`L_2`$, and observes the sky with a compound spin and precession that covers the full sky every six months. The differential data are processed on the ground to produce full sky maps of the CMB anisotropy (Hinshaw et al., 2003b).

Full sky maps provide the smallest record of the CMB anisotropy without loss of information. They permit a wide variety of statistics to be computed from the data – one of the most fundamental is the angular power spectrum of the CMB. Indeed, if the temperature fluctuations are Gaussian, with random phase, then the angular power spectrum provides a complete description of the statistical properties of the CMB. Komatsu et al. (2003) have analyzed the first-year WMAP sky maps to search for evidence of non-Gaussianity and find none, aside from a modest level of point source contamination which we account for in this paper. Thus, the measured power spectrum may be compared to predictions of cosmological models to develop constraints on model parameters.

This paper presents the angular power spectrum obtained from the first-year WMAP sky maps. Companion papers present the maps and an overview of the basic results (Bennett et al., 2003b), and describe the foreground removal process that precedes the power spectrum analysis (Bennett et al., 2003c). Spergel et al. (2003), Verde et al. (2003), Peiris et al. (2003), Kogut et al. (2003), and Page et al. (2003b) discuss the implications of the WMAP power spectrum for cosmological parameters and carry out a joint analysis of the WMAP spectrum together with other CMB data and data from large-scale structure probes. Hinshaw et al. (2003b), Jarosik et al. (2003b), Page et al. (2003a), Barnes et al. (2003), and Limon et al. (2003) discuss the data processing, the radiometer performance, the instrument beam characteristics and the spacecraft in-orbit performance, respectively.

A sky map $`\mathrm{\Delta }T(𝐧)`$ defined over the full sky can be decomposed in spherical harmonics

$$\mathrm{\Delta }T(𝐧)=\underset{l>0}{}\underset{m=l}{\overset{l}{}}a_{lm}Y_{lm}(𝐧)$$
(1)
with

$$a_{lm}=𝑑\mathrm{\Omega }_𝐧\mathrm{\Delta }T(𝐧)Y_{lm}^{}(𝐧),$$
(2)
where $`𝐧`$ is a unit direction vector, and $`Y_{lm}(𝐧)`$ are the spherical harmonic functions evaluated in the direction $`𝐧`$.

If the CMB temperature fluctuation $`\mathrm{\Delta }T`$ is Gaussian distributed, then each $`a_{lm}`$ is an independent Gaussian deviate with

$$a_{lm}=0,$$
(3)
and

$$a_{lm}a_{l^{}m^{}}^{}=\delta _{ll^{}}\delta _{mm^{}}C_l,$$
(4)
where $`C_l`$ is the ensemble average power spectrum predicted by models, and $`\delta `$ is the Kronecker symbol. The actual power spectrum realized in our sky is

$$C_l^{\mathrm{sky}}=\frac{1}{2l+1}\underset{m=l}{\overset{l}{}}\left|a_{lm}\right|^2.$$
(5)
In the absence of noise, and with full sky coverage, the right hand side of equation (5) provides an unbiased estimate of the underlying theoretical power spectrum, which is limited only by cosmic variance. However, realistic CMB anisotropy measurements contain noise and other sources of error that cause the quadratic estimator in equation (5) to be biased. In addition, while WMAP measures the anisotropy over the full sky, the data near the Galactic plane are sufficiently contaminated by foreground emission that only a portion of the sky ($``$85%) can be used for CMB power spectrum estimation. Thus, the integral in equation (2) cannot be evaluated as such, and other methods must be found to estimate $`C_l`$.

In Appendix A, we review two methods that have appeared in the literature for estimating the angular power spectrum in the presence of instrument noise and sky cuts. The first (Hivon et al., 2002) is a quadratic estimator that evaluates equation (2) on the cut sky yielding a “pseudo power spectrum” $`\stackrel{~}{C}_l`$. The ensemble average of this quantity is related to the true power spectrum, $`C_l`$ by means of a mode coupling matrix $`G_{ll^{}}`$ (Hauser & Peebles, 1973). The second method (Oh et al., 1999) uses a maximum likelihood approach optimized for fast evaluation with WMAP-like data. In Appendix A we demonstrate that the two methods produce consistent results.

A quadratic estimator offers the possibility of computing both an “auto-power” spectrum, proportional to $`_m|a_{lm}|^2`$, and a “cross-power” spectrum, proportional to $`_ma_{lm}^ia_{lm}^j`$, where the $`a_{lm}`$ coefficients are estimated from two independent CMB maps, $`i`$ and $`j`$. This latter form has the advantage that, if the noise in the two maps is uncorrelated, the quadratic estimator is not biased by the noise. For all cosmological analyses, we use only the cross-power spectra between statistically independent channels. As a result, the angular spectra are, for all intents and purposes, independent of the noise properties of an individual radiometer. This is analogous to interferometric data which exhibits a high degree of immunity to systematic errors. The precise form of the estimator we use is given in Appendix A.

The plan of this paper is as follows. In §2 we review the properties of the WMAP instrument and how they affect the derived power spectrum. In §3 we present results for the angular power spectra obtained from individual pairs of radiometers, the cross-power spectra, and examine numerous consistency tests. In §4, in preparation for generating a final combined power spectrum, we present the full covariance matrix of the cross-power spectra. In §5 we present the methodology used to produce the final combined power spectrum and its covariance matrix. In §6 we compare the WMAP first-year power spectrum to a compilation of previous CMB measurements and to a prediction based on a combination of previous CMB data and the 2dFGRS data. We summarize our results in §7 and outline the power spectrum data products being made available through the Legacy Archive for Microwave Background Data Analysis (LAMBDA). Appendix A reviews two methods used to estimate the angular power spectrum from CMB maps. Appendix B describes how we account for point source contamination. Appendix C presents our approach to combining multi-channel data. Appendix D describes how the foreground mask correlates multipole moments in the Fisher matrix, and Appendix E collects some useful properties of the spherical harmonics.

## 2 INSTRUMENTAL PROPERTIES

The WMAP instrument is composed of 10 “differencing assemblies” (DAs) spanning 5 frequencies from 23 to 94 GHz (Bennett et al., 2003a). The 2 lowest frequency bands (K and Ka) are primarily Galactic foreground monitors, while the 3 highest (Q, V, and W) are primarily cosmological bands. There are 8 high frequency differencing assemblies: Q1, Q2, V1, V2, and W1 through W4. Each DA is formed from two differential radiometers which are sensitive to orthogonal linear polarization modes; the radiometers are designated 1 or 2 (e.g., V11 or W12) depending on which polarization mode is being sensed.

The temperature measured on the sky is modified by the properties of the instrument. The most important properties that affect the angular power spectrum are finite resolution and instrument noise. Let $`C_l^{ii^{}}`$ denote the auto or cross-power spectrum evaluated from two sky maps, $`i`$ and $`i^{}`$, where $`i`$ is a DA index. Further, define the shorthand $`𝐢(i,i^{})`$ to denote a pair of indices, e.g., (Q1,V2). This spectrum will have the form

$$C_l^𝐢=w_l^𝐢C_l^{\mathrm{sky}}+N_l^𝐢,$$
(6)
where $`w_l^𝐢b_l^ib_l^i^{}p_l^2`$ is the window function that describes the combined smoothing effects of the beam and the finite sky map pixel size. Here $`b_l^i`$ is the beam transfer function for DA $`i`$, given by Page et al. (2003a) \[note that they reserve the term “beam window function” for $`(b_l^i)^2`$\], and $`p_l`$ is the pixel transfer function supplied with the HEALPix package (Górski et al., 1998). $`N_l^𝐢`$ is the noise spectrum realized in this particular measurement. On average, the observed spectrum estimates the underlying power spectrum, $`C_l`$,

$$C_l^𝐢=w_l^𝐢C_l+N_l^i\delta _{ii^{}},$$
(7)
where $`N_l^i`$ is the average noise power spectrum for differencing assembly $`i`$, and the Kronecker symbol indicates that the noise is uncorrelated between differencing assemblies. To estimate the underlying power spectrum on the sky, $`C_l`$, the effects of the noise bias and beam convolution must be removed. The determination of transfer functions and noise properties are thus critical components of any CMB experiment.

In §2.1 we summarize the results of Page et al. (2003a) on the WMAP window functions and their uncertainties. We propagate these uncertainties through to the final Fisher matrix for the angular power spectrum. In §2.2 we present a model of the WMAP noise properties appropriate to power spectrum evaluation. For cross power spectra ($`ii^{}`$ above), the noise bias term drops out of equation (7) if the noise between the two DAs is uncorrelated. These cross-power spectra provide a nearly optimal estimate of the true power spectrum, essentially independent of errors in the noise model, thus we use them exclusively in our final power spectrum estimate. The noise model is primarily used to propagate noise errors through the analysis, and to test a variety of different power spectrum estimates for consistency with the combined cross-power spectrum.

### 2.1 Window Functions

As discussed in Page et al. (2003a), the instrument beam response was mapped in flight using observations of the planet Jupiter. The signal to noise ratio is such that the response, relative to the peak of the beam, is measured to approximately $`35`$ dB in W band, the band with the highest angular resolution. The beam widths, measured in flight, range from $`0\text{.}^{}82`$ at K band down to $`0\text{.}^{}20`$ in some of the W band channels (FWHM). Maps of the full two-dimensional beam response are presented in Page et al. (2003a), and are available with the WMAP first-year data release. The radial beam profiles obtained from these maps have been fit to a model consisting of a sum of Hermite polynomials that accurately characterize the main Gaussian lobe and small deviations from it. The model profiles are then Legendre transformed to obtain the beam transfer functions $`b_l^i`$ for each DA $`i`$. Full details of this procedure are presented in Page et al. (2003a), and the resulting transfer functions are also provided in the first-year data release. We have chosen to normalize the transfer function to 1 at $`l=1`$ because WMAP calibrates its intensity response using the modulation of the CMB dipole ($`l=1`$). This effectively partitions calibration uncertainty from window function uncertainties.

The beam processing described above provides a straightforward means of propagating the noise uncertainty directly from the time-ordered data through to the final transfer functions. The result is the covariance matrix $`\mathrm{\Sigma }_{b,ll^{}}^i`$ for the normalized transfer function. Plots of the diagonal elements of $`\mathrm{\Sigma }_{b,ll^{}}^i`$ are presented in Page et al. (2003a). The fractional error in the transfer functions $`b_l^i`$ are typically 1-2% in amplitude. In the end, these window function uncertainties dominate the small off-diagonal elements of the final covariance matrix for the combined power spectrum (see §4). Additional observations of Jupiter will reduce these uncertainties.

An additional source of error in our treatment of the beam response arises from non-circularity of the main beam. The effects of this non-circularity are mitigated by WMAP’s scan strategy which results in most sky pixels being observed over a wide range of azimuth angles. The effective beam response on the sky is thus largely symmetrized. We estimate that the effects of imperfect symmetrization produce window function errors of $`<1`$% relative to a perfectly symmetrized beam window function (Page et al., 2003a; Hinshaw et al., 2003b). This error is well within the formal uncertainty given in $`\mathrm{\Sigma }_{b,ll^{}}^i`$. In §5 we infer the optimal power spectrum by combining the 28 cross-power spectra measured by the 8 high frequency DAs Q1 through W4. As part of this process we marginalize over the window function uncertainty, which automatically propagates these errors into the final covariance matrix for the combined power spectrum. Both the combined power spectrum and the corresponding Fisher matrix are part of the first-year data release.

### 2.2 Instrument Noise Properties

The noise bias term in equation (7) is the noise per $`a_{lm}`$ mode on the sky. If auto-power spectra are used in the final power spectrum estimate, the noise bias term must be known very accurately because it exponentially dominates the convolved power spectrum at high $`l`$. If only cross correlations are used, the noise bias is only required for propagating errors. Our final best spectrum is based only on cross correlations, and is independent of this term. However, as an independent check of our results, we evaluate the maximum likelihood spectrum based on a combined Q+V+W sky map. The noise bias must be estimated accurately for this application.

In the limit that the time-ordered instrument noise is white, the noise bias is a constant, independent of $`l`$. If the time-ordered noise has a $`1/f`$ component, the bias term will rise at low $`l`$. In this subsection we estimate the noise bias properties for each of the high frequency WMAP radiometers based on the time-ordered noise properties presented in Jarosik et al. (2003b). While the WMAP radiometer noise is nearly white by design (Jarosik et al., 2003a) with $`1/f`$ knee frequencies of less than 10 mHz for 9 out of 10 differencing assemblies, one of the radiometers (W41) has a $`1/f`$ knee frequency of $``$45 mHz. The latter is large enough that the deviations of $`N_l^i`$ from a constant must be accounted for.

The most reliable way to estimate the effects of $`1/f`$ noise on the measured power spectra is by Monte Carlo simulation. Using the pipeline simulator discussed in Hinshaw et al. (2003b) we have generated a library of noise maps with flight-like properties. Specifically we have included flight-like $`1/f`$ noise in the simulated time-ordered data, and have run each full-year realization through the map-making pipeline, including the baseline pre-whitening discussed in Hinshaw et al. (2003b). We evaluate the power spectra of these maps using the quadratic estimator described in Appendix A with 3 different pixel weighting schemes. (See Appendix A.1.2 for definitions of the weights, and the $`l`$ range in which each is used.) We define the effective noise as a function of $`l`$ based on fits to these Monte Carlo noise spectra. For the analyses in this paper, we fit the spectra to a model of the form

$$\mathrm{ln}N_l^i=\left[\underset{n=0}{\overset{n_{\mathrm{max}}}{}}c_n^i(\mathrm{ln}l)^n\right]^1$$
(8)
where the $`c_n^i`$ are fit coefficients given in Table 1, with $`n_{\mathrm{max}}=3`$ for $`l<200`$, and $`n_{\mathrm{max}}=1`$ for $`l>200`$.

Figure 1 shows the noise spectrum derived from the simulations for each of the 8 high frequency DAs, using uniform weighting over the entire $`l`$ range. For comparison, we also plot an estimate of the CMB power spectrum from §5 in grey. Note that the W4 spectrum is the only one of this set to exhibit deviations from white noise in an $`l`$ range where the signal-to-noise is relatively low, and we believe this simulation slightly over-estimates the $`1/f`$ noise in the flight W4 differencing assembly (Hinshaw et al., 2003a).

### 2.3 Systematic Errors

Hinshaw et al. (2003a) present limits on systematic errors in the first-year sky maps. They consider the effects of absolute and relative calibration errors, artifacts induced by environmental disturbances (thermal and electrical), errors from the map-making process, pointing errors, and other miscellaneous effects. The combined errors due to relative calibration errors, environmental effects, and map-making errors are limited to $`<15`$ $`\mu `$K2 ($`2\sigma `$) in the quadrupole moment $`C_2`$ in any of the 8 high-frequency DAs. Tighter limits are placed on higher-order moments. We conservatively estimate the absolute calibration uncertainty in the first-year WMAP data to be 0.5%.

Random pointing errors are accounted for in the beam mapping procedure; the beam transfer functions presented by Page et al. (2003a) incorporate random pointing errors automatically. A systematic pointing error of $``$10$`^{\prime \prime }`$at the spin period is suspected in the quaternion solution that defines the spacecraft pointing. This is much smaller than the smallest beam width ($``$12$`^{}`$at W band), and we estimate that it would produce $`<`$1% error in the angular power spectrum at $`l=1000`$, thus we do not attempt to correct for this effect. Barnes et al. (2003) place limits on spurious contributions due to stray light pickup through the far sidelobes of the instrument. They place limits of $`<`$10 $`\mu `$K2 on spurious contributions to $`C_l`$, at Q through W band, due to far sidelobe pickup.

A detailed model of Galactic foreground emission based on the first-year WMAP data is presented by Bennett et al. (2003c) and is summarized in §3.1. We show that diffuse foreground emission is a modest source of contamination at large angular scales ($`\genfrac{}{}{0pt}{}{_>}{^{}}`$2). Systematic errors on these angular scales are negligible compared to the (modest) level of foreground emission. On smaller angular scale ($`\genfrac{}{}{0pt}{}{_<}{^{}}`$2), the 1-3% uncertainty in the individual beam transfer functions is the largest source of uncertainty, while for multipole moments greater than $``$600, random white noise from the instrument is the largest source of uncertainty.

## 3 THE DATA

Figure 2 shows the cross-power spectra obtained from all 28 combinations of the 8 differencing assemblies Q1 through W4 using the quadratic estimator described in Appendix A.1. These spectra have been evaluated with the Kp2 sky cut described in Bennett et al. (2003c). The spectra are color coded by effective frequency, $`\sqrt{\nu ^i\nu ^i^{}}`$, where $`\nu ^i`$ is the frequency of differencing assembly $`i`$. The low frequency (41 GHz) data are shown in red, the high frequency (94 GHz) data in blue, with intermediate frequencies following the colors of the rainbow. The top panel shows $`l(l+1)C_l/2\pi `$ in $`\mu `$K2, while the bottom panel plots the ratio of each channel to our final combined spectrum presented in §5. The top panel shows a very robust measurement of the first acoustic peak with a maximum near $`l220`$ and a shape that is consistent with the predictions of adiabatic fluctuation models. There is also a clear indication of the rise to a second peak at $`l540`$. See Page et al. (2003b) for an analysis and discussion of the peaks and troughs in the first-year WMAP power spectrum.

The red data in the top panel show very clearly that the low frequency data are contaminated by diffuse Galactic emission at low $`l`$ and by point sources at higher $`l`$. The higher frequency data show less apparent contamination, consistent with the foreground emission being dominated by radio emission, rather than thermal dust emission, as expected in this frequency range.

### 3.1 Galactic and Extragalactic Foregrounds

Bennett et al. (2003c) present a detailed model of the Galactic foreground emission based on a Maximum Entropy analysis of all 5 WMAP frequency bands, in combination with external tracer templates. They demonstrate that the emission is well modeled by three distinct emission components. 1) Synchrotron emission from cosmic ray electrons, with a steeply falling spectrum in the WMAP frequency range: $`T_A(\nu )\nu ^\beta `$ with $`\beta <3`$, steepening with increasing frequency. 2) Free-free emission from the ionized interstellar medium that is well traced by H$`\alpha `$ emission in regions where the dust extinction is low. 3) Thermal emission from interstellar dust grains with an emissivity index $`\beta 2.2`$. The model has a Galactic signal minimum between V and W band.

In principal we could subtract the above model from each WMAP channel and recompute the power spectrum. However, since the model is based on WMAP data that have been smoothed to an angular resolution of $`1\text{.}^{}0`$, the resulting maps would have complicated noise properties. For the purposes of power spectrum analysis, we adopt a more straightforward approach based on fitting foreground tracer templates to the Q, V, and W band data. The details of this procedure, the resulting fit coefficients, and a comparison of the fits to the Maximum Entropy model are given in Bennett et al. (2003c). They estimate that the template model removes $``$85% of the foreground emission in Q, V, and W bands and that the remaining emission constitutes less than $``$2% of the CMB variance (up to $`l=200`$) in Q band, and $`\genfrac{}{}{0pt}{}{_<}{^{}}`$1% of the CMB variance in V and W bands.

The contribution from extragalactic radio sources has been analyzed in three separate ways. Bennett et al. (2003c) directly fit for sources in the sky maps. The result of this analysis is that we have identified 208 sources in the WMAP data with sufficient signal to noise ratio to pass the detection criterion (we estimate that $``$5 of these are likely to be spurious). The derived source count law is consistent with the following model for the power spectrum of the unresolved sources

$$C_l^{𝐢,src}=A\left(\frac{\nu _i}{\nu _0}\right)^\beta \left(\frac{\nu _i^{}}{\nu _0}\right)^\beta w_l^𝐢,$$
(9)
with $`A=0.015`$ $`\mu `$K2 sr (measured in thermodynamic temperature), $`\beta =2.0`$, and $`\nu _045`$ GHz. Komatsu et al. (2003) evaluate the bispectrum of the WMAP data and are able to fit a non-Gaussian source component based on a particular configuration of the bispectrum. They find the same source model, equation (9), fits the bispectrum data. For the remainder of this section, we adopt this model for correcting the cross-power spectra. At Q band (41 GHz) the correction to $`l(l+1)C_l/2\pi `$ is 868 and 3468 $`\mu `$K2 at $`l`$=500 and 1000, respectively. At W band (94 GHz), the correction is only 31 and 126 $`\mu `$K2 at the same $`l`$ values. For comparison, the CMB power in this $`l`$ range is $``$ 2000 $`\mu `$K2. Later, when we derive a final combined spectrum from the multi-frequency data, we adopt equation (9) as a model with $`A`$ as a free parameter. We simultaneously fit for a combined CMB spectrum and source amplitude and marginalize over the residual uncertainty in A. The best-fit source amplitude from this process is consistent with the other two methods.

Figure 3 shows the cross-power spectra obtained from the same 28 combinations as in Figure 2, this time with the Galactic template model and source model subtracted. The bottom panel of the Figure shows the ratio of the 28 channels to the combined spectrum obtained in §5. The 28 cross-power spectra are consistent with each other at the 5 to 20% level over the $`l`$ range 2 – 500. Similar scatter is seen in Monte Carlo simulations of an ensemble of 28 cross-power spectra with WMAP’s beam and noise properties. The only significant deviation lies in the Q band data at low $`l`$ which is $``$10% below the higher frequency bands at $`l<20`$. This is consistent with the accuracy estimated above for the Galactic template model, see also Figure 11 of Bennett et al. (2003c) for images of the maps after Galactic template subtraction. Since the WMAP data are not noise limited at low $`l`$, we use only V and W band data in the final combined spectrum for $`l<100`$.

The subtraction of the source model, equation (9), brings the Q band spectrum into good agreement with the other cross-power spectra up to $`l500`$. At higher $`l`$, the Q band data contributes very little to the final combined spectrum because the (normalized) Q band window function has dropped to less than 5% (Page et al., 2003a). As discussed in §5, we marginalize over the source amplitude uncertainty, $`\delta A`$, when obtaining the final power spectrum estimate and associated covariance matrix. Thus the uncertainty is also accounted for in subsequent cosmological parameter fits (Verde et al., 2003; Spergel et al., 2003; Peiris et al., 2003).

Figure 4 shows a close-up of the 28 cross-power spectra in Figure 3 up to $`l=100`$. The top panel shows the raw (unbinned) data which has correlations of $`<`$2% between neighboring points is this $`l`$ range (see §4). These data are strikingly consistent with each other and support the conclusion that systematic errors at low $`l`$ are insignificant. To assess the level of scatter that does exist between the spectra, we have generated a Monte Carlo simulation in which we compute the rms scatter among the 28 spectra at each $`l`$, relative to the measured power. The bottom panel of Figure 4 shows the results of this simulation, averaged over 1000 realizations, compared to the relative rms scatter in the data. The agreement is excellent, indicating that the uncertainty in the measured power spectrum in this $`l`$ range is a few percent and is consistent with a combination of instrument noise and mode coupling due to the 15% sky cut.

Another striking feature is the low amplitude of the observed quadrupole, and the sharp rise in power, almost linear in $`l`$, to $`l=5`$. Bennett et al. (2003b) quote a value for the rms quadrupole amplitude, $`Q_{rms}=\sqrt{(5/4\pi )C_2}=8\pm 2`$ $`\mu `$K, where the uncertainty is largely due to Galactic model uncertainty. This is consistent with the amplitude measured by the COBE-DMR experiment, $`Q_{rms}=10_4^{+7}`$ $`\mu `$K (Bennett et al., 1996). The fast, nearly linear rise to $`l=5`$ produces an angular correlation function with essentially no power on angular scales $`\genfrac{}{}{0pt}{}{_>}{^{}}`$60, again in excellent agreement with the COBE-DMR correlation function (Bennett et al., 2003b; Hinshaw et al., 1996). In the context of a standard $`\mathrm{\Lambda }`$CDM model, the probability of observing this little power on scales greater than 60 is $`2\times 10^3`$ (Bennett et al., 2003b).

## 4 THE FULL COVARIANCE MATRIX

In §5 we derive our best estimate of the angular power spectrum by optimally combining the 28 cross-power spectra discussed above. The procedure for combining spectra requires the full covariance matrix of the individual cross-power spectra – in this section we outline the salient features of this matrix. There are six principal sources of variance for the measured spectra, $`C_l^𝐢`$: cosmic variance, instrument noise, mode coupling due to the foreground mask, point source subtraction errors, uncertainty in the beam window functions, and an overall calibration uncertainty. We ignore uncertainties in the diffuse foreground correction since they are everywhere sub-dominant to the cosmic variance uncertainty (see §3.1).

We may write the covariance matrix as

$$(\mathrm{\Sigma }_{\mathrm{full}})_{ll^{}}^{\mathrm{𝐢𝐣}}=[C_l^𝐢(C_l+AS^𝐢)w_l^𝐢][C_l^{}^𝐣(C_l^{}+AS^𝐣)w_l^{}^𝐣],$$
(10)
where the angle brackets represent an ensemble average, $`C_l`$ is the true underlying power spectrum, $`w_l^𝐢`$ is the window function of spectrum $`𝐢`$, and $`AS^𝐢`$ is the point source contribution to spectrum $`𝐢`$. Here we have defined a point source spectral function, $`S^𝐢`$, as

$$S^𝐢\left(\frac{\nu _i}{\nu _0}\right)^\beta \left(\frac{\nu _i^{}}{\nu _0}\right)^\beta ,$$
(11)
where $`\nu _i`$, $`\nu _0`$, and $`\beta `$ are as defined after equation (9). Note that $`(\mathrm{\Sigma }_{\mathrm{full}})_{ll^{}}^{\mathrm{𝐢𝐣}}`$ is symmetric in both $`(ll^{})`$ and $`(\mathrm{𝐢𝐣})`$.

In the process of forming the combined spectrum we will estimate a best-fit point source amplitude, $`\overline{A}`$, and subtract the corresponding source contribution from each spectrum $`𝐢`$. We thus rewrite $`\mathrm{\Sigma }_{\mathrm{full}}`$ as

$$(\mathrm{\Sigma }_{\mathrm{full}})_{ll^{}}^{\mathrm{𝐢𝐣}}=[\overline{C}_l^𝐢(C_l+\delta AS^𝐢)w_l^𝐢][\overline{C}_l^{}^𝐣(C_l^{}+\delta AS^𝐣)w_l^{}^𝐣],$$
(12)
where $`\overline{C}_l^𝐢C_l^𝐢\overline{A}S^𝐢w_l^𝐢`$ is the source-subtracted spectrum, and $`\delta AA\overline{A}`$ is the residual source amplitude, which we will marginalize over as a nuisance parameter.

We may expand the covariance matrix as

$$\mathrm{\Sigma }_{\mathrm{full}}=\mathrm{\Sigma }_{\mathrm{cv}}+\mathrm{\Sigma }_{\mathrm{mask}}+\mathrm{\Sigma }_{\mathrm{src}}+\mathrm{\Sigma }_{\mathrm{win}}.$$
(13)
We discuss each of these contributions in more detail below.

Cosmic Variance, Noise, and Mode Coupling.— The first two terms, $`\mathrm{\Sigma }_{\mathrm{cv}}+\mathrm{\Sigma }_{\mathrm{mask}}`$, incorporate the combined uncertainty due to cosmic variance, instrument noise, and mode coupling due to the foreground mask,

$$(\mathrm{\Sigma }_{\mathrm{cv}}+\mathrm{\Sigma }_{\mathrm{mask}})_{ll^{}}^{\mathrm{𝐢𝐣}}=[\overline{C}_l^𝐢C_lw_l^𝐢][\overline{C}_l^{}^𝐣C_l^{}w_l^{}^𝐣],$$
(14)
where $`w_l^𝐢`$ is fixed at its measured value. We have split this contribution into two pieces to mimic the procedure we actually use to compute the covariance matrix. As outlined in §5.1, we start with $`\mathrm{\Sigma }_{\mathrm{cv}}`$, then incorporate the effects of point source error and window function error. We do not add the effects of mode coupling, $`\mathrm{\Sigma }_{\mathrm{mask}}`$, until the very end of the computation. We consider this term in more detail in Appendix D.

Point Source Subtraction Errors.— The third term, $`\mathrm{\Sigma }_{\mathrm{src}}`$, is due to uncertainty in the point source amplitude determination

$$(\mathrm{\Sigma }_{\mathrm{src}})_{ll^{}}^{\mathrm{𝐢𝐣}}=S^𝐢w_l^𝐢\sigma _{\mathrm{src}}^2w_l^{}^𝐣S^𝐣,$$
(15)
where $`\sigma _{\mathrm{src}}^2(A\overline{A})^2`$ is the variance in the best-fit amplitude $`\overline{A}`$, and we assume that the frequency dependence, $`S^𝐢`$, is perfectly known. In practice, we do not explicitly evaluate $`\mathrm{\Sigma }_{\mathrm{src}}`$ as given above, rather we employ a method based on marginalizing a Gaussian likelihood function, $`(C_l^𝐢,A)`$, over a nuisance parameter $`A`$. This process, which is discussed in Appendix B, effectively yields $`\mathrm{\Sigma }_{\mathrm{cv}}+\mathrm{\Sigma }_{\mathrm{src}}`$.

Window Function and Calibration Errors.— The fourth term, $`\mathrm{\Sigma }_{\mathrm{win}}`$, is due to uncertainty in the beam window function, $`w_l^𝐢`$. This term arises from fluctuations in the window function which cause the measured spectrum, $`C_l^𝐢`$ to differ from our estimate of the convolved spectrum, $`C_lw_l^𝐢`$, where $`w_l^𝐢`$ is the estimated window function. This contribution has the form

$$(\mathrm{\Sigma }_{\mathrm{win}})_{ll^{}}^{\mathrm{𝐢𝐣}}=C_l\mathrm{\Delta }w_l^𝐢\mathrm{\Delta }w_l^{}^𝐣C_l^{}.$$
(16)
Recall from §2 that $`w_l^𝐢=b_l^ib_l^i^{}p_l^2`$ where $`b_l^i`$ is the beam transfer function for DA $`i`$ and $`p_l^2`$ is the pixel window function. Define $`u_l^i\mathrm{\Delta }b_l^i/b_l^i`$ to be the fractional error in $`b_l^i`$, then to first order in $`u_l^i`$ we have

$$\mathrm{\Delta }w_l^𝐢\mathrm{\Delta }w_l^{}^𝐣=w_l^𝐢u_l^iu_l^{}^j+u_l^iu_l^{}^j^{}+u_l^i^{}u_l^{}^j+u_l^i^{}u_l^{}^j^{}w_l^{}^𝐣.$$
(17)
For WMAP the uncertainty in $`b_l^i`$ is uncorrelated between DAs, thus the above expression reduces to

$$\mathrm{\Delta }w_l^𝐢\mathrm{\Delta }w_l^{}^𝐣=w_l^𝐢\left[\mathrm{\Sigma }_{\mathrm{u},ll^{}}^i(\delta _{ij}+\delta _{ij^{}})+\mathrm{\Sigma }_{\mathrm{u},ll^{}}^i^{}(\delta _{i^{}j}+\delta _{i^{}j^{}})\right]w_l^{}^𝐣,$$
(18)
where we have defined

$$\mathrm{\Sigma }_{\mathrm{u},ll^{}}^iu_l^iu_l^{}^i=(b_l^ib_l^{}^i)^1\mathrm{\Sigma }_{\mathrm{b},ll^{}}^i,$$
(19)
and where $`\mathrm{\Sigma }_{\mathrm{b},ll^{}}^i`$ is the covariance matrix of the beam transfer function for DA $`i`$ given by Page et al. (2003a). When generating the combined spectrum in the next section, we add $`\mathrm{\Sigma }_{\mathrm{win}}`$ to the above contributions, giving $`\mathrm{\Sigma }_{\mathrm{cv}}+\mathrm{\Sigma }_{\mathrm{src}}+\mathrm{\Sigma }_{\mathrm{win}}`$.

The WMAP absolute calibration uncertainty is 0.5%. We do not explicitly incorporate this contribution in the covariance matrix. Instead, we propagate a 0.5% uncertainty into the normalization of the final power spectrum amplitude (Spergel et al., 2003).

## 5 THE COMBINED POWER SPECTRUM

In §3 we demonstrated that the three high frequency bands of WMAP data produced consistent estimates of the angular power spectrum, after a modest correction for diffuse Galactic emission and extragalactic point sources. It is therefore justifiable to combine these data into a single “optimal” estimate of the angular power spectrum of the CMB. In this section, we provide an overview of two methods we use to generate a single combined spectrum. The first is a multi-step process that simultaneously fits the 28 cross-power spectra presented above to a single CMB power spectrum and a point source model, equation (9), while correctly propagating beam and residual point source uncertainties through to a final Fisher matrix. This spectrum constitutes our best estimate of the CMB power spectrum from the first-year WMAP data. The second spectrum, which serves as a cross check of the first, is based on forming a single co-added sky map from the Q1 through W4 maps, and using the quadratic estimator with noise bias subtraction. We compare the two spectra in §5.3.

### 5.1 Method I – Optimal Combination of Cross-Power Spectra

Since this method is relatively complicated, we outline the basic procedure here and relegate the details to Appendices, as indicated. We present the result in §5.3. The steps are as follows.

1. Subtract best-fit Galactic foreground templates from each of the maps Q1 through W4, using the coefficients given in Table 3 of Bennett et al. (2003c).
2. Evaluate the 28 cross-power spectra from the maps Q1 through W4, where each spectrum has been evaluated using the quadratic estimator of Appendix A.1 with the weighting scheme defined in Appendix A.1.2.
3. Collect the noise bias estimate, $`N_l^i`$, for each DA from §2.2. These estimates are used in the calculation of the covariance matrix for the combined spectrum, and in setting the relative weight of each cross-power spectrum in the final combined spectrum.
4. Apply the procedure presented in Appendix B.1 to obtain an estimate for the point source amplitude. The value obtained is $`\overline{A}=0.0155\pm 0.0017`$, roughly independent of $`l_{\mathrm{max}}`$ in equation (B6). This value is in good agreement with an estimate based on the bispectrum (Komatsu et al., 2003), and on an extrapolation of point source counts (Bennett et al., 2003b). Subtract the point source contribution from the cross-power spectra: $`\overline{C}_l^𝐢=C_l^𝐢\overline{A}S^𝐢w_l^𝐢`$.
5. Compute an approximate form of the full covariance matrix discussed in §4, $`\stackrel{~}{\mathrm{\Sigma }}_{\mathrm{full}}`$. The procedure we use produces a covariance matrix which includes cosmic variance, instrument noise, source subtraction uncertainties, and window function uncertainties. At this stage of the process, it does not yet include the effects of mode coupling. More details are given in §4 and Appendix B.
6. Invert the approximate covariance matrix for use in computing the optimal spectrum. This is the most computationally intensive step in the process.
7. Compute the final combined spectrum from the 28 $`C_l^𝐢`$ as per the procedure given in Appendix C. In particular, assume a fiducial cosmological model (as specified in the Appendix), and use equation (C6), with $`\mathrm{\Sigma }_{\mathrm{full}ll^{}}^{\mathrm{𝐢𝐣}}=\stackrel{~}{\mathrm{\Sigma }}_{\mathrm{full}ll^{}}^{\mathrm{𝐢𝐣}}`$. This produces a final spectrum which is very nearly optimal.

Note: for $`l<100`$ we use a surrogate procedure for computing the combined spectrum. In order to minimize Galactic foreground contamination, we use only V and W band data. Moreover, because the statistics of the $`C_l`$ are mildly non-Gaussian, and because point source contamination and window function uncertainties are small, the “optimal” machinery developed above is unnecessary. Rather we simply form a weighted average spectrum from the V and W band $`C_l^𝐢`$.
8. Compute the approximate inverse-covariance matrix, $`\stackrel{~}{Q}`$, for the combined spectrum using equation (C7). This matrix is approximate in two ways. 1) It does not yet incorporate the effects of mode coupling – this is added below. 2) It has been evaluated for a fiducial cosmological model, $`\stackrel{~}{Q}(C_l^{\mathrm{fid}})`$, while in a likelihood application, we need to evaluate $`\stackrel{~}{Q}(C_l^{\mathrm{th}})`$ for an arbitrary model, $`C_l^{\mathrm{th}}`$ – we add this next.
9. Introduce the dependence on cosmological model into $`\stackrel{~}{Q}`$ as follows. Invert $`\stackrel{~}{Q}`$ to obtain the approximate covariance matrix of the combined spectrum, $`\stackrel{~}{\mathrm{\Sigma }}`$. The off-diagonal terms of $`\stackrel{~}{\mathrm{\Sigma }}`$ are small and weakly dependent on cosmological model, so we expand $`\stackrel{~}{\mathrm{\Sigma }}`$ as

$$\stackrel{~}{\mathrm{\Sigma }}_{ll^{}}=D_l\delta _{ll^{}}+ϵ_{ll^{}},$$
(20)
where $`ϵ_{ll^{}}`$ encodes the mode coupling due to window function and source subtraction uncertainties. This relation defines $`ϵ_{ll^{}}`$, which we take to be zero on the diagonal. $`D_l`$ is dominated by cosmic variance and noise; in order to separate the two contributions, we introduce an ansatz of the form

$$D_l\frac{2}{2l+1}\frac{1}{f_{\mathrm{sky}}}(C_l^{\mathrm{fid}}+N_l^{\mathrm{eff}})^2,$$
(21)
where $`C_l^{\mathrm{fid}}`$ is the fiducial model used to generate the combined spectrum (see Appendix C), and $`N_l^{\mathrm{eff}}`$ is the effective noise in the combined spectrum, which is defined by this equation. The dependence on cosmological model is introduced in the covariance matrix by re-computing $`D_l`$ with $`C_l^{\mathrm{fid}}C_l^{\mathrm{th}}`$, leaving $`N_l^{\mathrm{eff}}`$ fixed.
10. Estimate the coupling induced by the foreground mask as described in Appendix D. The effect of the mask on the off-diagonal elements of the Fisher (inverse-covariance) matrix can be written as

$$(\mathrm{\Sigma }_{\mathrm{mask}}^1)_{ll^{}}r_{ll^{}}(D_l^{}D_l^{}^{})^{1/2},$$
(22)
where $`D_l^{}D_lf_{\mathrm{sky}}^1`$. This expression parameterizes, and defines, the off-diagonal mode coupling in the form of a correlation matrix, $`r_{ll^{}}`$. Note that $`\stackrel{~}{F}_{ll^{}}`$, defined in Appendix D.2, is related to $`r_{ll^{}}`$ by $`\stackrel{~}{F}_{ll^{}}\delta _{ll^{}}+r_{ll^{}}`$.
11. The final Fisher (curvature) matrix, $`Q_{ll^{}}`$, is obtained using

$$Q_{ll^{}}=(D_l^{})^1\delta _{ll^{}}ϵ_{ll^{}}(D_l^{}D_l^{}^{})^{1/2}+r_{ll^{}}(D_l^{}D_l^{}^{})^{1/2},$$
(23)

We further calibrate $`D_l^{}`$ with Monte Carlo simulations. This calibration process, and a description of how the curvature matrix is used in a maximum likelihood determination of cosmological parameters, is given in §2 of Verde et al. (2003). As part of the first-year data release, we provide a Fortran 90 subroutine to evaluate the likelihood of a given cosmological model, $`C_l^{\mathrm{th}}`$, given the WMAP data (supplied in the routine). The code also optionally returns the Fisher (inverse-covariance) matrix for the combined spectrum.

### 5.2 Method II – Combined Sky Map

Our second approach is to form a single co-added map from the Q1 through W4 maps,

$$T=\frac{_{i=3}^{10}T_i^{}/\sigma _{0,i}^2}{_{i=3}^{10}1/\sigma _{0,i}^2},$$
(24)
where $`T_i^{}`$ is the sky map for DA $`i`$ with the best-fit Galactic template model subtracted, and $`\sigma _{0,i}^2`$ is the noise per observation for DA $`i`$, given by Bennett et al. (2003b). We evaluate the power spectrum of this map on the Kp2 cut sky using the quadratic estimator in Appendix A.1. An effective noise model is obtained by using the same approach as described in §2.2: we generate co-added noise maps from the library of end-to-end simulations, evaluate their average spectra, then fit a noise model. The noise bias model is then subtracted from the power spectrum of the combined temperature map. We have performed this analysis with 3 distinct pixel weighting schemes (see Appendix A.1.2) and 3 corresponding noise models. The results are shown in Figure 5 where it is seen that the three cases are virtually indistinguishable.

In effect, this analysis uses both the auto- and cross-power spectra. We view this as a useful check of the more sophisticated procedure described in §5.1, but we do not rely on it for a final result. Uncertainties in the noise model only effect the fourth moment of the cross-power spectra, but they effect the second moment of the auto-power spectra and potentially bias the final result. The $``$6% sensitivity advantage gained by including auto-power spectra was not deemed worth the effort required to guarantee that the final result was not biased.

### 5.3 Comparison of Results

Figure 6 compares the power spectra obtained from methods I and II above. The combined cross-power spectrum from §5.1 is shown in black, the auto-power spectrum obtained in §5.2 from the co-added map is shown in grey. The two methods agree extremely well with the only notable deviation being at the highest $`l`$ range probed by the first-year data. This is the regime where the auto-power spectrum will be most sensitive to the noise bias subtraction. As can be seen in the error estimates shown in Figure 8, the deviation is less than 1$`\sigma `$.

A separate test of robustness is to compute the angular power spectrum in separate regions of the sky to see if the spectrum changes. We have computed the power spectrum in two subsets of the sky – the ecliptic poles, and the ecliptic plane, using the quadratic estimator with the combined sky map. The results are shown in Figure 7 where the pole data is shown in grey and the plane data in black. The two spectra are very consistent overall, but some of the features that appear in the combined spectrum, such as the “peak” at $`l40`$ and the “bite” at $`l210`$, are not robust to this test, thus we consider these features to be of marginal significance. There is also no evidence that beam ellipticity, which would be more manifest in the plane than in the poles, systematically biases the spectrum. This is consistent with estimates of the effect given by Page et al. (2003a).

## 6 DISCUSSION

Our best estimate of the angular power spectrum of the CMB is shown in Figure 8. Also shown is the best-fit $`\mathrm{\Lambda }`$CDM model from Spergel et al. (2003) which is based on a fit to the this spectrum plus a compilation of additional CMB and large-scale structure data. The WMAP data points are plotted with measurement errors based on the diagonal elements of the Fisher matrix presented in Appendix D. The cosmic variance errors, which include the effects of the sky cut, are plotted as a 1$`\sigma `$ band around the best-fit model. As discussed in Spergel et al. (2003), the model is an excellent fit to the data. The combined spectrum provides a definitive measurement of the CMB power spectrum with uncertainties limited by cosmic variance up to $`l350`$. The spectrum clearly exhibits a first acoustic peak at $`l=220.1\pm 0.8`$ and a second acoustic peak at $`l=546\pm 10`$. Page et al. (2003b) present an analysis and interpretation of the peaks and troughs in the first-year WMAP power spectrum.

Figure 9 compares the first-year WMAP spectrum to a compilation of recent balloon and ground-based measurements. In order to make this Figure meaningful, we plot the best-fit model spectrum to represent the WMAP results. The data points are plotted with errors that include both measurement uncertainty and cosmic variance, so no error band is included with the model curve. (Since individual groups report band power measurement with different band-widths, it is not possible to represent a single cosmic variance band that applies to all data sets.) The model spectrum fit to WMAP agrees very well with the ensemble of previous observations.

Wang et al. (2002a) have recently distilled a CMB power spectrum from an optimal combination of the extent pre-WMAP data. In Figure 10 we plot their derived band power points alongside the WMAP data. To make this comparison meaningful, we plot the WMAP data with cosmic variance plus measurement errors and omit the error band from the model spectrum. The distilled spectrum is notably lower than the WMAP data in the vicinity of the first acoustic peak. In a previous version of this work (Wang et al., 2002b) the authors noted that the first peak of their combined spectrum was lower than a significant fraction of their input data. They attribute this to their formalism allowing for a renormalization of individual experiments within their respective calibration uncertainties. Figure 1 in Bennett et al. (2003a) presents a similarly distilled spectrum from the data extent in late 2001 and found a first peak amplitude that was more intuitively consistent with the bulk of the input data, and which is now seen to be consistent with the WMAP power spectrum.

Figure 11 shows the WMAP combined power spectrum compared to the locus of predicted spectra, in red, based on a joint analysis of pre-WMAP CMB data and 2dFGRS large-scale structure data (Percival et al., 2002). As in Figure 8, the WMAP data are plotted with measurement uncertainties, and the best-fit $`\mathrm{\Lambda }`$CDM model (Spergel et al., 2003) is plotted with a 1$`\sigma `$ cosmic variance error band. Percival et al. (2002) predict the location of the first peak should occur at $`l=221.8\pm 2.4`$, which is quite consistent with the value reported by Page et al. (2003b) of $`l=220.1\pm 0.8`$. The height of the first peak was predicted to be in the range $`4920\pm 170`$ $`\mu `$K2, while Page et al. (2003b) report a measured height of $`5580\pm 75`$ $`\mu `$K2, about 13% higher. Unlike the position, the amplitude of the first peak has a complicated dependence on cosmological parameters. Percival et al. (2002) report best-fit parameters for a $`\mathrm{\Lambda }`$CDM model that are mostly consistent with those reported by Spergel et al. (2003) for the same class of models. The only mildly disparate comparison lies in the combination of normalization, $`\sigma _8`$, and optical depth, $`\tau `$. Percival et al. (2002) report the product $`\sigma _8e^\tau =0.72\pm 0.03\pm 0.02`$, where the first error is a “theory” error and the second is measurement error. While Spergel et al. (2003) does not report a maximum likelihood range for this explicit parameter combination, the product of their maximum likelihood values for $`\sigma _8`$ and $`\tau `$ yields $`\sigma _8e^\tau =0.74`$, which is consistent with Percival et al. (2002), but would make the first peak a few percent higher. Small differences in $`\mathrm{\Omega }_bh^2`$, $`\mathrm{\Omega }_mh^2`$, and $`n_s`$, may also contribute to the difference.

## 7 CONCLUSIONS

We present measurements of the angular power spectrum of the cosmic microwave background from the first-year WMAP data. The eight high-frequency sky maps from DAs Q1 through W4 were used to estimate 28 cross-power spectra, which are largely independent of the noise properties of the experiment. These data were tested for consistency in §3, then used in §5 as input to a final combined spectrum, discussed in §6. The procedure for estimating the uncertainties in the final combined spectrum were discussed in §4 and in numerous Appendices.

The combined spectrum provides a definitive measurement of the CMB power spectrum, with uncertainties limited by cosmic variance up to $`l350`$, and a signal to noise per mode $`>1`$ up to $`l650`$. The spectrum clearly exhibits a first acoustic peak at $`l=220.1\pm 0.8`$ and a second acoustic peak at $`l=546\pm 10`$. Page et al. (2003b) present an analysis and interpretation of the peaks and troughs in the first-year WMAP power spectrum. Spergel et al. (2003), Verde et al. (2003), and Peiris et al. (2003) analyze the combined spectrum in the context of cosmological models. They conclude that the data provide strong support adiabatic initial conditions, and they give precise measurements of a number of cosmological parameters. Kogut et al. (2003) analyze the correlation between WMAP’s temperature and polarization signals, the $`C_l^{TE}`$ spectrum, and present evidence for a relatively high optical depth, and an early period of cosmic reionization. Among other things, this result implies that the temperature power spectrum is suppressed by $``$30% on degree angular scales, due to secondary scattering.

A variety of first-year WMAP data products are being made available by NASA’s new Legacy Archive for Microwave Background Data Analysis (LAMBDA). In addition to the sky maps and calibrated time-ordered data, we are providing the 28 cross power spectra used in this paper (with diffuse foregrounds subtracted), the combined spectrum from §5.1, and a Fortran 90 subroutine to compute the likelihood of a given cosmological model, (the code will also optionally return the Fisher (inverse-covariance) matrix for the combined spectrum.) The LAMBDA URL is `http://lambda.gsfc.nasa.gov/`.

The WMAP mission is made possible by the support of the Office of Space Sciences at NASA Headquarters and by the hard and capable work of scores of scientists, engineers, technicians, machinists, data analysts, budget analysts, managers, administrative staff, and reviewers. We thank Mike Nolta for helpful comments on an earlier draft of this manuscript. LV is supported by NASA through Chandra Fellowship PF2-30022 issued by the Chandra X-ray Observatory center, which is operated by the Smithsonian Astrophysical Observatory for an on behalf of NASA under contract NAS8-39073. We acknowledge use of the HEALPix package.

## Appendix A POWER SPECTRUM ESTIMATION METHODS

For the analysis of WMAP’s first-year data, we have chosen two distinct methods for inferring the power spectrum. The first is a fast and accurate quadratic method for estimating the power spectrum of a partial sky map (Hivon et al., 2002). We summarize the basic approach here, highlighting the aspects of the method that are especially pertinent to WMAP, and refer the reader to Hivon et al. (2002) for details. The second is a maximum likelihood method that provides an independent estimate of the spectrum measured by WMAP (Oh et al., 1999).

We have applied both of these methods to the WMAP data. The results are shown in Figure 12, which shows spectra estimated from the V band map, up to $`l=200`$, for the two methods. The maximum likelihood estimate has slightly lower uncertainties at low $`l`$, because the method optimally weights the data with a pixel-pixel covariance $`𝐂=𝐒+𝐍𝐒`$ where $`𝐒`$ is the covariance of the CMB signal and $`𝐍`$ is the covariance of the noise (see Appendix A.2). Our quadratic estimator uses uniform pixel weights at low $`l`$ (see Appendix A.1.2) though it is clear from the Figure that the difference is not significant. At high $`l`$, where the WMAP data are noise dominated, the two estimators give essentially identical results because they effectively weight the data in the same way.

To obtain our “best” estimate of the WMAP power spectrum, we adopt the quadratic estimator because it can be readily applied to pairs of WMAP radiometers in a way that is nearly independent of the properties of the instrument noise. In §5 we discuss our methodology for combining spectra from pairs of radiometers and present the final combined spectrum.

### A.1 Quadratic Estimation

Hivon et al. (2002) start with the full-sky estimator in equation (2), add a position dependent weight, $`W(𝐧)`$, and set $`W`$ to zero in the regions where the sky is contaminated. In other regions, $`W`$ is chosen to optimize the sensitivity of the estimator. A temperature map $`\mathrm{\Delta }T(𝐧)`$ on which a weight $`W(𝐧)`$ is applied can be decomposed in spherical harmonic coefficients as

$`\stackrel{~}{a}_{lm}`$ $`=`$ $`{\displaystyle 𝑑\mathrm{\Omega }_𝐧\mathrm{\Delta }T(𝐧)W(𝐧)Y_{lm}^{}(𝐧)}`$ (A1)
$``$ $`\mathrm{\Omega }_p{\displaystyle \underset{p}{}}\mathrm{\Delta }T(p)W(p)Y_{lm}^{}(p),`$ (A2)
where the integral over the sky is approximated by a discrete sum over map pixels, each of which subtends solid angle $`\mathrm{\Omega }_p`$. Hivon et al. (2002) then define the “pseudo power spectrum” $`\stackrel{~}{C}_l`$ as

$$\stackrel{~}{C}_l=\frac{1}{2l+1}\underset{m=l}{\overset{l}{}}\left|\stackrel{~}{a}_{lm}\right|^2.$$
(A3)
The pseudo power spectrum $`\stackrel{~}{C}_l`$, given by the weighted spherical harmonic transform of a map, is clearly different from the full sky angular spectrum, $`C_l^{\mathrm{sky}}`$, but the ensemble averages of the two spectra can be related by

$$\stackrel{~}{C}_l=\underset{l^{}}{}G_{ll^{}}C_l^{}^{\mathrm{sky}}$$
(A4)
where $`G_{ll^{}}`$ describes the mode coupling resulting from the weight function $`W(𝐧)`$ (Hauser & Peebles, 1973). Hivon et al. (2002) give the following expression for the coupling matrix, which depends only on the geometry of the weight function $`W(𝐧)`$

$$G_{l_1l_2}=\frac{2l_2+1}{4\pi }\underset{l_3}{}(2l_3+1)𝒲_{l_3}\left(\begin{array}{ccc}l_1\hfill & l_2& \hfill l_3\\ 0\hfill & 0& \hfill 0\end{array}\right)^2,$$
(A5)
where the final term in parentheses is the Wigner 3-$`j`$ symbol, and $`𝒲_l`$ is the angular power spectrum of the weight function

$$𝒲_l=\frac{1}{2l+1}\underset{m}{}\left|w_{lm}\right|^2,$$
(A6)
where

$$w_{lm}=𝑑\mathrm{\Omega }_𝐧W(𝐧)Y_{lm}^{}(𝐧).$$
(A7)
Upon inverting the coupling matrix $`G_{ll^{}}`$ and making the identification $`C_l^{\mathrm{sky}}=C_l`$, we obtain the following estimator of the power spectrum

$$C_l^{\mathrm{obs}}=\underset{l^{}}{}G_{ll^{}}^1\stackrel{~}{C}_l^{}.$$
(A8)
The computation of equation (A2) for each $`(l,m)`$ up to $`l=l_{\mathrm{max}}`$ would scale as $`N_{\mathrm{pix}}l_{\mathrm{max}}^2`$ if performed on an arbitrary pixelization of the sphere, where $`N_{\mathrm{pix}}`$ is the number of sky map pixels. However, for a pixelization scheme with iso-latitude pixel centers, fast FFT methods may be employed to speed up the evaluation, so it scales like $`N_{\mathrm{pix}}^{1/2}l_{\mathrm{max}}^2`$ (Muciaccia et al., 1997). The WMAP sky maps have been produced using the HEALPix layout (Górski et al., 1998) which supports such fast spherical harmonic transforms. In particular, the HEALPix routine `map2alm` evaluates equation (A2).

#### A.1.1 Auto- and Cross-Power Spectra from the WMAP Data

For a multi-channel experiment like WMAP it is quite powerful to evaluate the power spectra from different maps and compare results. In particular, the quadratic estimator described above may be used on 1 or 2 maps at a time by generalizing the expression for the pseudo power spectrum equation (A3) as

$$\stackrel{~}{C}_l=\frac{1}{2l+1}\underset{m=l}{\overset{l}{}}\stackrel{~}{a}_{lm}^i\stackrel{~}{a}_{lm}^j$$
(A9)
where $`\stackrel{~}{a}_{lm}^i`$ refers to the transform of map $`i`$ and $`\stackrel{~}{a}_{lm}^j`$ refers to map $`j`$, which needn’t be the same as map $`i`$. As discussed in §2, if $`ij`$ and the noise in the two maps is uncorrelated, the estimator equation (A9) provides an unbiased estimate of the underlying power spectrum.

We have tested the auto- and cross-power estimator extensively with Monte Carlo simulations of the first-year WMAP data. Selected results from this testing are shown in Figure 13. The auto- and cross-power spectra that obtain from the flight data are presented in detail in §3. The cross-power spectra form the basis for our final combined spectrum, presented in §5.

#### A.1.2 Choice of Weighting

We seek a weighting scheme that mimics the maximum likelihood estimation (Appendix A.2), which effectively weights the data by the full inverse covariance matrix $`𝐂^1=(𝐒+𝐍)^1`$. For the combined spectrum presented in §5, we use three distinct weighting functions in three separate $`l`$ ranges.

1. For $`l<200`$ we give equal weight to all un-cut pixels,

$$W(p)=M(p)$$
(A10)
where $`M(p)`$ is the Kp2 sky mask, defined by Bennett et al. (2003c). It takes values of 0 within the mask and 1 otherwise.
2. For $`l>450`$ we use inverse-noise weighting,

$$W(p)=M(p)N_{\mathrm{obs}}(p)$$
(A11)
where $`N_{\mathrm{obs}}(p)`$ is the number of observations of pixel $`p`$.
3. For $`200<l<450`$ we use a transitional weighting,

$$W(p)=\frac{M(p)}{1/N_{\mathrm{obs}}+1/N_{\mathrm{obs}}(p)}$$
(A12)
where $`N_{\mathrm{obs}}`$ is the mean of $`N_{\mathrm{obs}}`$ evaluated over the cut sky.

Verde et al. (2003) discuss the choice of weighting in more detail.

### A.2 Maximum Likelihood Estimation

This Appendix provides a summary of the maximum likelihood approach to power spectrum estimation originally presented by Oh et al. (1999). If the temperature fluctuations are Gaussian, and the a priori probability of a given set of cosmological parameters is uniform, then the power spectrum may be estimated by maximizing the multi-variate Gaussian likelihood function

$$(C_l|𝐦)=\frac{\mathrm{exp}(\frac{1}{2}𝐦^T𝐂^1𝐦)}{\sqrt{det𝐂}}$$
(A13)
where $`𝐦`$ is a data vector (see below) and $`𝐂`$ is the covariance matrix of the data, which has contributions from both the signal and the instrument noise, $`𝐂=𝐒+𝐍`$. We can work in whatever basis is most convenient. In the pixel basis the data are the sky map pixel temperatures, and in the spherical harmonic basis the data are the $`a_{lm}`$ coefficients of the map. In the former basis the noise covariance is nearly diagonal, while in the latter, the signal covariance is.

$$\begin{array}{ccccccc}& & \text{Pixel basis:}\hfill & & & & \text{Spherical harmonic basis:}\hfill \\ 𝐦\hfill & \hfill & T_i\hfill & & 𝐦\hfill & \hfill & a_{lm}\hfill \\ 𝐒\hfill & \hfill & _l\frac{(2l+1)}{4\pi }C_lP_l(\mathrm{cos}\theta _{ij})\hfill & & 𝐒\hfill & \hfill & \mathrm{diag}(C_2,C_2,\mathrm{}C_3,C_3,\mathrm{})\hfill \\ 𝐍\hfill & \hfill & \sigma _i^2\delta _{ij}\hfill & & 𝐍\hfill & \hfill & N_{(lm)(lm)^{}}\text{(see below)}.\hfill \end{array}$$

For WMAP, the length of the data vector, $`N_{\mathrm{data}}`$, is 2,672,361, the number of $`7^{}`$ sky map pixels (HEALPix $`N_{\mathrm{side}}=512`$) that survive the Galaxy cut, so it is necessary to find methods for evaluating $``$ that do not require a full inversion of the covariance matrix $`𝐂`$, which requires $`O\left(N_{\mathrm{data}}^3\right)`$ operations. We use an iterative method for evaluating the likelihood that exploits the ability to find an approximate inverse $`\stackrel{~}{𝐂}^1`$. The most important features in the data that make this possible are that WMAP observes the full sky and the Galaxy cut is predominantly azimuthally symmetric in Galactic coordinates (Bennett et al., 2003b). Of secondary importance for this pre-conditioner is that WMAP’s noise per pixel does not vary strongly across the sky (Bennett et al., 2003a). We discuss the pre-conditioner in more detail below.

Defining $`f2\mathrm{ln}`$ and $`𝐏_l\frac{𝐂}{C_l}`$, we maximize the likelihood by solving

$$\frac{f}{C_l}=0=𝐦^T𝐂^1𝐏_l𝐂^1𝐦+\mathrm{tr}(𝐂^1𝐏_l)$$
(A14)
using a Newton-Raphson root finding method that generates an iterative estimate of the angular power spectrum

$$C_l^{(n+1)}=C_l^{(n)}\frac{1}{2}\underset{l^{}}{}F_{ll^{}}\frac{f}{C_l^{}}$$
(A15)
at each step. Here $`F_{ll^{}}`$ is the Fisher matrix

$$F_{ll^{}}=^2\mathrm{ln}/C_lC_l^{}=\frac{1}{2}\mathrm{tr}(𝐂^1𝐏_l𝐂^1𝐏_l^{}).$$
(A16)
To implement the solution in equation (A15) we need a fast way to evaluate the following three components of $`_l^{}F_{ll^{}}\frac{f}{C_l^{}}`$:

1. $`𝐦^T𝐂^1𝐏_l𝐂^1𝐦`$
2. $`\mathrm{tr}(𝐂^1𝐏_l)`$
3. $`\mathrm{tr}(𝐂^1𝐏_l𝐂^1𝐏_l^{})`$.

We use the spherical harmonic basis in which the data vector consists of the $`a_{lm}`$ coefficients of the map obtained by least squares fitting on the cut sky. The signal covariance is diagonal in this basis, while the noise matrix is obtained from the normal equations for the $`a_{lm}`$ fit

$$\underset{(lm)^{}}{}N_{(lm)(lm)^{}}^1a_{(lm)^{}}=y_{(lm)}$$
(A17)
where

$`N_{(lm)(lm)^{}}^1`$ $``$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{Y_{(lm)}(\widehat{n}_i)Y_{(lm)^{}}(\widehat{n}_i)}{\sigma _i^2}}`$ (A18)
$`y_{(lm)}`$ $``$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{T_iY_{(lm)}(\widehat{n}_i)}{\sigma _i^2}}.`$ (A19)
The sums are over all sky map pixels that survive the Galaxy cut, and we have assumed that the noise is uncorrelated from pixel to pixel.

#### A.2.1 Evaluation of $`𝐂^1𝐦`$

The term $`𝐂^1𝐦`$ appears repeatedly in the evaluation of equation (A15). We compute this by solving $`𝐂𝐳=(𝐒+𝐍)𝐳=𝐦`$ for $`𝐳`$. A more numerically tractable system is obtained by multiplying both sides by $`𝐒^{\frac{1}{2}}𝐍^1`$, so

$$(𝐈+𝐒^{\frac{1}{2}}𝐍^1𝐒^{\frac{1}{2}})𝐒^{\frac{1}{2}}𝐳=𝐒^{\frac{1}{2}}𝐍^1𝐦=𝐒^{\frac{1}{2}}𝐲$$
(A20)
where $`𝐲`$ is the spherical harmonic transform of the map, defined in equation (A19). Note that $`𝐲`$ can be quickly computed in any pixelization scheme that has iso-latitude pixel centers with fixed longitude spacing. We then solve equation (A20) using an iterative conjugate gradient method with a pre-conditioner for the matrix $`𝐀(𝐈+𝐒^{\frac{1}{2}}𝐍^1𝐒^{\frac{1}{2}})`$. We find the following block-diagonal form of $`𝐀`$ to be a good starting point

$$\stackrel{~}{𝐀}=\left(\begin{array}{cc}𝐈+𝐒^{\frac{1}{2}}\stackrel{~}{𝐍}^1𝐒^{\frac{1}{2}}& 0\\ 0& \text{diag}(𝐈+𝐒^{\frac{1}{2}}\stackrel{~}{𝐍}^1𝐒^{\frac{1}{2}})\end{array}\right)$$
(A21)
where $`\stackrel{~}{𝐍}^1`$ is a block-diagonal approximation of the noise matrix discussed below. The lower-right block of $`\stackrel{~}{𝐀}`$ occupies the high $`l`$ portion of the matrix where the signal to noise ratio $`𝐒^{\frac{1}{2}}𝐍^1𝐒^{\frac{1}{2}}`$ is low, so a diagonal approximation is adequate. The upper-left block occupies the low $`l`$ portion of the matrix where the signal dominates the noise, so we need a better estimate of $`𝐍^1`$. In practice we find this split works well at $`l=512`$ for the WMAP noise levels. As for the approximate form of $`𝐍^1`$, defined in equation (A19), note that the dominant off-diagonal contributions arise from the azimuthally symmetric Galaxy cut, which couples different $`l`$ modes, but not $`m`$ modes. Thus $`𝐍^1`$ is approximately block diagonal, with perturbations induced by the non-uniform sky coverage of WMAP. We therefore use a block diagonal form of $`𝐍^1`$ as the pre-conditioner,

$$\stackrel{~}{N}_{(lm)(lm)^{}}^1=N_{(lm)(lm)^{}}^1\delta _{mm^{}}.$$
(A22)
Using the pre-conditioner equation (A21) we find that the conjugate gradient solution of equation (A20) converges in approximately six iterations and requires only cpu-minutes of processing on an SGI Origin 2000.

## Appendix B POINT SOURCE SUBTRACTION

In this Appendix we describe the procedure we use to estimate and subtract the point source contribution directly from the multi-frequency cross-power spectra. We then show how we incorporate the source model uncertainty into the covariance matrix of the source-subtracted spectra by marginalizing a Gaussian likelihood function over the source model amplitude parameter.

This marks the starting point of the multi-frequency analysis which will lead to the combined power spectrum, discussed in §5.1. In order to generate the combined spectrum, we need the full covariance matrix of the cross-power spectra (see §4). Our approach to generating the full covariance is to start with the ideal, full-sky form, which only includes cosmic variance and instrument noise, then we incorporate additional effects step by step, as outlined in §5.1 and in these Appendices. For an ideal experiment with full sky coverage, no point source contamination, and no beam uncertainty the covariance matrix is

$$\mathrm{\Sigma }_{ll}^{ii^{}jj^{}}=\frac{1}{(2l+1)}\left[\left(C_l^{\mathrm{th}}w_l^{ij}+n^in^j\delta _{ij}\right)\left(C_l^{\mathrm{th}}w_l^{i^{}j^{}}+n^i^{}n^j^{}\delta _{i^{}j^{}}\right)+\left(C_l^{\mathrm{th}}w_l^{ij^{}}+n^in^j^{}\delta _{ij^{}}\right)\left(C_l^{\mathrm{th}}w_l^{i^{}j}+n^i^{}n^j\delta _{i^{}j}\right)\right]$$
(B1)
where $`\delta _{ij}`$ denotes the Kronecker symbol, $`w_l^{ij}b_l^ib_l^jp_l^2`$ is the window function, and $`n^in^i^{}\delta _{ii^{}}=N^𝐢`$.

### B.1 Estimating the Point Source Amplitude

We start by assuming a Gaussian likelihood for the sky model, given the WMAP data

$$2\mathrm{ln}(A,C_l|C_l^𝐢)=\underset{\mathrm{𝐢𝐣}l}{}[C_l^𝐢(C_l+AS^𝐢)w_l^𝐢](\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}[C_l^𝐣(C_l+AS^𝐣)w_l^𝐣],$$
(B2)
where $`C_l^𝐢`$ is cross-power spectrum $`𝐢`$, $`w_l^𝐢`$ is the window function for spectrum $`𝐢`$, $`C_l`$ is the true CMB power spectrum, and $`AS^𝐢`$ is the source model defined in equations (9) and (11). Here we assume the diagonal form of the covariance matrix in equation (B1).

To determine the best-fit source amplitude, $`A`$, we marginalize this likelihood over the CMB spectrum, $`C_l`$. First we expand equation (B2) as

$`2\mathrm{ln}`$ $`=`$ $`{\displaystyle \underset{l}{}}{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}(C_l^𝐢AS^𝐢w_l^𝐢)(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}(C_l^𝐣AS^𝐣w_l^𝐣)`$ (B3)
$``$ $`{\displaystyle \underset{l}{}}2C_l{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}(C_l^𝐢AS^𝐢w_l^𝐢)(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}w_l^𝐣`$
$`+`$ $`{\displaystyle \underset{l}{}}C_l^2{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}w_l^𝐢(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}w_l^𝐣,`$
which is a quadratic form in $`C_l`$, $`_l\mathrm{exp}[\frac{1}{2}(aC_l^22bC_l+c)]`$, with

$`a`$ $``$ $`{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}w_l^𝐢(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}w_l^𝐣,`$
$`b`$ $``$ $`{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}(C_l^𝐢AS^𝐢w_l^𝐢)(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}w_l^𝐣,`$
$`c`$ $``$ $`{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}(C_l^𝐢AS^𝐢w_l^𝐢)(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}(C_l^𝐣AS^𝐣w_l^𝐣).`$ (B4)
We wish to evaluate the marginalized likelihood, $`_{C_l}(A)=𝑑C_l(A,C_l)`$. This is readily evaluated using the substitution $`C_lC_lb/a`$, giving $`_{C_l}_l\mathrm{exp}[\frac{1}{2}(cb^2/a)]`$, where we drop a multiplicative term proportional to $`a`$, which is independent of $`A`$. The marginalized likelihood function is thus

$`2\mathrm{ln}_{C_l}`$ $`=`$ $`{\displaystyle \underset{l}{}}{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}(C_l^𝐢AS^𝐢w_l^𝐢)(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}(C_l^𝐣AS^𝐣w_l^𝐣)`$ (B5)
$``$ $`{\displaystyle \underset{l}{}}{\displaystyle \frac{1}{_{𝐢^{}𝐣^{}}w_l^𝐢^{}(\mathrm{\Sigma }^1)_l^{𝐢^{}𝐣^{}}w_l^𝐣^{}}}\left[{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}(C_l^𝐢AS^𝐢w_l^𝐢)(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}w_l^𝐣\right]^2.`$
Setting $`_{C_l}(A)/A=0`$ gives the most likely point sources amplitude

$$\overline{A}=\frac{_{\mathrm{𝐢𝐣}l}C_l^𝐣(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}H_l^𝐢w_l^𝐢}{_{\mathrm{𝐢𝐣}l}S^𝐣w_l^𝐣(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}H_l^𝐢w_l^𝐢},$$
(B6)
where

$$H_l^𝐢=S^𝐢\frac{_{𝐢^{}𝐣^{}}S^𝐢^{}(\mathrm{\Sigma }^1)_l^{𝐢^{}𝐣^{}}w_l^𝐣^{}}{_{𝐢^{}𝐣^{}}w_l^𝐢^{}(\mathrm{\Sigma }^1)_l^{𝐢^{}𝐣^{}}w_l^𝐣^{}}.$$
(B7)
The standard error on the best fit value of $`A`$ is

$$\sigma _{\mathrm{src}}^2=\left[\underset{\mathrm{𝐢𝐣}l}{}H_l^𝐢w_l^𝐢(\mathrm{\Sigma }^1)_l^{\mathrm{𝐢𝐣}}S_l^𝐣w_l^𝐣\right]^1.$$
(B8)
For WMAP, the off-diagonal terms in the covariance matrix are small. Here, neglecting them changes the inferred point source amplitude by less than 1%.

### B.2 Marginalizing over Point Source Amplitude

The source subtraction procedure discussed above is uncertain. In this section we incorporate this uncertainty into the full covariance matrix for the cross-power spectra. We again assume a Gaussian likelihood function of the form

$$2\mathrm{ln}=\underset{\mathrm{𝐢𝐣}ll^{}}{}[\overline{C}_l^𝐢(C_l^{\mathrm{th}}+\delta AS_l^𝐢)w_l^𝐢](\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}[\overline{C}_l^{}^𝐣(C_l^{}^{\mathrm{th}}+\delta AS_l^{}^𝐣)w_l^{}^𝐣]$$
(B9)
where $`\overline{C}_l^𝐢`$ is the source-subtracted cross-power spectrum for DA pair $`𝐢`$, obtained above, $`w_l^𝐢`$ is the window function for spectrum $`𝐢`$, $`C_l^{\mathrm{th}}`$ is now a fixed CMB model spectrum, and $`\delta A`$ is the residual source amplitude.

We can marginalize the likelihood function over the residual point source amplitude as follows. Expand equation (B9) as

$`2\mathrm{ln}`$ $`=`$ $`{\displaystyle \underset{\mathrm{𝐢𝐣}ll^{}}{}}(\overline{C}_l^𝐢C_l^{\mathrm{th}}w_l^𝐢)(\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}(\overline{C}_l^{}^𝐣C_l^{}^{\mathrm{th}}w_l^{}^𝐣)`$ (B10)
$``$ $`2(\delta A){\displaystyle \underset{\mathrm{𝐢𝐣}ll^{}}{}}(\overline{C}_l^𝐢C_l^{\mathrm{th}}w_l^𝐢)(\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}S^𝐣w_l^{}^𝐣`$
$`+`$ $`(\delta A)^2{\displaystyle \underset{\mathrm{𝐢𝐣}ll^{}}{}}S^𝐢w_l^𝐢(\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}S^𝐣w_l^{}^𝐣,`$
which is a quadratic form in $`\delta A`$, $`\mathrm{exp}[\frac{1}{2}(a(\delta A)^22b(\delta A)+c)]`$, with

$`a`$ $``$ $`{\displaystyle \underset{\mathrm{𝐢𝐣}ll^{}}{}}S^𝐢w_l^𝐢(\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}S^𝐣w_l^{}^𝐣,`$
$`b`$ $``$ $`{\displaystyle \underset{\mathrm{𝐢𝐣}ll^{}}{}}(\overline{C}_l^𝐢C_l^{\mathrm{th}}w_l^𝐢)(\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}S^𝐣w_l^{}^𝐣,`$
$`c`$ $``$ $`{\displaystyle \underset{\mathrm{𝐢𝐣}ll^{}}{}}(\overline{C}_l^𝐢C_l^{\mathrm{th}}w_l^𝐢)(\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}(\overline{C}_l^{}^𝐣C_l^{}^{\mathrm{th}}w_l^{}^𝐣).`$ (B11)
We wish to evaluate the marginalized likelihood, $`_A=d(\delta A)`$. This is readily evaluated using the substitution $`\delta A\delta Ab/a`$, giving $`_A\mathrm{exp}[\frac{1}{2}(cb^2/a)]`$, where we drop a multiplicative term proportional to $`a`$, which is independent of $`\overline{C}_l^𝐢`$. The marginalized likelihood function is thus

$`2\mathrm{ln}_A`$ $`=`$ $`{\displaystyle \underset{\mathrm{𝐢𝐣}ll^{}}{}}(\overline{C}_l^𝐢C_l^{\mathrm{th}}w_l^𝐢)(\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}(\overline{C}_l^{}^𝐣C_l^{}^{\mathrm{th}}w_l^{}^𝐣)`$ (B12)
$``$ $`{\displaystyle \frac{1}{a}}{\displaystyle \underset{\mathrm{𝐢𝐣}ll^{}}{}}(\overline{C}_l^𝐢C_l^{\mathrm{th}}w_l^𝐢)(\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}S^𝐣w_l^{}^𝐣{\displaystyle \underset{𝐢^{}𝐣^{}l^{\prime \prime }l^{\prime \prime \prime }}{}}(\overline{C}_{l^{\prime \prime }}^𝐢^{}C_{l^{\prime \prime }}^{\mathrm{th}}w_{l^{\prime \prime }}^𝐢^{})(\mathrm{\Sigma }^1)_{l^{\prime \prime }l^{\prime \prime \prime }}^{𝐢^{}𝐣^{}}S^𝐣^{}w_{l^{\prime \prime \prime }}^𝐣^{}.`$ (B13)
This expression may be recast in the form

$$2\mathrm{ln}_A=\underset{\mathrm{𝐢𝐣}ll^{}}{}(\overline{C}_l^𝐢C_l^{\mathrm{th}}w_l^𝐢)(F^{\mathrm{src}})_{ll^{}}^{\mathrm{𝐢𝐣}}(\overline{C}_l^{}^𝐣C_l^{}^{\mathrm{th}}w_l^{}^𝐣),$$
(B14)
where $`(F^{\mathrm{src}})_{ll^{}}^{\mathrm{𝐢𝐣}}`$ is

$$(F^{\mathrm{src}})_{ll^{}}^{\mathrm{𝐢𝐣}}(\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}\frac{1}{a}B_{ll^{}}^{\mathrm{𝐢𝐣}},$$
(B15)
with

$$B_{ll^{}}^{\mathrm{𝐢𝐣}}\underset{𝐢^{}𝐣^{}l^{\prime \prime }l^{\prime \prime \prime }}{}(\mathrm{\Sigma }^1)_{ll^{\prime \prime }}^{\mathrm{𝐢𝐢}^{}}S^𝐢^{}w_{l^{\prime \prime }}^𝐢^{}(\mathrm{\Sigma }^1)_{l^{}l^{\prime \prime \prime }}^{\mathrm{𝐣𝐣}^{}}S^𝐣^{}w_{l^{\prime \prime \prime }}^𝐣^{}.$$
(B16)
The superscript “src” indicates that the Fisher matrix so obtained includes point source subtraction uncertainty, in addition to whatever effects have been included in $`(\mathrm{\Sigma }^1)_{ll^{}}^{\mathrm{𝐢𝐣}}`$ already, in this case only cosmic variance and noise. Note that equation (B15) neglects a term proportional to $`\mathrm{ln}deta`$ which has a weak dependence on cosmological parameters. In the actual calculation, as in the previous section, we assume the diagonal form of the covariance matrix in equation (B1).

## Appendix C OPTIMAL WEIGHTING OF MULTI-CHANNEL SPECTRA

We use the 8 high-frequency differencing assemblies Q1 through W4 to generate the final combined spectrum. In this Appendix we show how we combine these spectra into a single estimate of the angular power spectrum.

The ultimate goal of the WMAP analysis is to produce a likelihood function for a set of cosmological parameters, $`\{\stackrel{}{\alpha }\}`$, given the data, $`C_l^𝐢`$. Specifically

$$𝒫(\stackrel{}{\alpha }|C_l^𝐢)=(C_l^𝐢|C_l^{\mathrm{th}}(\stackrel{}{\alpha }))𝒫(\stackrel{}{\alpha }),$$
(C1)
where $`𝒫(\stackrel{}{\alpha }|C_l^𝐢)`$ is the probability of $`\{\stackrel{}{\alpha }\}`$ given the data, $`(C_l^𝐢|C_l^{\mathrm{th}}(\stackrel{}{\alpha }))`$ is the likelihood of the data given the model, $`C_l^{\mathrm{th}}(\stackrel{}{\alpha })`$, and $`𝒫(\stackrel{}{\alpha })`$ is the prior probability of the parameter set (Verde et al., 2003). To this end, we seek a combined spectrum, $`\widehat{C}_l`$, that estimates the power spectrum in our sky, $`C_l^{\mathrm{sky}}`$, with the property that $`𝒫(\stackrel{}{\alpha }|\widehat{C}_l)=𝒫(\stackrel{}{\alpha }|C_l^𝐢)`$, and hence $`(\widehat{C}_l|C_l^{\mathrm{th}}(\stackrel{}{\alpha }))=(C_l^𝐢|C_l^{\mathrm{th}}(\stackrel{}{\alpha }))`$.

To estimate the combined spectrum, we approximate the likelihood function for the cross-power spectra as Gaussian

$$2\mathrm{ln}(\overline{C}_l^𝐢|C_l^{\mathrm{th}})=\underset{\mathrm{𝐢𝐣}ll^{}}{}(\overline{C}_l^𝐢C_l^{\mathrm{th}}w_l^𝐢)(\mathrm{\Sigma }_{\mathrm{full}}^1)_{ll^{}}^{\mathrm{𝐢𝐣}}(\overline{C}_l^{}^𝐣C_l^{}^{\mathrm{th}}w_l^{}^𝐣),$$
(C2)
where $`\overline{C}_l^𝐢`$ is the spectrum with the best-fit source model subtracted, defined after equation (12), $`w_l^𝐢`$ is the window function of spectrum $`𝐢`$, defined after equation (6), and $`\mathrm{\Sigma }_{\mathrm{full}}`$ is the covariance matrix of the 28 cross-power spectra. Note that the treatment in the remainder of this section does not depend on any specific property of the covariance, so we use generic notation for readability. However, when we generate the WMAP first-year combined spectrum, the actual form of the covariance used at this step is $`(\stackrel{~}{\mathrm{\Sigma }}_{\mathrm{full}})_{ll^{}}^{\mathrm{𝐢𝐣}}`$, where $`\stackrel{~}{\mathrm{\Sigma }}`$ indicates the approximate covariance, obtained in §4, that has not yet had the effects of the foreground mask accounted for.

We seek a spectrum $`\widehat{C}_l`$ such that

$$2\mathrm{ln}(\widehat{C}_l|C_l^{\mathrm{th}})\underset{ll^{}}{}(\widehat{C}_lC_l^{\mathrm{th}})Q_{ll^{}}(\widehat{C}_l^{}C_l^{}^{\mathrm{th}})=2\mathrm{ln}(\overline{C}_l^𝐢|C_l^{\mathrm{th}}),$$
(C3)
where $`Q_{ll^{}}`$ is the inverse-covariance matrix of the combined spectrum which comes with the estimate of $`\widehat{C}_l`$. Suppose, for simplicity, that $`\mathrm{\Sigma }_{\mathrm{full}}`$ is diagonal, $`(\mathrm{\Sigma }_{\mathrm{full}})_{ll^{}}^{\mathrm{𝐢𝐣}}=(\mathrm{\Sigma }_{\mathrm{full}})_l^{\mathrm{𝐢𝐣}}\delta _{ll^{}}`$, then it is straightforward to show that the deconvolved, weighted-average spectrum

$$\widehat{C}_l=\frac{_{\mathrm{𝐢𝐣}}\overline{C}_l^𝐢(\mathrm{\Sigma }_{\mathrm{full}}^1)_l^{\mathrm{𝐢𝐣}}w_l^𝐣}{Q_l},$$
(C4)
with

$$Q_l=\underset{\mathrm{𝐢𝐣}}{}w_l^𝐢(\mathrm{\Sigma }_{\mathrm{full}}^1)_l^{\mathrm{𝐢𝐣}}w_l^𝐣,$$
(C5)
is the desired spectrum. Substituting equations (C4) and (C5) into equation (C3) produces equation (C2) up to a term which has a weak dependence on comological model, which we ignore. This combined spectrum is equivalent to the result we would obtain using the “optimal data compression” of Tegmark et al. (1997).

If the inverse covariance matrix is not diagonal in $`l`$, it can be shown that the optimal combined spectrum is given by

$$\widehat{C}_l=\frac{_{\mathrm{𝐢𝐣}l^{}}\overline{C}_l^𝐢(\mathrm{\Sigma }_{\mathrm{full}}^1)_{ll^{}}^{\mathrm{𝐢𝐣}}C_l^{}^{\mathrm{fid}}w_l^{}^𝐣}{_l^{}Q_{ll^{}}C_l^{}^{\mathrm{fid}}},$$
(C6)
with

$$Q_{ll^{}}=\underset{\mathrm{𝐢𝐣}}{}w_l^𝐢(\mathrm{\Sigma }_{\mathrm{full}}^1)_{ll^{}}^{\mathrm{𝐢𝐣}}w_l^{}^𝐣,$$
(C7)
and where $`C_l^{\mathrm{fid}}`$ is a fiducial cosmological model which we take to be a flat $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_bh^2=0.021`$, $`\mathrm{\Omega }_ch^2=0.129`$, $`\mathrm{\Omega }_{\mathrm{tot}}=1`$, $`h=0.68`$, $`n_s=1.2`$, and $`\tau =0.2`$. While this model has parameters values that are substantially different than the best-fit WMAP model obtained (afterwards) by Spergel et al. (2003), the parameter degeneracies are such that $`C_l^{\mathrm{fid}}`$ is close to the best-fit model for $`l>100`$ where this estimator is actually used (see §5.1). The combined spectrum is optimal if the fiducial model chosen is the correct one; otherwise it is still unbiased but slightly sub-optimal (Gupta & Heavens, 2002).

## Appendix D CUT-SKY FISHER MATRIX

The WMAP sky maps have nearly diagonal pixel-pixel noise covariance (Hinshaw et al., 2003a) which greatly simplifies the properties of the power spectrum Fisher matrix. In this Appendix we present an analytic derivation of the effect of a sky cut and non-uniform pixel weighting. In D.1, we assume that we have an optimal estimator of the power spectrum. In the noise dominated limit, we can obtain an exact expression, while in the signal dominated limit, we need to approximate the signal correlation matrix to obtain an analytic expression. In D.2, we interpolate the Fisher matrix between the signal and noise-dominated regimes and show that it agrees with numerical estimates. In D.3, we estimate the power spectrum covariance matrix from Monte Carlo simulations of the sky and calibrate the interpolation formula.

### D.1 Cut-Sky Fisher Matrix: Analytic Evaluation

We equate the Fisher matrix of the power spectrum to the curvature of the likelihood function, equation (A16), then develop an approximate form for the covariance matrix, $`𝐂=𝐒+𝐍`$, where $`𝐒`$ is the signal matrix and $`𝐍`$ is the noise matrix. We split the noise matrix into two pieces: a weight term and a mask term. In the limit that the pixel noise is diagonal, the weight term has the form

$$\left(𝐍_𝐰^\mathrm{𝟏}\right)_{ij}=\frac{n_i}{\sigma _0^2}\delta _{ij}w_i\delta _{ij}$$
(D1)
where $`i`$ and $`j`$ are pixel indices, $`n_i`$ is the number of observations of pixel $`i`$, $`\sigma _0`$ is the rms noise of a single observation, and, by definition, $`w_i`$ is the weight of pixel $`i`$. The mask, $`M_i`$, is defined so that $`M_i`$ equals 0 in pixels that are not used due to foreground contamination, and equals 1 otherwise. The noise matrix can then be written as the product of the two terms

$$𝐍^\mathrm{𝟏}=𝐍_𝐰^\mathrm{𝟏}𝐌=\mathrm{𝐌𝐍}_𝐰^\mathrm{𝟏}=\frac{\widehat{n}_i}{\sigma _0^2}\delta _{ij}=\widehat{w}_i\delta _{ij}$$
(D2)
where $`\widehat{n}_i=n_i`$ and $`\widehat{w}_i=w_i`$ in the unmasked pixels, and are set to 0 otherwise. We thus define the covariance matrix over the full sky, which allows us to exploit the orthogonality properties of the spherical harmonics. Note that $`𝐌^2=𝐌`$.

The inverse of the full covariance matrix can now be written as

$`𝐂^\mathrm{𝟏}=(𝐒+𝐍)^\mathrm{𝟏}`$ $`=`$ $`𝐍^\mathrm{𝟏}(\mathrm{𝐒𝐍}^\mathrm{𝟏}+𝐈)^\mathrm{𝟏}`$ (D3)
$`=`$ $`𝐌^\mathrm{𝟏}(\mathrm{𝐒𝐌}^\mathrm{𝟏}+𝐍_𝐰)^\mathrm{𝟏}.`$ (D4)
The covariance matrix has two limits. In the noise dominated limit, $`\mathrm{𝐒𝐍}^\mathrm{𝟏}𝐈`$,

$$𝐂^\mathrm{𝟏}𝐍^\mathrm{𝟏}.$$
(D5)
In the signal dominated limit, $`\mathrm{𝐒𝐌}^\mathrm{𝟏}𝐍_𝐰`$, only the mask alters the covariance matrix, so

$$𝐂^\mathrm{𝟏}\mathrm{𝐌𝐒}^\mathrm{𝟏}𝐌,$$
(D6)
where we have set the inverse of the mask to zero where there is no data, i.e. $`𝐌^\mathrm{𝟏}𝐌`$.

#### D.1.1 Cut-Sky Fisher Matrix in the Noise Dominated Limit

In the noise dominated limit, the Fisher matrix, equation (A16), is approximately

$$F_{ll^{}}^\mathrm{N}=\frac{1}{2}\mathrm{tr}\left(𝐂^\mathrm{𝟏}𝐏_l𝐂^\mathrm{𝟏}𝐏_l^{}\right)\frac{1}{2}\mathrm{tr}\left(𝐍^\mathrm{𝟏}𝐏_l𝐍^\mathrm{𝟏}𝐏_l^{}\right).$$
(D7)
Using

$$\left(𝐏_l\right)_{ij}\frac{𝐂}{C_l}=\frac{(2l+1)}{4\pi }w_lP_l(\mathrm{cos}\theta _{ij})=w_l\underset{m}{}Y_{lm}^{}(p_i)Y_{lm}(p_j)$$
(D8)
where $`\theta _{ij}`$ is the angle between pixels $`i`$ and $`j`$, $`Y_{lm}(p_i)`$ is a spherical harmonic evaluated in the direction of pixel $`i`$, and we have employed the addition theorem for spherical harmonics in the last equality. Then the Fisher matrix takes the form

$$F_{ll^{}}^\mathrm{N}=\frac{1}{2}w_lw_l^{}\underset{mm^{}}{}\underset{ij}{}\widehat{w}_iY_{lm}^{}(p_i)Y_{lm}(p_j)\widehat{w}_jY_{l^{}m^{}}^{}(p_j)Y_{l^{}m^{}}(p_i).$$
(D9)
Now expand the weight array as

$$\widehat{w}_i=\underset{lm}{}\widehat{w}_{lm}Y_{lm}(p_i)$$
(D10)
and substitute this into equation (D9) to yield

$$F_{ll^{}}^\mathrm{N}=\frac{1}{2}w_lw_l^{}\underset{mm^{}}{}\underset{l^{\prime \prime }m^{\prime \prime }}{}\underset{l^{\prime \prime \prime }m^{\prime \prime \prime }}{}\widehat{w}_{l^{\prime \prime }m^{\prime \prime }}\widehat{w}_{l^{\prime \prime \prime }m^{\prime \prime \prime }}\underset{ij}{}Y_{l^{\prime \prime }m^{\prime \prime }}(p_i)Y_{lm}^{}(p_i)Y_{lm}(p_j)Y_{l^{\prime \prime \prime }m^{\prime \prime \prime }}(p_j)Y_{l^{}m^{}}^{}(p_j)Y_{l^{}m^{}}(p_i).$$
(D11)
Since $`\widehat{w}_i`$ was defined over the full sky, the sum over pixels may be expressed in terms of products of Wigner 3-$`j`$ symbols. As shown in Appendix E, we may use an orthogonality property of the 3-$`j`$ symbols to reduce the expression for the Fisher matrix to

$$F_{ll^{}}^\mathrm{N}=\frac{1}{2}w_lw_l^{}\frac{(2l+1)(2l^{}+1)}{4\pi \mathrm{\Omega }_p^2}\underset{l^{\prime \prime }}{}\widehat{W}_{l^{\prime \prime }}\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime }\\ 0\hfill & 0& \hfill 0\end{array}\right)^2$$
(D12)
where $`\widehat{W}_l_m|\widehat{w}_{lm}|^2`$, and $`\mathrm{\Omega }_p`$ is the solid angle per pixel.

#### D.1.2 Cut-Sky Fisher Matrix in the Signal Dominated Limit

In the signal dominated limit, the Fisher matrix, equation (A16), is approximately

$$F_{ll^{}}^\mathrm{S}=\frac{1}{2}\mathrm{tr}\left(𝐂^\mathrm{𝟏}𝐏_l𝐂^\mathrm{𝟏}𝐏_l^{}\right)\frac{1}{2}\mathrm{tr}\left(\widehat{𝐒}^\mathrm{𝟏}𝐏_l\widehat{𝐒}^\mathrm{𝟏}𝐏_l^{}\right),$$
(D13)
where $`\widehat{𝐒}^\mathrm{𝟏}M_iS_{ij}^1M_j`$. Using equation (D8) we can write $`F_{ll^{}}^\mathrm{S}`$ in the form

$$F_{ll^{}}^\mathrm{S}=\frac{1}{2}w_lw_l^{}\underset{mm^{}}{}\underset{ijk}{}M_iM_jS_{ij}^1Y_{lm}^{}(p_j)Y_{lm}(p_k)\underset{kp}{}M_kM_pS_{kp}^1Y_{l^{}m^{}}^{}(p_p)Y_{l^{}m^{}}(p_i).$$
(D14)
Substituting a normalized expression for $`𝐒^\mathrm{𝟏}`$ in the pixel basis,

$$S_{ij}^1=\mathrm{\Omega }_p^2\underset{lm}{}\frac{1}{C_lw_l}Y_{lm}^{}(p_i)Y_{lm}(p_j),$$
(D15)
into equation (D14), one obtains a term

$$\mathrm{\Omega }_p^2\underset{l^{\prime \prime }}{}\frac{1}{C_{l^{\prime \prime }}w_{l^{\prime \prime }}}\underset{m^{\prime \prime }}{}\underset{i}{}M_iY_{l^{}m^{}}(p_i)Y_{l^{\prime \prime }m^{\prime \prime }}^{}(p_i)\underset{j}{}M_jY_{l^{\prime \prime }m^{\prime \prime }}(p_j)Y_{lm}^{}(p_j).$$
(D16)
Since our mask cuts only $``$15% of the sky, the coupling sum, $`_iM_iY_{l^{}m^{}}(p_i)Y_{l^{\prime \prime }m^{\prime \prime }}^{}(p_i)`$, peaks very sharply at $`|l^{}l^{\prime \prime }|l^{}`$. Therefore, one may approximate equation (D16) with

$`{\displaystyle \frac{\mathrm{\Omega }_p^2}{(C_lw_lC_l^{}w_l^{})^{1/2}}}{\displaystyle \underset{l^{\prime \prime }m^{\prime \prime }}{}}{\displaystyle \underset{i}{}}M_iY_{l^{}m^{}}(p_i)Y_{l^{\prime \prime }m^{\prime \prime }}^{}(p_i){\displaystyle \underset{j}{}}M_jY_{l^{\prime \prime }m^{\prime \prime }}(p_j)Y_{lm}^{}(p_j)`$ (D17)
$`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_p}{(C_lw_lC_l^{}w_l^{})^{1/2}}}{\displaystyle \underset{i}{}}M_iY_{l^{}m^{}}(p_i)Y_{lm}^{}(p_i),`$
where, in the last equality, we have used the completeness relation for the spherical harmonics

$$\underset{lm}{}Y_{lm}^{}(p_i)Y_{lm}(p_j)=\mathrm{\Omega }_p^1\delta _{ij}.$$
(D18)
With this approximation equation (D14) reduces to

$$F_{ll^{}}^\mathrm{S}\frac{1}{2}\frac{\mathrm{\Omega }_p^2}{C_lC_l^{}}\underset{mm^{}}{}\underset{ij}{}M_iY_{lm}^{}(p_i)Y_{lm}(p_j)M_jY_{l^{}m^{}}^{}(p_j)Y_{l^{}m^{}}(p_i).$$
(D19)
Following the same steps outlined above for the derivation of equation (D12) yields

$$F_{ll^{}}^\mathrm{S}=\frac{1}{2}\frac{1}{C_lC_l^{}}\frac{(2l+1)(2l^{}+1)}{4\pi }\underset{l^{\prime \prime }}{}M_{l^{\prime \prime }}\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime }\\ 0\hfill & 0& \hfill 0\end{array}\right)^2$$
(D20)
where $`M_l_m|m_{lm}|^2`$, and $`m_{lm}`$ is the spherical harmonic transform of the mask. Note that $`M_0=4\pi f_{\mathrm{sky}}^2`$, where $`f_{\mathrm{sky}}`$ is the fraction of the sky that survives the mask.

### D.2 Interpolating the Cut-Sky Fisher Matrix

We can combine the two limiting cases of the Fisher matrix, obtained in the previous section, into a single expression

$$F_{ll^{}}^{\mathrm{mask}}=\frac{\sqrt{(2l+1)(2l^{}+1)}f_{\mathrm{sky}}^2}{2(C_l+N_l)(C_l^{}+N_l^{})}\stackrel{~}{F}_{ll^{}},$$
(D21)
where $`N_l\mathrm{\Omega }_p\sigma _0^2/(\overline{n}_{obs}w_l)`$ is the deconvolved noise power spectrum. Using equations (D12) and (D20) derived in the previous section, $`\stackrel{~}{F}_{ll^{}}`$ can be expressed in the low $`l`$ (signal-dominated) and high $`l`$ (noise-dominated) limits as

$$\stackrel{~}{F}_{ll^{}}^\mathrm{S}=\frac{\sqrt{(2l+1)(2l^{}+1)}}{4\pi f_{\mathrm{sky}}^2}\underset{l^{\prime \prime }}{}M_{l^{\prime \prime }}\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime }\\ 0\hfill & 0& \hfill 0\end{array}\right)^2,$$
(D22)
and

$$\stackrel{~}{F}_{ll^{}}^\mathrm{N}=\frac{\sqrt{(2l+1)(2l^{}+1)}}{4\pi f_{\mathrm{sky}}^2\overline{n}_{obs}^2/\sigma _0^4}\underset{l^{\prime \prime }}{}\widehat{W}_{^{\prime \prime }}\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime }\\ 0\hfill & 0& \hfill 0\end{array}\right)^2.$$
(D23)
By construction, these matrices are normalized to 1 on the diagonal, $`\stackrel{~}{F}_{ll}^\mathrm{S}=\stackrel{~}{F}_{ll}^\mathrm{N}=1`$, since $`M_0=4\pi f_{\mathrm{sky}}^2`$ and $`\widehat{W}_0=4\pi f_{\mathrm{sky}}^2\overline{n}_{obs}^2/\sigma _0^4`$.

We have computed $`F_{ll^{}}^{\mathrm{mask}}`$ directly, for selected $`l`$ values, by evaluating $`Y_{lm}𝐂^1Y_{l^{}m^{}}`$ using the pre-conditioner code described by Oh et al. (1999). We find that deviations between the numerical and analytical results are consistent with numerical noise in the Fisher matrix estimate. Figure 14 shows the estimated Fisher matrix for the Kp2 cut in the noise and the signal dominated limits.

We interpolate between the two regimes with the following expression:

$$F_{ll^{}}^{\mathrm{mask}}=\frac{\left(C_l\sqrt{F_{ll^{}}^\mathrm{S}}+N_l\sqrt{F_{ll^{}}^\mathrm{N}}\right)\left(C_l^{}\sqrt{F_{ll^{}}^\mathrm{S}}+N_l^{}\sqrt{F_{ll^{}}^\mathrm{N}}\right)}{(C_l+N_l)(C_l^{}+N_l^{})}.$$
(D24)
Note that the form of the Fisher matrix primarily depends on $`\mathrm{\Delta }l`$, but is weakly dependent on $`l`$ – there is more coupling in the noise dominated limit.

### D.3 Power Spectrum Covariance: Monte Carlo Evaluation

As discussed in Appendix A.1.2 we compute the $`C_l`$ using three different pixel weightings. The uniform weighting and the inverse-noise inverse weighting are optimal in the signal-dominated regime and the noise-dominated regime, respectively. In between these limits the transitional weighting performs better. In order to determine which ranges in $`l`$ correspond to which regimes, and to calibrate our ansatz for the covariance matrix, equation (D21), we proceed as follows. Using 100,000 Monte Carlo simulations of signal plus noise (with WMAP noise levels and symmetrized beams), we compute the diagonal elements of the covariance matrix for the three different weighting schemes, evaluated with the Kp2 sky cut. Denote these estimates $`D_l^{\mathrm{sim}}`$

We find that the uniform weighting produces the smallest $`D_l^{\mathrm{sim}}`$ below $`l=200`$. Inverse-noise weighting is the best above $`l=450`$, and the transitional weighting produces the lowest variance in between. In each of these regimes we use the resulting $`D_l^{\mathrm{sim}}`$ to calibrate our ansatz for the diagonal elements of the covariance matrix

$$D_l=\frac{2}{2l+1}\frac{1}{f_{\mathrm{sky}}^{\mathrm{eff}}}(C_l+N_l^{\mathrm{eff}})^2,$$
(D25)
as illustrated in Verde et al. (2003). These calibrations produce a smooth correction to equation (D24) of at most 6%. No correction at all is required in the signal dominated regime.

## Appendix E SOME USEFUL PROPERTIES OF SPHERICAL HARMONICS

The derivation of the form of the Fisher matrix in the signal and noise dominated limits led to expressions which included a term

$$\underset{mm^{}}{}\underset{ij}{}Y_{l^{\prime \prime }m^{\prime \prime }}(p_i)Y_{lm}^{}(p_i)Y_{lm}(p_j)Y_{l^{\prime \prime \prime }m^{\prime \prime \prime }}(p_j)Y_{l^{}m^{}}^{}(p_j)Y_{l^{}m^{}}(p_i),$$
(E1)
where the sum is effectively a double integral over the full sky. This can be evaluated in terms of the Wigner 3-$`j`$ symbols, defined as

$$𝑑\mathrm{\Omega }_𝐧Y_{lm}(𝐧)Y_{l^{}m^{}}^{}(𝐧)Y_{l^{\prime \prime }m^{\prime \prime }}(𝐧)\left(1\right)^m^{}\sqrt{\frac{(2l+1)(2l^{}+1)(2l^{\prime \prime }+1)}{4\pi }}\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime }\\ 0\hfill & 0& \hfill 0\end{array}\right)\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime }\\ m\hfill & m^{}& \hfill m^{\prime \prime }\end{array}\right).$$
(E2)
Substituting this into equation (E1) gives

$`{\displaystyle \underset{mm^{}}{}}\left(1\right)^{m+m^{}}{\displaystyle \frac{(2l+1)(2l^{}+1)}{4\pi }}\sqrt{(2l^{\prime \prime }+1)(2l^{\prime \prime \prime }+1)}\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime }\\ 0\hfill & 0& \hfill 0\end{array}\right)\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime \prime }\\ 0\hfill & 0& \hfill 0\end{array}\right)`$ (E7)
$`\times \left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime }\\ m\hfill & m^{}& \hfill m^{\prime \prime }\end{array}\right)\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime \prime }\\ m\hfill & m^{}& \hfill m^{\prime \prime \prime }\end{array}\right).`$ (E12)
We simplify this using the Wigner 3-$`j`$ orthogonality condition

$$\underset{mm^{}}{}(1)^{m+m^{}}\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime }\\ m\hfill & m^{}& \hfill m^{\prime \prime }\end{array}\right)\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime \prime }\\ m\hfill & m^{}& \hfill m^{\prime \prime \prime }\end{array}\right)=\frac{1}{2l^{\prime \prime }+1}\delta _{l^{\prime \prime }l^{\prime \prime \prime }}\delta _{m^{\prime \prime }m^{\prime \prime \prime }}\delta (l,l^{},l^{\prime \prime }),$$
(E13)
where $`\delta (l,l^{},l^{\prime \prime })=1`$ for $`|ll^{}|l^{\prime \prime }l+l^{}`$ and is 0 otherwise. This reduces equation (E1) to

$$\frac{(2l+1)(2l^{}+1)}{4\pi \mathrm{\Omega }_p^2}\left(\begin{array}{ccc}l\hfill & l^{}& \hfill l^{\prime \prime }\\ 0\hfill & 0& \hfill 0\end{array}\right)^2,$$
(E14)
where the factor of $`\mathrm{\Omega }_p^2`$ accounts for the fact that equation (E1) is a double sum over pixels, instead of a double integral.