# Three-Year Wilkinson Microwave Anisotropy Probe
(WMAP^{1}^{1}affiliation: \map is the result of a partnership between Princeton
University and NASA’s Goddard Space Flight Center. Scientific guidance is
provided by the \map Science Team. ) Observations:

Temperature Analysis

###### Abstract

We present new full-sky temperature maps in five frequency bands from 23 to 94 GHz, based on data from the first three years of the \map sky survey. The new maps are consistent with the first-year maps and are more sensitive. The three-year maps incorporate several improvements in data processing made possible by the additional years of data and by a more complete analysis of the polarization signal. These include several new consistency tests as well as refinements in the gain calibration and beam response models (Jarosik et al., 2006).

We employ two forms of multi-frequency analysis to separate astrophysical foreground signals from the CMB, each of which improves on our first-year analyses. First, we form an improved “Internal Linear Combination” (ILC) map, based solely on \map data, by adding a bias correction step and by quantifying residual uncertainties in the resulting map. Second, we fit and subtract new spatial templates that trace Galactic emission; in particular, we now use low-frequency \map data to trace synchrotron emission instead of the 408 MHz sky survey. The \map point source catalog is updated to include 115 new sources whose detection is made possible by the improved sky map sensitivity.

We derive the angular power spectrum of the temperature anisotropy using a hybrid approach that combines a maximum likelihood estimate at low (large angular scales) with a quadratic cross-power estimate for . The resulting multi-frequency spectra are analyzed for residual point source contamination. At 94 GHz the unmasked sources contribute K to at . After subtracting this contribution, our best estimate of the CMB power spectrum is derived by averaging cross-power spectra from 153 statistically independent channel pairs. The combined spectrum is cosmic variance limited to , and the signal-to-noise ratio per -mode exceeds unity up to . For bins of width %, the signal-to-noise ratio exceeds unity up to . The first two acoustic peaks are seen at and , respectively, while the first two troughs are seen at and , respectively. The rise to the third peak is unambiguous; when the \map data are combined with higher resolution CMB measurements, the existence of a third acoustic peak is well established.

Spergel et al. (2006) use the three-year temperature and polarization data to constrain cosmological model parameters. A simple six parameter CDM model continues to fit CMB data and other measures of large scale structure remarkably well. The new polarization data (Page et al., 2006) produce a better measurement of the optical depth to re-ionization, . This new and tighter constraint on helps break a degeneracy with the scalar spectral index which is now found to be . If additional cosmological data sets are included in the analysis, the spectral index is found to be .

^{†}

^{†}slugcomment: ApJ, in press, January 5, 2007

## 1 Introduction

The Wilkinson Microwave Anisotropy Probe (\map) is a Medium-class Explorer (MIDEX) mission designed to elucidate cosmology by producing full-sky maps of the cosmic microwave background (CMB) anisotropy. Results from the first year of \map observations were reported in a suite of papers published in the Astrophysical Journal Supplement Series in September 2003 (Bennett et al., 2003b; Jarosik et al., 2003a; Page et al., 2003a; Barnes et al., 2003; Hinshaw et al., 2003a; Bennett et al., 2003c; Komatsu et al., 2003; Hinshaw et al., 2003b; Kogut et al., 2003; Spergel et al., 2003; Verde et al., 2003; Peiris et al., 2003; Page et al., 2003c; Bennett et al., 2003a; Page et al., 2003b; Barnes et al., 2002; Jarosik et al., 2003b; Nolta et al., 2004). The data were made available to the research community via the Legacy Archive for Microwave Background Data Analysis (LAMBDA), NASA’s CMB Thematic Data Center, and were described in detail in the \map Explanatory Supplement (Limon et al., 2003).

Papers based on the first-year \map results cover a wide range of topics, including: constraints on inflation, the nature of the dark energy, the dark matter density, implications for supersymmetry, the CMB and \map as the premier baryometer, intriguing features in the large-scale data, the topology of the universe, deviations from Gaussian statistics, time-variable cosmic parameters, the Galactic interstellar medium, microwave point sources, the Sunyaev-Zeldovich effect, and the ionization history of the universe. The \map data has also been used to establish the calibration of other CMB data sets.

Our analysis of the first three years of \map data is now complete and the results are presented here and in companion papers (Jarosik et al., 2006; Page et al., 2006; Spergel et al., 2006). The three-year \map results improve upon the first-year set in many ways, the most important of which are the following. (1) A thorough analysis of the polarization data has produced full-sky polarization maps and power spectra, and an improved understanding of many aspects of the data. (2) Additional data reduces the instrument noise, producing power spectra that are 3 times more sensitive in the noise limited regime. (3) Independent years of data enable cross-checks that were not previously possible. (4) The instrument calibration and beam response have been better characterized.

This paper presents the analysis of the three-year temperature data, focusing on foreground modeling and removal, evaluation of the angular power spectrum, and selected topics beyond the power spectrum. Companion papers present the new polarization maps and polarization-specific scientific results (Page et al., 2006), and discuss the cosmological implications of the three-year \map data (Spergel et al., 2006). Jarosik et al. (2006) present our new data processing methods and place systematic error limits on the maps.

In §2 we summarize the major changes we have made to the data processing since the first-year analysis, and §3 presents a synopsis of the three-year temperature maps. In §4 we discuss Galactic foreground emission and our attempts to separate the emission components using a Maximum Entropy Method (MEM) analysis. §5 illustrates two methods we employ to remove Galactic foreground emission from the maps in preparation for CMB analysis. §6 updates the \map point source catalog and presents a search for the Sunyaev-Zeldovich effect in the three-year maps. §7 evaluates the angular power spectrum and compares it to the previous \map spectrum and to other contemporary CMB results. In §8 we survey the claims that have been made regarding odd features in the \map first-year sky maps, and we offer conclusions in §9.

## 2 Changes in the Three-Year Data Analysis

The first-year data analysis was described in detail in the suite of first-year \map papers listed above. In large part, the three-year analysis employs the same methods, with the following exceptions.

In the first-year analysis we subtracted the COBE dipole from the time-ordered data to minimize the effect of signal aliasing that arises from pixelizing a signal with a steep gradient. Since the \map gain calibration procedure uses the Doppler effect induced by \map’s velocity with respect to the Sun to establish the absolute calibration scale, \map data may be used independently to determine the CMB dipole. Consequently, we subtract the \map first-year dipole (Bennett et al., 2003b) from the time-ordered data in the present analysis.

A small temperature dependent pointing error ( arcmin) was found during the course of the first-year analysis. The effect is caused by thermal stresses on the spacecraft structure that induce slight movement of the star tracker with respect to the instrument. While the error was small enough to ignore in the first-year data, it is now corrected with a temperature dependent model of the relative motion (Jarosik et al., 2006).

The radiometer gain model described by Jarosik et al. (2003b) has been updated to include a dependence on the temperature of the warm-stage (RXB) amplifiers. While this term was not required by the first-year data, it is required now for the model to fit the full three-year data with a single parameterization. The new model, and its residual errors, are discussed by Jarosik et al. (2006).

The \map beam response has now been measured with six independent “seasons” of Jupiter observations. In addition, we have now produced a physical model of one side of our symmetric optical system, the A-side, based on simultaneous fits to all 10 A-side beam pattern measurements (Jarosik et al., 2006). We use this model to augment the beam response data at very low signal-to-noise ratio, which in turn allows us to determine better the total solid angle and window function of each beam.

The far sidelobe response of the beam was determined from a combination of ground measurements and in-flight lunar data taken early in the mission (Barnes et al., 2003). In the first-year processing we applied a small far-sidelobe correction to the K-band sky map. For the current analysis, we have implemented a new far sidelobe correction and gain re-calibration that operates on the time-ordered data (Jarosik et al., 2006). These corrections have now been applied to data from all 10 differencing assemblies.

When producing polarization maps, we account for differences in the frequency pass-band between the two linear polarization channels in a differencing assembly (Page et al., 2006). If this difference is not accounted for, Galactic foreground signals would alias into linear polarization signals.

Due to a combination of 1/f noise and observing strategy, the noise in the \map sky maps is correlated from pixel to pixel. This results in certain low- modes on the sky being less well measured than others. This effect can be completely ignored for temperature analysis since the low- signal-to-noise ratio is so high, and the effect is not important at high- (§7.1.2). However, it is very important for polarization analysis because the signal-to-noise ratio is so much lower. In order to handle this complexity, the map-making procedure has been overhauled to produce genuine maximum likelihood solutions that employ optimal filtering of the time-ordered data and a conjugate-gradient algorithm to solve the linear map-making equations (Jarosik et al., 2006). In conjunction with this we have written code to evaluate the full pixel-to-pixel weight (inverse covariance) matrix at low pixel resolution. (The HEALPix convention is to denote pixel resolution by the parameter , with (Gorski et al., 2005). We define a resolution parameter such that . The weight matrices have been evaluated at resolution r4, , .) The full noise covariance information is propagated through the power spectrum analysis (Page et al., 2006).

When performing template-based Galactic foreground subtraction, we now use templates based on \map K- and Ka-band data in place of the 408 MHz synchrotron map (Haslam et al., 1981). As discussed in §5.3, this substitution reduces errors caused by spectral index variations that change the spatial morphology of the synchrotron emission as a function of frequency. A similar model is used for subtracting polarized synchrotron emission from the polarization maps (Page et al., 2006).

We have performed an error analysis of the internal linear combination (ILC) map and have now implemented a bias correction as part of the algorithm. We believe the map is now suitable for use in low- CMB signal characterization, though we have not performed a full battery of non-Gaussian tests on this map, so we must still advise users to exercise caution. Accordingly, we present full-sky multipole moments for , derived from the three-year ILC map.

We have improved the final temperature power spectrum () by using a maximum likelihood estimate for low- and a pseudo- estimate for (see §7). The pseudo- estimate is simplified by using only V- and W-band data, and by reducing the number of pixel weighting schemes to two, “uniform” and “” (§7.5). With three individual years of data and six V- and W-band differencing assemblies (DAs) to choose from, we can now form individual cross-power spectra from 15 DA pairs within a year and from 36 DA pairs across 3 year pairs, for a total of 153 independent cross-power spectra. In the first-year spectrum we included Q-band data, which gave us 8 DAs and 28 independent cross-power spectra. The arguments for dropping Q-band from the three-year spectrum are given in §7.2.

We have developed methods for estimating the polarization power spectra ( for XX = TE, TB, EE, EB, BB) from temperature and polarization maps. The main technical hurdle we had to overcome in the process was the proper handling of low signal-to-noise ratio data with complex noise properties (Page et al., 2006). This step, in conjunction with the development of the new map-making process, was by far the most time consuming aspect of the three-year analysis.

We have improved the form of the likelihood function used to infer cosmological parameters from the Monte Carlo Markov Chains (Spergel et al., 2006). In addition to using an exact maximum likelihood form for the low- TT data, we have developed a method to self-consistently evaluate the joint likelihood of temperature and polarization data given a theoretical model (described in Appendix D of Page et al. (2006)). We also now account for Sunyaev-Zeldovich (SZ) fluctuations when estimating parameters. Within the \map frequency range, it is difficult to distinguish between a primordial CMB spectrum and a thermal SZ spectrum, so we adopt the Komatsu & Seljak (2002) model for the SZ power spectrum and marginalize over the amplitude as a nuisance parameter.

## 3 Observations and Maps

The three-year \map data encompass the period from 00:00:00 UT, 10 August 2001 (day number 222) to 00:00:00 UT, 9 August 2004 (day number 222). The observing efficiency during this time is roughly 99%; Table 2 lists the fraction of data that was lost or flagged as suspect. The Table also gives the fraction of data that is flagged due to potential contamination by thermal emission from Mars, Jupiter, Saturn, Uranus, and Neptune. These data are not used in map-making, but are useful for in-flight beam mapping (Limon et al., 2006).

Sky maps are created from the time-ordered data using the procedure described by Jarosik et al. (2006). For several reasons, we produce single-year maps for each year of the three-year observing period (after performing an end-to-end analysis of the instrument calibration). We produce three-year maps by averaging the annual maps. Figure 1 shows the three-year maps at each of the five \map observing frequencies: 23, 33, 41, 61, and 94 GHz. The number of independent observations per pixel, , is displayed in Figure 2. The noise per pixel, , is given by , where is the noise per observation, given in Table 1. To a very good approximation, the noise per pixel in the three-year maps is a factor of times lower than in the one-year maps. The noise properties of the data are discussed in more detail in Jarosik et al. (2006).

The three-year maps are compared to the previously released maps in Figure 3. Both set of maps have been smoothed to 1 resolution to minimize the noise difference between them. When viewed side by side they look indistinguishable. The right column of Figure 3 shows the difference of the maps at each frequency on a scale of 30 K. Aside from the noise reduction and a few bright variable quasars, such as 3C279, the main difference between the maps is in the large-scale (low-) emission. This is largely due to improvements in our model of the instrument gain as a function of time, which is made possible by having a longer time span with which to fit the model (Jarosik et al., 2006). In the specific case of K-band, the improved far-sidelobe pickup correction produced an effective change in the absolute calibration scale by 1%. This, in turn, is responsible for the difference seen in the bright Galactic plane signal in K-band (Jarosik et al., 2006). We discuss the low- emission in detail in §7.4 and §8, but we stress here that the changes shown in Figure 3 are small, even compared to the low quadrupole moment seen in the first-year maps. Table 3 gives the amplitude of the dipole, quadrupole, and octupole moments in these difference maps. For comparison, we estimate the CMB power at to be = 236 and 1053 K, respectively (§7.4).

As discussed in §8, several authors have noted unusual features in the large-scale signal recorded in the first-year maps. We have not attempted to reproduce the analyses presented in those papers, but based on the small fractional difference in the large-scale signal, we anticipate that most of the previously reported results will persist when the three-year maps are analyzed.

## 4 Galactic Foreground Analysis

The CMB signal in the \map sky maps is contaminated by microwave emission from the Milky Way Galaxy and from extragalactic sources. In order to use the maps reliably for cosmological studies, the foreground signals must be understood and removed from the maps. In this section we present an overview of the mechanisms that produce significant diffuse microwave emission in the Milky Way and we assess what can be learned about them using a Maximum Entropy Method (MEM) analysis of the \map data. We discuss foreground removal in §5.

### 4.1 Free-Free Emission

Free-free emission arises from electron-ion scattering which produces microwaves with a brightness spectrum for frequencies GHz, where is the emission measure, , and we assume an electron gas temperature K. As discussed in Bennett et al. (2003c), high-resolution maps of H emission (Dennison et al., 1998; Haffner et al., 2003; Reynolds et al., 2002; Gaustad et al., 2001) can serve as approximate tracers of free-free emission. The intensity of H emission is given by

(1) |

where is in Rayleighs (1 R ergs cms sr at the H wavelength of 0.6563 m), the helium contribution is assumed to be small, and is an extinction factor that depends on the dust optical depth, , at the wavelength of H. If the emitting gas is co-extensive with dust, then . H is in R-band, where the extinction is 0.75 times visible, ; thus, , and . Finkbeiner (2003) assembled a full-sky H map using data from several surveys: the Wisconsin H-Alpha Mapper (WHAM), the Virginia Tech Spectral-Line Survey (VTSS), and the Southern H-Alpha Sky Survey Atlas (SHASSA). We use this map, together with the Schlegel et al. (1998) (SFD) extinction map, to predict a map of free-free emission in regions where , under the assumption that the dust and ionized gas are co-extensive. As discussed in Bennett et al. (2003c), this template has known sources of uncertainty and error. We use it as a prior estimate in the MEM analysis (§4.5), and as a free-free estimate in the template-based foreground removal (§5.3).

### 4.2 Synchrotron Emission

Synchrotron emission arises from the acceleration of cosmic ray electrons in magnetic fields. In our Galaxy, discrete supernova remnants contribute only 10% of the total synchrotron emission at 1.5 GHz (Lisenfeld & Völk, 2000; Biermann, 1976; Ulvestad, 1982), while 90% of the observed emission arises from a diffuse component. Hummel et al. (1991) present maps of synchrotron emission at 610 MHz and 1.49 GHz from the edge-on spiral galaxy NGC891. They find that the synchrotron spectral index varies from in most of the galactic plane to in the halo. Similar spectral index variations are seen in the Milky Way at 1 GHz, where the synchrotron signal is complex. Variations of the synchrotron spectral index are both expected and observed. Moreover, the emission is dominated at low frequencies by components with steep spectra, whereas at higher frequencies it is dominated by components with flatter spectra, usually with a different spatial distribution. As a result, great care must be taken when using low frequency maps, like the 408 MHz map of Haslam et al. (1981), as tracers of the synchrotron emission at microwave frequencies.

Synchrotron emission can be highly polarized. Theoretically, the linear polarization fraction can be as high as 75%, though values % are more typically observed. See Page et al. (2006) for a discussion of the new full-sky observations of polarized synchrotron emission in the three-year \map data.

### 4.3 Thermal Dust Emission

Thermal dust emission has been mapped over the full sky in several infrared bands by the IRAS and COBE missions. Schlegel et al. (1998) combined data from both missions to produce an absolutely calibrated full-sky map of the thermal dust emission. Finkbeiner et al. (1999) (FDS) extended this work to far-infrared and microwave frequencies using the COBE-FIRAS and COBE-DMR data to constrain the low-frequency dust spectrum. They fit the data to a particular two-component model that gives power-law emissivity indices and , and temperatures of K and K. The fraction of power emitted by each component is and , and the relative ratio of IR thermal emission to optical opacity of the two components is . The cold component is potentially identified as emission from amorphous silicate grains while the warm component is plausibly carbon based. Independent of the physical interpretation of the model, FDS found that it fit the data moderately well, with for 118 degrees of freedom. Bennett et al. (2003a) noted that this model, call “Model 8” by FDS, did well predicting the first-year \map dust emission.

It is reasonable to assume that the Milky Way is like other spiral galaxies and that the microwave properties of external galaxies should help to inform our understanding of the global properties of the Milky Way. It has long been known that a remarkably tight correlation exists between the broadband far-infrared and broadband synchrotron emission in external galaxies. This relation has been extensively studied and modeled (Dickey & Salpeter, 1984; de Jong et al., 1985; Helou et al., 1985; Sanders & Mirabel, 1985; Gavazzi et al., 1986; Hummel, 1986; Wunderlich et al., 1987; Wunderlich & Klein, 1988; Beck & Golla, 1988; Fitt et al., 1988; Hummel et al., 1988; Mirabel & Sanders, 1988; Bicay et al., 1989; Devereux & Eales, 1989; Unger et al., 1989; Voelk, 1989; Chi & Wolfendale, 1990; Wunderlich & Klein, 1991; Condon, 1992; Bressan et al., 2002). All theories attempting to explain this tight correlation are tied to the level of the star formation activity. During this cycle, stars form, heat, and destroy dust grains; create magnetic fields and relativistic electrons; and create the O- and B-stars that ionize the surrounding gas. However, it is not clear what these models predict on a “microscopic” (cloud by cloud) level within a galaxy.

Bennett et al. (2003c) showed that the synchrotron and dust emission in our own Galaxy are spatially correlated at \map frequencies. Many authors have argued that this correlation is actually due to radio emission from dust grains themselves, rather than from a tight dust-synchrotron correlation. We review the evidence for this more fully in the next section.

### 4.4 “Anomalous” Microwave Emission from Dust?

With the advent of high-quality diffuse microwave emission maps in the early 1990s, it became possible to study the high-frequency tail of the synchrotron spectrum and the low-frequency tail of the interstellar dust spectrum. Kogut et al. (1996a, b) analyzed foreground emission in the COBE-DMR maps and reported a signal that was significantly correlated with 240 m dust emission (Arendt et al., 1998) but not with 408 MHz synchrotron emission (Haslam et al., 1981). The correlated signal was notably brighter at 31 GHz than at 53 GHz (), hence they concluded it was consistent with free-free emission that was spatially correlated with dust. The same conclusion was reached by de Oliveira-Costa et al. (1997), who found the Saskatoon 40 GHz data to be correlated with infrared dust, but not with radio synchrotron emission.

Leitch et al. (1997); Leitch (1998), and Leitch et al. (2000) analyzed data from the “RING5m” experiment. A likelihood fit to their 14.5 GHz and 31.7 GHz data, assuming CMB anisotropy and a single foreground component, produced a foreground spectral index of . The data would have preferred a steeper value had it not been for an assumed prior limit of . This signal was fully consistent with synchrotron emission. However, a puzzle arose in comparing the RING5m data with a the Westerbork Northern Sky Survey (WENSS) (Rengelink et al., 1997) at 325 MHz: the WENSS data showed no detectable signal in the vicinity of the RING5m field. The 325 MHz limit implies that and rules out conventional synchrotron emission as the dominant foreground. (Since the WENSS is an interferometric survey primarily designed to study discrete sources, the data are insensitive to zero-point flux from extended emission. It is not clear how much this affects the above conclusion.) As with the DMR and Saskatoon data, the 14.5 GHz foreground emission was correlated with dust, but it was difficult to attribute it to spatially correlated free-free emission because there was negligible H emission in the vicinity. To reconcile this, a gas temperature in excess of a million degrees would be needed to suppress the H. Flat spectrum synchrotron was also suggested as a possible source; it had been previously observed in other sky regions and it would obviate the need for such a high temperature and pressure.

Draine & Lazarian (1998) dismissed the hot ionized gas explanation on energetic grounds and instead suggested that the emission (which they described as “anomalous”) be attributed to electric dipole rotational emission from very small dust grains – a mechanism first proposed by Erickson (1957) in a different astrophysical context. One of the hallmarks of this mechanism is that it produces a frequency spectrum that peaks in the 10-60 GHz range and falls off fairly steeply on either side.

de Oliveira-Costa et al. (1998) analyzed the nearly full-sky 19 GHz sky map (Boughn et al., 1992) and found some correlation with the 408 MHz synchrotron emission, but found a stronger correlation with the COBE-DIRBE 240 m dust emission. They concluded the 19 GHz data were consistent with either free-free or spinning dust emission.

Leitch (1998) commented that the preferred model of Draine & Lazarian (1998) could produce the RING5m foreground component at 31.7 GHz, but that it only accounted for at most 30% of the 14.5 GHz emission, even when adopting unlikely values of the grain dipole moment. Since Leitch et al. (2000) were primarily interested in studying the CMB anisotropy, they considered using the IRAS 100 m map as a foreground template to remove the “anomalous” emission, regardless of its physical origin. They found, however, that fitting only a CMB component and a dust-correlated component produced an unacceptably high per degree of freedom. Thus, while the radio foreground morphology correlates with dust, the correlation is not perfect.

Finkbeiner et al. (2002) used the Green Bank 140 foot telescope to search for spinning dust emission in a set of dusty sources selected to be promising for detection. Ten infrared-selected dust clouds were observed at 5, 8, and 10 GHz. Eight of the ten sources yielded negative results, one was marginal, and one (the only diffuse H II region of the ten, LPH 201.6) was claimed as a tentative detection based on its spectral index of . Recognizing that this spectrum does not necessarily imply spinning dust emission, Finkbeiner et al. (2002) offer three additional requirements to convincingly demonstrate the detection of spinning dust, and concluded that none of the three requirements was met by the existing data. The absence of a rising spectrum in most of the sources may be taken as evidence that spinning dust emission is not typically dominant in this spectral region, at least for this type of infrared-selected cloud. The tentative detection in LPH 201.6 has met with three criticisms: (1) lack of evidence for the premise that its radio emission is proportional to its far-infrared dust emission (Casassus et al., 2004), (2) the putative spinning dust emission is stronger than theory predicts (McCullough & Chen, 2002), and (3) the positive spectral index may be accounted for by unresolved optically thick emission (McCullough & Chen, 2002). On the latter point, follow-up observations failed to identify a compact HII region candidate (McCullough, private communication).

Bennett et al. (2003c) fit the first-year \map foreground data to within 1% using a Maximum Entropy Method (MEM) analysis (see also §4.5). As with the above-cited results, \map found that the 22 GHz to 33 GHz foregrounds are dominated by a component with a synchrotron-like spectrum, but a dust-like spatial morphology. Bennett et al. (2003c) suggested that this may be due to spatially varying synchrotron spectral indices acting over a large frequency range, significantly altering the synchrotron morphology with frequency. A spinning dust component with a thermal dust morphology and the Draine and Lazarian spectrum could not account for more than % of the emission at 33 GHz. Of course, the \map fit did not rule out spinning dust as a sub-dominant emission source (as it surely must be at some level), nor did it rule out spinning dust models with other spectra or spatial morphologies.

Casassus et al. (2004) report evidence for anomalous microwave emission in the Helix planetary nebula at 31 GHz, where at least 20% of the emission is correlated with 100 m dust emission. They rejected several explanations. The observed features are not seen in H, ruling out free-free emission as the source. Cold grains are also ruled out as the source by the absence of 250 GHz continuum emission. Very small grains are not expected to survive in planetary nebulae, and none have been detected in the Helix, but Fe is strongly depleted in the gas. Instead, Casassus et al. (2004) favor the notion of magnetic dipole emission (produced by variations in grain magnetization) from hot ferromagnetic classical grains (Draine & Lazarian, 1999). Although the derived emissivity per nucleon in the Helix is a factor of larger than the highest end of the range predicted by Draine & Lazarian (1999), this excess could be explained by a high dust temperature, since Draine & Lazarian (1999) assume an ISM temperature of 18 K instead of a typical planetary nebula dust temperature of K. The fraction of 31 GHz Helix emission attributable to free-free is estimated to be in the range %. This low level of free-free emission implies an electron temperature of K, which is much lower than the value K based on collisionally excited lines. This discrepancy may be due to strong temperature variations within the nebula. Casassus et al. (2004) suggest that the Finkbeiner et al. (2002) measurement of LPH 201.6 may also be produced by magnetic dipole emission from classical dust grains.

de Oliveira-Costa et al. (2004) correlate the Tenerife 10 and 15 GHz data (Gutiérrez et al., 2000) with the \map non-thermal (“synchrotron”) map that was produced as part of the first-year Maximum Entropy Method (MEM) analysis of the \map foreground signal. They detect a low frequency roll-off in the correlated emission, as shown in Figure 6. We consider this result after discussing the three-year MEM analysis in the next section.

Fernandez-Cerezo et al. (2006) report new measurements with the COSMOSOMAS experiment covering 9000 square degrees with 1 resolution at frequencies of 12.7, 14.7, and 16.3 GHz. In addition to CMB signal and what is interpreted to be a population of unresolved radio sources, they find evidence for diffuse emission that is correlated with the DIRBE 100 m and 240 m bands. As with many of the above-mentioned results, they find that the correlated signal amplitude rises from 22 GHz (using \map data) to 16.3 GHz, and that it shows signs of flattening below 16.3 GHz “compatible with predictions for anomalous microwave emission related to spinning dust.”

The topic of anomalous dust emission remains unsettled, and is likely to remain so until high-quality diffuse measurements are available over a modest fraction of sky in the 5-15 GHz frequency range. We offer some further comments after presenting the three-year MEM results in the following section.

### 4.5 Maximum Entropy Method (MEM) Foreground Analysis

Bennett et al. (2003c) described a MEM-based approach to modeling the multifrequency \map sky maps on a pixel-by-pixel basis. Since the method is non-linear, it produces maps with complicated noise properties that are difficult to propagate in cosmological analyses. As a result, it is not a promising method for foreground removal. However, the method is quite useful in helping to separate Galactic foreground components by emission mechanism, which in turn informs our understanding of the foregrounds.

We model the temperature map at each frequency, , and pixel, , as

(2) |

where the subscripts cmb, s, ff, and d denote the CMB, synchrotron (including any anomalous dust component), free-free, and thermal dust components, respectively. is the spatial distribution of emission component in pixel , and is the spectrum of emission , which is not assumed to be uniform across the sky. We normalize the spectra () at K-band for the synchrotron and free-free components, and at W-band for the dust component.

The model is fit in each pixel by minimizing the functional (Press et al., 1992), where , is the standard of the model fit, and is the MEM functional (see below). The parameter controls the relative weight of (the data) and (the prior information) in the fit. In the functional , the sum over is restricted to Galactic emission components, and is a prior estimate of . The form of ensures the positivity of the solution for the Galactic components, which greatly alleviates degeneracy between the foreground components.

Throughout the MEM analysis, we smooth all maps to a uniform angular resolution. To improve our ability to constrain and understand the foreground components, we first subtract a prior estimate of the CMB signal from the data rather than fit for it. We use the ILC map described in §5.2 for this purpose and subtract it from all 5 frequency band maps. Since \map employs differential receivers, the zero level of each temperature map is unspecified. For the MEM analysis we adopt the following convention. In the limit that the Galactic emission is described by a plane-parallel slab, we have , where is the temperature at the Galactic pole. For each of the 5 frequency band maps, we assign the zero level such that a fit of the form , over the range , yields .

We construct a prior estimate for dust emission, , using Model 8 of Finkbeiner et al. (1999), evaluated at 94 GHz. The dust spectrum is modeled as a straight power law, . For free-free emission, we estimate the prior, , using the extinction-corrected H map (Finkbeiner, 2003). This is converted to a free-free signal using a conversion factor of 11.4 K R (units of antenna temperature at K-band). We model the spectrum as a straight power law, . As noted in §4.1, the main source of uncertainty in this free-free estimate is the level of extinction correction (in addition to any H photometry errors). We reduce in regions of high dust optical depth to minimize the effect of errors in the prior.

For the synchrotron emission, we construct a prior estimate, , by subtracting an extragalactic brightness of 5.9 K from the Haslam 408 MHz map (Lawson et al., 1987) and scaling the result to K-band assuming . Since the synchrotron spectrum varies with position on the sky, this prior estimate is expected to be imperfect. We account for this in choosing , as described below. We construct an initial spectral model for the synchrotron, , using the spectral index map . Specifically, we form for Ka-band, and for Q-, V-, and W-band ( for K-band). This allows for a -dependent steepening of the synchrotron spectrum at microwave frequencies. In our first-year results, use of this initial spectral model produced solutions with zero synchrotron signal, , in a few low latitude pixels for which the K-Ka spectrum is flatter than free-free emission. For the three-year analysis, this problem is handled by setting the initial model, , to be flatter than free-free emission in these pixels.

For all three emission components, the priors and the spectra are fixed during each minimization of . As described below, we iteratively improve the synchrotron spectrum model based on the residuals of the fit.

The parameter controls the degree to which the solution follows the prior. In regions where the signal is strong, the data alone should constrain the model without the need for prior information, so can be small (though we maintain to naturally impose positivity on and reduce degeneracy among the emission components). As in the first-year analysis, we base on the foreground signal strength: , where is the K-band map (with the ILC map subtracted) in units of mK, antenna temperature. This gives

After each minimization of , we update for K- through V-band according to

(3) |

where is the solution residual, and is a gain factor that we set to 0.5. For W-band we update by power-law extrapolation using the inferred Q,V-band spectral index. For some low signal-to-noise ratio pixels, we occasionally find . In a change from our first-year analysis, we now extrapolate the synchrotron spectrum from V- to W-band using . This change affects only 7% of pixels.

We iterate the minimization of and the update of eleven times. At the end of this cycle, the residual, , is less than 1% of the total signal . However, there are still several sources of potential error in the component decomposition, including:

(1) Zero-level uncertainty. As noted above, we use a plane-parallel slab model to assign the sky map zero point at each frequency. We estimate the uncertainty in this convention by fitting the model separately in the northern and southern Galactic hemispheres. These separate fits change the zero levels by as much as K. When these differences are propagated into , the output maps, , change by %, %, and % for the free-free, synchrotron, and dust components, respectively, at low Galactic latitudes.

(2) Dust spectrum uncertainty. Changing by 0.2 changes the component maps by %, %, and % for the free-free, synchrotron, and dust components, respectively, at low Galactic latitudes.

(3) Dipole subtraction uncertainty. A 0.5% dipole error would result in a systematic 0.5% gain error in all bands. The dominant effect would be to rescale each component map by 0.5%.

(4) CMB signal subtraction uncertainty. Errors in the ILC estimate of the CMB signal will produce errors in the corrected Galactic data, , used in the MEM analysis. We quantify ILC errors in §5.2, but note here that they are small compared to the total Galactic signal at low latitudes. However, they plausibly dominate the total uncertainty at high Galactic latitudes, where the \map data are (fortunately) dominated by CMB signal. Fractional uncertainties of 100% are not unlikely in the foreground model at high latitudes.

(5) While the total model residuals are small, there is still potentially significant uncertainty in the individual foreground components. The program outlined above produces three component maps, , , and , and the synchrotron spectrum model . To illustrate the degeneracy between these outputs, we consider a simplified single-pixel model of the form

(4) |

where the angle brackets indicate a full-sky average, and we explicitly evaluate the average MEM functional,

(5) |

for selected pairs of parameters, while marginalizing over the rest. Contour plots of are shown in Figure 4. For the panels that explore the shape of the synchrotron spectrum, , we parameterize it as a power law with a steepening parameter, and evaluate as a function of , while marginalizing over . For the rest, we iteratively update as per equation (3). For the most part, the output parameters are only weakly correlated; the most notable degeneracy is between the free-free amplitude and the synchrotron amplitude. \map data tightly constrain the sum of the two, but their difference is determined by the relative amplitude of the prior estimates for these two components, and by the initial synchrotron spectrum model. The synchrotron spectrum is found with modest significance to be steepening with increasing frequency. However, the dust index, (not shown), is poorly constrained; thus our conclusion in the first-year analysis that was not well founded upon further investigation.

Figure 5 shows the three input prior maps, , and the corresponding output component maps, , obtained from the three-year data. These maps are available on LAMBDA as part of the three-year data release. The maps are displayed using a logarithmic color stretch to highlight a range of intensity levels. The morphology and amplitude of the thermal dust emission are well predicted by the prior (FDS) dust map (see also §5.3). The free-free emission is generally over-predicted by the prior discussed above, especially in regions of high extinction. But even in regions where the extinction is low, we find the mean free-free to H ratio at K-band to be closer to 8 K R than the value of 11.4 K R assumed in generating the prior. Moreover, we find considerable variation in this ratio (a factor of 2) depending on location.

The most notable discrepancy between prior and output maps is seen in the synchrotron emission. Specifically, the K-band signal has a much more extended Galactic longitude distribution than does the 408 MHz emission, and it is remarkably well correlated with the thermal dust emission. Is this K-band non-thermal component due to anomalous dust emission or to mostly flat-spectrum synchrotron emission that dominates at microwave frequencies and is well correlated with dusty active star-forming regions? We cannot answer this question with \map data alone because the frequency range of \map does not extend low enough to see the predicted rollover in the low-frequency anomalous dust spectrum. However, we note the following points.

In Figure 6 we show the mean spectra of the three Galactic emission components observed by \map, in addition to the sum of the three. For comparison, we also show the signal observed at 408 MHz and we infer the signals at 10 and 15 GHz based on correlation analyses of the Tenerife CMB data by de Oliveira-Costa et al. (2004c) and de Oliveira-Costa et al. (1999). The curves show the mean signal in the range , as computed from the output MEM component maps: blue is dust, green is free-free, red is the non-thermal signal, and black is the sum of the three. The total intensity of the 408 MHz emission is remarkably well matched to a simple power-law extrapolation of the total \map signal measured from K-band to V-band. The spectral index of the dashed black curve is 2.65 between 408 MHz and K-band, and 2.69 between 408 MHz and Q-band. If we interpret the non-thermal emission as synchrotron, the implied spectral index between 408 MHz and the K-band non-thermal component is 2.73.

As noted in §4.4, de Oliveira-Costa et al. (2004) correlate the Tenerife 10 and 15 GHz data with the first-year \map non-thermal MEM map. We infer the non-thermal signal at 10 and 15 GHz (shown in red) by scaling the three-year non-thermal map using the reported correlation coefficients. We find a frequency rollover consistent with de Oliveira-Costa et al. (2004). Using the same scaling method, we also infer the 10 and 15 GHz signals derived from correlations with H emission (green arrows) and the Haslam 408 MHz (synchrotron) map (grey points). The derived free-free emission is lower than the extrapolation of the free-free emission inferred from \map – which is not physically tenable. Similarly, the Haslam-correlated emission at 15 GHz is substantially lower than it is at either 10 or 23 GHz. Thus, the sum of all correlated components in the 10 and 15 GHz data requires a substantial dip in the spectrum of the total emission (solid black points), which is not hinted at by the Haslam or \map data (dashed black points). This may be a sign of substantial anomalous emission, but one must be cautious interpreting the spectra of correlation coefficients.

Page et al. (2006) have used the three-year \map polarization data to construct a novel decomposition of the intensity foreground signal at each \map frequency band. In brief, they predict a synchrotron polarization fraction from a model of the Galaxy’s magnetic field strength and electron density. This fraction is de-rated by an empirical factor to account for missing structure in the model, then a “high-polarization-fraction” component of the intensity signal is formed as , where is the polarization intensity at frequency , in mK, and is the model polarization fraction. After subtracting an estimate of the CMB and free-free signal from the intensity maps, the remaining non-thermal signal is attributed to a “low-polarization fraction” component, , where is the temperature map corrected for CMB and free-free emission. The high-fraction component has a morphology similar to the 408 MHz emission (Figure 9 in Page et al. (2006)), while low-fraction component has a very dust-like morphology (Figure 13 in Page et al. (2006)). While the accuracy of this decomposition depends on the model fraction, , the basic picture should inform future studies of anomalous emission.

## 5 Galactic Foreground Removal

The primary goal of foreground removal is to provide a clean map of the CMB for cosmological analysis; an improved understanding of foreground astrophysics is a secondary goal. Removal techniques typically rely on the fact that the foreground signals have quite different spatial and spectral distributions than the CMB. In this section we describe two-and-a-half approaches to foreground removal that use complementary information. The first is an update of the Internal Linear Combination (ILC) method we employed in the first-year analysis (Bennett et al., 2003c). The second is an updated approach to fitting Galactic emission templates to each \map frequency band map. The remaining strategy is to mask regions of the sky that are too contaminated to be reliably cleaned. Extragalactic sources are treated in §6. For our primary CMB results, we analyze the masked, template-subtracted maps, but for some low- applications we also analyze the ILC map as a consistency check.

### 5.1 Temperature Masks

Many regions of the sky are so strongly contaminated by foreground signals that reliable cleaning cannot be assured. These regions are masked for cosmological analysis, though the extent of the masking required depends on the type of analysis being done. Bennett et al. (2003c) defined a set of pixel masks based on the first-year K-band temperature map. Since these masks were based on high signal-to-noise ratio Galactic signal contours in the K-band data, we have not modified the diffuse emission masks for the three-year analysis.

In addition to diffuse Galactic emission, point sources also contaminate the \map data. A point source mask was constructed for the first-year analysis that included all of the sources from Stickel et al. (1994); sources with 22 GHz fluxes Jy from Hirabayashi et al. (2000); flat spectrum objects from Teräsranta et al. (2001); and sources from the blazar survey of Perlman et al. (1998) and Landt et al. (2001). The mask contained nearly 700 objects, including all 208 of the sources directly detected by \map in the first-year data. Each source was masked to a radius of . For the three-year analysis, we have supplemented the source mask with objects from the three-year \map source catalog discussed in §6.1. Of the newly detected sources, 81 were not included in the previous mask and have been added to the three-year mask. Weaker, undetected sources still contribute to the high- angular power spectrum. As discussed in §7.2, we fit and subtract a residual source contribution to the multi-frequency power spectrum data.

Figure 7 gives an overview of the microwave sky and indicates the extent of the various foreground masks. The yellow, salmon, and red shaded bands indicate the diffuse masks defined in Bennett et al. (2003c). The violet shading shows the “P06” polarization analysis mask described in Page et al. (2006). The small blue dots indicate point sources detected by \map (to alleviate crowding, the full source mask described above is not shown). In addition, some well-known sources and regions are specifically called out.

### 5.2 The Internal Linear Combination (ILC) Method

Linear combinations of the multi-frequency \map sky maps can be formed using coefficients that approximately cancel Galactic signals while preserving the CMB signal. This approach exploits the fact that the frequency spectrum of foreground emission is different from that of the CMB. The method is “internal” in that it relies only on \map data, so the calibration and systematic errors of other experiments do not enter. There are a number of ways the coefficients can be determined, some of which require only minimal assumptions about the nature of the foreground signals. In the first year \map papers we introduced a method in which the coefficients were determined by minimizing the variance of the resulting map subject to the constraint that the coefficients sum to unity, in order to preserve the CMB signal. We called the resulting map the “ILC” map. In this section we elaborate on the strengths and limitations of the ILC method and quantify the uncertainties in the ILC map.

Eriksen et al. (2004) have also analyzed the method as an approach to foreground removal. They devised an approach to variance minimization that employed a Lagrange multiplier to linearize the problem and dubbed the resulting map the “LILC” map, where the first L denotes Lagrange. They found their LILC map differed somewhat from the ILC map in certain regions of the sky. We have since verified that the two minimization methods produce identical results for a given set of inputs and that the differences were due to an ambiguity in the way the regions were defined in the original ILC description. Because the linearized algorithm is considerably faster than our original nonlinear minimization, we have adopted it in the present work.

#### 5.2.1 Uniform Foreground Spectra

In order to better understand how errors arise in the ILC map, we first consider a simple scenario in which instrument noise is negligible and the spectrum of the foreground emission is uniform across the sky, or within a defined region of the sky. In this case, a frequency map, , may be written as a superposition of a CMB term, , and a foreground term, , where describes the composite frequency spectrum of the foreground emission, and describes the spatial distribution, so that . A linear combination map has the form

(6) |

where we have imposed the constraint , and have defined .

Suppose we choose to determine the coefficients by minimizing the variance of . Then,

(7) | |||||

(8) | |||||

(9) |

where the angle brackets indicate an average over pixels, and we have defined the variance and covariance in terms of these averages. Note that this expression would still hold if we added an arbitrary constant to each frequency map, . The ILC variance will be minimized when

(10) |

Thus the coefficients that minimize give , and in the absence of noise, the corresponding ILC solution is

(11) |

with

(12) |

In this ideal case, the frequency maps combine in such a way as to maximize the cancellation between CMB signal and foreground signal, producing a biased CMB map with . We have tested this result with ideal simulations in which we generate 5 frequency maps, , which include a Galaxy signal with a constant spectrum, , and random realizations of CMB signal and instrument noise. We then generate ILC maps from each realization and compare the residual map, , to the bias prediction, . The results confirm that the above description is correct, and that instrument noise is not a significant concern in this situation. The level of the bias is typically 10 K in the Galactic plane.

#### 5.2.2 Non-uniform Foreground Spectra

To minimize the anti-correlation bias we should choose regions that minimize the covariance between the CMB and the foreground, . However, in the previous analysis we assumed that the spectra of the foreground signals were constant over the sky. In reality these will vary as the ratio of synchrotron, free-free, and dust emission varies across the sky (and as the intrinsic synchrotron and dust spectra vary). In this case, the bias analysis becomes more complex. Specifically, the foreground component at each frequency may be written as , and the ILC map takes the form

(13) |

where . The ILC variance then generalizes to

(14) |

Using the same reasoning that led to equation (10), we obtain the following result for the minimum variance solution

(15) |

This has the same interpretation as equation (10), in the sense that it relates the foreground variance to the CMB-foreground covariance. We can solve this equation for by noting that , so that

(16) |

Now define and , whereby

(17) |

which is the multi-frequency analog of equation (11). Once again though, the bias in the ILC solution is proportional to (minus) the CMB-foreground covariance.

We have tested this expression with simulations like the ones described above, except this time we employ a three-component Galaxy model with variable spectra, , based on the first-year MEM model. As we discuss in more detail below, the simulations verify that the output ILC map is biased, , with as given in equation (17). Unfortunately, we do not know and a priori, and it has proven difficult to relate this bias expression to the frequency band maps, , in a way the can be used to minimize the bias. As a result, we have primarily resorted to correcting the bias with Monte Carlo simulations, as we describe below.

Given the results above, and our previous experience with the ILC method, it is clear that one should subdivide the sky into regions selected by foreground spectra, in order to reduce bias prior to correcting it. We have carried out such a program for the three-year analysis and have found it very difficult to improve on the region selection made in the first-year analysis (Bennett et al., 2003c). Nonetheless, we have adopted a few changes: (1) we eliminate the Taurus A region, as it is too small to ensure a reliable CMB-foreground separation (Eriksen et al., 2004); (2) we add a new region to minimize dust residuals in the Galactic plane. This region is based on a color selection and encompasses the outer Galactic plane within the Kp2 cut. The new ILC region map is shown in Figure 8. The region designated 1, shown in red, replaces the old Taurus A region; the remaining regions are unchanged. This map is available on LAMBDA as part of the three-year data release.

For each region , we determine a set of band weights, , by minimizing the variance of the linear combination map in that region, subject to the constraint . There are two exceptions to note. The coefficients for region 0 were derived from a subset of the data in that region, specifically pixels inside the Kp2 cut with . The coefficients for region 1 were derived from a slightly larger region of data, specifically pixels inside the Kp2 cut with , that pass the color cut. To ensure uniformity, the band maps have been smoothed to a common resolution of FWHM. The coefficients are given in Table 5.

We form a full-sky map by combining the region maps, ; but to minimize edge effects, we blend the region maps as follows. We create a set of full-sky weight maps, such that for and otherwise. We smooth these maps (which contain only ones and zeros) with a smoothing kernel, to get smoothed weight maps, . The final full-sky map is then given by

(18) |

where the sum is over the twelve sky regions.

In order to obtain the final bias correction, we generate multi-region Monte Carlo simulations, using the variable spectrum, MEM-based Galaxy model as input. We evaluate the error, , for each realization and compute a bias from the mean error averaged over 100 realizations. The result, shown in Figure 8, is roughly 20-30 K in the Galactic plane, but substantially less off the plane. This map was used to correct the three-year ILC map, which is shown in the middle panel of Figure 9. This Figure also shows the first-year ILC map (top panel) for comparison. The difference between the two (bottom panel) is primarily due to the new bias correction, but a small quadrupole difference, due to the changes noted in Figure 3, is also visible.

Based on the Monte Carlo simulations carried out for this ILC study, we estimate that residual Galactic removal errors in the three-year ILC map are less than 5 K on angular scales greater than . But we caution that on smaller scales, there is significant structure in the bias correction map that is still uncertain. On larger scales, we believe the three-year ILC map provides a reliable estimate of the CMB signal, with negligible instrument noise, over the full sky. We analyze the low- multipole moments of this map in §7.3.

### 5.3 Foreground Template Subtraction

The ILC method discussed above produces a CMB map with complicated noise properties, while the MEM method discussed in §4.5 is primarily used to identify and separate foreground components from each other. For most cosmological analyses one must retain the well-defined noise properties of the \map frequency band maps. To achieve this we form low-noise model templates of each foreground emission component and fit them to the \map sky maps at each frequency. After subtracting the best-fit model we mask regions that cannot be reliably cleaned because of limitations in the template models. In this section we describe the model templates we use for synchrotron, free-free, and dust emission and we estimate the residual foreground uncertainties that remain after these templates have been fit and subtracted. The \map band maps are calibrated in thermodynamic temperature units; where appropriate, we convert Galactic signals to units of antenna temperature using the factors given in Table 1.

In our first-year model we used the Haslam 408 MHz map as a template for synchrotron emission. We now use the \map K- and Ka-band data to provide a synchrotron template, as described below. This is preferable because: (1) the intrinsic systematic measurement errors are smaller in the \map data than in the Haslam data, and (2) the non-uniform synchrotron spectrum produces morphological changes in the brightness as a function of frequency (Bennett et al., 2003c), so that the low frequency Haslam map is less reliable at tracing microwave synchrotron emission than the \map data.

There are two potential pitfalls associated with using the K- and Ka-band data for cleaning: (1) the data are somewhat noisy, and since the template subtraction will be common to all cleaned channels, there can be a noise bias introduced in the inferred angular power spectrum. (But note that we use separate templates for each year of data, so the correlation only acts across frequency bands within a single year.) (2) Since the K- and Ka- band data are contaminated with point sources, this signal could interfere with the primary goal of cleaning the diffuse emission. Using the fitting coefficients obtained below and the known noise properties of the K- and Ka-band data, we estimate the noise bias in the final power spectrum to be 5 K near the 1st acoustic peak (0.1% of the CMB signal), and even smaller at lower and higher multipoles. Further, assuming the point source model given in equation (43), and the fact that the template has been smoothed to an effective resolution of FWHM, we estimate that sources contribute 1 K to the power spectrum at in the template, and thus may be safely ignored. In the end, these pitfalls are not a source of concern for the three-year analysis.

The difference map , in thermodynamic units, cancels CMB signal while it contains a specific linear combination of synchrotron and free-free emission (and a minimal level of thermal dust emission). We use this map as the first template in the model. For the second template we use the full-sky H map compiled by Finkbeiner (2003) with a correction for dust extinction (Bennett et al., 2003c). This template independently traces free-free emission, allowing the model to produce an arbitrary ratio of synchrotron to free-free emission at a given frequency (the limitations of H as a proxy for free-free are discussed below). For dust emission, we adopt “Model 8” from the Finkbeiner et al. (1999) analysis of IRAS and COBE data, evaluated at 94 GHz (see §4.3). The full model has the form

(19) |

where are the fit coefficients for each template at frequency , and is the dust map. As discussed below, this model is simultaneously fit to the Q-, V-, and W-band maps, and we constrain the coefficients and to follow specified frequency spectra to minimize component degeneracy.

To clarify the physical interpretation of and we first note that may be rewritten in terms of synchrotron and free-free emission as

(20) |

where and are the synchrotron and free-free maps in antenna temperature at K-band, is the surviving fraction of emission component (synchrotron or free-free) in , and is the spectrum of component , in antenna temperature, relative to K-band. To a very good approximation, the spectrum of free-free emission is (§4.1), so that . For synchrotron emission, variations in the spectrum as a function of position will produce variations in . For spectral indices in the range we have . In the remainder of this section we assume and neglect the 2.5% error introduced by spectral index variations between K and Ka bands. Note that this value of is only used to estimate the level of synchrotron emission in the template ; we do not constrain the fit coefficients to follow a specified frequency spectrum.

By adding the H map to the model, we allow the synchrotron to free-free ratio to vary as a function of frequency, but we must be cognizant of potential errors introduced by the use of H as a proxy for the free-free emission. Nominally we have , where is the H intensity in Rayleighs and is the free-free to H ratio. At K-band is predicted to be 11.4 K R (Bennett et al., 2003c) but the actual ratio is both uncertain and dependent locally on extinction and reflection. These effects make the H proxy unacceptable in the Galactic plane, and force one to mask these regions for CMB analysis. Outside the masked regions the variations in are primarily due to residual extinction and to variations in the temperature of the emitting gas. Here higher fractional errors can be tolerated because the total free-free signal is fainter.

Due to the uncertainties in the free-free to H ratio, and the fact that contains a mixture of synchrotron and free-free emission, care must be taken to interpret the model correctly. Let the combined synchrotron and free-free emission in the data at frequency be

(21) |

where the terms are as defined above. The synchrotron and free-free terms in the model may be written as

(22) | |||||

(23) |

Comparing the synchrotron terms in equations (23) and (21) we can infer the mean synchrotron spectral index returned by the fit

(24) |

Comparing the free-free terms, and assuming is known, we can solve for the free-free to H ratio at K-band

(25) |

We fit the template model, equation (19), simultaneously to each of the eight Q- through W-band differencing assembly (DA) maps (the three-year maps smoothed to 1 resolution) by minimizing

(26) |

where is the \map sky map from DA (in thermodynamic units), is the mean noise variance per pixel for DA , and the second sum is over pixels outside the Kp2 sky cut. To regularize the model, we impose the following constraints on the fit coefficients: (1) all coefficients must be positive-definite, (2) the dust coefficients must follow a spectrum , and (3) the free-free coefficients for each DA must follow a free-free spectrum, which leads to the following form

(27) |

The synchrotron coefficients are fit separately for each differencing assembly. Given the 10 coefficients from the three-year fit, we subtract the model from each single-year DA map to produce a set of cleaned maps. In doing so, we form separate single-year maps of to maintain rigorously independent noise between separate years of data.

The fit coefficients are given in Table 4 along with derived values for and . To facilitate model subtraction, we tabulate values for and for each DA using the above constraints. Note that the FDS dust model, which predates \map by a few years, predicts the 94 GHz dust signal remarkably well. The synchrotron emission shows a steady steepening with increasing frequency, as seen in the first-year data (Bennett et al., 2003c). Also, the free-free to H ratio is seen to be 6.5 K R, which is roughly half of the 11.4 K R prediction. Taken together, this fit finds a remarkably low total Galactic foreground amplitude at V-band.

Figure 10 shows the three-year band maps before and after subtracting the above model. For comparison, the Figure also shows the same three-year maps after subtracting the first-year template-based model (Bennett et al., 2003c). In all panels an estimate of the CMB signal (the ILC map) has been subtracted to better show residual foreground errors. The main visible difference between the first-year and three-year residual maps is the synchrotron subtraction error in the first-year model due to the use of the Haslam 408 MHz map. This is especially visible in the region of the North Polar Spur and around the inner Galaxy. Note also the significant model errors visible inside the Kp2 sky cut. This is presumably caused by a combination of synchrotron spectral index variations and errors in the extinction correction applied to the H template. Indeed, errors of up to 30 K are also clearly visible in isolated regions outside the cut, especially in the vicinity of the Gum Nebula and the Ophiuchus complex. In §7.2 we compare the foreground signal to the CMB and assess the degree to which these residual errors contaminate the CMB power spectrum.

## 6 Extragalactic Foregrounds

### 6.1 Point Sources

Extragalactic point sources contaminate the \map anisotropy data and a few hundred of them are strong enough that they should be masked and discarded prior to undertaking any CMB analysis. In this section we describe a new direct search for sources in the three-year \map band maps. Based on this search, we update the source mask that was used in the first-year analysis. In §7.2 we describe our approach to fitting and subtracting residual sources in the data. Page et al. (2006) discuss the treatment of polarized sources.

For the first-year analysis, we constructed a catalog of sources surveyed at 4.85 GHz using the northern hemisphere GB6 catalog (Gregory et al., 1996) and the southern hemisphere PMN catalog (Griffith et al., 1994, 1995; Wright et al., 1994, 1996). The GB6 catalog covers the declination range to a flux limit of 18 mJy, while the PMN catalog covers to a flux limit between 20 and 72 mJy. Combined, these catalogs contain 119,619 sources, with 93,799 in the region . We have examined the three-year \map sky maps for evidence of these sources as follows: we bin the catalog by source brightness and, for each bin, we cull the corresponding sky map pixels that contain those sources. The data show a clear correlation between source strength and mean sky map temperature that disappears if the sky map pixels are randomized. The multi-frequency \map data suggest that the detected sources are primarily flat-spectrum, with .

In the first-year analysis, we produced a catalog of bright point sources in the \map sky maps, independent of their presence in external surveys. This process has been repeated with the three-year maps as follows. We filter the weighted maps, ( is the number of observations per pixel) in harmonic space by , (Tegmark & de Oliveira-Costa, 1998; Refregier et al., 2000), where is the transfer function of the \map beam response (Page et al., 2003a; Jarosik et al., 2006), is the CMB angular power spectrum, and is the noise power. Peaks that are in the filtered maps are fit in the unfiltered maps to a Gaussian profile plus a planar baseline. The Gaussian amplitude is converted to a source flux density using the conversion factors given in Table 5 of Page et al. (2003a). When a source is identified with confidence in any band, the flux densities for other bands are given if they are and the fit source width is within a factor of 2 of the true beam width. We cross-correlate detected sources with the GB6, PMN, and Kühr et al. (1981) catalogs to identify 5 GHz counterparts. If a 5 GHz source is within of the \map source position (the \map source position uncertainty is ) we tag the \map source. When more than one source lies within the cutoff radius the brightest one is assumed to be the \map counterpart.

The catalog of 323 sources obtained from the three-year maps is listed in Table 9. In the first-year catalog, source ID numbers were assigned on the basis of position (sorted by galactic longitude). Now, rather than assigning new numbers to the newly detected sources, we recommend that WMAP sources be referred to by their coordinates, e.g., WMAP J0006-0622. For reference, we give the first-year source ID in column 3 of Table 9. The 5 GHz IDs are given in the last column. The source count distribution, , obtained from the Q-band data is shown in Figure 11.

The first-year catalog contained 208 sources. Given the increased sensitivity in the three-year maps, the number of new sources detected is consistent with expectations based on differential source count models. By the same token, 6 sources from the first-year catalog were dropped from the three-year list (numbers 15, 31, 61, 96, 156, and 168). Simulations of the first-year catalog suggested that it contained false detections, so the number of dropped sources is consistent with expectations. Five of the 208 sources in the first-year catalog could not be identified with 5 GHz counterparts; now 6 out of 323 can not be. Of the 6 sources that dropped off the first-year catalog, source numbers 15 and 61 did not have 5 GHz identifications. The remaining 4 may have dropped off because of variable (declining) flux density.

Trushkin (2003) has compiled multifrequency radio spectra and high resolution radio maps of the sources in the first-year \map catalog. Reliable identifications are claimed for 205 of the 208 sources. Of the 203 sources with optical identifications, Trushkin (2003) finds 141 quasars, 29 galaxies, 19 active galactic nuclei, 19 BL Lac-type objects and one planetary nebula, IC418. Forty percent of the sources are identified as having flat and inverted radio spectra, 13% have GHz-peaked spectra, 8% are classical power-law sources, and 7% have a classical low frequency power-law combined with a flat or inverted spectrum component (like 3C84). Trushkin (2003) suggests that \map source number 116 is likely to be spurious and, for source 61 no radio component was found.

### 6.2 Sunyaev-Zeldovich (SZ) Effect

Hot gas in galaxy clusters produces secondary anisotropy in the CMB via the Sunyaev-Zeldovich effect: a systematic frequency shift of CMB photons produced by Compton scattering off hot electrons in the cluster gas. At frequencies less than 217 GHz, this produces a temperature decrement in the CMB along sight lines that pass through clusters. The effect is relatively small for \map due to its moderate angular resolution; nevertheless, nearby clusters are large enough to produce a measurable effect in the sky maps. In our first year analysis we reported detections of an SZ signal from the Coma cluster and from an aggregate sample of X-ray selected Abell clusters, the XBAC catalog. Since then, there have been numerous additional reports of SZ signal detection in the first-year \map data (Fosalba & Gaztañaga, 2004; Fosalba et al., 2003; Myers et al., 2004; Afshordi et al., 2004; Hernández-Monteagudo et al., 2004; Hernández-Monteagudo & Rubiño-Martín, 2004; Afshordi et al., 2005; Hansen et al., 2005; Atrio-Barandela & Muecket, 2006). In this section, we update our results for the Coma cluster and the XBAC catalog; in subsequent cosmological studies, we mask these clusters from the data.

The brightest SZ source is the Coma cluster. For SZ analysis we model its temperature profile with an isothermal -model of the form

(28) |

where is the angular distance from the core, located at ; Briel et al. (1992) give and . Using these values, we convolve the profile with the \map beam response at each frequency (Page et al., 2003a; Jarosik et al., 2006) and fit for by minimizing

(29) |

where is the temperature in pixel , is a temperature baseline, is the angular distance of pixel from the cluster core, and is the pixel-pixel covariance matrix of the sky map. We model as , where is the 2-point correlation function of the CMB and is the variance due to instrument noise in pixel .

For W-band we find a Coma decrement of mK, with a reduced of 1.05 for 196 degrees of freedom. In V-band the decrement is mK. Herbig et al. (1995) measured a Rayleigh-Jeans decrement of mK at 32 GHz. Given the SZ frequency spectrum, their result predicts mK at V-band and mK at W-band, consistent with what we observe.

The XBAC catalog of X-ray clusters produced by Ebeling et al. (1996) is an essentially complete flux-limited sample of 242 Abell clusters selected from the ROSAT All-Sky Survey data. We cross-correlate this catalog with the \map 94 GHz map. Treating the XBAC clusters as point sources, we estimate their flux density at 94 GHz using the expression from Refregier et al. (2000)

(30) |

where is the angular diameter distance to the cluster, is the gas mass fraction, and is the electron temperature. The overall normalization of this relation is uncertain due to our ignorance of the correct gas fraction and cluster virialization state. We fix . Extended clusters (more than one-third the extent of Coma) are omitted; the remaining fluxes are all Jy. A template map is constructed by convolving the clusters with the \map W-band beam response (Page et al., 2003a) and fitting the result to the \map 94 GHz map. We find a template normalization of for the three-year data (compared to in the first year data). Since the fluxes used to construct the template are positive, the negative scaling is consistent with a 2.5 SZ decrement.

CMB photons that travel to us through the plane of the Milky Way undergo an SZ distortion of , where is the Thomson scattering cross-section and is the electron-pressure-weighted path length through our Galaxy. Taking K cm and kpc we find ; thus our Galaxy does not significantly affect CMB photons. The SZ effect can safely be ignored as a diffuse contaminating foreground signal.

## 7 The Angular Power Spectrum

Full-sky maps provide the most compact record of CMB anisotropy without loss of information. They permit a wide variety of statistics to be computed from the data, the most fundamental of which is the CMB angular power spectrum. Indeed, if the temperature fluctuations are Gaussian distributed, with random phase, the angular power spectrum provides a complete description of the statistical properties of the CMB. While there have been a number of papers based on the first-year data that claim evidence of non-Gaussianity and/or non-random phases in the fluctuations (see below) it is clear that the fluctuations are not strongly deviant from Gaussian, random phase. Thus, the measured power spectrum does provide the primary point of contact between data and cosmological parameters. This section presents the angular power spectrum obtained from the first three years of \map observations.

A sky map defined over the full sky may be decomposed in spherical harmonics as

(31) |

with

(32) |

where is a unit direction vector. If the CMB anisotropy is Gaussian distributed with random phases, then each is an independent Gaussian deviate with , and

(33) |

where is the angular power spectrum and is the Kronecker symbol. is the mean variance per that would be observed by a hypothetical ensemble of observers distributed throughout the universe. The actual power spectrum observed in our sky, by a supposedly typical member of this ensemble, would be

(34) |

If we had noiseless CMB data over the full sky, equation (32) could be evaluated exactly and equation (34) would give an unbiased estimate of the true power spectrum, in the sense that when averaged over the ensemble. However, we do not have ideal data (see below), and even if we did, since we only measure modes per (per sky), the above estimate of the variance has an intrinsic uncertainty (or “cosmic variance”) of

(35) |

Beyond cosmic variance, there are two effects that preclude using equations (32) and (34) to estimate the power spectrum; (1) real CMB data contains noise and other sources of error that cause the quadratic expression in equation (34) to be biased; and (2) data near the Galactic plane must be masked because Galactic emission cannot be reliably cleaned there. Masking precludes the proper evaluation of the integral in equation (32), so other methods must be found to estimate and .

In the first-year analysis, we addressed problem (1) by adopting a “cross-power” estimator in which we replace by , where and denote data channels with uncorrelated noise (Hinshaw et al., 2003b), see also Polenta et al. (2005); Patanchon (2003); Tristram et al. (2004). This form removes the noise bias from the estimate of . We addressed problem (2) by using a “pseudo-” estimator (Peebles & Hauser, 1974; Hivon et al., 2002) that statistically corrects for the aliasing introduced by non-uniform pixel weights, of which the sky cut is an extreme case.

An alternative approach is to invoke a maximum likelihood estimator that employs Bayes’ theorem to relate the likelihood of the power spectrum to the likelihood of the data,

(36) |

Here is the likelihood of the model (the underlying power spectrum, ) given the data (usually a sky map, or its transform, the coefficients), is the likelihood of the data given the model, and is the prior probability of the power spectrum. It is standard to assume that the data follow a multivariate Gaussian distribution (as predicted by most inflationary models, for example) in which case the likelihood takes the form

(37) |

where is the covariance matrix of the data, including contributions from both signal and noise, . Since the covariance matrix includes both terms, the confidence regions on deduced from automatically incorporate the uncertainty from both cosmic variance and instrument noise. In addition, if the data are restricted to lie outside a Galaxy cut, the confidence regions on will incorporate the uncertainty due to aliasing. Finally, it is also standard to assume a uniform prior distribution for the model, , though given sufficiently robust data, this choice will be unimportant. There are many technical challenges one must overcome when evaluating equation (37) and most of these have been well covered in the literature (Peebles, 1973; Peebles & Hauser, 1974; Hamilton, 1997a, b; Tegmark, 1997; Bond et al., 1998; Borrill, 1999; Oh et al., 1999). See Page et al. (2006) for an extension of this methodology to a simultaneous analysis of the temperature and polarization data.

Since the \map first-year data release, several authors have re-derived the angular power spectrum from the sky maps using a variety of other estimators. Some of these analyses cover the full range to which \map is sensitive, while others restrict attention to the low regime where certain technical issues matter most, and where hints of unusual behavior have been suggested. Of the analyses covering the full range, Fosalba & Szapudi (2004) employ a pixel-space estimator called SpICE that produces a spectrum that agrees remarkably well with the released \map spectrum. These authors do find small discrepancies at the very highest values (where the \map data is noise-dominated) but these differences do not significantly affect parameter determination. O’Dwyer et al. (2004) and Eriksen et al. (2004b) use a novel method based on Gibbs Sampling to undertake a full Bayesian analysis of the power spectrum. This method also produces results that are in good agreement with the \map spectrum. In addition, the method is able to generate a full conditional likelihood for the spectrum that allows one to rigorously evaluate cosmological models without resorting to likelihood approximations.

Efstathiou (2004b) has surveyed a variety of methods for estimating the power spectrum of large data sets such as \map and concluded that the best approach is a hybrid one which employs a maximum likelihood estimate at low and a pseudo- based estimate at high . We agree with this approach and adopt it for the three-year analysis, making the transition at . This gives results that agree well with the first-year spectrum for , and reasonably well for low-, as detailed below. In the original version of this paper, we chose to make the transition at , but a detailed study of the three-year spectrum by Eriksen et al. (2006) suggested that the transition needed to be at to avoid biasing cosmological parameter estimates, particularly the scalar spectral index, . Eriksen et al. (2006) and Spergel et al. (2006) (Appendix A) discuss the effect of this choice on the spectral index.

In the following sections we discuss the instrumental properties that are important to know for accurate power spectrum estimation, and assess the ability of the new Galactic foreground models to clean the power spectrum data. We then analyze the low and high- power spectrum in detail and compare the result with the previous \map spectrum and with other contemporary CMB data.

### 7.1 Instrumental properties

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 denote the auto or cross-power spectrum evaluated from two sky maps, and , where is a composite DA/year index ( denotes the DA, the year). Further, define the shorthand to denote a pair of composite indices, e.g., (Q1/yr-1,V2/yr-1), (Q1/yr-1,Q1/yr-2), etc. This spectrum will have the form

(38) |

where is the window function that describes the combined smoothing effects of the beam and the finite sky map pixel size. Here is the beam transfer function for DA , defined by Page et al. (2003a) and updated by Jarosik et al. (2006) [note that we reserve the term “beam window function” for ], and is the pixel transfer function supplied with the HEALPix package (Gorski et al., 2005). is the noise spectrum realized in this particular measurement. On average, the observed spectrum estimates the underlying power spectrum, ,

(39) |

where is the noise bias for differencing assembly , and the Kronecker symbol indicates that the noise is uncorrelated between differencing assemblies and between years of data. To estimate the underlying power spectrum on the sky, , the effects of the noise bias and beam convolution must be removed. The ability to determine these terms accurately is a critical element of any CMB experiment design.

In §7.1.1 we summarize the results of Jarosik et al. (2006) on the \map window functions and their uncertainties. We propagate these uncertainties through to the final Fisher matrix for the angular power spectrum. In §7.1.2 we present a model of the \map noise properties appropriate to power spectrum evaluation. For cross power spectra ( above), the noise bias term drops out of the signal estimate if the noise between the two DAs (or between years of data from a single DA) is uncorrelated. However, the noise bias term still enters into the error estimate, even for cross spectra. Therefore, the noise model is used primarily in error propagation.

#### 7.1.1 Window Functions

The instrument beam response was mapped in flight using observations of the planet Jupiter (Jarosik et al., 2006; Page et al., 2003a). The beam widths, measured in flight, range from (FWHM) at K band down to in some of the W band channels. The signal-to-noise ratio is such that the response, relative to the peak of the beam, is measured to approximately dB in W band. As part of the three-year analysis, we have produced a physical model of the A-side optics based on simultaneous fits to all 10 A-side beam pattern measurements (Jarosik et al., 2006). We use this model to augment the beam response data at very low signal-to-noise ratio ( to dB, depending on frequency band) which in turn allows us to better determine the total solid angle and window function of each A-side beam. For the B-side response we scale the A-side model by fitting it to the B-side Jupiter measurements in the high signal-to-noise regime. We then form similar hybrid response maps by augmenting the B-side data with the scaled model in the low signal-to-noise regime. (We plan to update the A-side model and produce an independent B-side model in the near future.) These hybrid beam response maps are available on LAMBDA as part of the \map three-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 for each DA (Jarosik et al., 2006). These transfer functions are also provided with the three-year data release.

The constraints imposed by the physical optics model have allowed us to extend the beam analysis out to much larger radii than was possible with the first-year analysis. As a result, we have determined that the beam solid angles were being systematically underestimated by 1%. Since we normalize the transfer functions to 1 at this meant that the normalized functions were being systematically over-estimated by a comparable amount for . Consequently, as we discuss in §7.6, the final three-year power spectrum is 1-2% higher than the first-year spectrum for . We believe this new determination of the beam response is more accurate than that given in the first-year analysis, but until we can complete the B-side model analysis, we have held the window function uncertainties at the first-year level (Page et al., 2003a). The first-year uncertainties are large enough to encompass the differences between the previous and current estimates of . The propagation of random beam errors through to the power spectrum Fisher matrix is discussed in Appendix A.

An additional source of systematic error in our treatment of the beam response arises from non-circularity of the main beam (Wu et al., 2001; Souradeep & Ratra, 2001; Mitra et al., 2004). The effects of non-circularity are mitigated by \map’s scan strategy, which causes most sky pixels to be observed over a wide range of azimuth angles. However, residual asymmetry in the effective beam response on the sky will, in general, bias estimates of the power spectrum at high . In Appendix B we compute the bias induced in a cross-power spectrum due to residual beam asymmetry. In principle the formalism is exact, but in order to make the calculation tractable, we make two approximations: (1) that \map’s scan coverage is independent of ecliptic longitude, and (2) that the azimuthal structure in the \map beams is adequately described by modes up to azimuthal quantum number . We define the bias as the ratio of t