# Three-Year Wilkinson Microwave Anisotropy Probe (Wmap)^{1}^{1}affiliation: Wmap 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:
Beam Profiles, Data Processing, Radiometer Characterization and Systematic Error Limits

###### Abstract

The WMAP satellite has completed 3 years of observations of the cosmic microwave background radiation. The 3-year data products include several sets of full sky maps of the Stokes I, Q and U parameters in 5 frequency bands, spanning to , and supporting items, such as beam window functions and noise covariance matrices. The processing used to produce the current sky maps and supporting products represents a significant advancement over the first year analysis, and is described herein. Improvements to the pointing reconstruction, radiometer gain modeling, window function determination and radiometer spectral noise parametrization are presented. A detailed description of the updated data processing that produces maximum likelihood sky map estimates is presented, along with the methods used to produce reduced resolution maps and corresponding noise covariance matrices. Finally two methods used to evaluate the noise of the full resolution sky maps are presented along with several representative year-to-year null tests, demonstrating that sky maps produced from data from different observational epochs are consistent.

^{†}

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

## 1 Introduction

The Wilkinson Microwave Anisotropy Probe (WMAP) is a Medium-class Explorer mission designed to produce full sky maps of the cosmic microwave background (CMB) radiation in five frequency bands centered at , , , and GHz. WMAP was launched on 2001 June 30 and began taking survey data on 2001 August 10 with its 20 high electron mobility transistor (HEMT) based radiometers. A suite of papers (Bennett et al., 2003b; Jarosik et al., 2003b; Page et al., 2003; Barnes et al., 2003; Hinshaw et al., 2003a; Bennett et al., 2003a; 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., 2003; Bennett et al., 2003c; Page et al., 2003; Barnes et al., 2002; Jarosik et al., 2003a; Nolta et al., 2003) describing the design of the observatory and the results of the first year’s observations have been previously published. Analysis of the first three years of WMAP data has now been completed and the results are presented in companion papers (Page et al., 2006; Hinshaw et al., 2006; Spergel et al., 2006).

One of the major design goals of WMAP was the careful control of systematic errors to allow the production of high quality maps of the microwave sky. Systematic errors in the first year results were analyzed in several papers accompanying the data release. Analyses of the beam shapes, beam window functions and associated errors were presented in Page et al. (2003). Details of the data processing and related errors were discussed in Hinshaw et al. (2003a). Radiometer performance and systematic errors were analyzed in Jarosik et al. (2003b) and Hinshaw et al. (2003a) while errors related to sidelobe pickup were explored in Barnes et al. (2003).

This paper updates and extends the previous analyses through the use of three years of observational data and addresses additional issues, such as data processing specific to the production of polarization maps. The overall instrument performance and improvements to the instrument modeling are described in § 2. Included in this section are updates to the radiometer gain models and beam window function analysis. In § 3 changes and additions to the data processing procedures, including the technique used to produce the maximum likelihood sky maps for Stokes I, Q and U parameters, and the evaluation of the pixel-pixel inverse noise matrix are described. In § 4 tests on the sky maps, including year to year comparisons and evaluation of map noise levels, are discussed.

## 2 Instrument Performance and Modeling

### 2.1 Observatory Status and Observing Efficiency

The WMAP spacecraft and instrument continued to perform normally throughout its second and third years of science data collection, spanning 2002 August 10 through 2004 August 9. Two station-keeping maneuvers were performed during each of the second and third years of WMAP observations. Each of these maneuvers resulted in the loss of approximately 10 h of science data, approximately 0.5 h for the maneuver itself and the remainder for the recovery from the thermal disturbance to the observatory. The observatory also entered a safe-hold mode on 2003 10 Aug, believed to be the result of a cosmic ray event. This occurrence resulted in a loss of approximately 60 h of data, caused both by the time spent in safe hold mode itself, and the time required for observatory temperatures to re-stabilize after the associated thermal disturbance. These periods, plus a short ( h) data loss due to a data recorder overflow caused by ground station difficulties that occurred on 2003 May 25, resulted in the loss of a total of 103 h of science data, resulting in an overall observing efficiency of for WMAP’s second and third years of observations. A tabulation of the times of the data excluded from processing by the aforementioned events can be found in Limon et al. (2006).

During the second and third years of operation, 8 sudden jumps were observed in the outputs of the radiometers. These events, distributed among 5 of the 20 radiometers comprising WMAP, are believed to be the result of sudden releases of thermally induced mechanical stresses that slightly alter the properties of some radiometer components (Jarosik et al., 2003b). During the first year of operation, 21 such events were observed. It is believed that the reduction in the number of events is the result of a more stable thermal environment during the second and third years of operation. Although the events are of short duration, the data processing used to produce the sky maps requires that a 2-4 h segment of data for the affected radiometer be excluded from the final processing to avoid artifacts in the sky maps. All 8 events taken together caused a loss of of the science data from the second and third year of observations.

### 2.2 Instrument Thermal Environment

Providing a stable thermal environment for the WMAP instrument is a key component in the production of a high quality data set with well defined systematic error limits. Of particular importance are temperature variations synchronous with the 129 s spin period of the observatory. The methods used to design and predict the thermal environment of the instrument and how they affect the instrument performance are described in Bennett et al. (2003c) and Jarosik et al. (2003a). On orbit, 56 high-resolution platinum-resistance thermometers monitor the temperature of key instrument components. These data provide verification that the actual thermal stability meets design specifications and allow for estimates of spurious signal levels resulting from radiometer thermal fluctuations. Analysis of the flight thermal data indicates the expected annual temperature modulation arising from the eccentricity of WMAP’s orbit and the asymptotic gradual warming of the instrument components as the thermal blankets degrade from ultraviolet exposure. Values characterizing the annual temperature modulation and gradual warming of the instrument components are presented in Table 1. The slow warming of the instrument components has increased the radiometer noise levels by based on analysis of the noise in the sky maps as presented in § 4.3. Projections indicate that the WMAP observatory should be able to operate for years (with rare solar eclipse avoidance maneuvers). Plots of the instrument and observatory temperature histories are presented in Limon et al. (2006).

Significant spin-synchronous temperature variations are detected in 10 of the 13 thermometers attached to the thermal reflector system (TRS) with peak-to-peak variations ranging from . Similar variations were seen during the first year of observations (Page et al., 2003; Limon et al., 2003) and were identified as the result of impulsive heating by solar radiation scattered off the rim of the solar array. Figure 1 displays the largest of these signals for each year of observations. The measured temperature profile as a function of Solar azimuth with respect to the observatory follows the same pattern for all 3 years of data. The similarity of the temperature profiles indicates that no significant degradation has occurred in the solar screening properties of the solar arrays or insulation blankets. The predicted microwave emissivity of the reflector surface () yields an estimated spin synchronous peak-to-peak spurious signal of in the time ordered data (TOD).

No spin-synchronous temperature variations were detected in the
Focal Plane Assembly (FPA), which contains the cold radiometer components, the
Receiver Box (RXB), which contains the warm radiometer components and phase switch drivers,
the Analog Electronics Unit (AEU) containing the data collection electronics, or
the Power Distribution Unit (PDU) which contains the power supplies for the HEMT amplifiers.
The limits for such signals are presented in Table 1. The bounds
are limited by thermometer readout noise and are similar to the year-1 values.
Limits on possible spurious
signals arising from thermal fluctuations of the radiometers are obtained by
combining the on-orbit temperature variations measured during the second and
third years of WMAP observations with instrument thermal
susceptibility coefficients based on preflight measurements (Jarosik et al., 2003b)
and on-orbit measurements (Hinshaw et al., 2003a). Artifacts in the
TOD resulting from thermal fluctuations
of radiometer components are found to be less than for all 20 radiometers, consistent
with the year-1 results. *Given the low level of the expected spin synchronous signal, no corrections
for spin synchronous effects have been applied to the time ordered data in the 3-year analysis.*

### 2.3 Pointing Determination

The boresight directions of the WMAP beams relative to the observatory are measured by relating radiometric measurements of Jupiter to the attitude of the observatory as determined by two onboard star trackers (Hinshaw et al., 2003a). WMAP observes Jupiter days per year in two day “Jupiter seasons”. To test pointing stability, boresight directions are computed separately for each Jupiter season, and differences between the individual seasonal determinations are noted. During the first year of WMAP observations, the azimuthal beam positions found for the first two Jupiter seasons agreed to better than , but the elevation positions differed by . This small difference was consistent with expected error in the spacecraft quaternions and was treated as part of the error budget rather than being actively corrected.

The second year of WMAP operation provided two additional Jupiter observing seasons, 2002 Nov 01 - 2002 Dec 26 and 2003 Mar 13 - 2003 May 08. With the addition of the third season, it was found that, while azimuthal positions were stable, the apparent elevations of the beams relative to the spin axis of the satellite now differed by from the positions originally computed from the first Jupiter season. Jupiter observations from season 4 confirmed the systematic change in elevation. By comparing the attitude measurements of the two star trackers and the radiometric Jupiter observations, the systematic variation was traced to a temperature dependent flexure occurring in the spacecraft structure, which altered the pointing of the star trackers relative to the microwave telescope’s beams.

Processing of the 3-year WMAP data assumes a linear
dependence of apparent star tracker elevation with spacecraft temperature, and
corrects the spacecraft quaternions appropriately.
A new set of line-of-sight vectors
describing the telescope beam positions for use with the updated quaternions
is available with the data release^{1}^{1}1 All data products in the current release are
available through http://lambda.gsfc.nasa.gov.
The resulting pointing corrections are small compared to the size of the beams and therefore
produce negligible changes to signals in the maps, but will slightly change the noise
patterns since observations near pixel borders may have moved into adjacent pixels relative to the first year analysis.
Based on the data from the six Jupiter seasons, residual pointing errors after
application of the corrections are estimated to be . A detailed description of the
modeling of the pointing correction is contained in Limon et al. (2006).

### 2.4 Radiometer Gain Model

The year-1 analysis of the WMAP data employed a model that related small changes in the radiometer gains (typically 1 %) to values of radiometer housekeeping data based on a physical model of the radiometer performance (Jarosik et al., 2003b). Each detector was modeled separately and used the measured values of the FPA temperature, , the radio frequency (RF) bias powers on the detectors, , and three parameters determined by fitting the model to the hourly gain measurements, obtained from the CMB dipole signal. The form of the model used in the year-1 analysis was

(1) |

The parameter characterizes compression (non-linear response) of the amplifiers and detectors, characterizes the variation of the noise temperature of the input amplifiers with physical temperature, and normalizes the overall gain. An example of the gain model performance fit to the hourly dipole-based gain measurements can be found in Jarosik et al. (2003b).

When the gain models fit to the first year observation data were extended to the
3-year data significant deviations became evident. These deviations became apparent
as the input housekeeping data to the model ( and ) spanned
a larger range than was used to fit the model, due to the slow warming of the
observatory. Re-fitting the original gain model parameters ( and
) greatly improved the agreement between the model and dipole derived gain measurements,
but significant errors still were evident. An additional
term was added to the model to account for
the *variation* of the compression of the amplifiers and
square law diode detectors with their physical temperature by allowing
the parameter characterizing the nonlinearity to have a weak temperature dependence.
The resulting form of the gain model is

(2) |

where is the temperature (in Kelvin) of the RXB thermometer closest to the detector and is an additional parameter to be fit. Figure 2 shows the dipole derived gain measurements for the V223 detector and the improved gain model for this detector, fit to the entire 3 years of data. Also shown in this Figure is the original gain model, using model parameters fit to year-1 data only, but applied to the entire 3 years of data . The new model significantly improves the fit over the entire 3 year period. The first-year analysis implemented a small, time dependent weighting of the TOD using the denominator of equation (1) as a measure of the input referenced noise temperature of the HEMT amplifiers. The values of and are strongly coupled in the current gain model and are not independently determined with high accuracy. It is therefore not possible to use the denominator of equation (2) as a measure of the radiometer noise levels. The three year processing therefore assumes uniform radiometer noise within each year. We maintain the year-1 estimate of the absolute sky map calibration uncertainty of .

#### 2.4.1 The Gain Model’s Effect on Measurement of the Quadrupole

Comparison of the year-1 and 3-year gain models indicated that many of the year-1 gain models displayed small errors in the predicted gain roughly linear in time, with an error on the order of .

For each DA^{2}^{2}2 A Differencing Assembly (DA) is a pair of differential radiometers connected to
the two linear polarizations of a set of telescope feed horns., difference
maps were formed by subtracting maps processed with
the original gain model from those processed with the
improved gain model. These maps displayed a feature with a several micro-Kelvin quadrupolar component, an example of
which is presented in Figure 3.
This signal was found to have similar morphologies in many of the maps with its
amplitude correlated to the difference between the mean slopes
of original and improved gain models.
No other multipoles show significant correlation with the mean slopes of the gain models.
(Sinusoidal differences between the gain models with an annual period were also
examined and showed no correlations with any multipole of the difference maps.)
These small errors in the gain model resulted in errors in subtraction of the
dipole signal.
When processed through the map-making pipeline this residual dipole signal yielded a
systematic feature with a quadrupolar component common to
many of the sky maps.
Since the form of the residual gain error was common to many of the radiometers the
resultant spurious quadrupolar features were aligned, resulting in a biased
measurement of the quadrupole moments in both the auto and cross power
spectra (Hinshaw et al., 2003b). The alignment of this spurious feature was such that it partially canceled
the true sky quadrupole signal. Correction of this systematic error accounts for a part of the small increase in the
reported value of the quadrupole moment, , from in the year-1
release (Bennett et al., 2003b) to , the largest change being the result of the change of the estimator use
to evaluate the quadrupole (Hinshaw et al., 2006). It should be noted that the current value is still far smaller than the mean expected value of
obtained from the best fit -dominated cold dark matter models (Hinshaw et al., 2006).

Estimates of residual quadrupolar contamination remaining in the maps processed with the improved gain model are obtained by multiplying a coefficient relating the amplitude of the spurious quadrupolar signal to the slope error in the gain model, , by an estimate of the residual slope error in the improved gain model, . The value of was obtained by a linear fit to the data contained in Figure 4 constrained to pass through the origin. The residual slope error in the gain model, , was estimated as the sum of two terms: a fit to the residuals of the gain model to the dipole based gain measurements, and a term to allow for a possible bias in the slope of the dipole based gain determinations. The value of this second term was obtained from simulations of the complete calibration algorithm including effects of far sidelobe pickup of Galactic and dipole signals. These simulations indicate that the hourly dipole based gain measurements have slope biases ranging from to . The magnitudes of these biases were summed with the magnitudes of the residual errors between the dipole based gain determinations and the gain model to yield an estimate . Values of are all less than . The resulting estimates of the magnitude of possible spurious quadrupolar signals in the 3-year maps are plotted in red in Figure 4 . An estimate of spurious contributions to is obtained by calculating the quantity

(3) |

where the average is over the Ka-W band DAs, and the nominal value for the quadrupole moment is taken as . The resultant estimate of the uncertainty on the quadrupole moment arising from gain model errors, , is .

### 2.5 Time Ordered Data Power Spectra

The outputs of the WMAP radiometers exhibit excess power at low frequencies characterized by a “1/f knee frequency”, as described in Jarosik et al. (2003b). Production of sky maps involves application of a filter to the TOD, as described in §3.4. These filters are derived from fits to the autocorrelation function of the TOD after removal of an estimated sky signal based on preliminary sky maps. The measured autocorrelation functions for each radiometer are parameterized as

(4) |

where is the time lag between data points in units of samples, the parameters , , , and are determined by fitting to the autocorrelation data, and is the time lag at which the fit crosses zero, typically s. These fits were performed on a year-by-year basis to allow for gradual changes in the radiometer noise characteristics, even though the radiometer noise properties are very stable. An example of this procedure is presented in the top panel of Figure 5, which presents the measured autocorrelation data and the parameterized fit for the year-3 observations of the W11 radiometer.

Tabulated values of the parameterized filter coefficients for all years and DAs are available with the data release.

### 2.6 Transmission Imbalance Measurement

The year-1 WMAP analysis included a measurement of a set of transmission imbalance factors, , characterizing the different transmission of sky signals from the A side and B side optics into the radiometers. The response, , of the WMAP radiometers to sky signals and may be written

(5) |

where the factors parameterize a departure from ideal differential radiometer performance, specifically the level of response to common mode input signals (Jarosik et al., 2003b). For an ideal differential radiometer , and the radiometer exhibits no response to common mode signal.

The transmission imbalance factors are determined by fitting the raw TOD to templates composed of pure differential and pure common mode signals originating from the CMB dipole anisotropy (Jarosik et al., 2003b). The year-1 analysis was performed using 232 days of data, and did not remove an estimate of Galactic emission and CMB anisotropy before performing the fits, since reliable sky maps were not available at the time the analysis was performed. The 3-year analysis improves on the previous analysis in two respects; 1) estimates of the Galactic and CMB anisotropy signals are removed from the TOD, leaving only the dipole signal, before fitting to the aforementioned templates, and 2) three years of TOD are used in the analysis. A fit to all three years of data was used to determine the values of the transmission imbalance factors, which are presented in Table 2. These values agree with those from the year-1 analysis to the degree expected, but are more accurate due to the improvements in the processing.

Radiometer | Radiometer | ||
---|---|---|---|

K11 | K12 | ||

Ka11 | Ka12 | ||

Q11 | Q12 | ||

Q21 | Q22 | ||

V11 | V12 | ||

V21 | V22 | ||

W11 | W12 | ||

W21 | W22 | ||

W31 | W32 | ||

W41 | W42 |

Note. – Measurement of the fractional input transmission imbalance, , obtained from the WMAP 3-year data. They were obtained via measurements of the radiometer responses to the common mode signal arising from the CMB dipole. All the values are small, nevertheless corrections for this effect have been included in the map making algorithm. The values of the uncertainties are estimated from the variance of measurements from 3 single year analyses.

### 2.7 Beam and Window Function Determination

Knowledge of the optical beam shapes is of critical importance to the scientific interpretation of the data. The shape of the measured power spectrum, and through that the values of the cosmological parameters, depends on the beam window functions. In turn, the beam window functions are determined from the shapes of the beam profiles. Uncertainties in the beam profile at the -20 dBi level (30 to 40 dB down from the peak, near the noise level of the measurements) can influence the window function at the 1% level. This sensitivity has motivated an extensive investigation of the beams.

The beam modeling we present here is based on six seasons of Jupiter data acquired over three years of observations. The modeling has been significantly enhanced and the approximations made for the first year analysis have been reexamined. With the new models, the beam profiles can be extrapolated below the noise level of the maps, thereby obviating the cutoff radius, , introduced in Page et al. (2003). It is found that on average the beam solid angles are systematically larger by 1% than the year-1 estimates, leading to a systematic decrease in the window functions at , leading in turn to a systematic increase in the power spectrum. For example in the range, the new combined V and W band window functions are 1.5% lower () than the same combination for year one, justifying the assumptions and treatment in the year-1 release. We retain nearly the same uncertainties on the window function that were given for year-1.

The co- and cross-polarization beam profiles for the A side are shown in Figure 6. These update previous versions of the figures by including the cross-polar response and data further into the tail of the beams. The figures are based on pre-flight data, though the general agreement between the projections and measurements of Jupiter is excellent. The predicted and measured cross-polar patterns are in general agreement.

The cross polar response of the main beam can be parameterized as a combination of a rotation of the linear polarization reference axis around the line of sight, and a coupling to the Stokes V component. The relative magnitudes of these two components depends on the phase of the cross-polar coupling. The two outputs from each feed horn are coupled to different radiometers, so WMAP is not directly sensitive to this phase. While neither of these components can produce an errant polarization signal from an unpolarized input signal, they do affect the polarization measurement. The rotation of the polarization reference axis from their design directions was measured to be less then during ground testing. The effect of the coupling to the Stokes V component is a reduction in the sensitivity to the Stokes Q and U components, presuming no Stokes V signal is present on the sky. The largest cross polar response is in W-band (-22 dB), corresponding to, at most, a 0.6% reduction in polarization sensitivity relative to the temperature calibration. At the levels indicated, neither of these effects alter the polarization measurements significantly. More details can be found in Page et al. (2006).

The beams of the V2 DA are used to illustrate the salient aspects of the beam analysis. Figure 7 shows the beam profile for year-1 and the three years combined. It also shows the accumulated solid angle as a function of angle from the main beam axis. Figure 8 shows the window function derived from the profile. A number of features are evident. 1) The profile is stable and is the same between years. 2) There are contributions to the solid angle at the -20 dBi level (the forward gain is 55 dBi), near the noise floor of the measurement. 3) The new beam transforms are slightly different from the old ones.

The beams are modeled using a code based on the DADRA (Diffraction Analysis of a Dual Reflector Antenna) physical optics routines (Rahmat-Samii et al., 1995). The deviations from an ideal beam shape can be explained as deformations of the reflectors. The shape of the primary is parametrized with a set of 122 Fourier modes on a rectangular grid and the shape of the secondary with 30 modes. Based on pre-flight measurements of the cold optics, the surface correlation length is cm, which is just resolved by our basis functions. A conjugate-gradient least squares method is used to solve for the shape of the primary and secondary reflectors. The parameters of the fit are the amplitudes of the modes. The model simultaneously solves for all ten beam profiles and takes into account the measured passbands. For simplicity we average over polarizations. Though the input can be thought of as 10 Jupiter maps with net pixels, the fitting is done in the time ordered data to obviate complications associated with map pixelization. The actual fit is a multistep process. At first just a few modes on the primary are considered. After an approximate solution is found, more modes on the primary are added, and finally the modes on the secondary are added. The solution is annealed and the resolution of the model is adjusted as the fitting progresses. So far, we have only run the model on the A side.

The fit method works because of the wide range of frequencies, the high signal to noise ratio, the angular resolution of the measurement, and the large sampling of the focal plane. Related methods generally employ some sort of sampling of the phase of the wave front— e.g., Out of Focus (OOF) holography (B. Nikolic & Hills, 2002)— which of course cannot be done with WMAP.

The residuals of the model are shown in Figure 9. These should be compared to Figure 3 of Page et al. (2003) which shows an earlier version of the model for the A side. For the best fit, the reduced (with roughly pixels). In other words, the model fits to near the noise of the measurement with just 152 parameters. One measure of the accuracy of the model comes from the comparison between the modeled and the measured beam solid angle as shown in Table 3. Excluding W2, the rms deviation is 1.6%; the W2 solid angle is predicted to be 4.8% smaller than is measured. This gives us confidence that we understand the beams. However, as can be seen in the figure, there are small residual features that the model does not capture. These features occur mostly in the steep parts of the profiles corresponding to the deformations with the greatest length scales. Nevertheless, the model is successful at unveiling faint large scale components of the beam. For example, a -26 dB Q-band lobe located below Q2 was predicted by the model before it was found in the data. The model’s effectiveness lies in the fact that it uses measurements with a high signal-to-noise ratio at all frequencies to find the shape of the reflector and thereby predict the beam shape many beamwidths off the beam axis at the noise floor of the measurements.

Figure 7 shows a comparison between the measured and modeled beam. The model captures the main features of the beam though there are some discrepancies at the -20 dB level. Most importantly, the model reveals that the beam profile continues to drop with increasing radius. This is simply a result of the correlation length of the surface deformations. With the new model there is no longer a need to impose a cutoff radius as was done in Page et al. (2003).

The window functions are computed from the symmetrized beam profiles following the Hermite method in Page et al. (2003), although the method is modified to take advantage of the modeling. The range of spacecraft orientations corresponding to the observational data contained in each sky map pixel depends on the ecliptic latitude of the pixel. For pixels near the ecliptic poles the observations occur nearly uniformly for all rotations of the spacecraft around each beam line-of-sight, effectively symmetrizing the beam. At lower ecliptic latitude the range of rotations is reduced. The use of the symmetrized beam profile approximation is therefore very good near the ecliptic poles, but worsens at lower ecliptic latitude. The least symmetrization occurs in the ecliptic plane. For beams in the ecliptic plane window function errors arising from use of the symmetrized beam approximation are less than 1% for in V-band and W-band. In Q-band the errors are somewhat larger due to the larger beam ellipticities, with window functions error less than for . The uncertainties in the window function determinations arising from incomplete beam symmetrization are included in the final window function uncertainties. More details regarding beam symmetrization can be found in Page et al. (2003).

Figure 7 shows one such fit for the symmetrized radial profile of the V2 A-side beam with a Hermite expansion of order 170. To remove some of the sensitivity of the window functions to the noise tail, hybrid beams are constructed. The hybrid method retains the high signal to noise observations of the main beam and replaces the measurement with the model in regions of low signal. For the A side, the hybrid beams are constructed by (a) scaling the model to measured peak height, (b) choosing a threshold level, (Table 3), and (c) using the measurement above and the scaled model below it. The hybridization is done with the full two dimensional profile in the TOD before the symmetrization. Since a complete model of the B side does not yet exist, the A-side profile is translated and rotated to the B side focal plane, scaled to the peak of the B-side data, and interpolated to replace B side observations in the TOD when they are less than . For both A and B sides the Hermite fits of order are made to the symmetrized beam profiles. Analytic and numerical models show that the assumption of symmetrization is sufficient except in Q band where a 4% correction in made for (Hinshaw et al., 2006).

The window functions are available with the data release. In general they follow those in Figure 4 of Page et al. (2003). At low l, the windows for polarization are negligibly different than those for temperature. At , the polarization windows can differ from the temperature windows by a few percent because the effective central frequencies for polarization are different from those for temperature. Thus, the effective beams for polarization are different from those for temperature. We do not take this effect into account yet, as it is negligible at the current levels of sensitivity.

Band | (dBi) | Below peak (dB) | |
---|---|---|---|

K | 17 | -30 | 0.999 |

Ka | 17 | -32 | 0.982 |

Q | 18 | -33 | 0.971 |

V | 19 | -36 | 0.999 |

W | 20 | -38 | 1.018 |

Note. – “Below peak” is the forward gain minus . It is the level in dB below the peak at which the measurement is replaced by the model. is the measured beam solid angle ( from observations of Jupiter) divided by the modeled beam solid angle.

The uncertainty in the window function is ascertained by changing the criterion for merging the model and the measurements, changing the cutoff in the Hermite fit, comparing the Hermite and Legendre polynomial methods, and propagating the formal errors from the fit. The adopted uncertainties, along with a comparison to year-1, are shown in Figure 10. The new uncertainties are very close to the year-1 uncertainties except at low in K through Q bands where the new uncertainties are larger. A typical window function has an uncertainty of 2-3%. These uncertainties add in quadrature in the cosmological analysis. We anticipate that the uncertainty can be reduced a factor of two after the development of the B side beam model and with more beam maps of Jupiter.

In addition to the work described above, two end-to-end consistency checks are performed. 1) We check to ensure that within any frequency band, all beams yield the same value for the temperature of Jupiter. This tells us that the beam solid angles are consistently computed for both A and B sides. We postpone a full reanalysis of the Jupiter calibration and recommend using the values in Page et al. (2003). 2) We make separate maps of the sky for the A and B sides and compute their power spectra. We then show that the ratio of the power spectra is unity to within the noise limits. This shows that the A and B side window functions are consistent.

### 2.8 Sidelobe Corrections

All radio telescopes exhibit some degree of sensitivity to sources outside of the main beam. Such pickup is a particular concern for CMB experiments since the measured signals are small and foreground emission can be large. The total beam of a telescope is traditionally described as the sum of two components, the main beam and the sidelobes. Although the demarcation between the two components is somewhat arbitrary, the two components sum to the total beam response of the telescope. The WMAP raw time ordered data, and the maps reconstructed from this data, represent the sky signal convolved with the total beam response of the telescope.

The exponential factor in the Hermite expansion of the symmetrized beam profiles (see § 2.7) and the order of the polynomial fit determine the maximum radius for which the Hermite expansion accurately models the beam pattern. Outside of this radius the Hermite expansion quickly decays. The radius at which the Hermite expansion decays, , provides a natural point at which to separate the main beam and sidelobe response: the window function derived from the Hermite expansion represents the beam for and the sidelobes account for the remainder of the beams. Ideally for , the sidelobe map would contain residuals between the true beam pattern and that parameterized by the Hermite expansion. However, since the main beam is modeled to levels below the noise floor of the measurements we have no way to determine these residuals. We therefore zero the area of the sidelobe maps for . The values of are presented in Table 4. These sidelobe maps, when convolved with various input signal maps, are used to implement corrections for the sidelobe pickup in the TOD as described in § 3.3. These maps differ from those contained in the first year data release in two aspects: 1) The region within a radius of around each line of sight direction has been set to zero, and 2) The region of the W-band maps derived from the preflight sidelobe mapping has been scaled down by , based on a re-evaluation of the calibration of those measurements. These sidelobe maps are available with the data release.

## 3 Data Processing

The 3-year WMAP sky maps were processed as three independent data sets, each containing 1 year of observational data. The maps were later combined to produce various data products, including a three year combined map. The processing used to produce sky maps differs in several aspects from that used to produce the year-1 maps. The most significant change is that the 3-year maps, comprising Stokes I, Q and U components, are the maximum likelihood estimates of the sky map given the radiometer “1/f” noise characteristics. The year-1 maps were unbiased representations of the sky signal, but were sub-optimal in terms of noise since they were produced from a pre-whitened TOD archive (Hinshaw et al., 2003a). Other significant processing changes include applying corrections for sidelobe pickup to the TOD for all DAs, inclusion of an estimated polarization signal in the fitting of the hourly baseline solution, and elimination of a weighting factor based on a radiometer noise estimate originating from the radiometer gain model. ( see § 2.4 )

The current data processing flow is outlined in Figure 11. Note that in the
*Make Sky Maps* block three sets of maps are produced for each year of data. The first of these are full sky maps
produced at HEALPix r9. WMAP utilizes the HEALPix^{3}^{3}3see http://healpix.jpl.nasa.gov
pixelization scheme and designates map resolution with the notation r4, r5, r9 and r10, corresponding to
HEALPix parameter values of 16, 32, 512 and 1024, and approximate pixel side dimensions of ,
, and
respectively.

These full sky maps suffer from small processing artifacts associated with observations
in which one of the telescope beams is in a region of strong Galactic emission
while the other is in a low Galactic emission region.
These artifacts arise from slight errors in the
radiometer gain determination and pixelization effects, resulting in errant signals in some
high Galactic latitudes regions used in CMB analyses.
The year-1 map processing eliminated this problem by only updating the value in the pixel with the high Galactic emission
for such observations in the iterative map making procedure (Hinshaw et al., 2003a). Adapting this
asymmetric masking to the conjugate gradient
processing used to produce the current maps is not straightforward. Instead, a separate set of r9 maps,
termed *spm* maps, is produced using a symmetric processing mask. These spm maps omit all
observations for which either beam falls into a
high Galactic emission region. The processing mask used to identify the high emission regions excludes of the pixels, and is
available with the data release.
These spm maps do not suffer the aforementioned problem, but contain no data in the high Galactic emission regions.
Data from the full sky maps are used to fill the unobserved regions of the spm maps.
Details of this process are presented in § 3.4.2.

The last set of maps listed in the *Make Sky Maps* block of Figure 11
is produced at r10 and contains only spm intensity maps. These maps are used to evaluate the high-
temperature power spectra and are produced at higher resolution to minimize pixelization effects. To speed computation these
maps were produced using only intensity information, so no corresponding high resolution polarization maps are available.

The implementation of the entire data processing pipeline has been exhaustively verified through numerous simulations of TOD with the non-idealities described in § 3.4.1, both with and without inclusion of instrument noise. In the case excluding instrument noise the pipeline has been shown to reproduce the simulated input maps to sub-nanoKelvin accuracy. Simulations with instrument noise have been shown to produce unbiased estimates of the input simulated sky maps. The subsequent discussion follows the notation of Hinshaw et al. (2003a) and updates the description of the data processing therein.

### 3.1 SMOC processing and Archive Generation

Processing of the data in the Science and Mission Operations Center (SMOC) and Archive Generation steps (see Figure 11) is virtually unchanged from that used in year-1. The only substantive change is that corrections for the thermally induced star tracker position errors ( § 2.3) are applied to the star tracker data in calculation of the pointing quaternions included in the raw TOD archive.

### 3.2 Calibration and Baseline Fitting

The calibration and baseline fitting procedures, used to calibrate the radiometric data based on the CMB dipole signal, are similar to those used in the year-1 processing (Hinshaw et al., 2003a). This is an iterative process in which both calibration data and approximate sky map solutions are generated. Hourly estimates of the radiometer gain are made by fitting each 1 hour segment of the TOD, corrected for estimated sky signals (excluding the CMB dipole), to the sum of a template dipole signal and a baseline term that accounts for very slow drifts in the radiometer outputs. The template dipole signal is derived from WMAP velocity and attitude information and the measured barycenter dipole from the first-year WMAP results. Improved sky maps are then generated using the updated baseline and hourly gain solutions, which in turn are used to produce improved baseline and gain solutions. The 3-year processing is an improvement over the 1-year processing in that it corrects the TOD for estimated sky signals based on the Stokes I, Q and U components before fitting for the gain and baseline, whereas in the year-1 processing only the Stokes I signal was removed from the TOD before the baseline and gain solution were fit.

The WMAP scan strategy ensures that the Stokes I sky signal components average
very nearly to zero over a 1 hour precession period. However, the orientation of the polarization axis
of the feed horns and scan pattern can transform certain
sky *polarization* signals into TOD signals with very long periods, extending to many hours or days.
The baseline fitting routine forces the
mean of the signal corrected TOD to zero on time scales of one hour and longer. Removing the polarization signal before
performing the baseline fitting ensures that the polarization signals contained in the TOD remain unbiased.

After the hourly gain and baseline solutions are converged they are used to fit the parameters of the improved gain model described in § 2.4.

### 3.3 Calibrated TOD Archive Production

The production of a final calibrated TOD archive involves two major steps: first the gain and baseline are applied to the uncalibrated TOD using the parameters determined in the calibration and baseline fitting procedure. Two corrections related to sidelobe response are then applied to this initial calibrated archive.

The first sidelobe correction simply removes an estimate of the signal arising from sidelobe pickup. This signal is calculated by convolving the Stokes I sidelobe response maps (§ 2.8) for each DA with preliminary Stokes I sky maps using the technique of Wandelt & Górski (2001). The input sky maps contain three components: 1) An estimated CMB + Galactic signal obtained from a preliminary set of maps produced without sidelobe corrections, 2) a fixed barycenter CMB dipole signal, and 3) an annually modulated CMB dipole signal from Earth’s motion relative to the solar system barycenter. These corrections apply to the Stokes I signal only and are therefore applied equally to the data from both radiometers comprising each DA. This correction is applied to the TOD on a sample-by-sample basis. The effects of polarized sidelobe pickup are small (Barnes et al., 2003) and are not corrected in this data release.

The second sidelobe related correction compensates for a small calibration error introduced in the initial calibration process by multiplying the sidelobe corrected TOD archive by an overall scale factor. The uncalibrated TOD archive used to fit the gain and baseline solution contains signal arising from both the main beam and the sidelobe response of the telescope optics to the true sky signal. However, the template dipole signal used to fit the gain is based on ideal pencil beams in the boresight directions of each telescope beam. This approximation can lead to biased gain measurements since it ignores the component of the TOD arising from sidelobe pickup of the sky signal. This effect has been studied through simulations in which simulated TOD, based on the sidelobe response maps, is input to the hourly gain fitting algorithm to determine the level of bias originating from the sidelobe signals. The bias is found to be relatively small and varies only slightly over the course of a year as the precession axis of the WMAP scan pattern sweeps across the sky. Table 4 lists the values of these recalibration factors. A detailed description of this correction is presented in Appendix A.

### 3.4 Maximum Likelihood Map Solution

As described in Hinshaw et al. (2003a), the WMAP time ordered data may be written as

(6) |

where is the mapping matrix that transforms a sky map, , into discrete samples, , and is the noise added to each sample by the radiometer. The map making problem consists of generating an estimated sky map given and the statistical properties of the noise. Taking the noise description as

(7) | |||||

(8) |

it is well known that the maximum likelihood estimate of the sky map, , is

(9) |

In the case of a single differential radiometer and sky maps comprising only Stokes I, is a data vector comprised of time ordered data elements, is a sky map of pixels, is an matrix and is an element matrix. The matrix is sparse, each row corresponding to one observation, and each column corresponding to a map pixel. For an ideal differential radiometer each row contains a value of in the column corresponding to the pointing of the A-side beam for a given observation, and a value of in the column corresponding to the B-side beam. Given these matrices it is possible to produce maps from the WMAP data set using equation (9). The second term on the right hand side of equation (9) may be evaluated directly, while multiplication by the first term, , may be performed using a conjugate gradient iterative technique. It is straightforward to generalize this technique to process polarization maps.

As described in Bennett et al. (2003c), the signal from each telescope beam (A-side and B-side) is separated into two orthogonal linear polarizations by orthomode transducers at the base of the feed horns. Each pair of feed horns is associated with two differential radiometers comprising the DA. One linear polarization from each feed horn is fed into one differential radiometer, while the orthogonal set of linear polarizations are fed to the second radiometer comprising the DA. Details of the alignment of the polarization axis and beam boresights were presented in Page et al. (2003).

The outputs of radiometers “1” and “2”, and , for each observation correspond to sky signals

(10) |

and

(11) |

Here , and are the Stokes I, Q and U sky signals at sky position . The polarization angle of radiometer “1” with respect to the sky reference direction for sky position is as described in Hinshaw et al. (2003a). Variables with subscript B are the corresponding values for the sky position observed by the B-side beam.

The mapping matrix is generalized to include polarization processing by expanding it to matrix where the number of rows has been doubled since there are two data values from each observation, and the number of columns has been increased to accommodate 3 maps, corresponding to the three Stokes parameters. Each row of the polarization mapping matrix has 6 non-zero elements, with in columns corresponding to the observed pixels (A-side and B-side for each radiometer ) in the map and , and , in appropriate locations in the columns corresponding to the and maps. Each observation is associated with 12 non-zero values of the mapping matrix that are distributed in two rows. The noise matrix is similarly expanded to account for the signals from both radiometers. The noise of the two radiometers is uncorrelated, so it may be described by the relations,

(12) | |||

(13) | |||

(14) | |||

(15) |

The full noise covariance matrix is simply a block diagonal combination of and .

The noise from different radiometers was verified as uncorrelated (eq.[13]) by evaluation of the cross correlation coefficient

(16) |

where and were obtained from the TOD by subtracting an estimated signal based on a combined 3 year final sky map. For all DAs, except K1, the measured value of the magnitude of (at zero time lag) was . Similar results were obtained from simulations of sky signal plus noise which were constructed to have zero noise cross correlations. The small values of anti-correlation observed arise from noise in the common sky map used to remove the sky signal from the TOD of the two radiometers. For K-band, incomplete removal of the intense Galactic signal due to small gain errors () lead to a measured value of , similar to values obtained from the simulation which also included gain uncertainties.

#### 3.4.1 Radiometer Non-Idealities

If the WMAP radiometers were ideal, solutions based on the mapping matrix described in the previous section would reproduce the maximum likelihood estimate of the sky signal. There are two known instrumental non-idealities that affect the maps and are treated in the analysis: bandpass mismatch and input transmission imbalance.

Bandpass Mismatch
The construction of the mapping matrix presented in §3.4 implicitly
assumed that the microwave frequency response of the two radiometers comprising each DA were identical.
The WMAP radiometers have slightly differing frequency responses (Jarosik et al., 2003b) which can cause signals with spectra
different from that of the calibration source (the CMB dipole) to be aliased into polarization maps (Barnes et al., 2003). This
*spurious* signal has the characteristic that it appears as a difference between the response of the two
radiometers comprising a DA, , but is not modulated with the polarization angles
and as a true polarization signal would be, since it only depends on the intensity
and spectral index of the source region.
This effect can be treated in the map making procedure by considering the signal from the radiometers
to originate from 4 source maps, the original 3 Stokes parameters plus a spurious map, . Including this extra term, the TOD is
then described by the relations

(17) |

and

(18) |

Note that for the combination of radiometers outputs corresponding to Stokes I, , the spurious term cancels, and for the combination corresponding to polarization signals, , it appears as a scalar map, independent of polarization angles and . Equation (9) is easily generalized to handle these extra terms by expanding the mapping matrix to , where the extra columns now correspond to the spurious signal map. Each row of this matrix now has 8 non-zero entries, 2 corresponding to each type of map, and .

Input Transmission Imbalance An ideal differential radiometer only responds to differential input signals and completely rejects common mode signals. Due to several effects, the WMAP radiometers exhibit a small response to common mode signals (Jarosik et al., 2003b). This response is characterized by a transmission imbalance factor, , with the response of the radiometer, , to input signals and given by equation 5. The values of for the 3-year data are given in Table 2. The effect can easily be incorporated into the maximum likelihood maps solution (see eq.[9]) by modification of the mapping matrix, . The TOD, including both the spurious map term and the transmission imbalance terms is given by

(19) | |||

(20) |

and

(21) | |||

(22) |

The corresponding modification to is that each term involving A-side data is multiplied by a factor and those involving B-side data are multiplied by a factor .

#### 3.4.2 Map Evaluation

The maximum likelihood map solution (eq.[9]) is evaluated in two steps. First the product , the “iteration 0” maps, are formed. The maps are then multiplied by the term using an iterative technique.

Generation of the Maps The calibrated TOD, , represents the entire sky signal, including the CMB dipole components. To minimize numerical errors a TOD signal corresponding to a nominal CMB dipole is removed from the calibrated TOD before evaluation of the maps. The dipole signal removed from the TOD for each radiometer includes the effects of loss imbalance based on the values given in Table 2. A CMB dipole amplitude of 3.3463 mK (thermodynamic) in the direction is used to calculate the barycentric component. A CMB monopole temperature of 2.725 K (thermodynamic) (Mather et al., 1999) is used to calculate the annually varying component due to the spacecraft’s motion about the barycenter. The kinematic quadrupole in not removed from the TOD in the generation of the maps.

The dipole-subtracted TOD from each radiometer must be multiplied by a filter consisting of the inverse of each radiometer’s noise correlation matrix. This is accomplished by forming the inverse noise correlation function,

(23) |

where is the parametrized noise correlation function as described in § 2.5. Since is constructed to be zero at lags greater than , also falls to zero on the same time scale, so values of the filter, , are also set to zero for lags at which was set to zero. An example of the steps involved in forming this filter function for the W11 radiometer is shown in Figure 5. Next, a constant is added to the values of for lags less than to force the mean of this portion of the filter to zero. This ensures that very low frequencies, those with periods longer than the filter extent, are filtered out. Finally, the filter is normalized by an overall multiplicative constant, , such that the value at is unity. Note that none of these adjustments to biases the resultant sky maps since the same form of is used in both terms in equation (9). Tabulated versions of these filters are contained in the data release.

The radiometer noise properties are treated as stationary for each year of data, making the inverse noise matrix, , circulant. The matrix multiplication, , may therefore be implemented through use of the convolution , and is implemented using standard Fourier techniques. Missing samples in the TOD, due either to gaps in the data or masking when producing spm maps, are zero-padded to properly preserve the time relationship between data on either side of the gaps.

Given the filtered TOD, the maps are evaluated by accumulating the product of the corresponding elements of and the TOD sample by sample. Cut sky maps are formed by simply replacing the rows of with zeros for observations where the data are masked. These maps correspond to the sky maps weighted by the inverse of the pixel-pixel noise correlation matrix, ,

(24) |

where represents sky maps of the three Stokes parameters, I, Q and U, and a map corresponding to the spurious signal described in § 3.4.1. These maps are available with the data release. Solving for the sky maps simply requires multiplication by the pixel-pixel noise covariance matrix, .

Final Sky Map Production The final step in the sky map production, multiplication of the maps by the pixel-pixel noise correlation matrix, , is effected through a conjugate gradient iterative technique. This method allows solution of the symmetric positive definite system

(25) |

for vector given vector and the ability to multiply an arbitrary vector by the matrix . In this application vector is a map and matrix corresponds to the inverse of the corresponding pixel-pixel noise correlation matrix, . Multiplication of this matrix by a vector representing a set of sky maps is a straightforward operation. A simulated time stream is generated from the input maps by stepping through the archive evaluating the 16 non-zero elements of the mapping matrix (contained in two rows) corresponding to each observation (the operation). The resultant time stream is then filtered and re-accumulated into a new map (the ) operation using the same method previously described to produce the maps. The conjugate gradient method then iteratively improves estimates of , with the level of convergence being characterized by the quantity given by the relation

(26) |

where is the vector-norm operator. The iterative processing was stopped when the condition was satisfied. \subsubsubsectionPreconditioning The rate of convergence of the conjugate gradient method is quite sensitive to the properties of the matrix , with nearly diagonal forms preferable. If another tractable matrix multiplication operation can be found which approximates multiplication by , the rate of convergence can be greatly increased. Consider the case where the matrix meets these criteria, multiplication of equation (25) by yields

(27) |

which has the same solution as equation (25), but the matrix product is nearly diagonal, speeding convergence of the solution. This technique was used in production of the final sky maps. The preconditioner was implemented by separating the input map into high resolution (r9) and low resolution (r4) components. The low resolution map was formed by re-binning the high resolution map with uniform weights. This low resolution map was then subtracted from the high resolution map, leaving only small angular scale information in the high resolution map. The low resolution component was multiplied by the inverse of pixel-pixel noise correlation matrix (§ 3.5) and the high resolution component multiplied by the reciprocal of the diagonal components of the full resolution noise correlation matrix. The two components were then recombined to produce the preconditioned map. The map production algorithms were carefully checked using numerous simulations to verify that they converged to the correct solutions. It was also verified that this solution was independent of the exact form of the preconditioner employed, although as expected the convergence rates did vary. Using the preconditioner as described, the r9 sky map sets, comprising Stokes I, Q and U components, and the spurious map S, converged in iterations to the level described above.

#### 3.4.3 Combining Full Sky and Cut Sky Maps

Both full sky and cut sky versions of maps were produced for all three years of observations. The final sky maps for each year were produced by using data from the full sky maps to fill in the regions of missing data in the masked maps. Before combining the data from these maps the mean temperature of the full sky Stokes I and the S map required adjustment. For an ideal differential radiometer the mean value of the I and S maps would be undetermined, corresponding to singular modes in . However, inclusion of the non-zero valued transmission imbalance factors, , converts these previously singular modes into modes with very small eigenvalues, indicating that they are very poorly constrained by the measurement. While no physical significance is attached to the recovered values, allowing the mean to vary is required to achieve the level of convergence previously described with regard to the conjugate gradient solution. The fact that the map means are poorly constrained indicates that estimates from the two different map processings may differ significantly due to statistical fluctuations.

Failure to adjust the relative values of the means would introduce an obvious discontinuity at the border of the masked regions when the maps were combined. The procedure used to combine the maps is as follows. First the set of pixels which have the same number of observations in both versions of the sky maps were identified. A constant was then added to the full sky map such that the mean values of the previously identified sets of pixels in the full sky map equaled the mean value of the same set of pixels in each corresponding cut sky map. The pixels containing no observations in the cut sky map were then set to the values of the corresponding pixels in the adjusted full sky maps for the I and S maps. The mean values of the Q and U maps are well determined and were not adjusted. This procedure preserves the mean value for all four map components (I, Q, U and S) from the cut sky maps. Since the matrices were calculated using the cut sky coverage, cut sky maps with the appropriate mean value for use with these matrices may be obtained by simply masking the Galactic plane region.

### 3.5 Evaluation of the Inverse Pixel-Pixel Noise Matrix,

Several steps in the WMAP data processing and spectral processing require an explicit numerical representation of the inverse pixel-pixel noise correlation matrix, . Formally this matrix is described as

(28) |

where and are the inverse noise matrix and the mapping matrix as described in
§ 3.4.2 and § 3.4.1, and and are pixel indices spanning
the I, Q, U and S maps. The sums over and extend over all non-zero values of the
inverse noise matrix evaluated at time . These matrices were evaluated at r4 on a year-by year basis
for all 10 DAs, resulting in 30 matrices. One approximation was used in the evaluation of these matrices to speed processing.
Recall that each observation populates two rows of the mapping matrix, each row corresponding
the projection factors for each of the two radiometers comprising the DA.
Within each row the 8 non-zero elements correspond to
the two differential pixels (A-side and B-side) of the 4 maps I, Q, U and S.
The approximation consists of grouping time contiguous observations for which both the A-side beam and B-side beam remain in the same
r4 pixels. Each such group comprises rows of the mapping matrix (about 15 observations) whose non-zero elements reside
in the same 8 columns. For each contiguous group of observations, the averages of each of the columns containing non-zero values
is calculated for the rows associated with each radiometer. These averages are
tagged with the group’s starting and ending times, the pixel indices of the A-side and B-side beams and the radiometer identification.
Given the starting and stopping
times it is possible to calculate the appropriately weighted sum of corresponding to any pair of groups
and corresponding radiometer.
The sums over and are then calculated for all pairs of groups for which the weighted sum of is non-zero.
This is accomplished by multiplying this weighted sum by the average values of the appropriate columns of the
corresponding groups, and accumulating the values into pixel-pixel noise matrix elements based on the pixel indices
associated with each group. The result is a matrix, , which may be inverted
by standard numerical techniques to form the noise matrix . Using the group average values for the data in the columns
of rather than the exact row by row values in calculation of the inverse noise matrix
is a good approximation for the vast majority of the sky,
since the values being averaged are smooth and slowly varying. However these assumptions fail for pixels very close to the Galactic poles
since the polarizations angles and vary significantly *within* a r4 pixel. The effect of these
approximations has been investigated and is found to have negligible affect on the calculated power spectra. The r4 noise matrices
contained in the data release are calculated from mapping matrices corresponding to cut sky
coverage and are intended for use in regions outside of the
region excluded by the processing mask. Values corresponding to r4 pixels entirely contained within
the processing mask are set to zero.

#### 3.5.1 Projecting Transmission Imbalance Modes from the Matrices

One potential source of systematic artifacts in the sky maps arises from the error in the determination of the input transmission imbalance parameters, . These parameters are measured from the flight data and are used to calculate the estimated dipole signal that is removed from the TOD in preparation for the conjugate gradient map processing. Due to the large amplitude of the dipole signal, even small errors in the measured values of transmission imbalance parameters could introduce significant artifacts into the sky maps. This effect was studied through the use of simulations in which a simulated TOD archive was produced with a given set of values and was analyzed with the input value of as well as values higher and lower to simulate errors in the measurement of . Differences between the maps produced using the correct and incorrect displayed well defined spatial structure confined to low with small variations between DAs arising from the slightly different scan patterns. The geometry of these structures is completely determined by the scan patterns, while their amplitude in the final sky maps is determined by the difference between the value used to process the TOD and the “true” value of . The patterns from these simulation are therefore treated as templates, used to exclude the modes corresponding to the aforementioned spatial structures from subsequent analysis. This is accomplished through the use of modified versions of the inverse noise matrices, , that have these modes projected out. Each modified matrix is calculated as

(29) |

so that

(30) |

In the above expressions designates an outer-product and is the map mode to be removed, obtained from differences between the simulated maps made with the correct transmission imbalance factors and the incorrect factors. Using these forms of the matrices in the spectral analysis (Page et al., 2006) eliminates artifacts which might result from the errant signals cause by errors in the measured values of the transmission imbalance parameters. Figure 12 shows the effect of suppressing these modes on the low- polarization power spectra. This correction has a very small effect on the cosmological parameters, but has been included in the results presented in Page et al. (2006) for completeness.

Inverse noise *covariance* matrices based on these *correlation*
matrices are contained in the data release. The covariance matrices are the correlation matrices
scaled by a multiplicative factor, , where is the noise per observation corresponding to each DA.

### 3.6 Production of Low Resolution Maps

Evaluation of the low- power spectra with optimal signal-to-noise involves use of sky maps and the inverse pixel-pixel noise matrices, , which describe the maps’ noise properties. Since the matrices are produced at r4 due to computational constraints, corresponding low resolution sky maps must be produced that preserve the signal and noise properties as well as possible. There are several difficulties involved with production of low resolution maps. Simply binning maps to lower resolution, using either uniform or pixel-by-pixel (diagonal) inverse noise weighting will result in aliasing of high- signal and noise into lower- modes. This effect can be reduced by applying a low-pass filter to the map before degradation, but such a filter introduces additional noise correlations, not described by the inverse noise matrix as described in the previous section. Even without application of a filter, the omission of the inter-pixel noise correlations results in the degraded maps having noise properties not precisely described by the noise matrix as calculated in § 3.5.

It is possible to produce low resolution sky maps directly using the procedures described in
§ 3.4.2 by accumulating data directly into low resolution maps.
This corresponds to reducing the number of columns in the mapping matrix, , to reflect the reduced number of pixels
in the maps. Such maps will have noise properties accurately described by the pixel-pixel noise matrix, but will have
a *higher* degree of signal aliasing than maps initially produced at high resolution and subsequently degraded. This occurs since
in evaluation of the sky map (see eq.[9]) the low resolution form of the mapping matrix is applied to
the data multiple times, each time introducing aliased signal components, whereas degrading a map initially produced at high
resolution only introduces the aliased components once.

For the Stokes I maps, which are highly signal dominated, degraded res 9 maps are clearly preferable, since the small errors in the noise description are of no consequence and these maps contain lower levels of aliased signal. For the Stokes Q and U maps the signal to noise ratio in the low- regime is on order unity. Low resolution maps of Stokes Q and U were produced both directly at r4 and by degrading r9 maps without application of an anti-alias filter. Analysis of both sets of maps produced very similar low- polarization results for both simulated and flight data. Based on these results it was determined that the noise properties of the degraded form of low resolution polarization maps was adequately described by the pixel noise matrices, so the degraded form is used for all three Stokes parameters.

Although the inter-pixel noise correlations are ignored in the map degradation, the noise correlations between the
Stokes Q and U signals *within* each high resolution pixel are included. Two different versions of reduced
resolution maps are produced and are described in the following sections.

Production of Full Sky Reduced Resolution Maps
The full sky reduced resolution maps (r4) are comprised of two components. The first component, which excludes the Galactic plane region, is
obtained by degrading the r9 cut sky (spm) maps and has noise properties described by the inverse noise matrices
(§3.5). The second component, used in the Galactic plane region, results from the degradation of the
combined maps (§3.4.3) and has noise that is *not* described by the
inverse noise matrices. The noise in this portion of the maps is approximately described by the number of observations values, ,
contained in the sky map data products.

The first component of the reduced resolution map, comprising the high Galactic latitude regions, is generated as follows: During production of the high resolution spm maps, the total weights applied to each pixel in the I, Q, U and S maps, , , , and , are calculated. These are simply the sums of the squares of the corresponding columns of the mapping matrix, . The QU, QS and SU off-diagonal terms, , and , with each r9 pixel are also calculated. These are the sums of the products of the corresponding columns of the mapping matrix. Figure 13 shows maps of the , , , and patterns.

For each high resolution pixel, , an matrix is formed,

(31) |

which differs only by a multiplicative constant from the inverse noise matrix. Correlations between the Stokes I and the polarization related components are not used in the map degradation. The Q, U and S values in a low resolution map pixel, , are then simply the inverse noise weighted average of all the Q, U and S values of the high resolution pixels it contains,

(32) | |||||

(33) |

Here is the r9 processing mask (Hinshaw et al., 2006) containing values of and , and the summations are for all r9 pixels contained within each r4 pixel. Elements of the r4 maps containing one or more observations are retained. The Stokes I maps were degraded using simple weights.

Elements of the r4 maps described above with no observations are filled with values from degraded forms of the combined maps (§3.4.3). The combined maps are degraded following the same procedure, but the matrices are formed using weight values from the cut sky (spm) maps when outside the processing mask, and values from the full sky (fs) maps within the processing mask. Since the weights corresponding to the full sky maps differ from those of the spm maps, the inverse noise matrices, calculated from spm map weights, do not describe the noise in this portion of the map. All the rows and columns of the matrices corresponding to these pixels contain zeros.

Production of the Partial Sky, Foreground Cleaned Reduced Resolution Polarization Maps The reduced resolution sky maps used in the pixel based low- likelihood evaluation (Page et al., 2006) are formed from the foreground cleaned Q and U r9 spm maps, as described in Hinshaw et al. (2006) and Page et al. (2006). Including the spurious term in the map degradation and likelihood evaluation was found to have a negligible effect on the likelihood results (Page et al., 2006). Therefore, to simplify and speed calculation, only the Q and U terms were used in the the degradation and low- likelihood evaluation of the foreground cleaned maps (Page et al., 2006).

The foreground cleaned Q and U maps were degraded in a similar fashion to the full sky maps described above: For each high resolution pixel, , an matrix is formed,

(34) |

The values of the low resolution map pixels, , are the inverse noise weighted sum of the values of all the constituent high resolution pixels,

(35) | |||||

(36) |

In this case is the r9 mask (Page et al., 2006) containing values of and and is a low resolution form of this mask which contains zeros in elements for which more than half of the r9 pixels were masked and values of unity otherwise. These degraded maps are contained in the data release.

## 4 Tests on the Data Set

### 4.1 Limits on Spin-Synchronous Artifacts in the TOD

Among the most troublesome type of systematic errors in experiments such as WMAP are spurious signals which are synchronous with the scanning motion of the instrument. Such artifacts do not integrate to lower levels with the inclusion of additional data as does random noise. The most likely source of such driving signals is from varying insolation due to the observatory’s motion driving temperature variations in the instrument hardware, voltage fluctuation on the spacecraft power bus, or coupling directly into the microwave optics. For all of these effects the errant signal is expected to be synchronous with the 129 s motion of the sun relative to the spacecraft. Following the procedure used in the year-1 analysis (Jarosik et al., 2003b) the radiometer output signals were binned in a spacecraft fixed coordinate system based on the azimuth of the Sun about the spin axis of the spacecraft. When binned in such a manner any signal synchronous with the azimuthal position of the sun about WMAP should be preserved while asynchronous signal components should average away. Combinations of radiometer outputs corresponding to both temperature anisotropies, , and polarization signals, , were accumulated for the entire 3 year mission. The resultant signals were all found to be less that . We therefore believe that there is little contamination of the TOD arising from spin synchronous solar effects and no corrections are applied to the TOD for spin synchronous signals.

### 4.2 Year-to-Year Null Tests

One of the most fundamental tests of the WMAP maps is the year-to-year consistency between maps. Although such tests do not eliminate the possibility of processing or systematic errors that repeat year to year, they do test for spurious signals such as striping from radiometer noise and glitches. A comparison between the year-1 and 3-year processing of the first year data from the V2 DA was presented in § 2.4.1 and displayed a predominantly quadrupolar feature associated with the use of the improved radiometer gain model. Since both maps used in that comparison were generated from the same TOD, the white noise component largely canceled in the difference maps allowing this feature to be seen. In the case of year-to-year comparisons such white noise cancellations do not occur, so subtle effects may not be evident. Nevertheless two sets of year-by-year difference maps are presented to demonstrate the year-to-year consistency, one based on K bands with the largest foreground signal, and one based on V-band with the smallest foreground signal. The top four panels in Figure 14 display the 3 year average and year by year differences of the K-band full sky maps. The 3-year combined map is dominated by Galactic emission and also displays a number of point sources at high Galactic latitudes. The high level of Galactic emission makes K-band one of the more difficult bands to calibrate. While the year-by-year null maps are noise dominated some residual Galactic plane signal is evident. This signal arises from small calibration errors of and is consistent with the stated calibration uncertainty of . Also evident in these maps are some features of both signs at high latitude arising from known variable sources.

The lower four panels in Figure 14 display the 3-year average and year-by-year difference of the V-band data. In this case Galactic emission can be seen at low latitudes and CMB anisotropies at high latitude. The difference maps are again noise dominated with the spatially varying noise properties evident. The increased noise near the ecliptic plane is a consequence of the lower number of observations in these regions as determined by the scan strategy.

Simple power spectra of individual year-to-year null maps are dominated by the radiometer noise bias which grows as when
plotted with the customary scaling, and are not useful for assessing year-to-year map consistency. Figure 15 compares the
best fit WMAP CMB power spectrum with the year-1 to year-2 *c*ross power spectrum (Hinshaw et al., 2003b) of combined V and W band Stokes I sky maps,
averaged in bins of 100. For each year, the linear combination of sky maps (V1-V2) + (W1-W2) + (W3-W4) was formed.
This combination removes the true sky signals to the extent that the bandpasses and gains of each pair of
maps are matched, while the use of the year-to-year cross spectrum eliminates the noise bias since the
radiometer noise is uncorrelated year-to-year. The power spectrum of the null maps are very small relative to the measured power spectrum.

A set of null maps for the Stokes Q and U parameters is displayed in Figure 16 for K-band. Here in the year-combined maps the Galactic foreground emission, predominantly synchrotron, is clearly evident. The year-by-year difference maps are again noise dominated. The spatially varying noise structure for the Stokes Q and U maps is more complicated than that of the Stokes I maps due to the additional polarization projection factors (e.g. ) and the noise covariance between the Stoke Q and Stokes U maps of each DA. The faint cross-like features near the ecliptic poles are regions with large weights, and therefore low statistical noise, as can be seen by comparing these null maps with the maps in Figure 13. Simple null tests, as displayed in Figure 15 for Stokes I, are not useful in assessing year-to-year consistency of the polarization maps due to the much smaller level of the expected signals and more complicated noise properties of the polarization maps. Detailed analysis of the polarization results, including null tests, is presented in Page et al. (2006).

### 4.3 Determination of Map Noise Levels

Evaluation of the instrument noise level is an important aspect of
many analyses involving the WMAP sky maps. In principle, a single noise per observation parameter, ,
should be sufficient to characterize the noise level
for each set of maps (a set consisting of the I, Q and U maps for a given year of data and DA). Such a characterization
would require the use of full resolution pixel-pixel noise correlation matrices,
describing the noise correlations between all combinations of pixel indices and Stokes components.
For r9 sky maps, the size of such matrices precludes their use. However, it is possible to characterize the map noise levels
using only the diagonal components of the noise correlation matrices, , , and , provided the
*intra-pixel* components of noise correlation matrix, , , and (§ 3.6) are included.
This approximation uses different values of the noise per observation parameter for maps of the different Stokes components, parameterizing the
effects of the ignored, off diagonal components of the noise correlation matrices. The different values of these
parameters for the different map components result from the differences between the frequency distributions
of the intensity and polarization information in the TOD, and the radiometers’
frequency dependent, i.e. 1/f, noise properties. The noise levels of the radiometers are also found to be stable year-to-year, so
a single noise per observation value for each sky map/DA combination describes all three individual year sky maps.

#### 4.3.1 Noise Evaluation of the r9 Maps

The noise levels of the r9 maps in the current data release have been estimated using two different methods. In the first method the noise-per-observation, is calculated directly for each DA from a weighted average of the variance of year to year differences of sky maps. For the I maps the noise per observations based on year and year sky maps, , is given as

(37) | |||||

(38) |

where and are the Stokes I value and the number of observations for pixel of the year sky maps, and is the processing mask, used to exclude regions of high Galactic emission from the analysis. Extending this technique to the polarization maps requires a generalization to allow for the noise covariance between the Stokes Q and U maps and the spurious map, S. This is implemented using the relations

(39) | |||||

(40) | |||||

(41) |

where is a vector comprised of the differences between the year and year values of the Q, U and S maps
at pixel , and is the corresponding element matrix describing the
intra-pixel noise correlations between the Q, U and S maps for pixel and year .
Relative calibration errors
between the two maps and the year to year differences in the distributions of
observation pointings *within* pixels will allow some sky signal to leak
into the difference maps, biasing the resultant noise estimates high. Values of and obtained by
averaging the values of the three possible year-by-year difference combinations are presented in Table 5.

The second technique also measures the noise from the year to year difference of the sky maps. Difference maps of each Stokes parameter, I, Q and U, are formed from year and sky maps for each DA. The effective number of observations for each pixel of these difference maps are calculated. For the Stokes I the effective number of observation is given as

(42) |

and for the Stokes Q and U components it is given as

(43) |

and

(44) |

where and are the diagonal components of as defined in equation (41). For each Stokes component the following procedure is applied: 1) Pixels are divided into groups, , based on their values of , 2) The variance of the difference maps for the pixels in each group, , are calculated as well as the mean value of , , 3) These values are fit with two parameters, and in the form

(45) |

where is a parameter that allows for the presence of a statistically isotropic signal that does not scale with . Pixels in high Galactic emission regions, as determined by the processing mask, are excluded from the analysis. Fits of all three year-by-year difference maps are performed and the resultant values averaged together, these results are also presented in Table 5 and agree within with the values obtained using the first technique. Since there is no significant disagreement between the methods we use the results of the first technique for analysis. By comparing the noise level obtained from the year 3 - year 2 difference maps with those from the year 2 - year 1 difference maps we find that the sensitivity of the radiometers has decreased at most by over the 1 year period between the mean epoch of the data sets.

#### 4.3.2 Noise Evaluation the r4 Full Sky Maps

The noise of the full sky reduced resolution (r4) maps may be determined using the pixel-pixel noise matrices. This analysis is performed only on the high Galactic latitude regions, those pixels originating from degradation of the high resolution spm maps. For each DA, noise levels are evaluated for all three possible year-to-year difference map combinations. For combinations of years and the noise per observation, , is calculated as

(46) | |||||

(47) |

where and are reduced resolution sky map sets comprising I, Q, U and S components for years x and y respectively, and and are the pixel-pixel noise correlation matrices for the corresponding years. For each DA, the values from the three year-to-year combinations are averaged together. The ratio of these average values to the noise values of the r9 I maps are presented in Table 5. The differences between the two noise estimates arise from ignoring the pixel-pixel noise correlations in both the r9 - r4 map degradation, and the evaluation of the of the r9 maps.

## 5 Summary

WMAP continued to operate normally throughout its second and third years of science observations. The additional observations served not only to reduce the statistical noise, but also provided for improved understanding of the beam geometries and radiometer performance, important factors in interpretation and processing of the science data. Updated versions of the radiometer gain model, beam models and window functions were presented. New processing procedures that produce maximum likelihood sky maps for Stokes I, Q and U parameters have been described in detail as well as the calculation of the inverse noise matrix. The techniques used to quantify systematic errors were described along with details of the sky map processing used to correct for potential systematic artifacts. Sky maps produced from data from different epochs were used to display the year to year reproducibility of the signal and to measure the instrument noise levels. The performance of the WMAP instrument remains nominal with no significant decrease in sensitivity observed to date.

## Appendix A Sidelobe Correction and Rescaling

Following a notation similar to that of Hinshaw et al. (2003a) the *uncalibrated* time ordered
differential data for each detector, , may be represented as

(A1) |

Here is instrument responsivity (), is the mapping function which transforms spatial information into temporal information and is the sky signal being observed (). All the aforementioned quantities are time dependent as indicated by the use of boldface symbols.

The mapping function can be expressed in terms of the differential beam response of the instrument in instrument coordinates, , and , the time dependent rotation matrix that transforms sky coordinates, , into the satellite coordinates, .

(A2) |

The differential beam response is a signed quantity and is normalized such that

(A3) |

The fraction of the differential beam response in the sidelobes, , is obtained by integrating the magnitude of the differential beam response over the region of the map further than from the boresights of both main beams,

(A4) |

leaving the fraction of the gain in the main beams. As described in Hinshaw et al. (2003a), the calibration of WMAP is tied to the amplitude of the annual dipole signal arising from the motion of WMAP with respect to the solar system barycenter. If the annual dipole sky signal is represented as , the corresponding component of the uncalibrated radiometer signal is

(A5) |

Calibration is accomplished by performing a least-squares fit of the uncalibrated data to a predicted time ordered data set, , based on the known monopole temperature,

(A6) |

yielding an estimate of the radiometer gain, . The predicted time ordered data set is obtained using a mapping function assuming ideal ‘pencil’ beams, and the annual dipole sky signal,

(A7) |

The actual calibration procedure is iterative and provides for subtraction of better estimates of the fixed sky signal before each new gain determination. Once the iterative steps are completed the radiometer gain model (Sec. 2.4) is fit to smooth and interpolate the hourly values. Taken together, the above calibration procedure effectively measures by evaluating the mean of the ratio

(A8) |

Note that the mapping function in the numerator incorporates the full sky beam response. This function can be written as the sum of a main beam mapping function, , and a sidelobe mapping function, ,