text:
# First Year Wilkinson Microwave Anisotropy Probe (WMAPWMAP is the result of a partnership between Princeton University and NASA’s Goddard Space Flight Center. Scientific guidance is provided by the WMAP Science Team. ) Observations: Data Processing Methods and Systematic Errors Limits

## 1 INTRODUCTION

The Wilkinson Microwave Anisotropy Probe (WMAP) has produced full-sky maps of the cosmic microwave background (CMB) of unprecedented precision and accuracy. On angular scales larger than $`0\text{.}^{}5`$, the dominant uncertainty is not the instrument noise but rather the “cosmic variance” inherent when analyzing a single realization (the observable universe) of a random process. In the cosmic variance limit, no further improvement can be made by reducing instrument noise, placing even greater importance on the minimization of non-random instrumental effects in the data.

The WMAP design emphasizes control of systematic errors (Bennett et al., 2003a). The observatory was designed with a detailed systematic error budget in place, and a mature data analysis pipeline was written early to help inform many of the design decisions. Differential radiometers compare the temperature from widely-separated regions of the sky through back-to-back matched optics. Common-mode signals thus cancel before affecting the sky maps. The radiometer feed horns only illuminate a fraction of the primary mirrors, reducing the sidelobe response in the beam patterns. The instrument was designed to have minimal response to electrical or thermal perturbations and operates in an exceptionally stable environment at the second Earth-Sun Lagrange point. The observatory’s compound spin and precession allow rapid inter-comparison of different positions on the sky, greatly reducing the coupling of systematic error signals into the sky maps and effectively symmetrizing the beam response. WMAP data are calibrated in flight using the dipole modulation of the CMB from the observatory’s orbital motion around the Sun as a full-beam calibration source. We measure the beam pattern in flight using observations of the planet Jupiter.

We characterize or limit systematic errors in the WMAP first-year data using flight data supplemented where necessary with results from ground tests. Systematic errors may be classified into several broad categories including the following:

Calibration Errors.— The time-ordered data is simultaneously fit for the calibration and sky map. An iterative algorithm updates the calibration solution based on the previous iteration of the sky map solution, and vice-versa. The most important source of error is confusion between the dipole signal and higher-order sky signal, especially bright Galactic foreground emission in the low frequency WMAP bands. See §2.2.1.

Map-making errors.— These are due to poor convergence in the sky map solution, or to errors in the determination of the spacecraft pointing. See §2.1.1 and §3.3.

Beam Errors.— Instrument noise, background subtraction, and pointing errors can limit the in-flight measurement of the beam response from Jupiter. Although the beams are not symmetric, the observatory’s compound spin and precession effectively symmetrize the beam response. Uncertainties in both the beam solid angle and the window functions must be characterized. See Page et al. (2003a) for a complete discussion of beam mapping. We summarize and incorporate their results in §3.3.3.

Sidelobe Response.— Sidelobe pickup of bright sources (e.g. the Galactic plane) introduces an additive signal dependent on the orientation of the beams on the sky. Barnes et al. (2003) discuss the sidelobe response of each radiometer and estimate the effect on the first-year sky maps.

Baseline Errors.— Thermal or electrical perturbations can produce signals dominated by an additive term in the time domain. Slow drifts are removed as part of the calibration procedure, but signals near the spin period can couple to the sky maps with some efficiency. See §3.4.

Striping.— Correlations in the time-ordered data from sources not fixed on the sky (e.g. $`1/f`$ noise or post-detection filtering) introduce correlated noise in the sky maps. Application of a pre-whitening filter to the time-ordered data reduces this effect. See §2.3.2.

We have constructed a detailed model of the instrument that successfully reproduces all major aspects of the instrument performance. Software simulations using this model validate the map-making algorithm and allow us to assess the effect in the sky maps of various signals in the time-ordered data. Based on this, we apply a pre-whitening filter to the time-ordered data to remove a low level of $`1/f`$ noise. We infer and correct for a small ($``$1%) transmission imbalance between the two sky inputs to each differential radiometer, and we subtract a small sidelobe correction from the 23 GHz (K band) map prior to further analysis. No other systematic error corrections are applied to the data. Calibration and baseline artifacts, including the response to environmental perturbations, are negligible. Systematic uncertainties are comparable to statistical uncertainties in the characterization of the beam response. Both are accounted for in the covariance matrix of the window function and are propagated to uncertainties in the final power spectrum. We characterize the combined upper limits to residual systematic uncertainties through the pixel-pixel covariance

This paper is organized as follows. In §1.1 we define the terms and notation used throughout the paper. In §2 we discuss the iterative algorithm for making maps from time-ordered data, then generalize to the case of simultaneous calibration and sky map estimation. (Appendix C further generalizes to map-making with polarization data.) We also discuss the noise properties of the time-ordered data. In §3 we discuss combined systematic error limits from calibration and map-making. We also present the noise properties of the sky maps in terms of their pixel-pixel covariance. Finally, we derive upper limits to environmental perturbations and summarize the combined systematic error budget. In §4 we present our conclusions.

### 1.1 Notation and Overview

Throughout this paper, we denote vectors and scalars with bold and plain lowercase letters, respectively. Matrices and operators are denoted with uppercase bold letters. Following Stompor et al. (2002) we denote vector and matrix component indices in parentheses, saving subscripts and superscripts to further identify quantities. A summary of the most frequently used symbols is given in Table 1. Unless otherwise stated, all temperatures are specified in thermodynamic units.

WMAP measures the brightness temperature of the sky as a function of position, $`𝐭(\theta ,\varphi )𝐭(p)`$, where $`p`$ denotes the sky map pixel number, indexed from 0, in HEALPix nested format (Górski et al., 1998). To make this measurement, WMAP scans the sky and measures the temperature difference between two points at time $`t`$. The resulting time-ordered differential data (TOD) is denoted $`𝐝`$. The main goal of the map-making is to obtain the minimum variance estimate of the sky map, $`\stackrel{~}{𝐭}`$, by inverting the raw differential data. Note that $`\stackrel{~}{𝐭}`$ is the true sky temperature convolved with the nominal instrument beam, plus instrument noise. In the process of solving for the map, we calibrate the data by estimating the gain and baseline from the flight data itself; characterize the full instrument beam response function from observations of the planets; characterize the noise spectrum of the instrument, and place limits on residual systematic errors.

In order to produce stable data with a nearly-white noise spectrum WMAP employs 20 high-electron-mobility-transistor (HEMT) based differential radiometers. Each radiometer measures the brightness difference between two inputs, one fed by an A-side beam, the other by a B-side beam approximately 141apart. A detailed description of their design and fabrication may be found in Jarosik et al. (2003a); a summary of their in-flight performance is presented in Jarosik et al. (2003b). The 20 radiometers form 10 polarization-sensitive “differencing assemblies” (DA) which are designated based on their frequency or waveguide band: K1, Ka1, Q1, Q2, V1, V2, W1, W2, W3, W4. The two radiometers in a DA are sensitive to orthogonal linear polarization modes; the radiometers are designated 1 or 2 (e.g., K11 or K12) depending on which polarization mode is being sensed. Each of the 20 radiometers is intrinsically a 2-channel device, with channels designated 3 and 4 in the flight telemetry, e.g., K113 or K114. \[Channels 3 and 4 were designated left and right, respectively, in Jarosik et al. (2003a).\] There are 40 such data channels in the flight telemetry. As discussed below, each of the 40 channels is individually calibrated, then the 4 channels from a single differencing assembly are combined to form differential intensity and polarization signals as follows. Let $`𝐝_{ij}`$ be the calibrated differential signal from a single channel, $`j`$, of radiometer $`i`$. The differential intensity data (Stokes parameter $`I`$) is the average of all 4 channels

$$𝐝=\frac{1}{2}(𝐝_{13}+𝐝_{14})+\frac{1}{2}(𝐝_{23}+𝐝_{24}).$$
(1)
The differential polarization data is obtained by taking the difference between the two radiometer signals

$$𝐩=\frac{1}{2}(𝐝_{13}+𝐝_{14})\frac{1}{2}(𝐝_{23}+𝐝_{24}).$$
(2)
In Appendix C on polarization map making, we relate the differential polarization signal to the Stokes parameters $`Q`$ and $`U`$. Kogut et al. (2003) discuss additional aspects of polarization mapping and analyze the first-year temperature-polarization correlation data based on these maps. Note that we can also form null channels from the data by taking channel differences, $`(𝐝_{i3}𝐝_{i4})`$. As discussed in §3.2, these channel combinations provide valuable consistency tests for the final sky maps.

A single channel of uncalibrated differential data may be modeled as

$$𝐜=𝐠\left[𝐌(𝐭+𝐭_𝐬)+𝐧\right]+𝐛,$$
(3)
where each quantity is a function of time:

* $`𝐜(t)`$: raw differential data, in counts or digital units (“du”).
* $`𝐠(t)`$: instrument responsivity (here called gain), in du mK<sup>-1</sup>.
* $`𝐛(t)`$: instrument baseline, in du.
* $`𝐧(t)`$: instrument noise, in mK.
* $`𝐌(𝐭+𝐭_𝐬)𝚫𝐭(t)`$: differential sky signal from all sources, convolved with the instrument beam, in mK. This includes fixed sources, $`𝐭`$ (e.g., CMB and Galactic emission) and moving sources, $`𝐭_𝐬`$ (e.g., planets).

In practice, the differential signal is integrated for a fixed time $`\tau `$ and sampled at discrete times $`t_i`$, thus we may regard time series data as a vector with $`N_t`$ observations. The integration time per observation is 128.0 ms, 128.0 ms, 102.4 ms, 76.8 ms, and 51.2 ms for bands K through W, respectively.

The differential temperature at time $`t`$ is the convolution of a time-dependent mapping function, $`𝐌`$, with the sky signal at time $`t`$, $`𝐭+𝐭_𝐬(t)`$

$$𝚫𝐭(t)=𝑑\mathrm{\Omega }_𝐧^{}𝐌(𝐧(t),𝐧^{})(𝐭(𝐧^{})+𝐭_𝐬(𝐧^{},t)).$$
(4)
Here $`𝐭=𝐭_𝐜+𝐭_𝐠`$ is the fixed sky signal consisting of CMB anisotropy, $`𝐭_𝐜`$, and Galactic foreground signals, $`𝐭_𝐠`$, while $`𝐭_𝐬`$ represents all time dependent sources, especially the Sun, Earth, and Moon, which are potentially visible in the far sidelobes of the instrument. The operator $`𝐌`$ can be represented as an $`N_t\times N_p`$ matrix where each row is the normalized, full-sky beam response in sky-fixed coordinates as given by the scan pattern. Several features of the mapping function that pertain to the data processing are discussed in Appendix A. The main beam response is mapped in flight using observations of Jupiter as a far-field point source (Page et al., 2003a). An important aspect of the WMAP optical design (Page et al., 2003b) was to limit sidelobe pickup to negligible levels and to have the effective beam response in the final sky maps be approximately circularly symmetric. We discuss each of these topics in more detail in separate papers, (Barnes et al., 2003; Page et al., 2003a), while this paper summarizes the main results in terms of systematic error limits in the sky maps.

The instrument gain, baseline, and noise are determined from the flight data itself. This is an iterative process that we discuss in detail in §2.2.1. Here we provide a brief overview of our terminology in order to frame the following discussion of systematic errors. Let the gain and baseline measured in flight be $`\stackrel{~}{𝐠}`$ and $`\stackrel{~}{𝐛}`$, respectively. The calibrated differential signal is then

$$\stackrel{~}{𝐝}=\frac{(𝐜\stackrel{~}{𝐛})}{\stackrel{~}{𝐠}}=\frac{𝐠}{\stackrel{~}{𝐠}}𝚫𝐭+\frac{𝐠}{\stackrel{~}{𝐠}}𝐧+\frac{𝐛\stackrel{~}{𝐛}}{\stackrel{~}{𝐠}}.$$
(5)
With calibrated data available, the sky map is obtained by evaluating the linear equation

$$\stackrel{~}{𝐭}=𝐖\stackrel{~}{𝐝},$$
(6)
where $`𝐖`$ is a linear operator defined in §2.1. The properties of $`𝐖`$ are determined by the scan strategy of the observatory and the noise properties of the time-ordered data, $`𝐝`$.

The WMAP scan pattern is an integral part of the mission design (Bennett et al., 2003a). It consists of a compound spin and precession centered about the Sun–WMAP line, with parameters as given in Table 2. There are several aspects of this scan strategy that are important for high quality data: scans of a given pixel cross at many angles so that the effective beam response is symmetric; a given pixel is observed on many different time scales from minutes to months; the angular velocity of a given line of sight is nearly constant on the sky; the instrument observes a large fraction ($`>`$30%) of the sky each day; and the time-average of the differential data is approximately zero over an hourly calibration period, allowing for robust initial baseline estimation.

## 2 THE MAP-MAKING PIPELINE

A graphical overview of the WMAP data processing and analysis pipeline is shown in Figure 1 of Bennett et al. (2003b). The heart of the pipeline is a set of programs that bring science and housekeeping data from the Science and Mission Operations Center (SMOC) through to a set of calibrated full sky maps for each of the 10 WMAP differencing assemblies. Numerous additional programs are used to generate ancillary data products such as beam response maps, calibration files and analysis products. In this section we provide a high-level overview of each program in the map-making pipeline and provide references to later sections of this paper, or to companion papers, for further details.

Raw telemetry data from the satellite is transferred approximately once per day through NASA’s Deep Space Network (DSN), to the SMOC, located at the Goddard Space Flight Center. The SMOC monitors the basic health and safety of the Observatory, sends all commands, and requests re-transmissions of data that were missing from a previous transmission. The data are then “level-0” processed into a set of time-ordered, daily files which contain science data, instrument housekeeping data, spacecraft data (including attitude and ephemeris data), and event message files. These files are then transferred to the Science Team’s processing facility, also at the Goddard Space Flight Center.

Every time a new full day of data arrives, a series of automated procedures perform the following tasks: 1) Generate a standard set of daily plots that are archived and visually inspected; 2) Generate a reduced “trending archive”, which consists of subsets of the data sampled once every 10 minutes. In the case of the science data, we record the mean and rms of each channel in a 10 minute interval, whereas the housekeeping data are sub-sampled, once per 10 minutes; 3) Perform a series of data quality checks that search the data for violation of pre-set range limits or excessive time-gradients in the telemetry signals. Limit violations are logged and notification is sent to a member of the science team via e-mail. (Initial limit tests are performed at the SMOC as well.)

At selected time intervals, the level-0 telemetry files are collated by a pre-preprocessor into a master archive of raw (uncalibrated) data. The major functions of the pre-processor are to: 1) collate the science and housekeeping data into single daily files; 2) flag data that is suspected or known to be unusable; 3) interpolate the attitude and ephemeris data to times that are commensurate with the science data time stamps; and 4) apply a coarse flag to data that is within 7 of one of the outer planets (Mars through Neptune) so that it may be rejected from the initial sky maps, but identified for beam mapping (Page et al., 2003a). (Only Jupiter data is used for making the final beam maps.)

Initial sky maps and calibration data are generated from the raw archive using the iterative map-making algorithm first described in Wright et al. (1996), and further described in §2.1. As discussed in §2.2, the initial calibration is determined by fitting the raw time-ordered data to the known signal produced by the CMB dipole. Because the sky signal contains significant higher-order power ($`l>1`$), the calibration solution must be iteratively improved in concert with the initial sky map iterations. The convergence of this simultaneous fit has been demonstrated with end-to-end simulations, which are also described in §2.2.

The calibration solution converges more rapidly than the sky map does, so we freeze the initial calibration solution after $``$10 iterations before proceeding to convergence with the sky maps. With the initial dipole-based calibration data in place, a re-processor generates a refined gain and baseline solution and applies this to the data. The program also updates the data quality flags, as necessary, then writes the calibrated data to a new final time-ordered archive. The refinements to the dipole-based gain solution are discussed in Jarosik et al. (2003b) and in §2.2.1. The initial baseline solution is refined with a pre-whitening filter (Wright, 1996) which is presented in detail in §2.3.2.

The final sky maps are computed using the final calibrated time-ordered data as input. The first-year sky maps required 20 post-calibration iterations to be sufficiently converged. The map-making algorithm is fundamentally the same as is used in the initial estimates, but we add some refinements for this final stage of processing. These include: 1) correcting for a small ($`\genfrac{}{}{0pt}{}{_<}{^{}}`$1%) loss imbalance between the A and B-side sky beams. Jarosik et al. (2003b) demonstrate that this effect is nearly orthogonal to the gain solution, so its inclusion after the calibration processing does not invalidate the gain solution. 2) We weight individual time-ordered observations by their proper statistical weight to account for the small change in instrument noise ($`<`$1%) over the course of a year due to the 0.9 K seasonal temperature variation experienced by the instrument cold stage. 3) We compute the planet–boresight angle for each observation to minimize unnecessarily conservative data loss. The criterion used for the first-year maps is a cut of radius $`1\text{.}^{}5`$.

We also generate polarization maps using a generalization of the temperature algorithm, the details of which are presented in Appendix C. For the first-year data processing, we have generated maps of the Stokes parameters $`Q`$ and $`U`$, but we have not yet fully characterized all of the potential systematic errors in these maps. However, the temperature-polarization correlation data are much less prone to systematic errors than the polarization auto-correlation data. Kogut et al. (2003) have analyzed the temperature-polarization data and, supported by systematic error limits from Barnes et al. (2003), they find a significant correlation, including the signature of a relatively early epoch of cosmic reionization.

In the remainder of this paper, we present a more detailed description of each of the map-making and calibration procedures, including an assessment of their performance with the first-year WMAP data. We then derive detailed systematic error limits applicable to the time-ordered data, to the sky maps, and to the angular power spectrum.

### 2.1 Map-Making with Pre-Calibrated Data

While the process of generating the final sky maps from calibrated data comes last in the map-making pipeline, we discuss the algorithm first because the algebra of map-making is central to the entire data processing scheme, and it helps to guide the systematic error analysis.

We consider the problem of estimating a sky map, $`𝐭`$, from calibrated, differential time-ordered data, $`𝐝`$, which is a linear function of the sky map

$$𝐝=\mathrm{𝐌𝐭}+𝐧,$$
(7)
where $`𝐌`$ is the mapping function of the experiment, which has $`N_p`$ columns and $`N_t`$ rows, $`N_p`$ is the number of sky map pixels and $`N_t`$ is the total number of time-ordered observations. In its simplest form, each row (observation) of the scan matrix contains a $`+1`$ in the column (pixel) seen by the A-side beam, and a $`1`$ in the column (pixel) seen by the B-side beam. This matrix can be generalized to include the effects of beam convolution, but for WMAP these refinements are small and are being deferred to future processing. The effects of a differential loss imbalance between the A and B-side beams is readily accounted for by using values different from $`\pm `$1 in $`𝐌`$. An analysis of this effect in the WMAP radiometers is presented in Jarosik et al. (2003b). The details of how we account for it in the pipeline are given in Appendices A and C.

The noise $`𝐧`$ is assumed to have zero mean and covariance $`𝐍`$,

$`𝐧`$ $`=`$ $`0`$ (8)
$`\mathrm{𝐧𝐧}^T`$ $`=`$ $`𝐍.`$ (9)
We defer a detailed discussion of the WMAP noise properties to §2.3.2, but for most radiometers it is reasonable to approximate the noise covariance as diagonal, $`𝐍\sigma _0^2𝐈`$, where $`\sigma _0`$ is the rms noise per observation and $`𝐈`$ is the identity matrix, though this assumption is not required for the algorithm described below to converge.

The least-squares, or maximum-likelihood estimate of the sky map, $`\stackrel{~}{𝐭}`$, results from solving the normal equations

$$\stackrel{~}{𝐭}=(𝐌^T𝐍^1𝐌)^1(𝐌^T𝐍^1𝐝).$$
(10)
More generally, we obtain an unbiased estimate of the sky map by choosing any symmetric matrix $`𝐒`$ in place of $`𝐍^1`$ (Tegmark, 1997). To see this, substitute $`\mathrm{𝐌𝐭}+𝐧`$ in place of $`𝐝`$ in the above equation to get

$$\stackrel{~}{𝐭}=(𝐌^T\mathrm{𝐒𝐌})^1(𝐌^T𝐒[\mathrm{𝐌𝐭}+𝐧])=𝐭+(𝐌^T\mathrm{𝐒𝐌})^1(𝐌^T\mathrm{𝐒𝐧}).$$
(11)
Thus $`\stackrel{~}{𝐭}`$ reduces to $`𝐭`$ plus a noise term that is independent of $`𝐭`$ and has zero mean over an ensemble average. For the first-year processing we take $`𝐒=𝐈`$. To simplify notation, we define a matrix $`𝐖(𝐌^T𝐌)^1𝐌^T`$ in which case the map solution is $`\stackrel{~}{𝐭}=\mathrm{𝐖𝐝}`$.

The pixel-pixel noise covariance in the sky map solution is

$`𝚺`$ $`=`$ $`(\stackrel{~}{𝐭}𝐭)(\stackrel{~}{𝐭}𝐭)^T=(\mathrm{𝐖𝐧})(\mathrm{𝐖𝐧})^T=\mathrm{𝐖𝐍𝐖}^T`$ (12)
$`=`$ $`(𝐌^T𝐌)^1(𝐌^T\mathrm{𝐍𝐌})(𝐌^T𝐌)^1.`$
In the limit that $`𝐍`$ is diagonal and the rms noise per observation, $`\sigma _0`$, is constant, $`𝚺`$ reduces to $`\sigma _0^2(𝐌^T𝐌)^1`$. For the WMAP scan pattern, the matrix $`𝐌^T𝐌`$ is diagonally dominant with diagonal elements $`𝐧_{\mathrm{obs}}(p)`$, the number of observations of pixel $`p`$ by either the A or B-side beam. Thus, to a very good approximation, the pixel-pixel covariance matrix is diagonal

$$𝚺(p_i,p_j)\frac{\sigma _0^2}{𝐧_{\mathrm{obs}}(p_i)}\delta _{ij}.$$
(13)
Values for $`\sigma _0`$ are given by Bennett et al. (2003b).

The leading order off-diagonal terms occur at the beam separation angle ($`\theta _{\mathrm{beam}}141^{}`$), and are of order 0.3% of the diagonal elements. If the time-ordered noise $`𝐍`$ is not diagonal, then maps produced with the above algorithm will have correlated noise (stripes) along the scan paths defined by $`𝐌`$. This is a small, but not negligible, effect for some of the WMAP radiometers, and is entirely negligible for others. The noise properties of the time-ordered data and sky maps are further discussed in §2.3.2 and §3.2, respectively.

#### 2.1.1 Iterative Map Making

The evaluation of the sky map solution $`\mathrm{𝐖𝐝}`$ requires the inversion of the $`N_p\times N_p`$ matrix $`𝐃𝐌^T𝐌`$. We use the iterative approach introduced by Wright et al. (1996) to evaluate this expression. Briefly, suppose we have an initial guess for the sky map, $`𝐭_0`$, which differs from the true sky map, $`𝐭`$, by $`\delta 𝐭_0=𝐭_0𝐭`$. Then $`𝐃𝐭_0=𝐃(𝐭+\delta 𝐭_0)`$ can be recast as

$$𝐃\delta 𝐭_0=𝐃𝐭_0𝐌^T𝐝,$$
(14)
where we have used the fact that $`𝐃𝐭=𝐌^T𝐝`$. As noted above, $`𝐃`$ is diagonally dominant, so a good approximate inverse for $`𝐃`$ is

$$\stackrel{~}{𝐃}^1(p_i,p_j)\frac{1}{𝐧_{\mathrm{obs}}(p_i)}\delta _{ij}.$$
(15)
This leads to the approximate solution for the residual, $`\delta 𝐭_0\stackrel{~}{𝐃}^1[𝐃𝐭_0𝐌^T𝐝]`$, and suggests the following iterative solution

$`𝐭_{n+1}`$ $`=`$ $`𝐭_n\delta 𝐭_n`$ (16)
$`=`$ $`𝐭_n\stackrel{~}{𝐃}^1[𝐃𝐭_n𝐌^T𝐝]`$ (17)
$`=`$ $`(\stackrel{~}{𝐃}^1𝐌^T)𝐝+(𝐈\stackrel{~}{𝐃}^1𝐃)𝐭_n.`$ (18)
The interpretation of equation (18) is that for each pixel the new sky map temperature is the average of all differential observations of that pixel (accounting for the sign of the observing beam) corrected by an estimate of the signal in the paired beam, based on the previous sky map iteration. The expression in equation (18) can be efficiently evaluated because the sums can be accumulated by reading through the time-ordered data from disk, each iteration, and accumulating data into arrays of length $`N_p`$

$`𝐧_{\mathrm{obs}}𝐭_{n+1}(p_A)𝐧_{\mathrm{obs}}𝐭_{n+1}(p_A)+w_i[𝐝(t_i)+𝐭_n(p_B)]`$  $`𝐧_{\mathrm{obs}}(p_A)𝐧_{\mathrm{obs}}(p_A)+w_i`$
$`𝐧_{\mathrm{obs}}𝐭_{n+1}(p_B)𝐧_{\mathrm{obs}}𝐭_{n+1}(p_B)w_i[𝐝(t_i)𝐭_n(p_A)]`$  $`𝐧_{\mathrm{obs}}(p_B)𝐧_{\mathrm{obs}}(p_B)+w_i,`$ (19)
where $`w_i=1`$ in the initial sky map processing and is proportional to a noise weight (equation C31) in the final sky map processing. Note that it is never necessary to store or invert an $`N_p\times N_p`$ matrix.

We have tested this algorithm extensively with flight-like simulations. In this section we present results for an “ideal” noiseless instrument with circular beams and perfect calibration to isolate the performance of the map-making algorithm from other effects. More realistic data models are introduced to the simulation in subsequent sections. Figure 1 shows a sample residual map, $`𝐭_{out}𝐭_{in}`$, generated from a one-year simulation of Q2 data. The input sky map included realistic CMB signal with a peak-to-peak amplitude of $`\pm `$420 $`\mu `$K, and a Galactic signal with a peak brightness of $``$50 mK. The output sky map is recovered with an rms error of $`<0.2`$ $`\mu `$K, after 50 iterations. The dominant structure in the residual map is a mode aligned with the ecliptic plane. The power in this mode is concentrated in spherical harmonic mode $`l=4`$, due to a combination of the WMAP scan strategy and the beam separation angle. This is the mode on the sky that is least well measured by WMAP (except for the monopole!) and is thus the slowest to converge, though additional iterations would reduce its amplitude even further. The final first-year flight sky maps were effectively run for 80 iterations. Since the rms error associated with this term is very small, and since we build up a more realistic data model in subsequent simulations, we do not further quantify this contribution to the final systematic error budget. Rather, we subsume it into an overall map-making and calibration error budget that includes this and other effects together.

This iterative approach to map-making is readily generalized to polarization maps as well – the formalism is presented in Appendix C. We have tested that algorithm with the same simulations used above to test the temperature algorithm and find that the polarization maps converge even faster than the temperature maps. After 10 iterations, the map-making artifacts in a residual polarization map are $`<`$0.05 $`\mu `$K peak-peak.

### 2.2 Combined Calibration and Map Making

The processing algorithm described above assumes that the data have already been calibrated. In practice, we use the above algorithm in the second stage of map-making, after an initial stage in which we simultaneously solve for the radiometer calibration and the sky map. In the initial stage of map-making, we employ the same iterative algorithm to solve for the map, but rather than processing straight through the time-ordered data on each iteration, we process the data one hour at a time, pausing to solve for the calibration in each radiometer channel, before accumulating the calibrated data. The calibration solution then iteratively improves as the sky model improves. The following sub-sections lay out the procedure in detail, and present results for the flight data with an assessment of its precision and accuracy based on flight-like simulations.

#### 2.2.1 Instrument Calibration from the Dipole Modulation

For a sufficiently short period of time the instrument gain and baseline can be approximated as constant, $`𝐜g_k(𝚫𝐭+𝐧)+b_k`$, where $`g_k`$ and $`b_k`$ are the gain and baseline during the $`k^{\mathrm{th}}`$ 1-hr calibration period. Since the sky signal $`𝚫𝐭`$ is dominated by the CMB dipole measured by COBE, $`𝚫𝐭_𝐝`$, a single channel of raw data can be modeled as

$$𝐜_𝐦(g_k,b_k)=g_k(𝚫𝐭_𝐝+𝚫𝐭_𝐯)+b_k,$$
(20)
where $`𝚫𝐭_𝐯`$ is the additional dipole moment induced by the motion of WMAP relative to the solar system barycenter (the rest frame of the COBE dipole).

We fit for the gain and baseline in each calibration period, $`k`$, by minimizing

$$\chi ^2=\underset{ik}{}\frac{\left[𝐜(t_i)𝐜_𝐦(t_i|g_k,b_k)\right]^2}{\sigma _0^2},$$
(21)
where $`i`$ is a time-ordered datum index. We omit data that are flagged as unusable, and data when either the A or B-side beam points within a Galactic pixel mask. The mask used for this latter application is the Kp8 mask defined in Bennett et al. (2003c), without edge smoothing. This mask is used throughout the map-making pipeline. The fit is performed for each of the 40 WMAP channels independently. To minimize the covariance between the recovered gain and the baseline, it is necessary to have a scan strategy such that the time average of the sky signal, $`𝚫𝐭`$, is approximately zero in one calibration period. The combined spin and precession of WMAP is designed to produce time-ordered data that satisfies this requirement. For example, in K band, which has the largest sky signal, a 1-hr running mean of the differential sky signal has an rms fluctuation of 14 $`\mu `$K, compared to a dipole signal of greater than 3 mK. After each hour of data is processed for the calibration solution, the data are accumulated as per equation (19) to develop the sky map solution.

The largest source of error in the calibration fit is due to un-modeled sky signal from the CMB anisotropy and Galactic foreground emission, $`𝚫𝐭_𝐚𝚫𝐭𝚫𝐭_𝐝`$. This projects onto the dipole signal and, as shown below, causes errors in the gain solution as large as 5-10% in K band, where the Galactic signal is strongest. The calibration fit may be iteratively improved by subtracting an estimate of the anisotropy from the raw data prior to fitting. In particular, let $`g_k^{}`$ be the gain inferred for calibration period $`k`$ from the previous iteration of the calibration fit. Then minimize

$$\chi ^2=\underset{ik}{}\frac{\left[𝐜^{}(t_i)𝐜_𝐦(t_i|g_k,b_k)\right]^2}{\sigma _0^2},$$
(22)
where

$$𝐜^{}=𝐜g_k^{}𝚫𝐭_𝐚^{}$$
(23)
and $`𝚫𝐭_𝐚^{}`$ is the differential sky signal (less the dipole component) computed from the previous sky map iteration. This process is repeated until the calibration solution is sufficiently converged.

Note that the absolute calibration is tied to the time-dependent portion of the dipole signal, $`𝚫𝐭_𝐯`$; we use the fixed dipole as a short-term transfer standard only. In particular, when we update the sky model and apply the anisotropy correction in equation (23), any error in the fixed dipole moment, $`𝚫𝐭_𝐝`$, is assumed to be anisotropy, and is applied as a correction in the same way. For a data set of at least one year in length (one full cycle of $`𝚫𝐭_𝐯`$), the error in the absolute calibration will be essentially orthogonal to any error in the fixed dipole $`𝚫𝐭_𝐝`$.

#### 2.2.2 Performance of the Dipole-Based Gain Solution

As an illustration of the systematic gain error induced by higher-order ($`l>1`$) anisotropy, Figure 2 shows an example of the gain solution convergence from a one-year low-noise simulation. This simulation implements the simultaneous calibration and sky map estimation discussed above and was run for 30 iterations. The example shown is for one channel of K band data (the worst case) which exhibits $``$7% errors after one iteration, corresponding to a sky model that has only a dipole component. After 30 iterations, the residual errors are $`<`$0.1% over the entire year. Similar, or better, performance is achieved for all other WMAP channels.

In processing the final first-year maps, the combined calibration and map-making code was run for 10 iterations. However the initial sky model was based on an earlier “pathfinder” run of the pipeline that ran for a total of 30 iterations of combined calibration and map-making plus an additional 20 iterations of sky map convergence. Thus we conservatively estimate that the combined absolute and relative calibration errors due to incomplete calibration convergence to be $`<`$0.1%. We defer a discussion of the final calibration uncertainty to §3.1.

Figure 3 shows a sample of the converged gain solution from equation (22) for two WMAP channels, K113 and V113. Note that the radiometer gains are typically drifting by a few percent over the course of the first year. As we show below, the dipole-based fits easily track drifts at this level. The noise in the gain solution is typically a few percent per hourly calibration period, though, as is readily seen in the figure, the noise level changes with time of year as the scan pattern sweeps around the CMB dipole. The V113 gain exhibits an additional modulation that is clearly correlated with the physical temperature of the instrument. However, the time scale of the temperature change is slow enough that the corresponding gain changes are well tracked by the dipole fits. Quantitative limits on thermally induced gain and baseline errors are discussed in Jarosik et al. (2003b) and in §3.4.1. A summary of the gain statistics from the flight data is given in Table 3.

#### 2.2.3 The Initial Baseline Solution

The bottom two panels of Figure 3 show the converged baseline solution resulting from the fit in equation (22) for one year of K113 and V113 flight data. The fits have had a mean subtracted, and have been divided by the gain to convert to temperature. These plots, which are representative of all 40 channels, show that the offsets of the radiometers are typically stable to $`\pm `$5 mK over the course of the first year. Simulations demonstrate that the hourly baseline solution is unbiased. However, it is also clearly noisier than optimal, consistent with the flight measurements of the noise power spectral density (Jarosik et al., 2003b). In §2.3.2 we describe an improved baseline model that is based on the application of a pre-whitening filter tailored to the measured noise spectrum of each channel.

### 2.3 Improving the Calibration Model

The sky maps obtained with the hourly calibration are reasonable; however the noise in the calibration solution, particularly in the baseline, is significantly higher than optimal, and the use a piecewise continuous calibration in the final maps would introduce striping in its own right. In the following subsections, we present the steps undertaken to filter the gain and baseline solutions that enter into the final sky maps. Prior to generating the final maps, this refined calibration is applied to the data and written to disk as the final first year calibrated time-ordered archive.

#### 2.3.1 The Gain Model

Jarosik et al. (2003b) present a physical model for the gain that is based on the RF bias, or “total power” measured in each channel, and on the physical temperature of the instrument cold stage, which is monitored with high resolution platinum resistance thermistors (PRTs). Each of these quantities is recorded once every 23 s in the engineering telemetry with a relative noise that is substantially lower than the noise in the dipole-based gain solution. Thus, if the model fits the dipole-based data satisfactorily, it offers a means for measuring the gain with more precision, and on time scales shorter than the spin period. The model for the gain, $`𝐠(t)`$, has the form

$$𝐠(t)=g_0\frac{\overline{V}(t)V_0}{T_{\mathrm{FPA}}(t)T_0},$$
(24)
where $`\overline{V}`$ is the measured RF detector bias, $`T_{\mathrm{FPA}}`$ is the measured temperature of the FPA, and $`g_0`$, $`V_0`$, and $`T_0`$ are fit constants. See Jarosik et al. (2003b) for more detail.

Figure 8 in Jarosik et al. (2003b) shows the performance of the gain model when fit to the dipole-based gain solution (see also Figure 7 in this paper). In §3.1 we evaluate the overall performance of the gain model and the hourly gain solution in the context of an end-end simulation designed to place limits on combined calibration and map-making errors. We will conclude that the model provides an excellent description of the radiometer gain, and here we adopt it as the final gain solution for further processing. The gain model fits into the data processing sequence as follows. After we iterate the simultaneous calibration and sky map solution long enough for the calibration to converge (10 iterations when starting with a good sky model), we freeze the dipole-based calibration and fit the gain model parameters in equation (24). All subsequent data products are produced with data calibrated using this gain solution, including the time-ordered archive, the final sky maps, and the Jupiter beam maps.

#### 2.3.2 Baseline Filtering

The baseline that results from the initial calibration solution is not optimal. This is due to the fact that the initial baseline is sampled once per hour (0.28 mHz), while Jarosik et al. (2003b) show that the power spectral density of the noise has a $`1/f`$ knee frequency of a few mHz, typically. If the initial baseline estimate were used in the final sky maps, it would generate weak stripes of correlated noise along the scan paths, as per equation (12). Even so, it is important to note that $`1/f`$ effects are small relative to the white noise. In the worst WMAP radiometer, W41, the amplitude of the noise covariance $`𝐍`$ at small lag is $``$2% of the white noise variance. Thus we treat $`1/f`$ noise iteratively in the data processing by first ignoring it to obtain an estimate of the gain, baseline, and sky solution. Then we subtract the estimated sky signal from the time-ordered data, apply a pre-whitening baseline filter to the residual noise, add the sky signal back in, and write the data to a final calibrated, time-ordered data archive. The approach of first subtracting an estimated sky signal is designed to avoid biasing the gain solution and/or removing low-order power from the sky maps. The noise properties of maps constructed in this way must account for the filtering process. We discuss the map-making algebra appropriate to our filter implementation in Appendix B.

The steps we follow to define and apply the filter are as follows. We remove an estimate of the sky signal, in du, from the raw differential data using

$$𝐜^{}(t)=𝐜(t)𝐠(t)𝚫𝐭^{}(t)b_k,$$
(25)
where $`b_k`$ is the hourly baseline point appropriate to the current time, $`𝐠`$ is the final gain solution from equation (24) and $`𝚫𝐭^{}`$ is the differential sky signal computed from the initial sky map. We then evaluate the auto-correlation function of $`𝐜^{}(t)`$ to a lag of $`10^4`$ sec. Results for representative radiometers are shown in Figure 4. This range of lags is sufficient to account for both the long-range correlations due to $`1/f`$ noise and the correlation at a lag of 1 observation due to the low-pass post-demodulation filter in the Analog Electronics Unit (Bennett et al., 2003a). The baseline filters are then defined as follows (Wright, 1996)

1. Fit the auto-correlation function, $`C(\mathrm{\Delta }t)`$, to the model defined below.
2. Fourier transform the model correlation function to generate the power spectral density $`P(f)`$.
3. Compute the Fourier space filter $`w(f)=1/\sqrt{P(f)}`$ and set $`w(0)=0`$ to produce a zero mean output signal.
4. Fourier transform the filter $`w(f)`$ to generate the time domain filter $`w(t)`$, normalized to 1 at $`t=0`$.

By inspection, the auto-correlation functions are well modeled by a log-linear function

$$C(\mathrm{\Delta }t)/C(0)=\{\begin{array}{cc}C_1\hfill & \mathrm{\Delta }t=\tau \hfill \\ AB\mathrm{log}(\mathrm{\Delta }t/\text{1 s})\hfill & \tau <\mathrm{\Delta }t<10^{A/B}\text{ s}\hfill \\ 0\hfill & \mathrm{\Delta }t>10^{A/B}\text{ s},\hfill \end{array}$$
(26)
where $`\tau `$ is the integration time for a single observation, $`C_1`$ is the correlation at lag $`\tau `$, measured from the data, and $`A`$ and $`B`$ are fit parameters. Note that $`A`$ gives the typical fractional covariance at small lag, while the suppression of correlations at large lag ($``$2000 s) is dictated by the subtraction of the hourly baseline as a pre-filter. The best-fit parameters are given in Table 4, and fits for selected radiometers are shown in Figure 4.

The derived pre-whitening filters, $`w(f)`$, are plotted as a function of frequency in Figure 4. One point of particular interest is the filter response at the spin frequency, 7.7 mHz. As shown in Table 4, the best channels have a transmission of $`>`$ 95%, while the worst case, W41, is just above 35%. These values indicate the amount by which the dipole (calibration) signal would be suppressed if the filter were applied prior to calibration and sky signal subtraction. The convolution of $`𝐜^{}(t)`$ with $`w(t)`$ is performed in Fourier space using the Numerical Recipes routine `convlv` (Press et al., 1992). The number of data points convolved at any one time is chosen to be the smallest power of 2 such that the data span a full day with sufficient padding beyond the day to guarantee that wrap-around effects are negligible. This is $`2^{20}`$ for K–Q bands, and $`2^{21}`$ for V,W bands, which gives a minimum of 2.9 hours of padding on each end of a day. Sample auto-correlation functions obtained from the filtered data are shown in Figure 4. The filtering is clearly effective at removing low frequency noise in the time-ordered data. Another example of filtered data is seen in Figure 5, which shows 1 day of W42 data, one of the worst radiometers for $`1/f`$ noise, before and after filtering. These data are smoothed with a 46 s window to show structure in the unfiltered data since plots of unsmoothed data before and after filtering are virtually indistinguishable.

The above results are encouraging but not definitive, because the process of sky signal subtraction and re-addition could introduce correlated artifacts that these tests would miss. The ultimate test of a filter is its ability to “clean” the pixel-to-pixel covariance matrix of the final sky maps and the noise covariance of the angular power spectrum, without altering the underlying sky signal. The sky map noise properties are discussed in §3.2, while the noise properties of the power spectrum are quantified in Hinshaw et al. (2003b).

#### 2.3.3 Baseline Jumps

Limon et al. (2003) identify 21 instances of sudden baseline jumps, or “glitches”, during WMAP’s first year of operation. These events have been identified as small shifts in the properties of several microwave components resulting from sudden releases of internal mechanical stress, presumably from thermal changes. These events last for less than 1 s, and cause no discernible change in the radiometer gain or noise properties.

Care must be taken in the application of the baseline filter in the vicinity of these steps to avoid ringing in the filtered data. Each event is initially flagged by the pre-processor for at least 1.2 hours on either side of the event. Since the initial hourly baseline is derived entirely in the time domain, the $`\pm `$1.2 hr flagged interval ensures that this baseline estimator only “knows” about data on one side of the jump or the other. Prior to convolving the raw data with the pre-whitening filter, we subtract the initial hourly baseline from the data as a pre-conditioner. Thus all data that is input to the convolution routine has approximately zero mean. On output, the re-processor expands the flagged interval by 0.5 hr on either side of the event to ensure that no edge effects propagate into the usable data. In the first year of operation, a total of 0.13% of the data was lost to these steps. See Table 2 in Bennett et al. (2003b).

The threshold amplitude for jump detection by visual inspection is $``$0.05 du, which corresponds to a jump of $``$150 $`\mu `$K in the calibrated output of radiometer W12, the worst offender. To assess the effect of undetected baseline steps in the data, we have generated a test data set in which we take 24 hours of flight W12 data and insert a step of 0.05 du in each channel. We then run the data through the pre-whitening filter to see the magnitude of the transient response. The result is a transient baseline error with a peak magnitude of 80 $`\mu `$K, which lasts for less than one 2-minute spin period. The total time the baseline error exceeds 10 $`\mu `$K is 22 minutes, or approximately 11 spin periods. We pessimistically assume that there could be as many as 40 steps at or just below the threshold of detection, and that half of these are in W12. If we assume these occur at random times, and note that WMAP observes $``$30% of the sky in any given hour, then any given sky pixel is likely to “see” approximately $`11\times 20\times 0.366`$ data points with baseline errors greater than 10 $`\mu `$K. Since the sign of a given step is random, and since W11 data is combined with W12 in the sky maps, we estimate the residual systematic error in a given pixel of the W1 map is less than $`10\text{ }\mu \text{}/2/\sqrt{20}1`$ $`\mu `$K. We emphasize that no jumps have been observed in any other W band radiometer, thus DA-DA consistency is an excellent test of whether any statistic is sensitive to baseline errors of this nature. We have found no evidence that the W1 map is “out of family” (§3.2).

### 2.4 Final Sky Map Processing

Once the calibrated, time-ordered archive has been written, final sky map processing commences based on the algorithm presented in §2.1. At this stage in the sky map processing, we add a few features to the algorithm that, for simplicity, are not present in the combined calibration and map-making code. These include: 1) Weighting each datum by a true weight, $`1/\sigma _i^2`$ based on an estimate derived from the physical temperature of the instrument cold stage. This introduces $``$1% variations in the data weights over the year, since the instrument noise is a weak function of temperature, and the temperature varies by $``$1% over the year. 2) Accounting for loss imbalance, as discussed in Appendix A. In effect, we model each differential observation as $`𝚫𝐭=(1+x_{im})𝐭(p_A)(1x_{im})𝐭(p_B)`$, where $`x_{im}`$ is the small loss imbalance parameter given by Jarosik et al. (2003b). 3) Computing the planet avoidance flag at run time to reduce the amount of data lost. In the final sky maps, a total of 0.11% of the data was lost to planet avoidance. See Table 2 in Bennett et al. (2003b).

The final stage of sky map processing, based on the filtered data, was run for a total of 20 iterations. Convergence was determined by measuring the rms difference between pairs of iterations for a given differencing assembly. For example, the difference between the 10th and 20th iteration of the W2 sky map is 0.08 $`\mu `$K rms. We estimate that artifacts due to lack of solution convergence are $`<`$0.1 $`\mu `$K rms with all of the power being in the low multipoles, $`l<10`$. We present a final combined estimate of sky map artifacts due to calibration and map-making errors in §3.1. This estimate includes the convergence limits given above.

To assess the effect of the improved calibration model on the final sky maps, we form differences between the final post-filtered maps and the last iteration of the intial, pre-filtered maps. The results for DAs W3 and W4 are shown in Figure 6. The top panels in this Figure show the difference maps from a one-year simulation (§3.1) that included a realistic radiometer noise and gain model. The bottom two panels show differences from the flight data. Because these maps are largely based on the same data, most of the white noise drops out of these differences. The remaining “blobs” of white noise result from the change in the planet cut and can be ignored. The striking feature is the striping present in the W4 difference, but virtually absent in the W3 difference. As we show in Appendix B, this is the structure that has been removed from the data by the pre-whitening filter, an interpretation that is substantiated by the analysis of the simulation. The fact that the W3 difference is very small is an indication that the level of striping in the unfiltered W3 data was very small to start with. We estimate the level of residual striping in the final maps in §3.2. Images of the final maps at each frequency are presented by Bennett et al. (2003b).

## 3 SYSTEMATIC ERROR ANALYSIS

As discussed in the introduction, systematic errors may be classified by the nature of their source. In this section we place limits on the level of systematic errors in the final sky maps, using that classification to guide the analysis. In §3.1 we place limits on combined calibration and map-making artifacts, based largely on a detailed simulation of the first year of WMAP operation. In §3.2 we present null tests based on difference maps formed from a variety of data combinations, each of which should yield no sky signal. We use these maps to measure or place limits on correlated pixel noise (striping) in the final first-year maps. In §3.3 we discuss systematic errors relating to pointing and beam mapping errors. We conclude by placing stringent limits on residual errors due to environmental (thermal and electrical) and other miscellaneous sources.

### 3.1 Calibration and Map-Making Errors

To assess the combined errors from calibration and map-making artifacts, we have generated a high fidelity simulation that includes all of the effects we believe are important for calibration and map-making. In particular this simulation includes: 1) A sky model that closely mimics the statistical properties of the observed sky; 2) A realistic noise model for every channel, including $`1/f`$ noise \[see Jarosik et al. (2003b) for a tabulation of $`1/f`$ knee frequencies\]; 3) A model for the thermal drift of the gain, baseline and offset of each radiometer, based on measured susceptibility coefficients, and driven by the actual temperature profile measured in flight. This simulation generates the sky signal using a circular beam approximation. The effects of elliptical beams are treated in a separate, noiseless simulation in §3.3.4. We write simulated science data to files that mimic the raw telemetry, then process the data using the same pipeline as was used to process the flight data.

The top panel of Figure 7 shows the converged gain solution from the simulation for channel V113; the bottom panel shows the corresponding result from the first year of flight data. In both panels, the “noisy” black traces are the hourly gain data, binned in 24-hr samples to reduce the noise, and the green traces are the best-fit gain model (§2.3.1). For the simulation, the input gain used to generate the data is shown in grey. The absolute gain is recovered in the simulation to better than 0.1% for all 40 channels.

The dipole signal seen by an observer moving with speed $`v`$ relative to the rest frame of the CMB is $`T_0v/c`$, where $`T_0`$ is the absolute temperature of the CMB, and $`c`$ is the speed of light. Thus, additional sources of error that could affect the absolute calibration of the WMAP data include errors in the determination of WMAP’s velocity with respect to the solar system barycenter (the point of reference for the COBE dipole) and errors in the absolute temperature of the CMB. The velocity of WMAP is routinely measured with respect to geocentric inertial coordinates (GCI) with an accuracy of $`<`$1 cm s<sup>-1</sup>. The velocity of the Earth is determined from the JPL ephemeris with similar accuracy. The combined uncertainty from velocity errors is 0.1 nK. The uncertainty in the absolute temperature of the CMB is 0.1% (Mather et al., 1999). Combining these uncertainties with the results of the simulation, we conservatively estimate an absolute calibration error of 0.5% for the first-year WMAP data.

Errors in relative calibration can produce structure in the sky maps, beyond an overall normalization factor. The largest relative discrepancy between the dipole gain solution and the gain model in the flight data is $``$0.4% in K band, and $``$0.2% in the other bands. Similar deviations are seen in the simulation, thus we use the simulation as our primary tool for placing systematic error limits due to relative calibration and map-making errors. We have generated residual maps from the simulated data by subtracting the known input sky signal from the maps produced by the pipeline. These residual maps exhibit no visible structure aside from the pixel noise. In order to assess the errors due to map-making artifacts, we compute the angular power spectrum, $`C_l`$, of the residual maps and search for features in the spectra beyond a simple flat, white noise spectrum. The residual spectra for all 10 differencing assemblies are shown in Figure 8 and summarized in Table 5. In general, the spectra are consistent with white noise over a wide range of multipole moments, but clearly show the most variation at low $`l`$. Because of this, we specifically highlight these modes in Table 5, where we give $`C_2`$, $`C_l_{310}`$, and $`C_l_{11100}`$ for each of the DA’s. For combined systematic error limits due to calibration and map-making, we assign twice the excess variance in each $`l`$ range relative to the white noise plateau, $`\sigma ^{sys}2\left|C_l_{\mathrm{band}}C_l_{7001000}\right|`$. These values are also quoted in Table 5. For comparison, the average power in the CMB in each band is $`C_2130`$ $`\mu `$K<sup>2</sup>, $`C_l_{310}150`$ $`\mu `$K<sup>2</sup>, and $`C_l_{11100}6`$ $`\mu `$K<sup>2</sup>.

Because the simulation includes realistic models of $`1/f`$ noise and long-term thermal effects, these limits also implicitly limit artifacts at low $`l`$ due to these effects. As we demonstrate in subsequent section, we feel this simulation captures all of the important radiometric characteristics of the instrument. Potential artifacts due to optical effects, especially pickup through the far sidelobes, are treated in Barnes et al. (2003), and are summarized in §3.3.5.

### 3.2 Difference Maps and Noise Properties

Difference maps are combinations of the data that, ideally, should contain no sky signal. They provide insight to potential systematic errors and can be used to characterize the noise properties of the sky maps. The first set of difference maps we generate are between DA pairs with the same frequency and beam response, namely $`\frac{1}{2}`$(Q1$``$Q2), $`\frac{1}{2}`$(V1$``$V2), and $`\frac{1}{2}`$(W12$``$W34), where W12 = $`\frac{1}{2}`$(W1+W2), and W34 = $`\frac{1}{2}`$(W3+W4). Images of these difference maps are shown in Figure 9, along with low resolution versions of the sum (signal) maps to give a sense of the signal strength in each map. Aside from the pattern of the noise, which follows the sky coverage \[see Figure 3 of Bennett et al. (2003b)\], the only visible structure in these difference maps is in the Galactic plane, especially in V band. This is understood to be a result of a small difference in the effective center frequency of the V1 and V2 differencing assemblies (Jarosik et al., 2003a). In particular, the V1 map has an effective frequency approximately 1 GHz lower than V2. Since the spectrum of the Galaxy at V band follows $`T_A(\nu )\nu ^2`$ (Bennett et al., 2003c), we expect the Galactic signal to be $``$3% brighter in V1 than V2, which is consistent with the residual signal seen in the difference map. \[A complete tabulation of effective center frequencies, radiometer by radiometer, is given by Jarosik et al. (2003a) for diffuse sources, and by Page et al. (2003a) for point sources.\] Note that because the data are calibrated using the CMB dipole, there should be no residual CMB signal in such a difference map. A more sensitive comparison of the single DA maps is afforded by comparing their angular power spectra. In that case, it is easier to compare across frequencies because differences in beam response are readily accounted for by deconvolution. See Hinshaw et al. (2003a) for such a comparison.

We generate three additional sets of difference maps using different combinations of the 4 channels within a DA. Specifically, we form the differences in the time-ordered data then generate maps as follows

$`{\displaystyle \frac{1}{2}}(𝐝_{13}+𝐝_{14}){\displaystyle \frac{1}{2}}(𝐝_{23}+𝐝_{24})`$ $``$ $`𝚫_{12},`$
$`{\displaystyle \frac{1}{2}}(𝐝_{13}𝐝_{14})+{\displaystyle \frac{1}{2}}(𝐝_{23}𝐝_{24})`$ $``$ $`𝚫_{34},`$
$`{\displaystyle \frac{1}{2}}(𝐝_{13}𝐝_{14}){\displaystyle \frac{1}{2}}(𝐝_{23}𝐝_{24})`$ $``$ $`𝚫_{1234},`$ (27)
where $``$ indicates the map-making process. The $`𝚫_{12}`$ maps are based on the polarization data, but processed as temperature maps, i.e., without attempting to demodulate the polarization signal. Since the two radiometers within a DA have completely independent detection chains, and since the polarization signal is weak to begin with (and is further suppressed by the lack of demodulation) the noise properties of the $`𝚫_{12}`$ maps should be virtually identical to the nominal signal maps. The $`𝚫_{34}`$ and $`𝚫_{1234}`$ maps are based on channel differences, $`(𝐝_{i3}𝐝_{i4})`$, and since the two channels within a radiometer have partially correlated noise, the noise properties of these latter maps will be different than the maps based on $`(𝐝_{i3}+𝐝_{i4})`$. However, these maps do provide a check on the channel calibration, common-mode thermal effects and other potential artifacts.

For each difference map the two-point correlation function and the angular power spectrum are calculated. The results are shown in Figures 10 and 11 and summarized in Table 6. Figure 10 shows the two-point function computed from the $`𝚫_{12}`$ maps for Q2, V2, and W2. The most apparent feature in each of these functions is the slight bump at the beam separation angle of $`\theta _{\mathrm{beam}}141^{}`$, as expected (§2.1); the first data column of Table 6 gives $`C(\theta _{\mathrm{beam}})/C(0)`$, for each DA. Note that, with the exception of K band, the ratio is typically 0.3%. The larger K band values arise because Galactic leakage in these difference maps is most severe in this band. This is also the source of the weak large-scale feature in the V2 two-point function in Figure 10. While this residual signal is small compared to the temperature signal, it is a systematic error that must be accounted for in the analysis of polarization data (Kogut et al., 2003).

Figure 11 shows the angular power spectra of the difference maps for each of the 10 DAs, as well as for the final signal maps. Table 6 summarizes their statistics. Note that the white noise plateau in the 4 channel combinations per DA divide into two families, as noted above, due to the correlations between channels 3 and 4. As a result, the null combinations, $`𝚫_{34}`$ and $`𝚫_{1234}`$, cannot be used to estimate the white noise parameter $`\sigma _0`$ for the signal maps. However, the polarization channel is seen to be in excellent agreement with the temperature channel in the white noise tail, thus, to the extent that real polarization signals, or other systematics, such as bandpass mismatch are not important, these maps should provide an excellent noise model for the temperature data. Table 10 summarizes the shape of the angular power spectrum at low $`l`$ in the same way Table 5 did for the simulation. We find the spectra of these difference maps to be remarkably flat, with residual quadrupole moments of $`<`$4 $`\mu `$K<sup>2</sup> for all bands except K (in which the difference is dominated by bandpass mismatch) and a single combination of W3. This value is much smaller than the small quadrupole measured in our sky (Bennett et al., 2003b). The power in the other multipole ranges is very close to the white noise floor, as seen in the final columns of Table 6. Since the residual signals seen in the flight difference maps are somewhat lower than those seen in the simulation (Table 5), we adopt the more conservative limits from the simulation as systematic error limits for structure at low $`l`$. This allows for the possibility than some of the error seen in the simulation comes from, e.g., common-mode calibration errors that cancel in the difference maps.

The two-point correlation function of the $`𝚫_{12}`$ maps demonstrates that the angle-averaged off-diagonal terms of the pixel-pixel covariance matrix are less than $``$0.3%. However, the maps in Figure 6 indicate the potential for stripes along the scan directions for which the covariance can be locally higher than the angle- averaged value. In order to determine the magnitude of the covariance along the scan directions we perform the following computation. We form W band difference maps: W$`i`$W$`i^{}`$, where $`i`$ = 1-4, and W$`i^{}`$ is the average of the 3 other W band maps, e.g., W1$`{}_{}{}^{}=\frac{1}{3}`$(W2+W3+W4). We then form time-ordered data from this map using the pointing appropriate to DA W$`i`$ and compute the auto-correlation function, $`C(\mathrm{\Delta }t)`$, from 30 days of data. This provides a measure of the pixel-pixel covariance along a stripe. The results for W3 and W4 are shown in Figure 12; the top panels show the covariance, normalized to one at lag zero, computed from the unfiltered maps, while the bottom panels show the results for the filtered maps. In order to convert the time axis to angular displacement along a scan, recall that the 2$`\text{.}^{}`$784 s<sup>-1</sup>spin rate translates to a 2$`\text{.}^{}`$6 s<sup>-1</sup>rate for the beams in either focal plane (the second decimal place depends on position in the focal plane, and time in the precession cycle). The W3 result shows a hint of covariance ($``$0.1%) at lag 0.1 s, or 0$`\text{.}^{}`$26, prior to filtering, but none after ($`<`$0.05%). Prior to filtering, the W4 result shows clear covariance of up to 0.5% at small lag, decaying to $`<`$0.1% at lags of $`10^2`$ s, roughly one full spin. After filtering, the covariance is reduced by nearly a factor of two, but is still clearly detectable. This is understood to be residual covariance that survives the filtering process because of the fact that we subtract an estimated sky signal, based on the data, prior to filtering the noise, then add it back in to restore the signal. The algebra of this process is presented in Appendix B. We emphasize that W4 is the worst DA for $`1/f`$ stripes by at least a factor of 3 (Jarosik et al., 2003b), and we limit covariance along scan directions to be $`<`$0.1% for all other WMAP first-year sky maps.

### 3.3 Pointing and Beam Determination

#### 3.3.1 Spacecraft Attitude Control and Determination

The spacecraft attitude is determined from a combination of two autonomous star trackers (ASTs) with boresights perpendicular to the spin axis (along the spacecraft $`\pm y`$ axes), two rate gyroscopes, and two digital sun sensors. The sensor outputs are combined using a Kalman filter to determine the aspect solution. The sensor noise parameters and offsets were initially calibrated in flight during the in-orbit checkout (IOC) period in July 2001. By the end of IOC, the final tables were uploaded to the spacecraft.

Spacecraft quaternions output by the Kalman filter provide the definitive transformation from the spacecraft reference frame to the J2000 geocentric inertial (GCI) system. Errors in the attitude solution are estimated using the residuals of the individual sensor signals and propagated to the quaternions. After the final Kalman filter parameters were loaded, quaternion differences show a noise-like error with a $`10^{\prime \prime }`$ rms. In addition to the sensor noise, there is an apparent spin-synchronous error of $`10^{\prime \prime }`$ that is believed to be due to propagation errors in the Kalman filter. As discussed below, this error is apparently seen in the instrument boresight determination using Jupiter observations. Since the pointing performance exceeds the requirement of $`0\text{.}^{}9`$ (root-sum-square for three axes), no correction of the spacecraft quaternions is attempted for the first-year processing. Sufficient information exists in the raw telemetry to attempt a correction in the future, if it is warranted. Note that random quaternion errors are automatically accounted for in the flight beam response maps generated from the Jupiter observations (Page et al., 2003a).

#### 3.3.2 Instrument Boresight Determination

As mentioned above, the spacecraft quaternions provide the definitive reference frame for the spacecraft. The instrument boresights, 10 each on the A and B sides, are determined from the Jupiter beam maps (Page et al., 2003a), which are generated with respect to the spacecraft frame provided by the quaternions. The boresight is defined as the location of the peak of a circular Gaussian fit to the main beam. The results of this fitting are given in Table 7 as 20 unit vectors in spacecraft coordinates. These are the values used to determine instrument pointing in the first-year data processing. The uncertainty in the boresight position is $`2^{\prime \prime }`$ per beam in both spacecraft azimuth and elevation. In addition to statistical uncertainty in the boresight fits, there are two other potential sources of error in the boresight determination: changes (drifts) with time, and errors in the relative time-tagging of quaternion data and science data.

To test stability, we note that WMAP is in a position to see Jupiter twice per year for about 45 days each time. We refer to each $``$45 day period as a Jupiter “season”. During the first season of each year, the boresights scan across Jupiter from roughly ecliptic north to ecliptic south, and vice-versa in the second season. As a test of boresight stability, we have generated beam maps from each of the first two seasons of data separately, and have fit boresight directions to each. We find the azimuth positions agree to better than $`3^{\prime \prime }`$ on both the A and B sides, but the elevation positions differ by $`10^{\prime \prime }`$ on the A side, and a smaller amount on the B side. This difference between seasons is consistent with the $`10^{\prime \prime }`$ spin-synchronous error in the spacecraft quaternions discussed above. We ignore this small effect in the first year processing, and subsume the small systematic error that results into our error estimate for the beam transfer functions, as discussed in Page et al. (2003a).

The relative time-tag accuracy of telemetry packets was tested on the ground. A timing computer was set up to simultaneously receive test pulses from both the Attitude Control Electronics (ACE) box and the Digital Electronics Unit (DEU), the two computers that tag the attitude and science data packets, respectively. Each of these boxes in turn derives its time from the main “Mongoose” computer on WMAP (Bennett et al., 2003a). This test demonstrates a relative time-tag accuracy of 30 $`\mu `$s between the quaternion packets and the science packets. In observing mode, the boresights sweep the sky at a rate of $`2\text{.}^{}6`$ s<sup>-1</sup>, so a 30 $`\mu `$s time error produces a negligible pointing error of $`<0.3^{\prime \prime }`$.

Uncertainty in the spacecraft position is another potential source of boresight determination error. For the first year processing we use a predicted ephemeris that is uploaded to the spacecraft approximately weekly for on-board use by the Attitude Control System. This solution is returned in telemetry and is the basis for the ephemeris data supplied with the first-year release of time-ordered data. The uncertainty in these predictions are $`<7`$ km in position and $`<`$1 cm s<sup>-1</sup> in velocity (3 $`\sigma `$), relative to the Earth. An error of 7 km in WMAP’s position would result in error of $``$2 mas in the apparent position of Jupiter as seen from WMAP and is thus completely negligible.

#### 3.3.3 Beam and Window Function Determination

Along with gain calibration and noise properties, knowledge of the beam shapes and window functions are among the most important aspects of the instrument to characterize for accurate measurements of the CMB. Page et al. (2003a) describe in detail the process by which beam maps are formed from in-flight observations of Jupiter, and how those maps are transformed to determine the beam window functions. The primary result they derive are a set of 10 beam transfer functions, $`b_l`$, one per DA, based on azimuthally-averaged beam radial profiles. These transfer functions are included in the first-year data release (Limon et al., 2003). In addition, they derive a full covariance matrix for each transfer function, $`\mathrm{\Sigma }_{ll^{}}^b`$, which characterizes the uncertainty in $`b_l`$. Typically, the uncertainty for a single DA is about 1-2%, with moderate covariance in $`l`$. See Figure 5 of Page et al. (2003a). These estimates already include a systematic error allowance to bound the small range of results obtained from different analysis methods. As described in Hinshaw et al. (2003a), the window function covariance matrices are propagated into the Fisher matrix (inverse covariance matrix) for the final combined angular power spectrum. Thus the final power spectrum, and the parameter fits based on it, already include statistical and systematic window function uncertainties (Spergel et al., 2003; Verde et al., 2003; Peiris et al., 2003).

#### 3.3.4 Effects from Elliptical Beams

The WMAP beams are moderately elliptical, so the use of azimuthally-averaged radial profiles to describe the beam response is an approximation. This approximation is justified in Page et al. (2003a) by noting that the WMAP scan pattern produces excellent azimuthal averaging of the beam response in a large fraction of the sky. They have placed limits on the variation of the window function across the sky by comparing the effective window function in the ecliptic plane, based on a full two-dimensional transform of the beam response averaged over the flight range of scan angles, to the fully averaged transform, $`b_l`$. For the three highest frequency cosmology bands, these variations range from 2-3% at Q band to $``$1% in V and W bands. These variations are consistent with estimates of the angular power spectrum using data at high and low ecliptic latitudes separately (Hinshaw et al., 2003a). Since most of the statistical weight at high-$`l`$ resides in the V and W band data at high ecliptic latitudes, the use of fully averaged beam transforms is appropriate, and the systematic error estimate incorporated into $`\mathrm{\Sigma }_{ll^{}}^b`$ should encompass any error in this approximation.

Elliptical beams can also produce errors in the sky maps that are difficult to characterize in a simple way. We can define the sky map error due to non-circular beam response as

$$𝚫𝐭_{\mathrm{asym}}𝐭_{\mathrm{obs}}𝐭_{\mathrm{circ}},$$
(28)
where $`𝐭_{\mathrm{obs}}`$ is the hypothetical noise-free sky map obtained with the actual experimental beam and scan pattern, and $`𝐭_{\mathrm{circ}}`$ is the ideal sky map obtained by convolving the true sky with the averaged beam transform, $`b_l`$. For a differential experiment like WMAP, there are two effects that contribute to $`𝚫𝐭_{\mathrm{asym}}`$. The first, as noted, is incomplete azimuthal coverage in a given pixel, which gives rise to slightly elliptical peak structure at low ecliptic latitudes (see below), the second is a localized effect due to echoes from bright Galactic sources propagating to other pixels in the map. Specifically, as a bright source is observed in different orientations, the differential signal changes with orientation. Since the map-making algorithm must assign one average value to the pixel with the bright source, the ring of pair pixels at the beam separation will see an echo with a quadrupolar temperature distribution around the ring. We mitigate this effect by incorporating a bright source mask in the map-making algorithm, which is invoked as follows. If side A observes a pixel in the bright source mask, we only update the sky map accumulator for pixel A, but not for pixel B

$`𝐧_{\mathrm{obs}}𝐭_{n+1}(p_A)𝐧_{\mathrm{obs}}𝐭_{n+1}(p_A)+w_i[𝐝(t_i)+𝐭_n(p_B)]`$  $`𝐧_{\mathrm{obs}}(p_A)𝐧_{\mathrm{obs}}(p_A)+w_i`$
$`𝐧_{\mathrm{obs}}𝐭_{n+1}(p_B)\to ̸𝐧_{\mathrm{obs}}𝐭_{n+1}(p_B)w_i[𝐝(t_i)𝐭_n(p_A)]`$  $`𝐧_{\mathrm{obs}}(p_B)\to ̸𝐧_{\mathrm{obs}}(p_B)+w_i,`$ (29)
where the terms are as defined in after equation (19). In this way we obtain an estimate of $`𝐭(p_A)`$, but we do not propagate bright echoes to the ring of neighbor pixels, of which $`p_B`$ is one. The mask we use for assigning this cut is the same Kp8 processing mask we used for the calibration fits (§2.2.1).

We have generated a simulation to quantify the errors from both of these effects. Specifically, the simulation generates one year of noise-free differential sky signal which includes a model for the flight beam ellipticity. We run this data through the flight map-making pipeline to generate sky maps, $`𝐭_{\mathrm{obs}}`$. We also generate convolved maps $`𝐭_{\mathrm{circ}}`$ using the azimuthally-averaged beam transforms appropriate to the beam model. The residual map, $`𝚫𝐭_{\mathrm{asym}}`$, for DA K1 is shown in Figure 13. The K band radiometers have the largest beam ellipticity of all the DAs, so this represents a worst case result. The general “mottling” near the ecliptic plane results from the relatively limited azimuthal coverage in this region producing elliptical peaks and anti-peaks which, in $`𝚫𝐭_{\mathrm{asym}}`$, are differenced with circular counterparts. This is especially noticeable near bright Galactic sources. The rms amplitude of these fluctuations in the Kp2 cut sky is 2 $`\mu `$K in K band, and at least a factor of 2 lower in Q-W bands. The effect of this structure in the power spectrum is primarily represented as a variation in the window function across the sky, as discussed above, and in Page et al. (2003a) and Hinshaw et al. (2003b). However, this structure also contributes to the 4-point function of the data in the sense that it couples power from different $`l`$ ranges. This effect is potentially important for the interpretation of any gravitational lensing analysis of the WMAP data.

The figure also exhibits faint echoes of the brightest sources that evade the map-making cut discussed above. We limit localized features in the Kp2 cut sky to less than 10 $`\mu `$K in K band and less than 2 $`\mu `$K in Q-W bands due to a combination of dimmer sources and more circular beams. We estimate that such features occupy $`<`$0.1% of the Kp2 cut sky.

#### 3.3.5 Far Sidelobe Pick-up

The WMAP optical system was designed to produce minimal pick-up from signals entering the far sidelobes. Barnes et al. (2003) present a complete determination of the WMAP sidelobe response by combining measurements from a variety of ground-based sources with in-flight measurements of the Moon. They produce response maps covering $`4\pi `$ sr that are included as part of the first-year data release. They then use these response maps, with the first-year sky maps, to estimate the systematic artifacts remaining in the first-year maps based on the well-justified assumption that sidelobe artifacts are small relative to the sky signal. The K band data have the largest sidelobe signal due both to the largest sidelobe spill, and to the brightest Galactic signal. The signal was deemed to be large enough, and well enough characterized, to warrant a small post-processing correction to the first-year K band map. Limits on remaining sidelobe induced artifacts in all the bands are presented in Table 1 of Barnes et al. (2003).

### 3.4 Environmental Effects

#### 3.4.1 Thermal Effects

The radiometer gain and offset are dependent on temperature. There are several aspects of the WMAP design that are critical to mitigating this source of systematic error (Bennett et al., 2003a). The instrument is differential, so thermally induced gain changes act on a relatively small offset signal. The observatory environment was designed to be as stable as possible, consistent with other goals. For example, all nominal thermal control is passive to avoid heaters cycling on and off. The observatory is placed at the second Earth-Sun Lagrange point far from the Earth, and the solar panels maintain a fixed $`22\text{.}^{}5`$ angle with respect to the Sun during normal observing mode. The instrument temperature is monitored with precision platinum resistance thermistors (PRTs) to verify the degree to which thermal stability is in fact achieved.

Temperature variations at the spin period are the most critical since they can induce signals that couple relatively efficiently to the sky. However, owing to the relatively fast (129.3 s) spin period and the thermal mass of the instrument, any induced signals will have a very red spectrum and thus will couple only to the lowest few harmonic modes on the sky ($`l\genfrac{}{}{0pt}{}{_<}{^{}}10`$). In the analysis below we use flight data to estimate the susceptibility of the gain and baseline to temperature variations of the instrument. In turn we use limits on the instrument’s physical temperature variation at the spin period from Jarosik et al. (2003b) to put limits on thermally induced artifacts in the time-ordered data, and hence the sky maps.

Thermally induced signals can enter either through changes in the gain acting on the instrument offset or through changes in the offset itself. We show below that the latter are more significant for WMAP.

The radiometer gain model presented by Jarosik et al. (2003b) describes the gain in terms of the RF bias (“total power”) of the detector, and the temperature of the FPA. This model tracks thermal variations in the gain on the time scale of the RF bias readout (23.04 s), and the map-making algorithm updates the gain on this time scale. However, since this sample rate is only a few times per spin, it is possible that a systematic temperature variation at the spin period could induce gain changes that are not well sampled by this model. As a separate check of gain-induced artifacts, we infer the temperature susceptibility of the gain from data taken over a long time period where gain changes are measurable. Results are given in Table 8. We combine these measurements with the limits on temperature modulation at the spin period derived by Jarosik et al. (2003b) to place limits on gain induced artifacts, as shown in Table 9. We conclude that gain-induced signals at the spin period, which might be poorly tracked by the gain model, are $`<`$20 nK.

The instrument baseline is the product of the gain times the offset. As described in §2.2, we get an initial estimate of the baseline from the dipole calibration algorithm. This gives us an estimate of the instrument baseline once per hourly precession period. Sample hourly baselines for channels V113 and V114 are shown as a function of time over the first year of operation in Figure 14. Also shown is the temperature of the instrument FPA over the same time period; there is a clear temperature dependence in the baseline. We measure the baseline temperature susceptibility by fitting the hourly baseline estimates to a model of the form $`𝐛(t)=c_0+c_1t+c_2\mathrm{\Delta }T_{\mathrm{FPA}}(t)`$ where the $`c_i`$ are model coefficients and $`\mathrm{\Delta }T_{\mathrm{FPA}}=T_{\mathrm{FPA}}T_{\mathrm{FPA}}`$ is the deviation of the FPA temperature from its mean. The most robust susceptibility results come from fitting a portion of the data near the time of a partial battery cell failure which occurred on day 2002:054 (GMT) (Limon et al., 2003). In response to this event, the spacecraft bus voltage was autonomously commanded lower on day 2002:058 (GMT) causing the spacecraft to dissipate less power and thus cool slightly. The coefficients are given in Table 8. We have combined the results for the two channels in each radiometer because this is the combination that enters into the final sky maps. This has the effect of canceling some of the common-mode susceptibility measured in individual channels. As noted above, we combine these susceptibility measurements, taken over long time periods, with limits on the temperature variations at the spin period (Jarosik et al., 2003b) to place limits on induced signals at this time scale. The results are given in Table 9. We conclude that offset-induced signals at the spin period are $`<`$180 nK.

Slow drifts in the instrument temperature will be largely filtered out by the baseline pre-whitening discussed in §2.3.2. The steepest temperature gradient observed during the first year of observation occurred just after the above-mentioned battery cell failure. To assess the efficiency of the filtering process, we have analyzed the data during this period as follows. We applied the baseline pre-whitening filter to the temperature signal, $`T_{\mathrm{FPA}}(t)`$, to measure how much the cooling gradient was suppressed by the filter. The input temperature gradient on day 2002:058 (GMT) was $`1.7`$ mK hr<sup>-1</sup>. Applying the K11 filter to $`T_{\mathrm{FPA}}(t)`$ yielded an output gradient of $`0.1`$ $`\mu `$K hr<sup>-1</sup>, while the W41 filter yielded an upper limit 10 times smaller. We conservatively estimate upper limits on residual drift in the filtered baseline of $`<`$10 nK hr<sup>-1</sup>for the most susceptible channels.

#### 3.4.2 Electrical Effects

A variable electrical signal on board the observatory could induce an apparent signal in the radiometers. Sources of variable electrical signals include the reaction wheels, transponder, bus voltage fluctuations, and RF noise coupled to the instrument through the power bus. During the final observatory thermal vacuum test, in which the observatory was operating at temperatures close to those achieved in flight, searches for such electrically induced radiometric artifacts were conducted (Jarosik et al., 2003a). Upper limits on radiometer bus voltage susceptibility, based on ground tests, are given in Table 8. We combine these results with an upper limit on bus voltage variations of 3.0 mV rms, measured on-orbit, to conclude that electrically-induced signals at the spin period are $`<`$40 nK. See Table 9.

### 3.5 Miscellaneous Effects

#### 3.5.1 Radiometer Cross-talk

A large signal in one radiometer could induce an erroneous signal at the output of another radiometer due to electrical cross talk. Such cross talk is not expected but could arise from, e.g., non-ideal amplifier behavior, or other parasitic effects, such as pickup in the wiring harnesses.

A careful search was made for this effect during the instrument ground tests. Noise diodes were used to inject a large signal into one radiometer at a time while the input feeds of all other radiometers were covered by absorptive loads. The outputs of the 9 non-driven differencing assemblies were searched for any evidence of the injected signal. The tests were run with the amplifiers in the passive channels both on and off in order to distinguish pickup mechanisms. No pickup was found in any test. Table 10 gives $`2\sigma `$ upper limits on the pickup by any DA due to any of the 9 other DAs. The column labeled Electrical gives the results obtained from the test with the amplifiers turned off, and the column labeled Radiometric gives weaker limits from the test with the HEMTs turned on. The latter limits are weaker because the output noise levels are higher.

This level of pickup is quite small and could only be of potential concern when WMAP scans across Jupiter, the brightest source in the sky for WMAP at L2. The values in the Table 10 are thus referred to peak Jupiter signals in each band. For example, the first entry indicates that when Jupiter induces a signal of 185 mK in the W1 radiometer, the pickup in the K1 differencing assembly is $`<`$30 $`\mu `$K (95% confidence) which is $``$26.8 dB below the peak Jupiter signal of 14 mK in K band. This signal occurs when the beam of the pickup channel is within a few degrees of Jupiter, depending on channel separation in the focal plane, and in every case is less than direct radiometric detection of Jupiter in the near side lobes.

This limit on cross talk implies pickup is completely inconsequential in normal observing mode. Using the same Table example, a 100 $`\mu `$K signal in W1 could cause at most a 200 nK signal in the most susceptible the four K band differencing assemblies.

#### 3.5.2 Source Variability

Time variable objects are a potential source of contamination for observations of the CMB; see for example Sokasian et al. (2001) and references therein. One concern is that an object may grow in brightness over the course of WMAP’s observations, avoid detection during an initial source survey, and remain unmasked during subsequent data analysis. For example, blazars produce relatively rapid and large amplitude variability in all wavebands. Long term observations of such objects show that increases in flux by a factor of up to a few over a time scale of years can be anticipated (Flett & Henderson, 1983; Ennis et al., 1982; Stevens et al., 1994; Bower et al., 1997). Observations of Zw 2 by Falcke et al. (1999), provide an extreme example: a greater than 20-fold increase in brightness, from $``$0.1 Jy to $``$2 Jy, over a period of less than two years. While this object could produce a temperature response of a few hundred $`\mu `$K in the WMAP data, such objects are rare and, if left undetected, would have a minimal effect on cosmological inferences. Tests for point source contamination in the WMAP data are given by Bennett et al. (2003c) and Hinshaw et al. (2003b). These tests will need to be revisited on an annual basis.

Another source of concern is that a time variable source in the WMAP data has the effect of broadcasting noise to the ring of $``$1000 pixels which are separated from the variable source by the beam separation angle ($`\theta _{\mathrm{beam}}141^{}`$). The point source list derived from the WMAP first-year data is 98% reliable with $``$5 spurious detections at the $``$0.5 Jy flux limit of the survey (Bennett et al., 2003c). The nominal point source sensitivity of the WMAP telescope is $`\mathrm{\Gamma }200`$ $`\mu `$K/Jy, thus a noise level of $``$0.1 $`\mu `$K is expected from variable sources that evade detection.

## 4 CONCLUSIONS

The processing steps used to produce the first-year WMAP sky maps include an initial simultaneous estimate of the sky map and the instrument calibration. The instrument gain is then refined using a model based on engineering telemetry, and the instrument baseline is refined by the application of a pre-whitening filter. A final archive of calibrated data is produced and used to generate final sky maps using a slightly refined iterative algorithm.

We limit systematic artifacts due to calibration, map-making and environmental disturbances to less than $``$15 $`\mu `$K<sup>2</sup> in the quadrupole $`C_2`$, with tighter limits at higher multipole moments (Table 5). Beam transfer functions are measured for each beam with 1-3% over the entire range of multipole moments that WMAP is sensitive to (Page et al., 2003a). The covariance matrix of the beam transfer function is propagated through to the final power spectrum error matrix. We characterize pixel-pixel covariance matrix and place limits on residual stripes in the final maps.

All major data products from the first year of WMAP observations are being released through NASA’s new Legacy Archive for Microwave Background Data Analysis (LAMBDA) at `http://lambda.gsfc.nasa.gov/`.

The WMAP mission is made possible by the support of the Office of Space Sciences at NASA Headquarters and by the hard and capable work of scores of scientists, engineers, technicians, machinists, data analysts, budget analysts, managers, administrative staff, and reviewers. We acknowledge use of the HEALPix package.

## Appendix A THE MAPPING FUNCTION

Equation (4) defines the continuous form of the mapping function, which encodes both the scan strategy of an experiment, and convolution due to the beam response. We can relate this to the matrix form, in equation (7) as follows. The mapping function evaluated at time $`t_i`$ for a finite integration time $`\tau `$ may be written in terms of the beam response function as

$$𝐌(𝐧,t_i)=\frac{1}{\tau }_{t_i}^{t_i+\tau }𝑑t\left[\alpha B_A(𝐑(t)𝐧)\beta B_B(𝐑(t)𝐧)\right],$$
(A1)
where $`B_A(𝐧)`$ is the beam response of the A-side beam, in spacecraft coordinates, normalized to unit integral

$$𝑑\mathrm{\Omega }_𝐧B_A(𝐧)1,$$
(A2)
similarly for the B side, and $`𝐑(t)`$ is the rotation matrix from sky-fixed (Galactic) coordinates to spacecraft coordinates at time $`t`$. The terms $`\alpha `$ and $`\beta `$ in equation (A1) account for possible ohmic losses in the A and B-side optics that are not necessarily equal (Jarosik et al., 2003b). Since the data are calibrated using the modulation of the CMB dipole, we only need to parameterize the loss imbalance which, following Jarosik et al. (2003b), we parameterize as

$`\alpha `$ $``$ $`1+x_{\mathrm{im}},`$
$`\beta `$ $``$ $`1x_{\mathrm{im}}.`$ (A3)
Note that loss imbalance is separate from lossless differences in the beam response function, e.g. differences in the solid angle of the A and B-side beams. Once the calibration is applied, the differential sky signal is a measurement of the form

$$𝚫𝐭(t)=𝑑\mathrm{\Omega }_𝐧𝐭(𝐧)\left[(1+x_{\mathrm{im}})B_A(𝐑(t)𝐧)(1x_{\mathrm{im}})B_B(𝐑(t)𝐧)\right],$$
(A4)
which still includes the effects of any loss imbalance. We now separately consider how this calibrated differential data propagates into the sky maps and the Jupiter beam maps.

When making sky maps from the calibrated data, each datum is modeled simply as

$$𝚫𝐭(t_i)=(1+x_{\mathrm{im}})𝐭(p_A)(1x_{\mathrm{im}})𝐭(p_B),$$
(A5)
where $`p_A`$ is the pixel observed by the A-side beam at time $`t_i`$, and similarly for $`p_B`$. That is, each row of the mapping matrix in equation (7) has the form

$$𝐌(p,t_i)=[\mathrm{},(1+x_{\mathrm{im}}),\mathrm{},(1x_{\mathrm{im}}),\mathrm{}],$$
(A6)
with non-zero entries in pixel columns $`p_A`$ and $`p_B`$ only. Upon solving for the sky map, this ideally leads to an effective beam response of the form

$$B(𝐧)=\frac{1}{2}\left[\frac{B_A^{(s)}(𝐧)}{(1+x_{\mathrm{im}})}+\frac{B_B^{(s)}(𝐧)}{(1x_{\mathrm{im}})}\right],$$
(A7)
where $`B_A^{(s)}`$ is the symmetrized beam response for the A-side beam and similarly for the B-side. For this ideal case to obtain, the following must hold: 1) each pixel must be observed equally by the A and B-side beams, which is true to a very good approximation for WMAP, and 2) each pixel must be observed with uniform azimuthal coverage. Deviations from these assumptions are considered in the text.

The beam mapping data is compiled from calibrated observations of the bright source Jupiter. The calibrated data taken when side A is observing Jupiter has the form

$$𝚫𝐭=𝑑\mathrm{\Omega }_𝐧𝐭_J(𝐧)(1+x_{\mathrm{im}})B_A(𝐑𝐧)+𝚫𝐭_{\mathrm{sky}},$$
(A8)
where $`𝐭_J(𝐧)`$ is the brightness temperature of Jupiter in the direction $`𝐧`$ and $`𝚫𝐭_{\mathrm{sky}}`$ is the background sky temperature difference, which is subtracted during processing. An analogous equation holds for the B-side data. Assuming the beam response is constant over the extent of Jupiter, the integral reduces to

$$𝚫𝐭=T_J\mathrm{\Omega }_J(1+x_{\mathrm{im}})B_A(𝐑𝐧_J),+𝚫𝐭_{\mathrm{sky}},$$
(A9)
where $`T_J`$ is the disk brightness temperature of Jupiter, and $`\mathrm{\Omega }_J`$ is its solid angle. Beam maps are compiled by binning the corrected data $`𝚫𝐭𝚫𝐭_{\mathrm{sky}}`$ as a function of $`𝐧`$, in spacecraft coordinates. This produces maps proportional to the beam response

$$T_A(𝐧)=T_J\mathrm{\Omega }_J(1+x_{\mathrm{im}})B_A(𝐧),$$
(A10)
and similarly for the B-side.

Ultimately, we wish to compute the transfer function of the symmetrized beam response. This may be obtained from the symmetrized beam maps as

$$B(𝐧)=\frac{1}{2T_J\mathrm{\Omega }_J}\left[\frac{T_A^{(s)}(𝐧)}{(1+x_{\mathrm{im}})}+\frac{T_B^{(s)}(𝐧)}{(1x_{\mathrm{im}})}\right].$$
(A11)
We don’t know the brightness temperature of Jupiter a priori, but since this is an overall normalization factor, we are free to normalize the final transfer function to 1 at $`l=0`$.

## Appendix B MAP-MAKING WITH FILTERED DATA

§2.3.2 presents the filtering algorithm used to determine the final instrument baseline. This process included an estimated sky signal subtraction based on the initial sky maps produced with the hourly calibration. In this Appendix we derive the noise properties of sky maps produced with this filtered data. In the following, we assume the time-ordered data has a noise covariance $`𝐍=\mathrm{𝐧𝐧}^T`$ that includes a $`1/f`$ component, and that we have a pre-whitening filter $`𝐅`$ such that

$$(𝐅𝐧)(𝐅𝐧)^T𝐈.$$
(B1)
### B.1 Map-making with filtered signal + noise

We could filter the full data prior to any calibration or sky map estimation, then deconvolve the effects of the filter in the subsequent data processing. The input data would have the form

$$𝐝^{}=𝐅𝐝=𝐅\mathrm{𝐌𝐭}+𝐅𝐧.$$
(B2)
Then, in order to obtain an unbiased sky map estimate, we would need to evaluate the sky map estimator

$$𝐭^{}=(𝐌^T𝐅^T\mathrm{𝐅𝐌})^1(𝐌^T𝐅^T𝐝^{}),$$
(B3)
which deconvolves the action of the filter on the sky signal. Since the $`1/f`$ noise in the WMAP data is relatively small, implementing this estimator for the first-year sky maps was deemed unnecessary and would likely have delayed the release of the maps. The alternative is to filter only the noise by subtracting an estimate of the sky signal prior to filtering, then adding it back in to the time-ordered data prior to making new maps.

### B.2 Map-making with filtered noise

Let $`𝐭_0`$ be the sky map estimated from unfiltered data, using the hourly calibration. This is related to the true sky signal by

$`𝐭_0`$ $`=`$ $`\mathrm{𝐖𝐝}`$ (B4)
$`=`$ $`𝐖(\mathrm{𝐌𝐭}+𝐧)`$ (B5)
$`=`$ $`𝐭+\mathrm{𝐖𝐧}`$ (B6)
where $`𝐖=(𝐌^T𝐌)^1𝐌^T`$ is the map-making operator defined in §2.1, and we have used the fact that $`𝐖𝐌=1`$. We use this sky map to subtract a sky signal from the time-ordered data prior to filtering, then we add it back in after filtering. This produces a filtered data set

$`𝐝_1`$ $`=`$ $`𝐅(𝐝\mathrm{𝐌𝐭}_0)+\mathrm{𝐌𝐭}_0`$ (B7)
$`=`$ $`𝐅(\mathrm{𝐌𝐭}+𝐧\mathrm{𝐌𝐭}\mathrm{𝐌𝐖𝐧})+\mathrm{𝐌𝐭}+\mathrm{𝐌𝐖𝐧}`$ (B8)
$`=`$ $`\mathrm{𝐌𝐭}+\mathrm{𝐅𝐧}+(\mathrm{𝟏}𝐅)\mathrm{𝐌𝐖𝐧}.`$ (B9)
This time series data set consists of an unbiased sky signal $`\mathrm{𝐌𝐭}`$, a white noise term $`\mathrm{𝐅𝐧}`$, and a residual noise term $`(\mathrm{𝟏}𝐅)\mathrm{𝐌𝐖𝐧}`$ which is due to the off-diagonal “wings” of the filter $`(\mathrm{𝟏}𝐅)`$ acting on the noise from the initial sky map estimate, $`\mathrm{𝐌𝐖𝐧}`$.

We can make a map $`𝐭_1`$ from the data $`𝐝_1`$ using the algorithm of §2.1

$`𝐭_1`$ $`=`$ $`\mathrm{𝐖𝐝}_1`$ (B10)
$`=`$ $`𝐭+\mathrm{𝐖𝐅𝐧}+𝐖(\mathrm{𝟏}𝐅)\mathrm{𝐌𝐖𝐧}`$ (B11)
$`=`$ $`𝐭+\mathrm{𝐖𝐅𝐧}+(\mathrm{𝟏}\mathrm{𝐖𝐅𝐌})\mathrm{𝐖𝐧}`$ (B12)
$`=`$ $`𝐭+\mathrm{𝐖𝐅𝐧}+𝐑\mathrm{𝐖𝐧},`$ (B13)
where we have again used $`𝐖𝐌=1`$ and we have defined a “residual” operator $`𝐑(\mathrm{𝟏}\mathrm{𝐖𝐅𝐌})`$ which is small in the sense that only off-diagonal terms in $`𝐅`$ contribute to it. This is most easily seen if we recall that the filter operator is 1 on the diagonal and has small off-diagonal terms. We can then write $`𝐅\mathrm{𝟏}𝐄`$ from which it follows $`𝐑=\mathrm{𝐖𝐄𝐌}`$. It follows that $`𝐭_1`$ is an unbiased estimate of $`𝐭`$ that includes a white noise term $`\mathrm{𝐖𝐅𝐧}`$ (this noise still contains the small beam separation covariance), and a residual noise term, $`𝐑\mathrm{𝐖𝐧}`$, due to the noisy sky signal estimator used in the filtering process. This latter term is the “excess” noise seen in the W band single DA maps after filtering (Hinshaw et al., 2003a).

The residual noise term can be reduced somewhat by iterating the filter algorithm a second time, using $`𝐭_1`$ as a sky signal estimator. It is straightforward to show that, after some algebra, the resulting map is

$$𝐭_2=𝐭+\mathrm{𝐖𝐅𝐧}+𝐑\mathrm{𝐖𝐅𝐧}+𝒪(𝐑^2),$$
(B14)
where we have neglected a term of order $`𝐑^2\mathrm{𝐖𝐧}`$. The residual noise is reduced slightly since $`𝐑`$ is now acting on the white noise $`\mathrm{𝐖𝐅𝐧}`$ instead of the full noise $`\mathrm{𝐖𝐧}`$. But since the white noise dominates, this is a relatively insignificant improvement. It is also clear that subsequent iterations of the filtering only contribute higher order corrections which are negligible. This convergence has been verified with the flight data. Another feature we have verified with the flight data is the fact that the excess noise term decreases with time. The reason for this is simply that $`𝐑=(\mathrm{𝟏}\mathrm{𝐖𝐅𝐌})=\mathrm{𝐖𝐄𝐌}`$ gets smaller with additional data because the map-making operator $`𝐖`$ gets smaller as more observations accumulate in its “denominator”, $`(𝐌^T𝐌)^1`$. The first-year maps were only filtered once, as per equation (B13), because the improvement in noise properties was not deemed sufficient to warrant the additional processing time.

## Appendix C MAP-MAKING WITH POLARIZATION

WMAP observes the sky with two orthogonal linear polarization modes per feed, thus it is sensitive to the 3 Stokes parameters $`I`$, $`Q`$, and $`U`$. This Appendix outlines the algorithm with which these parameters can be mapped using the differential data from WMAP. The approach is an extension of the iterative method in §2.1.1 introduced by Wright et al. (1996).

### C.1 Polarization Mapping with Total Power Data

Suppose we observed the sky with a single beam, total power radiometer that is sensitive to a single linear polarization, denoted mode #1. In terms of the Stokes parameters, the temperature observed by the instrument at time $`t`$ in pixel $`p`$ would be

$$𝐝_1(t)=𝐢(p)+𝐪(p)\mathrm{cos}2\gamma +𝐮(p)\mathrm{sin}2\gamma ,$$
(C1)
where $`𝐢`$, $`𝐪`$, and $`𝐮`$ are Stokes parameter maps in units of temperature, and $`\gamma `$ is the angle between the polarization axis of the beam and the chosen reference direction for pixel $`p`$ (the choice of reference direction is discussed below). Note that we adopt the common convention that the instrument response reduces to the total temperature in the limit of an unpolarized source.

The signal in the orthogonal polarization channel, which is fed by the other port of the ortho-mode transducer (OMT) and denoted mode #2, is given by

$`𝐝_2(t)`$ $`=`$ $`𝐢(p)+𝐪(p)\mathrm{cos}2(\gamma +{\displaystyle \frac{\pi }{2}})+𝐮(p)\mathrm{sin}2(\gamma +{\displaystyle \frac{\pi }{2}})`$ (C2)
$`=`$ $`𝐢(p)𝐪(p)\mathrm{cos}2\gamma 𝐮(p)\mathrm{sin}2\gamma .`$ (C3)
We can isolate the intensity and polarization signals by taking sums and differences

$`𝐝(t)`$ $``$ $`{\displaystyle \frac{1}{2}}\left(𝐝_1+𝐝_2\right)=𝐢(p)`$ (C4)
$`𝐩(t)`$ $``$ $`{\displaystyle \frac{1}{2}}\left(𝐝_1𝐝_2\right)=𝐪(p)\mathrm{cos}2\gamma +𝐮(p)\mathrm{sin}2\gamma .`$ (C5)
Given noisy data, we can estimate the intensity $`𝐢(p)`$ by averaging all the data $`𝐝(t_i)`$. For the polarization we can only estimate a linear combination of $`𝐪`$ and $`𝐮`$ from a single observation. However, if we have several observations of pixel $`p`$ with a variety of polarization angles $`\gamma `$, we can estimate $`𝐪`$ and $`𝐮`$ in a given pixel by minimizing $`\chi ^2`$, defined as

$$\chi ^2\underset{ip}{}\frac{\left[𝐩(t_i)𝐪(p)\mathrm{cos}2\gamma _i𝐮(p)\mathrm{sin}2\gamma _i\right]^2}{\sigma _i^2},$$
(C6)
where $`i`$ is a time-ordered data index and the sum is over observations within pixel $`p`$, and $`\gamma _i`$ is the polarization angle for the $`i^{\mathrm{th}}`$ observation. The best-fit values for $`𝐪`$ and $`𝐮`$ are given by

$$\left(\begin{array}{c}𝐪(p)\\ 𝐮(p)\end{array}\right)=\frac{1}{\mathrm{\Delta }}\left(\begin{array}{cc}_is_i^2/\sigma _i^2& _ic_is_i/\sigma _i^2\\ _ic_is_i/\sigma _i^2& _ic_i^2/\sigma _i^2\end{array}\right)\left(\begin{array}{c}\frac{1}{2}_ic_i𝐩(t_i)/\sigma _i^2\\ \frac{1}{2}_is_i𝐩(t_i)/\sigma _i^2\end{array}\right).$$
(C7)
where $`c_i\mathrm{cos}2\gamma _i`$, $`s_i\mathrm{sin}2\gamma _i`$, and $`\mathrm{\Delta }_ic_i^2/\sigma _i^2_is_i^2/\sigma _i^2\left(_ic_is_i/\sigma _i^2\right)^2`$ is the determinant of the normal equations matrix. The standard errors for $`𝐪`$ and $`𝐮`$ are given by the inverse of the normal equations matrix

$`\sigma _𝐪^2`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }}}{\displaystyle \underset{i}{}}s_i^2/\sigma _i^2,`$ (C8)
$`\sigma _𝐮^2`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }}}{\displaystyle \underset{i}{}}c_i^2/\sigma _i^2,`$ (C9)
$`\sigma _{\mathrm{𝐪𝐮}}^2`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }}}{\displaystyle \underset{i}{}}c_is_i/\sigma _i^2.`$ (C10)
In the limit of uniform azimuthal coverage and constant noise per observation ($`\sigma _i=\sigma _0`$), the matrix elements in the linear system reduce to

$$\left(\begin{array}{cc}_ic_i^2& _ic_is_i\\ _ic_is_i& _is_i^2\end{array}\right)\left(\begin{array}{cc}N/2& 0\\ 0& N/2\end{array}\right),$$
(C11)
where $`N`$ is the number of observations of pixel $`p`$. In this limit, the noise in $`𝐪`$ and $`𝐮`$ is equal and uncorrelated and reduces to

$$\sigma _𝐪=\sigma _𝐮\sqrt{\frac{2}{N}}\sigma _0.$$
(C12)
Thus the noise in each polarization component is $`\sqrt{2}`$ times noisier than in the intensity $`𝐢`$.

### C.2 Polarization Mapping with Differential Data

We now generalize to the case of polarization mapping with differential input data. For clarity, we first consider the case where the loss in the two sides A and B are equal. We generalize to the case with unbalanced loss in the next subsection. In the case of WMAP, one radiometer in a differencing assembly (radiometer #1) is fed from the axial port of the OMT and the other (#2) is fed by the lateral port (Jarosik et al., 2003a). Following equations (C1,C3) the differential signals from radiometers 1 and 2 are

$`𝐝_1(t)={\displaystyle \frac{1}{2}}(𝐝_{13}+𝐝_{14})`$ $`=`$ $`𝐢(p_A)+𝐪(p_A)\mathrm{cos}2\gamma _A+𝐮(p_A)\mathrm{sin}2\gamma _A`$ (C13)
$``$ $`𝐢(p_B)𝐪(p_B)\mathrm{cos}2\gamma _B𝐮(p_B)\mathrm{sin}2\gamma _B,`$ (C14)
and

$`𝐝_2(t)={\displaystyle \frac{1}{2}}(𝐝_{23}+𝐝_{24})`$ $`=`$ $`𝐢(p_A)𝐪(p_A)\mathrm{cos}2\gamma _A𝐮(p_A)\mathrm{sin}2\gamma _A`$ (C15)
$``$ $`𝐢(p_B)+𝐪(p_B)\mathrm{cos}2\gamma _B+𝐮(p_B)\mathrm{sin}2\gamma _B,`$ (C16)
where $`\gamma _A`$ is the angle between the axial polarization plane and the reference direction in the pixel seen by the A beam, and similarly for $`\gamma _B`$. We take sums and differences of the two signals to isolate the unpolarized and polarized portions of the signal

$$𝐝(t)\frac{1}{2}(𝐝_1+𝐝_2)=𝐢(p_A)𝐢(p_B)$$
(C17)
$$𝐩(t)\frac{1}{2}(𝐝_1𝐝_2)=𝐪(p_A)\mathrm{cos}2\gamma _A+𝐮(p_A)\mathrm{sin}2\gamma _A𝐪(p_B)\mathrm{cos}2\gamma _B𝐮(p_B)\mathrm{sin}2\gamma _B.$$
(C18)
An iterative scheme for making maps of $`𝐪`$ and $`𝐮`$ follows the form used for intensity maps (§2.1.1). Let $`𝐪_n`$ and $`𝐮_n`$ be the $`n^{\mathrm{th}}`$ estimates of $`𝐪`$ and $`𝐮`$ respectively. Estimates for $`𝐪_{n+1}`$ and $`𝐮_{n+1}`$ are obtained by combining the per-pixel fitting algorithm in equation (C6) with the iterative map-making algorithm, as follows

$$\chi ^2\underset{ip}{}\frac{\left[𝐩^{}(t_i)𝐪_{n+1}(p)\mathrm{cos}2\gamma _i𝐮_{n+1}(p)\mathrm{sin}2\gamma _i\right]^2}{\sigma _i^2},$$
(C19)
where the sum is over all observations of pixel $`p`$ by either the A- or B-side beam, and $`𝐩^{}(t)`$ is the polarization data corrected with an estimate of the signal in the opposite beam

$$𝐩^{}(t_i)\{\begin{array}{cc}𝐩(t_i)+𝐪_n(p_B)\mathrm{cos}2\gamma _B+𝐮_n(p_B)\mathrm{sin}2\gamma _B\hfill & \text{beam }Ap\hfill \\ 𝐩(t_i)+𝐪_n(p_A)\mathrm{cos}2\gamma _A+𝐮_n(p_A)\mathrm{sin}2\gamma _A\hfill & \text{beam }Bp.\hfill \end{array}$$
(C20)
The best-fit solution for $`𝐪_{n+1}`$ and $`𝐮_{n+1}`$ is then

$$\left(\begin{array}{c}𝐪_{n+1}(p)\\ 𝐮_{n+1}(p)\end{array}\right)=\frac{1}{\mathrm{\Delta }}\left(\begin{array}{cc}_is_i^2/\sigma _i^2& _ic_is_i/\sigma _i^2\\ _ic_is_i/\sigma _i^2& _ic_i^2/\sigma _i^2\end{array}\right)\left(\begin{array}{c}_ic_i𝐩^{}(t_i)/\sigma _i^2\\ _is_i𝐩^{}(t_i)/\sigma _i^2\end{array}\right),$$
(C21)
where the sum on $`i`$ is as defined above. The uncertainties are as given for the total power case, where $`\sigma _i`$ is now the uncertainty per differential observation, $`𝐩(t_i)`$.

### C.3 Map Making with Unbalanced Differential Data

We now generalize to the case of map-making with unbalanced differential input data. In this case, losses in the A and B-side front ends are different and the differential signal is of the form

$`𝐝_1(t)`$ $`=`$ $`(1+x_{\mathrm{im},1})\left[𝐢(p_A)+𝐬(p_A,\gamma _A)\right](1x_{\mathrm{im},1})\left[𝐢(p_B)+𝐬(p_B,\gamma _B)\right]`$ (C22)
$`𝐝_2(t)`$ $`=`$ $`(1+x_{\mathrm{im},2})\left[𝐢(p_A)𝐬(p_A,\gamma _A)\right](1x_{\mathrm{im},2})\left[𝐢(p_B)𝐬(p_B,\gamma _B)\right],`$ (C23)
where $`x_{\mathrm{im},1},x_{\mathrm{im},2}`$ are the loss imbalance parameters in radiometers 1 and 2, respectively, as defined in Jarosik et al. (2003b), and $`𝐬(p,\gamma )`$ is short-hand for the linear combination of Stokes parameters

$$𝐬(p,\gamma )𝐪(p)\mathrm{cos}2\gamma +𝐮(p)\mathrm{sin}2\gamma .$$
(C24)
As before, we take sums and differences of the two signals to isolate the unpolarized and polarized portions of the signal. First, define the mean imbalance and the “imbalance in the imbalance” as

$`\overline{x}_{\mathrm{im}}`$ $``$ $`{\displaystyle \frac{1}{2}}(x_{\mathrm{im},1}+x_{\mathrm{im},2})`$ (C25)
$`\delta x_{\mathrm{im}}`$ $``$ $`{\displaystyle \frac{1}{2}}(x_{\mathrm{im},1}x_{\mathrm{im},2}),`$ (C26)
then

$`𝐝(t)={\displaystyle \frac{1}{2}}(𝐝_1+𝐝_2)`$ $`=`$ $`𝐢(p_A)𝐢(p_B)+\overline{x}_{\mathrm{im}}\left[𝐢(p_A)+𝐢(p_B)\right]+\delta x_{\mathrm{im}}\left[𝐬(p_A,\gamma _A)+𝐬(p_B,\gamma _B)\right]`$ (C27)
$`=`$ $`(1+\overline{x}_{\mathrm{im}})𝐢(p_A)(1\overline{x}_{\mathrm{im}})𝐢(p_B)+𝒪(\delta x_{\mathrm{im}})`$
$`𝐩(t)={\displaystyle \frac{1}{2}}(𝐝_1𝐝_2)`$ $`=`$ $`𝐬(p_A,\gamma _A)𝐬(p_B,\gamma _B)+\overline{x}_{\mathrm{im}}\left[𝐬(p_A,\gamma _B)+𝐬(p_B,\gamma _B)\right]+\delta x_{\mathrm{im}}\left[𝐢(p_A)+𝐢(p_B)\right]`$ (C28)
$`=`$ $`(1+\overline{x}_{\mathrm{im}})𝐬(p_A,\gamma _A)(1\overline{x}_{\mathrm{im}})𝐬(p_B,\gamma _B)+𝒪(\delta x_{\mathrm{im}}).`$
Note that the term of $`𝒪(\delta x_{\mathrm{im}})`$ in equation (C27) is negligible, because $`𝐢𝐬`$, but the term in equation (C28) must be considered more carefully. First note that $`\delta x_{\mathrm{im}}`$ is small – from Table 3 in Jarosik et al. (2003b), the largest value is 0.35% in W2, with most values being $``$0.1%. Second, the multiplier, $`(𝐢(p_A)+𝐢(p_B))`$ does not modulate with polarization angle, $`\gamma `$, thus it is effectively an offset term that is further suppressed by the map-making demodulation. Finally, the term is out of phase with the differential signal $`(𝐢(p_A)𝐢(p_B))`$ so it does not effectively couple to the sky. The effect of this term in the first-year data is further considered by Kogut et al. (2003).

We generalize the differential map-making algorithm to account for loss imbalance as follows. For intensity, let $`𝐢_n`$ be the $`n^{\mathrm{th}}`$ estimate of $`𝐢`$, then

$$𝐝^{}(t_i)\{\begin{array}{cc}\left[+𝐝(t_i)+(1\overline{x}_{\mathrm{im}})𝐢_n(p_B)\right]/(1+\overline{x}_{\mathrm{im}})\hfill & \text{beam A }p\hfill \\ \left[𝐝(t_i)+(1+\overline{x}_{\mathrm{im}})𝐢_n(p_A)\right]/(1\overline{x}_{\mathrm{im}})\hfill & \text{beam B }p.\hfill \end{array}$$
(C29)
The updated intensity map is then estimated by binning the corrected data

$$𝐢_{n+1}(p)=\frac{_iw_i𝐝^{}(t_i)}{_iw_i},$$
(C30)
where $`w_i`$ is the normalized statistical weight of each observation

$$w_i=\{\begin{array}{cc}(1+\overline{x}_{\mathrm{im}})^2\sigma _0^2/\sigma _i^2\hfill & \text{beam A }p\hfill \\ (1\overline{x}_{\mathrm{im}})^2\sigma _0^2/\sigma _i^2\hfill & \text{beam B }p.\hfill \end{array}$$
(C31)
For polarization, let $`𝐪_n`$ and $`𝐮_n`$ be the $`n^{\mathrm{th}}`$ estimates of $`𝐪`$ and $`𝐮`$ respectively, and define $`𝐬_n(p,\gamma )𝐪_n(p)\mathrm{cos}2\gamma +𝐮_n(p)\mathrm{sin}2\gamma `$. Estimates for $`𝐪_{n+1}`$ and $`𝐮_{n+1}`$ are obtained by combining the per-pixel $`𝐪`$ and $`𝐮`$ demodulation with the iterative map-making algorithm. We define $`\chi ^2`$ as follows

$$\chi ^2\underset{ip}{}\frac{\left[𝐩^{}(t_i)𝐪_{n+1}(p)\mathrm{cos}2\gamma _i𝐮_{n+1}(p)\mathrm{sin}2\gamma _i\right]^2}{\sigma _i^2},$$
(C32)
where the sum is over all observations of pixel $`p`$ by either the A- or B-side beam, and

$$𝐩^{}(t_i)\{\begin{array}{cc}\left[+𝐩(t_i)+(1\overline{x}_{\mathrm{im}})𝐬_n(p_B,\gamma _B)\right]/(1+\overline{x}_{\mathrm{im}})\hfill & \text{beam A }p\hfill \\ \left[𝐩(t_i)+(1+\overline{x}_{\mathrm{im}})𝐬_n(p_A,\gamma _A)\right]/(1\overline{x}_{\mathrm{im}})\hfill & \text{beam B }p.\hfill \end{array}$$
(C33)
The best-fit solution for $`𝐪_{n+1}`$ and $`𝐮_{n+1}`$ is then

$$\left(\begin{array}{c}𝐪_{n+1}(p)\\ 𝐮_{n+1}(p)\end{array}\right)=\frac{1}{\mathrm{\Delta }}\left(\begin{array}{cc}_iw_is_i^2& _iw_ic_is_i\\ _iw_ic_is_i& _iw_ic_i^2\end{array}\right)\left(\begin{array}{c}_iw_ic_i𝐩^{}(t_i)\\ _iw_is_i𝐩^{}(t_i)\end{array}\right),$$
(C34)
where the weight $`w_i`$ is the same as equation (C31).

### C.4 The Choice of Reference Direction

We adopt the sign convention of Zaldarriaga & Seljak (1997) in which the polarization components are defined in a right-handed coordinate system with the z-axis pointed outward towards the sky. The Stokes parameters are defined with respect to a fiducial direction in each pixel on the sky. We adopt the convention in which the reference direction is aligned with the local Galactic meridian, i.e., the great circle connecting a given point to the North Galactic Pole. The unit vector, $`\widehat{n}`$, tangent to this great circle, pointing in the direction of the North Pole, is given by

$$\widehat{n}=\widehat{l}\times \widehat{e}=\widehat{l}\times \frac{(\widehat{z}\times \widehat{l})}{\mathrm{sin}\theta }$$

where $`\widehat{l}`$ is the unit vector along the line of sight of the current observation, $`\widehat{e}`$ is a unit vector pointing east from $`\widehat{l}`$, $`\widehat{z}`$ is the unit vector to the North Galactic Pole, and $`\theta `$ is the polar angle (co-latitude) of $`\widehat{l}`$.

For reference, we give formulae for computing the factors $`\mathrm{cos}2\gamma `$ and $`\mathrm{sin}2\gamma `$ here. Let $`\widehat{l}`$ be the unit vector along the line of sight, $`\widehat{w}`$ be the unit vector pointing west from $`\widehat{l}`$, $`\widehat{n}`$ be the unit vector pointing north from $`\widehat{l}`$ (the polarization reference direction) and $`\widehat{p}`$ be the unit vector along the polarization plane defined by the axial port of the OMT. Then, for both the $`A`$ and $`B`$ sides, we have

$`\mathrm{cos}\gamma `$ $`=`$ $`\widehat{p}\widehat{n}`$
$`\mathrm{sin}\gamma `$ $`=`$ $`\widehat{p}\widehat{w}`$
$`\mathrm{cos}2\gamma `$ $`=`$ $`2\mathrm{cos}^2\gamma 1`$
$`\mathrm{sin}2\gamma `$ $`=`$ $`2\mathrm{sin}\gamma \mathrm{cos}\gamma .`$
Note that this defines a right-handed coordinate system with $`(x,y,z)`$ axes $`(\widehat{n},\widehat{w},\widehat{l})`$ whose $`z`$ axis is oriented outward following the sign conventions of Zaldarriaga & Seljak (1997).