SEVEN-YEAR WILKINSON MICROWAVE ANISOTROPY PROBE (WMAP^*) OBSERVATIONS: COSMOLOGICAL INTERPRETATION

Electron–photon scattering converts quadrupole temperature anisotropy in the CMB at the decoupling epoch, z = 1090, into linear polarization (Rees 1968; Basko & Polnarev 1980; Kaiser 1983; Bond & Efstathiou 1984; Polnarev 1985; Bond & Efstathiou 1987; Harari & Zaldarriaga 1993). This produces a correlation between the temperature pattern and the polarization pattern (Coulson et al. 1994; Crittenden et al. 1995). Different mechanisms for generating fluctuations produce distinctive correlated patterns in temperature and polarization:

• 1.  
  Adiabatic scalar fluctuations predict a radial polarization pattern around temperature cold spots and a tangential pattern around temperature hot spots on angular scales greater than the horizon size at the decoupling epoch, ≳2°. On angular scales smaller than the sound horizon size at the decoupling epoch, both radial and tangential patterns are formed around both hot and cold spots, as the acoustic oscillation of the CMB modulates the polarization pattern (Coulson et al. 1994). As we have not seen any evidence for non-adiabatic fluctuations (Komatsu et al. 2009a, see Section 4.4 for the seven-year limits), in this section we shall assume that fluctuations are purely adiabatic.
• 2.  
  Tensor fluctuations predict the opposite pattern: the temperature cold spots are surrounded by a tangential polarization pattern, while the hot spots are surrounded by a radial pattern (Crittenden et al. 1995). Since there is no acoustic oscillation for tensor modes, there is no modulation of polarization patterns around temperature spots on small angular scales. We do not expect this contribution to be visible in the WMAP data, given the noise level.
• 3.  
  Defect models predict that there should be minimal correlations between temperature and polarization on 2° ≲ θ ≲ 10° (Seljak et al. 1997). The detection of large-scale temperature polarization fluctuations rules out any causal models as the primary mechanism for generating the CMB fluctuations (Spergel & Zaldarriaga 1997). This implies that the fluctuations were either generated during an accelerating phase in the early universe or were present at the time of the initial singularity.

This section presents the first direct measurement of the predicted pattern of adiabatic scalar fluctuations in CMB polarization maps. We stack maps of Stokes Q and U around temperature hot and cold spots to show the expected polarization pattern at the statistical significance level of 8σ. While we have detected the TE correlations in the first year data (Kogut et al. 2003), we present here the direct real space pattern around hot and cold spots. In Section 2.5, we discuss the relationship between the two measurements.

We first identify temperature hot (or cold) spots, and then stack the polarization data (i.e., Stokes Q and U) on the locations of the spots. As we shall show below, the resulting polarization data are equivalent to the temperature peak–polarization correlation function which is similar to, but different in an important way from, the temperature–polarization correlation function.

2.2.1. Q_r and U_r: Transformed Stokes Parameters

Our definitions of Stokes Q and U follow that of Kogut et al. (2003): the polarization that is parallel to the Galactic meridian is Q>0 and U = 0. Starting from this, the polarization that is rotated by 45° from east to west (clockwise, as seen by an observer on Earth looking up at the sky) has Q = 0 and U>0, that perpendicular to the Galactic meridian has Q < 0 and U = 0, and that rotated further by 45° from east to west has Q = 0 and U < 0. With one more rotation we go back to Q>0 and U = 0. We show this in Figure 1.

Zoom In Zoom Out Reset image size

As the predicted polarization pattern around temperature spots is either radial or tangential, we find it most convenient to work with Q_r and U_r first introduced by Kamionkowski et al. (1997b):

$$\begin{eqnarray} Q_r({\btheta }) &=& -\!Q({\btheta })\cos (2\phi)-U({\btheta })\sin (2\phi), \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} U_r({\btheta }) &=& Q({\btheta })\sin (2\phi)-U({\btheta })\cos (2\phi). \end{eqnarray} \tag{ 2 }$$

These transformed Stokes parameters are defined with respect to the new coordinate system that is rotated by ϕ, and thus they are defined with respect to the line connecting the temperature spot at the center of the coordinate and the polarization at an angular distance θ from the center (also see Figure 1). Note that we have used the small-angle (flat-sky) approximation for simplicity of the algebra. This approximation is justified as we are interested in relatively small angular scales, θ < 5°.

The above definition of Q_r is equivalent to the so-called tangential shear statistic used by the weak gravitational lensing community. By following what has been already done for the tangential shear, we can find the necessary formulae for Q_r and U_r. Specifically, we shall follow the derivations given in Jeong et al. (2009).

With the small-angle approximation, Q and U are related to the E- and B-mode polarization in Fourier space (Seljak & Zaldarriaga 1997; Kamionkowski et al. 1997a) as

$$\begin{eqnarray} -\!Q({\btheta }) &=& \int \frac{d^2{\mathbf l}}{(2\pi)^2} \left[E_{\mathbf l}\cos (2\varphi)-B_{\mathbf l}\sin (2\varphi)\right]e^{i{\mathbf l}\cdot { {\btheta }}},\\ -U({\btheta }) &=& \int \frac{d^2{\mathbf l}}{(2\pi)^2} \left[E_{\mathbf l}\sin (2\varphi)+B_{\mathbf l}\cos (2\varphi)\right]e^{i{\mathbf l}\cdot { {\btheta }}}, \end{eqnarray} \tag{ 3 }$$

where φ is the angle between l and the line of Galactic latitude, l = (lcos φ, lsin φ). Note that we have included the negative signs on the left-hand side because our sign convention for the Stokes parameters is opposite of that used in Equation (38) of Zaldarriaga & Seljak (1997). The transformed Stokes parameters are given by

$$\begin{eqnarray} \nonumber - \! Q_r({\btheta }) &=& -\int \frac{d^2{\mathbf l}}{(2\pi)^2} \left\lbrace E_{\mathbf l}\cos [2(\phi -\varphi)]\right.\\ & & \left.+B_{\mathbf l}\sin [2(\phi -\varphi)]\right\rbrace e^{i{\mathbf l}\cdot { {\btheta }}},\\ \nonumber -U_r({\btheta }) &=& \int \frac{d^2{\mathbf l}}{(2\pi)^2} \left\lbrace E_{\mathbf l}\sin [2(\phi -\varphi)]\right.\\ & & \left.-B_{\mathbf l}\cos [2(\phi -\varphi)]\right\rbrace e^{i{\mathbf l}\cdot { {\btheta }}}. \end{eqnarray} \tag{ 5 }$$

The stacking of Q_r and U_r at the locations of temperature peaks can be written as

$$\begin{eqnarray} &&\langle Q_r\rangle ({\btheta }) = \frac{1}{N_{\rm pk}} \int d^2\hat{\mathbf n}M(\hat{\mathbf n})\langle n_{\rm pk}(\hat{\mathbf {n}})Q_r(\hat{\mathbf {n}}+\hat{ {\btheta }})\rangle,\\ &&\langle U_r\rangle ({\btheta }) = \frac{1}{N_{\rm pk}} \int d^2\hat{\mathbf n}M(\hat{\mathbf n})\langle n_{\rm pk}(\hat{\mathbf {n}})U_r(\hat{\mathbf {n}}+\hat{ {\btheta }})\rangle, \end{eqnarray} \tag{ 7 }$$

where the angle bracket, 〈...〉, denotes the average over the locations of peaks, $n_{\rm pk}(\hat{\mathbf {n}})$ is the surface number density of peaks (of the temperature fluctuation) at the location $\hat{\mathbf {n}}$, N_pk is the total number of temperature peaks used in the stacking analysis, and $M(\hat{\mathbf n})$ is equal to 0 at the masked pixels and 1 otherwise. Defining the density contrast of peaks, $\delta _{\rm pk}\equiv n_{\rm pk}/\bar{n}_{\rm pk}-1$, we find

$$\begin{eqnarray} &&\langle Q_r\rangle ({\btheta }) = \frac{1}{f_{\rm sky}} \int \frac{d^2\hat{\mathbf n}}{4\pi }M(\hat{\mathbf n})\langle \delta _{\rm pk}(\hat{\mathbf {n}})Q_r(\hat{\mathbf {n}}+\hat{ {\btheta }})\rangle, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &&\langle U_r\rangle ({\btheta }) = \frac{1}{f_{\rm sky}} \int \frac{d^2\hat{\mathbf n}}{4\pi }M(\hat{\mathbf n})\langle \delta _{\rm pk}(\hat{\mathbf {n}})U_r(\hat{\mathbf {n}}+\hat{ {\btheta }})\rangle, \end{eqnarray} \tag{ 10 }$$

where $f_{\rm sky}\equiv \int M(\hat{\mathbf n})d^2\hat{\mathbf n}/(4\pi)$ is the fraction of sky outside of the mask, and we have used $N_{\rm pk}=4\pi f_{\rm sky}\bar{n}_{\rm pk}$.

In Appendix B, we use the statistics of peaks of Gaussian random fields to relate 〈Q_r〉 to the temperature–E-mode polarization cross power spectrum C^TE_l, 〈U_r〉 to the temperature–B-mode polarization cross power spectrum C^TB_l, and the stacked temperature profile, 〈T〉, to the temperature power spectrum C^TT_l. We find

$$\begin{eqnarray} \langle Q_r\rangle (\theta) &=& -\int \frac{ldl}{2\pi } W_l^TW_l^P(\bar{b}_\nu +\bar{b}_\zeta l^2)C_l^{\rm TE}J_2(l\theta), \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \langle U_r\rangle (\theta) &=& -\int \frac{ldl}{2\pi } W_l^TW_l^P(\bar{b}_\nu +\bar{b}_\zeta l^2)C_l^{\rm TB}J_2(l\theta), \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \langle T\rangle (\theta) &=& \int \frac{ldl}{2\pi } (W_l^T)^2(\bar{b}_\nu +\bar{b}_\zeta l^2)C_l^{\rm TT}J_0(l\theta), \end{eqnarray} \tag{ 13 }$$

where W^T_l and W^P_l are the harmonic transform of window functions, which are a combination of the experimental beam, pixel window, and any other additional smoothing applied to the temperature and polarization data, respectively, and $\bar{b}_\nu +\bar{b}_\zeta l^2$ is the "scale-dependent bias" of peaks found by Desjacques (2008) averaged over peaks. See Appendix B for details.

2.2.2. Prediction and Physical Interpretation

What do 〈Q_r〉(θ) and 〈U_r〉(θ) look like? The Q_r map is expected to be non-zero for a cosmological signal, while the U_r map is expected to vanish in a parity-conserving universe unless some systematic error rotates the polarization plane uniformly.

To understand the shape of Q_r as well as its physical implications, let us begin by showing the smoothed C^TE_l spectra and the corresponding temperature–Q_r correlation functions, $C^{TQ_r}(\theta)$, in Figure 2. (Note that $C^{TQ_r}$ and $C^{TU_r}$ can be computed from Equations (11) and (12), respectively, with b_ν = 1 and b_ζ = 0.) This shows three distinct effects causing polarization of CMB (see Hu & White 1997, for a pedagogical review):

• 1.  
  θ ≳ 2θ_horizon, where θ_horizon is the angular size of the radius of the horizon size at the decoupling epoch. Using the comoving horizon size of r_horizon = 0.286 Gpc and the comoving angular diameter distance to the decoupling epoch of d_A = 14 Gpc as derived from the WMAP data, we find θ_horizon = 12. As this scale is so much greater than the sound horizon size (see below), only gravity affects the physics. Suppose that there is a Newtonian gravitational potential, Φ_N, at the center of a perturbation, θ = 0. If it is overdense at the center, Φ_N < 0, and thus it is a cold spot according to the Sachs–Wolfe formula (Sachs & Wolfe 1967), ΔT/T = Φ_N/3 < 0. The photon fluid in this region will flow into the gravitational potential well, creating a converging flow. Such a flow creates the quadrupole temperature anisotropy around an electron at θ ⩾ 2θ_horizon, producing polarization that is radial, i.e., Q_r>0. Since the temperature is negative, we obtain 〈TQ_r〉 < 0, i.e., anti-correlation (Coulson et al. 1994). On the other hand, if it is overdense at the center, then the photon fluid moves outward, producing polarization that is tangential, i.e., Q_r < 0. Since the temperature is positive, we obtain 〈TQ_r〉 < 0, i.e., anti-correlation. The anti-correlation at θ ⩾ 2θ_horizon is a smoking gun for the presence of super-horizon fluctuations at the decoupling epoch (Spergel & Zaldarriaga 1997), which has been confirmed by the WMAP data (Peiris et al. 2003).
• 2.  
  θ ≃ 2θ_A, where θ_A is the angular size of the radius of the sound horizon size at the decoupling epoch. Using the comoving sound horizon size of r_s = 0.147 Gpc and d_A = 14 Gpc as derived from the WMAP data, we find θ_A = 06. Again, consider a potential well with Φ_N < 0 at the center. As the photon fluid flows into the well, it compresses, increasing the temperature of the photons. Whether or not this increase can reverse the sign of the temperature fluctuation (from negative to positive) depends on whether the initial perturbation was adiabatic. If it was adiabatic, then the temperature would reverse sign at θ ≲ 2θ_horizon. Note that the photon fluid is still flowing in, and thus the polarization direction is radial, Q_r>0. However, now that the temperature is positive, the correlation reverses sign: 〈TQ_r〉>0. A similar argument (with the opposite sign) can be used to show the same result, 〈TQ_r〉>0, for Φ_N>0 at the center. As an aside, the temperature reverses sign on smaller angular scales for isocurvature fluctuations.
• 3.  
  θ ≃ θ_A. Again, consider a potential well with Φ_N < 0 at the center. At θ ≲ 2θ_A, the pressure of the photon fluid is so great that it can slow down the flow of the fluid. Eventually, at θ ∼ θ_A, the pressure becomes large enough to reverse the direction of the flow (i.e., the photon fluid expands). As a result the polarization direction becomes tangential, Q_r < 0; however, as the temperature is still positive, the correlation reverses sign again: 〈TQ_r〉 < 0.

Zoom In Zoom Out Reset image size

On even smaller scales, the correlation reverses sign again (see Figure 2 of Coulson et al. 1994) because the temperature gets too cold due to expansion. We do not see this effect in Figure 2 because of the smoothing. Lastly, there is no correlation between T and Q_r at θ = 0 because of symmetry.

These features are essentially preserved in the peak– polarization correlation as measured by the stacked polarization profiles. We show them in Figure 3 for various values of the threshold peak heights. The important difference is that, thanks to the scale-dependent bias ∝l², the small-scale trough at θ ≃ θ_A is enhanced, making it easier to observe. On the other hand, the large-scale anti-correlation is suppressed. We can therefore conclude that, with the WMAP data, we should be able to measure the compression phase at θ ≃ 2θ_A = 12, as well as the reversal phase at θ ≃ θ_A = 06. We also show the profiles calculated from numerical simulations (gray solid lines). The agreement with Equation (11) is excellent. We also show the predicted profiles of the stacked temperature data in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.3.1. Temperature Data

We use the foreground-reduced V + W temperature map at the HEALPix resolution of N_side = 512 to find temperature peaks. First, we smooth the foreground-reduced temperature maps in six differencing assemblies (DAs) (V1, V2, W1, W2, W3, W4) to a common resolution of 05 (FWHM) using

$$\begin{equation} \Delta T(\hat{\mathbf n}) = \sum _{lm} a_{lm}\frac{W_l^T}{b_l}Y_{lm}(\hat{\mathbf n}), \end{equation} \tag{ 14 }$$

where b_l is the appropriate beam transfer function for each DA (Jarosik et al. 2011), and W^T_l = p_lexp [ − l(l + 1)σ²_FWHM/(16ln 2)] is the pixel window function for N_side = 512, p_l, times the spherical harmonic transform of a Gaussian with σ_FWHM = 05. We then co-add the foreground-reduced V- and W-band maps with the inverse noise variance weighting, and remove the monopole from the region outside of the mask (which is already negligibly small, 1.07 × 10⁻⁴μK). For the mask, we combine the new seven-year KQ85 mask, KQ85y7 (defined in Gold et al. 2011; also see Section 3.1) and P06 masks, leaving 68.7% of the sky available for the analysis.

We find the locations of minima and maxima using the software " hotspot " in the HEALPix package (Gorski et al. 2005). Over the full sky (without the mask), we find 20953 maxima and 20974 minima. As the maxima and minima found by hotspot still contain negative and positive peaks, respectively, we further select the "hot spots" by removing all negative peaks from maxima, and the "cold spots" by removing all positive peaks from minima. This procedure corresponds to setting the threshold peak height to ν_t = 0; thus, our prediction for 〈Q_r〉(θ) is the top left panel of Figure 3.

Outside of the mask, we find 12,387 hot spots and 12,628 cold spots. The rms temperature fluctuation is σ₀ = 83.9μK. What does the theory predict? Using Equation (B15) with the power spectrum C^TT_l = (C^TT,signal_lp²_l + N^TT_l/b²_l)exp [ − l(l + 1) σ²_FWHM/(8ln 2)] where N^TT_l = 7.47 × 10⁻³μK²sr is the noise bias of the V+W map before Gaussian smoothing and C^TT,signal_l is the five-year best-fitting power-law ΛCDM temperature power spectrum, we find $4\pi f_{\rm sky}\bar{n}_{\rm pk}=12330$ for ν_t = 0 and f_sky = 0.687; thus, the number of observed hot and cold spots is consistent with the predicted number.¹⁵

2.3.2. Polarization Data

In Figures 5 and 6, we show the stacked images of T, Q, U, Q_r, and U_r around temperature cold spots and hot spots, respectively. The peak values of the stacked temperature profiles agree with the predictions (see the dashed line in the top left panel of Figure 4). A dip in temperature (for hot spots; a bump for cold spots) at θ ≃ 1° is clearly visible in the data. While the Stokes Q and U measured from the data exhibit the expected features, they are still fairly noisy. The most striking images are the stacked Q_r (and T). The predicted features are clearly visible, particularly the compression phase at 12 and the reversal phase at 06 in Q_r: the polarization directions around temperature cold spots are radial at θ ≃ 06 and tangential at θ ≃ 12, and those around temperature hot spots show the opposite patterns, as predicted.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

How significant are these features? Before performing the quantitative χ² analysis, we first compare Q_r and U_r using both the (V + W)/2 sum map (here, V + W refers to the inverse noise variance weighted average) as well as the (V − W)/2 difference map (bottom panels of Figures 5 and 6). The Q_r map (which is expected to be non-zero for a cosmological signal) shows clear differences between the sum and difference maps, while the U_r map (which is expected to vanish in a parity-conserving universe unless some systematic error rotates the polarization plane uniformly) is consistent with zero in both the sum and difference maps.

Next, we perform the standard χ² analysis. We summarize the results in Table 5. We report the values of χ² measured with respect to zero signal in the second column, where the number of degrees of freedom (dof) is 625. For each sum map combination, we fit the data to the predicted signal to find the best-fitting amplitude.

The largest improvement in χ² is observed for Q_r, as expected from the visual inspection of Figures 5 and 6: we find 0.82 ±  0.15 and 0.90 ± 0.15 for the stacking of Q_r around hot and cold spots, respectively. The improvement in χ² is Δχ² = −29.2 and −36.2, respectively; thus, we detect the expected polarization patterns around hot and cold spots at the level of 5.4σ and 6σ, respectively. The combined significance exceeds 8σ.

On the other hand, we do not find any evidence for U_r. The χ² values with respect to zero signal per dof are 629.2/625 (hot spots) and 657.8/625 (cold spots), and the probabilities of finding larger values of χ² are 44.5% and 18%, respectively. But, can we learn anything about cosmology from this result? While the standard model predicts C^TB_l = 0 and hence 〈U_r〉 = 0, models in which the global parity symmetry is violated can create C^TB_l = sin(2Δα)C^TE_l (Lue et al. 1999; Carroll 1998; Feng et al. 2005). Therefore, we fit the measured U_r to the predicted Q_r, finding a null result: sin(2Δα) = −0.13 ± 0.15 and 0.20 ± 0.15 (68% CL), or equivalently Δα = −37 ± 43 and 57 ± 43 (68% CL) for hot and cold spots, respectively. Averaging these numbers, we obtain Δα = 10 ± 30 (68% CL), which is consistent with (although not as stringent as) the limit we find from the full analysis presented in Section 4.5. Finally, all the χ² values measured from the difference maps are consistent with a null signal.

How do these results compare to the full analysis of the TE power spectrum? By fitting the seven-year C^TE_l data to the same power spectrum used above (five-year best-fitting power-law ΛCDM model from l = 24 to 800, i.e., dof=777), we find the best-fitting amplitude of 0.999 ± 0.048 and Δχ² = −434.5, i.e., a 21σ detection of the TE signal. This is reasonable, as we used only the V- and W-band data for the stacking analysis, while we used also the Q-band data for measuring the TE power spectrum; 〈Q_r〉(θ) is insensitive to information on θ ≳ 2° (see top left panel of Figure 3); and the smoothing suppresses the power at l ≳ 400 (see left panel of Figure 2). Nevertheless, there is probably a way to extract more information from 〈Q_r〉(θ) by, for example, combining data at different threshold peak heights and smoothing scales.

If the temperature fluctuations of the CMB obey Gaussian statistics and global parity symmetry is respected on cosmological scales, the temperature–E-mode polarization cross power spectrum, C^TE_l, contains all the information about the temperature–polarization correlation. In this sense, the stacked polarization images do not add any new information.

The detection and measurement of the temperature–E mode polarization cross-correlation power spectrum, C^TE_l (Kovac et al. 2002; Kogut et al. 2003; Spergel et al. 2003), can be regarded as equivalent to finding the predicted polarization patterns around hot and cold spots. While we have shown that one can write the stacked polarization profile around temperature spots in terms of an integral of C^TE_l, the formal equivalence between this new method and C^TE_l is valid only when temperature fluctuations obey Gaussian statistics, as the stacked Q and U maps measure correlations between temperature peaks and polarization. So far there is no convincing evidence for non-Gaussianity in the temperature fluctuations observed by WMAP (Komatsu et al. 2003, see Section 6 for the seven-year limits on primordial non-Gaussianity, and Bennett et al. 2011 for discussion on other non-Gaussian features).

Nevertheless, they provide striking confirmation of our understanding of the physics at the decoupling epoch in the form of radial and tangential polarization patterns at two characteristic angular scales that are important for the physics of acoustic oscillation: the compression phase at θ = 2θ_A and the reversal phase at θ = θ_A.

Also, this analysis does not require any analysis in harmonic space, nor decomposition to E and B modes. The analysis is so straightforward and intuitive that the method presented here would also be useful for null tests and systematic error checks. The stacked image of U_r should be particularly useful for systematic error checks.

Any experiments that measure both temperature and polarization should be able to produce the stacked images such as presented in Figures 5 and 6.

Foreground mask. The seven-year temperature analysis masks, KQ85y7 and KQ75y7, have been slightly enlarged to mask the regions that have excess foreground emission, particularly in the H ii regions Gum and Ophiuchus, identified in the difference between foreground-reduced maps at different frequencies. As a result, the new KQ85y7 and KQ75y7 masks eliminate an additional 3.4% and 1.0% of the sky, leaving 78.27% and 70.61% of the sky for the cosmological analyses, respectively. See Section 2.1 of Gold et al. (2011) for details. There is no change in the polarization P06 mask (see Section 4.2 of Page et al. 2007, for definition of this mask), which leaves 73.28% of the sky.

Point sources and the SZ effect. We continue to marginalize over a contribution from unresolved point sources, assuming that the antenna temperature of point sources declines with frequency as ν^−2.09 (see Equation (5) of Nolta et al. 2009). The five-year estimate of the power spectrum from unresolved point sources in Q band in units of antenna temperature, A_ps, was 10³A_ps = 11 ± 1μK²sr (Nolta et al. 2009), and we used this value and the error bar to marginalize over the power spectrum of residual point sources in the seven-year parameter estimation. The subsequent analysis showed that the seven-year estimate of the power spectrum is 10³A_ps = 9.0 ± 0.7μK²sr (Larson et al. 2011), which is somewhat lower than the five-year value because more sources are resolved by WMAP and included in the source mask. The difference in the diffuse mask (between KQ85y5 and KQ85y7) does not affect the value of A_ps very much: we find 9.3 instead of 9.0 if we use the five-year diffuse mask and the seven-year source mask. The source power spectrum is sub-dominant in the total power. We have checked that the parameter results are insensitive to the difference between the five-year and seven-year residual source estimates.

We continue to marginalize over a contribution from the SZ effect using the same template as for the 3- and five-year analyses (Komatsu & Seljak 2002). We assume a uniform prior on the amplitude of this template as 0 < A_SZ < 2, which is now justified by the latest limits from the SPT collaboration, A_SZ = 0.37 ± 0.17 (68% CL; Lueker et al. 2010), and the ACT Collaboration, A_SZ < 1.63 (95% CL; Fowler et al. 2010).

High-l temperature and polarization. We increase the multipole range of the power spectra used for the cosmological parameter estimation from 2–1000 to 2–1200 for the TT power spectrum, and from 2–450 to 2–800 for the TE power spectrum. We use the seven-year V- and W-band maps (Jarosik et al. 2011) to measure the high-l TT power spectrum in l = 33–1200. While we used only Q- and V-band maps to measure the high-l TE and TB power spectra for the five-year analysis (Nolta et al. 2009), we also include W-band maps in the seven-year high-l polarization analysis.

With these data, we now detect the high-l TE power spectrum at 21σ, compared to 13 σ for the five-year high-l TE data. This is a consequence of adding two more years of data and the W-band data. The TB data can be used to probe a rotation angle of the polarization plane, Δα, due to potential parity-violating effects or systematic effects. With the seven-year high-l TB data we find a limit Δα = −09 ± 14 (68% CL). For comparison, the limit from the five-year high-l TB power spectrum was Δα = −12 ± 22 (68% CL; Komatsu et al. 2009a). See Section 4.5 for the seven-year limit on Δα from the full analysis.

Low-l temperature and polarization. Except for using the seven-year maps and the new temperature KQ85y7 mask, there is no change in the analysis of the low-l temperature and polarization data: we use the internal linear combination map (Gold et al. 2011) to measure the low-l TT power spectrum in l = 2–32, and calculate the likelihood using the Gibbs sampling and Blackwell–Rao (BR) estimator (Jewell et al. 2004; Wandelt 2003; Wandelt et al. 2004; O'Dwyer et al. 2004; Eriksen et al. 2004, 2007a, 2007b; Chu et al. 2005; Larson et al. 2007). For the implementation of the BR estimator in the five-year analysis, see Section 2.1 of Dunkley et al. (2009). We use Ka-, Q-, and V-band maps for the low-l polarization analysis in l = 2–23, and evaluate the likelihood directly in pixel space as described in Appendix D of Page et al. (2007).

To get a feel for improvements in the low-l polarization data with two additional years of integration, we note that the seven-year limits on the optical depth, and the tensor-to-scalar ratio and rotation angle from the low-l polarization data alone, are τ = 0.088 ± 0.015 (68% CL; see Larson et al. 2011), r < 1.6 (95% CL; see Section 4.1), and Δα = −38 ± 52 (68% CL; see Section 4.5), respectively. The corresponding five-year limits were τ = 0.087 ± 0.017 (Dunkley et al. 2009), r < 2.7 (see Section 4.1), and Δα = −75 ± 73 (Komatsu et al. 2009a), respectively.

In Table 6, we summarize the improvements from the five-year data mentioned above.

Table 6. Polarization Data: Improvements from the Five-year data

l Range	Type	Seven Year	Five Year
High la	TE	Detected at 21σ	Detected at 13σ
	TB	Δα = −09 ± 14	Δα = −12 ± 22
Low lb	EE	τ = 0.088 ± 0.015	τ = 0.087 ± 0.017
	BB	r < 2.1 (95% CL)	r < 4.7 (95% CL)
	EE/BB	r < 1.6 (95% CL)	r < 2.7 (95% CL)
	TB/EB	Δα = −38 ± 52	Δα = −75 ± 73
All l	TE/EE/BB	r < 0.93 (95% CL)	r < 1.6 (95% CL)
	TB/EBc	Δα = −11 ± 14	Δα = −17 ± 21

Notes. ^al ⩾ 24. The Q-, V-, and W-band data are used for the seven-year analysis, whereas only the Q- and V-band data were used for the five-year analysis. ^b2 ⩽ l ⩽ 23. The Ka-, Q-, and V-band data are used for both the seven-year and five-year analyses. ^cThe quoted errors are statistical only and do not include the systematic error ±15 (see Section 4.5).

Download table as:  ASCIITypeset image

The WMAP data are statistically powerful enough to constrain six parameters of a flat ΛCDM model with a tilted spectrum. However, to constrain deviations from this minimal model, other CMB data probing smaller angular scales and astrophysical data probing the expansion rates, distances, and growth of structure are useful.

3.2.1. Small-scale CMB Data

The best limits on the primordial helium abundance, Y_p, are obtained when the WMAP data are combined with the power spectrum data from other CMB experiments probing smaller angular scales, l ≳ 1000.

We use the temperature power spectra from the Arcminute Cosmology Bolometer Array Receiver (ACBAR; Reichardt et al. 2009) and QUEST at DASI (QUaD) (Brown et al. 2009) experiments. For the former, we use the temperature power spectrum binned in 16 band powers in the multipole range 900 < l < 2000. For the latter, we use the temperature power spectrum binned in 13 band powers in 900 < l < 2000.

We marginalize over the beam and calibration errors of each experiment: for ACBAR, the beam error is 2.6% on a 5 arcmin (FWHM) Gaussian beam and the calibration error is 2.05% in temperature. For QUaD, the beam error combines a 2.5% error on 5.2 and 3.8 arcmin (FWHM) Gaussian beams at 100 GHz and 150 GHz, respectively, with an additional term accounting for the sidelobe uncertainty (see Appendix A of Brown et al. 2009, for details). The calibration error is 3.4% in temperature.

The ACBAR data are calibrated to the WMAP five-year temperature data, and the QUaD data are calibrated to the BOOMERanG data (Masi et al. 2006) which are, in turn, calibrated to the WMAP 1-year temperature data. (The QUaD team takes into account the change in the calibration from the 1-year to the five-year WMAP data.) The calibration errors quoted above are much greater than the calibration uncertainty of the WMAP five-year data (0.2%; Hinshaw et al. 2007). This is due to the noise of the ACBAR, QUaD, and BOOMERanG data. In other words, the above calibration errors are dominated by the statistical errors that are uncorrelated with the WMAP data. We thus treat the WMAP, ACBAR, and QUaD data as independent.

Figure 7 shows the WMAP seven-year temperature power spectrum (Larson et al. 2011) as well as the temperature power spectra from ACBAR and QUaD.

Zoom In Zoom Out Reset image size

We do not use the other, previous small-scale CMB data, as their statistical errors are much larger than those of ACBAR and QUaD, and thus adding them would not improve the constraints on the cosmological parameters significantly. The new power-spectrum data from the SPT (Lueker et al. 2010) and ACT (Fowler et al. 2010) Collaborations were not yet available at the time of our analysis.

3.2.2. Hubble Constant and Angular Diameter Distances

3.2.3. Power Spectrum of Luminous Red Galaxies

3.2.4. Luminosity Distances

3.2.5. Time-delay Distance

Can we measure angular diameter distances out to higher redshifts? Measurements of gravitational lensing time delays offer a way to determine absolute distance scales (Refsdal 1964). When a foreground galaxy lenses a background variable source (e.g., quasars) and produces multiple images of the source, changes of the source luminosity due to variability appear on multiple images at different times.

The time delay at a given image position ${\btheta }$ for a given source position ${\bbeta }$, $t({\btheta },{\bbeta })$, depends on the angular diameter distances as (see, e.g., Schneider et al. 2006, for a review)

$$\begin{equation} t({\btheta },{\bbeta }) = \frac{1+z_{\rm l}}{c}\frac{D_{\rm l}D_{\rm s}}{D_{\rm ls}}\phi _{\rm F}({\btheta },{\bbeta }), \end{equation} \tag{ 18 }$$

where D_l, D_s, and D_ls are the angular diameter distances out to a lens galaxy, to a source galaxy, and between them, respectively, and ϕ_F is the so-called Fermat potential, which depends on the path length of light rays and gravitational potential of the lens galaxy.

The biggest challenge for this method is to control systematic errors in our knowledge of ϕ_F, which requires a detailed modeling of mass distribution of the lens. One can, in principle, minimize this systematic error by finding a lens system where the mass distribution of lens is relatively simple.

The lens system B1608+656 is not a simple system, with two lens galaxies and dust extinction; however, it has one of the most precise time-delay measurements of quadruple lenses. The lens redshift of this system is relatively large, z_l = 0.6304 (Myers et al. 1995). The source redshift is z_s = 1.394 (Fassnacht et al. 1996). This system has been used to determine H₀ to 10% accuracy (Koopmans et al. 2003).

Suyu et al. (2009) have obtained more data from the deep HST Advanced Camera for Surveys (ACS) observations of the asymmetric and spatially extended lensed images, and constrained the slope of mass distribution of the lens galaxies. They also obtained ancillary data (for stellar dynamics and lens environment studies) to control the systematics, particularly the so-called "mass-sheet degeneracy," which the strong lensing data alone cannot break. By doing so, they were able to reduce the error in H₀ (including the systematic error) by a factor of two (Suyu et al. 2010). They find a constraint on the "time-delay distance," D_Δt, as

$$\begin{equation} D_{\Delta t} \equiv (1+z_{\rm l}) \frac{D_{\rm l}D_{\rm s}}{D_{\rm ls}} \simeq 5226\pm 206 \,{\rm Mpc}, \end{equation} \tag{ 19 }$$

where the number is found from a Gaussian fit to the likelihood of D_Δt¹⁶; however, the actual shape of the likelihood is slightly non-Gaussian. We thus use

• 1.  
  Likelihood of D_Δt out to the lens system B1608+656 given by Suyu et al. (2010),

  $$\begin{equation} P(D_{\Delta t}) = \frac{\exp [-(\ln (x-\lambda)-\mu)^2/(2\,\sigma ^2)]}{\sqrt{2\pi }(x-\lambda)\,\sigma }, \end{equation} \tag{ 20 }$$

  where x = D_Δt/(1 Mpc), λ = 4000, μ = 7.053, and  σ = 0.2282. This likelihood includes systematic errors due to the mass-sheet degeneracy, which dominates the total error budget (see Section 6 of Suyu et al. 2010, for more details). Note that this is the only lens system for which D_Δt (rather than H₀) has been constrained.¹⁷

When we evaluate the likelihood of external astrophysical data sets, we often need to compute the Hubble expansion rate, H(a). While we treated the effect of massive neutrinos on H(a) approximately for the five-year analysis of the external data sets presented in Komatsu et al. (2009a), we treat it exactly for the seven-year analysis, as described below.

The total energy density of massive neutrino species, ρ_ν, is given by (in natural units)

$$\begin{equation} \rho _\nu (a) = 2\int \frac{d^3p}{(2\pi)^3}\frac{1}{e^{p/T_\nu (a)}+1} \sum _{i} \sqrt{p^2+m^2_{\nu,i}}, \end{equation} \tag{ 21 }$$

where m_ν,i is the mass of each neutrino species. Using the comoving momentum, q ≡ pa, and the present-day neutrino temperature, T_ν0 = (4/11)^1/3T_cmb = 1.945 K, we write

$$\begin{equation} \rho _\nu (a) = \frac{1}{a^4}\int \frac{q^2dq}{\pi ^2}\frac{1}{e^{q/T_{\nu 0}}+1} \sum _{i} \sqrt{q^2+m^2_{\nu,i}a^2}. \end{equation} \tag{ 22 }$$

Throughout this paper, we shall assume that all massive neutrino species have the equal mass m_ν, i.e., m_ν,i = m_ν for all i.¹⁸

When neutrinos are relativistic, one may relate ρ_ν to the photon energy density, ρ_γ, as

$$\begin{equation} \rho _\nu (a)\rightarrow \frac{7}{8} \left(\frac{4}{11}\right)^{4/3} N_{\rm eff}\rho _\gamma (a) \simeq 0.2271 N_{\rm eff}\rho _\gamma (a), \end{equation} \tag{ 23 }$$

where N_eff is the effective number of neutrino species. Note that N_eff = 3.04 for the standard neutrino species.¹⁹ This motivates our writing (Equation (22)) as

$$\begin{equation} \rho _\nu (a) = 0.2271{N_{\rm eff}}\rho _\gamma (a) f(m_{\nu }a/T_{\nu 0}), \end{equation} \tag{ 24 }$$

where

$$\begin{equation} f(y)\equiv \frac{120}{7\pi ^4} \int _0^\infty dx\frac{x^2\sqrt{x^2+y^2}}{e^x+1}. \end{equation} \tag{ 25 }$$

The limits of this function are f(y) → 1 for y → 0, and $f(y)\rightarrow \frac{180\zeta (3)}{7\pi ^4}y$ for y → ∞, where ζ(3) ≃ 1.202 is the Riemann zeta function. We find that f(y) can be approximated by the following fitting formula:²⁰

$$\begin{equation} f(y)\approx [1+(Ay)^p]^{1/p}, \end{equation} \tag{ 26 }$$

where $A=\frac{180\zeta (3)}{7\pi ^4}\simeq 0.3173$ and p = 1.83. This fitting formula is constructed such that it reproduces the asymptotic limits in y → 0 and y → ∞ exactly. This fitting formula underestimates f(y) by 0.1% at y ≃  2.5 and overestimates by 0.35% at y ≃ 10. The errors are smaller than these values at other y's.

Using this result, we write the Hubble expansion rate as

$$\begin{eqnarray} \nonumber H(a) &=& H_0\left\lbrace \frac{\Omega _c+\Omega _b}{a^3}+\frac{\Omega _\gamma }{a^4} \left[1+0.2271N_{\rm eff} f(m_{\nu }a/T_{\nu 0})\right]\right.\\ & &\left.+ \, \frac{\Omega _k}{a^2} + \frac{\Omega _\Lambda }{a^{3(1+w_{\rm eff}(a))}} \right\rbrace ^{1/2}, \end{eqnarray} \tag{ 27 }$$

where Ω_γ = 2.469 × 10⁻⁵h⁻² for T_cmb = 2.725 K. Using the massive neutrino density parameter, Ω_νh² = ∑m_ν/(94eV), for the standard three neutrino species, we find

$$\begin{equation} \frac{m_\nu a}{T_{\nu 0}} = \frac{187}{1+z}\left(\frac{\Omega _\nu h^2}{10^{-3}}\right). \end{equation} \tag{ 28 }$$

One can check that (Ω_γ/a⁴)0.2271N_efff(m_νa/T_ν0) → Ω_ν/a³ for a → ∞. One may compare Equation (27), which is exact (if we compute f(y) exactly), to Equation (7) of Komatsu et al. (2009a), which is approximate.

Throughout this paper, we shall use Ω_Λ to denote the dark energy density parameter at present: Ω_Λ ≡ Ω_de(z = 0). The function w_eff(a) in Equation (28) is the effective equation of state of dark energy given by $w_{\rm eff}(a)\equiv \frac{1}{\ln a}\int _0^{\ln a} d\ln a^{\prime } w(a^{\prime })$, and w(a) is the usual dark energy equation of state, i.e., the dark energy pressure divided by the dark energy density: w(a) ≡ P_de(a)/ρ_de(a). For vacuum energy (cosmological constant), w does not depend on time, and w = −1.

The seven-year WMAP data combined with BAO and H₀ exclude the scale-invariant spectrum by 99.5% CL, if we ignore tensor modes (gravitational waves).

For a power-law spectrum of primordial curvature perturbations ${\cal R}_k$, i.e.,

$$\begin{equation} \Delta ^2_{\cal R}(k) = \frac{k^3\langle |{\cal R}_k|^2\rangle }{2\pi ^2} = \Delta ^2_{\cal R}(k_0)\left(\frac{k}{k_0}\right)^{n_s-1}, \end{equation} \tag{ 29 }$$

where k₀ = 0.002 Mpc⁻¹, we find

$$\begin{equation*} n_s=0.968\pm 0.012 (68\% {\rm CL}). \end{equation*}$$

For comparison, the WMAP data-only limit is n_s = 0.967 ± 0.014 (Larson et al. 2011), and the WMAP plus the small-scale CMB experiments ACBAR (Reichardt et al. 2009) and QUaD (Brown et al. 2009) is n_s = 0.966^+0.014_−0.013. As explained in Section 3.1.2 of Komatsu et al. (2009a), the small-scale CMB data do not reduce the error bar in n_s very much because of relatively large statistical errors, beam errors, and calibration errors.

How about tensor modes? While the B-mode polarization is a smoking gun for tensor modes (Seljak & Zaldarriaga 1997; Kamionkowski et al. 1997b), the WMAP data mainly constrain the amplitude of tensor modes by the low-l temperature power spectrum (see Section 3.2.3 of Komatsu et al. 2009a). Nevertheless, it is still useful to see how much constraint one can obtain from the seven-year polarization data.

We first fix the cosmological parameters at the five-year WMAP best-fit values of a power-law ΛCDM model. We then calculate the tensor mode contributions to the B-mode, E-mode, and TE power spectra as a function of one parameter: the amplitude, in the form of the tensor-to-scalar ratio, r, defined as

$$\begin{equation} r \equiv \frac{\Delta _h^2(k_0)}{\Delta _{\cal R}^2(k_0)}, \end{equation} \tag{ 30 }$$

where Δ²_h(k) is the power spectrum of tensor metric perturbations, h_k, given by

$$\begin{equation} \Delta _h^2(k) = \frac{4k^3\langle |h_k|^2\rangle }{2\pi ^2} = \Delta ^2_{h}(k_0)\left(\frac{k}{k_0}\right)^{n_t}. \end{equation} \tag{ 31 }$$

In Figure 8, we show the limits on r from the B-mode power spectrum only (r < 2.1, 95% CL), from the B- and E-mode power spectra combined (r < 1.6), and from the B-mode, E-mode, and TE power spectra combined (r < 0.93). These limits are significantly better than those from the five-year data (r < 4.7, 2.7, and 1.6, respectively), because of the smaller noise and shifts in the best-fitting values. For comparison, the B-mode power spectrum from the BICEP 2-year data gives r < 0.73 (95% CL; Chiang et al. 2010).

Zoom In Zoom Out Reset image size

If we add the temperature power spectrum, but still fix all the other cosmological parameters including n_s, then we find r < 0.15 (95% CL) from both five-year and seven-year data; however, due to a strong correlation between n_s and r, this would be an underestimate of the error. For a 7-parameter model (a flat ΛCDM model with a tilted spectrum, tensor modes, and n_t = −r/8), we find r < 0.36(95%CL) from the WMAP data alone (Larson et al. 2011), r < 0.33(95%CL) from WMAP plus ACBAR and QUaD,

$$$ r < 0.24\ \mbox{(95\% CL)} $$$

from WMAP+BAO+H₀, and r < 0.20(95%CL) from WMAP+BAO+SN, where "SN" is the Constitution samples compiled by Hicken et al. (2009a; see Section 3.2.4).

We give a summary of these numbers in Table 7.

Let us relax the assumption that the power spectrum is a pure power law, and add a "running index," dn_s/dln k as (Kosowsky & Turner 1995)

$$\begin{equation} \Delta ^2_{\cal R}(k) = \Delta ^2_{\cal R}(k_0)\left(\frac{k}{k_0}\right)^{n_s(k_0)-1+\frac{1}{2}\ln (k/k_0)dn_s/d\ln k}. \end{equation} \tag{ 32 }$$

Ignoring tensor modes again, we find

$$$ dn_s/d\ln {k} = -0.022\pm 0.020 (68\% {\rm CL}), $$$

from WMAP+BAO+H₀. For comparison, the WMAP data-only limit is dn_s/dln k = −0.034 ± 0.026 (Larson et al. 2011), and the WMAP+ACBAR+QUaD limit is dn_s/dln k = −0.041^+0.022_−0.023.

None of these data combinations require dn_s/dln k: improvements in a goodness-of-fit relative to a power-law model (Equation (29)) are Δχ² = −2ln(L_running/L_power-law) = −1.2, −2.6, and −0.72 for the WMAP-only, WMAP+ACBAR+QUaD, and WMAP+BAO+H₀, respectively. See Table 7 for the case where both r and dn_s/dln k are allowed to vary.

A simple power-law primordial power spectrum without tensor modes continues to be an excellent fit to the data. While we have not done a non-parametric study of the shape of the power spectrum, recent studies after the five-year data release continue to show that there is no convincing deviation from a simple power-law spectrum (Peiris & Verde 2010; Ichiki et al. 2010; Hamann et al. 2010).

While the WMAP data alone cannot constrain the spatial curvature parameter of the observable universe, Ω_k, very well, combining the WMAP data with other distance indicators such as H₀, BAO, or supernovae can constrain Ω_k (e.g., Spergel et al. 2007).

Assuming a ΛCDM model (w = −1), we find

$$$ -0.0133<\Omega _k<0.0084\ \mbox{(95\% CL)}, $$$

from WMAP+BAO+H₀.²¹ However, the limit weakens significantly if dark energy is allowed to be dynamical, w ≠ −1, as this data combination, WMAP+BAO+H₀, cannot constrain w very well. We need additional information from Type Ia supernovae to constrain w and Ω_k simultaneously (see Section 5.3 of Komatsu et al. 2009a). We shall explore this possibility in Section 5.

Non-adiabatic fluctuations are a powerful probe of the origin of matter and the physics of inflation. Following Section 3.6 of Komatsu et al. (2009a), we focus on two physically motivated models for non-adiabatic fluctuations: axion type (Seckel & Turner 1985; Linde 1985, 1991; Turner & Wilczek 1991) and curvaton type (Linde & Mukhanov 1997; Lyth & Wands 2003; Moroi & Takahashi 2001, 2002; Bartolo & Liddle 2002).

For both cases, we consider non-adiabatic fluctuations between photons and cold dark matter:

$$\begin{equation} {\cal S}=\frac{\delta \rho _c}{\rho _c}-\frac{3\delta \rho _\gamma }{4\rho _\gamma }, \end{equation} \tag{ 33 }$$

and report limits on the ratio of the power spectrum of ${\cal S}$ and that of the curvature perturbation ${\cal R}$ (e.g., Bean et al. 2006):

$$\begin{equation} \frac{\alpha (k_0)}{1-\alpha (k_0)} = \frac{P_{\cal S}(k_0)}{P_{\cal R}(k_0)}, \end{equation} \tag{ 34 }$$

where k₀ = 0.002 Mpc⁻¹. We denote the limits on axion-type and curvaton-type by α₀ and α₋₁, respectively.²²

We find no evidence for non-adiabatic fluctuations. The WMAP data-only limits are α₀ < 0.13(95%CL) and α₋₁ < 0.011 (95% CL; Larson et al. 2011). With WMAP+BAO+H₀, we find

$$$ \alpha _{0} < 0.077\ \mbox{(95\% CL)} \mbox{ and } \alpha _{-1} < 0.0047\ \mbox{(95\% CL)}, $$$

while with WMAP+BAO+SN, we find α₀ < 0.064(95%CL) and α₋₁ < 0.0037(95%CL).

The limit on α₀ has an important implication for axion dark matter. In particular, a limit on α₀ is related to a limit on the tensor-to-scalar ratio, r (Kain 2006; Beltran et al. 2007; Sikivie 2008; Kawasaki & Sekiguchi 2008). The explicit formula is given by Equation (48) of Komatsu et al. (2009a) as²³

$$\begin{eqnarray} r &=& \frac{4.7\times 10^{-12}}{\theta _a^{10/7}}\left(\frac{\Omega _ch^2}{\gamma }\right)^{12/7} \left(\frac{\Omega _c}{\Omega _a}\right)^{2/7} \frac{\alpha _{0}}{1-\alpha _{0}}, \end{eqnarray} \tag{ 36 }$$

where Ω_a ⩽ Ω_c is the axion density parameter, θ_a is the phase of the Peccei–Quinn field within our observable universe, and γ ⩽ 1 is a "dilution factor" representing the amount by which the axion density parameter, Ω_ah², would have been diluted due to a potential late-time entropy production by, e.g., decay of some (unspecified) heavy particles, between 200 MeV and the epoch of nucleosynthesis, 1 MeV.

Where does this formula come from? Within the context of the "misalignment" scenario of axion dark matter,²⁴ there are two observables one can use to place limits on the axion properties: the dark matter density and α₀. They are given by (e.g., Kawasaki & Sekiguchi 2008, and references therein)

$$\begin{eqnarray} \frac{\alpha _0(k)}{1-\alpha _0(k)} &=& \frac{\Omega ^2_a}{\Omega ^2_c}\frac{8\epsilon }{\theta ^2_a(f_a/M_{\rm pl})^2}, \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} \Omega _ah^2 &=& 1.0\times 10^{-3}\gamma \theta _a^2\left(\frac{f_a}{10^{10}\, {\rm GeV}}\right)^{7/6}, \end{eqnarray} \tag{ 38 }$$

where f_a is the axion decay constant, and $\epsilon =-\dot{H}_{\rm inf}/H^2_{\rm inf}$ is the so-called slow-roll parameter (where H_inf is the Hubble expansion rate during inflation). For single-field inflation models, is related to r as r = 16. By eliminating the axion decay constant, one obtains Equation (36).

In deriving the above formula for Ω_ah² (Equation (38)), we have assumed that the axion field began to oscillate before the QCD phase transition.²⁵ This is true when $f_a<{\cal O}(10^{-2})M_{\rm {pl}}$; however, when $f_a>{\cal O}(10^{-2})M_{\rm {pl}}$, the axions are so light that the axion field would not start oscillating after the QCD phase transition.²⁶ In this limit, the formula for Ω_ah² is given by

$$\begin{equation} \Omega _ah^2=1.6\times 10^5\gamma \theta _a^2\left(\frac{f_a}{10^{17} {\rm GeV}}\right)^{3/2}. \end{equation} \tag{ 39 }$$

By eliminating f_a from Equations (37) and (39), we obtain another formula for r:

$$\begin{equation} r=\frac{4.0\times 10^{-10}}{\theta _a^{2/3}} \left(\frac{\Omega _ch^2}{\gamma }\right)^{4/3} \left(\frac{\Omega _c}{\Omega _a}\right)^{2/3} \frac{\alpha _{0}}{1-\alpha _{0}}. \end{equation} \tag{ 40 }$$

Equations (36) and (40), combined with our limits on Ω_ch² and α₀, imply that the axion dark matter scenario in which axions account for most of the observed amount of dark matter, Ω_a ∼ Ω_c, must satisfy

$$\begin{eqnarray} r&<&\frac{7.6\times 10^{-15}}{\theta _a^{10/7}\gamma ^{12/7}} \quad \mbox{for $f_a<{\cal O}(10^{-2})M_{\rm {pl}}$},\\ r&<&\frac{1.5\times 10^{-12}}{\theta _a^{2/3}\gamma ^{4/3}} \quad \mbox{for $f_a>{\cal O}(10^{-2})M_{\rm {pl}}$}. \end{eqnarray} \tag{ 41 }$$

Alternatively, one can express this constraint as

$$\begin{eqnarray} \nonumber \theta _a\gamma ^{6/5}&<&3.3\times 10^{-9}\left(\frac{10^{-2}}{r}\right)^{7/10} \mbox{for $f_a<{\cal O}(10^{-2})M_{\rm {pl}}$},\\ \nonumber \theta _a\gamma ^{2}&<&1.8\times 10^{-15}\left(\frac{10^{-2}}{r}\right)^{3/2} \mbox{for $f_a>{\cal O}(10^{-2})M_{\rm {pl}}$}. \end{eqnarray}$$

Therefore, a future detection of tensor modes at the level of r = 10⁻² would imply a fine-tuning of θ_a or γ or both of these parameters (Komatsu et al. 2009a). If such fine-tunings are not permitted, axions cannot account for the observed abundance of dark matter (in the misalignment scenario that we have considered here).

Depending on one's interest, one may wish to eliminate the phase, leaving the axion decay constant in the formula (see Equation (B7) of Komatsu et al. 2009a):

$$\begin{equation} r = (1.6\times 10^{-12})\left(\frac{\Omega _ch^2}{\gamma }\right) \left(\frac{\Omega _c}{\Omega _a}\right)\left(\frac{f_a}{10^{12} {\rm GeV}}\right)^{5/6} \frac{\alpha _0}{1-\alpha _0}, \end{equation} \tag{ 43 }$$

for $f<{\cal O}(10^{-2})M_{\rm {pl}}$. This formula gives

$$\begin{equation} f_a>1.8\times 10^{26} {\rm GeV} \gamma ^{6/5}\left(\frac{r}{10^{-2}}\right)^{6/5}, \end{equation} \tag{ 44 }$$

which is inconsistent with the condition $f_a<{\cal O}(10^{-2})M_{\rm {pl}}$ (unless r is extremely small). The formula that is valid for $f>{\cal O}(10^{-2})M_{\rm {pl}}$ is

$$\begin{equation} r= (2.2\times 10^{-8})\left(\frac{\Omega _ch^2}{\gamma }\right) \left(\frac{\Omega _c}{\Omega _a}\right)\left(\frac{f_a}{10^{17} {\rm GeV}}\right)^{1/2} \frac{\alpha _0}{1-\alpha _0}, \end{equation} \tag{ 45 }$$

which gives

$$\begin{equation} f_a>3.2\times 10^{32} {\rm GeV} \gamma ^{2}\left(\frac{r}{10^{-2}}\right)^{2}. \end{equation} \tag{ 46 }$$

Requiring $f_a<M_{\rm {pl}}=2.4\times 10^{18}$ GeV, we obtain

$$\begin{equation} r < \frac{8.7\times 10^{-10}}{\gamma }. \end{equation} \tag{ 47 }$$

Thus, a future detection of tensor modes at the level of r = 10⁻² implies a significant amount of entropy production, γ ≪ 1, or a super-Planckian axion decay constant, $f_a\gg M_{\rm {pl}}$, or both. Also see Hertzberg et al. (2008), Mack (2009), and Mack & Steinhardt (2009) for similar studies.

For the implications of α₋₁ for curvaton dark matter, see Section 3.6.4 of Komatsu et al. (2009a).

While the TB and EB correlations vanish in a parity-conserving universe, they may not vanish when global parity symmetry is broken on cosmological scales (Lue et al. 1999; Carroll 1998). In pixel space, they would show up as a non-vanishing 〈U_r〉. As we showed already in Section 2.4, the WMAP seven-year 〈U_r〉 data are consistent with noise. What can we learn from this?

It is now a routine work of CMB experiments to deliver the TB and EB data, and constrain a rotation angle of the polarization plane due to a parity-violating effect (or a rotation due to some systematic error). Specifically, a rotation of the polarization plane by an angle Δα gives the following five transformations:

$$\begin{eqnarray} C_l^{\rm TE,obs} &=& C_l^{\rm TE}\cos (2\Delta \alpha),\\ C_l^{\rm TB,obs} &=& C_l^{\rm TE}\sin (2\Delta \alpha),\\ C_l^{\rm EE,obs} &=& C_l^{\rm EE}\cos ^2(2\Delta \alpha),\\ C_l^{\rm BB,obs} &=& C_l^{\rm EE}\sin ^2(2\Delta \alpha),\\ C_l^{\rm EB,obs} &=& \frac{1}{2}C_l^{\rm EE}\sin (4\Delta \alpha), \end{eqnarray} \tag{ 48 }$$

where C_l's on the right-hand side are the primordial power spectra in the absence of rotation, while C^obs_l's on the left-hand side are what we would observe in the presence of rotation.

Note that these equations are not exact but valid only when the primordial B-mode polarization is negligible compared to the E-mode polarization, i.e., $C_l^{\scriptsize\textit{BB}}\ll C_l^{\scriptsize\textit{EE}}$. For the full expression including $C_l^{\scriptsize\textit{BB}}$, see Lue et al. (1999) and Feng et al. (2005).

Roughly speaking, when the polarization data are still dominated by noise rather than by the cosmic signal, the uncertainty in Δα is given by a half of the inverse of the signal-to-noise ratio of TE or EE, i.e.,

$$\begin{eqnarray} {\rm Err}[\Delta \alpha ^{\rm TB}] \nonumber &\simeq & \frac{1}{2({S/N})^{\rm TE}},\\ \nonumber {\rm Err}[\Delta \alpha ^{\rm EB}] &\simeq & \frac{1}{2({S/N})^{\rm EE}}. \end{eqnarray}$$

(Note that we use the full likelihood code to find the best-fitting values and error bars. These equations should only be used to provide an intuitive feel of how the errors scale with signal-to-noise.) As we mentioned in the last paragraph of Section 2.4, with the seven-year polarization data we detect the TE power spectrum at 21σ from l = 24 to 800. We thus expect Err[Δα^TB] ≃ 1/42 ≃ 0.024rad ≃ 14, which is significantly better than the five-year value, 22 (Komatsu et al. 2009a). On the other hand, we detect the EE power spectrum at l ⩾ 24 only at a few σ level, and thus Err[Δα^EB] ≫ Err[Δα^TB], implying that we may ignore the high-l EB data.

The magnitude of polarization rotation angle, Δα, depends on the path length over which photons experienced a parity-violating interaction. As pointed out by Liu et al. (2006), this leads to the polarization angle that depends on l. We can divide this l-dependence in two regimes: (1) l ≲ 20: the polarization signal was generated during reionization (Zaldarriaga 1997). We are sensitive only to the polarization rotation between the reionization epoch and present epoch. (2) l ≳ 20: the polarization signal was generated at the decoupling epoch. We are sensitive to the polarization rotation between the decoupling epoch and present epoch; thus, we have the largest path length in this case.

Using the high-l TB data from l = 24 to 800, we find Δα = −09 ± 14, which is a significant improvement over the five-year high-l result, Δα = −12 ± 22 (Komatsu et al. 2009a).

Let us turn our attention to lower multipoles, l ⩽ 23. Here, with the seven-year polarization data, the EE power spectrum is detected at 5.1σ, whereas the TE power spectrum is only marginally seen (1.9σ). (The overall significance level of detection of the E-model polarization at l ⩽ 23, including EE and TE, is 5.5σ.) We therefore use both the TB and EB data at l ⩽ 23. We find Δα = −38 ± 52, which is also a good improvement over the five-year low-l value, Δα = −75 ± 73.

Combining the low-l TB/EB and high-l TB data, we find Δα = −11 ± 14 (the five-year combined limit was Δα = −17 ± 21), where the quoted error is purely statistical; however, the WMAP instrument can measure the polarization angle to within ±15 of the design orientation (Page et al. 2003, 2007). We thus add 15 as an estimate of a potential systematic error. Our final seven-year limit is

$$$ \Delta \alpha = -1\mbox{$.\!\!^\circ $}1 \pm 1\mbox{$.\!\!^\circ $}4 \mbox{({\rm stat.})} \pm 1\mbox{$.\!\!^\circ $}5 \mbox{({\rm syst.})} \mbox{(68\% CL)}, $$$

or −50 < Δα < 28 (95% CL), for which we have added the statistical and systematic errors in quadrature (which may be an underestimate of the total error). The statistical error and systematic error are now comparable.

Several research groups have obtained limits on Δα from various data sets (Feng et al. 2006; Kostelecký & Mewes 2007; Cabella et al. 2007; Xia et al. 2008a; Xia et al. 2008c; Wu et al. 2009; Gubitosi et al. 2009). Recently, the BOOMERanG Collaboration (Pagano et al. 2009) revisited a limit on Δα from their 2003 flight (B2K), taking into account the effect of systematic errors rotating the polarization angle by −09 ± 07. By removing this, they find Δα = −43 ± 41 (68% CL). The QUaD Collaboration used their final data set to find Δα = 064 ± 050(stat.) ± 050(syst.) (68% CL; Brown et al. 2009). Xia et al. (2010) used the BICEP 2-year data (Chiang et al. 2010) to find Δα = −26 ± 10 (68% CL statistical); however, a systematic error of ±07 needs to be added to this error budget (see "Polarization orientation uncertainty" in Table 3 of Takahashi et al. 2010). Therefore, basically the systematic errors in recent measurements of Δα from WMAP seven-year, QUaD final, and BICEP 2-year data are comparable to the statistical errors.

Adding the statistical and systematic errors in quadrature and averaging over WMAP, QUaD and BICEP with the inverse variance weighting, we find Δα = −025 ± 058 (68% CL), or −141 < Δα < 091 (95% CL). We therefore conclude that the microwave background data are comfortably consistent with a parity-conserving universe. See, e.g., Kostelecký & Mewes (2008), Arvanitaki et al. (2010), and references therein for implications of this result for potential violations of Lorentz invariance and CPT symmetry.

Following Section 6.1 of Komatsu et al. (2009a; also see references therein), we constrain the total mass of neutrinos, ∑m_ν = 94eV(Ω_νh²), mainly from the seven-year WMAP data combined with the distance information. A new component in the analysis is the exact treatment of massive neutrinos when calculating the likelihood of the BAO data, as described in Section 3.3 (also see Wright 2006).

For a flat ΛCDM model, i.e., w = −1 and Ω_k = 0, the WMAP-only limit is ∑m_ν < 1.3 eV(95%CL), while the WMAP+BAO+H₀ limit is

$$$ \sum m_\nu < 0.58\ \mbox{eV}\ \mbox{(95\% CL)} \quad (\mbox{for } w=-1). $$$

The latter is the best upper limit on ∑m_ν without information on the growth of structure, which is achieved by a better measurement of the early Integrated Sachs–Wolfe (ISW) effect through the third acoustic peak of the seven-year temperature power spectrum (Larson et al. 2011), as well as by a better determination of H₀ from Riess et al. (2009). For explanations of this effect, see Ichikawa et al. (2005) or Section 6.1.3 of Komatsu et al. (2009a).

Sekiguchi et al. (2010) combined the five-year version of WMAP+BAO+H₀ with the small-scale CMB data to find ∑m_ν < 0.66 eV (95% CL). Therefore, the improvement from this value to our seven-year limit, ∑m_ν < 0.58eV, indeed comes from a better determination of the amplitude of the third acoustic peak in the seven-year temperature data.

The limit improves when information on the growth of structure is added. For example, with WMAP+H₀ and the power spectrum of LRGs (Reid et al. 2010b; see Section 3.2.3) combined, we find ∑m_ν < 0.44 eV(95%CL) for w = −1.

The WMAP+BAO+H₀ limit on the neutrino mass weakens significantly to ∑m_ν < 1.3 eV(95%CL) for w ≠ −1 because we do not use information of Type Ia supernovae here to constrain w. This is driven by w being too negative: there is an anti-correlation between w and ∑m_ν (Hannestad 2005). The best-fitting value of w in this case is w = −1.44 ± 0.27 (68% CL).²⁷ For WMAP+LRG+H₀, we find ∑m_ν < 0.71eV (95% CL) for w ≠ −1. When the Constitution supernova data are included (WMAP+BAO+SN), we find ∑m_ν < 0.71²⁸ and 0.91 eV (95% CL) for w = −1 and w ≠ −1, respectively.

Recent studies after the five-year data release combined the WMAP five-year data with information on the growth of structure to find various improved limits. Vikhlinin et al. (2009b) added the abundance of X-ray-selected clusters of galaxies, which were found in the ROSAT All Sky Survey and followed up by the Chandra X-ray Observatory (their cluster catalog is described in Vikhlinin et al. 2009a), to the WMAP five-year data, the BAO measurement from Eisenstein et al. (2005), and the Type Ia supernova data from Davis et al. (2007), to find ∑m_ν < 0.33 eV (95% CL) for w ≠ −1. Mantz et al. (2010b) added a different cluster catalog, also derived from the ROSAT All Sky Survey and followed up by the Chandra X-ray Observatory (their cluster catalog is described in Mantz et al. 2010a), and the measurement of the gas mass fraction of relaxed clusters (Allen et al. 2008) to the WMAP five-year data, the BAO measurement from Percival et al. (2007), and the "Union" Type Ia supernova samples from Kowalski et al. (2008) (all of which constitute the five-year "WMAP+BAO+SN" set in Komatsu et al. 2009a), to find ∑m_ν < 0.33 and 0.43 eV (95% CL) for w = −1 and w ≠ −1, respectively.

Reid et al. (2010a) added a prior on the amplitude of matter density fluctuations,  σ₈(Ω_m/0.25)^0.41 = 0.832 ± 0.033 (68% CL; Rozo et al. 2010), which was derived from the abundance of optically selected clusters of galaxies called the "maxBCG cluster catalog" (Koester et al. 2007), to the five-year WMAP+BAO+SN, and found ∑m_ν < 0.35 and 0.52 eV (95% CL) for w = −1 and w ≠ −1, respectively. Thomas et al. (2010) added the angular power spectra of photometrically selected samples of LRGs called "MegaZ" to the five-year WMAP+BAO+SN, and found ∑m_ν < 0.325 eV (95% CL) for w = −1. Wang et al. (2005) pointed out that the limit on ∑m_ν from galaxy clusters would improve significantly by not only using the abundance but also the power spectrum of clusters.

In order to exploit the full information contained in the growth of structure, it is essential to understand the effects of massive neutrinos on the nonlinear growth. All of the work to date (including WMAP+LRG+H₀ presented above) included the effects of massive neutrinos on the linear growth, while ignoring their nonlinear effects. The widely used phenomenological calculation of the nonlinear matter power spectrum called the HALOFIT (Smith et al. 2003) has not been calibrated for models with massive neutrinos. Consistent treatments of massive neutrinos in the nonlinear structure formation using cosmological perturbation theory (Saito et al. 2008, 2009; Wong 2008; Lesgourgues et al. 2009; Shoji & Komatsu 2009) and numerical simulations (Brandbyge et al. 2008; Brandbyge & Hannestad 2009) have just begun to be explored. More work along these lines would be necessary to exploit the information on the growth structure to constrain the mass of neutrinos.

How many relativistic species are there in the universe after the matter-radiation equality epoch? We parameterize the relativistic dof using the effective number of neutrino species, N_eff, given in Equation (23). This quantity can be written in terms of the matter density, Ω_mh², and the redshift of matter-radiation equality, z_eq, as (see Equation (84) of Komatsu et al. 2009a)

$$\begin{equation} N_{\rm eff} = 3.04 + 7.44 \left(\frac{\Omega _mh^2}{0.1308}\frac{3139}{1+z_{\rm eq}}-1\right). \end{equation} \tag{ 53 }$$

(Here, Ω_mh² = 0.1308 and z_eq = 3138 are the five-year maximum likelihood values from the simplest ΛCDM model.) This formula suggests that the variation in N_eff is given by

$$\begin{equation} \frac{\delta N_{\rm eff}}{N_{\rm eff}} \simeq 2.45 \frac{\delta (\Omega _mh^2)}{\Omega _mh^2} - 2.45\frac{\delta z_{\rm eq}}{1+z_{\rm eq}}. \end{equation} \tag{ 54 }$$

The equality redshift is one of the direct observables from the temperature power spectrum. The WMAP data constrain z_eq mainly from the ratio of the first peak to the third peak. As the seven-year temperature power spectrum has a better determination of the amplitude of the third peak (Larson et al. 2011), we expect a better limit on z_eq compared to the five-year one. For models where N_eff is different from 3.04, we find z_eq = 3145⁺¹⁴⁰₋₁₃₉ (68% CL) from the WMAP data only,²⁹ which is better than the five-year limit by more than 10% (see Table 8).

However, the fractional error in Ω_mh² is much larger, and thus we need to determine Ω_mh² using external data. The BAO data provide one constraint. We also find that Ω_mh² and H₀ are strongly correlated in the models with N_eff ≠ 3.04 (see Figure 9). Therefore, an improved measurement of H₀ from Riess et al. (2009) would help reduce the error in Ω_mh², thereby reducing the error in N_eff. The limit on Ω_mh² from the seven-year WMAP+BAO+H₀ combination is better than the five-year "WMAP+BAO+SN+HST" limit by 36%.

Zoom In Zoom Out Reset image size

We find that the WMAP+BAO+H₀ limit on N_eff is

$$$ N_{\rm eff} = 4.34^{+ 0.86}_{- 0.88} \mbox{(68\% CL)}, $$$

while the WMAP+LRG+H₀ limit is N_eff = 4.25^+0.76_−0.80 (68% CL), which are significantly better than the five-year WMAP+BAO+SN+HST limit, N_eff = 4.4 ± 1.5 (68% CL).

Reid et al. (2010a) added the maxBCG prior,  σ₈(Ω_m/0.25)^0.41 = 0.832 ± 0.033 (68% CL; Rozo et al. 2010), to the five-year WMAP+BAO+SN+HST, and found N_eff = 3.5 ± 0.9 (68% CL). They also added the above prior to the five-year version of WMAP+LRG+H₀, finding N_eff = 3.77 ± 0.67 (68% CL).

The constraint on N_eff can also be interpreted as an upper bound on the energy density in primordial gravitational waves with frequencies >10⁻¹⁵ Hz. Many cosmological mechanisms for the generation of stochastic gravitational waves exist, such as certain inflationary models, electroweak phase transitions, and cosmic strings. At low frequencies (10⁻¹⁷–10⁻¹⁶ Hz), the background is constrained by the limit on tensor fluctuations described in Section 4.1. Constraints at higher frequencies come from pulsar timing measurements at ∼10⁻⁸ Hz (Jenet et al. 2006), recent data from the Laser Interferometer Gravitational Wave Observatory (LIGO) at 100 Hz (with limits of Ω_gw < 6.9 × 10⁻⁶ Abbott et al. 2009), and at frequencies >10⁻¹⁰ Hz from measurements of light-element abundances. A large gravitational wave energy density at nucleosynthesis would alter the predicted abundances, and observations imply an upper bound of Ω_gwh² < 7.8 × 10⁻⁶ (Cyburt et al. 2005).

The CMB provides a limit that reaches down to 10⁻¹⁵ Hz, corresponding to the comoving horizon at recombination. The gravitational wave background within the horizon behaves as free-streaming massless particles, so affects the CMB and matter power spectra in the same way as massless neutrinos (Smith et al. 2006). The density contributed by N_gw massless neutrino species is Ω_gwh² = 5.6 × 10⁻⁶N_gw. Constraints have been found using the WMAP three-year data combined with additional cosmological probes by Smith et al. (2006), for both adiabatic and homogeneous initial conditions for the tensor perturbations. With the current WMAP+BAO+H₀ data combination, we define $N_{\rm {gw}} = N_{\rm {eff}}-3.04$, and find limits of

$$$ N_{{\rm gw}}< 2.85, \quad \Omega _{{\rm gw}}h^2 < 1.60 \times 10^{-5} (95\% {\rm CL}) $$$

for adiabatic initial conditions, imposing an $N_{\rm {eff}} \ge 3.04$ prior. Adiabatic conditions might be expected if the gravitational waves were generated by the appearance of cusps in cosmic strings (Damour & Vilenkin 2000, 2001; Siemens et al. 2006). For the WMAP+LRG+H₀ data, we find N_gw < 2.64, or Ω_gwh² < 1.48 × 10⁻⁵ at 95% CL. Given a particular string model, these bounds can be used to constrain the cosmic string tension (e.g., Siemens et al. 2007; Copeland & Kibble 2009).

A change in the primordial helium abundance affects the shape of the temperature power spectrum (Hu et al. 1995). The most dominant effect is a suppression of the power spectrum at l ≳ 500 due to an enhanced Silk damping effect.

For a given mass density of baryons (protons and helium nuclei), the number density of electrons, n_e, can be related to the primordial helium abundance. When both hydrogen and helium were ionized, n_e = (1 − Y_p/2)ρ_b/m_p. However, most of the helium recombines by z ∼ 1800 (see Switzer & Hirata 2008, and references therein), much earlier than the photon decoupling epoch, z = 1090. As a result, the number density of free electrons at around the decoupling epoch is given by n_e = (1 − Y_p)ρ_b/m_p ∝ (1 − Y_p)Ω_bh² (Hu et al. 1995). The larger Y_p is, the smaller n_e becomes. If the number of electrons is reduced, photons can free-stream longer (the mean free path of photons, 1/( σ_Tn_e), gets larger), wiping out more temperature anisotropy. Therefore, a larger Y_p results in a greater suppression of power on small angular scales.

Ichikawa et al. (2008; also see Ichikawa & Takahashi 2006) show that a 100% change in Y_p changes the heights of the second, third, and fourth peaks by ≈1%, 3%, and 3%, respectively. Therefore, one expects that a combination of the WMAP data and small-scale CMB experiments such as ACBAR (Reichardt et al. 2009) and QUaD (Brown et al. 2009) would be a powerful probe of the primordial helium abundance.

In Figure 10, we compare the WMAP, ACBAR, and QUaD data with the temperature power spectrum with the nominal value of the primordial helium abundance, Y_p = 0.24 (pink line), and that with a tiny amount of helium, Y_p = 0.01 (blue line). There is too much power in the case of Y_p = 0.01, making it possible to detect the primordial helium effect using the CMB data alone.

Zoom In Zoom Out Reset image size

However, one must be careful about a potential degeneracy between the effect of helium and those of the other cosmological parameters. First, as the number density of electrons is given by n_e = (1 − Y_p)n_b ∝ (1 − Y_p)Ω_bh², Y_p and Ω_bh² may be correlated. Second, a scale-dependent suppression of power such as this may be correlated with the effect of tilt, n_s (Trotta & Hansen 2004).

In the left panel of Figure 11, we show that Ω_bh² and Y_p are essentially uncorrelated: the baryon density is determined by the first-to-second peak ratio relative to the first-to-third peak ratio, which is now well measured by the WMAP data. Therefore, the current WMAP data allow Ω_bh² to be measured regardless of Y_p.

Zoom In Zoom Out Reset image size

In the middle panel of Figure 11, we show that there is a slight positive correlation between n_s and Y_p: an enhanced Silk damping produced by a larger Y_p can be partially canceled by a larger n_s (Trotta & Hansen 2004).

We find a 95% CL upper limit of Y_p < 0.51 from the WMAP data alone. When we add the ACBAR and QUaD data, we find a significant detection of the effect of primordial helium by more than 3σ (see the right panel of Figure 11),

$$$ Y_p=0.326\pm 0.075 \mbox{(68\% CL)}. $$$

The 95% CL limit is 0.16 < Y_p < 0.46. The 99% CL lower limit is Y_p>0.11. This value is broadly consistent with the helium abundances estimated from observations of low-metallicity extragalactic ionized (H ii) regions, Y_p ≃ 0.24–0.25 (Gruenwald et al. 2002; Izotov & Thuan 2004; Olive & Skillman 2004; Fukugita & Kawasaki 2006; Peimbert et al. 2007). See Steigman (2007) for a review.

We can improve this limit by imposing an upper limit on Y_p from these astrophysical measurements. As the helium is created by nuclear fusion in stars, the helium abundances measured from stars (e.g., the Sun; see Asplund et al. 2009, for a recent review) and H ii regions are, in general, larger than the primordial abundance. On the other hand, as we have just shown, the CMB data provide a lower limit on Y_p. Even with a very conservative hard prior, Y_p < 0.3, we find 0.23 < Y_p < 0.3 (68% CL)³⁰. Therefore, a combination of the CMB and the solar constraints on Y_p offers a new way for testing the predictions of theory of the big bang nucleosynthesis (BBN). For example, the BBN predicts that the helium abundance is related to the baryon-to-photon ratio, η, and the number of additional neutrino species (or any other additional relativistic dof) during the BBN epoch, ΔN_ν ≡ N_ν − 3, as (see Equation (11) of Steigman 2008)

$$\begin{equation} Y_p = 0.2485 + 0.0016[(\eta _{10} - 6) + 100(S - 1)], \end{equation} \tag{ 55 }$$

where $S\equiv \sqrt{1+(7/43)\Delta N_\nu }\simeq 1+0.081\Delta N_\nu$ and η₁₀ ≡ 10¹⁰η = 273.9(Ω_bh²) = 6.19 ± 0.15 (68% CL; WMAP+BAO+H₀). (See Simha & Steigman 2008, for more discussion on this method.) For ΔN_ν = 1, the helium abundance changes by ΔY_p = 0.013, which is much smaller than our error bar, but is comparable to the expected error bar from Planck (Ichikawa et al. 2008).

There have been several attempts to measure Y_p from the CMB data (Trotta & Hansen 2004; Huey et al. 2004; Ichikawa & Takahashi 2006; Ichikawa et al. 2008; Dunkley et al. 2009). The previous best-limit is Y_p = 0.25^+0.10(+0.15)_{−0.07(−0.17)} at 68% CL (95% CL), which was obtained by Ichikawa et al. (2008) from the WMAP five-year data combined with ACBAR (Reichardt et al. 2009), BOOMERanG (Jones et al. 2006; Piacentini et al. 2006; Montroy et al. 2006), and Cosmic Background Imager (CBI; Sievers et al. 2007). Note that the likelihood function of Y_p is non-Gaussian, with a tail extending to Y_p = 0; thus, the level of significance of detection was less than 3σ.

In a flat universe, Ω_k = 0, an accurate determination of H₀ helps improve a limit on a constant equation of state, w (Spergel et al. 2003; Hu 2005). Using WMAP+BAO+H₀, we find

$$$ w = -1.10\pm 0.14 \mbox{(68\% CL)}, $$$

which improves to w = −1.08 ± 0.13 (68% CL) if we add the time-delay distance out to the lens system B1608+656 (Suyu et al. 2010, see Section 3.2.5). These limits are independent of high-z Type Ia supernova data.

The high-z supernova data provide the most stringent limit on w. Using WMAP+BAO+SN, we find w = −0.980 ± 0.053 (68% CL). The error does not include systematic errors in supernovae, which are comparable to the statistical error (Kessler et al. 2009; Hicken et al. 2009a); thus, the error in w from WMAP+BAO+SN is about a half of that from WMAP+BAO+H₀ or WMAP+BAO+H₀+D_Δt.

The cluster abundance data are sensitive to w via the comoving volume element, angular diameter distance, and growth of matter density fluctuations (Haiman et al. 2001). By combining the cluster abundance data and the five-year WMAP data, Vikhlinin et al. (2009b) found w = −1.08 ± 0.15(stat) ± 0.025(syst) (68% CL) for a flat universe. By adding BAO of Eisenstein et al. (2005) and the supernova data of Davis et al. (2007), they found w = −0.991 ± 0.045(stat) ± 0.039(syst) (68% CL). These results using the cluster abundance data (also see Mantz et al. 2010c) agree well with our corresponding WMAP+BAO+H₀ and WMAP+BAO+SN limits.

When Ω_k ≠ 0, limits on w significantly weaken, with a tail extending to large negative values of w, unless supernova data are added.

In Figure 12, we show that WMAP+BAO+H₀ allows for w ≲ −2, which can be excluded by adding information on the time-delay distance. In both cases, the spatial curvature is well constrained: we find Ω_k = −0.0125^+0.0064_−0.0067 from WMAP+BAO+H₀, and −0.0111^+0.0060_−0.0063 (68% CL) from WMAP+BAO+H₀+D_Δt, whose errors are comparable to that of the WMAP+BAO+H₀ limit on Ω_k with w = −1, Ω_k = −0.0023^+0.0054_−0.0056 (68% CL; see Section 4.3).

Zoom In Zoom Out Reset image size

However, w is poorly constrained: we find w = −1.44 ±  0.27 from WMAP+BAO+H₀, and −1.40 ± 0.25 (68% CL) from WMAP+BAO+H₀+D_Δt.

Among the data combinations that do not use the information on the growth of structure, the most powerful combination for constraining Ω_k and w simultaneously is a combination of the WMAP data, BAO (or D_Δt), and supernovae, as WMAP+BAO (or D_Δt) primarily constrains Ω_k, and WMAP+SN primarily constrains w. With WMAP+BAO+SN, we find w = −0.999^+0.057_−0.056 and Ω_k = −0.0057^+0.0066_−0.0068 (68% CL). Note that the error in the curvature is essentially the same as that from WMAP+BAO+H₀, while the error in w is ∼4 times smaller.

Vikhlinin et al. (2009b) combined their cluster abundance data with the five-year WMAP+BAO+SN to find w = −1.03 ± 0.06 (68% CL) for a curved universe. Reid et al. (2010b) combined their LRG power spectrum with the five-year WMAP data and the Union supernova data to find w = −0.99 ± 0.11 and Ω_k = −0.0109 ± 0.0088 (68% CL). These results are in good agreement with our seven-year WMAP+BAO+SN limit.

As for a time-dependent equation of state, we shall find constraints on the present-day value of the equation of state and its derivative using a linear form, w(a) = w₀ + w_a(1 − a) (Chevallier & Polarski 2001; Linder 2003). We assume a flat universe, Ω_k = 0. (For recent limits on w(a) with Ω_k ≠ 0, see Wang 2009, and references therein.) While we have constrained this model using the WMAP distance prior in the five-year analysis (see Section 5.4.2 of Komatsu et al. 2009a), in the seven-year analysis we shall present the full Markov Chain Monte Carlo exploration of this model.

For a time-dependent equation of state, one must be careful about the treatment of perturbations in dark energy when w crosses −1. We use the "parameterized post-Friedmann" (PPF) approach, implemented in the CAMB code following Fang et al. (2008).³¹

In Figure 13, we show the seven-year constraints on w₀ and w_a from WMAP+H₀+SN (red), WMAP+BAO+H₀+SN (blue), and WMAP+BAO+H₀+D_Δt+SN (black). The angular diameter distances measured from BAO and D_Δt help exclude models with large negative values of w_a. We find that the current data are consistent with a cosmological constant, even when w is allowed to depend on time. However, a large range of values of (w₀, w_a) are still allowed by the data: we find

$$$ w_0=-0.93\pm 0.13 \ \mbox{and} \ w_a=-0.41^{+0.72}_{-0.71} \mbox{(68\% CL)}, $$$

from WMAP+BAO+H₀+SN. When the time-delay distance information is added, the limits improve to w₀ = −0.93 ± 0.12 and w_a = −0.38^+0.66_−0.65 (68% CL).

Zoom In Zoom Out Reset image size

Vikhlinin et al. (2009b) combined their cluster abundance data with the five-year WMAP+BAO+SN to find a limit on a linear combination of the parameters, w_a + 3.64(1 + w₀) = 0.05 ± 0.17 (68% CL). Our data combination is sensitive to a different linear combination: we find w_a + 5.14(1 + w₀) = −0.05 ± 0.32 (68% CL) for the seven-year WMAP+BAO+H₀+SN combination.

The current data are consistent with a flat universe dominated by a cosmological constant.

The growth of cosmological density fluctuations is a powerful probe of dark energy, modified gravity, and massive neutrinos. The WMAP data provide a useful normalization of the cosmological perturbation at the decoupling epoch, z = 1090. By comparing this normalization with the amplitude of matter density fluctuations in a low redshift universe, one may distinguish between dark energy and modified gravity (Ishak et al. 2006; Koyama & Maartens 2006; Amarzguioui et al. 2006; Doré et al. 2007; Linder & Cahn 2007; Upadhye 2007; Zhang et al. 2007; Yamamoto et al. 2007; Chiba & Takahashi 2007; Bean et al. 2007; Hu & Sawicki 2007; Song et al. 2007; Starobinsky 2007; Daniel et al. 2008; Jain & Zhang 2008; Bertschinger & Zukin 2008; Amin et al. 2008; Hu 2008) and determine the mass of neutrinos (Hu et al. 1998; Lesgourgues & Pastor 2006).

In Section 5.5 of Komatsu et al. (2009a), we provided a "WMAP normalization prior," which is a constraint on the power spectrum of curvature perturbation, $\Delta ^2_{\cal R}$. Vikhlinin et al. (2009b) combined this with the number density of clusters of galaxies to constrain the dark energy equation of state, w, and the amplitude of matter density fluctuations,  σ₈.

The matter density fluctuation in Fourier space, δ_m,k, is related to ${\cal R}_{\mathbf k}$ as $\delta _{m,{\mathbf k}}(z) = \frac{2k^3}{5H_0^2\Omega _m}{\cal R}_{\mathbf k}T(k)D(k,z)$, where D(k,z) and T(k) are the linear growth rate and the matter transfer function normalized such that T(k) → 1 as k → 0, and (1 + z)D(k, z) → 1 as k → 0 during the matter era, respectively. Ignoring the mass of neutrinos and modifications to gravity, one can obtain the growth rate by solving a single differential equation (Wang & Steinhardt 1998; Linder & Jenkins 2003).³²

The seven-year normalization prior is

$$$ \Delta _{\cal R}^2(k_{\scriptsize\textit{WMAP}}) = (2.208\pm 0.078)\times 10^{-9} \mbox{(68\% CL)}, $$$

where $k_{\scriptsize\textit{WMAP}}=0.027 \,{\rm Mpc}^{-1}$. For comparison, the five-year normalization prior was $\Delta _{\cal R}^2(0.02 \,{\rm Mpc}^{-1}) = (2.21\pm 0.09)\times 10^{-9}$. This normalization prior is valid for models with Ω_k ≠ 0, w ≠ −1, or m_ν>0. However, these normalizations cannot be used for the models that have the tensor modes, r>0, or the running index, dn_s/dln k ≠ 0.

The temperature power spectrum of CMB is sensitive to the physics at the decoupling epoch, z = 1090, as well as the physics between now and the decoupling epoch. The former primarily affects the amplitude of acoustic peaks, i.e., the ratios of the peak heights and the Silk damping. The latter changes the locations of peaks via the angular diameter distance out to the decoupling epoch. One can quantify this by (1) the "acoustic scale," l_A,

$$\begin{equation} l_A=(1+z_*)\frac{\pi D_A(z_*)}{r_s(z_*)}, \end{equation} \tag{ 56 }$$

where z_* is the redshift of decoupling, for which we use the fitting formula of Hu & Sugiyama (1996), as well as by (2) the "shift parameter," R (Bond et al. 1997),

$$\begin{equation} R =\frac{\sqrt{\Omega _mH_0^2}}{c}(1+z_*)D_A(z_*). \end{equation} \tag{ 57 }$$

These two parameters capture most of the constraining power of the WMAP data for dark energy properties (Wang & Mukherjee 2007; Wright 2007; Elgarøy & Multamäki 2007; Corasaniti & Melchiorri 2008), with one important difference. The distance prior does not capture the information on the growth of structure probed by the late-time ISW effect. As a result, the dark energy constraints derived from the distance prior are similar to, but weaker than, those derived from the full analysis (Komatsu et al. 2009a; Li et al. 2008).

We must understand the limitation of this method. Namely, the distance prior is applicable only when the model in question is based on

• 1.  
  the standard Friedmann–Lemaitre–Robertson–Walker universe with matter, radiation, dark energy, as well as spatial curvature,
• 2.  
  neutrinos with the effective number of neutrinos equal to 3.04, and the minimal mass (m_ν ∼ 0.05 eV), and
• 3.  
  nearly power-law primordial power spectrum of curvature perturbations, |dn_s/dln k| ≪ 0.01, negligible primordial gravitational waves relative to the curvature perturbations, r ≪ 0.1, and negligible entropy fluctuations relative to the curvature perturbations, α ≪ 0.1.

In Tables 9 and 10, we provide the seven-year distance prior. The errors in l_A, R, and z_* have improved from the five-year values by 12%, 5%, and 2%, respectively. To compute the likelihood, use

$$\begin{equation} -2\ln L = \sum _{ij}(x_i-d_i)(C^{-1})_{ij}(x_j-d_j), \end{equation} \tag{ 58 }$$

where x_i = (l_A, R, z_*) is the values predicted by a model in question, $d_i=(l_A^{\scriptsize\textit{WMAP}},R^{\scriptsize\textit{WMAP}},z_*^{\scriptsize\textit{WMAP}})$ is the data given in Table 9, and C⁻¹_ij is the inverse covariance matrix given in Table 10. Also see Section 5.4.1 of Komatsu et al. (2009a) for more information.

During the period of cosmic inflation (Starobinskiiˇ 1979; Starobinsky 1982; Guth 1981; Sato 1981; Linde 1982; Albrecht & Steinhardt 1982), quantum fluctuations were generated and became the seeds for the cosmic structures that we observe today (Mukhanov & Chibisov 1981; Hawking 1982; Starobinsky 1982; Guth & Pi 1982; Bardeen et al. 1983). (Also see Linde 2008, 1990, Mukhanov et al. 1992, Liddle & Lyth 2000, 2009, Bassett et al. 2006 for reviews.)

Inflation predicts that the statistical distribution of primordial fluctuations is nearly a Gaussian distribution with random phases. Measuring deviations from a Gaussian distribution, i.e., non-Gaussian correlations in primordial fluctuations, is a powerful test of inflation, as how precisely the distribution is (non-)Gaussian depends on the detailed physics of inflation (see Bartolo et al. 2004; Komatsu et al. 2009b, for reviews).

In this paper, we constrain the amplitude of non-Gaussian correlations using the angular bispectrum of CMB temperature anisotropy, the harmonic transform of the three-point correlation function (see Komatsu 2001, for a review). The observed angular bispectrum is related to the three-dimensional bispectrum of primordial curvature perturbations, $\langle \zeta _{{\mathbf k}_1}\zeta _{{\mathbf k}_2}\zeta _{{\mathbf k}_3}\rangle =(2\pi)^3\delta ^D({\mathbf k}_1+{\mathbf k}_2+{\mathbf k}_3)B_\zeta (k_1,k_2,k_3)$. In the linear order, the primordial curvature perturbation is related to Bardeen's curvature perturbation (Bardeen 1980) in the matter-dominated era, Φ, by $\zeta =\frac{5}{3}\Phi$ (e.g., Kodama & Sasaki 1984). The CMB temperature anisotropy in the Sachs–Wolfe limit (Sachs & Wolfe 1967) is given by $\Delta T/T=-\frac{1}{3}\Phi =-\frac{1}{5}\zeta$. We write the bispectrum of Φ as

$$\begin{equation} \langle \Phi _{{\mathbf k}_1}\Phi _{{\mathbf k}_2}\Phi _{{\mathbf k}_3}\rangle =(2\pi)^3\delta ^D({\mathbf k}_1+{\mathbf k}_2+{\mathbf k}_3)F(k_1,k_2,k_3). \end{equation} \tag{ 59 }$$

We shall explore three different shapes of the primordial bispectrum: "local," "equilateral," and "orthogonal." They are defined as follows:

• 1.  
  Local form. The local form bispectrum is given by (Gangui et al. 1994; Verde et al. 2000; Komatsu & Spergel 2001)

  $$\begin{eqnarray} \nonumber F_{\rm local}(k_1,k_2,k_3) &=& 2f_{\scriptsize\textit{NL}}^{\rm local}[P_\Phi (k_1)P_\Phi (k_2)+P_\Phi (k_2)P_\Phi (k_3)\\ \nonumber & &+P_\Phi (k_3)P_\Phi (k_1)]\\ &=& 2A^2f_{\scriptsize\textit{NL}}^{\rm local}\left[\frac{1}{k^{4-n_s}_1k^{4-n_s}_2}+(2 \mbox{perm}.)\right],\nonumber\\ \end{eqnarray} \tag{ 60 }$$

  where $P_\Phi =A/k^{4-n_s}$ is the power spectrum of Φ with a normalization factor A. This form is called the local form, as this bispectrum can arise from the curvature perturbation in the form of $\Phi =\Phi _L+f_{\scriptsize\textit{NL}}^{\rm local}\Phi ^2_L$, where both sides are evaluated at the same location in space (Φ_L is a linear Gaussian fluctuation).³³ Equation (60) peaks at the so-called squeezed triangle for which k₃ ≪ k₂ ≈ k₁ (Babich et al. 2004). In this limit, we obtain

  $$\begin{equation} F_{\rm local}(k_1,k_1,k_3\rightarrow 0)=4f_{\scriptsize\textit{NL}}^{\rm local}P_\Phi (k_1)P_\Phi (k_3). \end{equation} \tag{ 61 }$$

  How large is $f_{\scriptsize\textit{NL}}^{\rm local}$ from inflation? The earlier calculations showed that $f_{\scriptsize\textit{NL}}^{\rm local}$ from single-field slow-roll inflation would be of order the slow-roll parameter, ∼ 10⁻² (Salopek & Bond 1990; Falk et al. 1993; Gangui et al. 1994). More recently, Maldacena (2003) and Acquaviva et al. (2003) found that the coefficient of P_Φ(k₁)P_Φ(k₃) from the simplest single-field slow-roll inflation with the canonical kinetic term in the squeezed limit is given by

  $$\begin{equation} F_{\rm local}(k_1,k_1,k_3\rightarrow 0)=\frac{5}{3}(1-n_s) P_\Phi (k_1)P_\Phi (k_3). \end{equation} \tag{ 62 }$$

  Comparing this result with the form predicted by the $f_{\scriptsize\textit{NL}}^{\rm local}$ model, one obtains $f_{\scriptsize\textit{NL}}^{\rm local}=(5/12)(1-n_s)$, which gives $f_{\scriptsize\textit{NL}}^{\rm local}=0.015$ for n_s = 0.963.
• 2.  
  Equilateral form. The equilateral form bispectrum is given by (Creminelli et al. 2006)

  $$\begin{eqnarray} \nonumber && F_{\rm equil}(k_1,k_2,k_3)= 6A^2f_{\scriptsize\textit{NL}}^{\rm equil}\left\lbrace -\frac{1}{k^{4-n_s}_1k^{4-n_s}_2}-\frac{1}{k^{4-n_s}_2k^{4-n_s}_3} \right.\\ \nonumber &&\quad \left. - \frac{1}{k^{4-n_s}_3k^{4-n_s}_1} -\frac{2}{(k_1k_2k_3)^{2(4-n_s)/3}} \right.\\ &&\quad \left. + \left[\frac{1}{k^{(4-n_s)/3}_1k^{2(4-n_s)/3}_2k^{4-n_s}_3} +\mbox{(5 perm.)}\right]\right\rbrace. \end{eqnarray} \tag{ 63 }$$

  This function approximates the bispectrum forms that arise from a class of inflation models in which scalar fields have non-canonical kinetic terms. One example is the so-called Dirac–Born–Infeld inflation (Silverstein & Tong 2004; Alishahiha et al. 2004), which gives $f_{\scriptsize\textit{NL}}^{\rm equil}\propto -1/c_s^2$ in the limit of c_s ≪ 1, where c_s is the effective sound speed at which scalar field fluctuations propagate. There are various other models that can produce $f_{\scriptsize\textit{NL}}^{\rm equil}$ (Arkani-Hamed et al. 2004; Seery & Lidsey 2005; Chen et al. 2007; Cheung et al. 2008; Li et al. 2008). The local and equilateral forms are nearly orthogonal to each other, which means that both can be measured nearly independently.
• 3.  
  Orthogonal form. The orthogonal form, which is constructed such that it is nearly orthogonal to both the local and equilateral forms, is given by (Senatore et al. 2010)

  $$\begin{eqnarray} \nonumber && F_{\rm orthog}(k_1,k_2,k_3)\,{=}\, 6A^2f_{\scriptsize\textit{NL}}^{\rm orthog}\!\left\lbrace\! -\frac{3}{k^{4-n_s}_1k^{4-n_s}_2}\,{-}\,\frac{3}{k^{4-n_s}_2k^{4-n_s}_3} \right.\\ \nonumber & &\quad \left. -\frac{3}{k^{4-n_s}_3k^{4-n_s}_1} -\frac{8}{(k_1k_2k_3)^{2(4-n_s)/3}} \right.\\ & &\quad \left. +\left[\frac{3}{k^{(4-n_s)/3}_1k^{2(4-n_s)/3}_2k^{4-n_s}_3} +\mbox{(5 perm.)}\right]\right\rbrace. \end{eqnarray} \tag{ 64 }$$

  This form approximates the forms that arise from a linear combination of higher-derivative scalar-field interaction terms, each of which yields forms similar to the equilateral shape. Senatore et al. (2010) found that, using the "effective field theory of inflation" approach (Cheung et al. 2008), a certain linear combination of similarly equilateral shapes can yield a distinct shape which is orthogonal to both the local and equilateral forms.

Note that these are not the most general forms one can write down, and there are other forms which would probe different aspects of the physics of inflation (Moss & Xiong 2007; Moss & Graham 2007; Chen et al. 2007; Holman & Tolley 2008; Chen & Wang 2010; Chen & Wang 2010).

Of these forms, the local form bispectrum has special significance. Creminelli & Zaldarriaga (2004) showed that not only models with the canonical kinetic term, but all single-inflation models predict the bispectrum in the squeezed limit given by Equation (62), regardless of the form of potential, kinetic term, slow-roll, or initial vacuum state (also see Seery & Lidsey 2005; Chen et al. 2007; Cheung et al. 2008). This means that a convincing detection of $f_{\scriptsize\textit{NL}}^{\rm local}$ would rule out all single-field inflation models.

The first limit on $f_{\scriptsize\textit{NL}}^{\rm local}$ was obtained from the COBE 4-year data (Bennett et al. 1996) by Komatsu et al. (2002), using the angular bispectrum. The limit was improved by an order of magnitude when the WMAP first year data were used to constrain $f_{\scriptsize\textit{NL}}^{\rm local}$ (Komatsu et al. 2003). Since then the limits have improved steadily as WMAP collect more years of data and the bispectrum method for estimating $f_{\scriptsize\textit{NL}}^{\rm local}$ has improved (Komatsu et al. 2005; Creminelli et al. 2006, 2007; Spergel et al. 2007; Yadav & Wandelt 2008; Komatsu et al. 2009a; Smith et al. 2009).³⁴

In this paper, we shall adopt the optimal estimator (developed by Babich 2005; Creminelli et al. 2006, 2007; Smith & Zaldarriaga 2006; Yadav et al. 2008), which builds on and significantly improves the original bispectrum estimator proposed by Komatsu et al. (2005), especially when the spatial distribution of instrumental noise is not uniform. For details of the method, see Appendix A of Smith et al. (2009) for $f_{\scriptsize\textit{NL}}^{\rm local}$, and Section 4.1 of Senatore et al. (2010) for $f_{\scriptsize\textit{NL}}^{\rm equil}$ and $f_{\scriptsize\textit{NL}}^{\rm orthog}$. To construct the optimal estimators, we need to specify the cosmological parameters. We use the five-year ΛCDM parameters from WMAP+BAO+SN, for which n_s = 0.96.

We also constrain the bispectrum due to residual (unresolved) point sources, b_src. The optimal estimator for b_src is constructed by replacing a_lm/C_l in Equation (A24) of Komatsu et al. (2009a) with (C⁻¹a)_lm, and using their Equations (A17) and (A5). The C⁻¹ matrix is computed by the multigrid-based algorithm of Smith et al. (2007).

We use the V- and W-band maps at the HEALPix resolution N_side = 1024. As the optimal estimator weights the data optimally at all multipoles, we no longer need to choose the maximum multipole used in the analysis, i.e., we use all the data. We use both the raw maps (before cleaning foreground) and foreground-reduced (clean) maps to quantify the foreground contamination of f_NL parameters. For all cases, we find the best limits on f_NL parameters by combining the V- and W-band maps, and marginalizing over the synchrotron, free–free, and dust foreground templates (Gold et al. 2011). As for the mask, we always use the KQ75y7 mask (Gold et al. 2011).

In Table 11, we summarize our results.

• 1.  
  Local form results. The seven-year best estimate of $f_{\scriptsize\textit{NL}}^{\rm local}$ is

  $$$ f_{\scriptsize\textit{NL}}^{\rm local}=32\pm 21 \mbox{(68\% CL)}. $$$

  The 95% limit is $-10<f_{\scriptsize\textit{NL}}^{\rm local}<74$. When the raw maps are used, we find $f_{\scriptsize\textit{NL}}^{\rm local}=59\pm 21$ (68% CL). When the clean maps are used, but foreground templates are not marginalized over, we find $f_{\scriptsize\textit{NL}}^{\rm local}=42\pm 21$ (68% CL). These results (in particular the clean-map versus the foreground marginalized) indicate that the foreground emission makes a difference at the level of $\Delta f_{\scriptsize\textit{NL}}^{\rm local}\sim 10$.³⁵ We find that the V + W result is lower than the V-band or W-band results. This is possible, as the V + W result contains contributions from the cross-correlations of V and W such as 〈VVW〉 and 〈VWW〉.
• 2.  
  Equilateral form results. The seven-year best estimate of $f_{\scriptsize\textit{NL}}^{\rm equil}$ is

  $$$ f_{\scriptsize\textit{NL}}^{\rm equil}=26\pm 140 \mbox{(68\% CL)}. $$$

  The 95% limit is $-214<f_{\scriptsize\textit{NL}}^{\rm equil}<266$. For $f_{\scriptsize\textit{NL}}^{\rm equil}$, the foreground marginalization does not shift the central values very much, $\Delta f_{\scriptsize\textit{NL}}^{\rm equil}=-3$. This makes sense, as the equilateral bispectrum does not couple small-scale modes to very large-scale modes l ≲ 10, which are sensitive to the foreground emission. On the other hand, the local form bispectrum is dominated by the squeezed triangles, which do couple large- and small-scale modes.
• 3.  
  Orthogonal form results. The seven-year best estimate of $f_{\scriptsize\textit{NL}}^{\rm orthog}$ is

  $$$ f_{\scriptsize\textit{NL}}^{\rm orthog}=-202\pm 104 \mbox{(68\% CL)}. $$$

  The 95% limit is $-410<f_{\scriptsize\textit{NL}}^{\rm orthog}<6$. The foreground marginalization has little effect, $\Delta f_{\scriptsize\textit{NL}}^{\rm orthog}=-4$.

Table 11. Estimates^a and the Corresponding 68% Intervals of the Primordial non-Gaussianity Parameters ($f_{\scriptsize\textit{NL}}^{\rm local}$, $f_{\scriptsize\textit{NL}}^{\rm equil}$, $f_{\scriptsize\textit{NL}}^{\rm orthog}$) and the Point-source Bispectrum Amplitude, b_src (in units of 10⁻⁵μK³sr²), from the WMAP Seven-year Temperature Maps

Band	Foregroundb	$f_{\scriptsize\textit{NL}}^{\rm local}$	$f_{\scriptsize\textit{NL}}^{\rm equil}$	$f_{\scriptsize\textit{NL}}^{\rm orthog}$	bsrc
V + W	Raw	59 ± 21	33 ± 140	−199 ± 104	N/A
V + W	Clean	42 ± 21	29 ± 140	−198 ± 104	N/A
V + W	Marg.c	32 ± 21	26 ± 140	−202 ± 104	−0.08 ± 0.12
V	Marg.	43 ± 24	64 ± 150	−98 ± 115	0.32 ± 0.23
W	Marg.	39 ± 24	36 ± 154	−257 ± 117	−0.13 ± 0.19

Notes. ^aThe values quoted for "V + W" and "Marg." are our best estimates from the WMAP seven-year data. In all cases, the full-resolution temperature maps at HEALPix N_side = 1024 are used. ^bIn all cases, the KQ75y7 mask is used. ^c"Marg." means that the foreground templates (synchrotron, free–free, and dust) have been marginalized over. When the foreground templates are marginalized over, the raw and clean maps yield the same f_NL values.

Download table as:  ASCIITypeset image

As for the point-source bispectrum, we do not detect b_src in V, W, or V + W. In Komatsu et al. (2009a), we estimated that the residual sources could bias $f_{\scriptsize\textit{NL}}^{\rm local}$ by a small positive amount and applied corrections using Monte Carlo simulations. In this paper, we do not attempt to make such corrections, but we note that sources could give $\Delta f_{\scriptsize\textit{NL}}^{\rm local}\sim 2$ (note that the simulations used by Komatsu et al. (2009a) likely overestimated the effect of sources by a factor of two). As the estimator has changed from that used by Komatsu et al. (2009a), extrapolating the previous results is not trivial. Source corrections to $f_{\scriptsize\textit{NL}}^{\rm equil}$ and $f_{\scriptsize\textit{NL}}^{\rm orthog}$ could be larger (Komatsu et al. 2009a), but we have not estimated the magnitude of the effect for the seven-year data.

We used the linear perturbation theory to calculate the angular bispectrum of primordial non-Gaussianity (Komatsu & Spergel 2001). Second-order effects (Pyne & Carroll 1996; Mollerach & Matarrese 1997; Bartolo et al. 2006, 2007; Pitrou 2009, 2010) are expected to give $f_{\scriptsize\textit{NL}}^{\rm local}\sim 1$ (Nitta et al. 2009; Senatore et al. 2009a, 2009b; Khatri & Wandelt 2009, 2010; Boubekeur et al. 2009; Pitrou et al. 2008) and are negligible given the noise level of the WMAP seven-year data.

Among various sources of secondary non-Gaussianities which might contaminate measurements of primordial non-Gaussianity (in particular $f_{\scriptsize\textit{NL}}^{\rm local}$), a coupling between the ISW effect and the weak gravitational lensing is the most dominant source of confusion for $f_{\scriptsize\textit{NL}}^{\rm local}$ (Goldberg & Spergel 1999; Verde & Spergel 2002; Smith & Zaldarriaga 2006; Serra & Cooray 2008; Hanson et al. 2009; Mangilli & Verde 2009). While this contribution is expected to be detectable and bias the measurement of $f_{\scriptsize\textit{NL}}^{\rm local}$ for Planck, it is expected to be negligible for WMAP: using the method of Hanson et al. (2009), we estimate that the expected signal-to-noise ratio of this term in the WMAP seven-year data is about 0.8. We also estimate that this term can give $f_{\scriptsize\textit{NL}}^{\rm local}$ a potential positive bias of $\Delta f_{\scriptsize\textit{NL}}^{\rm local}\sim 2.7$. Calabrese et al. (2010) used the skewness power spectrum method of Munshi et al. (2009) to search for this term in the WMAP five-year data and found a null result. If we subtract $\Delta f_{\scriptsize\textit{NL}}^{\rm local}$ estimated above (for the residual source and the ISW-lensing coupling) from the measured value, $\Delta f_{\scriptsize\textit{NL}}^{\rm local}$ becomes more consistent with zero.

From these results, we conclude that the WMAP seven-year data are consistent with Gaussian primordial fluctuations to within 95% CL. When combined with the limit on $f_{\scriptsize\textit{NL}}^{\rm local}$ from SDSS, $-29<f_{\scriptsize\textit{NL}}^{\rm local}<70$ (95% CL Slosar et al. 2008), we find $-5<f_{\scriptsize\textit{NL}}^{\rm local}<59$ (95% CL).

When CMB photons encounter hot electrons in clusters of galaxies, the temperature of CMB changes due to the inverse Compton scattering by these electrons. This effect, known as the thermal SZ effect (Zel'dovich & Sunyaev 1969; Sunyaev & Zel'dovich 1972), is a source of significant additional (secondary) anisotropies in the microwave sky (see Rephaeli 1995; Birkinshaw 1999; Carlstrom et al. 2002, for reviews).

The temperature change due to the SZ effect in units of thermodynamic temperature, ΔT_SZ, depends on frequency, ν, and is given by (for a spherically symmetric distribution of gas):

$$\begin{equation} \frac{\Delta T_{\rm SZ}(\theta)}{T_{\rm cmb}} = g_\nu \frac{\,\sigma _T}{m_ec^2}\int _{-l_{\rm out}}^{l_{\rm out}} dl P_e\big(\sqrt{l^2+\theta ^2D_A^2}\big), \end{equation} \tag{ 65 }$$

where θ is the angular distance from the center of a cluster of galaxies on the sky, D_A is the proper (not comoving) angular diameter distance to the cluster center, l is the radial coordinates from the cluster center along the line of sight, P_e(r) is the electron pressure profile,  σ_T is the Thomson cross section, m_e is the electron mass, c is the speed of light, and g_ν is the spectral function given by

$$\begin{equation} g_\nu \equiv x\coth \left(\frac{x}{2}\right)-4, \end{equation} \tag{ 66 }$$

where x ≡ hν/(k_BT_cmb) ≃ ν/(56.78GHz) for T_cmb = 2.725 K. In the Rayleigh–Jeans limit, ν → 0, one finds g_ν → −2. At the WMAP frequencies, g_ν = −1.97, −1.94, −1.91, −1.81, and −1.56 at 23, 33, 41, 61, and 94 GHz, respectively. The integration boundary, l_out, will be given later.

The thermal SZ effect (when relativistic corrections are ignored) vanishes at ≃217 GHz. One then finds g_ν>0 at higher frequencies; thus, the thermal SZ effect produces a temperature decrement at ν < 217 GHz, vanishes at 217 GHz, and produces a temperature increment at ν>217 GHz.

The angular power spectrum of temperature anisotropy caused by the SZ effect is sensitive to both the gas distribution in clusters (Atrio-Barandela & Mücket 1999; Komatsu & Kitayama 1999) and the amplitude of matter density fluctuations, i.e.,  σ₈ (Komatsu & Kitayama 1999; Komatsu & Seljak 2002; Bond et al. 2005). While we have not detected the SZ power spectrum in the WMAP data, we have detected the SZ signal from the Coma cluster (Abell 1656) in the one-year (Bennett et al. 2003c) and three-year (Hinshaw et al. 2007) data.

We have also made a statistical detection of the SZ effect by cross-correlating the WMAP data with the locations of known clusters in the X-ray Brightest Abell-type Cluster (XBAC; Ebeling et al. 1996) catalog (Bennett et al. 2003c; Hinshaw et al. 2007). In addition, there have been a number of statistical detections of the SZ effect reported by many groups using various methods (Fosalba et al. 2003; Hernández-Monteagudo & Rubiño-Martín 2004; Hernández-Monteagudo et al. 2004; Myers et al. 2004; Afshordi et al. 2005; Lieu et al. 2006; Bielby & Shanks 2007; Afshordi et al. 2007; Atrio-Barandela et al. 2008; Kashlinsky et al. 2008; Diego & Partridge 2010; Melin et al. 2010).

The Coma cluster (Abell 1656) is a nearby (z = 0.0231) massive cluster located near the north Galactic pole (l, b)=(5675, 8805). The angular diameter distance to Coma, calculated from z = 0.0231 and (Ω_m, Ω_Λ) = (0.277, 0.723), is D_A = 67h⁻¹ Mpc; thus, 10 arcmin on the sky corresponds to the physical distance of 0.195h⁻¹ Mpc at the redshift of Coma.

To extract the SZ signal from the WMAP temperature map, we use the optimal method described in Appendix C: we write down the likelihood function that contains CMB, noise, and the SZ effect, and marginalize it over CMB. From the resulting likelihood function for the SZ effect, which is given by Equation (C7), we find the optimal estimator for the SZ effect in a given angular bin α, $\hat{p}_\alpha$, as

$$\begin{equation} \hat{p}_\alpha = F^{-1}_{\alpha \beta }(t_\beta)_{\nu ^{\prime } p^{\prime }} [N_{\rm pix}+\tilde{C}]^{-1}_{\nu ^{\prime }p^{\prime },\nu p} d_{\nu p}, \end{equation} \tag{ 67 }$$

where the repeated symbols are summed. Here, d_νp is the measured temperature at a pixel p in a frequency band ν, (t_α)_νp is a map of an annulus corresponding to a given angular bin α, which has been convolved with the beam and scaled by the frequency dependence of the SZ effect, $N_{{\rm pix},\nu p,\nu ^{\prime }p^{\prime }}$ is the noise covariance matrix (which is taken to be diagonal in pixel space and ν, i.e., $N_{{\rm pix},\nu p,\nu ^{\prime }p^{\prime }}=\,\sigma ^2_{\nu p}\delta _{\nu \nu ^{\prime }}\delta _{pp^{\prime }}$), and $\tilde{C}_{\nu p,\nu ^{\prime } p^{\prime }}\equiv \sum _{lm}C_lb_{\nu l}b_{\nu ^{\prime } l}Y_{lm,p}Y_{lm,p^{\prime }}^*$ is the signal covariance matrix of CMB convolved with the beam (C_l and b_νl are the CMB power spectrum and the beam transfer function, respectively). A matrix F_αβ gives the 1σ error of $\hat{p}_\alpha$ as $\sqrt{(F^{-1})_{\alpha \alpha }}$, and is given by

$$\begin{equation} F_{\alpha \beta } = (t_\alpha)_{\nu ^{\prime } p^{\prime }} [N_{\rm pix}+\tilde{C}]^{-1}_{\nu ^{\prime }p^{\prime },\nu p} (t_\beta)_{\nu p}. \end{equation} \tag{ 68 }$$

For d_νp, we use the foreground-cleaned V- and W-band temperature maps at the HEALPix resolution of N_side = 1024, masked by the KQ75y7 mask. Note that the KQ75y7 mask includes the seven-year source mask, which removes a potential bias in the reconstructed profile due to any sources which are bright enough to be resolved by WMAP, as well as the sources found by other surveys. Specifically, the seven-year point-source mask includes sources in the seven-year WMAP source catalog (Gold et al. 2011); sources from Stickel et al. (1994); sources with 22 GHz fluxes ⩾0.5 Jy from Hirabayashi et al. (2000); flat spectrum objects from Teräsranta et al. (2001); and sources from the blazar survey of Perlman et al. (1998) and Landt et al. (2001).

In Figure 14, we show the measured angular radial profiles of Coma in 16 angular bins (separated by Δθ = 20 arcmin), in units of the Rayleigh–Jeans temperature, for the V- and W-band data, as well as for the V + W combined data. The error bar at a given angular bin is given by $\sqrt{(F^{-1})_{\alpha \alpha }}$.

Zoom In Zoom Out Reset image size

We find that all of these measurements agree well at θ ⩾ 130 arcmin; however, at smaller angular scales, θ ⩽ 110 arcmin, the V + W result shows less SZ than both the V- and W-only results. Does this make sense? As described in Appendix C, our optimal estimator uses both the C⁻¹-weighted V + W map and the N⁻¹-weighted V − W map. While the latter map vanishes for CMB, it does not vanish for the SZ effect. Therefore, the latter map can be used to separate CMB and SZ effectively.

This explains why the V + W result and the other results are different only at small angular scales: at θ ⩾ 130 arcmin, the measured signal is |ΔT| ≲ 50μK. If this was due to SZ, the difference map, V − W, would give |ΔT| ≲ (1–1.56/1.81) × 50μK ≃ 7μK, which is smaller than the noise level in the difference map, and thus would not show up. In other words, our estimator cannot distinguish between CMB and SZ at θ ⩾ 130 arcmin.

On the other hand, at θ ⩽ 110 arcmin, each of the V- and W-band data shows much bigger signals, |ΔT| ≳ 100 μK. If this was due to SZ, the difference map would give |ΔT| ≳ 14μK, which is comparable to or greater than the noise level in the difference map, and thus would be visible. We find that the difference map does not detect signals in 50 ⩽ θ ⩽ 110 arcmin, which suggests that the measured signal, −100 μK, is not due to SZ, but due to CMB. As a result, the V + W result shows less SZ than the V- and W-only results.

In order to quantify a statistical significance of detection and interpret the result, we model the SZ profile using a spherical β model (Cavaliere & Fusco-Femiano 1976):

$$\begin{equation} \Delta T_{\rm SZ}(\theta) = \Delta T_{\rm SZ}(0) [1+(\theta /\theta _c)^2]^{(1-3\beta)/2}. \end{equation} \tag{ 69 }$$

To make our analysis consistent with previous measurements described later, we fix the core radius, θ_c, and the slope parameter, β, at θ_c = 10.5 arcmin and β = 0.75 (Briel et al. 1992), and vary only the central decrement, ΔT_SZ(0). In this case, the optimal estimator is

$$\begin{equation} \Delta T_{\rm SZ}(0) = \frac{1}{F}t_{\nu ^{\prime } p^{\prime }} [N_{\rm pix}+\tilde{C}]^{-1}_{\nu ^{\prime }p^{\prime },\nu p} d_{\nu p}, \end{equation} \tag{ 70 }$$

where t_νp is a map of the above β model with ΔT_SZ(0) = 1, and

$$\begin{equation} F = t_{\nu ^{\prime } p^{\prime }}[N_{\rm pix}+\tilde{C}]^{-1}_{\nu ^{\prime }p^{\prime },\nu p}t_{\nu p}, \end{equation} \tag{ 71 }$$

gives the 1σ error as $1/\sqrt{F}$.

For V + W, we find

$$$ \Delta T_{\rm SZ,RJ}(0)=-377\pm 105\, \mu {\rm K} \mbox{(68\% CL)}, $$$

which is a 3.6σ measurement of the SZ effect toward Coma. In terms of the Compton y-parameter at the center, we find

$$\begin{eqnarray*} y_{\scriptsize\textit{WMAP}}(0)&=&\frac{-1}{2}\frac{\Delta T_{\rm SZ,RJ}(0)}{T_{\rm cmb}}\\ &=&(6.9\pm 1.9)\times 10^{-5} \mbox{(68\% CL)}. \end{eqnarray*}$$

Let us compare this measurement with the previous measurements. Herbig et al. (1995) used the 5.5 m telescope at the Owens Valley Radio Observatory (OVRO) to observe Coma at 32 GHz. Using the same θ_c and β as above, they found the central decrement of ΔT_SZ,RJ(0) = −505 ± 92μK (68% CL), after subtracting 38μK due to point sources (5C4.81 and 5C4.85). These sources have been masked by our point-source mask, and thus we do not need to correct for point sources.

While our estimate of ΔT_SZ,RJ(0) is different from that of Herbig et al. (1995) only by 1.2σ, and thus is statistically consistent, we note that Herbig et al. (1995) did not correct for the CMB fluctuation in the direction of Coma. As the above results indicate that the CMB fluctuation in the direction of Coma is on the order of −100 μK, it is plausible that the OVRO measurement implies ΔT_SZ,RJ(0) ∼ −400 K, which is an excellent agreement with the WMAP measurement.

The Coma cluster has been observed also by the Millimetre and Infrared Testagrigia Observatory (MITO) experiment (De Petris et al. 2002). Using the same θ_c and β as above, Battistelli et al. (2003) found ΔT_SZ(0) = −184 ± 39, −32 ± 79, and +172 ± 36μK (68% CL) at 143, 214, and 272 GHz, respectively, in units of thermodynamic temperatures. As MITO has three frequencies, they were able to separate SZ, CMB, and the atmospheric fluctuation. By fitting these three data points to the SZ spectrum, ΔT_SZ/T_cmb = g_νy, we find y_MITO(0) = (6.8 ±  1.0 ±  0.7) × 10⁻⁵, which is an excellent agreement with the WMAP measurement. The first error is statistical and the second error is systematic due to 10% calibration error of MITO. The calibration error of the WMAP data (0.2%; Jarosik et al. 2011) is negligible.

Finally, one may try to fit the multi-wavelength data of T_SZ(0) to separate the SZ effect and CMB. For this purpose, we fit the WMAP data in V- and W-band to the β model without correcting for the CMB fluctuation. We find −381 ± 126 μK and −523 ± 127 μK in thermodynamic units (68% CL). The OVRO measurement, T_SZ,RJ(0) = −505 ± 92μK (Herbig et al. 1995), has been scaled to the Rayleigh–Jeans temperature with the SZ spectral dependence correction, and thus we use this measurement at ν = 0. Fitting the WMAP and OVRO data to the SZ effect plus CMB, and the MITO data only to the SZ effect (because the CMB was already removed from MITO using their multi-band data), we find y(0) = (6.8 ± 1.0) × 10⁻⁵ and ΔT_cmb(0) = −136 ± 82 μK (68% CL). This result is consistent with our interpretation that the y-parameter of the center of Coma is 7 × 10⁻⁵ and the CMB fluctuation is on the order of −100 μK.

The analysis presented here shows that our optimal estimator is an excellent tool for extracting the SZ effect from multi-frequency data.

The Coma cluster is the brightest SZ cluster on the sky. There are other clusters that are bright enough to be seen by WMAP.

7.3.1. Sample of Nearby (z < 0.1) Clusters

7.3.2. WMAP Versus X-ray: Cluster-by-cluster Comparison

In Figure 15, we show the measured SZ effect in the symbols with error bars, as well as the expected SZ from the X-ray data in the dashed lines.

Zoom In Zoom Out Reset image size

To compute the expected SZ, we use Equation (65) with P_e = n_ek_BT_e, where n_e and T_e are fits to the X-ray data. Specifically, we use (see Equations (3) and (8) of Vikhlinin et al. 2006)³⁷

$$\begin{eqnarray} n_e^2(r) &=& n_0^2\frac{(r/r_c)^{-\alpha }}{\big(1+r^2/r_c^2\big)^{3\beta -\alpha /2}}\frac{1}{(1+r^\gamma /r_s^\gamma)^{\epsilon /\gamma }}\nonumber\\ && +\frac{n_{02}^2}{\big(1+r^2/r_{c2}^2\big)^{3\beta _2}}, \end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray} \frac{T_e(r)}{T_{mg}}&=&1.35\frac{(x/0.045)^{1.9}+0.45}{(x/0.045)^{1.9}+1}\frac{1}{[1+(x/0.6)^2]^{0.45}}, \end{eqnarray} \tag{ 73 }$$

where x ≡ r/r₅₀₀. The parameters in the above equations are found from the Chandra X-ray data, and kindly made available to us by A. Vikhlinin.

For a given pressure profile, P_e(r), we compute the SZ temperature profile as

$$\begin{eqnarray} \nonumber \Delta T_{\rm SZ}(\theta)&=& g_\nu T_{\rm cmb} \frac{\,\sigma _T}{m_ec^2}P^{\rm 2d}_e(\theta)\\ &\simeq & 273 \mu {\rm K} g_\nu \left[\frac{P^{\rm 2d}_{e}(\theta)}{25 {\rm eV cm^{-3} Mpc}}\right], \end{eqnarray} \tag{ 74 }$$

where P^2d_e(θ) is the projected electron pressure profile on the sky:

$$\begin{eqnarray} &&P^{\rm 2d}_{e}(\theta) = \int _{-\sqrt{r_{\rm out}^2-\theta ^2D_A^2}}^{\sqrt{r_{\rm out}^2-\theta ^2D_A^2}} dl P_e\big(\sqrt{l^2+\theta ^2D_A^2}\big). \end{eqnarray} \tag{ 75 }$$

Here, we truncate the pressure profile at r_out. We take this to be r_out = 6r₅₀₀. While the choice of the boundary is somewhat arbitrary, the results are not sensitive to the exact value because the pressure profile declines fast enough.

We find a good agreement between the measured and expected SZ signals (see Figure 15), except for A754: A754 is a merging cluster with a highly disturbed X-ray morphology, and thus the expected SZ profile, which is derived assuming spherical symmetry (Equation (74)), may be different from the observed one.

To make the comparison quantitative, we select clusters within a given mass bin, and fit the expected SZ profiles to the WMAP data with a single amplitude, a, treated as a free parameter. The optimal estimator for the normalization of pressure, a, is

$$\begin{equation} a = \frac{1}{F}t_{\nu ^{\prime } p^{\prime }} [N_{\rm pix}+\tilde{C}]^{-1}_{\nu ^{\prime }p^{\prime },\nu p} d_{\nu p}, \end{equation} \tag{ 76 }$$

where t_νp is a map containing the predicted SZ profiles around clusters, and the 1σ error is $1/\sqrt{F}$ where F is given by Equation (71).

We summarize the results in the second column of Table 12. We find that the amplitudes of all mass bins are consistent with unity (a = 1) to within 2σ (except for the "non-cooling flow" case, for which a is less than unity at 2.2σ; we shall come back to this important point in the next section). The agreement is especially good for the highest mass bin (M₅₀₀ ⩾ 6 × 10¹⁴ h⁻¹ M_☉), a = 0.90 ± 0.16 (68% CL).

Table 12. Best-fitting Amplitude for the SZ Effect in the WMAP Seven-year data

Mass Rangea	Number of Clusters	Vikhlinin et al.b	Arnaud et al.c
6 ⩽ M500 < 9	5	0.90 ± 0.16	0.73 ± 0.13
4 ⩽ M500 < 6	6	0.73 ± 0.21	0.60 ± 0.17
2 ⩽ M500 < 4	9	0.71 ± 0.31	0.53 ± 0.25
1 ⩽ M500 < 2	9	−0.15 ± 0.55	−0.12 ± 0.47
4 ⩽ M500 < 9	11	0.84 ± 0.13	0.68 ± 0.10
1 ⩽ M500 < 4	18	0.50 ± 0.27	0.39 ± 0.22
4 ⩽ M500 < 9
Cooling flowd	5	1.06 ± 0.18	0.89 ± 0.15
Non-cooling flowe	6	0.61 ± 0.18	0.48 ± 0.15
2 ⩽ M500 < 9	20	0.82 ± 0.12	0.660 ± 0.095
1 ⩽ M500 < 9	29	0.78 ± 0.12	0.629 ± 0.094

Notes. ^aIn units of 10¹⁴ h⁻¹ M_☉. Coma is not included. The masses are derived from the mass–Y_X relation, and are given in the sixth column of Table 2 in Vikhlinin et al. (2009a), times h_vikhlinin = 0.72. ^bDerived from the X-ray data on the individual clusters (Vikhlinin et al. 2009a). ^cThe "universal pressure profile" given by Arnaud et al. (2010). ^dDefinition of "cooling flow" follows that of Vikhlinin et al. (2007). All of cooling-flow clusters here are also "relaxed," according to the criterion of Vikhlinin et al. (2009a). ^eDefinition of "non-cooling flow" follows that of Vikhlinin et al. (2007). All of non-cooling-flow clusters here are also "non-relaxed" (or mergers or morphologically disturbed), according to the criterion of Vikhlinin et al. (2009a).

Download table as:  ASCIITypeset image

Note that this is a 5.6σ detection of the SZ effect, just from stacking five clusters. By stacking 11 clusters with M₅₀₀ ⩾ 4 × 10¹⁴ h⁻¹ M_☉ (i.e., all clusters in Figure 15 but Coma), we find a = 0.84 ± 0.13 (68% CL), a 6.5σ detection. In other words, one does not need to stack many tens or hundreds of clusters to see the SZ effect in the WMAP data, contrary to what is commonly done in the literature (Fosalba et al. 2003; Hernández-Monteagudo & Rubiño-Martín 2004; Hernández-Monteagudo et al. 2004; Myers et al. 2004; Afshordi et al. 2005; Lieu et al. 2006; Bielby & Shanks 2007; Afshordi et al. 2007; Atrio-Barandela et al. 2008; Kashlinsky et al. 2008; Diego & Partridge 2010; Melin et al. 2010).

From this study, we conclude that the WMAP data and the expectation from the X-ray data are in good agreement.

7.3.3. WMAP Versus a "Universal Pressure Profile" of Arnaud et al.: Effect of Recent Mergers

To explore the SZ effect in a large number of clusters, we use a galaxy cluster catalog consisting of the ROSAT-ESO flux-limited X-ray (REFLEX) galaxy cluster survey (Böhringer et al. 2004) in the southern hemisphere above the Galactic plane (δ < 25 and |b|>20°) and the extended Brightest Cluster Sample (eBCS; Ebeling et al. 1998, 2000) in the northern hemisphere above the Galactic plane (δ>0° and |b|>20°). Some clusters are contained in both samples. Eliminating the overlap, this catalog contains 742 clusters of galaxies. Of these, 400, 228, and 114 clusters lie in the redshift ranges of z ⩽ 0.1, 0.1 < z ⩽ 0.2, and 0.2 < z ⩽ 0.45, respectively.

We use the foreground-reduced V- and W-band maps at the HEALPix resolution of N_side = 1024, masked by the KQ75y7 mask, which eliminates the entire Virgo cluster. Note that this mask also includes the point-source mask, which masks sources at the locations of some clusters (such as Coma). After applying the mask, we have 361, 214, and 109 clusters in z ⩽ 0.1, 0.1 < z ⩽ 0.2, and 0.2 < z ⩽ 0.45, respectively.

We again use Equation (67) to find the angular radial profile in four angular bins. For this analysis, (t_α)_νp is a map containing many annuli (one annulus around each cluster) corresponding to a given angular bin α, convolved with the beam and scaled by the frequency dependence of the SZ effect.

We show the measured profile in the top panel of Figure 16. We have done this analysis using three different choices of the maximum redshift, z_max, to select clusters: z_max = 0.1, 0.2, and 0.45. We find that the results are not sensitive to z_max. As expected, the results for z_max = 0.1 have the largest error bars. The error bars for z_max = 0.2 and 0.45 are similar, indicating that we do not gain much more information from z>0.2. The error bars have contributions from instrumental noise and CMB fluctuations. The latter contribution correlates the errors at different angular bins.

Zoom In Zoom Out Reset image size

The top panel shows a decrement of −3.6 ± 1.4 μK at a very large angular distance from the center, θ = 63 arcmin, for z_max = 0.2. As we do not expect to have such an extended gas distribution around clusters, one may wonder if this result implies that we have a bias in the zero level. In order to check for a potential systematic bias, we perform the following null test: instead of measuring the SZ signals from the locations of clusters, we measure them from random locations in the WMAP data. In the middle panel of Figure 16, we show that our method passes a null test. We find that the measured profiles are consistent with zero; thus, our method does not introduce a bias.

Is this signal at a degree scale real? For example, are there nearby massive clusters (such as Coma) which give a significant SZ effect at a degree scale? While the Virgo cluster has the largest angular size on the sky, the KQ75y7 mask eliminates Virgo. In order to see if other nearby clusters give significant contributions, we remove all clusters at z ⩽ 0.03 (where there are 57 clusters) and remeasure the SZ profile. We find that the changes are small, less than 1 μK at all angular bins. At θ = 63 arcmin, the change is especially small, ∼0.1 μK, and thus nearby clusters do not make much contribution to this bin.

The apparent decrement at θ = 63 arcmin is probably due to a statistical fluctuation. The angular bins are correlated with the following correlation matrix:

$$\begin{equation} \left(\begin{array}{@{}cccc@{}} 1 & 0.5242 & 0.0552 & 0.0446 \\ 0.5242 & 1 & 0.4170 & 0.0638 \\ 0.0552 & 0.4170 & 1 & 0.4402 \\ 0.0446 & 0.0638 & 0.4402 & 1 \end{array} \right), \end{equation} \tag{ 77 }$$

where the columns correspond to θ = 7, 35, 63, and 91 arcmin, respectively. The decrements at the first two bins (at θ = 7 and 35 arcmin) can drive the third bin at θ = 63 arcmin to be more negative. Note also that one of the realizations shown in the bottom panel ("Random 1" in the middle panel of Figure 16) shows ∼−3.5 μK at θ = 63 arcmin. The second bin is also negative with a similar amplitude. On the other hand, "Random 2" shows both positive temperatures at the second and third bins, which is also consistent with a positive correlation between these bins.

Finally, in the bottom panel of Figure 16, we compare the measured SZ profile with the expected profiles from various cluster gas models (described in Section 7.5). None of them show a significant signal at θ = 63 arcmin, which is also consistent with our interpretation that it is a statistical fluctuation.

7.5.1. General Idea

In order to interpret the measured SZ profile, we need a model for the electron pressure profile, P_e(r) (see Equation (65)). For fully ionized gas, the electron pressure is related to the gas (baryonic) pressure, P_gas(r), by

$$\begin{equation} P_e(r) = \left(\frac{2+2X}{3+5X}\right)P_{\rm gas}(r), \end{equation} \tag{ 78 }$$

where X is the abundance of hydrogen in clusters. For X = 0.76, one finds P_e(r) = 0.518P_gas(r).

We explore three possibilities: (1) Arnaud et al.'s profile that we have used in Section 7.3, (2) theoretical profiles derived by assuming that the gas pressure is in hydrostatic equilibrium with gravitational potential given by a Navarro–Frenk–White (NFW; Navarro et al. 1997) mass density profile (Komatsu & Seljak 2001), and (3) theoretical profiles from hydrodynamical simulations of clusters of galaxies with and without gas cooling and star formation (Nagai et al. 2007).

Case (2) is relevant because this profile is used in the calculation of the SZ power spectrum (Komatsu & Seljak 2002) that has been used as a template to marginalize over in the cosmological parameter estimation since the three-year analysis (Spergel et al. 2007; Dunkley et al. 2009; Larson et al. 2011). Analytical models and hydrodynamical simulations for the SZ signal are also the basis for planned efforts to use the SZ signal to constrain cosmological models.

As we have shown in the previous section using 29 nearby clusters, Arnaud et al.'s pressure profile overpredicts the SZ effect in the WMAP data by ∼30%. An interesting question is whether this trend extends to a larger number of clusters.

7.5.2. Komatsu–Seljak Profile

The normalization of the KS profile has been fixed by assuming that the gas density at the virial radius is equal to the cosmic mean baryon fraction, Ω_b/Ω_m, times the total mass density at the virial radius. This is an upper limit: for example, star formation turns gas into stars, reducing the amount of gas. KS also assumes that the gas is virialized and in thermal equilibrium (i.e., electrons and protons share the same temperature) everywhere in a cluster, that virialization converts potential energy of the cluster into thermal energy only, and that the pressure contributed by bulk flows, cosmic rays, and magnetic fields are unimportant.

We give details of the gas pressure profiles in Appendix D. In the top left panel of Figure 17, we show Arnaud et al.'s pressure profiles (see Appendix D.1) in the solid lines, and the KS profiles (see Appendix D.2) in the dotted and dashed lines. One of the inputs for the KS profile is the so-called concentration parameter of the NFW profile. The dotted line is for the concentration parameter of c = 10(M_vir/3.42 × 10¹² h⁻¹ M_☉)^−0.2/(1 + z) (Seljak 2000), which was used by Komatsu & Seljak (2002) for their calculation of the SZ power spectrum. Here, M_vir is the virial mass, i.e., mass enclosed within the virial radius. The dashed line is for c = 7.85(M_vir/2 × 10¹² h⁻¹ M_☉)^−0.081/(1 + z)^0.71, which was found from recent N-body simulations with the WMAP five-year cosmological parameters (Duffy et al. 2008).

Zoom In Zoom Out Reset image size

We find that the KS profiles and Arnaud et al.'s profiles generally agree. The agreement is quite good especially for the KS profile with the concentration parameter of Duffy et al. (2008). The KS profiles tend to overestimate the gas pressure relative to Arnaud et al.'s one for low-mass clusters (M_☉ ≲ 10¹⁴ h⁻¹ M_☉). Can we explain this trend by a smaller gas mass fraction in clusters than the cosmic mean? To answer this, we compute the gas mass fraction by integrating the gas density profile:

$$\begin{equation} f_{\rm gas}\equiv \frac{M_{\rm gas,500}}{M_{500}} = \frac{4\pi \int _0^{r_{500}} r^2dr \rho _{\rm gas}(r)}{M_{500}}, \end{equation} \tag{ 79 }$$

where M₅₀₀ and M_gas,500 are the total mass and gas mass contained within r₅₀₀, respectively.

In Figure 18, we show f_gas from X-ray observations (Vikhlinin et al. 2009a):

$$\begin{equation} f_{\rm gas}(h/0.72)^{3/2}=0.125+0.037\log _{10}(M_{500}/10^{15}\, h^{-1}\, M_\odot), \end{equation} \tag{ 80 }$$

for h = 0.7, and f_gas from the KS profiles with the concentration parameters of Seljak (2000) and Duffy et al. (2008). We find that the KS predictions, f_gas ≃ 0.12, are always much smaller than the cosmic mean baryon fraction, Ω_b/Ω_m = 0.167, and are nearly independent of mass. A slight dependence on mass is due to the dependence of the concentration parameters on mass. While the KS profile is normalized such that the gas density at the virial radius is Ω_b/Ω_m times the total mass density, the gas mass within r₅₀₀ is much smaller than Ω_b/Ω_m times M₅₀₀, as the gas density and total matter density profiles are very different near the center: while the gas density profile has a constant-density core, the total matter density, which is dominated by dark matter, increases as ρ_m ∝ 1/r near the center.

Zoom In Zoom Out Reset image size

However, the behavior of f_gas measured from X-ray observations is very different. It has a much steeper dependence on mass than predicted by KS. The reason for such a steep dependence on mass is not yet understood. It could be due to star formation occurring more effectively in lower mass clusters. In any case, for M₅₀₀ = 3 × 10¹⁴ h⁻¹ M_☉, the observed gas mass fraction is f_gas ≃ 0.11, which is only 10% smaller than the KS value, 0.12. For M₅₀₀ = 3 × 10¹³ h⁻¹ M_☉, the observed gas mass fraction, 0.08, is about 30% smaller than the KS value. This is consistent with the difference between the KS and Arnaud et al.'s pressure profiles that we see in Figure 17; thus, once the observed mass dependence of f_gas is taken into account, these profiles agree well.

To calibrate the amplitude of gas pressure, we shall use the KS pressure profile (without any modification to f_gas) as a template, and find its normalization, a, from the WMAP data using the estimator given in Equation (76). We shall present the results for hydrodynamical simulations later.

For a given gas pressure profile, P_gas(r), we compute the electron pressure as P_e = 0.518P_gas (see Equation (78)). We then use Equation (74) to calculate the expected SZ profile, ΔT_SZ(θ). We take the outer boundary of the pressure to be three times the virial radius, r_out = 3r_vir, which is the same as the parameter used by Komatsu & Seljak (2002). In the right panels of Figure 17, we show the predicted ΔT_SZ(θ), which will be used as templates, i.e., t_νp.

7.5.3. Luminosity–Size Relation

Now, in order to compute the expected pressure profiles from each cluster in the catalog, we need to know r₅₀₀. We calculate r₅₀₀ from the observed X-ray luminosity in ROSAT's 0.1–2.4 keV band, L_X, as

$$\begin{eqnarray} \nonumber r_{500} &=& \frac{(0.753\pm 0.063)\, h^{-1}\, {\rm Mpc}}{E(z)}\\ & &\times \left(\frac{L_{\rm X}}{10^{44}\, h^{-2}\, {\rm erg\, s^{-1}}}\right)^{0.228\pm 0.015}, \end{eqnarray} \tag{ 81 }$$

where E(z) ≡ H(z)/H₀ = [Ω_m(1 + z)³ + Ω_Λ]^1/2 for a ΛCDM model. This is an empirical relation found from X-ray observations (see Equation (2) of Böhringer et al. 2007) based upon the temperature–L_X relation from Ikebe et al. (2002) and the r₅₀₀–temperature relation from Arnaud et al. (2005). The error bars have been calculated by propagating the errors in the temperature–L_X and r₅₀₀–temperature relations. Admittedly, there is a significant scatter around this relation, which is the most dominant source of systematic error in this type of analysis. (The results presented in Section 7.3 do not suffer from this systematic error, as they do not rely on L_X–r₅₀₀ relations.) As M₅₀₀ ∝ r³₅₀₀, a ≈10% error in the predicted values of r₅₀₀ gives the mass calibration error of ≈30%. Moreover, the SZ effect is given by M₅₀₀ times the gas temperature, the latter being proportional to M^2/3₅₀₀ according to the virial theorem. Therefore, the total calibration error can be as big as ≈50%.

In order to quantify this systematic error, we repeat our analysis for three different size–luminosity relations: (1) the central values, (2) the normalization and slope shifted up by 1σ to 0.816 and 0.243, and (3) the normalization and slope shifted down by 1σ to 0.690 and 0.213. We adopt this as an estimate for the systematic error in our results due to the size–luminosity calibration error. For how this error would affect our conclusions, see Section 7.7.

Note that this estimate of the systematic error is conservative, as we allowed all clusters to deviate from the best-fit scaling relation at once by ±1σ. In reality, the nature of this error is random, and thus the actual error caused by the scatter in the scaling relation would probably be smaller. Melin et al. (2010) performed such an analysis and found that the systematic error is sub-dominant compared to the statistical error.

Nevertheless, we shall adopt our conservative estimate of the systematic error, as the mean scaling relation also varies from authors to authors. The mean scaling relations used by Melin et al. (2010) are within the error bar of the scaling relation that we use (Equation (81)).

In Figure 19, we show the distribution of M₅₀₀ estimated from clusters in the catalog using the measured values of L_X and Equations (81) and (D2). The distribution peaks at M₅₀₀ ∼ 3 × 10¹⁴ h⁻¹ M_☉ for z_max = 0.2 and 0.45, while it peaks at M₅₀₀ ∼ 1.5 × 10¹⁴ h⁻¹ M_☉ for z_max = 0.1.

Zoom In Zoom Out Reset image size

7.5.4. Results: Arnaud et al.'s Profile

7.5.5. Results: KS Profile and Hydrodynamical Simulation

To see the dependence of the best-fitting normalization on X-ray luminosities (hence M₅₀₀), we divide the cluster samples into three luminosity bins: (1) "High L_X" with 4.5 < L_X/(10⁴⁴ergs^-1) ⩽ 45, (2) "Low L_X" with 0.45 < L_X/(10⁴⁴ergs^-1) ⩽ 4.5, and (3) clusters fainter than (2). There are 82, 417, and 129 clusters in (1), (2), and (3), respectively. In Table 13, we show that we detect significant SZ signals in (1) and (2), despite the smaller number of clusters used in each luminosity bin. We do not have a statistically significant detection in (3).

The high L_X clusters have M ≳ 4 × 10¹⁴ h⁻¹ M_☉. For these clusters, the agreement between the WMAP data and the expected SZ signals is much better. In particular, for the REXCESS scaling relation, we find a = 0.90 ± 0.12, which is consistent with unity within the 1σ statistical error. This implies that, at least for high X-ray luminosity clusters, our results and the results of Melin et al. (2010) agree within the statistical uncertainty.

On the other hand, we find that less luminous clusters tend to have significantly lower best-fitting amplitudes for all models of gas-pressure profiles and scaling relations that we have explored. This trend is consistent with, for example, the gas mass fraction being lower for lower mass clusters. It is also consistent with radio point sources filling some of the SZ effect seen in the WMAP data. For the point-source contamination, see Section 7.7.

The best-fitting amplitudes may be shifted up and down by ≈50% due to the calibration error in the size–luminosity relation (Equation (81)). As we have shown already, the best-fitting amplitudes for the KS profiles can be shifted up to 0.77 and 0.91 for the concentration parameters of Seljak (2000) and Duffy et al. (2008), respectively. Similarly, the amplitude for Arnaud et al.'s profile can be shifted up to 0.97. As this calibration error shifts all amplitudes given in Table 13 by the same amount, it does not affect our conclusion that all of the gas pressure profiles considered above yield similar results.

This type of systematic error can be reduced by using a subset of clusters of galaxies for which the scaling relations are more tightly constrained (see, e.g., Pratt et al. 2009; Vikhlinin et al. 2009a; Mantz et al. 2010a); however, reducing the number of samples increases the statistical error. Indeed, the analysis presented in Section 7.3 does not suffer from the ambiguity in the scaling relations.

How important are radio point sources? While we have not attempted to correct for potential contamination from unresolved radio point sources, we estimate the magnitude of effects here. If, on average, each cluster has an F_src = 10 mJy source, then the corresponding temperatures,

$$\begin{equation} \Delta T_{\rm src} = 40.34 \mu {\rm K} \left[\frac{\sinh ^2(x/2)}{x^4}\frac{F_{\rm src}}{10 {\rm mJy}}\frac{10^{-5} {\rm sr}}{\Omega _{\rm beam}}\right], \end{equation} \tag{ 86 }$$

are 2.24, 2.29, and 2.19 μK in Q, V, and W bands, respectively. Here, x = ν/(56.78GHz), and Ω_beam = 9.0 × 10⁻⁵, 4.2 × 10⁻⁵, and 2.1 × 10⁻⁵ sr are the solid angles of beams in Q, V, and W bands, respectively (Jarosik et al. 2011). Using the radio sources observed in clusters of galaxies by Lin et al. (2009), Diego & Partridge (2010) estimated that the mean flux of sources in Q band is 10.4 mJy, and that at 90 GHz (which is close to 94 GHz of W band) is ≈4–6 mJy. Using these estimates, we expect the source contamination at the level of ≈1–2 μK in V and W bands, which is ≈5%–10% of the measured SZ temperature. Therefore, the best-fitting amplitudes reported in Table 13 could be underestimated by ≈5%–10%.

The gas pressure profile is not the only factor that determines the SZ power spectrum. The other important factor is the mass function, $\textit{dn}/\textit{dM}$:

$$\begin{equation} C_l\propto \int dz\frac{\textit{dV}}{\textit{dz}}\int \textit{dM}\frac{\textit{dn}}{\textit{dM}} \big|\tilde{P}^{\rm 2d}_l\big|^2, \end{equation} \tag{ 87 }$$

where V(z) is the comoving volume of the universe and $\tilde{P}^{\rm 2d}_l$ is the two-dimensional Fourier transform of P^2d(θ). Therefore, a lower-than-expected A_SZ may imply either a lower-than-expected amplitude of matter density fluctuations, i.e.,  σ₈, or a lower-than-expected gas pressure, or both.

As the predictions for the SZ power spectrum available today (see, e.g., Shaw et al. 2009; Sehgal et al. 2010, and references therein) are similar to the prediction of Komatsu & Seljak (2002) (for example, Lueker et al. 2010 found A_SZ = 0.55 ± 0.21 for the prediction of Sehgal et al. 2010, which is based on the gas model of Bode et al. 2009), a plausible explanation for a lower-than-expected A_SZ is a lower-than-expected gas pressure.

Arnaud et al. (2007) find that the X-ray observed integrated pressure enclosed within r₅₀₀, Y_X ≡ M_gas,500T_X, for a given M₅₀₀ is about a factor of 0.75 times the prediction from the Cooling+SF simulation of Nagai et al. (2007). This is in good agreement with our corresponding result for the "High L_X" samples, 0.79 ± 0.10 (68% CL; statistical error only).

While the KS profile is generally in good agreement with Arnaud et al.'s profile, the former is more extended than the latter (see Figure 17), which makes the KS prediction for the projected SZ profiles bigger. Note, however, that the outer slope of the fitting formula given by Arnaud et al. (2010) (Equation (D3)) has been forced to match that from hydrodynamical simulations of Nagai et al. (2007) in r ⩾ r₅₀₀. See the bottom panels of Figure 17. The steepness of the profile at r ≳ r₅₀₀ from the simulation may be attributed to a significant non-thermal pressure support from ρv², which makes it possible to balance gravity by less thermal pressure at larger radii. In other words, the total pressure (i.e., thermal plus ρv²) profile would probably be closer to the KS prediction, but the thermal pressure would decline more rapidly than the total pressure would.

If the SZ effect seen in the WMAP data is less than theoretically expected, what would be the implications? One possibility is that protons and electrons do not share the same temperature. The electron–proton equilibration time is longer than the Hubble time at the virial radius, so that the electron temperature may be lower than the proton temperature in the outer regions of clusters which contribute a significant fraction of the predicted SZ flux (Rudd & Nagai 2009; Wong & Sarazin 2009). The other sources of non-thermal pressure support in outskirts of the cluster (turbulence, magnetic field, and cosmic rays) would reduce the thermal SZ effect relative to the expectation, if these effects are not taken into account in modeling the intracluster medium. Heat conduction may also play some role in suppressing the gas pressure (Loeb 2002, 2007).

In order to explore the impact of gas pressure at r>r₅₀₀, we cut the pressure profile at r_out = r₅₀₀ (instead of 6r₅₀₀) and repeat the analysis. We find a = 0.74 ± 0.09 and 0.44 ± 0.14 for high and low L_X clusters, respectively. (We found a = 0.67 ±  0.09 and 0.43 ± 0.12 for r_out = 6r₅₀₀. See Table 13.) These results are somewhat puzzling—the X-ray observations directly measure gas out to r₅₀₀, and thus we would expect to find a ≈ 1 at least out to r₅₀₀. This result may suggest that, as we have shown in Section 7.3, the problem is not with the outskirts of the cluster, but with the inner parts where the cooling flow has the largest effect.

The relative amplitudes between high and low L_X clusters suggest that a significant amount of pressure is missing in low-mass (M₅₀₀ ≲ 4 × 10¹⁴ h⁻¹ M_☉) clusters, even if we scale all the results such that high-mass clusters are forced to have a = 1. A similar trend is also seen in Figure 3 of Melin et al. (2010). This interpretation is consistent with the SZ power spectrum being lower than theoretically expected. The SPT measures the SZ power spectrum at l ≳ 3000. At such high multipoles, the contributions to the SZ power spectrum are dominated by relatively low-mass clusters, M₅₀₀ ≲ 4 × 10¹⁴ h⁻¹ M_☉ (see Figure 6 of Komatsu & Seljak 2002). Therefore, a plausible explanation for the lower-than-expected SZ power spectrum is a missing pressure (relative to theory) in lower mass clusters.

Scaling relations, gas pressure, and entropy of low-mass clusters and groups have been studied in the literature.⁴⁰Leauthaud et al. (2010) obtained a relation between L_X of 206 X-ray-selected galaxy groups and the mass (M₂₀₀) derived from the stacking analysis of weak lensing measurements. Converting their best-fitting relation to r₂₀₀–L_X relation, we find $r_{200}=\frac{1.26 \,h^{-1}\, {\rm Mpc}}{E^{0.89}(z)} [L_X/(10^{44} \,h^{-2} \;{\rm erg\; s^{-1}})]^{0.22}$. (Note that the pivot luminosity of the original scaling relation is 2.6 × 10⁴² h⁻² erg s^-1.) As r₅₀₀ ≈ 0.65r₂₀₀, their relation is ≈1σ higher than the fiducial scaling relation that we adopted (Equation (81)). Had we used their scaling relation, we would find even lower normalizations.

The next generation of simulations or analytical calculations of the SZ effect should be focused more on understanding the gas pressure profiles, both the amplitude and the shape, especially in low-mass clusters. New measurements of the SZ effect toward many individual clusters with unprecedented sensitivity are now becoming available (Staniszewski et al. 2009; Hincks et al. 2009; Plagge et al. 2010). These new measurements would be important for understanding the gas pressure in low-mass clusters.

In order to calculate the stacked profiles of temperature and polarization of the CMB at the locations of temperature peaks, we need to relate the peak number density contrast, δ_pk, to the underlying temperature fluctuation, ΔT.

One often encounters a similar problem in the large-scale structure of the universe: how can we relate the number density contrast of galaxies to the underlying matter density fluctuation? It is often assumed that the number density contrast of peaks with a given peak height ν is simply proportional to the underlying density field. If one adopted such a linear and scale-independent bias prescription, one would find⁴²

$$\begin{equation} \delta _{\rm pk}(\hat{\mathbf n}) = b_\nu \Delta T(\hat{\mathbf n}). \end{equation} \tag{ B1 }$$

However, our numerical simulations show that the linear bias does not give an accurate description of 〈Q_r〉 or 〈T_r〉. In fact, breakdown of the linear bias is precisely what is expected from the statistics of peaks. From detailed investigations of the statistics of peaks, Desjacques (2008) found the following scale-dependent bias:

$$\begin{equation} \delta _{\rm pk}(\hat{\mathbf n}) = \left[b_\nu - b_\zeta \big(\partial _1^2+\partial _2^2\big)\right] \Delta T(\hat{\mathbf n}). \end{equation} \tag{ B2 }$$

While the first, constant term b_ν has been known for a long time (Kaiser 1984; Bardeen et al. 1986), the second term b_ζ has been recognized only recently. The presence of b_ζ is expected because, to define peaks, one needs to use the information on the first and second derivatives of ΔT. As the first derivative must vanish at the locations of peaks, the above equation does not contain the first derivative.

Desjacques (2008) has derived the explicit forms of b_ν and b_ζ:

$$\begin{eqnarray} &&b_\nu = \frac{1}{\,\sigma _0}\frac{\nu -\gamma \bar{u}}{1-\gamma ^2},\qquad b_\zeta = \frac{1}{\,\sigma _2}\frac{\bar{u}-\gamma \nu }{1-\gamma ^2}, \end{eqnarray} \tag{ B3 }$$

where ν ≡ ΔT/ σ₀, γ ≡  σ²₁/( σ₀ σ₂),  σ_j is the rms of jth derivatives of the temperature fluctuation:

$$\begin{equation} \,\sigma ^2_j = \frac{1}{4\pi }\sum _l (2l+1)[l(l+1)]^jC_l^{\rm TT}\big(W_l^T\big)^2, \end{equation} \tag{ B4 }$$

and W^T_l is the harmonic transform of a window function (which is a combination of the experimental beam, pixel window, and any other additional smoothing applied to the temperature data). The quantity $\bar{u}$ is called the "mean curvature," and is given by $\bar{u}\equiv G_1(\gamma,\gamma \nu)/G_0(\gamma,\gamma \nu)$, where

$$\begin{equation} G_n(\gamma,x_*) \equiv \int _0^\infty dx x^nf(x) \frac{\exp \big[-\!\frac{(x-x_*)^2}{2(1-\gamma ^2)}\big]}{\sqrt{2\pi (1-\gamma ^2)}}. \end{equation} \tag{ B5 }$$

While Desjacques (2008) applied this formalism to a three-dimensional Gaussian random field, it is straightforward to generalize his results to a two-dimensional case, for which f(x) is given by (Bond & Efstathiou 1987),

$$\begin{equation} f(x) = x^2-1+\exp (-x^2). \end{equation} \tag{ B6 }$$

With the bias given by Equation (B2), we find

$$\begin{eqnarray} \nonumber & &\langle \delta _{\rm pk}(\hat{\mathbf {n}})Q_r(\hat{\mathbf {n}}+\hat{ {\btheta }})\rangle = \int \frac{d^2{\mathbf l}}{(2\pi)^2} W_l^TW_l^P(b_\nu +b_\zeta l^2)\\ & &\quad\times \left\lbrace C_l^{\rm TE}\cos [2(\phi -\varphi)] +C_l^{\rm TB}\sin [2(\phi -\varphi)]\right\rbrace e^{i{\mathbf l}\cdot { {\btheta }}},\qquad \end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray} \nonumber & &\langle \delta _{\rm pk}(\hat{\mathbf {n}})U_r(\hat{\mathbf {n}}+\hat{ {\btheta }})\rangle = -\int \frac{d^2{\mathbf l}}{(2\pi)^2} W_l^TW_l^P(b_\nu +b_\zeta l^2)\\ & &\quad\times \left\lbrace C_l^{\rm TE}\sin [2(\phi -\varphi)] -C_l^{\rm TB} \cos [2(\phi -\varphi)]\right\rbrace e^{i{\mathbf l}\cdot { {\btheta }}},\qquad \end{eqnarray} \tag{ B8 }$$

where W^T_l and W^P_l are spherical harmonic transforms of the smoothing functions applied to the temperature and polarization data, respectively. Recalling ${\mathbf l}\cdot { {\btheta }}=l\theta \cos (\phi -\varphi)$, ∫^2π₀dφsin [2(ϕ − φ)]e^{ixcos(ϕ−φ)} = 0, and

$$\begin{eqnarray} J_m(x)&=&\int _\alpha ^{2\pi +\alpha }\frac{d\psi }{2\pi }e^{i(m\psi -x\sin \psi)}, \end{eqnarray} \tag{ B9 }$$

with m = 2, ψ = φ − ϕ − π/2 and α = −ϕ − π/2, we find

$$\begin{eqnarray} & &\langle \delta _{\rm pk}(\hat{\mathbf {n}})Q_r(\hat{\mathbf {n}}+\hat{ {\btheta }})\rangle =-\int \frac{ldl}{2\pi } W_l^TW_l^P(b_\nu +b_\zeta l^2)C_l^{\rm TE}J_2(l\theta),\nonumber\\ \end{eqnarray} \tag{ B10 }$$

$$\begin{eqnarray} & &\langle \delta _{\rm pk}(\hat{\mathbf {n}})U_r(\hat{\mathbf {n}}+\hat{ {\btheta }})\rangle = -\int \frac{ldl}{2\pi } W_l^TW_l^P(b_\nu +b_\zeta l^2)C_l^{\rm TB}J_2(l\theta).\nonumber\\ \end{eqnarray} \tag{ B11 }$$

Using these results in Equations (9) and (10), we finally obtain the desired formulae for the stacked polarization profiles:

$$\begin{eqnarray} \langle Q_r\rangle (\theta) &=& -\int \frac{ldl}{2\pi } W_l^TW_l^P(b_\nu +b_\zeta l^2)C_l^{\rm TE}J_2(l\theta),\quad \end{eqnarray} \tag{ B12 }$$

$$\begin{eqnarray} \langle U_r\rangle (\theta) &=& -\int \frac{ldl}{2\pi } W_l^TW_l^P(b_\nu +b_\zeta l^2)C_l^{\rm TB}J_2(l\theta).\quad \end{eqnarray} \tag{ B13 }$$

Incidentally, the stacked profile of the temperature fluctuation can also be calculated in a similar way:

$$\begin{equation} \langle T\rangle (\theta) = \int \frac{ldl}{2\pi } \big(W_l^T\big)^2(b_\nu +b_\zeta l^2)C_l^{\rm TT}J_0(l\theta). \end{equation} \tag{ B14 }$$

How can we evaluate Equations (B12)–(B14)? One may follow the following steps:

• 1.  
  Compute  σ₀,  σ₁, and  σ₂ from Equation (B4). For example, the WMAP five-year best-fitting temperature power spectrum for a power-law ΛCDM model (Dunkley et al. 2009),⁴³ multiplied by a Gaussian smoothing of 05 full-width-at-half-maximum (FWHM) and the pixel window function for the HEALPix resolution of N_side = 512, gives  σ₀ = 87.9μK,  σ₁ = 1.16 × 10⁴μK, and  σ₂ = 2.89 × 10⁶μK.
• 2.  
  Compute γ =  σ²₁/( σ₀ σ₂). For the above example, we find γ = 0.5306.
• 3.  
  As we need to integrate over peak heights ν, we need to compute the functions G₀(γ, γν) and G₁(γ, γν) for various values of ν. The former function, G₀(γ, γν), can be found analytically (see Equation (A1.9) of Bond & Efstathiou 1987). For G₁, we need to integrate Equation (B5) numerically.
• 4.  
  Compute $\bar{u}=G_1/G_0$. For the above example, we find $\bar{u}=1.596$, 1.831, 3.206, and 5.579 for ν = 0, 1, 5, and 10, respectively.
• 5.  
  Choose a threshold peak height ν_t, and compute the mean surface number density of peaks, $\bar{n}_{\rm pk}$, from Equation (A1.9) of Bond & Efstathiou (1987):

  $$\begin{equation} \bar{n}_{\rm pk}(\nu _t) =\frac{\,\sigma _2^2}{(2\pi)^{3/2}(2\,\sigma _1^2)}\int d\nu e^{-\nu ^2/2}G_0(\gamma,\gamma \nu). \end{equation} \tag{ B15 }$$

  The integration boundary is taken from ν_t to +∞ for temperature hot spots, and from −∞ to −|ν_t| for temperature cold spots. For the above example, we find $4\pi \bar{n}_{\rm pk}=15354.5$, 8741.5, 2348.9, and 247.5 for ν_t = 0, 1, 2, and 3, respectively.
• 6.  
  Compute b_ν and b_ζ from Equation (B3) for various values of ν.
• 7.  
  Average b_ν and b_ζ over ν. We calculate the averaged bias parameters, $\bar{b}_\nu$ and $\bar{b}_\zeta$, by integrating b_ν and b_ζ multiplied by the number density of peaks for a given ν. We then divide the integral by the mean number density of peaks, $\bar{n}_{\rm pk}$, to find

  $$\begin{eqnarray} \bar{b}_\nu &=& \frac{1}{\bar{n}_{\rm pk}(\nu _t)} \frac{\,\sigma _2^2}{(2\pi)^{3/2}\big(2\,\sigma _1^2\big)}\int d\nu e^{-\nu ^2/2}G_0(\gamma,\gamma \nu)b_\nu,\nonumber\\ \end{eqnarray} \tag{ B16 }$$

  $$\begin{eqnarray} \bar{b}_\zeta &=& \frac{1}{\bar{n}_{\rm pk}(\nu _t)} \frac{\,\sigma _2^2}{(2\pi)^{3/2}\big(2\,\sigma _1^2\big)}\int d\nu e^{-\nu ^2/2}G_0(\gamma,\gamma \nu)b_\zeta.\nonumber\\ \end{eqnarray} \tag{ B17 }$$

  The integration boundary is taken from ν_t to +∞ for temperature hot spots, and from −∞ to −|ν_t| for temperature cold spots. For the above example, we find $(\bar{b}_\nu,\bar{b}_\zeta)= (3.289\times 10^{-3}$, 6.039 × 10⁻⁷), (1.018 × 10⁻², 5.393 × 10⁻⁷), (2.006 × 10⁻², 4.569 × 10⁻⁷), and (3.128 × 10⁻², 3.772 × 10⁻⁷) for ν_t = 0, 1, 2, and 3, respectively (all in units of μK⁻¹). The larger the peak height is, the larger the linear bias and the smaller the scale-dependent bias becomes.
• 8.  
  Use $\bar{b}_\nu$ and $\bar{b}_\zeta$ in Equations (B12) and (B13) to compute 〈Q_r〉(θ) and 〈U_r〉(θ) for a given set of C^TE_l and C^TB_l, respectively.