FIVE-YEAR WILKINSON MICROWAVE ANISOTROPY PROBE^* OBSERVATIONS: COSMOLOGICAL INTERPRETATION

The shape of the power spectrum of primordial curvature perturbations, $P_{\cal R}(k)$, is one of the most powerful and practical tool for distinguishing among inflation models. Inflation models with featureless scalar-field potentials usually predict that $P_{\cal R}(k)$ is nearly a power law (Kosowsky & Turner 1995)

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

Here, $\Delta ^2_{\cal R}(k)$ is a useful quantity, which gives an approximate contribution of ${\cal R}$ at a given scale per logarithmic interval in k to the total variance of ${\cal R}$, as $\langle {\cal R}^2({\mathbf{x}})\rangle = \int d\ln k \Delta ^2_{\cal R}(k)$. It is clear that the special case with n_s = 1 and dn_s/dln k = 0 results in the "scale-invariant spectrum," in which the contributions of ${\cal R}$ at any scales per logarithmic interval in k to the total variance are equal (hence, the term "scale invariance"). Following the usual terminology, we shall call n_s and dn_s/dln k the tilt of the spectrum and the running index, respectively. We shall take k₀ to be 0.002 Mpc⁻¹.

The significance of n_s and dn_s/dln k is that different inflation models motivated by different physics make specific, testable predictions for the values of n_s and dn_s/dln k. For a given shape of the scalar field potential, V(ϕ), of a single-field model, for instance, one finds that n_s is given by a combination of the first derivative and second derivative of the potential, 1 − n_s = 3M²_pl(V'/V)² − 2M²_pl(V''/V) (where M²_pl = 1/(8πG) is the reduced Planck mass), and dn_s/dln k is given by a combination of V'/V, V''/V, and V‴/V (see Liddle & Lyth 2000, for a review).

This means that one can reconstruct the shape of V(ϕ) up to the first three derivatives in this way. As the expansion rate squared is proportional to V(ϕ) via the Friedmann equation, H² = V/(3M²_pl), one can reconstruct the expansion history during inflation by measuring the shape of the primordial power spectrum.

How generic are n_s and dn_s/dln k? They are physically motivated by the fact that most inflation models satisfy the slow-roll conditions, and thus deviations from a pure power-law, scale-invariant spectrum, n_s = 1 and dn_s/dln k = 0, are expected to be small, and the higher-order derivative terms such as V'''' and higher are negligible. However, there is always a danger of missing some important effects, such as sharp features, when one relies too much on a simple parametrization like this. Therefore, a number of people have investigated various, more general ways of reconstructing the shape of $P_{\cal R}(k)$ (Matsumiya et al. 2002, 2003; Mukherjee & Wang 2003b, 2003a, 2003c; Bridle et al. 2003; Kogo et al. 2004, 2005; Hu & Okamoto 2004; Hannestad 2004; Shafieloo & Souradeep 2004, 2008; Sealfon et al. 2005; Tocchini-Valentini et al. 2005, 2006; Spergel et al. 2007; Verde & Peiris 2008) and V(ϕ) (Lidsey et al. 1997; Grivell & Liddle 2000; Kadota et al. 2005; Covi et al. 2006; Lesgourgues & Valkenburg 2007; Powell & Kinney 2007).

These studies have indicated that the parametrized form (Equation (15)) is basically a good fit, and no significant features were detected. Therefore, we do not repeat this type of analysis in this paper, but focus on constraining the parametrized form given by Equation (15).

Finally, let us comment on the choice of priors. We impose uniform priors on n_s and dn_s/dln k, but there are other possibilities for the choice of priors. For example, one may impose uniform priors on the slow-roll parameters, = (M²_pl/2)(V'/V)², η = M²_pl(V''/V) and ξ = M⁴_pl(V'V‴/V²), and on the number of e-foldings, N, rather than on n_s and dn_s/dln k (Peiris & Easther 2006a, 2006b; Easther & Peiris 2006). It has been found that, as long as one imposes a reasonable lower bound on N, N > 30, both approaches yield similar results.

To constrain n_s and dn_s/dln k, we shall use the WMAP 5-year temperature and polarization data, the small-scale CMB data, and/or BAO and SN distance measurements. In Table 4, we summarize our results presented in Sections 3.1.2, 3.1.3, and 3.2.4.

First, we test the inflation models with dn_s/dln k = 0 and negligible gravitational waves. The WMAP 5-year temperature and polarization data yield n_s = 0.963^+0.014_−0.015, which is slightly above the 3-year value with a smaller uncertainty, n_s(3yr) = 0.958 ± 0.016 (Spergel et al. 2007). We shall provide the reason for this small upward shift in Section 3.1.3.

The scale-invariant, Harrison–Zel'dovich–Peebles spectrum, n_s = 1, is at 2.5 standard deviations away from the mean of the likelihood for the WMAP-only analysis. The significance increases to 3.1 standard deviations for WMAP+BAO+SN. Looking at the two-dimensional constraints that include n_s, we find that the most dominant correlation that is still left is the correlation between n_s and Ω_bh² (see Figure 1). The larger the Ω_bh² is, the smaller the second peak becomes, and the larger the n_s is required to compensate it. Also, the larger the Ω_bh² is, the larger the Silk damping (diffusion damping) becomes, and the larger the n_s is required to compensate it.

This argument suggests that the constraint on n_s should improve as we add more information from the small-scale CMB measurements that probe the Silk damping scales; however, the current data do not improve the constraint very much yet: n_s = 0.960 ± 0.014 from WMAP+CMB, where "CMB" includes the small-scale CMB measurements from CBI (Mason et al. 2003; Sievers et al. 2003, 2007; Pearson et al. 2003; Readhead et al. 2004), VSA (Dickinson et al. 2004), ACBAR (Kuo et al. 2004, 2007), and BOOMERanG (Ruhl et al. 2003; Montroy et al. 2006; Piacentini et al. 2006), all of which go well beyond the WMAP angular resolution, so that their small-scale data are statistically independent of the WMAP data.

We find that the small-scale CMB data do not improve the determination of n_s because of their relatively large statistical errors. We also find that the calibration and beam errors are important. Let us examine this using the latest ACBAR data (Reichardt et al. 2008). The WMAP+ACBAR yields 0.964 ± 0.014. When the beam error of ACBAR is ignored, we find n_s = 0.964 ± 0.013. When the calibration error is ignored, we find n_s = 0.962 ± 0.013. Therefore, both the beam and calibration error are important in the error budget of the ACBAR data.

The Big Bang nucleosynthesis (BBN), combined with measurements of the deuterium-to-hydrogen ratio, D/H, from quasar absorption systems, has been extensively used for determining Ω_bh², independent of any other cosmological parameters (see Steigman 2007, for a recent summary). The precision of the latest determination of Ω_bh² from BBN (Pettini et al. 2008) is comparable to that of the WMAP data-only analysis. More precise measurements of D/H will help reduce the correlation between n_s and Ω_bh², yielding a better determination of n_s.

There is still a bit of correlation left between n_s and the electron-scattering optical depth, τ (see Figure 1). While a better measurement of τ from future WMAP observations as well as the Planck satellite should help reduce the uncertainty in n_s via a better measurement of τ, the effect of Ω_bh² is much larger.

We find that the other datasets, such as BAO, SN, and the shape of galaxy power spectrum from SDSS or 2dFGRS, do not improve our constraints on n_s, as these datasets are not sensitive to Ω_bh² or τ; however, this will change when we include the running index, dn_s/dln k (Section 3.1.3) and/or the tensor-to-scalar ratio, r (Section 3.2.4).

Next, we explore more general models in which a sizable running index may be generated (we still do not include gravitational waves; see Section 3.2 for the analysis that includes gravitational waves). We find no evidence for dn_s/dln k from WMAP only, −0.090 < dn_s/dln k < 0.019(95%CL), or WMAP+BAO+SN −0.068 < dn_s/dln k < 0.012(95%CL). The improvement from WMAP-only to WMAP+BAO+SN is only modest.

We find a slight upward shift from the 3-year result, dn_s/dln k = −0.055^+0.030_−0.031 (68% CL; Spergel et al. 2003), to the 5-year result, dn_s/dln k = −0.037 ± 0.028 (68% CL; WMAP only). This is caused by a combination of three effects:

• 1.  
  The 3-year number for dn_s/dln k was derived from an older analysis pipeline for the temperature data, namely the resolution 3 (instead of 4) pixel-based low-l temperature likelihood and a higher point-source amplitude, A_ps = 0.017 μK²sr.
• 2.  
  With 2 years of more integration, we have a better signal-to-noise near the third acoustic peak, whose amplitude is slightly higher than the 3-year determination (Nolta et al. 2009).
• 3.  
  With the improved beam model (Hill et al. 2009), the temperature power spectrum at l ≳ 200 has been raised nearly uniformly by ∼2.5% (Hill et al. 2009; Nolta et al. 2009).

All of these effects go in the same direction, that is, to increase the power at high multipoles and reduce a negative running index. We find that these factors contribute to the upward shift in dn_s/dln k at comparable levels.

Note that an upward shift in n_s for a power-law model, 0.958 to 0.963 (Section 3.1.2), is not subject to (1) because the 3-year number for n_s was derived from an updated analysis pipeline using the resolution 4 low-l code and A_ps = 0.014 μK² sr. We find that (2) and (3) contribute to the shift in n_s at comparable levels. An upward shift in σ₈ from the 3-year value, 0.761, to the 5-year value, 0.796, can be similarly explained.

We do not find any significant evidence for the running index when the WMAP data and small-scale CMB data (CBI, VSA, ACBAR07, BOOMERanG) are combined, −0.1002 < dn_s/dln k < −0.0037(95%CL), or the WMAP data and the latest results from the analysis of the complete ACBAR data (Reichardt et al. 2008) are combined, −0.082 < dn_s/dln k < 0.015(95%CL).

Our best 68% CL constraint from WMAP+BAO+SN shows no evidence for the running index, dn_s/dln k = −0.028 ± 0.020. In order to improve upon the limit on dn_s/dln k, one needs to determine n_s at small scales, as dn_s/dln k is simply given by the difference between n_s's measured at two different scales, divided by the logarithmic separation between two scales. The Lyα forest provides such information (see Section 7; also Table 4).

The presence of primordial gravitational waves is a robust prediction of inflation models, as the same mechanism that generated primordial density fluctuations should also generate primordial gravitational waves (Grishchuk 1975; Starobinsky 1979). The amplitude of gravitational waves relative to that of density fluctuations is model-dependent.

The primordial gravitational waves leave their signatures imprinted on the CMB temperature anisotropy (Rubakov et al. 1982; Fabbri & Pollock 1983; Abbott & Wise 1984; Starobinsky 1985; Crittenden et al. 1993), as well as on polarization (Basko & Polnarev 1980; Bond & Efstathiou 1984; Polnarev 1985; Crittenden et al. 1993, 1995; Coulson et al. 1994).²⁰ The spin-2 nature of gravitational waves leads to two types of a polarization pattern on the sky (Zaldarriaga & Seljak 1997; Kamionkowski et al. 1997a): (1) the curl-free mode (E mode), in which the polarization directions are either purely radial or purely tangential to hot/cold spots in temperature, and (2) the divergence-free mode (B mode), in which the pattern formed by polarization directions around hot/cold spots possess nonzero vorticity.

In the usual gravitational instability picture, in the linear regime (before shell crossing of fluid elements), velocity perturbations can be written in terms of solely a gradient of a scalar velocity potential u, $\vec{v}=\vec{\nabla }u$. This means that no vorticity would arise, and, therefore, no B mode polarization can be generated from density or velocity perturbations in the linear regime. However, primordial gravitational waves can generate both E and B mode polarization; thus, the B mode polarization offers a smoking-gun signature for the presence of primordial gravitational waves (Seljak & Zaldarriaga 1997; Kamionkowski et al. 1997a). This gives us a strong motivation to look for signatures of the primordial gravitational waves in CMB.

In Table 4, we summarize our constraints on the amplitude of gravitational waves, expressed in terms of the tensor-to-scalar ratio, r, defined by Equation (20).

We quantify the amplitude and spectrum of primordial gravitational waves in the following form:

$$\begin{equation} \Delta ^2_h(k) \equiv \frac{k^3P_h(k)}{2\pi ^2}= \Delta ^2_h(k_0) \left(\frac{k}{k_0}\right)^{n_t}, \end{equation} \tag{ 16 }$$

where we have ignored a possible scale dependence of n_t(k), as the current data cannot constrain it. Here, by P_h(k), we mean

$$\begin{equation} \langle \tilde{h}_{ij}({\mathbf{k}})\tilde{h}^{ij}({\mathbf{k}^{\prime }})\rangle = (2\pi)^3P_h(k)\delta ^3({\mathbf{k}}-{\mathbf{k}}^{\prime }), \end{equation} \tag{ 17 }$$

where $\tilde{h}_{ij}({\bf k})$ is the Fourier transform of the tensor metric perturbations, g_ij = a²(δ_ij + h_ij), which can be further decomposed into the two polarization states (+ and ×) with the appropriate polarization tensor, e^(+,×)_ij, as

$$\begin{equation} \tilde{h}_{ij}({\mathbf{k}}) = \tilde{h}_+({\mathbf{k}})e_{ij}^+({\mathbf{k}}) + \tilde{h}_\times ({\mathbf{k}})e_{ij}^{\times }({\mathbf{k}}), \end{equation} \tag{ 18 }$$

with the normalization that e⁺_ije^+,ij = e^×_ije^×,ij = 2 and e⁺_ije^×,ij = 0. Unless there was a parity-violating interaction term such as $f(\phi) R_{\mu \nu \rho \sigma }\tilde{R}^{\mu \nu \rho \sigma }$, where f(ϕ) is an arbitrary function of a scalar field, R_μνρσ is the Riemann tensor, and $\tilde{R}^{\mu \nu \rho \sigma }$ is a dual tensor (Lue et al. 1999), both polarization states are statistically independent with the same amplitude, meaning

$$\begin{equation} \langle |\tilde{h}_+|^2\rangle = \langle |\tilde{h}_\times |^2\rangle \equiv \langle |\tilde{h}^2|\rangle,\qquad \langle \tilde{h}_\times \tilde{h}_+^*\rangle =0. \end{equation} \tag{ 19 }$$

This implies that parity-violating correlations, such as the TB and EB correlations, must vanish. We shall explore such parity-violating correlations in Section 4 in a slightly different context. For limits on the difference between $\langle |\tilde{h}_+|^2\rangle$ and $\langle |\tilde{h}_\times |^2\rangle$ from, respectively, the TB and EB spectra of the WMAP 3-year data, see Saito et al. (2007).

In any case, this definition suggests that P_h(k) is given by $P_h(k) = 4\langle |\tilde{h}|^2\rangle$. Notice a factor of 4. This is the definition of P_h(k) that we have been consistently using since the first year release (Peiris et al. 2003; Spergel et al. 2003, 2007; Page et al. 2007). We continue to use this definition.

With this definition of Δ²_h(k) (Equation (16)), we define the tensor-to-scalar ratio, r, at k = k₀, as

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

where $\Delta _{\cal R}(k)$ is the curvature perturbation spectrum given by Equation (15). We shall take k₀ to be 0.002 Mpc⁻¹. In this paper, we sometimes loosely call this quantity the "amplitude of gravitational waves."

What about the tensor spectral tilt, n_t? In single-field inflation models, there exists the so-called consistency relation between r and n_t (see Liddle & Lyth 2000, for a review)

$$\begin{equation} n_t = -\frac{r}{8}. \end{equation} \tag{ 21 }$$

In order to reduce the number of parameters, we continue to impose this relation at k = k₀ = 0.002 Mpc⁻¹. For a discussion on how to impose this constraint in a more self-consistent way, see Peiris & Easther (2006a).

To constrain r, we shall use the WMAP 5-year temperature and polarization data, the small-scale CMB data, and/or BAO and SN distance measurements.

Let us show how the gravitational wave contribution is constrained by the WMAP data (see Figure 3). In this pedagogical analysis, we vary only r and τ, while adjusting the overall amplitude of fluctuations such that the height of the first peak of the temperature power spectrum is always held fixed. All the other cosmological parameters are fixed at Ω_k = 0, Ω_bh² = 0.02265, Ω_ch² = 0.1143, H₀ = 70.1 km s^-1 Mpc⁻¹, and n_s = 0.960. Note that the limit on r from this analysis should not be taken as our best limit, as this analysis ignores the correlation between r and the other cosmological parameters, especially n_s. The limit on r from the full analysis will be given in Section 3.2.4.

• 1.  
  (The gray contours in the left panel and the upper right of the right panel of Figure 3.) The low-l polarization data (TE/EE/BB) at l ≲ 10 are unable to place meaningful limits on r. A large r can be compensated by a small τ, producing nearly the same EE power spectrum at l ≲ 10. (Recall that the gravitational waves also contribute to EE.) As a result, r that is as large as 10 is still allowed within 68% CL.²¹
• 2.  
  (The red contours in the left panel and the lower left of the right panel of Figure 3.) Such a high value of r, however, produces too negative a TE correlation between 30 ≲ l ≲ 150. Therefore, we can improve the limit on r significantly—by nearly an order of magnitude—by simply using the high-l TE data. The 95% upper bound at this point is still as large as r ∼ 2.²²
• 3.  
  (The blue contours in the left panel and the upper left of the right panel of Figure 3.) Finally, the low-l temperature data at l ≲ 30 severely limit the excess low-l power due to gravitational waves, bringing the upper bound down to r ∼ 0.2. Note that this bound is about a half of what we actually obtain from the full MCMC of the WMAP-only analysis, r < 0.43(95%CL). This is because we have fixed n_s, and thus ignored the correlation between r and n_s shown by Figure 2.

Zoom In Zoom Out Reset image size

Having understood which parts of the temperature and polarization spectra constrain r, we obtain the upper limit on r from the full exploration of the likelihood space using the MCMC. Figure 2 shows the one-dimensional constraint on r and the two-dimensional constraint on r and n_s, assuming a negligible running index. With the WMAP 5-year data alone, we find r < 0.43(95%CL). Since the B-mode contributes little here, and most of the information essentially comes from TT and TE, our limit on r is highly degenerate with n_s, and thus we can obtain a better limit on r only when we have a better limit on n_s.

When we add BAO and SN data, the limit improves significantly to r < 0.22(95%CL). This is because the distance information from BAO and SN reduces the uncertainty in the matter density, and thus it helps to determine n_s better because n_s is also degenerate with the matter density. This "chain of correlations" helped us improve our limit on r significantly from the previous results. This limit, r < 0.22(95%CL), is the best limit on r to date.²³ With the new data, we are able to get more stringent limits than our earlier analyses (Spergel et al. 2007) that combined the WMAP data with measurements of the galaxy power spectrum and found r < 0.30 (95% CL).

A dramatic reduction in the uncertainty in r has an important implication for n_s as well. Previously, n_s > 1 was within the 95% CL when the gravitational wave contribution was allowed, owing to the correlation between n_s and r. Now, we are beginning to disfavor n_s > 1 even when r is nonzero: with WMAP+BAO+SN, we find −0.0022 < 1 − n_s < 0.0589(95%CL).²⁴

However, these stringent limits on r and n_s weaken to −0.246 < 1 − n_s < 0.034(95%CL) and r < 0.55(95%CL) when a sizable running index is allowed. The BAO and SN data helped reduce the uncertainty in dn_s/dln k (Figure 4), but not enough to improve on the other parameters compared to the WMAP-only constraints. The Lyα forest data can improve the limit on dn_s/dln k even when r is present (see Section 7; also Table 4).

Zoom In Zoom Out Reset image size

The flatness of the observable universe is one of the predictions of conventional inflation models. How much curvature can we expect from inflation? The common view is that inflation naturally produces the spatial curvature parameter, Ω_k, on the order of the magnitude of quantum fluctuations, that is, Ω_k ∼ 10⁻⁵. However, the current limit on Ω_k is of order 10⁻²; thus, the current data are not capable of reaching the level of Ω_k that is predicted by the common view.

Would a detection of Ω_k rule out inflation? It is possible that the value of Ω_k is just below our current detection limit, even within the context of inflation: inflation may not have lasted for so long, and the curvature radius of our universe may just be large enough for us not to see the evidence for curvature within our measurement accuracy, yet. While this sounds like fine-tuning, it is a possibility.

This is something we can test by constraining Ω_k better. There is also a revived (and growing) interest in measurements of Ω_k, as Ω_k is degenerate with the equation of state of dark energy, w. Therefore, a better determination of Ω_k has an important implication for our ability to constrain the nature of dark energy.

Measurements of the CMB power spectrum alone do not strongly constrain Ω_k. More precisely, any experiments that measure the angular diameter or luminosity distance to a single redshift are not able to constrain Ω_k uniquely, as the distance depends not only on Ω_k, but also on the expansion history of the universe. For a universe containing matter and vacuum energy, it is essential to combine at least two absolute distance indicators, or the expansion rates, out to different redshifts, in order to constrain the spatial curvature well. Note that CMB is also sensitive to Ω_Λ, via the late-time integrated Sachs–Wolfe (ISW) effect, and to Ω_m, via the signatures of gravitational lensing in the CMB power spectrum. These properties can be used to reduce the correlation between Ω_k and Ω_m (Stompor & Efstathiou 1999) or Ω_Λ (Ho et al. 2008; Giannantonio et al. 2008).

It has been pointed out by a number of people (e.g., Eisenstein et al. 2005) that a combination of distance measurements from BAO and CMB is a powerful way to constrain Ω_k. One needs more distances, if dark energy is not a constant but dynamical.

In this section, we shall make one important assumption that the dark energy component is vacuum energy, that is, a cosmological constant. In Section 5, we shall study the case in which the equation of state, w, and Ω_k are varied simultaneously.

Figure 6 shows the limits on Ω_Λ and Ω_k. While the WMAP data alone cannot constrain Ω_k (see the left panel), the WMAP data combined with the HST's constraint on H₀ tighten the constraint significantly, to −0.052 < Ω_k < 0.013(95%CL). The WMAP data combined with SN yield ∼50% better limits, −0.0316 < Ω_k < 0.0078(95%CL), compared to WMAP+HST. Finally, the WMAP+BAO yields the smallest statistical uncertainty, −0.0170 < Ω_k < 0.0068(95%CL), which is a factor of 2.6 and 1.7 better than WMAP+HST and WMAP+SN, respectively. This shows how powerful the BAO is in terms of constraining the spatial curvature of the universe; however, this statement needs to be re-evaluated when dynamical dark energy is considered, for example, w ≠ −1. We shall come back to this point in Section 5.

Zoom In Zoom Out Reset image size

Finally, when WMAP, BAO, and SN are combined, we find −0.0178 < Ω_k < 0.0066(95%CL). As one can see from the right panel of Figure 6, the constraint on Ω_k is totally dominated by that from WMAP+BAO; thus, the size of the uncertainty does not change very much from WMAP+BAO to WMAP+BAO+SN. Note that the above result indicates that we have reached 1.3% accuracy (95% CL) in determining Ω_k, which is rather good. The future BAO surveys at z ∼ 3 are expected to yield an order of magnitude better determination, that is, 0.1% level, of Ω_k (Knox 2006).

It is instructive to convert our limit on Ω_k to the limits on the curvature radius of the universe. As Ω_k is defined as Ω_k = −kc²/(H²₀R²_curv), where R_curv is the present-day curvature radius, one can convert the upper bounds on Ω_k into the lower bounds on R_curv, as $R_{\rm curv}=(c/H_0)/\sqrt{\left|\Omega _k\right|}=3/\sqrt{\left|\Omega _k\right|}\,h^{-1}\,{\rm Gpc}$. For negatively curved universes, we find R_curv > 37 h⁻¹ Gpc, whereas for positively curved universes, R_curv > 22 h⁻¹ Gpc. Incidentally these values are greater than the particle horizon at present, 9.7 h⁻¹ Gpc (computed for the same model).

The 68% limits from the 3-year data (Spergel et al. 2007) were Ω_k = −0.012 ± 0.010 from WMAP-3 yr+BAO (where BAO is from the SDSS LRG of Eisenstein et al. (2005)), and Ω_k = −0.011 ± 0.011 from WMAP-3 yr+SN (where SN is from the SNLS data of Astier et al. 2006). The 68% limit from WMAP-5 yr+BAO+SN (where both BAO and SN have more data than for the 3-year analysis) is Ω_k = −0.0050^+0.0061_−0.0060. A significant improvement in the constraint is due to a combination of the better WMAP, BAO, and SN data.

We conclude that, if dark energy is vacuum energy (cosmological constant) with w = −1, we do not find any deviation from a spatially flat universe.

What does this imply for inflation? Since we do not detect any finite curvature radius, inflation had to last for a long enough period in order to make the observable universe sufficiently flat within the observational limits. This argument allows us to find a lower bound on the total number of e-foldings of the expansion factor during inflation, from the beginning to the end, N_tot ≡ ln(a_end/a_begin) (see also Section 4.1 of Weinberg 2008).

When the curvature parameter, Ω_k, is much smaller than unity, it evolves with the expansion factor, a, as Ω_k ∝ a⁻², a², and a during inflation, radiation, and matter era, respectively. Therefore, the observed Ω_k is related to Ω_k at the beginning of inflation as³⁰

$$\begin{equation} \frac{\Omega _k^{\rm obs}}{\Omega _k^{\rm begin}} = \left(\frac{a_{\rm today}}{a_{\rm eq}}\right)\left(\frac{a_{\rm eq}}{a_{\rm end}}\right)^2\left(\frac{a_{\rm begin}}{a_{\rm end}}\right)^2 \end{equation} \tag{ 28 }$$

$$\begin{equation} = (1+z_{\rm eq})\left(\frac{T_{\rm end}g_{*,\rm end}^{1/3}}{T_{\rm eq}g_{*,\rm eq}^{1/3}}\right)^2e^{-2N_{\rm tot}}, \end{equation} \tag{ 29 }$$

where g_* is the effective number of relativistic dof contributing to entropy, z_eq is the matter–radiation equality redshift, and T_end and T_eq are the reheating temperature of the universe at the end of inflation³¹ and the temperature at the equality epoch, respectively. To within 10% accuracy, we take 1 + z_eq = 3200, T_eq = 0.75 eV, and g_*,eq = 3.9. We find

$$\begin{equation} N_{\rm tot} = 47 -\frac{1}{2}\ln \frac{\Omega _k^{\rm obs}/0.01}{\Omega _k^{\rm begin}} +\ln \frac{T_{\rm end}}{10^8 {\rm GeV}}+\frac{1}{3}\ln \frac{g_{*,\rm end}}{200}. \end{equation} \tag{ 30 }$$

Here, it is plausible that g_*,end ∼ 100 in the standard model of elementary particles, and ∼200 when the supersymmetric partners are included. The difference between these two cases gives the error of only ΔN_tot = −0.2, and thus can be ignored.

The curvature parameter at the beginning of inflation must be below of order unity, as inflation would not begin otherwise. However, it is plausible that Ω^begin_k was not too much smaller than 1; otherwise, we have to explain why it was so small before inflation, and probably we would have to explain it by inflation before inflation. In that case, N_tot would refer to the sum of the number of e-foldings from two periods of inflation. From this argument, we shall take Ω^begin_k ∼ 1.

The reheating temperature can be anywhere between 1 MeV and 10¹⁶ GeV. It is more likely that it is between 1 TeV and 10⁸ GeV for various reasons, but the allowed region is still large. If we scale the result to a reasonably conservative lower limit on the reheating temperature, T_end ∼ 1 TeV, then we find, from our limit on the curvature of the universe,

$$\begin{equation} N_{\rm tot}>36 + \ln \frac{T_{\rm end}}{1 {\rm TeV}}. \end{equation} \tag{ 31 }$$

A factor of 10 improvement in the upper limit on |Ω^begin_k| will raise this limit by ΔN_tot = 1.2.

Again, N_tot here refers to the total number of e-foldings of inflation. In Section 3.3, we use N ≡ ln(a_end/a_WMAP), which is the number of e-foldings between the end of inflation and the epoch when the wavelength of fluctuations that we probe with WMAP leaves the horizon during inflation. Therefore, by definition, N is less than N_tot.

In the simplest model of inflation, the distribution of primordial fluctuations is close to a Gaussian with random phases. The level of deviation from a Gaussian distribution and random phases, called non-Gaussianity, predicted by the simplest model of inflation is well below the current limit of measurement. Thus, any detection of non-Gaussianity would be a significant challenge to the currently favored models of the early universe.

The assumption of Gaussianity is motivated by the following view: the probability distribution of quantum fluctuations, P(φ), of free scalar fields in the ground state of the Bunch–Davies vacuum, φ, is a Gaussian distribution; thus, the probability distribution of primordial curvature perturbations (in the comoving gauge), ${\cal R}$, generated from φ (in the flat gauge) as ${\cal R}=-[H(\phi)/\dot{\phi _0}]\varphi$ (Mukhanov & Chibisov 1981; Hawking 1982; Starobinsky 1982; Guth & Pi 1982; Bardeen et al. 1983), would also be a Gaussian distribution. Here, H(ϕ) is the expansion rate during inflation and ϕ₀ is the mean field, that is, ϕ = ϕ₀ + φ.

This argument suggests that non-Gaussianity can be generated when (a) scalar fields are not free, but have some interactions, (b) there are nonlinear corrections to the relation between ${\cal R}$ and φ, and (c) the initial state is not in the Bunch–Davies vacuum.

For (a), one can think of expanding a general scalar field potential V(ϕ) to the cubic order or higher, $V(\phi)=\bar{V}+V^{\prime }\varphi +(1/2)V^{\prime \prime }\varphi ^2+(1/6)V^{\prime \prime \prime }\varphi ^3+\cdots$. The cubic (or higher-order) interaction terms can yield non-Gaussianity in φ (Falk et al. 1993). When perturbations in gravitational fields are included, there are many more interaction terms that arise from expanding the Ricci scalar to the cubic order, with coefficients containing derivatives of V and ϕ₀, such as $\dot{\phi _0}V^{\prime \prime }, \dot{\phi _0}^3/H$, etc. (Maldacena 2003).

For (b), one can think of this relation, ${\cal R}=-[H(\phi)/\dot{\phi _0}]\varphi$, as the leading-order term of a Taylor series expansion of the underlying nonlinear (gauge) transformation law between ${\cal R}$ and φ. Salopek & Bond (1990) showed that, in the single-field models, ${\cal R}=4\pi G\int _{\phi _0}^{\phi _0+\varphi }d\phi \left(\partial \ln H/\partial \phi \right)^{-1}$. Therefore, even if φ is precisely Gaussian, ${\cal R}$ can be non-Gaussian due to nonlinear terms such as φ² in a Taylor series expansion of this relation. One can write this relation in the following form, up to second order in ${\cal R}$,

$$\begin{equation} {\cal R} = {\cal R}_{\rm L} - \frac{1}{8\pi G}\left(\frac{\partial ^2\ln H}{\partial \phi ^2}\right){\cal R}_{\rm L}^2, \end{equation} \tag{ 32 }$$

where ${\cal R}_{\rm L}$ is a linear part of the curvature perturbation. We thus find that the second term makes ${\cal R}$ non-Gaussian, even when ${\cal R}_{\rm L}$ is precisely Gaussian. This formula has also been found independently by other researchers, extended to multifield cases, and often referred to as the "δN formalism" (Sasaki & Stewart 1996; Lyth et al. 2005; Lyth & Rodriguez 2005).

The observers like us, however, do not measure the primordial curvature perturbations, ${\cal R}$, directly. A more observationally-relevant quantity is the curvature perturbation during the matter era, Φ. At the linear order, these quantities are related by $\Phi =(3/5){\cal R}_{\rm L}$ (e.g., Kodama & Sasaki 1984), but the actual relation is more complicated at the nonlinear order (see Bartolo et al. 2004, for a review). In any case, this argument has motivated our defining the "local nonlinear coupling parameter," f^local_NL, as (Komatsu & Spergel 2001)³²

$$\begin{equation} \Phi = \Phi _{\rm L} + f_{\rm NL}^{\rm local}\Phi ^2_{\rm L}. \end{equation} \tag{ 33 }$$

If we take equation (32), for example, we find f^local_NL = −(5/24πG)(∂²ln H/∂ϕ²) (Komatsu 2001). Here, we have followed the terminology proposed by Babich et al. (2004) and called f^local_NL the "local" parameter, as both sides of Equation (33) are evaluated at the same location in space. (Hence, the term "local.")

Let us comment on the magnitude of the second term in Equation (33). Since Φ ∼ 10⁻⁵, the second term is smaller than the first term by 10⁻⁵f^local_NL; thus, the second term is only 0.1% of the first term for f^local_NL ∼ 10². As we shall see below, the existing limits on f^local_NL have already reached this level of Gaussianity, and thus it is clear that we are already talking about a tiny deviation from Gaussian fluctuations. This limit is actually better than the current limit on the spatial curvature, which is only on the order of 1%. Therefore, Gaussianity tests offer a stringent test of early universe models.

In the context of single-field inflation in which the scalar field slowly rolls down the potential, the quantities, H, V, and ϕ, change slowly. Therefore, one generically expects that f^local_NL is small, on the order of the so-called slow-roll parameters and η, which are typically of order 10⁻² or smaller. In this sense, the single-field, slow-roll inflation models are expected to result in a tiny amount of non-Gaussianity (Salopek & Bond 1990; Falk et al. 1993; Gangui et al. 1994; Maldacena 2003; Acquaviva et al. 2003). These contributions from the epoch of inflation are much smaller than those from the ubiquitous, second-order cosmological perturbations, that is, the nonlinear corrections to the relation between Φ and ${\cal R}$, which result in f^local_NL of order unity (Liguori et al. 2006; Smith & Zaldarriaga 2006). Also see Bartolo et al. (2004) for a review on this subject.

One can use the cosmological observations, such as the CMB data, to constrain f^local_NL (Verde et al. 2000; Komatsu & Spergel 2001). While the temperature anisotropy, ΔT/T, is related to Φ via the Sachs–Wolfe formula as ΔT/T = −Φ/3 at the linear order on very large angular scales (Sachs & Wolfe 1967), there are nonlinear corrections (nonlinear Sachs–Wolfe effect, nonlinear Integrated Sachs–Wolfe effect, gravitational lensing, etc.) to this relation, which add terms of order unity to f^local_NL by the time we observe it in CMB (Pyne & Carroll 1996; Mollerach & Matarrese 1997). On smaller angular scales, one must include the effects of acoustic oscillations of photon-baryon plasma by solving the Boltzmann equations. The second-order corrections to the Boltzmann equations can also yield f^local_NL of order unity (Bartolo et al. 2006, 2007).

Any detection of f^local_NL at the level that is currently accessible would have a profound implication for the physics of inflation. How can a large f^local_NL be generated? We essentially need to break either (a) single field or (b) slow-roll. For example, a multifield model, known as the curvaton scenario, can result in much larger values of f^local_NL (Linde & Mukhanov 1997; Lyth et al. 2003), and so can the models with field-dependent (variable) decay widths for reheating of the universe after inflation (Dvali et al. 2004b, 2004a). A more violent, nonlinear reheating process called "preheating" can give rise to a large f^local_NL (Enqvist et al. 2005; Jokinen & Mazumdar 2006; Chambers & Rajantie 2008).

Although breaking of slow-roll usually results in a premature termination of inflation, it is possible to break it temporarily for a brief period, without terminating inflation, by some features (steps, dips, etc.) in the shape of the potential. In such a scenario, a large non-Gaussianity may be generated at a certain limited scale at which the feature exists (Kofman et al. 1991; Wang & Kamionkowski 2000; Komatsu et al. 2003). The structure of non-Gaussianity from features is much more complex and model-dependent than f^local_NL (Chen et al. 2007a, 2008).

There is also a possibility that non-Gaussianity can be used to test alternatives to inflation. In a collapsing universe followed by a bounce (e.g., new Ekpyrotic scenario), f^local_NL is given by the inverse (as well as inverse-squared) of slow-roll parameters; thus, f^local_NL as large as of order 10–10² is a fairly generic prediction of this class of models (Creminelli & Senatore 2007; Koyama et al. 2007; Buchbinder et al. 2008; Lehners & Steinhardt 2008a, 2008b).

Using the angular bispectrum,³³ the harmonic transform of the angular three-point correlation function, Komatsu et al. (2002) have obtained the first observational limit on f^local_NL from the COBE 4-year data (Bennett et al. 1996), finding −3500 < f^local_NL < 2000 (95% CL). The uncertainty was large due to a relatively large beam size of COBE, which allowed us to go only to the maximum multipole of l_max = 20. Since the signal-to-noise ratio of f^local_NL is proportional to l_max, it was expected that the WMAP data would yield a factor of ∼50 improvement over the COBE data (Komatsu & Spergel 2001).

The full bispectrum analysis was not feasible with the WMAP data, as the computational cost scales as N^5/2_pix, where N_pix is the number of pixels, which is on the order of millions for the WMAP data. The "KSW" estimator (Komatsu et al. 2005) has solved this problem by inventing a cubic statistic that combines the triangle configurations of the bispectrum optimally so that it is maximally sensitive to f^local_NL.³⁴ The computational cost of the KSW estimator scales as N^3/2_pix. We give a detailed description of the method that we use in this paper in Appendix A.

We have applied this technique to the WMAP 1-year and 3-year data, and found −58 < f^local_NL < 134 (l = 265; Komatsu et al. 2003) and −54 < f^local_NL < 114 (l_max = 350; Spergel et al. 2007), respectively, at 95% CL. Creminelli et al. performed an independent analysis of the WMAP data and found similar limits: −27 < f^local_NL < 121 (l_max = 335; Creminelli et al. 2006) and −36 < f^local_NL < 100 (l = 370; Creminelli et al. 2007) for the 1-year and 3-year data, respectively. These constraints are slightly better than the WMAP team's, as their estimator for f^local_NL was improved from the original KSW estimator.

While these constraints are obtained from the KSW-like fast bispectrum statistics, many groups have used the WMAP data to measure f^local_NL using various other statistics, such as the Minkowski functionals (Komatsu et al. 2003; Spergel et al. 2007; Gott et al. 2007; Hikage et al. 2008), real-space three-point function (Gaztañaga & Wagg 2003; Chen & Szapudi 2005), integrated bispectrum (Cabella et al. 2006), 2-1 cumulant correlator power spectrum (Chen & Szapudi 2006), local curvature (Cabella et al. 2004), and spherical Mexican hat wavelet (Mukherjee & Wang 2004). The suborbital CMB experiments have also yielded constraints on f^local_NL: MAXIMA (Santos et al. 2003), VSA (Smith et al. 2004), Archeops (Curto et al. 2007), and BOOMERanG (De Troia et al. 2007).

We stress that it is important to use different statistical tools to measure f^local_NL if any signal is found, as different tools are sensitive to different systematics. The analytical predictions for the Minkowski functionals (Hikage et al. 2006) and the angular trispectrum (the harmonic transform of the angular 4-point correlation function; Okamoto & Hu 2002; Kogo & Komatsu 2006) as a function of f^local_NL are now available. Studies on the forms of the trispectrum from inflation models have just begun, and some important insights have been obtained (Boubekeur & Lyth 2006; Huang & Shiu 2006; Byrnes et al. 2006; Seery & Lidsey 2007; Seery et al. 2007; Arroja & Koyama 2008). It is now understood that the trispectrum is at least as important as the bispectrum in discriminating inflation models: some models do not produce any bispectra but produce significant trispectra, and other models produce similar amplitudes of the bispectra but produce very different trispectra (Huang & Shiu 2006; Buchbinder et al. 2008).

In this paper, we shall use the estimator that further improves upon Creminelli et al. (2006) by correcting an inadvertent numerical error of a factor of 2 in their derivation (Yadav et al. 2008). Yadav & Wandelt (2008) used this estimator to measure f^local_NL from the WMAP 3-year data. We shall also use the Minkowski functionals to find a limit on f^local_NL.

In addition to f^local_NL, we shall also estimate the "equilateral nonlinear coupling parameter," f^equil_NL, which characterizes the amplitude of the three-point function (i.e., the bispectrum) of the equilateral configurations, in which the lengths of all the three wave vectors forming a triangle in Fourier space are equal. This parameter is useful and highly complementary to the local one: while f^local_NL mainly characterizes the amplitude of the bispectrum of the squeezed configurations, in which two wave vectors are large and nearly equal and the other wave vector is small, and thus it is fairly insensitive to the equilateral configurations, f^equil_NL is mainly sensitive to the equilateral configurations with little sensitivity to the squeezed configurations. In other words, it is possible that one may detect f^local_NL without any detection of f^equil_NL and vice versa.

These two parameters cover a fairly large class of models. For example, f^equil_NL can be generated from inflation models in which the scalar field takes on the nonstandard (noncanonical) kinetic form, such as ${\cal L}=P(X,\phi)$, where X = (∂ϕ)². In this class of models, the effective sound speed of ϕ can be smaller than the speed of light, c²_s = [1 + 2X(∂²P/∂X²)/(∂P/∂X)]⁻¹ < 1. While the sign of f^equil_NL is negative for the DBI inflation, f^equil_NL ∼ −1/c²_s < 0 in the limit of c_s ≪ 1, it can be positive or negative for more general models (Seery & Lidsey 2005; Chen et al. 2007b; Cheung et al. 2008; Li et al. 2008). Such models can be realized in the context of String Theory via the noncanonical kinetic action called the DBI form (Alishahiha et al. 2004), and in the context of an IR modification of gravity called the ghost condensation (Arkani-Hamed et al. 2004).

The observational limits on f^equil_NL have been obtained from the WMAP 1-year and 3-year data as −366 < f^equil_NL < 238 (l = 405; Creminelli et al. 2006) and −256 < f^equil_NL < 332 (l_max = 475; Creminelli et al. 2007), respectively.

There are other forms, too. Warm inflation might produce a different form of f_NL(Moss & Xiong 2007; Moss & Graham 2007). Also, the presence of particles at the beginning of inflation, that is, a departure of the initial state of quantum fluctuations from the Bunch–Davies vacuum, can result in an enhanced non-Gaussianity in the "flattened" triangle configurations (Chen et al. 2007b; Holman & Tolley 2008). We do not consider these forms of non-Gaussianity in this paper.

In this paper, we do not discuss the non-Gaussian signatures that cannot be characterized by f^local_NL, f^equil_NL, or b_src (the point-source bispectrum amplitude). There have been many studies on non-Gaussian signatures in the WMAP data in various forms (Chiang et al. 2003, 2007; Naselsky et al. 2007; Park 2004; de Oliveira-Costa et al. 2004; Tegmark et al. 2003; Larson & Wandelt 2004; Eriksen et al. 2004a, 2004b, 2004c, 2007a; Copi et al. 2004, 2006, 2007; Schwarz et al. 2004; Gordon et al. 2005; Bielewicz et al. 2005; Jaffe et al. 2005, 2006; Vielva et al. 2004; Cruz et al. 2005, 2006, 2007b, 2007a; Cayón et al. 2005; Bridges et al. 2008; Wiaux et al. 2008; Räth et al. 2007; Land & Magueijo 2005b, 2005a, 2007; Rakić & Schwarz 2007; Park et al. 2007; Bernui et al. 2007; Hajian & Souradeep 2003, 2006; Hajian et al. 2005; Prunet et al. 2005; Hansen et al. 2004a, 2004b), many of which are related to the large-scale features at l ≲ 20. We expect these features to be present in the WMAP 5-year temperature map, as the structure of CMB anisotropy in the WMAP data on such large angular scales has not changed very much since the 3-year data.

The largest concern in measuring primordial non-Gaussianity from the CMB data is the potential contamination from the Galactic diffuse foreground emission. To test how much the results would be affected by this, we measure f_NL parameters from the raw temperature maps and from the foreground-reduced maps.

We shall mainly use the KQ75 mask, the new mask that is recommended for tests of Gaussianity (Gold et al. 2009). The important difference between the new mask and the previous Kp0 mask (Bennett et al. 2003c) is that the new mask is defined by the difference between the K-band map and the ILC map, and that between the Q band and ILC. Therefore, the CMB signal was absent when the mask was defined, which removes any concerns regarding a potential bias in the distribution of CMB on the masked sky.³⁵

To carry out tests of Gaussianity, one should use the KQ75 mask, which is slightly more conservative than Kp0, as the KQ75 mask cuts slightly more sky: we retain 71.8% of the sky with KQ75, while 76.5% with Kp0. To see how sensitive we are to the details of the mask, we also tried Kp0 as well as the new mask that is recommended for the power spectrum analysis, KQ85, which retains 81.7% of the sky. The previous mask that corresponds to KQ85 is the Kp2 mask, which retains 84.6% of the sky.

In addition, we use the KQ75p1 mask, which replaces the point-source mask of KQ75 with the one that does not mask the sources identified in the WMAP K-band data. Our point-source selection at K band removes more sources and sky in regions with higher CMB flux. We estimate the amplitude of this bias by using the KQ75p1 mask, which does not use any WMAP data for the point-source identification. The small change in f^local_NL shows that this is a small bias.

The unresolved extra-galactic point sources also contribute to the bispectrum (Refregier et al. 2000; Komatsu & Spergel 2001; Argüeso et al. 2003; Serra & Cooray 2008), and they can bias our estimates of primordial non-Gaussianity parameters such as f^local_NL and f^equil_NL. We estimate the bias by measuring f^local_NL and f^equil_NL from Monte Carlo simulations of point sources, and list them as Δf^local_NL and Δf^equil_NL in Tables 5 and 7, respectively. As the errors in these estimates of the bias are limited by the number of Monte Carlo realizations (which is 300), one may obtain a better estimate of the bias using more realizations.

We give a detailed description of our estimators for f^local_NL, f^equil_NL, and b_src, the amplitude of the point-source bispectrum, as well as of Monte Carlo simulations in Appendix A.

In Table 5, we show our measurement of f^local_NL from the template-cleaned V+W map (Gold et al. 2009) with four different masks, KQ85, Kp0, KQ75p1, and KQ75, in the increasing order of the size of the mask. For KQ85 and KQ75, we show the results from different maximum multipoles used in the analysis, l_max = 400, 500, 600, and 700. The WMAP 5-year temperature data are limited by cosmic variance to l ∼ 500.

We find that both KQ85 and Kp0 for l_max = 500 show evidence for f^local_NL > 0 at more than 95% CL, 9 < f^local_NL < 113 (95% CL), before the point-source bias correction, and 6.5 < f^local_NL < 110.5 (95% CL) after the correction. For a higher l_max, l_max = 700, we still find evidence for f^local_NL > 0, 1.5 < f^local_NL < 125.5 (95% CL), after the correction.³⁶

This evidence is, however, reduced when we use larger masks, KQ75p1 and KQ75. For the latter, we find −5 < f^local_NL < 115 (95% CL) before the source bias correction, and −9 < f^local_NL < 111 (95% CL) after the correction, which we take as our best estimate. This estimate improves upon our previous estimate from the 3-year data, −54 < f^local_NL < 114 (95% CL; Spergel et al. 2007, for l = 350), by cutting much of the allowed region for f^local_NL < 0. To test whether the evidence for f^local_NL > 0 can also be seen with the KQ75 mask, we need more years of WMAP observations.

Let us study the effect of mask further. We find that the central value of f^local_NL (without the source correction) changes from 61 for KQ85 to 55 for KQ75 at l_max = 500. Is this change expected? To study this, we have computed f^local_NL from each of the Monte Carlo realizations using KQ85 and KQ75. We find the root mean square (rms) scatter of 〈(f^local_NL^KQ85 − f^local_NL^KQ75)²〉^1/2_MC = 13, 12, 15, and 15, for l_max = 400, 500, 600, and 700, respectively. Therefore, the change in f^local_NL measured from the WMAP data is consistent with a statistical fluctuation. For the other masks at l_max = 500, we find 〈(f^local_NL^Kp0 − f^local_NL^KQ75)²〉^1/2_MC = 9.7 and 〈(f^local_NL^KQ75p1 − f^local_NL^KQ75)²〉^1/2_MC = 4.0.

In Table 6, we summarize the results from various tests. As a null test, we have measured f^local_NL from the difference maps such as Q−W and V−W, which are sensitive to non-Gaussianity in noise and the residual foreground emission, respectively. Since the difference maps do not contain the CMB signal, which is a source of a large cosmic variance in the estimation of f^local_NL, the errors in the estimated f^local_NL are much smaller. Before the foreground cleaning ("Raw" in the second column), we see negative values of f^local_NL, which is consistent with the foreground emission having positively skewed temperature distribution and f^local_NL > 0 mainly generating negative skewness. We do not find any significant signal of f^local_NL at more than 99% CL for raw maps, or at more than 68% CL for cleaned maps, which indicates that the temperature maps are quite clean outside of the KQ75 mask.

From the results presented in Table 6, we find that the raw-map results yield more scatter in f^local_NL estimated from various data combinations than the clean-map results. From these studies, we conclude that the clean-map results are robust against the data combinations, as long as we use only the V- and W-band data.

In Table 7, we show the equilateral bispectrum, f^equil_NL, from the template-cleaned V+W map with the KQ75 mask. We find that the point-source bias is much more significant for f^equil_NL: we detect the bias in f^equil_NL at more than the 5σ level for l_max = 600 and 700. After correcting for the bias, we find −151 < f^equil_NL < 253 (95% CL; l_max = 700) as our best estimate. Our estimate improves upon the previous one, −256 < f^equil_NL < 332 (95% CL; Creminelli et al. 2006, for l = 475), by reducing the allowed region from both above and below by a similar amount.

Finally, the bispectrum from very high multipoles, for example, l_max = 900, can be used to estimate the amplitude of residual point-source contamination. One can use this information to check for a consistency between the estimate of the residual point sources from the power spectrum and that from the bispectrum. In Table 8, we list our estimates of b_src. The raw maps and cleaned maps yield somewhat different values, indicating a possible leakage from the diffuse foreground to an estimate of b_src. Our best estimate in the Q band is b_src = 4.3 ± 1.3 μK³ sr² (68% CL). See Nolta et al. (2009) for the comparison between b_src, C_ps, and the point-source counts.

Incidentally, we also list b_src from the KQ75p1 mask, whose source mask is exactly the same as we used for the first-year analysis. We find b_src = 8.7 ± 1.3 μK³ sr² in the Q band, which is in an excellent agreement with the first-year result, b_src = 9.5 ± 4.4 μK³ sr² (Komatsu et al. 2003).

For the analysis of the Minkowski functionals, we follow the method described in Komatsu et al. (2003) and Spergel et al. (2007). In Figure 7, we show all of the three Minkowski functionals (Gott et al. 1990; Mecke et al. 1994; Schmalzing & Buchert 1997; Schmalzing & Gorski 1998; Winitzki & Kosowsky 1998) that one can define on a two-dimensional sphere: the cumulative surface area (bottom), the contour length (middle), and the Euler characteristics (which is also known as the genus; top), as a function the "threshold," ν, which is the number of σ of hot and cold spots, defined by

$$\begin{equation} \nu \equiv \frac{\Delta T}{\sigma _0}, \end{equation} \tag{ 34 }$$

where σ₀ is the standard deviation of the temperature data (which includes both signal and noise) at a given resolution of the map that one works with. We compare the Minkowski functionals measured from the WMAP data with the mean and dispersion of Gaussian realizations that include CMB signal and noise. We use the KQ75 mask and the V+W-band map.

Zoom In Zoom Out Reset image size

While Figure 7 shows the results at resolution 7 (N_side = 128), we have carried out Gaussianity tests using the Minkowski functionals at six different resolutions from resolution 3 (N_side = 8) to resolution 8 (N_side = 256). We find no evidence for departures from Gaussianity at any resolutions, as summarized in Table 9; in this table, we list the values of χ² of the Minkowski functionals relative to the Gaussian predictions:

$$\begin{eqnarray} \nonumber \chi ^2_{\it WMAP} &=& \sum _{ij}\sum _{\nu _1\nu _2} \left[F^i_{\it WMAP}-\langle F_{\rm sim}^i\rangle \right]_{\nu _1}\\ &&\times (\Sigma ^{-1})_{\nu _1\nu _2}^{ij} \big[F^j_{\it WMAP}-\langle F_{\rm sim}^j\rangle \big]_{\nu _2}, \end{eqnarray} \tag{ 35 }$$

where Fⁱ_WMAP and Fⁱ_sim are the ith Minkowski functionals measured from the WMAP data and Gaussian simulations, respectively, the angular bracket denotes the average over realizations, and $\Sigma _{\nu _1\nu _2}^{ij}$ is the covariance matrix estimated from the simulations. We use 15 different thresholds from ν = −3.5 to ν = +3.5, as indicated by the symbols in Figure 7, and thus the number of dof in the fit is 15 × 3 = 45. We show the values of χ²_WMAP and the dof in the second column, and the probability of having χ² that is larger than the measured value, F(>χ²_WMAP), in the third column. The smallest probability is 0.1 (at N_side = 8), and thus we conclude that the Minkowski functionals measured from the WMAP 5-year data are fully consistent with Gaussianity.

What do these results imply for f^local_NL? We find that the absence of non-Gaussianity at N_side = 128 and 64 gives the 68% limits on f^local_NL as f^local_NL = −57 ± 60 and −68 ± 69, respectively. The 95% limit from N_side = 128 is −178 < f^local_NL < 64. The errors are larger than those from the bispectrum analysis given in Section 3.5.3 by a factor of 2, which is partly because we have not used the Minkowski functional at all six resolutions to constrain f^local_NL. For a combined analysis of the WMAP 3-year data, see Hikage et al. (2008).

It is intriguing that the Minkowski functionals prefer a negative value of f^local_NL, f^local_NL ∼ −60, whereas the bispectrum prefers a positive value, f^local_NL ∼ 60. In the limit that non-Gaussianity is weak, the Minkowski functionals are sensitive to three "skewness parameters": (1) 〈(ΔT)³〉, (2) 〈(ΔT)²[∂²(ΔT)]〉, and (3) 〈[∂(ΔT)]²[∂²(ΔT)]〉, all of which can be written in terms of the weighted sum of the bispectrum; thus, the Minkowski functionals are sensitive to some selected configurations of the bispectrum (Hikage et al. 2006). It would be important to study where an apparent "tension" between the Minkowski functionals and the KSW estimator comes from. This example shows how important it is to use different statistical tools to identify the origin of non-Gaussian signals on the sky.

"Adiabaticity" of primordial fluctuations offers important tests of inflation as well as clues to the origin of matter in the universe. The negative correlation between the temperature and E-mode polarization (TE) at l ∼ 100 is a generic signature of adiabatic superhorizon fluctuations (Spergel & Zaldarriaga 1997; Peiris et al. 2003). The improved measurement of the TE power spectrum and the temperature power spectrum from the WMAP 5-year data, combined with the distance information from BAO and SN, now provide tight limits on deviations of primordial fluctuations from adiabaticity.

Adiabaticity may be loosely defined as the following relation between fluctuations in radiation density and those in matter density:

$$\begin{equation} \frac{3\delta \rho _r}{4\rho _r} = \frac{\delta \rho _m}{\rho _m}. \end{equation} \tag{ 36 }$$

This version³⁷ of the condition guarantees that the entropy density (dominated by radiation, s_r ∝ ρ^3/4_r) per matter particle is unperturbed, that is, δ(s_r/n_m) = 0.

There are two situations in which the adiabatic condition may be satisfied: (1) there is only one degree of freedom in the system, foe example, both radiation and matter were created from decay products of a single scalar field that was solely responsible for generating fluctuations, and (2) matter and radiation were in thermal equilibrium before any nonzero conserving quantum number (such as baryon number minus lepton number, B − L) was created (e.g., Weinberg 2004).

Therefore, detection of any nonadiabatic fluctuations, that is, any deviation from the adiabatic condition (Equation (36)), would imply that there were multiple scalar fields during inflation, and either matter (baryon or dark matter) was never in thermal equilibrium with radiation, or a nonzero conserving quantum number associated with matter was created well before the era of thermal equilibrium. In any case, the detection of nonadiabatic fluctuations between matter and radiation has a profound implication for the physics of inflation and, perhaps more importantly, the origin of matter.

For example, axions, a good candidate for dark matter, generate nonadiabatic fluctuations between dark matter and photons, as axion density fluctuations could be produced during inflation independent of curvature perturbations (which were generated from inflaton fields, and responsible for CMB anisotropies that we observe today), and were not in thermal equilibrium with radiation in the early universe (see Kolb & Turner 1990; Sikivie 2008, for reviews). We can, therefore, place stringent limits on the properties of axions by looking at a signature of deviation from the adiabatic relation in the CMB temperature and polarization anisotropies.

In this paper, we focus on the nonadiabatic perturbations between CDM and CMB photons. Nonadiabatic perturbations between baryons and photons are exactly the same as those between CDM and photons, up to an overall constant; thus, we shall not consider them separately in this paper. For neutrinos and photons, we consider only adiabatic perturbations. In other words, we consider only three standard neutrino species (i.e., no sterile neutrinos) and assume that the neutrinos were in thermal equilibrium before the lepton number was generated.

The basic idea behind this study is not new, and adiabaticity has been extensively constrained using the WMAP data since the first-year release, including general (phenomenological) studies without references to specific models (Peiris et al. 2003; Crotty et al. 2003a; Bucher et al. 2004; Moodley et al. 2004; Lazarides et al. 2004; Kurki-Suonio et al. 2005; Beltrán et al. 2005a; Dunkley et al. 2005; Bean et al. 2006; Trotta 2007; Keskitalo et al. 2007), as well as constraints on specific models such as double inflation (Silk & Turner 1987; Polarski & Starobinsky 1992, 1994), axion (Weinberg 1978; Wilczek 1978; Seckel & Turner 1985; Linde 1985, 1991; Turner & Wilczek 1991), and curvaton (Lyth & Wands 2003; Moroi & Takahashi 2001, 2002; Bartolo & Liddle 2002), all of which can be constrained from the limits on nonadiabatic fluctuations (Gordon & Lewis 2003; Gordon & Malik 2004; Beltrán et al. 2004, 2005b, 2007; Lazarides 2005; Parkinson et al. 2005; Kawasaki & Sekiguchi 2008).

We shall use the WMAP 5-year data, combined with the distance information from BAO and SN, to place more stringent limits on two types of nonadiabatic CDM fluctuations: (1) axion-type and (2) curvaton-type. Our study given below is similar to that by Kawasaki & Sekiguchi (2008) for the WMAP 3-year data.

We define the nonadiabatic, or entropic, perturbation between the CDM and photons, ${\cal S}_{c,\gamma }$, as

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

and report on the limits on the ratio of the power spectrum of ${\cal S}_{c,\gamma }, P_{\cal S}(k)$, to the curvature perturbation, $P_{\cal R}(k)$, at a given pivot wavenumber, k₀, given by (e.g., Bean et al. 2006)

$$\begin{equation} \frac{\alpha (k_0)}{1-\alpha (k_0)}\equiv \frac{P_{\cal S}(k_0)}{P_{\cal R}(k_0)}. \end{equation} \tag{ 38 }$$

We shall take k₀ to be 0.002 Mpc⁻¹.

While α parametrizes the ratio of the entropy power spectrum to the curvature power spectrum, it may be more informative to quantify "how much the adiabatic relation (Equation (36)) can be violated." To quantify this, we introduce the adiabaticity deviation parameter, δ_adi, given by

$$\begin{equation} \delta ^{(c,\gamma)}_{\rm adi}\equiv \frac{\delta \rho _c/\rho _c-3\delta \rho _\gamma /(4\rho _\gamma)}{\frac{1}{2}[\delta \rho _c/\rho _c+3\delta \rho _\gamma /(4\rho _\gamma)]}, \end{equation} \tag{ 39 }$$

which can be used to say, "the deviation from the adiabatic relation between dark matter and photons must be less than 100δ^(c,γ)_adi%." The numerator is just the definition of the entropy perturbation, ${\cal S}_{c,\gamma }$, whereas the denominator is given by

$$\begin{equation} \frac{1}{2}\left(\frac{\delta \rho _c}{\rho _c}+\frac{3\delta \rho _\gamma }{4\rho _\gamma }\right) = 3{\cal R} + O({\cal S}). \end{equation} \tag{ 40 }$$

Therefore, we find, up to the first order in ${\cal S}/{\cal R}$,

$$\begin{equation} \delta ^{(c,\gamma)}_{\rm adi}\approx \frac{\cal S}{3\cal R}\approx \frac{\sqrt{\alpha }}{3}, \end{equation} \tag{ 41 }$$

for α ≪ 1.

There could be a significant correlation between ${\cal S}_{c,\gamma }$ and ${\cal R}$ (Langlois 1999; Langlois & Riazuelo 2000; Gordon et al. 2001). We take this into account by introducing the cross-correlation coefficient, as

$$\begin{equation} -\beta (k_0)\equiv \frac{P_{S,\cal R}(k_0)}{\sqrt{P_{\cal S}(k_0)P_{\cal R}(k_0)}}, \end{equation} \tag{ 42 }$$

where $P_{{\cal S},\cal R}(k)$ is the cross-correlation power spectrum.

Here, we have a negative sign on the left-hand side because of the following reason. In our notation that we used for Gaussianity analysis, the sign convention of the curvature perturbation is such that it gives temperature anisotropy on large scales (the Sachs–Wolfe limit) as $\Delta T/T=-(1/5){\cal R}$. However, those who investigate correlations between ${\cal S}$ and ${\cal R}$ usually use an opposite sign convention for the curvature perturbation, such that the temperature anisotropy is given by

$$\begin{equation} \frac{\Delta T}{T} = \frac{1}{5}{\tilde{\cal R}} - \frac{2}{5}{\cal S} \end{equation} \tag{ 43 }$$

on large angular scales (Langlois 1999), where $\tilde{\cal R}\equiv -{\cal R}$, and define the correlation coefficient by

$$\begin{equation} \beta (k_0)= \frac{P_{S,\tilde{\cal R}}(k_0)}{\sqrt{P_{\cal S}(k_0)P_{\tilde{\cal R}}(k_0)}}. \end{equation} \tag{ 44 }$$

Therefore, in order to use the same sign convention for β as most of the previous work, we shall use Equation (42), and call β = +1 "totally correlated" and β = −1 "totally anticorrelated," entropy perturbations.

It is also useful to understand how the correlation or anticorrelation affects the CMB power spectrum at low multipoles. By squaring Equation (43) and taking the average, we obtain

$$\begin{equation} \frac{\langle (\Delta T)^2\rangle }{T^2} = \frac{1}{25}\left(P_{\tilde{\cal R}}+4P_{\cal S}-4\beta \sqrt{P_{\tilde{\cal R}}P_{\cal S}(k)}\right). \end{equation} \tag{ 45 }$$

Therefore, the "correlation," β>0, reduces the temperature power spectrum on low multipoles, whereas the "anticorrelation," β < 0, increases the power. This point will become important when we interpret our results: namely, models with β < 0 will result in a positive correlation between α and n_s (Gordon & Lewis 2003). Note that this property is similar to that of the tensor mode: as the tensor mode adds a significant power only to l ≲ 50, the tensor-to-scalar ratio, r, is degenerate with n_s (see Figure 2).

Finally, we specify the power spectrum of S as a pure power law,

$$\begin{equation} P_{\cal S}(k)\propto k^{m-4},\qquad P_{{\cal S},{\cal R}}(k)\propto k^{(m+n_s)/2-4}, \end{equation} \tag{ 46 }$$

in analogy to the curvature power spectrum, $P_{\cal R}(k)\propto k^{n_s-4}$. Note that β does not depend on k for this choice of $P_{{\cal S},{\cal R}}(k)$. With this parametrization, it is straightforward to compute the angular power spectra of the temperature and polarization of CMB.

In this paper, we shall pay attention to two limiting cases that are physically motivated: totally uncorrelated (β = 0) entropy perturbations, axion-type (Seckel & Turner 1985; Linde 1985, 1991; Turner & Wilczek 1991), and totally anticorrelated (β = −1) entropy perturbations, curvaton-type (Linde & Mukhanov 1997; Lyth & Wands 2003; Moroi & Takahashi 2001, 2002; Bartolo & Liddle 2002). Then, we shall use α₀ to denote α for β = 0 and α₋₁ for β = −1.

First, let us consider the axion case in which ${\cal S}$ and ${\cal R}$ are totally uncorrelated, that is, β = 0. This case represents the entropy perturbation between photons and axions, with axions accounting for some fraction of dark matter in the universe. For simplicity, we take the axion perturbations to be scale invariant, that is, m = 1. In Appendix B, we show that this choice corresponds to taking one of the slow-roll parameters, , to be less than 10⁻² or adding a tiny amount of gravitational waves, r ≪ 0.1, which justifies our ignoring gravitational waves in the analysis.

The left panel of Figure 8 shows that we do not find any evidence for the axion entropy perturbations. The limits are α₀ < 0.16(95%CL) and α₀ < 0.072(95%CL) for the WMAP-only analysis and WMAP+BAO+SN, respectively. The latter limit is the most stringent to date, from which we find the adiabaticity deviation parameter of δ^c,γ_adi < 0.089 (Equation (39)); thus, we conclude that the axion dark matter and photons should obey the adiabatic relation (Equation (36)) to 8.9%, at the 95% CL.

Zoom In Zoom Out Reset image size

We find that n_s and α₀ are strongly degenerate (see the middle panel of Figure 8). It is easy to understand the direction of correlation. As the entropy perturbation with a scale-invariant spectrum adds power to the temperature anisotropy on large angular scales only, the curvature perturbation tries to compensate it by reducing power on large scales with a larger tilt, n_s. However, since a larger n_s produces too much power on small angular scales, the fitting tries to increase Ω_bh² to suppress the second peak and reduce Ω_ch² to suppress the third peak. Overall, Ω_mh² needs to be reduced to compensate an increase in n_s, as shown in the right panel of Figure 8.

Adding the distance information from the BAO and SN helps to break the correlation between Ω_mh² and n_s by constraining Ω_mh², independent of n_s. Therefore, with WMAP+BAO+SN we find an impressive, factor of 2.2 improvement in the constraint on α₀.

What does this imply for the axions? It has been shown that the limit on the axion entropy perturbation can be used to place a constraint on the energy scale of inflation which, in turn, leads to a stringent constraint on the tensor-to-scalar ratio, r (Kain 2006; Beltrán et al. 2007; Sikivie 2008; Kawasaki & Sekiguchi 2008).

In Appendix B, we study a particular axion cosmology called the "misalignment angle scenario," in which the Pecci–Quinn symmetry breaking occurred during inflation and was never restored after inflation. In other words, we assume that the Pecci–Quinn symmetry breaking scale set by the axion decay constant, f_a, which has been constrained to be greater than 10¹⁰ GeV from the SN 1987 A (Yao et al. 2006), is at least greater than the reheating temperature of the universe after inflation. This is a rather reasonable assumption, as the reheating temperature is usually taken to be as low as 10⁸ GeV in order to avoid overproduction of unwanted relics (Pagels & Primack 1982; Coughlan et al. 1983; Ellis et al. 1986). Such a low reheating temperature is natural also because a coupling between inflaton and matter had to be weak; otherwise, it would terminate inflation prematurely.

There is another constraint. The Hubble parameter during inflation needs to be smaller than f_a, that is, H_inf ≲ f_a; otherwise, the Pecci–Quinn symmetry would be restored by quantum fluctuations (Lyth & Stewart 1992).

In this scenario, axions acquired quantum fluctuations during inflation, in the same way that inflaton fields would acquire fluctuations. These fluctuations were then converted to mass density fluctuations when axions acquired mass at the QCD phase transition at ∼200 MeV. We observe a signature of the axion mass density fluctuations via CDM-photon entropy perturbations imprinted in the CMB temperature and polarization anisotropies.

We find that the tensor-to-scalar ratio, r, the axion density, Ω_a, the CDM density, Ω_c, the phase of the Pecci–Quinn field within our observable universe, θ_a, and α₀ are related as (for an alternative expression that has f_a left instead of θ_a, see Equation (B7))

$$\begin{equation} 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{equation} \tag{ 47 }$$

$$\begin{equation} <\frac{(0.99\times 10^{-13})\alpha _{0}}{\theta _a^{10/7}\gamma ^{12/7}} \left(\frac{\Omega _c}{\Omega _a}\right)^{2/7} \end{equation} \tag{ 48 }$$

where γ ⩽ 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, for example, decay of some (unspecified) heavy particles, between 200 MeV and the epoch of nucleosynthesis, 1 MeV. Here, we have used the limit on the CDM density parameter, Ω_ch², from the axion entropy perturbation model that we consider here, Ω_ch² = 0.1052^+0.0068_−0.0070, as well as the observational fact that α₀ ≪ 1.

With our limit, α₀ < 0.072(95%CL), we find a limit on r within this scenario as

$$\begin{equation} r<\frac{6.6\times 10^{-15}}{\theta _a^{10/7}\gamma ^{12/7}}\left(\frac{\Omega _c}{\Omega _a}\right)^{2/7}. \end{equation} \tag{ 49 }$$

Therefore, in order for the axion dark matter scenario that we have considered here to be compatible with Ω_c ∼ Ω_a and the limits on the nonadiabaticity and Ω_ch², the energy scale of inflation should be low, and hence the gravitational waves are predicted to be negligible, unless the axion density was diluted severely by a late-time entropy production, γ ∼ 0.8 × 10⁻⁷ (for θ_a ∼ 1), the axion phase (or the misalignment angle) within our observable universe was close to zero, θ_a ∼ 3 × 10⁻⁹ (for γ ∼ 1), or both γ and θ_a were close to zero with lesser degree. All of these possibilities would give r ∼ 0.01, a value that could be barely detectable in the foreseeable future. One can also reverse Equation (49) to obtain

$$\begin{equation} \frac{\Omega _a}{\Omega _c}<\frac{3.0\times 10^{-39}}{\theta _a^5\gamma ^6}\left(\frac{0.01}{r}\right)^{7/2}. \end{equation} \tag{ 50 }$$

Therefore, the axion density would be negligible for the detectable r, unless θ_a, or γ, or both are tuned to be small.

Whether such an extreme production of entropy is highly unlikely, or such a tiny angle is an undesirable fine-tuning, can be debated. In any case, it is clear that the cosmological observations, such as the CDM density, entropy perturbations, and gravitational waves, can be used to place a rather stringent limit on the axion cosmology based upon the misalignment scenario, one of the most popular scenarios for axions to become a dominant dark matter component in the universe.

Next, let us consider one of the curvaton models in which ${\cal S}$ and $\tilde{\cal R}$ are totally anticorrelated, that is, β = −1 (Lyth & Wands 2003; Moroi & Takahashi 2001, 2002; Bartolo & Liddle 2002). One can also write ${\cal S}$ as ${\cal S}=B{\cal R}=-B\tilde{\cal R}$ where B > 0; thus, B² = α₋₁/(1 − α₋₁).³⁸ We take the spectral index of the curvaton entropy perturbation, m, to be the same as that of the adiabatic perturbation, n_s, that is, n_s = m.

The left panel of Figure 9 shows that we do not find any evidence for the curvaton entropy perturbations, either. The limits, α₋₁ < 0.011(95%CL) and α₋₁ < 0.0041(95%CL) for the WMAP-only analysis and WMAP+BAO+SN, respectively, are more than a factor of 10 better than those for the axion perturbations. The WMAP-only limit is better than the previous limit by a factor of 4 (Bean et al. 2006). From the WMAP+BAO+SN limit, we find the adiabaticity deviation parameter of δ^(c,γ)_adi < 0.021 (Equation (39)); thus, we conclude that the curvaton dark matter and photons should obey the adiabatic relation (Equation (36)) to 2.1% at the 95% CL.

Zoom In Zoom Out Reset image size

Once again, adding the distance information from the BAO and SN helps to reduce the correlation between n_s and Ω_mh² (see the right panel of Figure 9) and reduces the correlation between n_s and α₋₁. The directions in which these parameters are degenerate are similar to those for the axion case (see Figure 8), as the entropy perturbation with β = −1 also increases the CMB temperature power spectrum on large angular scales, as we described in Section 3.6.2.

What is the implication for this type of curvaton scenario, in which β = −1? This scenario would arise when CDM was created from the decay products of the curvaton field. One then finds a prediction (Lyth et al. 2003)

$$\begin{equation} \frac{\alpha _{-1}}{1-\alpha _{-1}}\approx 9\left(\frac{1-\rho _{\rm curvaton}/{\rho _{\rm total}}}{\rho _{\rm curvaton}/\rho _{\rm total}}\right)^2, \end{equation} \tag{ 51 }$$

where ρ_curvaton and ρ_total are the curvaton density and total density at the curvaton decay, respectively. Note that there would be no entropy perturbation if curvaton dominated the energy density of the universe completely at the decay. The reason is simple: in such a case, all of the curvaton perturbation would become the adiabatic perturbation, so would the CDM perturbation. Our limit, α₋₁ < 0.0041(95%CL), indicates that ρ_curvaton/ρ_total is close to unity, which simplifies the relation (Equation (51)) to give

$$\begin{equation} \frac{\rho _{\rm curvaton}}{\rho _{\rm total}} \approx 1 - \frac{\sqrt{\alpha _{-1}}}{3} = 1-\delta _{\rm adi}^{(c,\gamma)}. \end{equation} \tag{ 52 }$$

Note that it is the adiabaticity deviation parameter given by Equation (39) that gives the deviation of ρ_curvaton/ρ_total from unity. From this result, we find

$$\begin{equation} 1 \ge \frac{\rho _{\rm curvaton}}{\rho _{\rm total}}\gtrsim 0.98\qquad \mbox{(95\% CL)}. \end{equation} \tag{ 53 }$$

As we mentioned in Section 3.5.1, the curvaton scenario is capable of producting the local form of non-Gaussianity, and f^local_NL is given by (Lyth & Rodriguez 2005, and references therein³⁹)

$$\begin{equation} f_{\rm NL}^{\rm local}= \frac{5\rho _{\rm total}}{4\rho _{\rm curvaton}}-\frac{5}{3} -\frac{5\rho _{\rm curvaton}}{6\rho _{\rm total}}, \end{equation} \tag{ 54 }$$

which gives −1.25 ⩽ f^local_NL(curvaton) ≲ −1.21, for α₋₁ < 0.0041(95%CL). While we need to add additional contributions from postinflationary, nonlinear gravitational perturbations of order unity to this value in order to compare with what we measure from CMB, the limit from the curvaton entropy perturbation is consistent with the limit from the measured f^local_NL (see Section 3.5.3).

However, should the future data reveal f^local_NL ≫ 1, then either this scenario would be ruled out (Beltrán 2008; Li et al. 2008b) or the curvaton dark matter must have been in thermal equilibrium with photons.

For the other possibilities, including possible baryon entropy perturbations, see Gordon & Lewis (2003).

Dark energy influences the distance scales as well as the growth of structure. The CMB power spectra are sensitive to both, although sensitivity to the growth of structure is fairly limited, as it influences the CMB power spectrum via the ISW effect at low multipoles (l ≲ 10), whose precise measurement is hampered by a large cosmic variance.

However, CMB is sensitive to the distance to the decoupling epoch via the locations of peaks and troughs of the acoustic oscillations, which can be measured precisely. More specifically, CMB measures two distance ratios: (1) the angular diameter distance to the decoupling epoch divided by the sound horizon size at the decoupling epoch, D_A(z_*)/r_s(z_*), and (2) the angular diameter distance to the decoupling epoch divided by the Hubble horizon size at the decoupling epoch, D_A(z_*)H(z_*)/c. This consideration motivates our using these two distance ratios to constrain various dark energy models, in the presence of the spatial curvature, on the basis of distance information (Wang & Mukherjee 2007; Wright 2007).

We shall quantify the first distance ratio, D_A(z_*)/r_s(z_*), by the "acoustic scale," l_A, defined by

$$\begin{equation} l_{\rm A}\equiv (1+z_*)\frac{\pi D_A(z_*)}{r_s(z_*)}, \end{equation} \tag{ 65 }$$

where a factor of (1 + z_*) arises because D_A(z_*) is the proper (physical) angular diameter distance (Equation (2)), whereas r_s(z_*) is the comoving sound horizon at z_* (Equation (6)). Here, we shall use the fitting function of z_* proposed by Hu & Sugiyama (1996):

$$\begin{equation} z_*=1048[1+0.00124(\Omega _{b}h^2)^{-0.738}] [1+g_1(\Omega _{m}h^2)^{g_2}], \end{equation} \tag{ 66 }$$

where

$$\begin{equation} g_1=\frac{0.0783(\Omega _bh^2)^{-0.238}}{1+39.5(\Omega _bh^2)^{0.763}}, \end{equation} \tag{ 67 }$$

$$\begin{equation} g_2=\frac{0.560}{1+21.1(\Omega _bh^2)^{1.81}}. \end{equation} \tag{ 68 }$$

Note that one could also use the peak of the probability of last scattering of photons, that is, the peak of the visibility function, to define the decoupling epoch, which we denote as z_dec. Both quantities yield similar values. We shall adopt z_* here because it is easier to compute, and, therefore, it allows one to implement the WMAP distance priors in a straightforward manner.

The second distance ratio, D_A(z_*)H(z_*)/c, is often called the "shift parameter," R, given by (Bond et al. 1997)

$$\begin{equation} R(z_*)\equiv \frac{\sqrt{\Omega _{m}H_0^2}}{c}(1+z_*)D_A(z_*). \end{equation} \tag{ 69 }$$

This quantity is different from D_A(z_*)H(z_*)/c by a factor of $\sqrt{1+z_*}$, and also ignores the contributions from radiation, curvature, or dark energy to H(z_*). Nevertheless, we shall use R to follow the convention in the literature.

We give the 5-year WMAP constraints on l_A, R, and z_* that we recommend as the WMAP distance priors for constraining dark energy models. However, we note an important caveat. As pointed out by Elgarøy & Multamäki (2007) and Corasaniti & Melchiorri (2008), the derivation of the WMAP distance priors requires us to assume the underlying cosmology first, as all of these quantities are derived parameters from fitting the CMB power spectra. Therefore, one must be careful about which model one is testing. Here, we give the WMAP distance priors, assuming the following model.

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

In Figure 13, we show the constraints on w and Ω_k from the WMAP distance priors (combined with BAO and SN). We find a good agreement with the full MCMC results (compare the middle and right panels of Figure 12 with the left and right panels of Figure 13, respectively). The constraints from the WMAP distance priors are slightly weaker than the full MCMC, as the distance priors use only a part of the information contained in the WMAP data.

Zoom In Zoom Out Reset image size

Of course, the agreement between Figures 12 and 13 does not guarantee that these priors yield good results for the other, more complex dark energy models with a time-dependent w; however, the previous studies indicate that, under the assumptions given above, these priors can be used to constrain a wide variety of dark energy models (Wang & Mukherjee 2007; Elgarøy & Multamäki 2007; Corasaniti & Melchiorri 2008). Also see Li et al. (2008a) for the latest comparison between the WMAP distance priors and the full analysis.

Here is the prescription for using the WMAP distance priors.

• 1.  
  For a given Ω_bh² and Ω_mh², compute z_* from Equation (66).
• 2.  
  For a given H₀, Ω_mh², Ω_rh² (which includes N_eff = 3.04), Ω_Λ, and w(z), compute the expansion rate, H(z), from Equation (7), as well as the comoving sound horizon size at z_*, r_s(z_*), from Equation (6).
• 3.  
  For a given Ω_k and H(z) from the previous step, compute the proper angular diameter distance, D_A(z), from Equation (2).
• 4.  
  Use Equations (65) and (69) to compute l_A(z_*) and R(z_*), respectively.
• 5.  
  Form a vector containing x_i = (l_A, R, z_*) in this order.
• 6.  
  Use Table 10 for the data vector, d_i = (l^WMAP_A, R^WMAP, z^WMAP_*). We recommend the maximum likelihood (ML) values.
• 7.  
  Use Table 11 for the inverse covariance matrix, (C⁻¹)_ij.
• 8.  
  Compute the likelihood, L, as χ²_WMAP ≡ −2ln L = (x_i − d_i)(C⁻¹)_ij(x_j − d_j).
• 9.  
  Add this to the favorite combination of the cosmological datasets. In this paper we add χ²_WMAP to the BAO and SN data, that is, χ²_total = χ²_WMAP + χ²_BAO + χ²_SN.
• 10.  
  Marginalize the posterior distribution over Ω_bh², Ω_mh², and H₀ with uniform priors. Since the WMAP distance priors combined with the BAO and SN data provide tight constraints on these parameters, the posterior distribution of these parameters is close to a Gaussian distribution. Therefore, the marginalization is equivalent to minimizing χ²_total with respect to Ω_bh², Ω_mh², and H₀ (also see Cash 1976; Wright 2007). We use a downhill simplex method for minimization (amoeba routine in Numerical Recipes; Press et al. 1992). The marginalization over Ω_bh², Ω_mh², and H₀ leaves us with the the marginalized posterior distribution of the dark energy function, w(a), and the curvature parameter, Ω_k.

Note that this prescription eliminates the need for running the MCMC entirely, and thus the computational cost for evaluating the posterior distribution of w(a) and Ω_k is not demanding. In Section 5.4.2, we shall apply the WMAP distance priors to constrain the dark energy equation of state that depends on redshifts, w = w(z).

For those who wish to include an additional prior on Ω_bh², we give the inverse covariance matrix for the "extended" WMAP distance priors: (l_A(z_*), R(z_*), z_*, 100Ω_bh²), as well as the maximum likelihood value of 100Ω_bh², in Table 12. We note, however, that it is sufficient to use the reduced set, (l_A(z_*), R(z_*), z_*), as the extended WMAP distance priors give very similar constraints on dark energy (see Wang 2008).

In this subsection, we explore a time-dependent equation of state of dark energy, w(z). We use the following parametrized form:

$$\begin{equation} w(a) = \frac{a\tilde{w}(a)}{a+a_{\rm trans}} - \frac{a_{\rm trans}}{a+a_{\rm trans}}, \end{equation} \tag{ 70 }$$

where $\tilde{w}(a)=\tilde{w}_0+(1-a)\tilde{w}_a$. We give a motivation, derivation, and detailed discussion on this form of w(a) in Appendix C. This form has a number of desirable properties, as follows.

• 1.  
  w(a) approaches −1 at early times, a < a_trans, where a_trans = 1/(1 + z_trans) is the "transition epoch" and z_trans is the transition redshift. Therefore, the dark energy density tends to a constant value at a < a_trans.
• 2.  
  The dark energy density remains totally subdominant relative to the matter density at the decoupling epoch.
• 3.  
  We recover the widely-used linear form, w(a) = w₀ + (1 − a)w_a (Chevallier & Polarski 2001; Linder 2003), at late times, a > a_trans.
• 4.  
  The early-time behavior is consistent with some of scalar field models classified as the "thawing models" (Caldwell & Linder 2005) in which a scalar field was moving very slowly at early times and then began to move faster at recent times.
• 5.  
  Since the late-time form of w(a) allows w(a) to go below −1, our form is more general than models based upon a single scalar field.
• 6.  
  The form is simple enough to give a closed, analytical form of the effective equation of state, w_eff(a) = (ln a)⁻¹∫^ln a₀dln a'w(a') (Equation C6), which determines the evolution of the dark energy density, $\rho _{\rm de}(a)=\rho _{\rm de}(0)a^{-3[1+w_{\rm eff}(a)]}$; thus, it allows one to easily compute the evolution of the expansion rate and cosmological distances.

While this form contains three free parameters, $\tilde{w}_0,\tilde{w}_a$, and z_trans, we shall give constraints on the present-day value of w, w₀ ≡ w(a = 1), and the first derivative of w at present, w' ≡.dw/dz|_{z = 0}, instead of $\tilde{w}_0$ and $\tilde{w}_a$, as they can be compared to the previous results in the literature more directly. We find that the results are not sensitive to the exact values of z_trans.

In Figure 14, we present the constraint on w₀ and w' that we have derived from the WMAP distance priors (l_A, R, z_*), combined with the BAO and SN data. Note that we have assumed a flat universe in this analysis, although it is straightforward to include the spatial curvature. Wang & Mukherjee (2007) and Wright (2007) showed that the two-dimensional distribution extends more towards southeast, that is, w > − 1 and w' < 0, when the spatial curvature is allowed.

Zoom In Zoom Out Reset image size

The 95% limit on w₀ for z_trans = 10 is −0.33 < 1 + w₀ < 0.21.⁴⁷ Our results are consistent with the previous work using the WMAP 3-year data (see Zhao et al. 2007; Wang & Mukherjee 2007; Wright 2007; Lazkoz et al. 2008, for recent work and references therein). The WMAP 5-year data help tighten the upper limit on w' and the lower limit on w₀, whereas the lower limit on w' and the upper limit on w₀ come mainly from the Type Ia SN data. As a result, the lower limit on w' and the upper limit on w₀ are sensitive to whether we include the systematic errors in the SN data. For this investigation, see Appendix D.

Alternatively, one may take the linear form, w(a) = w₀ + (1 − a)w_a, literally and extend it to an arbitrarily high redshift. This can result in an undesirable situation in which the dark energy is as important as the radiation density at the epoch of the BBN; however, one can severely constrain such a scenario by using the limit on the expansion rate from BBN (Steigman 2007). We follow Wright (2007) to adopt a Gaussian prior on

$$\begin{eqnarray} \nonumber & & \sqrt{1+\frac{\Omega _\Lambda (1+z_{\rm BBN})^{3[1+w_{\rm eff}(z_{\rm BBN})]}}{\Omega _m(1+z_{\rm BBN})^3+\Omega _r(1+z_{\rm BBN})^4+\Omega _k(1+z_{\rm BBN})^2}}\\ & &\ \ \ \quad= 0.942 \pm 0.030, \end{eqnarray} \tag{ 71 }$$

where we have kept Ω_m and Ω_k for definiteness, but they are entirely negligible compared to the radiation density at the redshift of BBN, z_BBN = 10⁹. Figure 15 shows the constraint on w₀ and w' for the linear evolution model derived from the WMAP distance priors, the BAO and SN data, and the BBN prior. The 95% limit on w₀ is −0.29 < 1 + w₀ < 0.21,⁴⁸ which is similar to what we have obtained above.

Zoom In Zoom Out Reset image size

The existence of nonzero neutrino masses has been firmly established by the experiments detecting atmospheric neutrinos (Hirata et al. 1992; Fukuda et al. 1994, 1998; Allison et al. 1999; Ambrosio et al. 2001), solar neutrinos (Davis et al. 1968; Cleveland et al. 1998; Hampel et al. 1999; Abdurashitov et al. 1999; Fukuda et al. 2001b, 2001a; Ahmad et al. 2002; Ahmed et al. 2004), reactor neutrinos (Eguchi et al. 2003; Araki et al. 2005), and accelerator beam neutrinos (Ahn et al. 2003; Michael et al. 2006). These experiments have placed stringent limits on the squared mass differences between the neutrino mass eigenstates, Δm²₂₁ ≃ 8 × 10⁻⁵ eV², from the solar and reactor experiments, and Δm²₃₂ ≃ 3 × 10⁻³ eV², from the atmospheric and accelerator beam experiments.

One needs different experiments to measure the absolute masses. The next-generation tritium β-decay experiment, KATRIN,⁴⁹ is expected to reach the electron neutrino mass of as small as ∼0.2 eV. Cosmology has also been providing useful limits on the mass of neutrinos (see Dolgov 2002; Elgarøy & Lahav 2005; Tegmark 2005; Lesgourgues & Pastor 2006; Fukugita 2006; Hannestad 2006b, for reviews). Since the determination of the neutrino mass is of fundamental importance in physics, there is enough motivation to pursue cosmological constraints on the neutrino mass.

How well can CMB constrain the mass of neutrinos? We do not expect massive neutrinos to affect the CMB power spectra very much (except through the gravitational lensing effect), if they were still relativistic at the decoupling epoch, z ≃ 1090. This means that, for massive neutrinos to affect the CMB power spectra, at least one of the neutrino masses must be greater than the mean energy of relativistic neutrinos per particle at z ≃ 1090 when the photon temperature of the universe is T_γ ≃ 3000 K ≃ 0.26 eV. Since the mean energy of relativistic neutrinos is given by 〈E〉 = (7π⁴T_ν)/(180ζ(3)) ≃ 3.15T_ν = 3.15(4/11)^1/3T_γ, we need at least one neutrino species whose mass satisfies m_ν > 3.15(4/11)^1/3T_γ ≃ 0.58 eV; thus, it would not be possible to constrain the neutrino mass using the CMB data alone, if the mass of the heaviest neutrino species is below this value.

If the neutrino mass eigenstates are degenerate with the effective number of species equal to 3.04, this argument suggests that ∑m_ν ∼ 1.8 eV would be the limit to which the CMB data are sensitive. Ichikawa et al. (2005) argued that ∑m_ν ∼ 1.5 eV would be the limit for the CMB data alone, which is fairly close to the value given above.

In order to go beyond ∼1.5 eV, therefore, one needs to combine the CMB data with the other cosmological probes. We shall combine the WMAP data with the distance information from BAO and SN to place a limit on the neutrino mass. We shall not use the galaxy power spectrum in this paper, and, therefore, our limit on the neutrino mass is free from the uncertainty in the galaxy bias. We discuss this in more detail in Section 6.1.2.

We assume that, for definiteness, the neutrino mass eigenstates are degenerate, which means that all of the three neutrino species have equal masses.⁵⁰ We measure the neutrino mass density parameter, Ω_νh², and convert it to the total mass, ∑m_ν, via

$$\begin{equation} \sum m_\nu \equiv N_\nu m_\nu = 94\,{\rm eV}(\Omega _\nu h^2), \end{equation} \tag{ 81 }$$

where N_ν is the number of massive neutrino species. We take it to be 3.04. Note that in this case, the mass density parameter is the sum of baryons, CDM, and neutrinos: Ω_m = Ω_b + Ω_c + Ω_ν.

Since the release of the 1 year (Spergel et al. 2003) and 3 year (Spergel et al. 2007) results on the cosmological analysis of the WMAP data, there have been a number of studies with regard to the limits on the mass of neutrinos (Hannestad 2003; Elgarøy & Lahav 2003; Allen et al. 2003; Tegmark et al. 2004a; Barger et al. 2004; Hannestad & Raffelt 2004; Crotty et al. 2004; Seljak et al. 2005a, 2005b; Ichikawa et al. 2005; Hannestad 2005; Lattanzi et al. 2005; Hannestad & Raffelt 2006; Goobar et al. 2006; Feng et al. 2006; Lesgourgues et al. 2007). These analyses reached different limits depending upon (1) the choice of datasets and (2) the parameters in the cosmological model.

The strongest limits quoted on neutrino masses come from combining CMB measurements with measurements of the amplitude of density fluctuations in the recent universe. Clustering of galaxies and Lyα forest observations have been used to obtain some of the strongest limits on neutrino masses (Seljak et al. 2005b, 2006; Viel et al. 2006). As the neutrino mass increases, the amplitude of mass fluctuations on small scales decreases (Bond et al. 1980; Bond & Szalay 1983; Ma 1996; also see the middle panel of Figure 17), which can be used to weigh neutrinos in the universe (Hu et al. 1998; Lesgourgues & Pastor 2006).

These analyses are sensitive to correctly calculating the relationship between the level of observed fluctuations in galaxies (or gas) and the mass fluctuation with the strongest limits coming from the smallest scales in the analyses. With the new WMAP data, these limits are potentially even stronger. There are, however, several potential concerns with this limit: there is already "tension" between the high level of fluctuations seen in the Lyman alpha forest and the amplitude of mass fluctuations inferred from WMAP (Lesgourgues et al. 2007), the relationship between gas temperature and density appears to be more complicated than assumed in the previous Lyα forest analysis (Kim et al. 2007; Bolton et al. 2008), and additional astrophysics could potentially be the source of some of the small-scale fluctuations seen in the Lyα forest data. Given the power of the Lyα forest data, it is important to address these issues; however, they are beyond the scope of this paper.

In this paper, we take the more conservative approach and use only the WMAP data and the distance measures to place limits on the neutrino masses. Our approach is more conservative than Goobar et al. (2006), who have found a limit of 0.62 eV on the sum of the neutrino mass from the WMAP 3-year data, the SDSS–LRG BAO measurement, the SNLS SN data, and the shape of the galaxy power spectra from the SDSS main sample and 2dFGRS. While we use the WMAP 5-year data, the SDSS+2dFGRS BAO measurements, and the union SN data, we do not use the shape of the galaxy power spectra. See Section 2.3 in this paper or Dunkley et al. (2009) for more detail on this choice.

In summary, we do not use the amplitude or shape of the matter power spectrum, but exclusively rely on the CMB data and the distance measurements. As a result, our limits are weaker than the strongest limits in the literature.

Next, let us discuss (2), the choice of parameters. A few correlations between the neutrino mass and other cosmological parameters have been identified. Hannestad (2005) has found that the limit on the neutrino mass degrades significantly when the dark energy equation of state, w, is allowed to vary (also see Figure 18 of Spergel et al. 2007). This correlation would arise only when the amplitude of the galaxy or Lyα forest power spectrum was included, as the dark energy equation of state influences the growth rate of the structure formation. Since we do not include them, our limit on the neutrino mass is not degenerate with w. We shall come back to this point later in the Section 6.1.3. Incidentally, our limit is not degenerate with the running index, dn_s/dln k, or the tensor-to-scalar ratio, r.

Figure 17 summarizes our limits on the sum of neutrino masses, ∑m_ν.

With the WMAP data alone we find ∑m_ν < 1.3 eV(95%CL) for the ΛCDM model in which w = −1, and ∑m_ν < 1.5 eV(95%CL) for w ≠ −1. We assume a flat universe in both cases. These constraints are very similar, which means that w and ∑m_ν are not degenerate. We show this more explicitly on the right panel of Figure 17.

When the BAO and SN data are added, our limits improve significantly, by a factor 2, to ∑m_ν < 0.67 eV(95%CL) for w = −1, and ∑m_ν < 0.80 eV(95%CL) for w ≠ −1. Again, we do not observe much correlation between w and ∑m_ν. While the distances out to either BAO or SN cannot reduce correlation between Ω_m (or H₀) and w, a combination of the two can reduce this correlation effectively, leaving little correlation left on the right panel of Figure 17.

What information do BAO and SN add to improve the limit on ∑m_ν? It's the Hubble constant, H₀, as shown in the left panel of Figure 17. This effect has been explained by Ichikawa et al. (2005) as follows.

The massive neutrinos modify the CMB power spectrum by their changing the matter-to-radiation ratio at the decoupling epoch. If the sum of degenerate neutrino masses is below 1.8 eV, the neutrinos were still relativistic at the decoupling epoch. However, they are definitely nonrelativistic at the present epoch, as the neutrino oscillation experiments have shown that at least one neutrino species is heavier than 0.05 eV. This means that the Ω_m that we measure must be the sum of Ω_b, Ω_c, and Ω_ν; however, at the decoupling epoch, neutrinos were still relativistic, and thus the matter density at the decoupling epoch was actually smaller than a naive extrapolation from the present value.

As the matter-to-radiation ratio was smaller than one would naively expect, it would accelerate the decay of gravitational potential around the decoupling epoch. This leads to an enhancement in the so-called early integrated Sachs–Wolfe (ISW) effect. The larger ∑m_ν is, the larger early ISW becomes, as long as the neutrinos were still relativistic at the decoupling epoch, that is, ∑m_ν ≲ 1.8 eV.

The large ISW causes the first peak position to shift to lower multipoles by adding power at l ∼ 200; however, this shift can be absorbed by a reduction in the value of H₀.⁵¹ This is why ∑m_ν and H₀ are anticorrelated (see Ichikawa et al. 2005, for a further discussion on this subject).

It is the BAO distance that provides a better limit on H₀, as BAO is an absolute distance indicator. The SN is totally insensitive to H₀, as their absolute magnitudes have been marginalized over (SN is a relative distance indicator); however, the SN data do help reduce the correlation between w and H₀ when w is allowed to vary. As a result, we have equally tight limits on ∑m_ν regardless of w.

Our limit, ∑m_ν < 0.67 eV(95%CL) (for w = −1), is weaker than the best limit quoted in the literature, as we have not used the information on the amplitude of fluctuations traced by the large-scale structure. The middle panel of Figure 17 shows how the WMAP data combined with BAO and SN predict the present-day amplitude of matter fluctuations, σ₈, as a function of ∑m_ν. From this, it is clear that an accurate, independent measurement of σ₈ will reduce the correlation between σ₈ and ∑m_ν, and provide a significant improvement in the limit on ∑m_ν.

Improving upon our understanding of nonlinear astrophysical effects, such as those raised by Bolton et al. (2008) for the Lyα forest data and Sánchez & Cole (2008) for the SDSS and 2dFGRS data, is a promising way to improve upon the numerical value of the limit, as well as the robustness of the limit, on the mass of neutrinos.

While the absolute mass of neutrinos is unknown, the number of neutrino species is well known: it is 3. The high precision measurement of the decay width of Z into neutrinos (the total decay width minus the decay width to quarks and charged leptons), carried out by LEP using the production of Z in e⁺e⁻ collisions, has yielded N_ν = 2.984 ± 0.008 (Yao et al. 2006). However, are there any other particles that we do not know yet, and that are relativistic at the photon decoupling epoch?

Such extra relativistic particle species can change the expansion rate of the universe during the radiation era. As a result, they change the predictions from the BBN for the abundance of light elements such as helium and deuterium (Steigman et al. 1977). One can use this property to place a tight bound on the relativistic dof, expressed in terms of the "effective number of neutrino species," N_eff (see Equation (84) for the precise definition). As the BBN occurred at the energy of ∼0.1 MeV, which is later than the neutrino decoupling epoch immediately followed by e⁺e⁻ annihilation, the value of N_eff for three neutrino species is slightly larger than 3. With other subtle corrections included, the current standard value is N^standard_eff = 3.04 (Dicus et al. 1982; Gnedin & Gnedin 1998; Dolgov et al. 1999; Mangano et al. 2002). The 2σ interval for N_eff from the observed helium abundance, Y_P = 0.240 ± 0.006, is 1.61 < N_eff < 3.30 (see Steigman 2007, for a recent summary).

Many people have been trying to find evidence for the extra relativistic dof in the universe, using the cosmological probes such as CMB and the large-scale structure (Pierpaoli 2003; Hannestad 2003; Crotty et al. 2003b, 2004; Barger et al. 2003; Trotta & Melchiorri 2005; Lattanzi et al. 2005; Cirelli & Strumia 2006; Hannestad 2006a; Ichikawa et al. 2007; Mangano et al. 2007; Hamann et al. 2007; de Bernardis et al. 2008). There is a strong motivation to seek the answer for the following question, "can we detect the cosmic neutrino background, and confirm that the signal is consistent with the expected number of neutrino species that we know?" Although we cannot detect the cosmic neutrino background directly yet, there is a possibility that we can detect it indirectly by looking for the signatures of neutrinos in the CMB power spectrum.

In this section, we shall revisit this classical problem by using the WMAP 5-year data as well as the distance information from BAO and SN, and Hubble's constant measured by HST.

It is common to write the energy density of neutrinos (including antineutrinos), when they were still relativistic, in terms of the effective number of neutrino species, N_eff, as

$$\begin{equation} \rho _\nu = N_{\rm eff}\frac{7\pi ^2}{120}T_\nu ^4, \end{equation} \tag{ 82 }$$

where T_ν is the temperature of neutrinos. How do we measure N_eff from CMB?

The way that we use CMB to determine N_eff is relatively simple. The relativistic particles that stream freely influence CMB in two ways: (1) their energy density changing the matter-radiation equality epoch and (2) their anisotropic stress acting as an additional source for the gravitational potential via Einstein's equations. Incidentally, the relativistic particles that do not stream freely, but interact with matter frequently, do not have a significant anisotropic stress because they isotropize themselves via interactions with matter; thus, anisotropic stress of photons before the decoupling epoch was very small. Neutrinos, on the other hand, decoupled from matter much earlier (∼2 MeV), and thus their anisotropic stress was significant at the decoupling epoch.

The amount of the early ISW effect changes as the equality redshift changes. The earlier the equality epoch is, the more the ISW effect CMB photons receive. This effect can be measured via the height of the third acoustic peak relative to the first peak. Therefore, the equality redshift, z_eq, is one of the fundamental observables that one can extract from the CMB power spectrum.

One usually uses z_eq to determine Ω_mh² from the CMB power spectrum, without noticing that it is actually z_eq that they are measuring. However, the conversion from z_eq to Ω_mh² is automatic only when one knows the radiation content of the universe exactly—in other words, when one knows N_eff exactly:

$$\begin{equation} 1+z_{\rm eq} = \frac{\Omega _m}{\Omega _{r}} =\frac{\Omega _mh^2}{\Omega _\gamma h^2}\frac{1}{1+0.2271N_{\rm eff}}, \end{equation} \tag{ 83 }$$

where Ω_γh² = 2.469 × 10⁻⁵ is the present-day photon energy density parameter for T_cmb = 2.725 K. Here, we have used the standard relation between the photon temperature and neutrino temperature, T_ν = (4/11)^1/3T_γ, derived from the entropy conservation across the electron–positron annihilation (see, e.g., Weinberg 1972; Kolb & Turner 1990).

However, if we do not know N_eff precisely, it is not possible to use z_eq to measure Ω_mh². In fact, we lose our ability to measure Ω_mh² from CMB almost completely if we do not know N_eff. Likewise, if we do not know Ω_mh² precisely, it is not possible to use z_eq to measure N_eff. As a result, N_eff and Ω_mh² are linearly correlated (degenerate), with the width of the degeneracy line given by the uncertainty in our determination of z_eq.

The distance information from BAO and SN provides us with an independent constraint on Ω_mh², which helps to reduce the degeneracy between z_eq and Ω_mh².

The anisotropic stress of neutrinos also leaves distinct signatures in the CMB power spectrum, which is not degenerate with Ω_mh² (Hu et al. 1995; Bashinsky & Seljak 2004). Trotta & Melchiorri (2005; also see Melchiorri & Serra 2006) have reported on evidence for the neutrino anisotropic stress at slightly more than 95% CL. They have parametrized the anisotropic stress by the viscosity parameter, c²_vis (Hu 1998), and found c²_vis > 0.12 (95% CL). However, they had to combine the WMAP 1-year data with the SDSS data to see the evidence for nonzero c²_vis.

In Dunkley et al. (2009), we reported on the lower limit to N_eff solely from the WMAP 5-year data. In this paper, we shall combine the WMAP data with the distance information from BAO and SN as well as Hubble's constant from HST to find the best-fitting value of N_eff.

Figure 18 shows our constraint on N_eff. The contours in the left panel lie on the expected linear correlation between Ω_mh² and N_eff given by

$$\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{ 84 }$$

which follows from Equation (83). Here, Ω_mh² = 0.1308 and z_eq = 3138 are the maximum likelihood values from the simplest ΛCDM model. The width of the degeneracy line is given by the accuracy of our determination of z_eq, which is given by z_eq = 3141⁺¹⁵⁴₋₁₅₇ (WMAP-only) for this model. Note that the mean value of z_eq for the simplest ΛCDM model with N_eff = 3.04 is z_eq = 3176⁺¹⁵¹₋₁₅₀, which is close. This confirms that z_eq is one of the fundamental observables, and N_eff is merely a secondary parameter that can be derived from z_eq. The middle panel of Figure 18 shows this clearly: z_eq is determined independently of N_eff. For each value of N_eff along a constant z_eq line, there is a corresponding Ω_mh² that gives the same value of z_eq along the line.

Zoom In Zoom Out Reset image size

However, the contours do not extend all the way down to N_eff = 0, although Equation (84) predicts that N_eff should go to zero when Ω_mh² is sufficiently small. This indicates that we see the effect of the neutrino anisotropic stress at a high significance. While we need to repeat the analysis of Trotta & Melchiorri (2005) in order to prove that our finding of N_eff > 0 comes from the neutrino anisotropic stress, we believe that there is a strong evidence that we see nonzero N_eff via the effect of neutrino anisotropic stress, rather than via z_eq.

While the WMAP data alone can give a lower limit on N_eff (Dunkley et al. 2009), they cannot give an upper limit owing to the strong degeneracy with Ω_mh². Therefore, we use the BAO, SN, and HST data to break the degeneracy. We find N_eff = 4.4 ± 1.5 (68%) from WMAP+BAO+SN+HST, which is fully consistent with the standard value, 3.04 (see the right panel of Figure 18).