BAYESIAN ANALYSIS OF SPARSE ANISOTROPIC UNIVERSE MODELS AND APPLICATION TO THE FIVE-YEAR WMAP DATA

From Equation (2) and the definition of ξ (see ACW), it is clear that Δ also contributes to the diagonal of the signal covariance matrix, and therefore affects the total angular power spectrum, not only the correlations among a_ℓm's. This introduces a strong degeneracy between g_* and the amplitude of the power spectrum of scalar perturbations, A_s or σ₈. Unless one attempts to estimate the standard ΛCDM parameters jointly with the new anisotropy parameters, one must therefore ensure that a given choice of g_* does not significantly affect the overall power spectrum, but only the anisotropic contribution.

The diagonal part of Δ, for which the integral over the transfer functions equals C_ℓ, is

$$\begin{equation} \Delta _{\ell m, \ell m} = g_* C_{\ell } \xi _{\ell m, \ell m}. \end{equation} \tag{ 5 }$$

The net extra power due to Δ is therefore

$$\begin{equation} D_{\ell } = \frac{g_*C_{\ell }}{2\ell +1} \sum _{m=-\ell }^{\ell } \xi _{\ell m, \ell m}. \end{equation} \tag{ 6 }$$

This may be greatly simplified by considering the detailed form of ξ_ℓm,ℓm,

$$\begin{equation} \xi _{\ell m, \ell m} = -2 n_+ n_- \frac{-1+\ell +\ell ^2+m^2}{(2\ell -1)(2\ell +3)} + n_0^2 \frac{2\ell ^2+2\ell -2m^2-1}{(2\ell -1)(2\ell +3)}, \end{equation} \tag{ 7 }$$

where

$$\begin{equation} n_+ = -\frac{n_x -in_y}{\sqrt{2}};\qquad n_- = \frac{n_x+in_y}{\sqrt{2}};\qquad n_0 = n_z. \end{equation} \tag{ 8 }$$

Averaging this expression over m, one finds that

$$\begin{equation} \frac{1}{2\ell +1} \sum _{m=-\ell }^{\ell } \xi _{\ell m, \ell m} = \frac{1}{3}, \end{equation} \tag{ 9 }$$

such that $D_{\ell } = \frac{1}{3} g_* C_{\ell }$. We therefore redefine the total signal covariance matrix to read

$$\begin{equation} S_{\ell m, \ell ^{\prime } m^{\prime }} = \frac{C_{\ell }\delta _{\ell \ell ^{\prime }} \delta _{m m^{\prime }} + \Delta _{\ell m, \ell ^{\prime } m^{\prime }} }{1+g_*/3}. \end{equation} \tag{ 10 }$$

With this definition, g_* is a direct measure of the anisotropic component of S, and does not directly depend on the power spectrum C_ℓ.

The effect of g_* on the power spectrum is demonstrated in Figure 3, where we plot the power spectra of a simulated anisotropic map with g_* = 3, with and without the above rescaling. Unless proper rescaling is performed, or some equivalent parametrization introduced, it is clear that the strongest constraints on g_* will come from the observed power spectrum rather than the correlations among a_ℓm's.

Zoom In Zoom Out Reset image size

The goal is now to estimate g_* and $\hat{\mathbf{n}}$ from observed CMB maps, by computing the posterior distribution $P(g_*, \hat{\mathbf{n}}|\mathbf{d})$, with d denoting the data. Because we assume that both the noise and CMB sky signal are Gaussian (but anisotropic) random fields, this distribution reads, by Bayes' theorem,

$$\begin{equation} P(g_*, \hat{\mathbf{n}}|\mathbf{d}) \propto \mathcal{L}(g_*, \hat{\mathbf{n}}) P(g_*, \hat{\mathbf{n}}), \end{equation} \tag{ 11 }$$

where $\mathcal{L}(g_*, \hat{\mathbf{n}}) = P(\mathbf{d}|g_*, \hat{\mathbf{n}})$ is the likelihood

$$\begin{equation} \mathcal{L}(g_*, \hat{\mathbf{n}}) \propto \frac{e^{-\frac{1}{2}\mathbf{d}^T \mathbf{C}^{-1}\mathbf{d}}}{\sqrt{|\mathbf{C}|}} \end{equation} \tag{ 12 }$$

and $P(g_*, \hat{\mathbf{n}})$ is a prior. Equation (12) can be evaluated in $\mathcal{O}\big(N_{{\rm pix}}^2\big)$ operations, as shown by Oh et al. (1999). In this expression, C is the signal-plus-noise covariance matrix. In principle, we could now simply map this distribution over a three-dimensional grid and our task would be completed. However, except for the special case of a data set with uniform noise and full-sky coverage, this is impossible in practice because C is a dense matrix, and inversion and matrix determinant therefore scales as $\mathcal{O}\big(N_{{\rm pix}}^3\big)$, N_pix being the number of pixels. For current and future data sets, one expects N_pix ∼ 10⁶ or more.

Fortunately, there is one specific feature of the ACW model that does make an exact analysis possible: although the full-sky CMB covariance matrix is nondiagonal, it is not dense. Rather, it has a well-defined shape in harmonic space (ℓ is coupled to ℓ' = {ℓ, ℓ ± 2} and m to m' = {m, m ± 1, m ± 2}) that allows for cheap matrix storage and fast Cholesky decomposition. This, combined with the development of the standard diagonal CMB Gibbs sampler mentioned in the introductory section (Jewell et al. 2004; Wandelt et al. 2004; Eriksen et al. 2004b), allows us to perform a full proper analysis, as explained in the next section.

Before describing this method, we consider first the special case of data having uniform noise and full-sky coverage, which is useful to illustrate the approach and highlight some particular issues. For this particular case, the full data covariance matrix, expressed in spherical harmonic space, has the same sparse filling pattern as the ACW covariance matrix, and direct evaluation is therefore possible using sparse matrix techniques (e.g., Davis 2005).

We simulated a single CMB realization from the ACW model, adopting a high anisotropy amplitude of g_* = 0.8 and a preferred direction (in Galactic longitude and latitude) of (l, b) = (57°, 33°), then convolved this realization with a 90' FWHM Gaussian beam, and projected it onto a HEALPix³ grid with resolution parameter N_side = 128. Finally, uniform, Gaussian noise with 10 μK rms was added to each pixel. This simulation was then analyzed by computing the raw likelihood over a three-dimensional grid, and finally marginalized likelihoods were produced by numerical integration. The results from this exercise are shown in Figure 4.

Zoom In Zoom Out Reset image size

As expected, the likelihood peaks close to the input values. However, there is also a second local maximum at g_* ∼ −0.5 with a direction of (l, b) ∼ (45°, − 50°), 90° with respect to the main axis. This maximum becomes visible only for large negative values of g_*. The existence of this maximum becomes intuitive when considering Figure 2: Flipping the sign of g_* and rotating the preferred axis by 90° leads to stripes in the same direction as the original parameters.

This is not a significant issue for a direct evaluation method, since the local maximum has a very small amplitude (note that the marginal likelihoods in Figure 4 are shown in logarithmic units). However, for Monte Carlo Markov chain (MCMC) methods it can cause problem in terms of burn-in; as explained in the next section, our method is based on the well-known MCMC and Gibbs sampling algorithms, and these essentially correspond to performing a random walk on the likelihood surface. Further, each chain is initialized randomly on the sphere. It is therefore a significant chance that a number of chains may get trapped in a local maximum, and thereby bias the final posterior. To avoid this, we impose a uniform prior of g_* ⩾ −0.2 in this paper, and a uniform prior on the sphere for $\hat{\mathbf{n}}$. If the final posteriors from the WMAP analysis actually happened to peak close to g_* = −0.2, we would have to reconsider this choice more carefully, but as we shall see, this is not the case.

We first review the CMB Gibbs sampler as previously described in literature. In any Bayesian analysis, the main goal is the posterior distribution P(θ|d), where θ is a set of parameters connected to some model and d is the observed data. For high-dimensional spaces, brute-force evaluations of the posterior are computationally infeasible, and one usually resorts to MCMC methods.

3.1.1. Notation and Data Model

We begin by defining a parametric model for the CMB observations. Given our current understanding of the CMB sky, the observed data may be accurately modeled as a sum of a CMB anisotropy term and a noise term,

$$\begin{equation} \mathbf{d} = \mathbf{A}\mathbf{s} + \mathbf{n}. \end{equation} \tag{ 13 }$$

Here d represents the observed data, A denotes convolution by an instrumental beam, s(θ, ϕ) = ∑_ℓ,ma_ℓmY_ℓm(θ, ϕ) is the CMB sky signal represented in either harmonic or real space, and n is instrumental noise.

Further, it is a good approximation to assume both the CMB and noise to be zero-mean Gaussian-distributed variates, with covariance matrices S and N, respectively. In harmonic space, the signal covariance matrix is defined by $\textbf S_{\ell m,\ell ^{\prime }m^{\prime }} = \langle a_{\ell m} a_{\ell ^{\prime }m^{\prime }}^*\rangle$, which may or may not be diagonal. The connection to cosmological parameters θ is made through this covariance matrix. Finally, for experiments such as WMAP, the noise is often assumed to be uncorrelated between pixels, N_ij = σ²_iδ_ij, for pixels i and j, and noise rms equals σ_i.

Our goal is now to compute the full joint posterior P(θ|d), which, as already mentioned, is given by $P(\theta |\mathbf{d}) \propto P(\mathbf{d}| \theta) P(\theta) = \mathcal{L}(\theta) P(\theta),$ where $\mathcal{L}(\theta)$ is the likelihood, and P(θ) is a prior. For a Gaussian data model, the likelihood is

$$\begin{equation} \mathcal{L}(\theta) \propto \frac{e^{-\frac{1}{2}\mathbf{d}^T \mathbf{C}^{-1}(\theta)\mathbf{d}}}{\sqrt{|\mathbf{C(\theta)}|}}. \end{equation} \tag{ 14 }$$

3.1.2. Posterior Mapping by Gibbs Sampling

When working with real-world CMB data, there are a number of issues that complicate the analysis. Two important examples are anisotropic noise and Galactic foregrounds. First, because of the scanning motion of a CMB satellite, the pixels in a given data set are observed by unequal amounts of time. This implies that the effective noise is a function of position on the sky. Second, large regions of the sky are obscured by Galactic foregrounds (e.g., synchrotron, free–free, and dust emission), and these regions must be rejected from the analysis by masking.

Because of these issues, the total data covariance matrix S + N is dense in both pixel and harmonic space. As a result, it is computationally difficult to evaluate the likelihood in Equation (14), since the computational cost of matrix inversion and determinant evaluation scale as $\mathcal{O}\big(N_{{\rm pix}}^3\big)$. Fortunately, this problem has already been solved for the CMB context, through the development of the CMB Gibbs sampler.

The idea behind the CMB Gibbs sampler is to estimate the CMB sky, s, together with the covariance parameters, by computing P(θ, s|d), and then subsequently marginalize over s. Specifically, the algorithm is the following. First choose any initial guess, (θ, s)⁰. Then alternately sample from each of the conditional distributions

$$\begin{equation} \theta ^{i+1} \leftarrow P(\theta | \mathbf{s}^i, \mathbf{d}), \end{equation} \tag{ 15 }$$

$$\begin{equation} \mathbf{s}^{i+1} \leftarrow P(\mathbf{s} | \theta ^{i+1}, \mathbf{d}). \end{equation} \tag{ 16 }$$

The theory of Gibbs sampling then guarantees that the joint samples (θ, s)ⁱ will, after some burn-in period, be drawn from the desired joint distribution. The remaining step is then simply to formulate sampling algorithms for each of the two conditionals, P(θ|s, d) and P(s|θ, d).

We first consider P(s|θ, d). This may, under the assumption of Gaussianity, be written as

$$\begin{equation} P(\mathbf{s} | \theta, \mathbf{d}) = P(\mathbf{d}|\mathbf{s}, \theta) P(\mathbf{s}|\theta) \end{equation} \tag{ 17 }$$

$$\begin{equation} \propto e^{-\frac{1}{2} (\mathbf{d}-\mathbf{s})^T \mathbf{N}^{-1} (\mathbf{d}-\mathbf{s})} e^{-\frac{1}{2}\mathbf{s}\mathbf{S}^{-1} \mathbf{s}} \end{equation} \tag{ 18 }$$

$$\begin{equation} = e^{-\frac{1}{2}(\mathbf{s}-\hat{\mathbf{s}})^T (\mathbf{S}^{-1} + \mathbf{N}^{-1}) (\mathbf{s}-\hat{\mathbf{s}})}, \end{equation} \tag{ 19 }$$

where we have defined the Wiener filtered map, $\hat{\mathbf{s}} = (\mathbf{S}^{-1} + \mathbf{N}^{-1})^{-1}\mathbf{N}^{-1}\mathbf{d}$. Thus, P(s|θ, d) is a Gaussian distribution with mean $\hat{\mathbf{s}}$ and covariance (S⁻¹ + N⁻¹)⁻¹.

Sampling from this distribution is straightforward, but implementationally somewhat involved. Draw two Gaussian random maps, η₀ and η₁, with zero mean and unit variance, and solve the following equation for $\bar{\mathbf{s}}$:

$$\begin{equation} (\mathbf{S}^{-1} + \mathbf{N}^{-1})\bar{\mathbf{s}} = \mathbf{N}^{-1}\mathbf{d} + \mathbf{L}^{-T} \eta _0 + \mathbf{N}^{-\frac{1}{2}}, \end{equation} \tag{ 20 }$$

where L is the Cholesky decomposition of S = LL^T. By multiplying both sides of this equation with (S⁻¹ + N⁻¹)⁻¹, one immediately sees that $\left<\bar{\mathbf{s}}\right> = \hat{\mathbf{s}}$, and a few more computations show that $\langle(\bar{\mathbf{s}}-\tilde{\mathbf{s}})(\bar{\mathbf{s}}-\tilde{\mathbf{s}})^T\rangle = (\mathbf{S}^{-1} + \mathbf{N}^{-1})^{-1}$, as required.

For improved numerical stability, this linear system is in practice rewritten into the following form:

$$\begin{equation} (\mathbf{1} + \mathbf{L}^{T}\mathbf{N}^{-1}\mathbf{L})(\mathbf{L}^{-1} \bar{\mathbf{s}}) = \mathbf{L}^{T}\mathbf{N}^{-1}\mathbf{d} + \eta _0 + \mathbf{L}^{T}\mathbf{N}^{-\frac{1}{2}}, \end{equation} \tag{ 21 }$$

which is first solved for $\mathbf{x} = \mathbf{L}^{-1} \bar{\mathbf{s}}$ by conjugate gradients, and then for $\bar{\mathbf{s}} = \mathbf{L} \mathbf{x}$. For further implementational details, see, e.g., Eriksen et al. (2008a). Note, however, that in previous papers Equation (21) was always written with symmetric signal covariance square roots, $\mathbf{S}^{\frac{1}{2}} = \big(\mathbf{S}^{\frac{1}{2}}\big)^T$. The current form is based on the Cholesky decomposition, which is computationally considerably cheaper than the symmetric form, especially for sparse matrices.

Finally, we need a sampling algorithm for P(θ|s, d). In previous publications, the main emphasis has been on covariance matrices parametrized by the angular CMB power spectrum, $C_{\ell m, \ell ^{\prime } m^{\prime }} = C_{\ell } \delta _{\ell \ell ^{\prime }} \delta _{m m^{\prime }}$. In this case, P(C_ℓ|s, d) reduces to a simple inverse gamma distribution, for which there is a simple textbook sampling algorithm available. As the details of this specific algorithm is of little use for the application presented here, we refer the interested reader to earlier papers for full details on this procedure, e.g., Wandelt et al. (2004) or Eriksen & Wehus (2008).

We now describe the two modifications to the CMB Gibbs sampler that allow us to analyze models with nondiagonal covariances. This involves adding support for nondiagonal covariance matrices for P(s|θ, d) and implementing a more general sampling algorithm for P(θ|s, d).

3.2.1. Sampling from P(s|θ, d)

We first consider sampling of sky maps, s, given a set of cosmological parameters, θ, and the associated covariance matrix S(θ). Formally, the sampling algorithm for P(s|θ, d) is identical to that given by Equation (21). However, in this case S is a nondiagonal matrix and the computational complexity is therefore greatly increased. Only special cases can be considered, for instance, models that predict a sparse covariance matrix. This is the case for the ACW model.

For general dense anisotropic covariance matrices, the memory requirements scale as $\mathcal{O}\big(\ell _{{\rm max}}^4\big)$, effectively rendering studies of anisotropic models where ℓ_max ≳ 100 impossible. However, working only with sparse matrices, the memory consumption scales as $\mathcal{O}\big(\ell _{{\rm max}}^2\big)$, enabling calculations of covariance matrices with ℓ_max well into the Planck regime (ℓ_max ∼ 2500).

To be able to handle sparse matrices efficiently, we have ported the LDL library of Davis (2005) to Fortran 90, and incorporated this into Commander. This library stores sparse matrices in a packed format, and supports fast Cholesky decomposition. Our F90 version of LDL may be obtained by sending an email to the authors, and it will be released publicly at a later time.

In the present paper, we are primarily concerned with the ACW model, and the corresponding covariance matrix exhibits correlations between ℓ and ℓ' = {ℓ, ℓ ± 2} and between m and m' = {m, m ± 1, m ± 2}. Thus, the number of elements up to ℓ_max is $\mathcal{O}\big(15\ell _{{\rm max}}^2\big)$. For example, for ℓ_max = 300 the memory requirements are ∼14 Mb with double precision complex numbers. Since the covariance matrix is very sparse, the CPU time required for Cholesky decomposition is nearly linear in ℓ²_max.

We define three different harmonic space limits in our code, namely ℓ_max, ℓ_low, and ℓ_high. The former denotes the maximum multipole moment of the full spherical harmonics composition used in the analysis, while the latter two denotes the range in which the anisotropic covariance matrix is used. In addition, we remove the monopole and dipole from the analysis. Thus, the total covariance matrix reads

$$\begin{eqnarray} &&S_{\ell m, \ell ^{\prime } m^{\prime }} \nonumber\\ &&=\left\lbrace \begin{array}{@{}ll@{}} 0 &\!\! \ell, \ell ^{\prime } \le 1 \\ C_{\ell }\delta _{\ell \ell ^{\prime }} \delta _{m m^{\prime }} &\!\! 2 \le \ell,\ell ^{\prime } < \ell _{{\rm low}} \\ (C_{\ell }\delta _{\ell \ell ^{\prime }} \delta _{m m^{\prime }} + \Delta _{\ell m, \ell ^{\prime } m^{\prime }}) / (1+g_*/3) &\!\! \ell _{{\rm low}} \le \ell, \ell ^{\prime } \le \ell _{{\rm high}} \\ C_{\ell }\delta _{\ell \ell ^{\prime }} \delta _{m m^{\prime }} & \!\!\ell _{{\rm high}} < \ell,\ell ^{\prime } \le \ell _{{\rm max.}} \end{array} \right.\nonumber\\ \end{eqnarray} \tag{ 22 }$$

The reason for defining ℓ_min and ℓ_max as free parameters is that it may be useful to study the dependence of θ on a particular ℓ-range. On the other hand, this implies that the model implemented in this paper is only an approximation to the full ACW model, for which the correlations extend over all ℓs, and is only exact when ℓ_high = ℓ_max.

With the sparse matrix operations implemented, the algorithm is precisely the same as for the diagonal case, and both rely on the solution of a linear system by conjugate gradients (CG; Eriksen et al. 2004b). In order to achieve an acceptable CG convergence rate, it is therefore necessary to establish a good preconditioner. However, as long as the off-diagonal elements remain small, the standard diagonal covariance matrix preconditioner performs reasonably even for the off-diagonal case. For the present paper, we therefore adopt the same preconditioner as described by Eriksen et al. (2004b), which consists of the directly inverted full matrix evaluated up to some ℓ_precond, and then a strictly diagonal matrix from ℓ_precond + 1 to ℓ_max.

The number of CG iterations per map making step is typically 70 for a WMAP-type run, and with a total CPU time per iteration of about 15 s, the total cost for a single sample is ∼20 CPU minutes. The average CPU time required to set up and perform a Cholesky decomposition of the corresponding covariance matrix for ℓ_max = 512, ℓ_low = 2, and ℓ_high = 300 is ∼20 s.

3.2.2. Sampling from P(θ|s, d)

Finally, we have to formulate a sampling algorithm for P(θ|s, d). Recall that for the diagonal power spectrum case, this step is typically performed by a standard inverse gamma distribution sampler (e.g., Gupta & Nagar 2000; Eriksen & Wehus 2008). For the general case considered here, we adopt a standard Metropolis MCMC sampler (e.g., Liu 2001).

First, note that P(θ|s, d) = P(θ|s); if we already know the CMB sky perfectly, no additional data can possibly tell us anything more about the anisotropy parameters θ. Second, although the CMB sky is now manifestly anisotropic, we still assume that it is Gaussian, and the target distribution therefore reads

$$\begin{equation} P(\theta |\mathbf{s}) \propto \frac{e^{-\frac{1}{2}\mathbf{s}^T \mathbf{S}^{-1}\mathbf{s}}}{\sqrt{|\mathbf{S}|}}. \end{equation} \tag{ 23 }$$

For sparse matrices, this may be directly evaluated by first computing the Cholesky decomposition of S = LL^t, and then, on the one hand, solve for x = Ls, and on the other hand, compute |S| = |L|².

We adopt a simple symmetric proposal rule for the Metropolis sampler, and the acceptance probability therefore simply reads

$$\begin{equation} p = \frac{P(\theta ^p|\mathbf{s})}{P(\theta ^i|\mathbf{s})}, \end{equation} \tag{ 24 }$$

where θ^p is the proposed sample and θⁱ is the current sample of the MCMC chain. Specifically, we adopt a Gaussian proposal density for g_* and a uniform proposal over a disk for $\hat{\mathbf{n}}$, centered on the current state. The proposal density is typically tuned by producing a short test chain before the main run, such that the final observed acceptance rate lies between 0.2 and 0.7.

Finally, because the computational cost is much lower for this step than for P(s|θ, d), we produce several θ samples per main Gibbs iteration, to improve the convergence properties of the chain. This essentially corresponds to performing a partial Rao–Blackwellization (Chu et al. 2005). A typical number of MCMC samples per main Gibbs iteration is 30.

To test our implementation and study the behavior of the algorithm in general, we simulate a few different maps from the ACW model, and analyze these maps with our modified Gibbs sampler. The CMB component of these maps is made by generating a random vector, η, of Gaussian uniform variates with zero mean and unit variance, and then computing s = Lη. This realization is then convolved with a beam function and the HEALPix pixel window, before it is projected on a HEALPix grid. Finally, Gaussian noise is added to each pixel.

The first two simulations have a resolution of N_side = 128, ℓ_max = 256, ℓ_low = 2, ℓ_high = 200, and a Gaussian beam of 90' FWHM. The noise rms is 10 μK uniformly over the full sky. The first of the two simulations has an anisotropy amplitude of g_* = 0.8 and a preferred direction toward (l, b) = (57°, 33°), and the other g_* = 0. These two simulations are primarily used to compare the Gibbs sampler with brute-force likelihood evaluation, which is only possible for uniform noise and full-sky coverage.

Second, we generate a full WMAP5-like simulation based on the V1 differencing assembly (DA), with g_* = 0.8, N_side = 512, ℓ_max = 600, ℓ_low = 2, ℓ_high = 300, and beam and noise properties appropriate for the V1 DA.⁴ In this case, we also apply the KQ85 sky cut (Gold et al. 2008), which removes 18% of the sky. This simulation is used to verify that correct results are obtained for realistic WMAP data, including anisotropic (but uncorrelated) noise and a sky cut.

We first consider the issue of burn-in and convergence, and analyze the simulation with g_* = 0.8, uniform noise, and full-sky coverage. In Figure 5, we show the first 6000 g_* samples produced by each of 14 chains. First, notice that the chains immediately divide into two classes, one which converges quickly toward g_* ∼ 0.8, and the other which hovers near the lower prior of g_* = −0.2. This is due to the fact that the chains are initialized randomly on the sphere, and those that happen to start close to the nonphysical local maximum (see Section 2.2) get temporarily trapped in this local maximum. However, as the chains explores the likelihood surface, they are able to converge into the right regime, and find the correct value. In this case, all chains have reached the equilibrium state after 1800 iterations. The pre–burn-in samples must be rejected from further analysis. For now, we inspect each chain individually, to make sure that they have all reached the common state.

Zoom In Zoom Out Reset image size

Note that there is a fundamental difference between low and high signal-to-noise cases in this respect: if g_* is low, the chains may jump between local maxima, while if g_* is high, some chains typically start out in the global maximum and stay there, while others start in the local maximum and eventually converge into the right regime. Which situation is relevant for a particular data set must be considered on a case-by-case basis, by checking whether the chains jump between states or if they stay in one place. It is also advisable to run many chains in parallel, randomly initialized over the full sphere, to understand how many local maxima the distribution has.

Second, once the chains have burned in, we must also ensure that they collectively have converged to the full posterior. One possible measure for this is the Gelman–Rubin R statistic (e.g., Gelman & Rubin 1992; Eriksen et al. 2006), which compares the variances within a single chain with the variance between chains. If the chains have converged properly, R should be close to unity. Typically, one recommends that R should be less than 1.1 or 1.2. For the chains shown in Figure 5, we find that R = 1.01 after rejecting the first 2000 burn-in samples, indicating very good convergence. Considering further subsets of these samples, we find that 5000 samples are sufficient to achieve R < 1.1 and 20,000 samples for R < 1.02. In the analyses presented later, we always have more than 20,000 samples.

We now analyze the two simulations with uniform noise and full-sky coverage, having g_* = 0 and g_* = 0.8, respectively. In addition to running the Gibbs sampler on these simulations, we also compute the full three-dimensional likelihood function over a grid in $(g_*, \hat{\mathbf{n}})$, and numerically integrate to produce brute-force marginal posteriors.

The resulting distributions are shown in Figure 6; the left column shows the g_* = 0.8 case, and the right column shows the g_* = 0 case. We see, as expected, that the two methods produce identical results, up to sampling uncertainty and grid resolution. Note that this holds both for high- and low-anisotropy amplitudes, indicating that the method is robust in all regimes.

Zoom In Zoom Out Reset image size

Next, we see that when the amplitude is large, there is only one visible preferred direction in the direction posterior; the secondary direction is too shallow to be seen. On the other hand, there are two "preferred" directions in the g_* = 0 case. However, the span in likelihood over the full sphere is in this case only a factor of 2 between the least and the most preferred directions, which essentially indicates a uniform distribution.

In Figure 7, we show similar plots for the WMAP simulation with uncorrelated noise, based on the V1 DA and g_* = 0.8. In this case it is not possible to evaluate the likelihood directly, since the noise is inhomogeneous and there is a sky cut. Still, we see that correct results are obtained. This concludes the verification of both the method and our implementation, and we are now ready to analyze the five-year WMAP temperature sky maps.

Zoom In Zoom Out Reset image size

Before turning to the analysis of the actual WMAP data, we compute the uncertainty in g_* as a function of ℓ_high for full-sky noiseless data (the lower limit is always kept at ℓ_low = 2). This topic was also considered by Pullen & Kamionkowski (2007), who presented both a more general formalism and forecasts for specific experiments. Note, however, that our parametrization is slightly different from theirs, as we introduce a rescaling of the covariance matrix to eliminate the power spectrum degeneracy (Section 2.1).

We carry out this analysis by simulating anisotropic ACW maps with g_* = 0 and different ℓ_high, and analyze these with the brute-force evaluation approach described above. No noise or beam effects are included. For each case, we marginalize over $\hat{\mathbf{n}}$ to obtain P(g_*|d), and compute the standard deviation, σ(g_*), from this distribution.

Figure 8 shows σ(g_*) as a function of ℓ_high. From this figure, we see that σ(g_*) is very close to a power law in ℓ_high, in good agreement with the arguments given by Pullen & Kamionkowski (2007). The best-fit power-law function is

$$\begin{equation} \sigma (\ell _{{\rm high}}; g_*) = 0.025\left(\frac{\ell _{{\rm high}}}{400}\right)^{-1.27}, \end{equation} \tag{ 25 }$$

and this can be used to produce rough forecasts for various experiments. For instance, in this paper we conservatively adopt ℓ_high = 400 for the WMAP analysis, to avoid possibly complicating high-ℓ issues, such as point-source confusion and noise misestimation. In that case, we expect an uncertainty of σ(g_*) = 0.025 before taking into account noise and sky cut. This is in excellent agreement with the σ(g_*) = 0.024 result of Pullen & Kamionkowski (2007), derived with slightly different data and model assumptions and a completely different approach.

Zoom In Zoom Out Reset image size

In this paper, we consider the five-year WMAP temperature sky maps (Hinshaw et al. 2008), and analyze the V and W bands (61 and 94 GHZ), which are believed to be the cleanest WMAP bands in terms of residual foregrounds. We adopt the template-corrected, foreground reduced maps recommended by the WMAP team for cosmological analysis, and impose both the KQ75 and KQ85 masks (Gold et al. 2008), which remove 28% and 18% of the sky, respectively. Point-source cuts are imposed in both cases.

We mainly analyze the data frequency-by-frequency, and consider the combinations V1+V2 and W1 through W4. In addition, we compute the posteriors for V1 and V2 separately. The noise rms patterns and beam profiles are taken into account for each DA individually. The noise is assumed uncorrelated between pixels and bands. For details on joint Gibbs analysis of multifrequency data, see Eriksen et al. (2004b). All data used in this analysis are available from Legacy Archive for Microwave Background Data Analysis (LAMBDA).

The angular resolutions of the V and W bands are 035 and 022, respectively, and the sky maps are pixelized at a HEALPix resolution of N_side = 512 with 7' pixels. We therefore adopt a harmonic space cutoff of ℓ_max = 700 and 800 for the two data sets, probing deeply into the noise-dominated regime. However, we never consider multipoles at ℓ>400 for the anisotropic part of the signal covariance matrix, in order to minimize the chance of systematic effects such as residual point-source contributions, beam uncertainties, or noise misestimation to affect our results. See Table 1 for a list of the specific ℓ-ranges considered.

We also note that the maps studied here are cleaned using external templates (Gold et al. 2008), which must be considered a fairly rough approach to foreground cleaning. A better approach is to use the joint foreground and the CMB Gibbs sampler (Eriksen et al. 2008a), which provides the user with a CMB map marginalized over very general foreground models. This work is currently underway for the five-year WMAP data, and the results will be reported elsewhere (C. Dickinson et al. 2008, in preparation). However, as an explicit foreground test, we also analyze the raw V-band data, from which no foreground templates have been subtracted, and find very consistent results.

We now present the marginal posteriors for the ACW model obtained from the five-year WMAP temperature sky maps, as computed with the method described in Section 3. First, in the top row of Figure 9 we show the marginal anisotropy amplitude posterior, P(g_*|d), for the V (left column) and W (right column) bands, and in the three bottom rows we show the preferred direction posteriors, $P(\hat{\mathbf{n}}|\mathbf{d})$. In Table 1, the full set of results are summarized quantitatively.

Zoom In Zoom Out Reset image size

First, we see that there is an apparently clear detection of g_* ≠ 0 when considering the full range of multipoles, ℓ = [2, 400]. The W-band posterior has g_* = 0.15 ± 0.04, nominally corresponding to a 3.8σ detection, and the V-band posterior has g_* = 0.10 ± 0.04, internally consistent with W band at ∼1σ. Second, the direction posteriors indicate a clearly preferred direction pointing toward (l, b) = (110°, 10°).

Further, this same direction is observed in both ℓ = [2 − 100] and ℓ = [100, 400], indicating that the structure is present over a large range of angular scales. The results are also stable with respect to sky cut, as the same pattern is seen with the KQ75 sky mask as with the KQ85 cut, removing an additional 10% of the sky.

Given the nominally strong results found in the previous section, it is imperative to search for possible systematic effects that might explain the observations. In particular, three major sources of uncertainty should be considered in detail, namely noncosmological foregrounds, correlated noise, and asymmetric beams.

First, residual Galactic foregrounds do not appear a priori as a particularly promising candidate, given that the results are robust with respect to both frequency and sky cut, and the preferred axis does not point toward any natural Galactic direction. Second, in Figure 10 we show the posteriors obtained from the raw V-band five-year sky maps. It should be clear that Galactic foregrounds have little impact on these results. Third, extragalactic point sources do not also seem to be likely candidate, because the signature is seen both at low and high ℓs, and we never consider multipoles above ℓ>400, precisely to avoid these type of concerns.

Zoom In Zoom Out Reset image size

The effect of correlated noise is harder to rule out. On the one hand, returning to the simulated anisotropic CMB realization in Figure 2, we see that the main signature of the ACW model is smoothed structures along the plane normal to the preferred direction, and essentially no modifications along the preferred direction. On the other hand, the main signature of correlated noise is striping along the scan direction.

Next, the ecliptic north pole has Galactic coordinates (l, b) = (96°, 30°), which is ∼24° (32°) away from the preferred direction for W band (V band) found in Section 5.2. The probability of obtaining such a close alignment by chance is ∼10% (∼16%), which is low enough for correlated noise to be considered relevant for this particular case.

To study the magnitude of this effect on g_*, we analyze realistic V- and W-band simulations with correlated noise. These noise simulations were produced and published by the WMAP team in their one-year data release. To mimic realistic five-year simulations, we co-add five independent realizations for each DA. These noise realizations are then added to an isotropic CMB sky realization, and the sum is analyzed using the same procedure as in Section 5.2. We also analyze two single one-year W4-band simulations, which serve as a worst-case scenario, as the knee frequency of this band is by far the highest of any WMAP DA, and the overall noise level is higher by a factor of ∼4.5 than the full five-year W-band data.

The results from these analyses are shown in Figure 11. The left and middle columns show the simulated five-year posteriors for V and W bands, respectively, and the right column shows the one-year W4 posterior. The top row shows P(g_*|d), and the bottom panel shows $P(\hat{\mathbf{n}}|\mathbf{d})$.

Zoom In Zoom Out Reset image size

First, note that with realistic five-year noise no detection is made in either the V or the W band. Further, the peak sky position is different in the V and W bands, and both have a very low significance. It therefore seems unlikely that correlated noise can explain the results found in Section 5.2.

Still, caution is warranted, as the one-year W4 posteriors do exhibit traces of anisotropic contributions, with a peak amplitude larger than the observed g_* in the actual five-year data. Yet, the match with the structures observed in the real data is less than striking. First, the correlated noise simulations show two independent peaks in the directional posterior, while the WMAP data show one. Second, a detailed study of the joint posteriors for the simulations show that the peak along the ecliptic north pole corresponds to the negative peak in P(g_*|d), while the WMAP data have a positive g_* along their preferred direction. Third, the preferred axis found in the WMAP data are further away from the ecliptic pole than the corresponding peak in the simulation posteriors.

Finally, in order to make a complete analysis, we should also consider the impact of asymmetric beams. Ideally, one would prefer to address this issue in the same manner as correlated noise, by analyzing simulated CMB realizations with asymmetric beams. Unfortunately, we do not have access to such simulations at this time, and it is difficult to do a rigorous analysis. However, there are some arguments against the asymmetric beams hypothesis. First, the effect is observed both at low and high ℓ's, with very consistent positions. Second, the observed preferred axis is ∼25–30 degrees away from the ecliptic pole, and the posterior ratio of the ecliptic poles to the maximum posterior is low. Finally, similar signatures are observed in both the V and W bands, and in V1 and V2, which all have slightly different beam patterns.

Nevertheless, at this point it would be unwise to make strong claims concerning a possible cosmological interpretation of the signature found in Section 5.2. Proper analysis of fully realistic five-year WMAP simulations is required before one can attach cosmological significance to these findings.