NINE-YEAR WILKINSON MICROWAVE ANISOTROPY PROBE (WMAP) OBSERVATIONS: FINAL MAPS AND RESULTS

The full nine-year WMAP archive of nominal survey data covers 00:00:00 UT 2001 August 10 (day number 222) to 00:00:00 UT 2010 August 10 (day number 222). Individual year demarcations begin at 00:00:00 UT on day number 222 of a year and end at 23:59:59 UT on day 221 of the following year. In addition to processing improvements, the WMAP nine-year release includes new data accumulated during mission years 8 and 9. Flight operations during those final two years included five scheduled station-keeping maneuvers, a lunar shadow passage, and special commanding procedures invoked within the last mission year to accommodate a compromised battery and transmitter. Overall, WMAP achieved a total mission observing efficiency of roughly 98.4%. The bulk of data excluded from science analysis use are dominated by time intervals that do not exhibit sufficient thermal stability.

The WMAP solar arrays were exposed to constant sunlight so the battery was trickle charged for almost a decade. This activated an internal battery design imperfection and caused battery voltage fluctuations in the final months of the mission (Greason et al. 2012). The resulting thermal variations were beyond what had been experienced earlier in the mission. A detailed analysis of time-ordered data (TOD) with sky signal subtracted showed no detectable dependence on thermal variations associated with battery events, and thus preservation of data was preferred to excision. Out of an abundance of caution, time sequences that contained some of the more egregious temperature excursions were flagged as suspect and omitted from use in the nine-year data processing even though there was no specific evidence of adverse effects.

For each observation, sky pointings of individual WMAP feed horns are computed using boresight vectors in spacecraft body coordinates coupled with the spacecraft attitude solution provided by on-board star trackers. After the first mission year, it was discovered that the apparent attitude computed by the trackers includes small errors induced by thermal flexure of the tracker mounting structure, as described by Jarosik et al. (2007). The amplitude of the flexure is time-dependent and driven by spacecraft temperature gradients. The spacecraft temperature responds both to solar heating and internal power dissipation, and is monitored by thermistors mounted at different locations on the spacecraft (Greason et al. 2012).

Telemetered spacecraft quaternions from the star trackers are corrected for this thermal effect at the very beginning of ground processing, when the raw science archive is created. Originally, we adopted a simple linear model, assuming a fixed angular rate of elevation change in units of arcsec per unit temperature change. As the mission progressed and additional data was used to improve the accumulated thermal profile history, the model has evolved to include angular corrections both in elevation (the dominant term) and azimuth. The nine-year quaternion correction model updates the rate coefficients in both azimuth and elevation, and uses readings from two separate thermistors to characterize the spacecraft temperature gradients. A more detailed description is provided by Greason et al. (2012). The residual pointing error after applying of the correction algorithm is computed using observations of Jupiter and Saturn. The upper limit of the estimated error is 10''.

Beam boresight vectors have been updated based on the full nine-year archive. The largest difference between the seven-year and nine-year line-of-sight (LOS) vectors is 3''. Both the calibrated and uncalibrated WMAP archive data products include documentation of these LOS vectors.

Calibration of TOD from each WMAP radiometer channel requires the derivation of time-dependent gains (responsivity, in units of counts mK⁻¹) and baselines (in units of counts) that are used to convert raw differential data into temperature units. Algorithmic details and underlying concepts are set forth in Hinshaw et al. (2007). Jarosik et al. (2011) outline the calibration process as consisting of two general steps. The first step determines baselines and preliminary gains on an hourly or daily basis via an iterative process that combines a skymap estimation with a calibration solution that updates with each iteration. Baselines and gains are computed by fitting sky-subtracted TOD to the dipole anisotropy induced by the motion of the WMAP spacecraft with respect to the CMB rest frame. The second calibration step determines absolute gain and fits a parameterized gain model to the dipole gains derived in the first step.

The form of the parameterized gain model is based on a physical understanding of radiometer performance, and uses telemetered measures of instrument temperatures and the radio frequency (RF) biases. The model provides a smooth characterization of the responsivity with time and allows higher time resolution than provided by the dipole-fit gains. For the nine-year analysis, we augment the gain model by adding a time-dependent linear trend term, mΔt + c, to the parameterized form presented in Jarosik et al. (2007). Here Δt is an elapsed mission time in days, and m, c are additional fit parameters. Physically, the linear trend can be thought of as a radiometer aging term. Without the addition of this term, model fits to the nine-year dipole gain measurements exhibited small systematic deviations from zero-mean residuals for nine of the 40 WMAP channels. The four Ka1 channels were most affected; the inclusion of the gain model aging term prevents an induced total gain error of about 0.1% in this band. Of the 40 WMAP radiometer channels, W323 alone has shown poor convergence in the iterative procedure that determines dipole-fit gains. Upon investigation we found that this problem is peculiar to the iterative algorithm and not the data itself. The W323 calibration has not been substantially affected in previous releases, but for the nine-year analysis the diverging mode was identified and we disallowed it in the gain model fit.

We continue to conservatively estimate an absolute calibration uncertainty of 0.2% (1σ), based on end-to-end gain recovery simulations. The overall change in calibration for the nine-year processing relative to the seven-year release is −0.031%, +0.048%, −0.005%, +0.041%, and +0.025% for K-, Ka-, Q-, V-, and W-bands respectively; a positive change indicates that features in the nine-year maps are slightly larger than those in the equivalent seven-year maps (i.e., a slight decrease in nine-year absolute gain compared to seven-year).

The transmission efficiencies of sky signals through the A-side and B-side optical systems into each WMAP radiometer differ slightly from one another. This deviation from ideal behavior is characterized in map-making and data analysis through the use of time-independent transmission imbalance factors. The method by which these factors are determined from the WMAP data was described by Jarosik et al. (2007). The determination improves with additional data. These factors have been updated for the nine-year analysis and are presented in Table 1. The nine-year values compare well against the previously published seven-year values (Jarosik et al. 2011) within the quoted uncertainties.

The standard WMAP map-making procedure is unchanged from the previous release and the resulting maps are used for the core cosmological analyses. Progress has been made on the algorithm for estimating the noise properties of the maps. The Stokes I noise levels (σ₀) are now more self-consistent between maps at angular resolution r9 and r10¹⁹ than they had been previously. Another difference from previous analyses is that this procedure now determines the noise in the polarized maps from the Stokes Q and U year-to-year differences while including a spurious ("S") map term, and a mean monopole is subtracted from each S map, as is done separately for Stokes I in the temperature map analysis. A detailed discussion is in Section 4.1.

Data are masked in the map-making process when one feed observes bright foregrounds (e.g., in the Galactic plane) while the corresponding differencing feed observes a far fainter sky. This masking prevents the contamination of faint pixels. Previous WMAP data analysis efforts used a single processing mask, based on the K-band temperature maps, to define which pixel-pairs to mask for all of the frequency bands. In the current processing we have changed to masking based on the brightness in each individual band.

The standard maps are used to subtract the background from Jupiter observations to create beam maps, as has been done in previous processing. We correct three seasons of Jupiter maps in the latter part of the mission for the proximity of Uranus and Neptune to Jupiter. Two-dimensional profiles from the newly updated beam map data are now also used as inputs for the new beam-symmetrized map-making procedure, described below.

In addition to the standard map-making, a new map-making procedure, described in Section 4.2, effectively deconvolves the beam sidelobes to produce maps with the true sky signal convolved by symmetrized beams. As a result of this new procedure, the previously reported map power asymmetry, which we speculated was due to the asymmetric beams and not cosmology (Bennett et al. 2011) has indeed been mitigated in the new beam-symmetrized maps.

In this paper we use the beam-symmetrized maps for diffuse foreground analysis (Section 5.3), but not for estimating the angular power spectrum and cosmological parameters. This is because the deconvolution process introduces correlations in the pixel noise on the beam scale and it is impractical to track these correlations at the full pixel resolution. Diffuse foreground analyses, on the other hand, used maps smoothed to a 1° scale. Appendix B of Hinshaw et al. (2007) demonstrated that the cosmological power spectrum, C_l, is insensitive to beam asymmetry at WMAP's sensitivity level. (It is the 4-point bipolar power spectrum, not the 2-point angular power spectrum, that is sensitive to beam asymmetry.) Use of the beam-symmetrized maps for high-l angular power spectrum estimation would invoke the need for high resolution noise covariance matrices, along with far greater computational and storage demands than are now feasible. Given that dense r9 noise covariance matrices are computationally undesirable and the cosmological power spectrum is insensitive to beam asymmetry, we do not use beam-symmetrized maps for cosmology.

The algorithm used to reconstruct sky maps from differential data masks selected observations to minimize artifacts associated with regions of high foreground intensity. (Jarosik et al. 2011). Observations for which one of the telescope beams is in a region of high foreground intensity gradients while the other is in a low gradient region are only applied to the pixel in the high foreground region as the map solutions are generated. This "asymmetric" masking suppresses map reconstruction artifacts in the low foreground emission regions used for CMB analysis. These artifacts arise from small variations in the power sampled by the telescope beams for different observations that fall within the same map pixel. The variations result from a combination of the finite pixel size and beam ellipticity that both couple to spatial intensity gradients. A processing mask is used to delineate the regions of high foreground intensity gradients. Previous data releases used a common processing mask for all frequency bands based on the K-band temperature maps, even though the foreground intensities vary greatly by band. The current release uses different masks for each frequency band and therefore utilizes the data more efficiently.

Masks for each frequency band are generated using an algorithm that estimates the magnitude of processing artifacts in each r4 pixel given the WMAP scan pattern, a candidate processing mask and the seven-year map of the sky temperature in that band. The magnitude of artifacts, ξ, in a resolution r4 pixel, p₄, is modeled as proportional to the mean magnitude of the temperature gradients within all the reference pixels used in the observations contributing to the original pixel,

$$\begin{eqnarray} \ \ \xi (p_4, n) &\simeq & \frac{\alpha }{N_{{\rm tot}}(p_4, n)}\left[ \sum _{ p_A(i)=p_4}w_n(p_B(i))|\nabla T(p_B(i))| \right.\nonumber\\ &&\left.+\, \sum _{ p_B(i)=p_4}w_n(p_A(i))|\nabla T(p_A(i))| \right], \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &&N_{{\rm tot}}(p_4, n) = \sum _{ p_A(i)=p_4}w_n(p_B(i)) + \sum _{ p_B(i)=p_4}w_n(p_A(i)). \qquad\ \ \end{eqnarray} \tag{ 2 }$$

Here p_A(i) and p_B(i) are the r4 pixel indices for the A- and B-side beams for TOD observation i, w_n represents a candidate processing mask with n pixels masked, and the sums are over observations for which the A-side beam and B-side beam point to pixel p₄. The proportionality constant α was evaluated as the amplitude of the response for each telescope beam as it was rotated about its axis while viewing a uniform temperature gradient, yielding values from 0032 to 0087 for the different beams. The magnitude of the temperature gradient in each r4 pixel is approximated as the standard deviation of the r9 pixels comprising each r4 pixel

$$\begin{eqnarray} \qquad |\nabla T(p_4)| &\simeq & \beta \cdot [{\rm var}(p_9 \in p_4) - \sigma ^2(p_4)]^{1/2}, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \sigma ^2(p_4) &=& \frac{ \sum _{p_9 \in p_4} \sigma _0^2\left/\right. N_{{\rm obs}}(p_9)}{ \sum _{p_9 \in p_4}1},\quad \end{eqnarray} \tag{ 4 }$$

where the last term in Equation (3) removes the bias introduced by the radiometer noise, σ₀ is the noise for one observation and N_obs(p₉) is the number of observations in r9 pixel p₉. The constant β ≃ 1.1 deg⁻¹ for r4 pixels.

Figure 1 shows a map of ξ(p₄, 0) for the Ka1 DA with no pixels masked in the candidate processing mask (n = 0). The highest value areas in this map correspond to regions that are ≈140° from the Galactic center corresponding to the spacing between the WMAP A-side and B-side telescope beams. Processing masks for each frequency band are generated iteratively starting from an empty mask, n = 0. The r4 pixel added to the candidate mask w_n at each step is that which produces the greatest reduction in the mean value of ξ(p₄, n) for the current value of n. The value of ξ is then recalculated with the updated candidate mask, w_{n + 1}, and the process repeated. Figure 2 displays how the maximum and mean value of ξ(p₄, n) vary as pixels are added to the mask. The mean and maximum values decrease rapidly as n increases and asymptotes to an approximately constant value for large n. The maximum value of ξ in the asymptotic region is calculated as

$$\begin{eqnarray} \xi ^{{\rm max}}(n) &=& \max _{p_4} \xi (p_4, n),\qquad\quad \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \qquad \xi ^{{\rm max}}_{{\rm sat}} &=& \overline{ \xi ^{{\rm max}}(n)},\ 180 \ge n \ge 360. \end{eqnarray} \tag{ 6 }$$

These steps are executed for each DA and masks for the Ka-, Q-, V-, and W-band DAs are selected by choosing the smallest value of n for which $\xi ^{{\rm max}}(n) < 1.15 \xi ^{{\rm max}}_{{\rm sat}}$. This criterion was selected by requiring that $\overline{\xi } < 5\,\mu {\rm K}$ for Ka-, Q-, V-, and W-bands and that the resulting Q-band mask have approximately the same number of excluded pixels as the mask used in earlier data releases. The mask created in this manner for the Ka1 DA is the final processing mask. Masks for frequency bands with multiple DAs are formed by merging the individual DA masks such that if a pixel was masked in either of the DA masks it is masked in the combined mask. The K-band processing mask requires special treatment due to the brightness of the foregrounds. Applying the criterion above yields a very large sky mask that leaves many pixels with few or no observations causing convergence problems in the conjugate gradient map solution. The adopted K-band processing mask is the largest w_n formed with K-band inputs for which the sky map solution converges for all years except year 2. Year 2 is particularly problematic due to the location of Jupiter. Achieving convergence requires selection of the w₂₇₀ mask and reduction of the Jupiter exclusion radius to 25. Even with these special considerations the size of the processing mask is still substantially larger than used in previous data releases and should result in reduced artifacts. Table 4 summarizes the mask sizes and planet exclusion radii for the nine-year maps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Map Generation Masking Parameters

	Masked Pixels	$\overline{\xi }$	Planet Exclusion Radii (in °)
Band	(of 3072 Total)	(μK)	Mars	Jupiter	Saturn	Uranus	Neptune
K (yr ≠ 2)	312	7.12	2.0	3.0	2.0	2.0	2.0
K (yr = 2)	270	7.59	2.0	2.5	2.0	2.0	2.0
Ka	212	4.46	1.5	2.5	1.5	1.5	1.5
Q	201	4.31	1.5	2.5	1.5	1.5	1.5
V	125	3.78	1.5	2.2	1.5	1.5	1.5
W	98	3.66	1.5	2.0	1.5	1.5	1.5

Download table as:  ASCIITypeset image

The TOD, d, for a differential radiometer sensitive only to a Stokes I signal may be written as

$$\begin{equation} {\bf d} = {\bf Mt} + {\bf n}. \end{equation} \tag{ 7 }$$

Here M is a sparse N_t × N_p matrix that contains information about the scan pattern and transforms the input sky signal array, t, into a sequence of differential observations, d. The number of time-ordered observations is given by N_t, the number of pixels in array t is N_p, and n is an N_t element array representing the radiometer noise. For the standard map processing each row of M contains two non-zero elements representing the weights given to the input map pixels nearest the A- and B-side telescope LOSs. The first step in generating a sky map is evaluation of the "iteration zero" map,

$$\begin{equation} {\bf \tilde{t}_0} = {\bf M}^{\rm T}_{\rm am} {\bf N}^{-1}{\bf d}, \end{equation} \tag{ 8 }$$

where ${\bf M_{\rm am}^{\rm T}}$ is the transpose of a masked version of the observation matrix, and N⁻¹ is the inverse of the radiometer noise covariance matrix,

$$\begin{equation} {\bf N}^{-1} = \langle {\bf nn}^{\rm T}\rangle ^{-1}, \end{equation} \tag{ 9 }$$

with the angle brackets representing an average. The masking contained in ${\bf M}_{\rm am}^{\rm T}$ prevents contamination of regions of the map with low foreground emission that can occur when one of the telescope beams is in a region of high foreground emission. (See Section 4.1.1.) The reconstructed sky map, ${\bf \tilde{t}}$, is then calculated by solving

$$\begin{equation} {\bf \tilde{t}} = \left({\bf M}_{\rm am}^{\rm T} {\bf N}^{-1}{\bf M}\right)^{-1} {\bf \tilde{t}_0}. \end{equation} \tag{ 10 }$$

The form of matrix M described above ignores the effects of the finite WMAP beam sizes since each observation is modeled using only the value of the input sky signal nearest the LOS direction. The actual radiometric data is an average of the input sky signal spatially weighted by the beam response. Each row of M should therefore contain additional non-zero elements describing the signal contribution from the off-axis beam response. If the beam response was the same for the A- and B-side beams and azimuthally symmetric about the LOS, the observation matrix including the off-axis signal contributions, M_s, could be written in the form

$$\begin{equation} {\bf M_{\rm s}} = {\bf M C}, \end{equation} \tag{ 11 }$$

where C is an N_p × N_p element matrix that performs a convolution by the symmetric beam pattern. Substituting M_s for M in Equation (7) shows that in this limit the sky map reconstructed using Equation (10) is the input map convolved by the symmetric beam pattern, ${\bf \tilde{t}}_{\rm c} = {\bf C t}$.

Following the approach discussed above, we present the nine-year temperature (Stokes I) full sky maps in Figure 3. The corresponding Stokes Q and Stokes U full sky maps are shown in Figures 4 and 5, respectively. Figure 6 shows the nine-year polarized intensity maps of P = (Q² + U²)^0.5 with superposed polarization angle line segments where the signal-to-noise ratio exceeds unity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The noise in the r9 and r10 maps is described assuming the radiometer noise distribution is stationary, has a white spectrum and is normally distributed. With these assumptions it can be shown that the noise component of a Stokes I sky map, t_n, is given by (Jarosik et al. 2011)

$$\begin{equation} \tilde{\bf t}_{\rm n} = ({\bf M^TM})^{-1}\cdot {\bf M^Tn}, \end{equation} \tag{ 12 }$$

where M is the mapping matrix as described in Section 4.1.2 and n is a vector of normally distributed random numbers that characterizes the radiometer noise,

$$\begin{equation} \langle {\bf n} \rangle = 0,\quad \langle {\bf n n^T} \rangle = \sigma _0^2 {\bf I}, \end{equation} \tag{ 13 }$$

where the brackets indicate an ensemble average and σ₀ describes the noise amplitude. The pixel-pixel noise correlation matrix is then

$$\begin{equation} {\bf \Sigma } = \frac{\left\langle {\bf t}_{\rm n} {\bf t}_{\rm n}^{T} \right\rangle }{\sigma _0^2} = ({\bf M^T M})^{-1}. \end{equation} \tag{ 14 }$$

Ideally the value of σ₀ is obtained by evaluating

$$\begin{equation} \sigma _0^2 N_{{\rm pix}}= \left\langle {\bf t}_n^T {\bf \Sigma }^{-1} {\bf t}_n \right\rangle , \end{equation} \tag{ 15 }$$

where N_pix is the number of map pixels, but such a calculation is intractable with high resolution maps. In practice only the diagonal elements of Equation (15) are evaluated. Since

$$\begin{equation} {\bf \Sigma }^{-1} = {\bf M^T M} \end{equation} \tag{ 16 }$$

the diagonal elements of Σ⁻¹ are simply the number of observations²¹ of each pixel, N_obs. Each data sample from a WMAP differential radiometer is a measure of the temperature difference between the sky locations at the A- and B-side telescope boresights. Reconstructing a map from differential data involves two different pixels for each observation, a pixel that is being updated and a reference pixel. The noise in each pixel therefore has contributions both from the noise in the radiometric data for each sample and noise in the value of the reference pixel. If σ₀ represents the radiometer noise for an individual sample, the noise contribution from the reference pixel is approximately $\sigma _0/\sqrt{N_{{\rm obs}}(p)}$, where N_obs(p) is the number of observations used to calculate the value of the reference pixel, p. As the map resolution increases the mean value of N_obs decreases, making the reference pixel noise more significant relative to the radiometer noise. The omission of the off-diagonal terms in the evaluation of Equation (15) ignores the contribution to the noise from the reference beam pixels in the evaluation of σ₀. This effect is evident when the σ₀ values for r9 and r10 versions of the Stokes I sky maps are compared. The σ₀ values from the r10 maps have values from 0.3% (W-band) to 1.5% (K-band) higher than those obtained form the corresponding r9 sky maps. The low sampling rate of the K-band radiometer results in lower N_obs values and hence the largest effect.

A more accurate determination of σ₀ can be made by equating the diagonal elements of Equation (14) since these quantities are directly measurable from the sky maps. The diagonal elements of Σ may be calculated relatively simply given two well justified assumptions: (1) The sky map noise is uncorrelated between pixels; and (2) The reference pixels associated with each main pixel are distinct. With these assumptions diagonal elements of Σ are estimated as

$$\begin{eqnarray} {\bf \Sigma }_{y,y} &=& \left[ \sum _{i, p_A(i) = y}w(p_B(i))\frac{(1 + x_{\rm im})^2}{1 + 1/N_{\rm obs}(p_B(i))}\right.\nonumber\\ &&\left. +\, \sum _{i, p_B(i) = y}w(p_A(i))\frac{(1 - x_{\rm im})^2}{1 + 1/N_{\rm obs}(p_A(i))} \right]^{-1}, \qquad \end{eqnarray} \tag{ 17 }$$

where i is a sample index of the TOD and the sums are over observations for which the A-side and B-side beams observe pixel y. The processing mask is represented by w, which has value zero in masked pixels and unity elsewhere. The 1 ± x_im factors are corrections arising from the transmission imbalance factors and N_obs represents the number of observations contained in the reference beam pixel of the sky map. The 1/N_obs terms in the denominators increase the value of Σ_{y, y} to account for the additional noise arising from the reference beam pixels. In the limit where N_obs is very large for all observations the value of Σ_{y, y} is simply $1/N_{{\rm obs}}(y) = 1/{\bf \Sigma }^{-1}_{y,y}$. The values of σ₀ obtained from r9 and r10 Stokes I maps, evaluated using diagonal elements of Equation (14), agree to ≈0.05% with the worst discrepancy being ≈0.1%. This is a significant improvement relative to the former method.

The N_obs fields of the nine-year r9 and r10 sky maps now contain the reciprocals of the diagonal element of the Σ matrix as it is now considered a more accurate measure of the pixel noise. This change allows the map noise in each pixel to still be calculated as $N = \sigma _0/\sqrt{N_{{\rm obs}}}$ providing that the values of σ₀ published with the nine-year analysis are used. Because the σ₀ values computed from r10 maps differ by less than 0.1% from those computed from r9 maps, the r9 values are adopted for all WMAP nine-year analysis.

These methods have been extended and applied to the Stokes Q and U maps and the spurious response map S. This change improved the agreement between the σ₀ values for the temperatures and polarization maps to ≈0.5% from ≈1.1% in earlier data releases. Table 5 gives the nine-year σ₀ values by DA for temperature (Stokes I) and polarization (Stokes Q, Stokes U), and spurious response S.

The beam-symmetrized maps are generated by solving Equation (21) iteratively using a stabilized bi-conjugate gradient method (Barrett et al. 1994). In this procedure the product

$$\begin{equation} {\bf M}_{\rm am}^{\rm T} {\bf N}^{-1} ({\bf M } + {\bf M}_{\rm n} {\bf C}^{-1})\cdot {\bf \tilde{t}}_{{\rm c},i} \end{equation} \tag{ 22 }$$

is evaluated repeatedly and the solution ${\bf \tilde{t}}_{{\rm c}, i}$ updated after each iteration, i, driving the value of this expression to ${\bf \tilde{t}}_0$. The product (22) is evaluated by looping through the TOD; each observation corresponds to multiplying one row of M + M_nC⁻¹ by the current iteration of the solution, ${\bf \tilde{t}}_{{\rm c},i}$. The first term in each multiplication, ${\bf M} {\bf \tilde{t}}_{{\rm c},i}$, is the weighted sum of the map pixels values nearest the LOS directions of the two beams, corresponding to the differential sky signal smoothed by the axisymmetric beam response. Each row of the matrix M contains two non-zero elements with values (1 + x_im) and (− 1 + x_im), the weight factors for the A- and B-side beams. (The x_im term (|x_im| ≪ 1) accounts for a small imbalance in radiometer response to beam filling signals from the A and B sides.)

The second term in the product of Equation (22), ${\bf M}_{\rm n} {\bf C}^{-1}{\bf \tilde{t}}_{{\rm c},i}$, corresponds to the differential signal from the non-axisymmetric beam response for the current LOS and azimuthal beam orientations. The nonzero elements in each row of M_n are the pixel weights of the non-axisymmetric beam response of the two beams, also weighted by the (± 1 + x_im) factors. To keep the computation time tractable only contributions within a radius r_sl (30 mrad for K-, Ka-band, 26 mrad for Q-, V-, and W-band) of the LOS of each beam are used. The circular regions contributing to the signal for the A and B beams do not overlap, so their contributions may be calculated separately then summed.

The matrix C⁻¹ performs a deconvolution by the symmetrized beam pattern. It is therefore rotationally symmetric and the product M_nC⁻¹ may be evaluated once for each beam, forming convolution kernels K_A and K_B. The contribution of ${\bf M}_{\rm n} {\bf C}^{-1}{\bf \tilde{t}}_{{\rm c},i}$ for each beam is then evaluated by mapping these kernels to the corresponding pixels of ${\bf \tilde{t}}_{{\rm c},i}$ for the LOS and azimuthal orientation for each observation and summing their products.

Figure 8 illustrates the steps used in forming the kernel for the Q1 A-side beam. First (in panel a) a map of the non-axisymmetric beam response, M_n, is formed on a Cartesian grid by subtracting the best-fit symmetrized beam profile from the total beam profile in Equation (19). Next the product M_nC⁻¹, is evaluated by performing a Wiener deconvolution of M_n. A Wiener deconvolution is used to minimize the impact of noise on the deconvolved map. (In performing the Wiener weighting the signal component of the result was assumed to be proportional to the input, M_n, while the noise was assumed to be white and its magnitude obtained from portions of the beam map far from the LOS direction.) Even using the Wiener weighting, some noise remains in the deconvolved maps at relatively large radii from the LOS direction. A cosine apodization function is therefore introduced to smoothly taper the value of the kernel to zero at radial distance r_sl from the beam LOS. This procedure eliminates artifacts in the maps that would be caused by a sharp cutoff of the kernel noise at the radius r_sl. The fidelity of the kernel is demonstrated in Figures 8(e) and (f) that show the kernel re-convolved with the symmetrized beam. After re-convolution the majority of the non-axisymmetric beam response is recovered without the introduction of excessive noise.

Zoom In Zoom Out Reset image size

Ideally the kernel weights representing the non-axisymmetric beam response sum to zero for each observation. This is only approximately true in practice since the HEALPix pixelization used for the solution ${\bf \tilde{t}}_{{\rm c},1}$ and the Cartesian grid of the kernel are incommensurate, resulting in slightly different combinations of weights being used for different LOS directions and azimuthal beam orientations. This results in small variations of the total weight for observation of different points on the sky.

The mean value of a map generated by ideal differential data is unconstrained. The non-idealities in the radiometers parameterized by the transmission imbalance factors, x_im, weakly constrain the mean value of the maps, but occasionally maps solutions with relative large mean values are generated. The spatially varying total weights described above can couple to these mean values resulting in small spurious map features. This problem is remedied by subtracting the sum of the kernel weights used for each observation from the value in M corresponding to the weight of the LOS pixel, resulting in a uniform weight for each observation. This choice insures that the total weight of the A- and B-side observations are (1 + x_im) and (− 1 + x_im) respectively, guaranteeing that the beam-symmetrized maps agree with the normal maps at angular scales larger than the characteristic size of the convolution kernels.

Figure 9 displays the ratio of the TT power spectra of the beam-symmetrized maps to those of the normally processed maps and ratios as predicted in Hinshaw et al. (2007). The spectra from the different map processings agree exactly at low l as expected and agree with the predictions within 2% in regions of adequate signal-to-noise ratios.

Zoom In Zoom Out Reset image size

A detailed analysis of WMAP seven-year observations of planets and selected celestial calibrators is given by Weiland et al. (2011), including intercomparisons with relevant results in the literature. Here we concentrate on updated nine-year WMAP results for some of these sources.

5.2.1.1. Jupiter. Mean nine-year Jupiter temperatures are derived from the l = 0 component of the unnormalized beam transfer functions B_l. The symmetrized beam response to Jupiter, T_pkΩ_beam, may be directly derived from B₀. As described in Weiland et al. (2011), all Jupiter observations have been corrected to a fiducial solid angle $\Omega _{{\rm Jup}}^{{\rm ref}} = 2.481\times 10^{-8}\,{\rm sr}$. Mean Jupiter temperatures T_Jup are thus computed using the relation $T_{{\rm Jup}} = T_{pk}\Omega _{{\rm beam}}/\Omega _{{\rm Jup}}^{{\rm ref}}$. These temperatures are presented in Table 6. Quoted uncertainties are a quadrature sum of estimated beam solid angle errors from Table 3 and the uncertainty in the absolute calibration. The mean Jupiter temperatures derived from the five-year, seven-year, and nine-year data releases are consistent with each other within the quoted uncertainties.

Table 6. Nine-year Mean Jupiter Temperatures

Band	$\nu _e^\mathrm{RJ}$ a	λ b	Tc	σ(T)d
	(GHz)	(mm)	(K)	(K)
Per DA
K1	22.82	13.1	136.1	0.75
Ka1	33.07	9.1	147.1	0.68
Q1	40.88	7.3	153.9	0.78
Q2	40.67	7.4	154.7	0.76
V1	60.37	5.0	164.9	0.71
V2	61.24	4.9	165.9	0.68
W1	93.25	3.2	172.5	0.84
W2	93.73	3.2	173.4	0.85
W3	92.72	3.2	173.1	0.87
W4	93.57	3.2	172.3	0.86
Per band
K	22.82	13.1	136.1	0.75
Ka	33.07	9.1	147.1	0.68
Q	40.78	7.3	154.3	0.59
V	60.81	4.9	165.4	0.54
W	93.32	3.2	172.8	0.52

Notes. ^aNine-year values; see Appendix A. ^b$\lambda = c/\nu _e^{\mathrm{RJ}}$. ^cBrightness temperature calculated for a solid angle Ω_ref = 2.481 × 10⁻⁸ sr at a fiducial distance of 5.2 AU. Temperature is with respect to blank sky: absolute brightness temperature is obtained by adding 2.2, 2.0, 1.9, 1.5, and 1.1 K in bands K, Ka, Q, V, and W respectively (Page et al. 2003a). Jupiter temperatures are uncorrected for a small synchrotron emission component (see Weiland et al. 2011). ^dComputed from errors in Ω_B (Table 3) summed in quadrature with absolute calibration error of 0.2%.

Download table as:  ASCIITypeset image

The stability of Jupiter emission over the nine-year baseline is evaluated by computing seasonal temperatures per DA and comparing them to their nine-year means. We compute ΔT/T as the mean deviation of all DAs from their nine-year mean values, and include a 1σ standard deviation as a measure of coherency. These results are listed in Table 7. From the seven-year analysis, Weiland et al. (2011) placed an upper limit on variability of 0.2% ± 0.4%. Although consistent with this value, the Jupiter observations from the last two seasons introduce the statistically weak (probability to exceed, PTE = 14%) possibility of a decreasing trend in temperature with time. Given our measurement uncertainties, a constant temperature is a very good fit to the data and that is what we use in our analysis.

Out of caution, we examined the hypothesis that there might be instrumental or calibration issues contributing to slightly lower Jupiter temperatures computed for the last few seasons of data. To determine if there might be a systematic calibration error within the last two years of the mission, yearly flux values for celestial sources Cas A, Cyg A, and Tau A were computed and compared against seven-year trends; no evidence for any calibration inconsistency was found. Since Jupiter is not a steep-spectrum source, bandpass center frequency variations are also not an important factor; we expect an effect of less than ±0.05% over the 9 years in the K- through V-bands. In terms of Jupiter itself, there is no clear temperature trend with Sun-Jupiter distance or sub-WMAP latitude.

5.2.1.2. Saturn. As seen by the WMAP satellite, the spatially unresolved microwave brightness of Saturn varies with orbital phase as the projected area of the ring system and oblate planetary spheroid changes aspect. Weiland et al. (2011) developed an empirical, geometrically motivated model to predict Saturn's apparent brightness at WMAP frequencies, based on the first seven years (14 seasons) of observations. The available range of observable ring opening angles during this seven year interval falls in the range −28° < B < −6°. Weiland et al. (2011) found that parameter covariance and potential systematics in their model fit permitted a determination of Saturn's disk temperature to within roughly 3–4 K, but noted that the inclusion of lower inclination observations in the fit should decrease the uncertainty in the derived model parameters. WMAP observations from the last two mission years include four new Saturn observing seasons, numbered 15 through 18. Since the Saturn ring system presented an "edge-on" configuration in early 2009, these four new seasons span the cross-over from viewing the rings from below (negative B) to viewing them from above (positive B) as seen in Table 8. These new observations at low B provide the opportunity to better constrain the predictive model for WMAP frequencies.

We apply the analysis methods of Weiland et al. (2011) to the nine-year compendium of Saturn observations to derive mean apparent temperatures of the Saturn system per DA per observing season, presented in Table 8. The analysis can be summarized as a three-step process. First, a time-ordered archive of Saturn observations is created, and sky signals arising from the Galaxy and CMB are removed, either through use of sky subtraction or masking. Second, the individual observations from this background subtracted archive are binned to form mean radial Saturn response profiles for each season and DA. Finally, the WMAP beam radial profile per DA (as determined from Jupiter observations) is fit to the Saturn radial response for that DA and an apparent temperature is derived. Temperature entries for the first 14 seasons listed in Table 8 may be directly compared against those in Table 9 of Weiland et al. (2011). There are small differences of order 0.5 to 1 σ between some of entries in common between the seven-year analysis and the nine-year analysis presented here. Differences of this nature are expected and can be traced to small variations in calibration, beam characterization and data masking between the seven-year and nine-year processing.

The temperatures in Table 8 may be fit with an empirical model that predicts Saturn's unresolved microwave brightness T as a function of ring opening angle and frequency. We adopt the same model formulation as in the seven-year analysis of Weiland et al. (2011), which employs a simple geometrical summation of emission from the unobscured planetary disk, emission from the ring system and emission from those portions of the disk obscured by the rings:

$$\begin{eqnarray} T(\nu ,B) &=& T_{\mathrm{disk}}(\nu) \left[A_{\mathrm{ud}} + \sum _{i=1}^7 e^{-\tau _{0,i}|\csc B|} A_{\mathrm{od},i}\right]\nonumber\\ && + T_{\mathrm{ring}}(\nu) \sum _{i=1}^7 A_{r,i}. \end{eqnarray} \tag{ 23 }$$

At a given frequency ν, a single temperature is assumed for the planetary disk, T_disk(ν). The model allows for seven radially concentric ring divisions. All rings are characterized by the same temperature T_ring(ν), but each of the seven ring sectors has its own ring-normal optical depth τ_{0, i}, with 1 ⩽ i ⩽ 7. Each τ_{0, i} is assumed to be both constant within its ring and frequency independent. A_ud, A_{od, i} and A_{r, i} are the projected areas of the unobscured disk, the portion of the disk that is obscured by ring i, and ith ring, respectively. These areas are normalized to the total (obscured+unobscured) disk area. Model fit parameters are the five Saturn uniform disk temperatures and five mean ring temperatures (one for each WMAP frequency). The geometrical ring boundaries and relative ratios τ_{0, i}/τ_{0, max} are constrained as per Table 10 of Weiland et al. (2011), where τ_{0, max} is the ring-normal optical depth for the most optically thick ring (ring 3, i.e., the outer B ring). For the nine-year fit, the value of τ_{0, max} was also allowed to be a fit parameter, although in practice its inclusion makes very little difference in the fit results.

The nine-year model fit returns a reduced χ² of ∼1.04 for ∼150 degrees of freedom; the model fit and residuals per WMAP frequency are shown in Figure 13. On average, the rms of the residuals is ∼1% per frequency; the value for Q-band is somewhat higher (1.3%) and that for V-band is lowest (0.7%). Model parameters and their formal errors σ_fit are presented in Table 9. By construction, the T_disk and T_ring model parameters are anti-correlated. The covariance between these parameters allows the possibility of systematic errors not accounted for in the fitting formalism. Although the mean disk temperature is reasonably well constrained by the new WMAP observations from seasons 15–18, hemispheric temperature gradients or local hot spots would negate the assumed symmetry of the empirical model, and would affect the derived mean ring temperatures. The nine-year baseline unfortunately does not extend far enough toward positive B to assess the limits of the symmetry assumption. Additionally, the model's assumed ring optical depth profile may not be accurate. As with the seven-year analysis, we use a model variant to estimate systematic differences between models which return similar values of χ². Our worst case estimate allows for differences of 0.9 K in T_disk and 0.7 K in T_ring; we add these to the formal fitting errors in Table 9 to produce the tabulated adopted error, σ_adopted. The T_disk and T_ring parameters are plotted along with their adopted errors in Figure 14. Within the conservative adopted errors, the nine-year derived disk and ring temperatures are in agreement with those from the seven-year fit; the nine-year adopted errors for T_disk are roughly half those quoted for the seven-year fit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As for the seven-year analysis, two separate methods have been used for the identification of point sources from WMAP maps and two separate point source tables have been produced. Both methods are largely unchanged from the seven-year analysis (Gold et al. 2011). Since the use of beam-symmetrized maps would result in only minor changes to the recovered source fluxes and since there is benefit to continuity with previous WMAP point source analyses, we have generated the source catalogs from maps that are not deconvolved. The first method searches for point sources in each of the five WMAP wavelength bands. The nine-year signal-to-noise ratio map in each band is filtered in harmonic space by $b_l/(b^2_lC^{\rm cmb}_l + C^{\rm noise}_l)$, where b_l is the transfer function of the WMAP beam response, $C^{\rm cmb}_l$ is the CMB angular power spectrum, and $C^{\rm noise}_l$ is the noise power (Tegmark & de Oliveira-Costa 1998; Refregier et al. 2000). The filtering suppresses CMB and Galactic foreground fluctuations relative to point sources. For each peak in the filtered maps that is >5σ in any band, the unfiltered temperature map in each band is fit with the sum of a planar base level and a beam template formed by convolving an azimuthally symmetrized beam profile with a skymap pixel. (This method was previously used by Weiland et al. (2011) for selected celestial calibration sources and is more accurate than the Gaussian fitting that was used for the seven-year and earlier point source analyses.) The peak temperature from each beam template fit is converted to a source flux density using the conversion factor Γ given in Table 3. The flux density uncertainty is calculated from the 1σ uncertainty in the peak temperature, and does not include any additional uncertainty due to Eddington bias. Uncertainty due to beam asymmetry effects has been found to be negligible, about 0.1% or less, by comparing results from beam template fits to the normally processed K-band map with those to the beam-symmetrized K-band map for Tau A, Cas A, and Cyg A. Flux density values are entered into the catalog for bands where they exceed 2σ and where the source width from an initial Gaussian fit is within a factor of two of the beam width. A point source catalog mask is used to exclude sources in the Galactic plane and Magellanic cloud regions. This mask has changed from the seven-year analysis in accordance with changes in the KQ85 temperature analysis mask. A map pixel is outside of the nine-year point source catalog mask if it is either outside of the diffuse component of the nine-year KQ85 temperature analysis mask or outside of the seven-year point source catalog mask. The new catalog mask admits 83% of the sky.

The second method of point source identification is the CMB-free method originally applied to one-year and three-year V- and W-band maps by Chen & Wright (2008) and to five-year V- and W-band maps by Wright et al. (2009). The method used here is that applied to five-year Q-, V-, and W-band maps by Chen & Wright (2009) and to seven-year Q-, V-, and W-band maps by Gold et al. (2011). The V- and W-band maps are smoothed to Q-band resolution. An ILC map (see Section 5.3.3) is then formed from the three maps using weights such that CMB fluctuations are removed, flat-spectrum point sources are retained with fluxes normalized to Q-band, and the variance of the ILC map is minimized. The ILC map is filtered to reduce noise and suppress large angular scale structure. Peaks in the filtered map that are >5σ and outside of the nine-year point source catalog mask are identified as point sources, and source positions are obtained by fitting the beam profile plus a baseline to the filtered map for each source. For the nine-year analysis, the position of the brightest pixel is adopted instead of the fit position in rare instances where they differ by >01. Source fluxes are estimated by integrating the Q, V, and W temperature maps within 125 of each source position, with a weighting function to enhance the contrast of the point source relative to background fluctuations, and applying a correction for Eddington bias due to noise (sometimes called "deboosting").

We identify possible 5 GHz counterparts to the WMAP sources found by both methods by cross-correlating with the GB6 (Gregory et al. 1996), PMN (Griffith et al. 1994, 1995; Wright et al. 1994, 1996), Kühr et al. (1981), and Healey et al. (2009) catalogs. A 5 GHz source is identified as a counterpart if it lies within 11' of the WMAP source position (the mean WMAP source position uncertainty is 4'). When two or more 5 GHz sources are within 11', the brightest is assumed to be the counterpart and a multiple identification flag is entered in the catalog.

The nine-year five-band point source catalog is presented in Appendix B and the nine-year QVW point source catalog is presented in Appendix C. The five-band catalog contains 501 sources, the QVW catalog contains 502 sources, and the two catalogs have 387 sources in common. As noted by Gold et al. (2011), differences in the source populations detected by the two search methods are largely caused by Eddington bias in the five-band source detections due to CMB fluctuations and noise. At low flux levels, the five-band method tends to detect point sources located on positive CMB fluctuations and to overestimate their fluxes, and it tends to miss sources located in negative CMB fluctuations. Other point source detection methods have been applied to WMAP data and have identified sources not found by our methods (e.g., Scodeller et al. 2012; Lanz 2012; Ramos et al. 2011, and references therein).

In this section we evaluate the diffuse foreground emission both for the purpose of separation from the CMB anisotropy and for characterizing the nature of the foreground components. As a prelude to our cosmological analyses we fit and remove external foreground template map data from the WMAP maps and we mask remaining regions estimated to be significantly contaminated. We discuss this temperature and polarization cleaning, and the masks, below. To elucidate the characteristics and nature of the diffuse foregrounds we implement four techniques: ILC technique; MEM; MCMC fits; and χ² fits.

Our analysis of the diffuse foregrounds generally uses the five bands of WMAP data in conjunction with other data sets. WMAP was designed to observe in the spectral region where the ratio of the CMB to the foregrounds is at its maximum. This minimizes the amplitude of contamination and needed corrections or masking, which is good for cosmology. To achieve an improved signal-to-noise ratio of the foregrounds themselves, it is sometimes useful to use external data where the foreground emission is weak.

Foreground analyses are done using 1° smoothed beam-symmetrized nine-year temperature maps in the five WMAP bands. As in our previous foreground studies, the zero level of each map is set such that a fit to the ILC-subtracted map of the form $T(\vert b \vert) = T_p \, \csc \vert b \vert + c$, over the range −90° < b < −15°, yields c = 0. This assumes a plane-parallel slab model for the Galactic emission. Formal 1σ uncertainties in the map zero levels (calculated as the quadrature sum of (1) the uncertainty in the fit intercept c and (2) the difference in intercepts from southern and northern Galactic hemisphere fits) are 7.2, 5.9, 3.6, 1.8, and 0.76 μK in thermodynamic units for K-, Ka-, Q-, V-, and W-bands respectively. The South Galactic pole brightness T_p from the fitting is 77.9 ± 1.5, 30.1 ± 0.6, 17.7 ± 0.4, 8.6 ± 0.2, and 9.4 ± 0.3 μK in thermodynamic units for K-, Ka-, Q-, V-, and W-bands respectively. The Stokes Q and U maps have well-defined zero levels and no monopole corrections are applied to them.

Previous WMAP team analyses have used the Finkbeiner (2003) Hα map corrected for extinction as a template for free–free emission (Bennett et al. 2003c). The Finkbeiner map is a composite of the Virginia Tech Spectral line Survey (Dennison et al. 1998), the Southern H-Alpha Sky Survey Atlas (Gaustad et al. 2001), and the Wisconsin H-Alpha Mapper survey (Haffner et al. 2003). The extinction correction assumes that Hα emission and extinction are uniformly mixed along each LOS,

$$\begin{equation} I(\rm H\alpha)_{\rm extinction\hbox{-}corrected} = \it I(\rm H\alpha) \: \tau /(1-e^{-\tau }). \end{equation} \tag{ 24 }$$

Here τ is the dust optical depth at the wavelength of Hα and was calculated from the E(B − V) map of Schlegel et al. (1998) as

$$\begin{equation} \tau = 2.2 \thinspace E(B-V), \end{equation} \tag{ 25 }$$

which assumes an extinction law for R_V = 3.1, characteristic of the diffuse interstellar medium.

Recent studies of selected dust clouds at 20° < |b| < 40° have shown that scattered Hα can make a significant contribution to the observed Hα brightness for some LOSs (Mattila et al. 2007; Lehtinen et al. 2010; Witt et al. 2010). Here we apply an approximate correction to our previous Hα template for the contribution of scattered Hα, based on correlations between Hα and 100 μm emission found by Witt et al. (2010) for four selected clouds and by Brandt & Draine (2012) for Sloan Digital Sky Survey blank sky regions at intermediate to high Galactic latitudes. Brandt and Draine noted that I(100 μm) varies in proportion to the product of the dust column density and the radiation field that heats the dust. If the spatial variation of the illuminating Hα radiation field in the Galaxy is similar to that of the radiation responsible for dust heating, I(100 μm) may be a good tracer of scattered Hα. The scattering correction we adopt is

$$\begin{equation} I({\rm H}\alpha)_{\rm scattering\scriptsize {\hbox{-}}corrected} = I({\rm H}\alpha)_{\rm extinction\scriptsize {\hbox{-}}corrected} - 0.11 \, I(100\,\mu {\rm m}), \end{equation} \tag{ 26 }$$

where I(Hα) is in Rayleighs, I(100 μm) is the Schlegel et al. (1998) 100 μm map in MJy sr⁻¹, and the I(100 μm) coefficient is a mean of the values of 0.129 ± 0.015 R/(MJy sr⁻¹) found by Witt et al. (2010) and 0.090 ± 0.017 R/(MJy sr⁻¹) found by Brandt & Draine (2012). These correlation slopes were measured for regions with τ < 1, but we apply Equation (26) over the entire sky. This assumes that the Equation (24) extinction correction is valid for the scattered component (i.e., the scattered Hα emissivity and the dust extinction are uniformly mixed along each LOS) and it neglects effects of multiple scattering that may be important for LOSs with high optical depth. The Hα template is made by applying the corrections for extinction and scattering to version 1.1 of the Finkbeiner Hα map, smoothing from 6' FWHM to 1° FWHM, and setting a small number of negative pixels to zero. The resulting Hα-based microwave template is shown in Figure 15 as the "Free–Free Template."

Zoom In Zoom Out Reset image size

Uncertainties in both the extinction correction and the scattering correction are large for high τ, but we find that results of our analyses using the template are not sensitive to these uncertainties. For the foreground cleaning of the temperature maps, the mask used in template fitting is chosen to minimize the combined effects of template error and foreground-CMB covariance (Section 5.3.2). For the MEM foreground fitting, the free–free prior is formed from the Hα template, but for high τ LOSs the observed brightness in the WMAP bands is great enough that the MEM results are not strongly affected by the free–free prior.

Prior to the nine-year analysis, the Haslam map used in the MCMC fitting and as a prior in the MEM fitting was the Fourier-filtered version available from LAMBDA. This version mitigates scan striping in the Haslam map, but also removes many strong point sources. Removal of the point sources affected the quality of some foreground fits for pixels in the Galactic plane. For this reason, the unfiltered Haslam map (also available on LAMBDA) is now used for these applications and its projection to K-band is shown in Figure 15.

All-sky templates of Galactic foregrounds or combinations of foregrounds which are "CMB-free" are fit in a least-squares sense to the WMAP sky maps to construct a foreground model at each frequency. The foreground model is subtracted from the WMAP sky maps to produce reduced foreground, or "cleaned" maps, which are used in turn for power spectrum analysis. The cleaning is applied to sky maps from the standard map-making procedure, not to beam-symmetrized sky maps. Cleaning of temperature and polarization maps is treated independently.

5.3.2.1. Temperature cleaning

A limited set of all-sky foreground templates is available for use in modeling potential contributions from synchrotron, free–free and dust emission. After testing a number of different template combinations, we continue to adopt a foreground model map, M(ν, p), of the form

$$\begin{eqnarray} M(\nu ,p) &=& c_\mathrm{1}(\nu)[T_K(p) - T_{Ka}(p)] + c_\mathrm{2}(\nu)I_\mathrm{H\alpha }(p)\nonumber\\ && + c_\mathrm{3}M_\mathrm{dust}(p), \end{eqnarray} \tag{ 27 }$$

where p indicates the pixel, the frequency dependence is entirely contained in the coefficients c_i, and the spatial templates are the WMAP K−Ka temperature difference map in thermodynamic mK (T_K − T_Ka), an Hα map (I_Hα) in units of Rayleighs, and dust model 8 from Finkbeiner et al. (1999) evaluated at 94 GHz in units of mK antenna temperature (M_dust). The K−Ka template is formed using standard (not beam symmetrized) maps. The values of the coefficients c are such that the model map M(ν, p) is in thermodynamic mK.

However, although the form of the model is the same as that used in previous WMAP analyses, there are modifications in the details of its application. As described in Section 5.3.1, the nine-year extinction corrected Hα template incorporates a scattering correction, a refinement not present in the seven-year analysis. Also, in recognition of the possible contribution of spinning dust to the Galactic emission and the uncertain synchrotron behavior with frequency, spectral and coefficient positivity constraints are no longer imposed in the template fitting. This allows maximum freedom in the fit, but makes physical interpretation of model coefficients more difficult.

There has also been a change in the portion of sky used in computing the foreground model fit. Derived model coefficients are dependent on the fraction of the sky which is fit: a full sky fit minimizes the covariance of the templates with the CMB signature in the WMAP data, but maximizes potential template cleaning residuals (bias) by including sky regions where the templates are more uncertain (generally close to the Galactic plane). For example, the extinction correction applied to the Hα map is only approximate and this template is an imperfect tracer of free–free emission in optically thick regions. In general, as more sky is excluded from the fit, CMB-template covariance increases, while template cleaning bias decreases. The "optimal" sky cut for template fitting may be determined by examining these two competing errors as a function of sky cut, and choosing the mask for which the sum of the two errors is a minimum. For this purpose, several simulated five-band Galaxy models of differing complexity were constructed. Each model is added to a CMB realization, and then cleaned using the algorithm in Equation (27) and a chosen sky cut. This is performed for 100 CMB realizations per sky cut; the mean bias is the template cleaning error and the variance is the CMB covariance. We have used the "KpX" series of Galactic masks, described by Bennett et al. (2003a) as a graduated set of sky cuts. The masking in the "KpX" series is based on K-band intensity: higher values of X indicate a smaller portion of bright sky is cut. For each simulation, the sum of both errors were plotted as a function of sky cut and a rough minimum chosen. Prior to the nine-year analysis, we had used the Kp2 mask for template fitting. However, the simulations indicated a more conservative choice would employ a smaller sky cut. The Kp8 mask was adopted for the nine-year cleaning.

Template cleaning coefficients derived using the updated procedure are shown in Table 10 for the Q, V, and W DAs. As noted previously, the ability of the fit to trade freely among the three templates makes physical interpretation difficult. Monte Carlo simulations have shown that the negative coefficients c₁ derived for W-band result from template covariance with the CMB. The change of template cleaning method from the seven-year to the nine-year analysis has little effect on power spectrum analysis. There is a slight change in the evaluated low-l power spectrum. For 2 ⩽ l ⩽ 16, using the Monte Carlo Apodised Spherical Transform EstimatoR (MASTER) method with the KQ85y9 mask, the absolute value of the change in l(l + 1)/(2π)C_l due to the change in template cleaning is typically 4% of cosmic variance per l.

Table 10. Template Cleaning Temperature Coefficients

DAa	c1b	c2b	c3b
		(μK/R−1)
Q1	0.284	0.890	0.231
Q2	0.284	0.898	0.226
V1	0.0630	0.554	0.686
V2	0.0567	0.541	0.716
W1	−0.0179	0.351	1.609
W2	−0.0182	0.349	1.617
W3	−0.0146	0.342	1.587
W4	−0.0153	0.345	1.594

Notes. ^aWMAP has two differencing assemblies (DAs) for the Q- and V-bands, and four for the W-band; the high signal-to-noise total intensity allows each DA to be fit independently. ^bThe c_i coefficients produce model maps in thermodynamic mK.

Download table as:  ASCIITypeset image

5.3.2.2. Polarization cleaning

The polarization cleaning method is unchanged from the seven-year analysis. The nine-year Stokes Q and U maps are degraded to low resolution (N_side = 16) and the data for pixels outside of the Q-band processing mask are fit to a linear combination of low resolution templates. The fit has the form

$$\begin{equation} [Q(\nu), U(\nu)] = a_1(\nu) \thinspace [Q, U]_K + a_2(\nu) \thinspace [Q, U]_\mathrm{dust}. \end{equation} \tag{ 28 }$$

The template used for synchrotron is the nine-year WMAP K-band polarization, [Q, U]_K. The template for dust, [Q, U]_dust, is constructed from the Finkbeiner et al. (1999, hereafter FDS) dust model 8 evaluated at 94 GHz together with a polarization direction map derived from starlight measurements and a geometric suppression map to account for the magnetic field geometry, as described in Page et al. (2007). The coefficients of the fit to the nine-year data are listed in Table 11 and plotted against frequency in Figure 16.

Zoom In Zoom Out Reset image size

Table 11. Template Cleaning Polarization Coefficients

Band	a1a	βs(νK, ν)b	a2a	βd(ν, νW)b
Ka	0.3204	−3.13	0.0145	1.41
Q	0.1682	−3.13	0.0182	1.50
V	0.0594	−2.97	0.0364	1.41
W	0.0398	−2.43	0.0758	⋅⋅⋅

Notes. ^aThe a_i coefficients are dimensionless and produce model maps in thermodynamic mK. ^bThe spectral indices refer to antenna temperature.

Download table as:  ASCIITypeset image

Full-resolution (N_side = 512) foreground-reduced Stokes Q and U maps were produced using the coefficients in Table 11 with full-resolution versions of the K-band and dust polarization templates smoothed to 1° FWHM. In making the full resolution dust template, the starlight polarization map used to determine polarization direction was upgraded to full resolution using nearest neighbor sampling. Smoothing of the templates to 1° FWHM potentially leaves artifacts in the foreground-reduced maps due to small-scale power or beam asymmetries, but previous analyses have found no sign of these effects (Gold et al. 2011). Data sets for all templates are available on the LAMBDA website.

The spectrum of K-band polarization template coefficients flattens significantly with increasing frequency, which is unexpected for synchrotron emission. This flattening can be understood if, due to shortcomings of the dust polarization template, some fraction of the dust polarization is traced by the K-band template. We illustrate this using a simple model. The polarization maps are modeled as a sum of synchrotron and thermal dust components,

$$\begin{equation} [Q(\nu), U(\nu)] = [Q(\nu), U(\nu)]_\mathrm{synch} + [Q(\nu), U(\nu)]_\mathrm{dust}. \end{equation} \tag{ 29 }$$

Assuming the synchrotron polarization has a power law spectrum and is traced exactly in all bands by the K-band polarization template, the synchrotron component is

$$\begin{equation} [Q(\nu), U(\nu)]_\mathrm{synch} = \frac{g(\nu)}{g(\nu _K)} \left(\frac{\nu }{\nu _K} \right)^{\beta _{{\rm synch}}} [Q, U]_K, \end{equation} \tag{ 30 }$$

where the antenna temperature to thermodynamic temperature conversion factors g are needed because the polarization maps and K-band template are in thermodynamic units. Assuming the dust polarization has a power law spectrum and is traced by a combination of the dust polarization template and the K-band polarization template, with the relative contributions of the two templates independent of frequency, the dust component is

$$\begin{eqnarray} [Q(\nu), U(\nu)]_\mathrm{dust} &=& \frac{g(\nu)}{g(\nu _W)} \!\left(\frac{\nu }{\nu _W} \right)^{\beta _{{\rm dust}}}\nonumber\\ &&\times (f_1 \thinspace [Q,U]_\mathrm{dust} \,{+}\, f_2 \thinspace [Q, U]_K), \end{eqnarray} \tag{ 31 }$$

where f₁ and f₂ are constants. Inserting Equations (30) and (31) in Equation (29) and comparing with Equation (28) gives expressions for the template fit coefficients,

$$\begin{equation} a_1(\nu) = \frac{g(\nu)}{g(\nu _K)} \left(\frac{\nu }{\nu _K} \right)^{\beta _{{\rm synch}}} + f_2 \thinspace \frac{g(\nu)}{g(\nu _W)} \left(\frac{\nu }{\nu _W} \right)^{\beta _{{\rm dust}}} \end{equation} \tag{ 32 }$$

and

$$\begin{equation} a_2(\nu) = f_1 \thinspace \frac{g(\nu)}{g(\nu _W)} \left(\frac{\nu }{\nu _W} \right)^{\beta _{{\rm dust}}}. \end{equation} \tag{ 33 }$$

Fitting these expressions to the a₁(ν) and a₂(ν) values in Table 11 gives β_synch = −3.13, β_dust = 1.44, f₁ = 0.076, and f₂ = 0.024. The fits are shown by the curves in Figure 16. They match the template coefficients very well with no need for an additional emission mechanism such as spinning dust or magnetic dust polarization. In this simple model, the K-band template component contributes about 1/3 of the rms dust polarization and the dust template component contributes about 2/3.

This suggests that there is room for improvement in the dust polarization template. Some alternate dust templates were tested in fitting the polarization maps, but none of them gave significant improvement in χ². These include a template based on K-band polarization directions instead of directions from starlight measurements, a template based on a geometric suppression map calculated from the ratio of observed K-band polarized intensity to K-band synchrotron total intensity from the seven-year MCMC shifted spinning dust model (Gold et al. 2011), and two templates from O'Dea et al. (2012) based on different Galactic magnetic field models.

5.3.2.3. Masks. Sky masks for CMB temperature analysis are generated as described by Gold et al. (2011). The process begins with K- and Q-band maps smoothed to $1\deg$ resolution, from which an estimate of the CMB is subtracted to leave maps that effectively consist of foreground emission alone. The CMB is estimated using the ILC method (Hinshaw et al. 2007). Both the K and the Q maps are masked at a flux contour that leaves either 75% or 85% of the sky unmasked. The K and Q-band sky masks for each cut level are combined so that any pixel excluded by either cut is excluded by the combination. The resulting combinations, dominated by diffuse Galactic emission, are called KQ75 and KQ85, labeled by the admitted sky fraction (f_sky) of the input masks.

These masks are intended primarily to be applied to the foreground-cleaned versions of the sky maps for power spectrum and non-Gaussian analysis. They are made more effective for this purpose by extending them to include regions where the cleaning algorithm is subject to possible systematic error. A χ² analysis is done using differences Q−V and V−W between cleaned band maps at a reduced HEALPix resolution of N_side = 32 (Gorski et al. 2005), or "res 5" in WMAP terminology. Regions of four or more contiguous pixels with χ² higher than four times that of the polar caps are used to define the mask extensions, after 3° smoothing and cleanup steps to remove small "islands" caused by noise.

A point source mask is added to each of the diffuse sky masks. The point source mask from the seven-year analysis is updated to include newly detected WMAP point sources and the 100 GHz sources in the Planck early release compact source catalog. An exclusion radius of 12 is used for sources stronger than 5 Jy in any band and an exclusion radius of 06 is used for weaker sources.

The nine-year versions of the final KQ85 and KQ75 sky masks, called KQ85y9 and KQ75y9, respectively, are compared to the seven-year versions in Figure 17. Changes in the foreground cleaning residuals have altered the morphology of the mask in the vicinity of the Gum Nebula, the Large Magellanic Cloud, and the Galactic center, with the largest relative changes occurring in the KQ85 mask. For KQ75, f_sky is decreased from 70.6% to 68.8%, a difference of 1.8% of the sky, and for KQ85, f_sky is decreased from 78.2% to 74.8%, a difference of 3.7% of the sky.

Zoom In Zoom Out Reset image size

The sky mask for CMB polarization analysis is generated using cuts in K-band polarized intensity and a model of polarized dust emission, together with masking of point sources, as described by Page et al. (2007) and Gold et al. (2009). The nine-year polarization mask is the same as the seven-year version except that three additional point sources were added using a 1° exclusion radius—Hydra A, HB89 1055+018, and BL Lac.

The ILC method seeks to produce a map of the CMB anisotropy from a linear combination of the five WMAP frequency bands. A first application of the method is described by Bennett et al. (2003c). The algorithm was revised slightly by Hinshaw et al. (2007); we refer to this version of the algorithm as the "classic ILC," it has remained unchanged throughout subsequent WMAP data releases. As described in Hinshaw et al. (2007), the algorithm divides the sky into 12 regions—a larger high latitude region and 11 smaller regions spread across the galactic plane. Use of the smaller regions along the plane allows for spatially varying foreground complexity. For each of these smaller regions, five band-weights are computed by minimizing the temperature variance in the region, under the constraint that common-mode CMB signal is unaffected. Weights for the larger high latitude region are computed in a similar manner, but using pixels from locations near the outer-Galactic plane. Exact definitions of these regions are provided on LAMBDA.

We compute the nine-year classic ILC using the coadded deconvolved band maps which have been smoothed to a FWHM of 1°. The weights applied to the 5 frequency maps for each of the 12 sky regions are shown in Table 12. Values for the weights change slightly compared to previous WMAP releases as pixel noise, calibration and beam profiles have been refined.

To the eye, the ILC presents a reasonably foreground-free image of the CMB anisotropy. The beauty of the algorithm is that it is relatively independent of assumptions about foregrounds. However, assessing the underlying uncertainty in the resultant anisotropy map is a difficult problem which heavily relies on knowledge of the Galactic foregrounds. In subsequent sections, we will discuss efforts to improve the classic ILC, as well as characterize the level to which foreground residuals remain.

A MEM-based approach originally developed by Bennett et al. (2003c) and Hinshaw et al. (2007) is used to model the Galactic foreground emission spectrum in the WMAP bands on a pixel-by-pixel basis. Spatial templates of different emission components from external data are used as priors, and the model is designed to revert to the priors in regions of low signal-to-noise ratio. Thus the analysis is of most interest for separating and characterizing the different emission components in high signal-to-noise regions. The model foreground maps that are produced have complicated noise properties so they are not useful for foreground removal in cosmological analyses.

The nine-year MEM analysis differs from previous analyses (Bennett et al. 2003c; Hinshaw et al. 2007; Gold et al. 2009; Gold et al. 2011) in that spinning dust emission is treated as a separate emission component. Previously, synchrotron emission and spinning dust emission were treated together as a single component and an iterative method was used to solve for the spectrum of this component for each pixel.

The analysis is done using 1° smoothed beam-symmetrized nine-year sky maps in the five WMAP bands, with the ILC map subtracted from each map and conversion to antenna temperature applied. The zero level of each map is set such that a $\csc |b|$ fit, for HEALPix N_side = 512 pixels at b < −15° and outside of the KQ85 mask, yields a value of zero for the intercept. The maps are degraded to HEALPix N_side = 128 pixelization, and a model is fit for each pixel p by minimizing the function

$$\begin{equation} H = \chi ^2 + \lambda (p) \sum _c T_c(p) \ln \left[\frac{T_c(p)}{eP_c(p)}\right]. \end{equation} \tag{ 34 }$$

Here T_c and P_c are the model brightness and template prior brightness for foreground component c (e is the base of natural logarithms). The form of the second term ensures positivity of the solution T_c for each component, which alleviates degeneracy between the components. The parameter λ controls the relative weight of the data and the priors in the fit. As in previous analyses, we base λ(p) on the foreground signal strength: λ(p) = 0.6 [T_K(p)]^−1.5, where T_K(p) is the K-band ILC-subtracted map in mK antenna temperature.

The MEM foreground model is a sum of synchrotron, free–free, spinning dust, and thermal dust components. The adopted spectra for synchrotron, free–free, and thermal dust emission are fixed power laws with β = −3.0, − 2.15, and +1.8, respectively. The adopted synchrotron spectral index is consistent with measurements of K- to Ka-band spectral index from WMAP polarization data, for which free–free and spinning dust contributions are expected to be negligible. For spinning dust emission, we adopt a spectral shape predicted by the model of Ali-Haïmoud et al. (2009) and Silsbee et al. (2011). The top panel of Figure 18 compares predictions of this model for different interstellar environments. We adopt the spectral shape for their nominal cold neutral medium conditions. The bottom panel shows that the predicted shape does not vary much for different conditions if a multiplicative frequency shift is allowed for. The MEM model includes a frequency scale factor for the spinning dust spectrum for pixels where the spinning dust prior is brighter than 0.1 mK. This is constrained such that the peak frequency is in the range from 10 to 30 GHz. For other pixels, the peak frequency is fixed at 14.4 GHz, a typical value found for the Galactic plane region.

Zoom In Zoom Out Reset image size

The adopted priors are shown in Figure 15. The synchrotron prior is based on the 408 MHz map of Haslam et al. (1982). We use the original version of this map; our previous MEM analyses used a filtered version in which striping and point sources are suppressed. We add a zero level offset of 3.9 K, as suggested by Tartari et al. (2008) based on absolute measurements of sky brightness at 600 and 820 MHz. We subtract the 2.725 K CMB monopole and an extragalactic contribution of 12.96 K, from the analysis of ARCADE 2 and other data by Fixsen et al. (2011). The 408 MHz map is then scaled to form the prior in K-band using a spectral index of −2.9. (The ARCADE 2 extragalactic background is used instead of a source count based value such as 2.6 K from Gervasi et al. (2008) because it gives a prior that is more consistent with the csc b normalized K-band map at high latitudes.) The free–free prior is the scattering-corrected, extinction-corrected Hα template described in Section 5.3.1, scaled to free–free brightness temperature in K-band using 11.4 μK R⁻¹ (Bennett et al. 2003c). The spinning dust prior is the temperature-corrected dust map of Schlegel et al. (1998), scaled to spinning dust brightness temperature in K-band using 9.5 μK MJy⁻¹ sr. This is the slope of the correlation between the dust map and a map of spinning dust brightness from fits to Haslam et al. (1982) 408 MHz, Duncan et al. (1995) 2.4 GHz, and ILC-subtracted WMAP data in the Galactic plane. The thermal dust prior is the prediction of model 8 of Finkbeiner et al. (1999) at 94 GHz. All of the prior maps have been smoothed to 1° FWHM.

The adopted model provides good fits to the data without iterative adjustment of the synchrotron component spectrum as used in previous analyses. For pixels at |b| < 5°, absolute residuals are typically less than 0.01%, 0.34%, 1.2, 2.1%, and 0.7% in K-, Ka-, Q-, V-, and W-bands, respectively. Maps of the foreground components and peak frequency of spinning dust from the MEM analysis are shown in Figure 19.

Zoom In Zoom Out Reset image size

We again perform a pixel-based MCMC fitting technique to the five bands of WMAP data. Our method is similar to that of Eriksen et al. (2007), but we focus more on Galactic foregrounds rather than CMB. The fit results of the prior releases have been reproduced, with the "base" model, which uses three power-law foregrounds, producing virtually the same reduced χ² per pixel. We have again also included the 408 MHz map of Haslam et al. (1981) with a zero-point determined using the same $\csc |b|$ method as for the WMAP data.

There are two main changes from the previous release. The first is that the MCMC fit now uses the spinning dust spectrum for grains in a "cold neutral medium" as computed by Silsbee et al. (2011), with an optional frequency shift parameter described below. This change was made so that the MCMC fit uses the same spinning dust spectrum as the rest of the nine-year analysis. The second significant change is that the spinning-dust model is now run with the synchrotron spectral index as a free parameter. This was done to improve the quality of the fit, also discussed below.

The MCMC fit is performed on 1° smoothed maps downgraded to HEALPix N_side = 64. A MCMC chain is run for each pixel, where the basic model is

$$\begin{equation} T(\nu) = T_{s} \left(\frac{\nu }{\nu _K} \right)^{\beta _s(\nu)} + T_{f} \left(\frac{\nu }{\nu _K} \right)^{\beta _f} + a(\nu) T_\mathrm{cmb} + T_{d} \left(\frac{\nu }{\nu _W}\right)^{\beta _d} \end{equation} \tag{ 35 }$$

for the antenna temperature. The subscripts s, f, d stand for synchrotron, free–free, and dust emission, ν_K and ν_W are the effective frequencies for K- and W-bands (22.5 and 93.5 GHz), and a(ν) accounts for the slight frequency dependence of a 2.725 K blackbody using the thermodynamic to antenna temperature conversion factors found in Bennett et al. (2003c). The fit always includes polarization data as well, where the model is

$$\begin{equation} Q(\nu) = Q_{s} \left(\frac{\nu }{\nu _K} \right)^{\beta _s(\nu)} + Q_{d} \left(\frac{\nu }{\nu _W}\right)^{\beta _d} + a(\nu) Q_\mathrm{cmb} \end{equation} \tag{ 36 }$$

$$\begin{equation} U(\nu) = U_{s} \left(\frac{\nu }{\nu _K} \right)^{\beta _s(\nu)} + U_{d} \left(\frac{\nu }{\nu _W}\right)^{\beta _d} + a(\nu) U_\mathrm{cmb} \end{equation} \tag{ 37 }$$

for Stokes Q and U parameters. Thus there are a total of 15 pieces of data for each pixel (T, Q, and U for five bands).

As for the previous two releases, the noise for each pixel at N_side = 64 is computed from maps of N_obs at N_side = 512. To account for the smoothing process, the noise is then rescaled by a factor calculated from simulated noise maps for each frequency band. The MCMC fit treats pixels as independent, and does not use pixel-pixel covariance, which leads to small correlations in χ² between neighboring pixels. This has negligible effect on results as long as goodness of fit is averaged over large enough regions.

K-band is used as a template for the polarization angle of synchrotron and dust emission, so U_s and U_d are not independent parameters, identical to the previous analyses. All models also fix the free–free spectral index to β_f = −2.16, the same as in the seven-year release.

Results from the models discussed below are listed in Table 13; see the LAMBDA Web site for further details. The "base" model uses three power-law foregrounds, where the synchrotron spectral index β_s(ν) is taken to be independent of frequency but may vary spatially, and the dust spectral index β_d is allowed to vary spatially. We assume the same spectral indices for polarized synchrotron and dust emission as for total intensity emission. This model has a total of 10 free parameters per pixel: T_s, T_f, T_d, T_cmb, β_s, β_d, Q_s, Q_d, Q_cmb, and U_cmb.

For models with a spinning dust component, another term is added to Equation (35)

$$\begin{equation} T_{sd}(\nu) = T_{sd}(\nu _\mathrm{K}) S_{sd}(\nu), \end{equation} \tag{ 38 }$$

Where S_sd(ν) parameterizes the shape of the spinning dust spectrum, and is interpolated from values for the "cold neutral medium" spectrum given by Silsbee et al. (2011). An optional shift parameter can be used to rescale the frequency dependence before interpolation. This shift parameter applies to the full sky and does not vary per pixel. After shifting and interpolation, the spectrum S_sd(ν) is normalized to unity at K-band, leaving T_sd(ν_K) as the only spinning dust parameter for each pixel. Independent fits were performed to determine the best-fit shift parameter, which for the averaged sky was found to be 0.84. Inside the Kp12 mask (within a few degrees of the galactic plane) the preferred shift parameter may be somewhat lower (0.77), but the evidence is not strong.

The spinning dust component is assumed to have negligible polarization, as theoretical expectations for the polarization fraction are low compared to synchrotron radiation (Lazarian & Draine 2000), and polarization data thus far show no evidence that such a component is necessary (Section 5.3.2; López-Caraballo et al. 2011; Dickinson et al. 2011; Rubiño-Martín et al. 2012). This model then has 11 free parameters per pixel: the 10 parameters of the base model, plus the spinning dust amplitude.

MCMC fits for the nine-year release were performed with the addition of the 408 MHz data compiled by Haslam et al. (1981). The error on the zero point for this data was estimated in that work to be ±3 K, with an overall calibration error of 10%. As the MCMC method treats all input maps equally, for consistency we estimate and subtract a nominal zero point offset of 7.4 K, as determined by the same $\csc |b|$ method we use for the WMAP sky maps. For comparison, Lawson et al. (1987) used a comparison with 404 MHz data to find a uniform (presumably extragalactic) component with a brightness of 5.9 K.

We find that to best fit the 408 MHz data, the spinning dust spectrum needs to have its peak frequency adjusted downward by approximately 15%, similar to the case in the previous release. We also find that a much better fit is achieved in the plane by varying the synchrotron spectral index, which for that region allows a $\chi ^2_\nu = 1.24$ versus $\chi ^2_\nu = 1.76$ with fixed index, for 8.5 effective degrees of freedom. The mean spinning dust fraction inside the KQ85 mask is somewhat lower than in the seven-year fit, at 10% of 22 GHz flux compared to 18% in the seven-year fit.

We also find that the fit is improved by taking into account some mild steepening of the synchrotron spectrum from 408 MHz to WMAP's frequency range. Strong et al. (2011) have compared mid-latitude synchrotron measurements and estimates from 22 MHz to 94 GHz with predictions of cosmic ray propagation models based on cosmic ray and gamma ray data. We adopted their best-fit pure diffusion model ("galdef_ID_54_z04LMPD_g0_1.3_withsecS") to compute an effective synchrotron spectral index between 408 MHz and 23 GHz (WMAP K-band), as well as the index from 23 GHz to 94 GHz over which range it remains nearly constant. We calculate the difference in these two indices to be −0.12. Our model g (hereafter MCMCg, and listed on the last line of Table 13) then uses this difference, so that while the model parameter β_s is used as the synchrotron index for the WMAP bands, the value β_s + 0.12 is used to extrapolate the synchrotron component down to 408 MHz for comparison to the map of Haslam et al. The parameters from this fit are shown in Figure 20.

Zoom In Zoom Out Reset image size

In this section we attempt to find a best-fit foreground model that is consistent with both the WMAP data and Haslam. This is intended to be a faster fit than was done with the MCMC method in Section 5.3.5, and so it allows us to experiment with models more rapidly. Because this method simply finds the maximum likelihood point of the foreground model, it does not provide errors bars as the MCMC method does. Also, we avoid priors in the form of foreground component sky maps, which were used in the MEM fitting in Section 5.3.4. The priors we use in this section are mostly in the form of the foreground spectral shapes (relative antenna temperature as a function of frequency) instead of in the form of sky maps. This is a complementary form of analysis to the MEM fitting.

5.3.6.1. Data and noise

Our data consists of maps smoothed to a common resolution of 1° FWHM, which we pixelize at r6. We use six maps: 408 MHz and the five WMAP bands. We use the original Haslam map (408 MHz) as in Section 5.3.4 with the same offsets, except in this case we do not use the ARCADE extragalactic background. Instead of subtracting 12.96 K, we subtract 2.6 K (Tartari et al. 2008), so the Haslam map used in this section is 10.36 K brighter in antenna temperature in all pixels. The rms noise in each pixel of the 408 MHz map is taken to be 10% of the antenna temperature, added in quadrature with a 0.6 K uncertainty in zero point (Haslam et al. 1982; Tartari et al. 2008).

We consider three noise components for the WMAP bands in this foreground fitting: the 0.2% overall gain uncertainty, the $\sigma _0/\sqrt{N_{\rm obs}}$ instrument noise, and the uncertainty in the diffuse foreground monopole corrected with the $\csc |b|$ offsets, discussed previously. Because our fitting is done on a per-pixel basis, we approximate these errors as uncorrelated between pixels, and we add them in quadrature.

The instrument noise can be treated carefully to account for the smoothing to 1° FWHM. Typically it is inaccuracies in the foreground model that cause χ² to be large and not the details of the noise. However, a detailed treatment of the noise smoothed to 1° in r6 pixels is given in Appendix D. Again, because we fit on a per-pixel basis, we ignore the correlations in noise between nearby pixels.

5.3.6.2. Foreground models. We start with a simple foreground model consisting of several simple power laws, and progressively add complexity to the model to improve the fit. The foreground model we use involves temperature only; we did not try to fit polarization. The sequence of foreground models we use is listed in Table 14, and details are discussed below.

Table 14. χ² Minimal Prior Fits of Foreground Models

Model	Synchrotrona	Δβsyncb	ffc	βdust	SDd	SD Peake	νf	χ2/pixelg	fbadh
1	Power	0	Yes	1.8	No	⋅⋅⋅	4	6.1	37%
2	Power	Vary	Yes	1.8	No	⋅⋅⋅	5	2.5	11%
3	Power	Vary	Yes	1.6–2.0	No	⋅⋅⋅	6	2.3	9.5%
4	Power	Vary	Yes	1.8	Yes	15.1	6	1.5	4.4%
5	Power	Vary	Yes	1.8	Yes	12.5–17.8	7	0.64	0.59%
6	Strong	0	Yes	1.8	Yes	15.1	5	5.4	30%
7	Strong	0	Yes	1.8	Yes	12.5–17.8	6	4.1	20%
8	Strong	Vary	Yes	1.8	Yes	15.1	6	1.2	2.1%
9	Strong	Vary	Yes	1.8	Yes	12.5–17.8	7	0.60	0.48%

Notes. ^aWhether the synchrotron is treated as a pure power law or modeled according to a model from Strong et al. (2011). ^bFor both power law and Strong et al. synchrotron models, we either set the spatial variation in spectral index Δβ_sync to zero or allow it to vary: −0.5 ⩽ Δβ_sync ⩽ 0.5. In the case of a power law, Δβ_sync is a perturbation added to β_sync = −3.0. ^cThe free–free spectrum is given by Oster (1961); we use an electron temperature of 8000 K. ^dWhether a spinning dust spectrum in the shape of the cold neutral medium is used. ^eRange of available peak frequencies for the spinning dust spectrum, in GHz. This is either fixed at 85% of the peak frequency 17.8 GHz for the cold neutral medium (which is 15.1 GHz), or allowed to be a range from 70% to 100% of the CNM peak frequency (which is 12.5 GHz to 17.8 GHz). ^fDegrees of freedom in the model. Most degrees of freedom are constrained: foreground amplitudes must all be positive, for example. The highly constrained CMB amplitude is included as a degree of freedom. ^gThe mean χ² per pixel, averaged over the whole sky (for temperature only, not polarization), where χ² values greater than 10 are set to exactly 10 so that a few extremely bad pixels do not throw off the whole fit. This χ² value includes deviations of the model from Haslam and WMAP bands, but not deviations from the ILC prior. Since there are six measurements in each pixel (and an ILC prior) and 4 ⩽ ν ⩽ 7 degrees of freedom in the model, we would expect χ²/pixel ≈ 6 − ν for a good fit if we had unconstrained variables. ^hThe fraction of the pixels where χ² > 10. This is an estimate of the sky fraction where the fit is bad. Again, the χ² used here includes the difference of the model from the six bands, but does not include deviations from the ILC prior.

Download table as:  ASCIITypeset image

The synchrotron spectrum is either taken to be a pure power law in antenna temperature, $T_A \propto \nu ^{\beta _{\rm sync}}$, or derived from assuming the spectral index curve from Strong et al. (2011), Figure 6, upper right corner. This is the curve for a low-energy electron injection index of 1.3 and is the same spectrum as used in the MCMC fitting. To this spectral index curve we optionally add an offset in β_sync, −0.5 ⩽ Δβ_sync ⩽ 0.5 independent of frequency. We numerically integrate this spectral index curve to obtain synchrotron antenna temperature.

The free–free spectrum is the slightly curved model given by Oster (1961) and rearranged for antenna temperature by Bennett et al. (2003c). This is

$$\begin{equation} T_A^{{\rm WMAP}}(\nu) \,{\propto}\, \frac{1 \,{+}\, 0.2218 \ln (T_e / 8000\,{\rm K}) \,{-}\, 0.1479 \ln (\nu / 41\,{\rm GHz})}{(\nu / 41\,{\rm GHz})^2 (T_e / 8000\,{\rm K})^{1/2}}. \end{equation} \tag{ 39 }$$

For simplicity we use an electron temperature of 8000 K. We expect variations in electron temperature, but these do not strongly affect the shape of the spectrum.

The dust spectrum is given by a pure power law, typically with a fixed spectral index of β_dust = 1.8.

Finally, we add a spinning dust component. This is an antenna temperature spectrum from the Silsbee et al. (2011) model prediction for cold neutral medium (CNM) conditions, with an optional frequency scale factor. If the spectrum is plotted as antenna temperature as a function of log frequency, the frequency scale factor simply shifts the spectrum left or right. However, instead of quoting the frequency scale factor, we instead quote the peak frequency, when the spectrum is measured in antenna temperature as a function of frequency. The peak frequency of the CNM spectrum is 17.8 GHz.

All of these foregrounds are assumed to have a positive scale factor associated with them. Synchrotron, free–free, and spinning dust are normalized to K-band antenna temperature, and dust is normalized to W-band antenna temperature.

The CMB is modeled as a blackbody with constant thermodynamic temperature. To make the CMB fit look statistically isotropic, we add a prior that the CMB must be within 5 μK rms of the nine-year ILC. Without this prior, the data do not constrain the CMB very tightly in the galactic plane, and we find the CMB preferring values lower than −250 μK.

To approximate the finite width of the WMAP bandpasses, we calculate these spectra at three frequencies per band and determine the WMAP response from a weighted average, as described in Appendix E.

5.3.6.3. Fitting code. Fitting the foregrounds is a least squares problem. However, we modify the simple linear least squares problem in two ways: we constrain the coefficients, and we allow nonlinear foreground spectra. Constraining the coefficients is essential, because we know the foregrounds are always positive. Unconstrained least squares fitting will frequently give a very negative and therefore unphysical foreground. Secondly, we allow nonlinear foregrounds, in the sense that the total foreground is not simply a linear combination of fixed foreground spectra. We allow the spectra to vary, for example by allowing the synchrotron spectral index to be a fit parameter, or by allowing the peak frequency of spinning dust to be a fit parameter.

There are several codes which can be used to solve this problem. We have not made a thorough search of all available software, and we only considered code in IDL since that is the language in which much of our other software is written. We have found two codes to be useful: a bound variable least squares routine and a Levenberg–Marquardt solver.

We found a bound variable least squares (BVLS) routine²² to be very fast, but it is restricted to linear foreground models and so it cannot solve for varying spectral indices or spinning dust frequency scale parameters. Because of this constraint we do not use it to report results in this paper. However, this code does have the advantage that parameters can be constrained to be positive, so it can provide physically reasonable fits.

For the results reported in this section (in Table 14) we use the mpfitfun.pro routine,²³ which uses the Levenberg–Marquardt algorithm and was written by Craig Marquardt, for the constrained nonlinear least squares fitting. This is somewhat slower than the BVLS code because it cannot use the assumption that the χ² function is precisely quadratic in all of the fit coefficients. The ability to calculate foreground spectra quickly is an important factor in the speed of these calculations. We discuss a useful rapid method of calculating the integral over the WMAP bandpasses in Appendix E.

5.3.6.4. Results. The results of this simple foreground fitting are shown in the last columns of Table 14. Additionally, maps from the Model 9 fit are shown in Figure 21. A set of three fixed power laws in Model 1 does not fit the data well. Allowing spatial variation of the synchrotron power-law spectral index in Model 2 substantially improves this, but 11% of the sky is still fit very poorly. Allowing spatial variation of the dust spectral index in Model 3 does not substantially improve the number of well fit pixels, so we fix the dust spectral index to β = 1.8. Adding a spinning dust component with peak frequency of 15.1 GHz (which is 0.85 times the CNM peak frequency of 17.8 GHz) does improve the fit, and allowing that peak frequency to vary between 12.5 GHz and 17.8 GHz helps even more. See Models 4 and 5.

Zoom In Zoom Out Reset image size

Because it is probable that the synchrotron is not a pure power law and because we use the Haslam data at 408 MHz, which is much lower in frequency than the WMAP data, we test a curved synchrotron model from Strong et al. (2011). If we do not allow the spectral index to vary, we again get bad fits in Models 6 and 7. However, a varying spectral index combined with a spinning dust component produces results that are fractionally better than a pure power law with the same spinning dust components, as can be seen by comparing models 5 and 9, and comparing models 4 and 8.

None of these fits is perfect. Even in Model 9, there remain a few pixels that are not fit well. These are primarily in Ophiuchus, the galactic plane, and the Gum nebula.

5.3.7.1. Cross-comparison of foreground fits

Maps of parameters from the MEM, MCMCg, and six-band χ² Model 9 fits are shown in Figures 19, 20, and 21. A summary of the parameter treatment for each of these three models is provided in Table 15, and a high-latitude consensus decomposition is in Figure 22.

Zoom In Zoom Out Reset image size

Table 15. Summary of Foreground Decomposition Model Assumptions

Parameter	MEM	MCMCg	χ2 Model 9
βsynca	−3.0, fixed	Strong, |Δβ| < 0.5	Strong, |Δβ| < 0.5
βdust	+1.8, fixed	Free	+1.8, fixed
βff	−2.15, fixed	−2.16, fixed	−2.09 to −2.17, fixedb
$\nu _{\mathrm peak}^{\mathrm sd}$c	10–30, constrained	14.95, fixed	12.5–17.8, constrained
CMB	ILC subtracted	Free	ILC prior
Polarization data fit	No	Yes	No
External foreground spatial priors	Haslam, SFD, FDS, Hαd	No	No

Notes. ^aSynchrotron is assumed to be a power law with a fixed spectral index, β_sync, or a variable power law based on a Strong et al. (2011) model, with a best-fit value of Δβ added to the spectral index. ^bThe free–free spectrum for the χ² Model 9 fit is given by Oster (1961) with a fixed electron temperature T_e = 8000 K. The spectral index, β_ff = −2.14 at K-band and −2.17 at W-band. It increases to −2.09 at 408 MHz. ^cA spinning dust cold neutral medium spectral shape is used with an allowed range of a peak frequency shift, specified in GHz. ^dHaslam: the 408 MHz survey of Haslam et al. (1982); SFD: the temperature-corrected dust map of Schlegel et al. (1998); FDS: thermal dust model 8 from Finkbeiner et al. (1999); Hα: Hα all-sky mosaic from Finkbeiner (2003).

Download table as:  ASCIITypeset image

Results from these three models are a sampling of the possible parameter space which can be used to produce a total foreground model in each WMAP band. Each of these models possesses strengths and weaknesses, which can be used to offset one another. Included in these considerations are the treatment of the CMB component, the use of spatial priors, and the use of spectral constraints.

Treatment of the CMB. Both the MEM and Model 9 make use of the ILC as the CMB estimator: the MEM subtracts the ILC from the WMAP data before fitting, and Model 9 uses the ILC as a strong prior. As discussed in Section 5.3.7.2, the ILC is an imperfect estimate of the true CMB, containing a residual foreground bias signal. MCMCg, on the other hand, treats the CMB as a free parameter in its fit solution. While this is a strength for MCMCg at high latitudes, CMB-foreground covariance is strongest in the Galactic plane, and MCMCg does not separate the CMB and foregrounds well there. Use of the ILC provides a better constraint in that case.

Use of spatial priors. The MEM uses spatial templates to constrain its fitting solution at high latitudes where signal-to-noise is lower than in the Galactic plane. This produces a less noisy parameter solution at high latitudes when compared to the MCMCg and Model 9 χ² fit. This is valuable to the extent that one trusts those priors, both in terms of zero levels and spatial structure.

Use of spectral constraints. The synchrotron spectral index β_s is a pivotal parameter in model fitting, since its behavior influences the model allocation between synchrotron, free–free and spinning dust. The MEM assumes a constant value of β_s = −3 at WMAP frequencies. Model 9 and MCMCg allow each pixel to fit for this parameter independently, within the constraints of a Strong et al. (2011) spectral dependence. Positional gradients, including a latitudinal gradient, are probably closer to physical reality than a constant value (Kogut et al. 2007). However, with this degree of freedom comes the possibility for degeneracies with the free–free and spinning dust parameters. In Figure 23 we show results from a foreground degeneracy analysis for a representative pixel in the six-band Model 9 fit. There are significant degeneracies between parameter pairs that include either synchrotron amplitude or β_s. (A similar result was presented by Gold et al. (2009) for the five-year MCMC analysis, although that lacked a spinning dust component). We believe degeneracies are a factor in the appearance of the MCMCg and Model 9 β_s maps, which show a strong latitudinal gradient and a dust-like morphology in some regions, e.g., extending south of the plane over 150° < l < 190° and in the north celestial pole H i loop that extends north of the plane over 120° < l < 150° (Meyerdierks et al. 1991). All three models share a common spectral shape for the spinning dust spectrum. This shape is allowed to shift peak frequencies for MEM and Model 9, while MCMCg adopts a fixed peak frequency. Although the use of a common shape seems well motivated (see Figure 18), there is no guarantee that it is correct for all pixels. This is an additional source of uncertainty in the fits, as observational deviations from this shape will be distributed primarily among free–free and synchrotron components. We note an apparent power deficit in the Model 9 free–free map, present to a lesser extent in the MCMCg result, which is dust-like in signature. Finally, we note that all models assume a fixed β_ff, and only MCMCg allows for a free β_dust. These are less uncertain values, but errors in fixed values can ripple into other components.

Zoom In Zoom Out Reset image size

It is nevertheless possible to find relative agreement among these models, especially at higher latitudes. The high latitude foreground spectral components in the WMAP bands are shown in Figure 22 and all of the fitting techniques support this spectral decomposition of the foregrounds.

The actual foregrounds, especially at low Galactic latitudes, are clearly more complex than our parameterizations allow, since variations in physical conditions exist along any LOS. There are some sky locations that were not well fit even with all of the degrees of freedom allowed by the χ² fitting, such as in Ophiuchus. Given the complexity of the foreground emission mechanisms sampled by the WMAP bands, separating the CMB from the total observed foreground is a more straightforward and reliable process than the decomposition of that total foreground into physical components. Although we have found imperfections in the dust and free–free templates we use for foreground cleaning, those imperfections are primarily confined to regions which are masked from use in the cosmological analysis, and the use of foreground cleaned maps in the power spectrum analysis is robust.

A remaining item of interest is the microwave "haze." The first claim of a haze (Finkbeiner 2004) suggested an excess of free–free emission compared to the expectation from Hα, and was dubbed a "free–free haze." No longer believed to be free–free emission, its exact shape and attribution has evolved in the literature. In general the haze is described as an excess extended diffuse emission near the Galactic center. This excess appeared as a residual from the decomposition of WMAP K, Ka, and Q maps using external templates (Finkbeiner 2004; Dobler & Finkbeiner 2008). The templates most often used for this purpose are the Haslam 408 MHz map, a de-extincted form of the Finkbeiner (2003) Hα all-sky mosaic and the Finkbeiner et al. (1999) thermal dust models.

While the excess compared to external templates is clear, the attribution to a physical mechanism associated with Galactic emission is not. One interesting possibility characterizes the haze as a separate hard spectrum synchrotron component associated with the gamma-ray bubbles (Planck Collaboration IX 2013; Dobler et al. 2010). Planck Collaboration IX (2013) uses a Gibbs sampler to fit a foreground model outside a Galactic mask that assumes separate hard and soft power-law spectra. The cut-sky maps with these spectra are further decomposed, using external templates, into emission components with a distinct residual identified as a β_s ∼ −2.55 synchrotron haze. It is also possible to find reasonable models which adequately describe the data without the invocation of a haze component, as in e.g., Dickinson et al. (2009). In these cases, the haze excess is absorbed and distributed amongst other low frequency Galactic components. For example, a typical K-band haze intensity at roughly ±20° latitude near the Galactic center is ∼100 μK (Planck Collaboration IX 2013), whereas K-band residuals in those locations for the MEM, MCMCg, and Model 9 models are roughly zero with a 1σ deviation of a few μK. Existence of the haze as a separate spatial component is model dependent. It depends on foreground spectral assumptions, which affect the emission allocation between the CMB and the decomposition of the Galactic foregrounds into physical components. Because the haze is easily absorbed into other model components if not explicitly accounted for, and a number of remaining uncertainties exist in the morphology and behavior of low-frequency emissions in general (e.g., spinning dust), we feel this is a topic which remains open. Additional observations would be beneficial, especially at frequencies below K-band.

Although the thermal dust and free–free parameter amplitudes differ between the models presented here in details, there are clear common-mode similarities when they are compared against their externally derived equivalents (which we have used in Section 5.3.2 for template cleaning). Figure 24 illustrates these common-mode features by taking the mean parameter amplitudes from three models presented in this paper (MCMCg, MEM and chi-square fitting Model 9), and differencing them against their template counterparts. On the left in Figure 24 is the mean thermal dust amplitude at W-band minus the 94 GHz estimate derived from IRAS and COBE data by Finkbeiner et al. (1999). We have chosen to difference against their model 8, but a similar result is obtained for their other two-component dust model, model 7. In the Galactic plane, all of the three WMAP models show more emission in the outer plane and less in the inner plane than that predicted from the FDS models. A more quantitative representation of the planar differences is shown in Figure 25. Correlations between MEM, MCMCg, and Model 9 have roughly unity slopes, whereas correlations against FDS model 8 indicate FDS is brighter by up to ∼20% in high intensity regions in the inner Galaxy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The right-hand image in Figure 24 shows the difference between a mean K-band free–free emission estimate from the same three models in this paper and that from scattering-corrected de-extincted Hα using a conversion factor of 11.4 μK R⁻¹. Scatter between models in the plane generally disallows a definitive free–free mapping there. However, differences between the free–free emission predicted from Hα and the free–free model estimates in this paper consistently indicate that the Hα prediction is higher by roughly 20%–30% in the Gum and Orion regions. Free–free differences for the Gum away from the plane, where the optical depth is <1, can be explained by a low electron temperature for this region (Dickinson et al. 2003; Woermann et al. 2000). Differences for other regions are most likely due to errors in the extinction correction, since the assumption of uniformly mixed dust and gas may not be valid. Although W-band Galactic emission is primarily either from thermal dust or free–free, linear combinations of the FDS dust model and Hα predicted free–free have consistently been unable to describe the WMAP data in the plane; these apparent errors in both templates are consistent with those fitting errors.

5.3.7.2. ILC errors

Here we consider two types of error in the ILC: error due to CMB-foreground covariance, and error due to an incorrect estimate of the bias. See for example Hinshaw et al. (2007). These are errors which leave residual foreground signatures in the ILC estimate of the CMB.²⁴

The bias correction is directly related to the foreground model. To determine the ILC bias, we take maps of our foreground-only estimate (without CMB) in each of the five WMAP bands and construct an ILC directly. The specific attribution of the foregrounds to individual components (synchrotron, free–free, etc.) is not needed in this step; we only require maps of the total foreground in each band. If the foregrounds are sufficiently complex (if they are not a linear combination of 4 or fewer spectra in each region), then there will be residuals in this foreground-only ILC, and this is the ILC bias. The ILC bias consists of foregrounds that cannot be removed by any set of ILC weights. With enough diversity in foreground spectral components, we can find a linear combination of foreground spectra that mimics the CMB, and we cannot remove the CMB signature from the ILC by construction, because the ILC weights must sum to 1. To deal with the ILC bias, we construct a foreground model, compute the ILC bias, and subtract it directly from the ILC. Inaccuracies in the foreground model will translate to an incorrect subtraction of the ILC bias.

An estimate of the ILC bias was computed by Hinshaw et al. (2007) from simulations and three-year data. We revisit the bias computation using the Galactic emission estimates in the five WMAP bands from Model 9, MEM, and MCMCg. If these models perfectly describe the total Galactic emission at WMAP frequencies, then a bias map can easily be constructed by applying the flight ILC weights (given in Table 12) to these foreground maps. Such an application is shown in Figure 26. For comparison, Figure 26 also shows the bias correction from the three-year analysis, which is non-zero within the Kp2 mask and zero everywhere outside the mask.

Zoom In Zoom Out Reset image size

Close to the Galactic plane, the bias computed from the MCMCg model is larger than that for the other two models. Removal of this bias from the uncorrected WMAP data ILC shows a clear negative residual in the plane for |l| < 120°, indicating over-correction. In addition, ILC regional weights computed for the MCMCg model are sufficiently different from flight data values to render the model "goodness" suspect near the plane within the Kp2 cut. This is in part due to poorly constrained apportionment between CMB and Galactic signals in the plane. In particular there is an inverse correlation between CMB and dust spectral index, resulting in higher fractional residuals in portions of the plane for the MCMCg fit to V-band. V-band typically has the highest ILC weight, so these residuals lead to a higher bias for this model. Within the Kp2 cut, both Model 9 and the MEM bias maps show similar behavior to the three-year bias map, although details vary. Both models also return foreground ILC regional weights similar to data values, with the MEM showing the closest correspondence. Bias levels within the Kp2 cut are estimated from these two models as near 20 μK or less. These levels are either of similar magnitude or smaller compared to those computed for the CMB-foreground covariance in the same location (see below).

Estimating the foreground bias at higher latitudes is more difficult than for the Galactic plane regions. Since classic ILC weights are primarily determined using sky pixels within the Kp2 cut (even for the high latitude region 0), correspondence between derived model and data weights is only a useful diagnostic for pixels within the Kp2 mask. In addition, both the MEM and Model 9 results are ILC dependent: MEM subtracts the ILC from the data as a prelude to foreground fitting, and the six-band χ² Model 9 fit relies on the ILC as a strong prior. Since the classic ILC algorithm applies no bias correction outside the Kp2 cut, it is possible for any existing high-latitude ILC foreground bias to either remove or add power to the high latitude sky which is being fit to a Galaxy model. Since Galactic signals are generally weaker here than in the plane, the fractional error is potentially higher. Here the MCMC method provides the most objective model for estimating high latitude bias, since the CMB contribution is determined independently as part of the fitting process. We have used an amalgam of the three model bias maps to construct a very crude estimate of ILC bias outside of the Kp2 cut, giving the most weight to the MCMCg result. All three bias maps show a common characteristic dust-like excess in the outer Galaxy near the edges of the Kp2 cut. Two of the three bias maps show a low-level inner Galaxy deficit with a synchrotron-like signature. Noise in the bias maps makes a clear determination of the morphology difficult; we have used templates to represent the spatial structure, but the fine structural detail of the templates should not be taken as truth. Our rough estimate of the high latitude ILC bias is shown at the bottom right of Figure 26. High-latitude ILC bias is estimated at 10 μK or less.

The CMB-foreground covariance was discussed in Hinshaw et al. (2007). Because the ILC weights are constructed by minimizing the variance in a region, the weights adjust to allow foreground fluctuations to cancel CMB fluctuations as much as possible. This is more of a problem for small regions. Because the total foreground level is well measured in the plane (even if we allow complete uncertainty in the CMB for an error term of σ ≈ 70 μK, the foregrounds are bright enough to make this term small), we can estimate how much the foregrounds could correlate with a random CMB sky with a given power spectrum. This estimate will not change substantially with different foreground models (different estimates of how much of the WMAP data is CMB and how much is foreground) because it only requires knowledge of the total foreground level, which is well constrained by the data. We can experimentally determine the CMB-foreground covariance by generating many CMB simulations, adding a foreground model to each CMB simulation, making a bias-subtracted ILC, and forming an error map by subtracting the true CMB from the ILC in each simulation. This gives us an ensemble of error maps, which span a 48 dimensional space. Since the CMB simulation is perfectly subtracted by any set of weights that add to 1, our error maps contain no CMB from the simulation. They only contain errors from residual foregrounds. Since there are 60 weights (going into the 12 regions of the ILC) and 12 constraints where sets of weights must add to 1, there are 48 degrees of freedom in the ILC error. As with the ILC bias, the results do depend on foreground model, but not nearly as strongly, as mentioned above.

We construct the 48 maps showing the ILC foreground-CMB covariance modes at res 6 as follows. We take the foreground Model 9 from Section 5.3.6 and prograde it directly to r9 (with no extra smoothing), where the ILC regions are defined. Then we form ILCs by the usual method, except that we do not smooth between regions as described in Equation (18) of Hinshaw et al. (2007) because we next degrade back to r6, which has a similar effect. We do this for 1000 CMB realizations, and form a 49152 × 1000 matrix of the maps, of which we take a singular value decomposition to determine the most common modes, taking care to normalize properly. There are only 48 singular values that are not effectively zero; we use the 1000 simulations to better sample these 48 modes and better determine their eigenvalues.

These modes provide the eigenvalues with nonzero eigenvectors of the foreground-CMB covariance error matrix. We compute the square root of the diagonal elements of this matrix to provide a visual estimate (that ignores correlations) of this error. The nine-year ILC map and this error map are shown in Figure 27.

Zoom In Zoom Out Reset image size

We demonstrate the use of this error description by propagating the foreground-CMB error to the quadrupole–octupole alignment, which we describe in Section 7.4.

5.3.7.3. ILC considerations

The primary difficulty with any method of extracting the CMB from the data is determining how much of the temperature in each pixel is foreground and how much is CMB. The data only constrain the sum of these two, and we must make other assumptions in order to separate them.

The ILC specifically assumes that the CMB has a blackbody spectrum while the foregrounds do not. In addition, the ILC assumes that while the foregrounds may change amplitude across a region, an individual foreground does not change its spectral shape (proportional to antenna temperature as a function of frequency), so that a set of ILC weights can null a given foreground everywhere in a region. Along with this, the ILC assumes that there are four or fewer foreground spectral shapes, since if there were more, we would not be able to remove them all with only the five bands of WMAP data. If there were five foreground spectra, some linear combination of them would be able to mimic a blackbody spectrum, which the ILC has been designed to keep.

Figure 28 is one way to visualize the foreground complexity of the WMAP data. It shows in color the regions that are approximate power laws, and it shows in grayscale regions that are not well fit by a single power law. The ILC methodology can handle more than a single power law foreground (it can remove up to four of them), so this is not directly a map of where the ILC will work well. However, this figure does show the varying nature of foreground spectra across the sky.

Zoom In Zoom Out Reset image size

Choosing the ILC region size is a trade-off between foreground complexity and foreground-CMB covariance. By choosing small regions, we give the foregrounds less chance to vary their shape over a region (such as by changing a synchrotron spectral index). But small regions are more susceptible to foreground-CMB covariance, as discussed in Hinshaw et al. (2007), which suppresses the variance of the ILC to the extent that the foregrounds and CMB correlate.

We could, for example, take minimum variance to be our figure of merit for an ILC map and allow arbitrary gerrymandering of the regions on a pixel-by-pixel basis. This could be done with a simulated annealing algorithm adjusting some small number of regions (e.g., 4) within a galactic mask. However, this would result in an ILC with variance inside the mask well below the expected CMB variance, because the regions optimize the foreground-CMB covariance to artificially suppress the ILC fluctuations. More knowledge than just the ILC variance is needed for intelligent region selection.

The foreground-CMB covariance can be estimated moderately well, since it only depends on an approximate foreground model and knowledge of the CMB power spectrum. We estimate this error in Section 5.3.7.2 and propagate it to the quadrupole–octupole alignment in Section 7.4. Other errors, such as those due to foregrounds changing spectral shape over a region or more than 4 foreground spectra in a region (these cause the ILC bias), are harder to estimate because they require an accurate separation of CMB from foregrounds in the first place. The demands on this foreground model accuracy depend on the amplitude of the foregrounds. For a pixel dominated by CMB, a slight foreground correction need not be extraordinarily accurate in a fractional sense. Yet for an extremely bright foreground location on the plane (say, a bright H ii region), the foreground model must have supreme fractional accuracy to distinguish meaningfully a tiny CMB contribution from the dominating foregrounds.

A more accurate ILC would require either a better bias subtraction or better region selection designed to minimize the needed bias correction; both of these require a highly accurate foreground model. A foreground model that separates out different components (such as synchrotron, free–free, etc.) is not needed, only a model that gives the total foreground in each band. The ILC bias can be directly calculated by making an ILC of this foreground-only data set, and regions could be selected to minimize the bias correction needed in each region. However, if we already have an accurate separation of the CMB from foregrounds, then the ILC method is no longer necessary, since we already have a map of the CMB.

In our power spectrum pipeline, we precompute three quantities: a transfer matrix $F_{ll^{\prime }}$ that represents the mean response of the unnormalized estimator at multipole l to CMB power at multipole l'; the bias of the power spectrum estimator due to unresolved point sources; and the noise bias, for the auto-correlation estimator ${\widehat{C}}_l$ (but not for the cross-correlation estimator ${\widehat{C}}_l^\times$). In Figure 29, we present end-to-end Monte Carlo tests of these precomputations using three simulated ensembles: CMB-only simulations, point source simulations, and noise-only simulations. In all cases the ratio of the recovered power spectrum (averaged over many Monte Carlo realizations) to the expected power spectrum is consistent with unity.

Zoom In Zoom Out Reset image size

Our pipeline uses interpolation in l to estimate transfer matrices, noise bias, and point source bias. We did an end-to-end test of the interpolation accuracy as follows. We reran the pipeline with half the interpolation step size, treated the difference between the two estimates as a power spectrum bias, and then we did a Fisher matrix forecast to determine whether the resulting bias was statistically significant. In all three cases, we found that the resulting bias is ≲ 0.02σ, i.e., much too small to be important.

We estimate the power spectrum covariance matrix $\mbox{Cov}({\widehat{C}}^\times _l, {\widehat{C}}^\times _{l^{\prime }})$ using Monte Carlo simulations. A direct Monte Carlo estimation of a 1200×1200 covariance matrix would require a prohibitive number of simulations, but this can be sped up using computational tricks: (1) the covariance $\mbox{Cov}({\widehat{C}}_l,{\widehat{C}}_{l^{\prime }})$ of the auto-estimator is equal to the inverse Fisher matrix $F_{ll^{\prime }}^{-1}$, so we only need Monte Carlos for the estimator difference $({\widehat{C}}_l^\times -{\widehat{C}}_l)$; (2) we only estimate variances and assume that off-diagonal covariances are given by appropriately rescaling Fisher matrix elements; and (3) we smooth the variance estimates in l. These tricks allow the covariance matrix to be accurately estimated from a small number of simulations. As an end-to-end convergence test, we compared covariance matrices C₂₅₆, C₅₁₂ constructed using 256 and 512 Monte Carlo simulations respectively. We found that all matrix entries were nearly identical in that all Karhunen–Loève eigenvalues of the matrix pair (C₂₅₆, C₅₁₂) are between 0.999 and 1.001.

In Figure 30, we show the binned power spectrum estimates for the seven-year WMAP data obtained using the optimal pipeline, described above, with the sub-optimal MASTER results used in the seven-year WMAP release (Larson et al. 2011) shown for comparison. The agreement is excellent; the two estimators agree to better than 1σ in every l-bin, as expected when comparing an optimal and near-optimal analysis of the same data.

Zoom In Zoom Out Reset image size

To compare the two estimators more closely, in the left panel of Figure 31 we show the difference between the optimal and sub-optimal estimators, before and after smoothing in l. No systematic trends are seen, as expected if the difference is pure statistical scatter. There is a small region near l = 50 where the optimal estimator fluctuates to a lower value of C_l than the sub-optimal estimator. This fluctuation slightly shifts the best-fit value of the spectral index n_s, as discussed by Hinshaw et al. (2013). This appears to be the most important difference between the two estimators for purposes of cosmological parameter estimation, aside from the effective sensitivity improvement discussed below.

Zoom In Zoom Out Reset image size

The right panel of Figure 31 shows the ratio between the power spectrum variance Var(C_l) obtained using the optimal and sub-optimal estimators. The optimal estimator improves the variance by 7%–17% depending on the value of l. This level of improvement is roughly comparable to the improvement in going from seven-year to nine-year data (which varies from no improvement at low l to a factor of 9/7 = 1.28 in C_l at high l).

We will limit our search for non-Gaussianity to the three-point function or bispectrum, and parameterize it by

$$\begin{eqnarray} \langle \zeta _{{\bf k}_1} \zeta _{{\bf k}_2} \zeta _{{\bf k}_3} \rangle &=& \left(f_{{\rm NL}}^{\rm loc}B_{\rm loc}(k_1,k_2,k_3) + f_{{\rm NL}}^{\rm eq}B_{\rm eq}(k_1,k_2,k_3) \right. \nonumber\\ &&\left. + f_{{\rm NL}}^{\rm orth}B_{\rm orth}(k_1,k_2,k_3) \right) (2\pi)^3 \delta ^3\left(\sum {\bf k}_i \right)\qquad \end{eqnarray} \tag{ 50 }$$

where $f_{{\rm NL}}^{\rm loc}$, $f_{{\rm NL}}^{\rm eq}$, $f_{{\rm NL}}^{\rm orth}$ are free parameters to be estimated, and the local, equilateral, and orthogonal template bispectra are defined by

$$\begin{eqnarray} B_{\rm loc}(k_1,k_2,k_3) &=& \frac{6}{5} (P_\zeta (k_1) P_\zeta (k_2) + \mbox{2 perm.} ) \end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray} &&B_{\rm eq}(k_1,k_2,k_3) = \frac{3}{5}(6 P_\zeta (k_1) P_\zeta (k_2)^{2/3} P_\zeta (k_3)^{1/3} \nonumber \\ &&\quad - 3 P_\zeta (k_1) P_\zeta (k_2) - 2 P_\zeta (k_1)^{2/3} P_\zeta (k_2)^{2/3} P_\zeta (k_3)^{2/3} + \mbox{5 perm.} )\nonumber\\ \end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray} &&B_{\rm orth}(k_1,k_2,k_3) = \frac{3}{5} (18 P_\zeta (k_1) P_\zeta (k_2)^{2/3} P_\zeta (k_3)^{1/3} \nonumber \\ &&\quad\! - 9 P_\zeta (k_1) P_\zeta (k_2) - 8 P_\zeta (k_1)^{2/3} P_\zeta (k_2)^{2/3} P_\zeta (k_3)^{2/3} + \mbox{5 perm.} ).\nonumber\\ \end{eqnarray} \tag{ 53 }$$

The $\lbrace f_{{\rm NL}}^{\rm loc}, f_{{\rm NL}}^{\rm eq}, f_{{\rm NL}}^{\rm orth}\rbrace$ basis for the three-point function is large enough to encompass a range of interesting models. Local-type non-Gaussianity is generic to some multi-field inflation models, for example curvaton models (Linde & Mukhanov 1997; Lyth et al. 2003) and variable reheating models (Dvali et al. 2004; Zaldarriaga 2004), and also to some alternatives to inflation, such as "new" ekpyrosis (Creminelli & Senatore 2007; Buchbinder et al. 2007) and cyclic (Lehners & Steinhardt 2008a, 2008b) models. Also, there is a theorem (Creminelli & Zaldarriaga 2004) that implies that no single-field model of inflation can generate detectable $f_{{\rm NL}}^{\rm loc}$. Equilateral-type and orthogonal-type non-Gaussianity can be generated in single-field models, and generically appear when there are non-negligible interaction terms in the inflationary Lagrangian.

We constrain the f_NL parameters using the optimal (i.e., minimum variance unbiased) bispectrum estimator implemented in Smith et al. (2009), which builds on previous work (Komatsu et al. 2005; Creminelli et al. 2006; Smith & Zaldarriaga 2011). The estimator optimally combines channels with different noise maps and beams by filtering the data with the inverse signal+noise covariance C⁻¹ = (S + N)⁻¹, and includes a one-point term (in addition to a three-point term) which reduces the variance. Unless otherwise specified, we use the V-band and W-band DAs from WMAP (six maps total), remove regions of high Galactic foreground and point source emission using the nine-year KQ75 mask, and marginalize three foreground templates corresponding to synchrotron, free–free, and dust emission. With foreground marginalization enabled, the same f_NL estimates are obtained on raw and template-cleaned maps.

Our "bottom line" constraints on non-Gaussianity are as follows:

$$\begin{eqnarray} f_{{\rm NL}}^{\rm loc}= 37.2 \pm 19.9 & \quad & \left({-}3 < f_{{\rm NL}}^{\rm loc}< 77 \mbox{ at 95% CL}\right) \nonumber \\ f_{{\rm NL}}^{\rm eq}= 51 \pm 136 & \quad & \left({-}221 < f_{{\rm NL}}^{\rm eq}< 323 \mbox{ at 95% CL}\right) \nonumber \\ f_{{\rm NL}}^{\rm orth}= -245 \pm 100 & \quad & \left({-}445 < f_{{\rm NL}}^{\rm orth}< -45 \mbox{ at 95% CL}\right). \nonumber\\ \end{eqnarray} \tag{ 54 }$$

The $f_{{\rm NL}}^{\rm loc}$ constraint includes a correction for the ISW-lensing contribution to the bispectrum, which arises from the large-scale correlation between the CMB temperature and the CMB lensing potential. We find that the ISW-lensing bispectrum biases the $f_{{\rm NL}}^{\rm loc}$ estimator by $\Delta f_{{\rm NL}}^{\rm loc}= 2.6$; this bias has been subtracted from the estimate in Equation (54). The ISW-lensing bias was computed using the Fisher matrix approximation, but this has been shown to be an excellent approximation to the exact result (Hanson et al. 2009; Lewis et al. 2011).

The constraint on each f_NL parameter in Equation (54) assumes that the other two f_NL parameters are zero. For a joint analysis of all three parameters, we need the bispectrum Fisher matrix:

$$\begin{equation} F = \left(\begin{array}{@{}ccc@{}}25.25 & 1.06 & -2.39 \\ 1.06 & 0.54 & 0.20 \\ -2.39 & 0.20 & 1.00 \end{array} \right) \times 10^{-4} \end{equation} \tag{ 55 }$$

where the ordering of the rows and columns is $f_{{\rm NL}}^{\rm loc}, f_{{\rm NL}}^{\rm eq}, f_{{\rm NL}}^{\rm orth}$. The statistical error on each f_NL parameter in Equation (54), with the other two f_NL parameters fixed to zero, is (F_ii)^−1/2, and the correlation between two estimators in Equation (54) is equal to the rescaled off-diagonal matrix element F_ij/(F_iiF_jj)^1/2.²⁵ An example of a two-parameter joint analysis is shown in Figure 42 below.

The most striking result in Equation (54) is the estimate for $f_{{\rm NL}}^{\rm orth}$, which is non-zero at 2.45σ. The (two-sided) probability of obtaining a value with this statistical significance in a Gaussian fiducial cosmology is 1.4%. This is not significant enough by itself to consider it a detection, but even further caution is required. When interpreting this probability, it must be kept in mind that we look for multiple deviations from the vanilla ΛCDM model,²⁶ so it is statistically unsurprising that one such deviation is at this significance level. The rest of this section will be devoted to consistency checks and interpretation of the $f_{{\rm NL}}^{\rm orth}$ result.

One possible source of systematic error is contamination by residual foregrounds. Since we marginalize over synchrotron, free–free and dust templates in our bispectrum estimator, any foreground contribution that is a linear combination of these spatial templates does not contribute to $f_{{\rm NL}}^{\rm orth}$. However, since the templates are not perfect, there will be residual contributions at some level. A simple procedure that gives the rough order of magnitude is to disable template marginalization in the estimator, and compute the foreground contribution to $f_{{\rm NL}}^{\rm orth}$ in an ensemble of simulated raw maps without any foreground cleaning. We simulate raw maps using random CMB and noise realizations, and a fixed dust realization given by model 8 of Finkbeiner et al. (1999). We do not include synchrotron and free–free foregrounds since dust dominates in W-band and is a significant fraction of the V-band foreground. In each simulation, we compute the difference $(\Delta f_{{\rm NL}}^{\rm orth})$ between the $f_{{\rm NL}}^{\rm orth}$ estimate obtained from the raw map, and the $f_{{\rm NL}}^{\rm orth}$ estimate that would be obtained from the CMB+noise contribution alone. We find that the mean value of $(\Delta f_{{\rm NL}}^{\rm orth})$ is 1.1 and the RMS scatter is 5.2. This presumably overestimates the dust contribution since we are not attempting to remove foregrounds at all. Since the shift $(\Delta f_{{\rm NL}}^{\rm orth})$ seen in these simple simulations is much smaller than the statistical error $\sigma (f_{{\rm NL}}^{\rm orth})$, we conclude that residual foregrounds are unlikely to be a significant contaminant.

As a first test for instrumental systematic effects, we check for consistency between different angular scales by splitting the $f_{{\rm NL}}^{\rm orth}$ estimator in l-bands. Our procedure is as follows: we write the $f_{{\rm NL}}^{\rm orth}$ estimator as a sum over triangles, restrict the sum to triangles whose maximum multipole max (l₁, l₂, l₃) is in a given bin (l_min, l_max), and then appropriately normalize so that the band-restricted sum is an unbiased estimator of $f_{{\rm NL}}^{\rm orth}$. This prescription for binning the $f_{{\rm NL}}^{\rm orth}$ estimator has the property that if we combine $f_{{\rm NL}}^{\rm orth}$ estimates in all bins up to some multipole l_max, the result agrees with simply rerunning the $f_{{\rm NL}}^{\rm orth}$ estimator with maximum multipole l_max. It also has the property that $f_{{\rm NL}}^{\rm orth}$ estimates in different l-bands are nearly uncorrelated.

In Figure 40, we show the $f_{{\rm NL}}^{\rm orth}$ estimate in l-bands, with the cumulative best-fit value $f_{{\rm NL}}^{\rm orth}=-245$ shown for comparison. Each bin is consistent with the cumulative best-fit value at 2σ, and the overall χ² of the fit to a constant $f_{{\rm NL}}^{\rm orth}$ value is good (χ² = 8.8 with seven degrees of freedom). We therefore conclude that there is no evidence for scale-dependent systematic contamination.

Zoom In Zoom Out Reset image size

As a second test for systematics, we can ask whether estimates of $f_{{\rm NL}}^{\rm orth}$ in different parts of the sky are consistent. The bispectrum estimator is naturally written as an integral over position on the sky, so a convenient way to visualize the position dependence is to simply plot the integrand as a skymap (Figure 41). This skymap is in units of "$f_{{\rm NL}}^{\rm orth}$ per steradian" and has the property that its integral over the whole sky is precisely equal to the estimated $f_{{\rm NL}}^{\rm orth}=-245$. If we restrict the integral to a subregion Ω of the sky, the value of the integral will roughly equal the value that would be obtained if we re-ran the estimator using masking to isolate the subregion Ω (appropriately rescaled by the area of Ω). Visual inspection of the skymap is a convenient way to look for an unexpected feature (e.g., a large contribution near the Galactic plane would suggest foreground contamination), although it might be difficult to assess the statistical significance of an a posteriori feature if found. Our interpretation of Figure 41 is that no visually striking features are seen; the skymap looks qualitatively similar to skymaps obtained from Gaussian simulations.

Zoom In Zoom Out Reset image size

As a more quantitative test for consistency between different parts of the sky, we estimated $f_{{\rm NL}}^{\rm orth}$ in the portions of the following regions that lie outside the KQ75 mask: the northern Galactic hemisphere, the southern Galactic hemisphere, within 30° of the ecliptic plane, and the ecliptic poles (>30° from the ecliptic plane). We find that for any pair of these regions, the estimated $f_{{\rm NL}}^{\rm orth}$ values are consistent at 2σ, relative to an ensemble of Monte Carlo simulations. The $f_{{\rm NL}}^{\rm orth}$ estimates in these four subregions are −139 ± 139, −361 ± 142, −132 ± 144, and −336 ± 138, respectively.

As a final test for systematics, we can compare $f_{{\rm NL}}^{\rm orth}$ estimates from different channels, or combinations of channels. In the first two columns of Table 16, we show the result of applying the $f_{{\rm NL}}^{\rm orth}$ estimator for several combinations of channels. To assess whether the $f_{{\rm NL}}^{\rm orth}$ estimates from a given pair of rows are statistically consistent, we subtract the two estimates, and compare the result to the same quantity (the difference of two $f_{{\rm NL}}^{\rm orth}$ estimates) evaluated in an ensemble of Monte Carlo simulations. This way of assessing consistency fairly incorporates the correlation between $f_{{\rm NL}}^{\rm orth}$ estimates that arises because the CMB realization (and the noise realizations, if the two rows have channels in common) is shared. The matrix in the rightmost columns of Table 16 shows the result of doing this consistency test for all pairs of rows in the table.

Table 16. A Test for Consistency between Channels

Channels	$f_{{\rm NL}}^{\rm orth}$	Discrepancy in "Sigmas"
		VW	V	W	V1	V2	W1	W2	W3	W4
All VW channels	−245.5 ± 99.6	⋅⋅⋅	2.2σ	1.5σ	1.5σ	2.1σ	0.7σ	1.4σ	1.1σ	2.2σ
All V-band channels	−125.9 ± 112.7	2.2σ	⋅⋅⋅	2.3σ	0.1σ	0.7σ	0.4σ	0.3σ	0.1σ	1.1σ
All W-band channels	−320.2 ± 112.1	1.5σ	2.3σ	⋅⋅⋅	2.1σ	2.5σ	1.7σ	2.2σ	2.0σ	3.2σ
V1 only	−119.3 ± 129.1	1.5σ	0.1σ	2.1σ	⋅⋅⋅	0.3σ	0.5σ	0.3σ	0.1σ	1.0σ
V2 only	−91.3 ± 124.2	2.1σ	0.7σ	2.5σ	0.3σ	⋅⋅⋅	0.8σ	0.0σ	0.2σ	0.8σ
W1 only	−172.1 ± 140.1	0.7σ	0.4σ	1.7σ	0.5σ	0.8σ	⋅⋅⋅	0.7σ	0.5σ	1.4σ
W2 only	−88.1 ± 152.2	1.4σ	0.3σ	2.2σ	0.3σ	0.0σ	0.7σ	⋅⋅⋅	0.2σ	0.7σ
W3 only	−111.0 ± 154.2	1.1σ	0.1σ	2.0σ	0.1σ	0.2σ	0.5σ	0.2σ	⋅⋅⋅	0.9σ
W4 only	−5.7 ± 147.7	2.2σ	1.1σ	3.2σ	1.0σ	0.8σ	1.4σ	0.7σ	0.9σ	⋅⋅⋅

Notes. The first two columns show $f_{{\rm NL}}^{\rm orth}$ estimates obtained from different subsets of WMAP channels. The matrix on the right shows the level of discrepancy between each pair of channel subsets, in "sigmas" after comparing to an ensemble of Monte Carlo simulations.

Download table as:  ASCIITypeset image

This "two-way" null test can be generalized to an N-way null test that tests mutual consistency between $f_{{\rm NL}}^{\rm orth}$ estimates obtained in all N rows of the table. We represent the $f_{{\rm NL}}^{\rm orth}$ estimates as a length-N vector f_i, and compute the N-by-N covariance matrix C_ij using Monte Carlo simulations with shared CMB and noise realizations. We then compute an overall best-fit $f_{{\rm NL}}^{\rm orth}$ value F which minimizes $\chi ^2 = (f_i-F) C^{-1}_{ij} (f_j-F)$. If the N estimates are mutually consistent, then the value of χ² at the minimum will be distributed as a χ² random variable with (N − 1) degrees of freedom.

We find that the channel-channel null tests are marginal. The N-way null test gives χ² = 16.3 with 8 degrees of freedom, corresponding to one-sided probability p = 0.038. The most discrepant pair of rows in Table 16 is (W, W4), which differ by 3.2σ relative to Monte Carlo simulations. This statistical significance should not be taken at face value since there are 36 matrix entries in Table 16, and we have chosen the most anomalous one. However, if we construct the same matrix for each member of an ensemble of simulations, we find that the probability that at least one pair of rows is discrepant by >3.2σ is 2.6%. Finally, we observe that the discrepancy between V-band and W-band channels, which is in some sense the most natural split, is 2.3σ, corresponding to probability p = 0.021.

We conclude that there is some tension in the channel–channel null tests, with p-value around a few percent depending on which test is chosen. Since we have also considered null tests that pass cleanly (i.e., the tests based on scale dependence and sky location), our interpretation is that one failure at the few-percent level does not indicate systematic contamination, although the discrepancy between V-band and W-band is of some concern. We therefore cautiously proceed to discuss the physical implications of the non-Gaussianity constraints.

We opt to work in the context of single-field inflation, and use the effective field theory developed in Cheung et al. (2008a, 2008b). The EFT provides a master Lagrangian which is general enough to describe almost all single-field models of inflation. See also Gruzinov (2005); Chen et al. (2007). The action consists of a standard kinetic term, plus small interaction terms whose coefficients parameterize allowed non-Gaussianity:

$$\begin{eqnarray} S &=& \int d^4x\, \sqrt{-g} \left[ -\frac{M_{\rm Pl}^2 \dot{H}}{c_s^2} \left(\dot{\pi }^2 - c_s^2 \frac{(\partial _i\pi)^2}{a^2} \right)\right.\nonumber\\ &&\left. +\, \left(M_{\rm Pl}^2 \dot{H}\right) \frac{1-c_s^2}{c_s^2} \left(\frac{\dot{\pi }(\partial _i\pi)^2}{a^2} + \frac{A}{c_s^2} \dot{\pi }^3 \right) + \cdots \right] \end{eqnarray} \tag{ 56 }$$

Non-Gaussianity is parameterized by a dimensionless sound speed c_s, and a dimensionless parameter A that represents the ratio between the coefficients the operators of $\dot{\pi }^3$ and $\dot{\pi }(\partial _i\pi)^2$. We treat c_s and A as free parameters, but specific models will make predictions. For example, in DBI inflation (Alishahiha et al. 2004), c_s is a free parameter (but related to the tensor-to-scalar ratio) and A = −1.

The coefficients in the action (56) can be related to the parameters $f_{{\rm NL}}^{\rm eq}$, $f_{{\rm NL}}^{\rm orth}$ by calculating the bispectra generated by the cubic operators $\dot{\pi }^3$ and $\dot{\pi }(\partial _i\pi)^2$, and projecting them onto the basis of template bispectra (Senatore et al. 2010). The result is:

$$\begin{eqnarray} f_{{\rm NL}}^{\rm eq}&=& \frac{1-c_s^2}{c_s^2} (-0.276 + 0.0785 A) \nonumber \\ f_{{\rm NL}}^{\rm orth}&=& \frac{1-c_s^2}{c_s^2} (0.0157 - 0.0163 A) \end{eqnarray} \tag{ 57 }$$

where the numerical coefficients are specific to the nine-year WMAP results and have been computed using the exact Fisher matrix, including CMB transfer functions and WMAP noise properties. For generic values of A, $f_{{\rm NL}}^{\rm eq}$ is larger than $f_{{\rm NL}}^{\rm orth}$ (by an order of magnitude) and equilateral non-Gaussianity is generated. However, there is an order-unity window of values (roughly 3.1 ≲ A ≲ 4.2) where $f_{{\rm NL}}^{\rm orth}$ is larger than $f_{{\rm NL}}^{\rm eq}$, and orthogonal non-Gaussianity is generated.

Since single-field models that produce $f_{{\rm NL}}^{\rm orth}$ are also expected to produce $f_{{\rm NL}}^{\rm eq}$ at some level, it is natural to analyze joint constraints in the two-parameter space $\lbrace f_{{\rm NL}}^{\rm eq}, f_{{\rm NL}}^{\rm orth}\rbrace$. To set up a joint analysis, we define notation as follows. Let $f_i = (f_{{\rm NL}}^{\rm eq},f_{{\rm NL}}^{\rm orth})$ be a two-component vector containing model parameters, let ${\hat{f}}_i = (51,-245)$ be the values of the associated estimators (i.e., the last two rows of Equation (54)), and let F_ij be the associated 2 × 2 Fisher matrix (i.e., the lower right corner of Equation (55). Then for given model parameters f_i, we define a χ² statistic,

$$\begin{equation} \chi ^2 = \sum _{ij} f_i F_{ij} f_j - 2 \sum _i F_{ii} f_i \hat{f}_i + \sum _{ij} {\hat{f}}_i F_{ii} F_{ij}^{-1} F_{jj} {\hat{f}}_j. \end{equation} \tag{ 58 }$$

We threshold this χ² to obtain confidence regions in the $(f_{{\rm NL}}^{\rm eq},f_{{\rm NL}}^{\rm orth})$ plane. These confidence regions are shown in the left panel of Figure 42. We note that the point $(f_{{\rm NL}}^{\rm eq},f_{{\rm NL}}^{\rm orth})=0$ is just outside the 2σ contour, which means that it is just barely a >2σ event when $f_{{\rm NL}}^{\rm eq}$ is included in the parameter space. More precisely, the relevant Δχ² is 7.16 with two degrees of freedom; the probability of getting a Δχ² this large in a Gaussian cosmology is 2.8%.

Zoom In Zoom Out Reset image size

In the right panel of Figure 42, we change variables to show confidence regions in the parameter space (c_s, A). These confidence regions were obtained under the assumption that the single-field bispectra are well-approximated by the equilateral and orthogonal template shapes. However, we have checked that nearly identical confidence regions are obtained if the exact tree-level bispectra for the operators ${\dot{\pi }}^3$ and $\dot{\pi }(\partial _i \pi)^2$ are used throughout the analysis.