Testing Gaussianity, Homogeneity, and Isotropy with the Cosmic Microwave Background

Abstract

We review the basic hypotheses which motivate the statistical framework used to analyze the cosmic microwave background, and how that framework can be enlarged as we relax those hypotheses. In particular, we try to separate as much as possible the questions of gaussianity, homogeneity, and isotropy from each other. We focus both on isotropic estimators of nongaussianity as well as statistically anisotropic estimators of gaussianity, giving particular emphasis on their signatures and the enhanced “cosmic variances” that become increasingly important as our putative Universe becomes less symmetric. After reviewing the formalism behind some simple model-independent tests, we discuss how these tests can be applied to CMB data when searching for large-scale “anomalies”.

1. Introduction

According to our current understanding of the Universe, the morphology of the cosmic microwave background (CMB) temperature field, as well as all cosmological structures that are now visible, like galaxies, clusters of galaxies, and the whole web of large-scale structure, are probably the descendants of quantum process that took place some 10⁻³⁵ seconds after the Big Bang. In the standard lore, the machinery responsible for these processes is termed cosmic inflation and, in general terms, what it means is that microscopic quantum fluctuations pervading the primordial Universe are stretched to what correspond, today, to cosmological scales (see [1–3] for comprehensive introductions to inflation.) These primordial perturbations serve as initial conditions for the process of structure formation, which enhance these initial perturbations through gravitational instability. The subsequent (classical) evolution of these instabilities preserves the main statistical features of these fluctuations that were inherited from their inflationary origin–provided, of course, that we restrain ourselves to linear perturbation theory.

However, given that matter has a natural tendency to cluster, and this inevitably leads to nonlinearities (not to mention the sorts of complications that come with baryonic physics), the structures which are visible today are far from ideal probes of those statistical properties. CMB photons, on the other hand, to an excellent approximation experience free streaming since the time of decoupling (𝑧≈1100) and are therefore exempt from these non-linearities (except, of course, for secondary anisotropies such as the Rees-Sciama effect or the Sunyaev-Zel'dovich effect), which implies that they constitute an ideal window to the physics of the early Universe—see, for example, [4–6]. In fact, we can determine the primary CMB anisotropies as well as most of the secondary anisotropies on large scales, such as the Integrated Sachs-Wolfe effect, completely in terms of the initial conditions by means of a linear kernel as follows:Θ(𝑛)̂≡Δ𝑇𝑛;̂𝜂⁰𝑇𝜂₀=𝑑³𝑥^^𝜂0₀𝑑𝜂^_𝑖𝐾_𝑖𝑥⃗^,𝜂^;𝑛̂𝑆^𝑖𝑥⃗^,𝜂^,(1) where 𝜂^ is conformal time, and 𝒮^𝑖 denote the initial conditions of all matter and metric fields (as well as their time derivatives, if the initial conditions are nonadiabatic). Here 𝐾_𝑖 is a linear kernel, or a retarded Green's function, that propagates the radiation field to the time and place of its detection, here on Earth. Since that kernel is insensitive to the statistical nature of the initial conditions (which can be thought of as constants which multiply the source terms), those properties are precisely transferred to the CMB temperature field Θ.

The statistical properties of the primordial fluctuations are, to lowest order in perturbation theory, quite simple; because the quantum fluctuations that get stretched and enhanced by inflation are basically harmonic oscillators in their ground state, the distribution of those fluctuations is Gaussian, with each mode an independent random variable. The Fourier modes of these fluctuations are characterized by random phases (corresponding to the random initial values of the oscillators), with zero mean, and variances which are given simply by the field mass and the mode's wavenumber 𝑘=2𝜋/𝜆. The presence of higher-order interactions (which exist even for free fields, because of gravity) changes this simple picture, introducing higher-order correlations which destroy gaussianity—even in the simplest scenario of inflation [7–9]. However, since these interactions are typically suppressed by powers of the factor 𝐺𝐻²≃10⁻¹², where 𝐺 is Newton's constant and 𝐻 the Hubble parameter during inflation, the corrections are small—but, at least in principle, detectable [10–12].

Since these statistical properties are a generic prediction of (essentially) all inflationary models, they can also be inferred from two ingredients that are usually assumed as a first approximation to our Universe. First, since inflation was designed to stretch our Universe until it became spatially homogeneous and isotropic, it is reasonable to expect that all statistical momenta of the CMB should be spatially homogeneous and rotationally invariant, regardless of their general form. Second, in linear perturbation theory [13] where we have a large number of cosmological fluctuations evolving independently, we can expect, based on the central limit theorem, that the Universe will obey a Gaussian distribution.

The power of this program lies, therefore, in its simplicity: if the Universe is indeed Gaussian, homogeneous, and statistically isotropic (SI), then essentially all the information about inflation and the linear (low redshift) evolution of the Universe is encoded in the variance, or two-point correlation function, of large-scale cosmological structures and/or the CMB. As it turns out, the five-year dataset from the Wilkinson Microwave anisotropy probe (WMAP) strongly supports these predictions [11, 14]. Moreover, the measurements of the CMB temperature power spectrum by the WMAP team, alongside measurements of the matter power spectrum from existing survey of galaxies [15, 16] and data from type Ia supernovae [17–19], have shown remarkable consistency with a concordance model (ΛCDM), in which the cosmos figures as a Gaussian, spatially flat, approximately homogeneous, and statistically isotropic web of structures composed mainly of baryons, dark matter, and dark energy.

However, while the detection of a nearly scale-invariant and Gaussian spectrum is a powerful boost to the idea of inflation, just knowing the variance of the primordial fluctuations is not sufficient to single out which particular inflationary model was realized in our Universe. For that, we will need not only the 2-point function, but the higher momenta of the distribution as well. Therefore, in order to break this model degeneracy, we must go beyond the framework of the ΛCDM, Gaussian, spatially homogeneous, and statistically isotropic Universe.

Reconstructing our cosmic history, however, is not the only reason to explore further the statistical properties of the CMB. The full-sky temperature maps by WMAP [11, 20] have revealed the existence of a series of large-angle anomalies—which, incidentally, were (on hindsight) already visible in the lower-resolution COBE data [21]. These anomalies suggest that at least one of our cherished hypotheses underlying the standard cosmological model might be wrong—even as a first-order approximation. Perhaps the most intriguing anomalies (described in more detail in other review papers in this volume) are the low value of the quadrupole and its alignment of the quadrupole (ℓ=2) with the octupole (ℓ=3) [22–27], the sphericity [26] (or lack of planarity [28]), of the multipole ℓ=5, and the north-south asymmetry [29–33]. In the framework of the standard cosmological model, these are very unlikely statistical events, and yet the evidence that they exist in the real data (and are not artifacts of poorly subtracted extended foregrounds—e.g., [34]) is strong.

Concerning theoretical explanations, even though we have by now an arsenal of ad hoc models designed to account for the existence of these anomalies, none has yet quite succeeded in explaining their origin. Nevertheless, they all share the point of view that the detected anomalies might be related to a deviation of gaussianity and/or statistical isotropy.

In this paper, we will describe, first, how to characterize, from the point of view of the underlying spacetime symmetries, both non-gaussianity and statistical anisotropy. We will adopt two guiding principles. The first is that gaussianity and SI, being completely different properties of a random variable, should be treated separately, whenever possible or practical. Second, since there is only one type of gaussianity and SI but virtually infinite ways away from them, it is important to try to measure these deviations without a particular model or anomaly in mind–although we may eventually appeal to particular models as illustrations or as a means of comparison. This approach is not new and, although not usually mentioned explicitly, it has been adopted in a number of recent papers [35, 36].

One of the main motivations for this model-independent approach is the difficult issue of aprioristic statistics; one can only test the random nature of a process if it can be repeated a very large (formally, infinite) number of times. Since the CMB only changes on a timescale of tens of millions of years, waiting for our surface of last scattering to probe a different region of the Universe is not a practical proposition. Instead, we are stuck with one dataset (a sequence of apparently random numbers), which we can subject to any number of tests. Clearly, by sheer chance, about 30% of the tests will give a positive detection with 70% confidence level (CL), 10% will give a positive detection with 90% CL, and so on. With enough time, anyone can come up with detections of arbitrarily high significance—and ingenuity will surely accelerate this process. Hence, it would be useful to have a few guiding principles to inform and motivate our statistical tests, so that we do not end up shooting blindly at a finite number of fish in a small wheelbarrow.

This paper is divided into two parts. We start Part I by reviewing the basic statistical framework behind linear perturbation theory (Section 2). This serves as a motivation for Section 3, where we discuss the formal aspects of non-Gaussian and statistically isotropic models (Section 2.1), as well as Gaussian models of statistical anisotropy (Section 2.2). Part II is devoted to a discussion on model-independent cosmological tests of non-gaussianity and statistical anisotropy and their application to CMB data. We focus on two particular tests, namely, the multipole vectors statistics (Section 3) and functional modifications of the two-point correlation function (Section 4). After discussing how such tests are usually carried out when searching for anomalies in CMB data (Section 6.1), we present a new formalism which generalizes the standard procedure by including the ergodicity of cosmological data as a possible source of errors (Section 6.2). This formalism is illustrated in Section 7, where we carry a search of planar-type deviations of isotropy in CMB data. We then conclude in Section 8.

Part I: The Linearized Universe

2. General Structure

We start by defining the temperature fluctuation field. Since the background radiation is known to have an average temperature of 2.725 K, we are interested only in deviations from this value at a given direction 𝑛̂ in the CMB sky. So let us consider the dimensionless function on 𝑆² as follows:Θ(𝑛)̂≡𝑇(𝑛)̂−𝑇⁰𝑇₀,(2) where 𝑇₀=2.725 K is the blackbody temperature of the mean photon energy distribution—which, if homogeneity holds, is also equal to the ensemble average of the temperature.

In full generality, the fluctuation field is not only a function of the position vector 𝑛⃗, but also of the time in which our measurements are taken. In practice, the time and displacement of measurements vary so slowly that we can ignore these dependences altogether. Therefore, we can equally well consider this function as one defined only on the unit radius sphere 𝑆², for which the following decomposition holds:Θ(𝑛)̂=_ℓ,𝑚𝑎_ℓ𝑚𝑌_ℓ𝑚(𝑛)̂.(3) Since the spherical harmonics 𝑌_ℓ𝑚(𝑛)̂ obey the symmetry 𝑌^∗ℓ𝑚(𝑛)̂=(−1)^𝑚𝑌_ℓ,−𝑚(𝑛)̂, the fact that the temperature field is a real function implies the identity 𝑎^∗ℓ𝑚=(−1)^𝑚𝑎_ℓ,−𝑚. This means that each temperature multipole ℓ is completely characterized by 2ℓ+1 real degrees of freedom.

2.1. From Inhomogeneities to Anisotropies: Linear Theory

The ultimate source of anisotropies in the Universe is the inhomogeneities in the baryon-photon fluid, as well as their associated spacetime metric fluctuations. If the photons were in perfect equilibrium with the baryons up to a sharply defined moment in time (the so-called instant recombination approximation), their distribution would have only one parameter (the equilibrium temperature at each point), so that photons flying off in any direction would have exactly the same energies. In that case, the photons we see today coming from a line-of-sight 𝑛̂ would reflect simply the density and gravitational potentials (the “sources”) at the position 𝑅𝑛̂, where 𝑅 is the radius to that (instantaneous) last scattering surface. Evidently, multiple scatterings at the epoch of recombination, combined with the fact that anisotropies themselves act as sources for more anisotropies, complicate this picture, and in general the relationship of the sources with the anisotropies must be calculated from either a set of Einstein-Boltzmann equations or, equivalently, from the line-of-sight integral equations coupled with the Einstein, continuity, and Euler equations [6].

Assuming for simplicity that recombination was instantaneous, at a time 𝜂_𝑅, the linear kernels of (1) reduce to 𝐾_𝑖(𝑥⃗^,𝜂^;𝑛)̂→𝛽_𝑖𝛿(𝜂^−𝜂_𝑅)𝛿(𝑥⃗^−𝑛𝑅̂), where 𝑅=𝜂₀−𝜂_𝑅 and 𝛽_𝑖 are constant coefficients. The photon distribution that we measure on Earth would therefore be given byΘ(𝑛)̂≈_𝑖0𝑥0200𝑑𝛽_𝑖𝑆^𝑖𝑥⃗^=𝑛𝑅̂,𝜂^=𝜂_𝑅.(4) We can also express this result in terms of the Fourier spectrum of the sources as follows:Θ(𝑛)̂≈_𝑖0𝑥0200𝑑𝛽_𝑖𝑑³𝑘(2𝜋)³𝑒^{𝑖⃗𝑘⋅}𝑛_𝑅̂𝑆^𝑖𝑘⃗,𝜂^𝑅.(5) Now we can use what is usually referred to as “Rayleigh's expansion” (though Watson, in his classic book on Bessel functions, attributes this to Bauer, J. f. Math. LVI, 1859) as follows:𝑒^{𝑖⃗𝑘⋅⃗𝑥}=4𝜋_ℓ𝑚𝑖^ℓ𝑗_ℓ(𝑘𝑥)𝑌^∗ℓ𝑚𝑘̂𝑌_ℓ𝑚(𝑥̂),(6) where 𝑗_ℓ(𝑧) are the spherical Bessel functions. Substituting (6) into (5) we obtain that𝑎_ℓ𝑚≈=𝑑𝑑²𝑛̂𝑌^∗ℓ𝑚(𝑛)̂Θ(𝑛)₃̂𝑘(2𝜋)³Θ𝑘⃗×4𝜋𝑖^ℓ𝑗_ℓ(𝑘𝑅)𝑌^∗ℓ𝑚𝑘̂,(7) where we have loosely collected the sources into the term Θ(𝑘⃗)≡∑_𝑖𝛽_𝑖𝑆^𝑖(𝑘⃗,𝜂_𝑅). This expression conveys well the simple relation between the Fourier modes and the spherical harmonic modes. Therefore, up to coefficients which are known given some background cosmology, the statistical properties of the harmonic coefficients 𝑎_ℓ𝑚 are inherited from those of the Fourier modes Θ(𝑘⃗) of the underlying matter and metric fields. Notice that the properties of the 𝑎_ℓ𝑚s under rotations, on the other hand, have nothing to do with the statistical properties of the fluctuations; they come directly from the spherical harmonic functions 𝑌_ℓ𝑚.

2.2. Statistics in Fourier Space

The characterization of the statistics of random variables is most commonly expressed in terms of the correlation functions. The two-point correlation function is the ensemble expectation value, 𝐶𝑘⃗,𝑘^⃗≡Θ𝑘⃗Θ𝑘^⃗.(8) In the absence of any symmetries, this would be a generic function of the arguments 𝑘⃗ and 𝑘^⃗, with only two constraints: first, because Θ(𝑥)⃗ is a real function, Θ^∗(𝑘⃗)=Θ(−𝑘⃗), hence, in our definition 𝐶^∗(𝑘⃗,𝑘^⃗)=𝐶(−𝑘⃗,−𝑘^⃗); second, due to the associative nature of the expectation value, 𝐶(𝑘⃗,𝑘^⃗)=𝐶(𝑘^⃗,𝑘⃗). It is obvious how to generalize this definition to 3, 4, or an arbitrary number of fields at different 𝑘⃗s (or “points”).

Let us first discuss the issue of gaussianity. If we say that the variables Θ(𝑘⃗) are Gaussian random numbers, then all the information that characterizes their distribution is contained in their two-point function 𝐶(𝑘⃗,𝑘^⃗). The probability distribution function (pdf) is then formally given by 𝑃Θ𝑘⃗,Θ𝑘^⃗∼exp⎡⎢⎢⎣−Θ𝑘⃗Θ𝑘^⃗2𝐶𝑘⃗,𝑘^⃗⎤⎥⎥⎦.(9) In this case, all higher-order correlation functions are either zero (for odd numbers of points) or they are simply connected to the two-point function by means of Wick's Theorem as follows: Θ𝑘¹Θ𝑘²⋯Θ𝑘^2𝑁_G=_{𝑁𝑖,𝑗}_𝛼=1𝐵^{𝛼𝑖,𝑗}Θ𝑘⃗^𝑖Θ𝑘^𝑗,(10) where the sum runs over all permutations of the pairs of wave vectors and 𝐵_𝑖,𝑗 are weights.

Second, let us consider the issue of homogeneity. A field is homogeneous if its expectation values (or averages) do not dependent on the spatial points where they are evaluated. In terms of the 𝑁-point functions in real space, we should have the following Θ𝑥₁Θ𝑥₂⋯Θ𝑥_𝑁_Homog.−−−−−→𝐶_𝑁𝑥₁−𝑥₂,…,𝑥_𝑁−1−𝑥_𝑁.(11) Writing this expression in terms of the Fourier modes, we get the followingΘ=𝑥₁𝑑Θ𝑥₂⋯Θ𝑥_𝑁₃𝑘₁𝑑³𝑘₁⋯𝑑³𝑘_𝑁(2𝜋)^3𝑁𝑒^{𝑘−𝑖1⋅𝑥1}𝑒^{𝑘−𝑖2⋅𝑥2}⋯𝑒^{𝑘−𝑖𝑁⋅𝑥𝑁}×Θ𝑘₁Θ𝑘²⋯Θ𝑘^𝑁.(12) Homogeneity demands that the expression in (12) is a function of the distances between spatial points only, not of the points themselves. Hence, the expectation value in Fourier space on the right-hand side of this expression must be proportional to 𝛿(𝑘₁+𝑘₂+⋯+𝑘_𝑁). In other words, the hypothesis of homogeneity constrains the 𝑁-point function in Fourier space to be of the following form: Θ𝑘¹Θ𝑘²⋯Θ𝑘^𝑁_H=×(2𝛿𝜋)³𝑘𝑁₁𝑘₁,𝑘₂,…,𝑘_𝑁+𝑘₂+⋯+𝑘_𝑁.(13) Notice that the “(𝑁−1)-spectrum” in Fourier space, 𝑁, can still be a function of the directions of the wavenumbers 𝑘⃗_𝑖 (it will be, in fact, a function of 𝑁−1 such vectors, due to the global momentum conservation expressed by the 𝛿-function.) Models which realize the general idea of (13) correspond to homogeneous but anisotropic universes [37–40].

There is a useful diagrammatic illustration for the 𝑁-point functions in Fourier space that enforce homogeneity. Notice that we could use the 𝛿-function in (13) to integrate out any one of the momenta 𝑘⃗_𝑖 in (12). Let us instead rewrite the 𝛿-functions in terms of triangles, so for the 4-point function we have 𝛿𝑘₁=+𝑘𝑑²₃+𝑘₃+𝑘₄𝑞𝛿𝑘₁+𝑘₂−𝑞⃗𝛿𝑘₃+𝑘₄+𝑞⃗,(14) whereas for the 5-point function we have 𝛿=𝑘₁+𝑑𝑘₂₃+𝑘₃+𝑘₄+𝑘₅×𝛿𝑞𝑞𝑑_⃗³𝑞^𝛿𝑘₁+𝑘₂−𝑞⃗𝛿𝑞⃗+𝑘₃−𝑞^⃗+𝑘₄+𝑘₅,(15) and so on, so that the 𝑁-point 𝛿-function is reduced to 𝑁−2 triangles with 𝑁−3 “internal momenta” (the idea is nicely illustrated in Figure 1.) Substituting the expression for the 𝑁-point 𝛿-function into (12) and integrating out all external momenta but the first (𝑘₁) and last (𝑘_𝑁), the result is as followsΘ=𝑥₁1Θ𝑥₂⋯Θ𝑥_𝑁(2𝜋)^3𝑁𝑑³𝑘₁𝑑³𝑞₁⋯𝑑³𝑞_𝑁−3𝑑³𝑘_𝑁×𝑒^{𝑖𝑘1⋅(𝑥1−𝑥2)}𝑒^{𝑖⃗𝑞1⋅(𝑥2−𝑥3)}⋯𝑒^{𝑖𝑞𝑁−3⋅(𝑥𝑁−2−𝑥𝑁−1)}×𝑒^{𝑖𝑘𝑁⋅(𝑥𝑁−1−𝑥𝑁)}Θ𝑘¹Θ𝑞¹⃗−𝑘₁⋯Θ𝑘^𝑁.(16) This expression shows explicitly that the real-space 𝑁-point function above does not depend on any particular spatial point, only on the intervals between points.

Figure 1 

Diagrammatic representation of the 2, 3, 4, and 5-point correlation functions in Fourier space. The dashed lines represent internal momenta.

Finally, what are the constraints imposed on the 𝑁-point functions that come from isotropy alone? Clearly, no dependence on the directions defined by the points, 𝑥⃗_𝑖−𝑥_𝑗, can arise in the final expression for the 𝑁-point functions in real space, so from (12) we see that the 𝑁-point function in Fourier space should depend only on the moduli of the wavenumbers—up to some momentum-conservation 𝛿-functions, which naturally carry vector degrees of freedom.

In this paper, we will mostly be concerned with tests of isotropy given homogeneity (but not necessarily Gaussianity), so in our case we will usually assume that the 𝑁-point function in Fourier space assumes the form given in (13).

2.3. Statistics in Harmonic Space

In the previous Section, we characterized the statistics of our field in Fourier space, which in most cases is most easily related to fundamental models such as inflation. Now we will change to harmonic representation, because that is what is most directly related to the observations of the CMB, Θ(𝑛)̂, which are taken over the unit sphere 𝑆². We will discuss mostly the two-point function here, and we defer a fuller discussion of 𝑁-point functions in harmonic space to Section 3.

From (7), we can start by taking the two-point function in harmonic space, and computing it in terms of the two-point function in Fourier space as follows 𝑎^ℓ𝑚𝑎^∗ℓ_′𝑚′=𝑑³𝑘𝑑³𝑘^(2𝜋)⁶×𝑌_{∗ℓ′𝑚′}(4𝜋)²𝑖^ℓ(−𝑖)^ℓ′𝑗_ℓ(𝑘𝑅)𝑗_ℓ^′𝑘^𝑅𝑌_ℓ𝑚𝑘̂𝑘̂^Θ𝑘⃗Θ^∗𝑘⃗^.(17) Under the hypothesis of homogeneity, this expression simplifies considerably, leading to 𝑎^ℓ𝑚𝑎^∗ℓ_′𝑚′_H=𝑑³𝑘2𝜋𝑖^ℓ(−𝑖)^ℓ′𝑗_ℓ(𝑘𝑅)𝑗_ℓ^′(𝑘𝑅)𝑌_ℓ𝑚𝑘̂𝑌^∗ℓ_′𝑚′𝑘̂×𝑁₂𝑘⃗.(18) If, in addition to homogeneity, we also assume isotropy, then 𝑁₂→𝑃(𝑘), and the integration over angles factors out, leading to the orthogonality condition for spherical harmonics as follows 𝑑²𝑘̂𝑌_ℓ𝑚𝑘̂𝑌^∗ℓ_′𝑚′𝑘̂=𝛿_ℓℓ^′𝛿_𝑚𝑚′,(19) and as a result the covariance of the 𝑎_ℓ𝑚s becomes diagonal as follows𝑎^ℓ𝑚𝑎^∗ℓ_′𝑚′_H,I==𝛿4𝜋_ℓℓ^′𝛿𝛿_{𝑚𝑚′ℓℓ′}𝑑𝑘𝑘𝑗^2ℓ(𝑘𝑅)𝜋2𝑘³𝑃(𝑘)𝛿_𝑚𝑚′𝑑log𝑘𝑗^2ℓ≡𝐶_ℓ(𝑘𝑅)Δ^2𝑇(𝑘)𝛿_ℓℓ^′𝛿_𝑚𝑚′,(20) where we have defined the usual temperature power spectrum Δ_𝑇(𝑘)=𝑘³𝑃(𝑘)/2𝜋² in the middle line, and the angular power spectrum 𝐶_ℓ in the last line of (20). As a pedagogical note, let us recall that the power spectrum basically expresses how much power the two-point correlation function has per unit log𝑘 as follows: Θ𝑥⃗Θ𝑥⃗^_H,𝐼=𝑑log𝑘sin𝑘||𝑥⃗−𝑥⃗^||𝑘||𝑥⃗−𝑥⃗^||Δ^2𝑇(𝑘).(21)

In an analogous manner to what was done above, we can also construct the angular two-point correlation function in harmonic space as follows:Θ(𝑛)̂Θ𝑛^̂=_ℓ𝑚_{ℓ′𝑚′}𝑎^ℓ𝑚𝑎^∗ℓ_′𝑚′𝑌^ℓ𝑚(𝑛)̂𝑌^∗ℓ𝑚𝑛^̂.(22) The hypothesis of homogeneity by itself does not lead to significant simplifications, but isotropy leads to a very intuitive expression for the angular two-point function as follows:Θ(𝑛)̂Θ𝑛^̂_H,I==_ℓ𝑚ℓ0𝑥0200𝑑_{ℓ′𝑚′}𝐶_ℓ𝛿_ℓℓ^′𝛿_𝑚𝑚′𝑌_ℓ𝑚(𝑛)̂𝑌^∗ℓ𝑚𝑛^̂𝐶_ℓ2ℓ+14𝜋𝑃^ℓ𝑛̂⋅𝑛^̂.(23) Clearly, not only is this expression the analogous in 𝑆² of (21), but in fact the Fourier power spectrum Δ^2𝑇(𝑘) and the angular power spectrum 𝐶_ℓ are defined in terms of each other as indicated in (23) as follows: 𝐶_ℓ=4𝜋𝑑log𝑘𝑗^2ℓ(𝑘𝑅)Δ^2𝑇(𝑘).(24) Now, using the facts that the spherical Bessel function of order ℓ peaks when its argument is approximately given by ℓ, and that ∫𝑑log𝑧𝑗^2ℓ(𝑧)=1/(2ℓ(ℓ+1)), we obtain the following (this is one type of what has become known in the literature as Limber's approximations )𝐶_ℓ≈2𝜋ℓ(ℓ+1)Δ^2𝑇𝑘=𝑅ℓ.(25) Incidentally, from this expression it is clear why it is customary to define 𝒞_ℓ≡ℓ(ℓ+1)2𝜋𝐶^ℓ≈Δ^2𝑇𝑘=ℓ𝑅.(26)

Using (12), we can easily generalize the results of this subsection to 𝑁-point functions in 𝑆² and in harmonic space, however, the assumption of isotropy alone does very little to simplify our life. The hypothesis of homogeneity, on the other hand, greatly simplifies the angular 𝑁-point functions, and most of the work in statistical anisotropy of the CMB that goes beyond the two-point function assumes that homogeneity holds. Notice that the issue of gaussianity is, as always, confined to the question of whether or not the two-point function holds all information about the distribution of the relevant variables and is therefore completely separated from questions about homogeneity and/or isotropy.

Also notice that the separable nature of the definition (22) implies here as well, like in Fourier space, a reciprocity relation for the correlation function 𝐶𝑛₁̂,𝑛₂̂=𝐶𝑛₂̂,𝑛₁̂.(27) This symmetry must hold regardless of underlying models and is important in order to analyze the symmetries of the correlation function, as we will see later.

Before we move on, it is perhaps important to mention that the decomposition (22) is not unique. In fact, instead of the angular momenta of the parts, (ℓ₁,𝑚₁;ℓ₂,𝑚₂), we could equally well have used the basis of total angular momentum (𝐿,𝑀;ℓ₁,ℓ₂) and decomposed that expression as 𝐶𝑛₁̂,𝑛₂̂=_𝐿,𝑀ℓ_1,ℓ2𝒜^ℓ𝐿𝑀_1ℓ2𝒴^ℓ𝐿𝑀_1ℓ2𝑛₁̂,𝑛₂̂,(28) where 𝒴^ℓ𝐿𝑀_1ℓ2 are known as the bipolar spherical harmonics, defined by [41]𝒴^ℓ𝐿𝑀_1ℓ2𝑛₁̂,𝑛₂̂=𝑌_ℓ1𝑛₁̂⊗𝑌_ℓ2𝑛₂̂_𝐿𝑀,(29) where 𝐿 and 𝑀=𝑚₁+𝑚₂ are the eigenvalues of the total and azimuthal angular momentum operators, respectively. This decomposition is completely equivalent to (22), and we can exchange from one decomposition to another by using the relation𝒜^ℓ𝐿𝑀_1ℓ2=_𝑚1𝑚2𝑎_ℓ1𝑚1𝑎_ℓ2𝑚2(−1)^{𝑀+ℓ1−ℓ2}√2𝐿+1ℓ¹𝑚₁ℓ₂𝐿𝑚₂−𝑀,(30) where the 3×2 matrices above are the well-known 3-j coefficients. At this point, it is only a matter of mathematical convenience whether we choose to decompose the correlation function as in (22) or as in (28). Although the bipolar harmonics behave similarly to the usual spherical harmonics in many aspects, the modulations of the correlation function as described in this basis have a peculiar interpretation. We will not go further into detail about this decomposition here, as it is discussed at length in another review article in this volume.

2.4. Estimators and Cosmic Variance

Returning to the covariance matrix (20), we see that, if we assume gaussianity of the 𝑎_ℓ𝑚s, then the angular power spectrum suffices to describe statistically how much the temperature fluctuates in any given angular scale; all we have to do is to calculate the average (20). This can be a problem, though, since we have only one Universe to measure, and therefore only one set of 𝑎_ℓ𝑚s. In other words, the average in (20) is poorly determined.

At this point, the hypothesis that our Universe is spatially homogeneous and isotropic at cosmological scales comes not only as simplifying assumption about the spacetime symmetries, but also as a remedy to this unavoidable smallness of the working cosmologist's sample space. If isotropy holds, different cosmological scales are statistically independent, which means that we can take advantage of the ergodic hypothesis and trade averaging over an ensemble for averaging over space. In other words, for a given ℓ we can consider each of the 2ℓ+1 real numbers in 𝑎_ℓ𝑚 as statistically independent Gaussian random variables, and define a statistical estimator for their variances as the average𝐶_ℓ≡12ℓ+1^ℓ_𝑚=−ℓ||𝑎_ℓ𝑚||².(31) The smaller the angular scales (ℓ bigger), the larger the number of independent patches that the CMB sky can be divided into. Therefore, in this limit we should havelim_ℓ→∞𝐶_ℓ=𝐶_ℓ.(32) On the other hand, for large angular scales (small ℓ's), the number of independent patches of our Universe becomes smaller, and (31) becomes a weak estimation of the 𝐶_ℓs. This means that any statistical analysis of the Universe on large scales will be plagued by this intrinsic cosmic sample variance. Notice that this is an unavoidable limit as long as we have only one observable Universe.

Finally, it is important to keep in mind the clear distinction between the angular power spectrum 𝐶_ℓ and its estimator (31). The former is a theoretical variable which can be calculated from first principles, as we have shown in Section 2.1. The latter, being a function of the data, is itself a random variable. In fact, if the 𝑎_ℓ𝑚s are Gaussian, then we can rewrite expression (31) as(2ℓ+1)𝐶_ℓ𝐶_ℓ=𝑋_ℓ,𝑋_ℓ=^ℓ_𝑚_=−ℓ||𝑎_ℓ𝑚||²𝐶_ℓ,(33) where 𝑋_ℓ is a chi-square random variable with 2ℓ+1 degrees of freedom. According to the central limit theorem, when ℓ→∞, 𝑋_ℓ approaches a standard normal variable (A standard normal variable is a Gaussian variable 𝑋 with zero mean and unit variance. Any other Gaussian variable 𝑌 with mean 𝜇 and variance 𝜎 can be obtained from 𝑋 through 𝑌=𝜎𝑋+𝜇.) which implies that 𝐶_ℓ will itself follow a Gaussian distribution. Its mean can be easily calculated using (20) and (31) and is of course given by𝐶^ℓ=𝐶^ℓ,(34) which shows that the 𝐶_ℓs are unbiased estimators of the 𝐶_ℓs. It is also straightforward to calculate its variance (valid for any ℓ)𝐶^ℓ−𝐶_ℓ𝐶^ℓ′−𝐶_ℓ^′=2ℓ2+1𝐶^2ℓ𝛿_ℓℓ^′.(35) Because this estimator does not couple different cosmological scales, it has the minimum cosmic variance we can expect from an estimator due to the finiteness of our sample—so it is optimal in that sense. 𝐶_ℓ is therefore the best estimator we can build to measure the statistical properties of the multipolar coefficients 𝑎_ℓ𝑚 when both statistical isotropy and gaussianity hold.

In later Sections, we will explore angular or harmonic 𝑁-point functions for which the assumption of isotropy does not hold. However, it is important to remember at all times that we have only one map, which means one set of 𝑎_ℓ𝑚s. The estimator for the angular power spectrum, 𝐶_ℓ, takes into account all the 𝑎_ℓ𝑚s by dividing them into the different ℓs and summing over all 𝑚∈(−ℓ,ℓ). Clearly, it will inherit a sample variance for small ℓs, when the 𝑎_ℓ𝑚s can only be divided into a few “independent parts”. As we try to estimate higher-order objects such as the 𝑁-point functions, we will have to subdivide the 𝑎_ℓ𝑚's into smaller and smaller subsamples, which are not necessarily independent (in the statistical sense) of each other. So, the price to pay for aiming at higher-order statistics is a worsening of the cosmic sample variance.

2.5. Correlation and Statistical Independence

The covariance given in (20) has two distinct, important properties. First, note that its diagonal entries, the 𝐶_ℓ's, are 𝑚-independent coefficients; this is crucial for having statistical isotropy, as we will show latter. Second, statistical isotropy at the Gaussian level implies that different cosmological “scales” (understood here as meaning the modes with total angular momentum ℓ and azimuthal momentum 𝑚) should be statistically independent of each other—and this is represented by the Kronecker deltas in (20).

In fact, statistical independence of cosmological scales is a particular property of Gaussian and statistically isotropic random fields and is not guaranteed to hold when gaussianity is relaxed. We will see in the next Section that the rotationally invariant 3-point correlation function (and in general any 𝑁>2 correlation function) couples to the three scales involved. In particular, if it happens that the Gaussian contribution of the temperature field is given by (20), but at least one of its non-Gaussian moments are nonzero, then the fact that a particular correlation is zero, like for example ⟨𝑎_2𝑚1𝑎^∗3𝑚₂⟩, does not imply that the scales ℓ=2 and ℓ=3 are (statistically) independent. This is just a restatement of the fact that, while statistical independence implies null correlation, the opposite is not necessarily true. This can be illustrated by the following example: consider a random variable 𝛼 distributed as𝑃(𝛼)=10,,𝛼oth∈er[w0i,s1e].,(36) Let us now define two other variables 𝑥=cos(2𝜋𝛼) and 𝑦=sin(2𝜋𝛼). From these definitions, it follows that 𝑥 and 𝑦 are statistically dependent variables, since knowledge of the mean/variance of 𝑥 automatically gives the mean/variance of 𝑦. However, these variables are clearly uncorrelated⟨𝑥𝑦⟩=12𝜋^02𝜋cos𝜂sin𝜂𝑑𝜂=0.(37) Although correlations are among cosmologist's most popular tools when analyzing CMB properties, statistical independence may turn out to be an important property as well, specially at large angular scales, where cosmic variance is more of a critical issue.

3. Beyond the Standard Statistical Model

Until now we have been analyzing the properties of Gaussian and statistically isotropic random temperature fluctuations. This gives us a fairly good statistical description of the Universe in its linear regime, as confirmed by the astonishing success of the ΛCDM model. This picture is incomplete though, and we have good reasons to search for deviations of either gaussianity and/or statistical isotropy. For example, the observed clustering of matter in galactic environments certainly goes beyond the linear regime where the central limit theorem can be applied, therefore leading to large deviations of gaussianity in the matter power spectrum statistics. Besides, deviations of the cosmological principle may leave an imprint in the statistical moments of cosmological observables, which can be tested by searching for spatial inhomogeneities [42] or directionalities [43].

But how do we plan to go beyond the standard model, given that there is only one Gaussian and statistically isotropic description of the Universe, but infinite possibilities otherwise? This is in fact an ambitious endeavor, which may strongly depend on observational and theoretical hints on the type of signatures we are looking for. In the absence of extra input, it is important to classify these signatures in a general scheme, differentiating those which are non-Gaussian from those which are anisotropic. Furthermore, given that the signatures of non-gaussianity may in principle be quite different from that of statistical anisotropy, such a classification is crucial for data analysis, which requires sophisticated tools capable of separating these two issues.(Although gaussianity and homogeneity/isotropy are mathematically distinct properties, it is possible for a Gaussian but inhomogeneous/anisotropic model to look like an isotropic and homogeneous non-gaussian model. See, e.g., [44].)

We therefore start Section 3.1 by analyzing deviations of gaussianity when statistical isotropy holds. In Section 3.2, we keep the hypothesis of gaussianity and analyze the consequences of breaking statistical rotational invariance.

3.1. Non-Gaussian and SI Models

3.1.1. Rotational Invariance of 𝑁-Point Correlation Functions

We turn now to the question of non-Gaussian but statistically isotropic probabilities distributions. We will keep working with the 𝑁-point correlation function defined in harmonic space,𝑎_ℓ1𝑚1𝑎_ℓ2𝑚2…𝑎_{ℓ𝑁𝑚𝑁},(38) since knowledge of these functions enables one to fully reconstruct the CMB temperature probability distribution. Specifically, we would like to know the form of any 𝑁-point correlation function which is invariant under arbitrary 3-dimensional spatial rotations. When rotated to a new (primed) coordinate system, the 𝑁-point correlation function transforms as𝑎=^ℓ₁^𝑚′1𝑎_ℓ2𝑚^′2_all𝑚…𝑎_ℓ𝑁𝑚^′𝑁𝑎_ℓ1𝑚1𝑎_ℓ2𝑚2…𝑎_{ℓ𝑁𝑚𝑁}𝐷^ℓ1_{𝑚′1𝑚1}𝐷^ℓ2_{𝑚′2𝑚2}⋯𝐷^ℓ𝑛_{𝑚′𝑁𝑚𝑁},(39) where the 𝐷^ℓ𝑚_𝑖𝑚′i(𝛼,𝛽,𝛾)s are the coefficients of the Wigner rotation-matrix, which depend on the three Euler-angles 𝛼, 𝛽, and 𝛾 characterizing the rotation. Notice that in this notation the primed (rotated) system is indicated by the primed 𝑚s. For the 2-point correlation function, we have already seen that the well-known expression𝑎_ℓ1𝑚1𝑎_ℓ2𝑚2=(−1)^𝑚2𝐶_ℓ1𝛿_ℓ1ℓ2𝛿_{𝑚1,−𝑚2}(40) does the job 𝑎^ℓ₁^𝑚′1𝑎_ℓ2𝑚^′2=𝐶^ℓ₁_𝑚1(−1)^𝑚1𝐷^ℓ1_{𝑚′1𝑚1}𝐷^ℓ1_{𝑚′2−𝑚1}𝛿^ℓ1ℓ2=(−1)^𝑚′2𝐶_ℓ1𝛿_ℓ1ℓ2𝛿_𝑚^′1_,−𝑚^′2.(41) Note the importance of the angular spectrum, 𝐶_ℓ, being an 𝑚-independent function.

What about the 3-point function? In this case, the invariant combination is found to be𝑎_ℓ1𝑚1𝑎_ℓ2𝑚2𝑎_ℓ3𝑚3=𝐵_ℓ1ℓ2ℓ3ℓ¹ℓ₂ℓ₃𝑚₁𝑚₂𝑚₃(42) which can be verified by straightforward calculations. Again, the nontrivial physical content of this statistical moment is contained in an arbitrary but otherwise 𝑚-independent function: the bispectrum 𝐵_ℓ1ℓ2ℓ3 [45–47]. As we anticipated in Section 2, rotational invariance of the 3-point correlation is not enough to guarantee statistical independence of the three cosmological scales involved in the bispectrum, although in principle a particular model could be formulated to ensure that 𝐵_ℓ1ℓ2ℓ3∝𝛿_ℓ1ℓ2𝛿_ℓ2ℓ3, at least for some subset of a general geometric configuration of the 3-point correlation function.

These general properties hold for all the 𝑁-point correlation function. For the 4-point correlation function, for example, Hu [48] has found the following rotationally invariant combination𝑎=_{ℓ1𝑚1𝐿𝑀}…𝑎_ℓ4𝑚4𝑄^ℓ1ℓ2_ℓ3ℓ4(𝐿)(−1)^𝑀ℓ¹𝑚₁ℓ₂𝐿𝑚₂−𝑀ℓ³𝑚₃ℓ₄𝐿𝑚₄𝑀,(43) where the 𝑄^ℓ1ℓ2_ℓ3ℓ4(𝐿) function is known as the trispectrum, and 𝐿 is an internal angular momentum needed to ensure parity invariance. In a likewise manner, it can be verified that the following expression 𝑎=_{ℓ1𝑚1𝐿𝑀}⋯𝑎_ℓ5𝑚50𝑥0200𝑑_{𝐿′𝑀′}𝑃^ℓ1ℓ2_ℓ3ℓ4ℓ5𝐿,𝐿^(−1)^𝑀+𝑀′ℓ¹𝑚₁ℓ₂𝐿×ℓ𝑚₂−𝑀₃ℓ₄𝐿^𝑚₃𝑚₄−𝑀^ℓ⁵𝐿𝐿^𝑚₅𝑀𝑀^(44) gives the rotationally invariant quadrispectrum 𝑃^ℓ1ℓ2_ℓ3ℓ4ℓ5(𝐿,𝐿^).

The examples above should be enough to show how the general structure of these functions emerges under SI; apart from an 𝑚-independent function, every pair of momenta ℓ_𝑖 in these functions are connected by a triangle, which in turn connects itself to other triangles through internal momenta when more than 3 scales are present. In Figure 2, we show some diagrams representing the functions above.

Figure 2 

Diagrammatic representation of the 2, 3, 4, and 5-point correlation functions in harmonic space. Here ℓ⃗ actually represents the pair (ℓ,𝑚).

Figure 3 

Geometrical representation of the 2pcf in terms of two unit vectors.

Although we have always shown 𝑁-point functions which are rotationally invariant, the procedure used for obtaining them was rather intuitive, and therefore does not offer a recipe for constructing general invariant correlation functions. Furthermore, it does not guarantee that this procedure can be extended for arbitrary 𝑁s. Here we will present a recipe for doing that, which also guarantees the uniqueness of the solution.

The general recipe for obtaining the rotationally invariant 𝑁-point function is as follows: from the expression (39) above, we start by contracting every pairs of Wigner functions, where by “contracting” we mean using the identity 𝐷^ℓ1_{𝑚′1𝑚1}(𝜔)𝐷^ℓ2_{𝑚′2𝑚2}(𝜔)=_{𝐿,𝑀,𝑀′}0𝑥0𝑒𝑓11ℓ¹𝑚_1ℓ₂𝐿𝑚^2−𝑀ℓ¹𝑚₁ℓ₂𝐿×(2𝐿+1)(−1)_𝑀+𝑀′𝑚₂−𝑀^𝐷^𝐿𝑀_′𝑀(𝜔),(45) and where 𝜔={𝛼,𝛽,𝛾} is a shortcut notation for the three Euler angles. Once this contraction is done, there will remain 𝑁/2𝐷-functions, which can again be contracted in pairs. This procedure should be repeated until there is only one Wigner function left, in which case we will have an expression of the following form:𝑎_ℓ1𝑚1⋯𝑎_{ℓ𝑁𝑚𝑁}=×_all𝑚′g𝑎eom^ℓ₁^𝑚′1e⋯tri𝑎calf^ℓ_𝑁^𝑚′𝑁actors×𝐷_{𝐿𝑀𝑀′}(𝜔).(46) Now, we see that the only way for this combination to be rotationally invariant is when the remaining 𝐷^𝐿𝑀𝑀_′ function above does not depend on 𝜔, that is, 𝐷^𝐿𝑀𝑀_′(𝜔)=𝛿_𝐿0𝛿_𝑀0𝛿_𝑀^′₀. Once this identity is applied to the geometrical factors, we are done, and the remaining terms inside the primed 𝑚-summation will give the rotationally invariant (𝑁−1)-spectrum.

As an illustration of this algorithm, let us construct the rotationally invariant spectrum and bispectrum. For the 2-point function, there is only one contraction to be done, and after we simplify the last Wigner function, we arrive at𝑎_ℓ1𝑚1𝑎_ℓ2𝑚2=⎡⎢⎢⎢⎢⎣_𝑚′1|||𝑎^ℓ1𝑚′1|||²2ℓ₁+1⎤⎥⎥⎥⎥⎦(−1)^𝑚2𝛿_ℓ1ℓ2𝛿_{𝑚1,−𝑚2},(47) where, of course𝐶_ℓ≡12ℓ+1_𝑚||𝑎^ℓ𝑚||²(48) is the well-known definition of the temperature angular spectrum. For the 3-point function, there are two contractions, and the simplification of the last Wigner function gives=𝑎⎡⎢⎢⎣_ℓ1𝑚1𝑎_{ℓ2𝑚2𝑚′1,𝑚′2,𝑚′3}𝑎_ℓ3𝑚3𝑎^ℓ₁^𝑚′1𝑎_ℓ2𝑚^′2𝑎_ℓ3𝑚^′3ℓ¹ℓ₂ℓ₃𝑚^1𝑚^2𝑚^3⎤⎥⎥⎦ℓ¹ℓ₂ℓ₃𝑚₁𝑚₂𝑚₃.(49) From this expression and the ortoghonality of the 3-js symbols (see the Appendix), we can immediately identify the definition of the bispectrum as follows:𝐵_ℓ1ℓ2ℓ3≡_{𝑚1,𝑚2,𝑚3}𝑎_ℓ1𝑚1𝑎_ℓ2𝑚2𝑎_ℓ3𝑚3ℓ¹ℓ₂ℓ₃𝑚₁𝑚₂𝑚₃.(50)

It should be mentioned that this recipe not only enables us to establish the rotational invariance of any 𝑁-point correlation function, but it also furnishes a straightforward definition of unbiased estimators for the 𝑁-point functions. All we have to do is to drop the ensemble average of the primed 𝑎_ℓ𝑚's. So, for example, for the 2- and 3-point functions above, the unbiased estimators are given, respectively, by 𝐶_ℓ=12ℓ+1_𝑚𝐵_ℓ1ℓ2ℓ3𝑎_ℓ𝑚𝑎^∗ℓ𝑚,=_{𝑚1,𝑚2,𝑚3}𝑎_ℓ1𝑚1𝑎_ℓ2𝑚2𝑎_ℓ3𝑚3ℓ¹ℓ₂ℓ₃𝑚₁𝑚₂𝑚₃.(51)

Notice that isotropy plays the same role, in 𝑆², that homogeneity plays in ℝ³. What enforces homogeneity in ℝ³ is the Fourier-space 𝛿-functions, as in the discussion around Figure 1. However, in 𝑆², the equivalents of the Fourier modes are the harmonic modes, for which there is only a discrete notion of orthogonality—and no Dirac 𝛿-function. What we found above is that the Wigner 3-j symbols play the same role as the Fourier space 𝛿-functions; they are the enforcers of isotropy (rotational invariance) for the 𝑁-point angular correlation function. Hence, the diagrammatic representations of the constituents of the 𝑁-point functions in Fourier (Figure 1) and in harmonic space (Figure 2) really do convey the same physical idea—one in ℝ³, the other in 𝑆².

3.2. Gaussian and Statistically Anisotropic Models

In the last Section, we have developed an algorithm which enables one to establish the rotational invariance of any 𝑁-point correlation function. As we have shown, this is also an algorithm for building unbiased estimators of non-Gaussian correlations. In this Section, we will change the perspective and analyze the case of Gaussian but statistically anisotropic models of the Universe.

There are many ways in which statistical anisotropy may be manifested in CMB. From a fundamental perspective, a short phase of inflation which produces just enough e-folds to solve the standard Big Bang problems may leave imprints on the largest scales of the Universe, provided that the spacetime is sufficiently anisotropic at the onset of inflation [39]. Another source of anisotropy may result from our inability to efficiently clean foreground contaminations from temperature maps. Usually, the cleaning procedure involves the application of a mask function in order to eliminate contaminations of the galactic plane from raw data. As a consequence, this procedure may either induce, as well as hide some anomalies in CMB maps [28].

It is important to mention that these two examples can be perfectly treated as Gaussian: in the first case, the anisotropy of the spacetime can be established in the linear regime of perturbation theory and therefore will not destroy gaussianity of the quantum modes, provided that they are initially Gaussian. In the second case, the mask acts linearly over the temperature maps, therefore preserving its probability distribution [49].

3.2.1. Primordial Anisotropy

Recently, there have been many attempts to test the isotropy of the primordial Universe through the signatures of an anisotropic inflationary phase [38–40, 50, 51]. A generic prediction of such models is the linear coupling of the scalar, vector, and tensor modes through the spatial shear, which is in turn induced by anisotropy of the spacetime [38]. Whenever that happens, the matter power spectrum, defined in a similar way as in (13), will acquire a directionality dependence due to this type of see-saw mechanism. This dependence can be accommodated in a harmonic expansion of the form𝑃𝑘⃗=_ℓ,𝑚𝑟_ℓ𝑚(𝑘)𝑌_ℓ𝑚𝑘̂,(52) where the reality of 𝑃(𝑘⃗) requires that 𝑟_ℓ𝑚(𝑘)=(−1)^𝑚𝑟^{∗ℓ,−𝑚}(𝑘). Given that temperature perturbations Θ(𝑥)⃗ are real, their Fourier components must satisfy the relation Θ(𝑘⃗)=Θ^∗(−𝑘⃗). This property taken together with the definition (13) implies that𝑃𝑘⃗=𝑃−𝑘⃗,(53) which in turn restricts the ℓs in (52) to even values. Also, note that by relaxing the assumption of spatial isotropy, we are only breaking a continuous spacetime symmetry, but discrete symmetries such as parity should still be present in this class of models. Indeed, by imposing invariance of the spectrum under the transformation 𝑧→−𝑧, we find that (−1)^ℓ−𝑚=1. Similarly, invariance under the transformations 𝑥→−𝑥 and 𝑦→−𝑦 implies the conditions 𝑟_ℓ𝑚=(−1)^𝑚𝑟_ℓ,−𝑚 and 𝑟_ℓ𝑚=𝑟_ℓ,−𝑚, respectively. Gathering all these constraints with the parity of ℓ, we conclude that𝑟_ℓ𝑚∈ℝ,ℓ,𝑚∈2ℕ.(54) That is, from the initial 2ℓ+1 degrees of freedom, only ℓ/2+1 of them contribute to the anisotropic spectrum [39].

3.2.2. Signatures of Statistical Anisotropies

The selection rules (54) are the most generic predictions we can expect from models with global break of anisotropy. We will now work out the consequences of these rules to the temperature power spectrum and check whether they can say something about the CMB large-scale anomalies.

From the expressions (17) and (52), we can immediately calculate the most general anisotropic covariance matrix [37] as follows𝑎^ℓ₁^𝑚₁𝑎^∗ℓ_2𝑚2=_ℓ3,𝑚3𝒢^ℓ1ℓ2ℓ3_{𝑚1𝑚2𝑚3}𝐻^ℓ3𝑚3_ℓ1ℓ2,(55) where𝒢^ℓ1ℓ2ℓ3_{𝑚1𝑚2𝑚3}=(−1)^𝑚12ℓ₁+12ℓ₂+12ℓ₃+1×ℓ4𝜋₁ℓ₂ℓ₃000ℓ¹ℓ₂ℓ₃−𝑚₁𝑚₂𝑚₃(56) are the Gaunt coefficients resulting from the integral of three spherical harmonics (see the Appendix). These coefficients are zero unless the following conditions are met:ℓ𝑚₁+ℓ₂+ℓ₃∈2ℕ,₁||ℓ+𝑚₂+𝑚₃=0,_𝑖−ℓ_𝑗||≤ℓ_𝑘≤ℓ_𝑖+ℓ_𝑗,∀𝑖,𝑗,𝑘∈{1,2,3}.(57) The remaining coefficients in (55) are given by𝐻^ℓ3𝑚3_ℓ1ℓ2=4𝜋𝑖^ℓ1−ℓ2^∞0𝑑log𝑘𝑟_ℓ3𝑚3(𝑘)𝑗_ℓ1(𝑘𝑅)𝑗_ℓ2(𝑘𝑅)(58) and correspond to the anisotropic generalization of temperature power spectrum (20).

The selection rules (57), taken together with (54), lead to important signatures in the CMB. In particular, since ℓ₃ is even, the quantity ℓ₁±ℓ₂ must also be even, that is, multipoles with different parity do not couple to each other in this class of models𝑎^ℓ₁^𝑚₁𝑎^∗ℓ_2𝑚2=0,ℓ¹+ℓ₂=odd.(59) This result is in fact expected on theoretical grounds, because by breaking a continuous symmetry (isotropy) we cannot expect to generate a discrete signature (parity) in CMB. However, notice that the absence of correlations between, say, the quadrupole and the octupole, does not imply that there will be no alignment between them. One example of this would be a covariance matrix of the form 𝐶_ℓ𝑚𝛿_ℓℓ^′𝛿_𝑚𝑚′. If the 𝐶_ℓ,0 happens to be zero, for example, then all multipoles will present a preferred direction (in this case, the 𝑧-axis).

3.2.3. Isotropic Signature of Statistical Anisotropy

We have just shown that a generic consequence of an early anisotropic phase of the Universe is the generation of even-parity signatures in CMB maps. Interestingly, these signatures may be present even in the isotropic angular spectrum, since the 𝐶_ℓs acquire some additional modulations in the presence of statistical anisotropies. In principle, these modulations could be constrained by measuring an effective angular spectrum of the form𝑎_ℓ𝑚𝑎^∗ℓ𝑚=𝐶_ℓ+𝜖_ℓ′>0𝒢^ℓℓℓ′_{𝑚,−𝑚,0}𝐻^{ℓ′0ℓℓ},(60) where we have introduced a small 𝜖 parameter to quantify the amount of primordial anisotropy.

In order to constrain these modulations, we have to build a statistical estimator for ⟨𝑎_ℓ𝑚𝑎^∗ℓ𝑚⟩. Since we are looking for the diagonal entries of the matrix (55), a first guess would be𝐶^eff_ℓ=12ℓ+1^ℓ_𝑚=−ℓ𝑎_ℓ𝑚𝑎^∗ℓ𝑚.(61) To check whether this is an unbiased estimator of the effective angular spectrum, we apply it to (55) and take its average𝐶^eff_ℓ=𝐶_ℓ+𝜖_{ℓ′>0,𝑚}𝒢^ℓℓℓ′_{𝑚,−𝑚,0}𝐻^{ℓ′0ℓℓ}.(62) Using the definition (56) of the Gaunt coefficients, the 𝑚-summation in the expression above becomes [41]_𝑚(−1)^ℓ−𝑚𝑚ℓ−ℓ𝑚ℓ0^=√2ℓ+1𝛿_ℓ^′_,0=0,(63) where the last equality follows because ℓ′>0. Consequently, we conclude that𝐶^effℓ=𝐶^ℓ≠𝑎_ℓ𝑚𝑎^∗ℓ𝑚.(64) At first sight, this result may seem innocuous, showing only that this is not an appropriate estimator for (60). Note, however, that (61) is in fact the estimator of the angular spectrum usually applied to CMB data under the assumption of statistical isotropy. In other words, by means of the usual procedure we may be neglecting important information about statistical anisotropy. Moreover, the cosmic variance induced by the application of this estimator on anisotropic CMB maps is small, because, as it can easily be checked,𝐶^effℓ−𝐶_ℓ𝐶^effℓ′−𝐶_ℓ^′=2ℓ2+1𝐶^2ℓ𝛿_ℓℓ^′+𝒪𝜖².(65)

This result shows that the construction of statistical estimators strongly depends on our prejudices about what non-gaussianity and statistical anisotropy should look like. Consequently, an estimator built to measure one particular property of the CMB may equally well hide other important signatures. One possible solution to this problem is to let the construction of our estimators be based on what the observations seem to tell us, as we will do in the next section.

Part II: Cosmological Tests

So much for mathematical formalism. We will now turn to the question of how the hypotheses of gaussianity and statistical isotropy of the Universe can be tested. Though we are primarily interested in applying these tests to CMB temperature maps, most of the tools we will be dealing with can be applied to polarization (𝐸- and 𝐵-mode maps) and to catalogs of large-scale structures as well.

Testing the gaussianity and SI of our Universe is a difficult task. Specially because, as we have seen, there is only one Gaussian and SI Universe, but infinitely many universes which are neither Gaussian nor isotropic. So what type of non-gaussianity and statistical anisotropy should we test for? In order to attack this problem, we can follow two different routes. In the bottom-up approach, models for the cosmic evolution are formulated in such a way as to account for some specific deviations from gaussianity and SI. These physical principles range from nontrivial cosmic topologies [52–54], primordial magnetic fields [55–59], and local [43, 60–62] and global [38–40, 50] manifestation of anisotropy, to nonminimal inflationary models [51, 63–67]. The main advantage of the bottom-up approach is that we know exactly what feature of the CMB is being measured. One of its drawbacks is the plethora of different models and mechanisms that can be tested.

The second possibility is the top-down, or model-independent approach. Here, we are not concerned with the mechanisms responsible for deviations of gaussianity or SI, but rather with the qualitative features of any such deviation. Once these features are understood, we can use them as a guide for model building. Examples here include constructs of a posteriori statistics [29, 35, 36] and functional modifications of the two-point correlation function [28, 37, 68–71].

In the next Section, we will explore two different model-independent tests: one based on functional modifications of the two-point correlation function, and another one based on the so-called Maxwell's multipole vectors.

4. Multipole Vectors

Multipole vectors were first introduced in cosmology by Copi et al. [23] as a new mathematical representation for the primordial temperature field, where each of its multipoles ℓ are represented by ℓ unit real vectors. Later, it was realized that this idea is in fact much older [72], being proposed originally by J. C. Maxwell in his Treatise on Electricity and Magnetism.

The power of this approach is that the multipole vectors can be entirely calculated in terms of a temperature map, without any reference to external reference frames. This makes them ideal tools to test the morphology of CMB maps, like the quadrupole-octupole alignment.

The purpose of the following presentation is only comprehensiveness. A mathematically rigorous introduction to the subject can be found in references [72–74], as well as other papers in this review volume.

4.1. Maxwell's Representation of Harmonic Functions

We start our presentation of the multipole vectors by recalling some terminology. A harmonic function in three dimensions is any (twice differentiable) function ℎ that satisfies Laplace's equation∇²ℎ=0,(66) where ∇² is the Laplace operator. In spherical coordinates, the formal solution to Laplace's equation which is regular at the origin (𝑟=0) isℎ=^∞_ℓ=0ℎ_ℓ(𝑟,𝜃,𝜑),ℎ_ℓ=^ℓ_𝑚_=−ℓ𝑎_ℓ𝑚𝑟^ℓ𝑌_ℓ𝑚(𝜃,𝜑).(67) The functions 𝑟^ℓ𝑌_ℓ𝑚 are known as the solid spherical harmonics [74]. Since they agree with the usual spherical harmonics on the unit sphere, it is sometimes stated in the literature that the latter form a set of harmonic functions. This is an abuse of nomenclature though, and the reader should be careful.

Given the scalar nature of Laplace's operator, it is possible to find solutions to (66) in terms of Cartesian coordinates. Such solutions can be constructed by combining homogeneous polynomials (i.e., sums of monomials of the same order) of order ℓ:ℎ=^∞_ℓ=0ℎ_ℓ(𝑥,𝑦,𝑧),ℎ_ℓ=^ℓ_𝑎𝑏𝑐𝜆_𝑎𝑏𝑐𝑥^𝑎𝑦^𝑏𝑧^𝑐,(𝑎+𝑏+𝑐=ℓ).(68)

In three dimensions, the most general homogeneous polynomial of order ℓ contains (ℓ+2)!/(2!ℓ!) independent coefficients. However, since each polynomial must independently satisfy (66), precisely ℓ!/(2!(ℓ−2)!) of these coefficients will depend on each other. This constraint leaves us with (ℓ+1)(ℓ+2)/2−ℓ(ℓ−1)/2=2ℓ+1 independent degrees of freedom in each multipole ℓ—which is, of course, the same number of independent degrees of freedom appearing in (67).

Based on this analysis, Maxwell introduced his own representation of harmonic functions. He noticed that by successively applying directional derivatives of the form 𝑣⃗⋅∇≡∇_⃗𝑣 over the monopole potential 1/𝑟, where 𝑟=√𝑥²+𝑦²+𝑧² and 𝑣⃗ is a unit vector, he could construct solutions of the form (68). That is,𝑓_ℓ(𝑥,𝑦,𝑧)=𝜆_ℓ∇_𝑣ℓ⋯∇_𝑣2∇_𝑣11𝑟|||_𝑟=1,(69) where 𝜆_ℓ are real constants. We are now going to show that this construction does indeed lead to solutions of the form (68). First, note that there is a pattern which emerges from successive application of directional derivatives over the monopole function𝑓₀=1𝑓₁𝑟,=∇_𝑣1𝑓₀=−𝑣¹⋅𝑟⃗𝑟³𝑓₂,=∇_𝑣2𝑓₁=3𝑣²⋅𝑟⃗𝑣₁⋅𝑟⃗−𝑟²𝑣₁⋅𝑣₂𝑟³,(70) and so on. By induction, one can show that the general expression will be given by𝑓_ℓ=(−1)^ℓ(2ℓ−1)!!∏^ℓ𝑖=1𝑣_𝑖⃗⋅𝑟⃗+𝑟²𝑄_ℓ−2𝑟^2ℓ+1,(71) where 𝑄_ℓ−2 is a homogeneous polynomial of order ℓ−2 which only involves combinations of the vectors 𝑣_𝑖⃗ and 𝑟⃗.

The numerator of the function 𝑓_ℓ given by (71) is clearly a homogeneous polynomial of order ℓ (as one can easily check for some ℓs and also prove by mathematical induction). 𝐴 not so obvious result is that this polynomial is also harmonic. To prove that, let us define𝑔_ℓ=(−1)^ℓ(2ℓ−1)!!^ℓ_𝑖=1𝑣_𝑖⃗⋅𝑟⃗+𝑟²𝑄_ℓ−2(72) and consider the application of the operator ∇² over the combination 𝑟^𝛼𝑔_ℓ. For the 𝑥-component we get𝜕^2𝑥𝑟^𝛼𝑔_ℓ=𝑟^𝛼𝜕^2𝑥𝑔_ℓ+2𝛼𝑟^𝛼−2𝑥𝜕_𝑥𝑔_ℓ+𝛼𝑟^𝛼−2+𝛼(𝛼−2)𝑥²𝑟^𝛼−4𝑔_ℓ.(73) Repeating this process for the 𝑦- and 𝑧-components and then adding the results, we find ∇²𝑟^𝛼𝑔_ℓ=+𝑟^𝛼𝛼∇(²𝛼𝑔_ℓ++1)2𝑟𝛼𝑟^𝛼−2_𝛼−2𝑥𝜕_𝑥𝑔_ℓ+𝑦𝜕_𝑦𝑔_ℓ+𝑧𝜕_𝑧𝑔_ℓ𝑔_ℓ=𝑟^𝛼∇²𝑔_ℓ+𝛼(𝛼+2ℓ+1)𝑟^𝛼−2𝑔_ℓ,(74) where in the last step we have used Euler's theorem on homogeneous functions, that is, 𝑟⃗⋅∇𝑔_ℓ=ℓ𝑔_ℓ. If we now choose 𝛼=−(2ℓ+1), we find immediately that∇²𝑔^ℓ𝑟^2ℓ+1=∇²𝑔_ℓ𝑟^2ℓ+1.(75)

By construction, the left-hand side of the above expression is equal to ∇²𝑓_ℓ. But according to the definition (69), this quantity is also zero, since Laplace's operator commutes with directional derivatives and ∇²(1/𝑟)=0 for 𝑟>0. Therefore, ∇²𝑔_ℓ=0 and 𝑔_ℓ is harmonic, which completes our proof.

In conclusion, Maxwell's construction of harmonic functions, (69), is completely equivalent to the standard representation in terms of spherical harmonics. More importantly, this gives a one-to-one relationship between temperature maps (given by the 𝑎_ℓ𝑚's) and ℓ unit vectors 𝑣_𝑖⃗. This means that the multipole vectors can be directly calculated from a CMB map, without any reference to external reference frames or additional geometrical constructs. The reader interested in algorithms to construct the vectors from CMB maps may check references [23, 72], as well as the other papers in this volume that review this approach.

4.2. Multipole Vectors Statistics

It should be clear from the discussion above that the multipole vectors give an intuitive way to discover, interpret, and visualize phase correlations between different multipoles in the CMB maps. But they also can reveal intramultipole features, such as planarity (when a given multipole presents a preferred plane). Some of the most conspicuous hints of statistical anisotropy in the CMB, like the quadrupole-octupole alignment, have indeed been first found with less mathematically elegant methods [22, 27], but are best described in terms of the multipole vector formalism [23, 24, 26, 72]. However, this case should also sound an alarm, because some feature of a (presumably) random realization was found, then further scrutinized with a certain test which was, to some extent (intentionally or not), tailored to single out that very feature.

The pitfalls of aprioristic approaches are sometimes unavoidable, and all we can do is to take a second look at our sample with a more generic set of tools, to try to assess how significant our result really is in the context of a larger set of statistical tests. Multipole vectors are, in fact, ideally suited to this, since they can be found for all multipoles, and it is relatively easy to construct scalar combinations of these vectors with the usual methods of linear algebra. Also convenient is the fact that simulating these vectors from maps, or even directly, is also relatively easy, so the standard model of a Gaussian random field can be easily translated into the pdfs for the tests constructed with the multipole vectors. An important drawback of the multipole vectors is that, because they are computed in terms of equations which are nonlinear in the temperature fields, the distribution functions of statistical tests involving these vectors are highly non-Gaussian [75]. This means that these statistics usually have to be estimated by means of simulations (usually assuming that the underlying temperature field is itself Gaussian). Another delicate issue is how stable the multipole vectors are to instrumental noise and sky cut—and here, again, we must rely on numerical simulations to compare the observations with theoretical models.

In order to see how one should go about constructing a general set of tests, it should be noted first that the multipole vectors define directions, but they have no sense (they are “headless vectors”). This follows from the fact that the sign of the constants 𝜆_ℓ is degenerated with the sign of the vectors 𝑣_ℓ,𝑝 (𝑝=1,…,ℓ). Hence, the first requirement is that our tests should be independent of this sign ambiguity. Notice that this ambiguity extends to ancillary constructs such as the normal vectors, defined as the vector product: 𝑤_{ℓ,𝑝𝑝′}≡𝑣_ℓ,𝑝⊗𝑣_ℓ,𝑝′.(76) Mathematically, the sign ambiguity implies that the multipole vectors belong to the quotient space 𝑆²/ℤ₂ (also known as ℝ𝑃²) while the normal vectors belong to ℝ³/ℤ₂.

In principle, any positive-definite scalar constructed through the multipole or the normal vectors is “fair game”, but there are some guidelines; for example, one should avoid double-counting the same degrees of freedom. Below we review some of these tests, based on the work of [35].

4.2.1. The 𝑅 Statistic

The first example of a test involving the multipole vectors would be a scalar product, such as 𝑣⃗⋅𝑣^⃗. The most natural test would involve asking whether the ℓ multipole vectors of the given multipole ℓ are especially aligned or not. This means computing𝑅_ℓℓ=2ℓ(ℓ−1)^ℓ_{𝑝,𝑝′>𝑝}||𝑣_ℓ,𝑝⋅𝑣_ℓ,𝑝′||,(77) where the normalization was introduced to make 0≤𝑅_ℓℓ≤1.

This idea could be generalized to test alignments between multipole vectors at different multipoles: 𝑅_ℓℓ^′=1ℓℓ^_𝑝,𝑝′||𝑣_ℓ,𝑝⋅𝑣_ℓ^′_,𝑝′||.(78) In fact, the quadrupole-octupole alignment can already be seen with this simple test: for essentially all CMB maps the significance of the alignment as measured by 𝑅₂₃ is of the order of 90%–95% CL [35].

4.2.2. The 𝑆 Statistic

The second most natural test does not involve directly the multipole vectors themselves, but the normal vectors that can be produced by taking the vector product between the multipole vectors. So, we take 𝑤_{ℓ,𝑝𝑝′}=𝑣_ℓ,𝑝⊗𝑣_ℓ^′_,𝑝′.(79) Notice that the number of normal vectors for a given multipole ℓ is 𝑙=ℓ(ℓ−1)/2—so the number of normal vectors grows rapidly for larger multipoles, making it harder to use and meaning that the same degrees of freedom may be overcounted, at least for ℓ>3.

Again, the best strategy is simplicity: we can ask whether the normal vectors are aligned, within and between multipoles. This means computing 𝑆_ℓℓ=2𝑙(𝑙−1)^𝑙_{𝑝,𝑝′>𝑝}||𝑤_ℓ,𝑝⋅𝑤_ℓ,𝑝′||,(80) where the normalization was introduced again to make 0≤𝑆_ℓℓ≤1.

And yet again this idea is easily generalized to test alignments between normal vectors at different multipoles: 𝑆_ℓℓ^′=1𝑙𝑙^_𝑝,𝑝′||𝑤_ℓ,𝑝⋅𝑤_ℓ^′_,𝑝′||.(81) With this test, the quadrupole-octupole alignment is much more significant: for essentially all CMB maps the test 𝑆₂₃ deviates from the expected range with 98% CL [35].

4.2.3. Other Tests with Multipole Vectors

One can go on and expand the types of tests using both multipole and normal vectors. One idea would be, for example, to disregard the moduli of the normal vectors: 𝐷_ℓℓ^′=1𝑙𝑙^_𝑝,𝑝′||𝑤_ℓ,𝑝⋅𝑤_ℓ^′_,𝑝′||.(82) This test is therefore insensitive to the relative angle between the multipole vectors that produce any given normal vectors. This idea is similar to the planar modulations that will be discussed in the next section. With this test, the quadrupole-octupole alignment is significant to about 95%–98% CL This test can also be generalized to a self-alignment test (ℓ=ℓ^), with just an adjustment to the normalization.

Another possibility would be to measure the alignment between multipole vectors and normal vectors: 𝐵_ℓℓ^′=1ℓ𝑙^_𝑝,𝑝′||𝑣_ℓ,𝑝⋅𝑤_ℓ^′_,𝑝′||.(83) Of course, this test cannot be easily generalized to a self-alignment. With the 𝐵₂₃ test, the quadrupole-octupole alignment is significant to about 95%–98% CL

We could go on here, but it should be clear that all the information (the 2ℓ+1 real degrees of freedom for each multipole ℓ) has already been exhausted in the tests above.

5. Temperature Correlation Function

Despite their strong cosmological appeal, the multipole vectors have some limitations. Not only their directions in the CMB sky are sometimes difficult to interpret physically [26], but they also have the additional drawback of mixing, in a nontrivial manner, information on both gaussianity and SI of the map being analyzed [26, 75].

Another way of quantifying deviations in the standard statistical framework of cosmology is through functional modifications of the two-point correlation function [37, 68]. Although this approach does not offer an optimal separation between gaussianity and SI (which is, by the way, an open problem in this field), working with the two-point correlation function makes it easier to test Gaussian models of statistical anisotropy [28, 71].

The most general 2-point correlation function (2pcf) of two independent unit vectors is a function 𝐶 of the form𝐶∶𝑆²×𝑆²⟶ℝ.(84) We have seen in Section 2.3 that if we choose spherical coordinates (𝜃_𝑖,𝜑_𝑖) to describe each vector 𝑛_𝑖̂, the function above can be decomposed either in terms of two spherical harmonics 𝑌_{ℓ𝑖𝑚𝑖} or in terms of the bipolar spherical harmonic 𝒴^ℓ𝐿𝑀_1ℓ2. In any case, the 2pcf will have the following functional dependence:𝐶=𝐶𝜃₁,𝜑₁,𝜃₂,𝜑₂.(85)

This function is absolutely general. If the Universe has any cosmological deviation of isotropy, whatever it is, it can be described by the function above (see also Figure 5).

Unfortunately, this function will be of limited theoretical interest unless we have some hints on how to select its relevant degrees of freedom. This difficulty is in fact a general characteristic of model-independent tools, which at some stages forces us to rely on our theoretical prejudices about the statistical nature of the Universe in order to construct estimators of non-Gaussianity and/or statistical anisotropy. Nonetheless, it is still possible to construct statistical estimators of anisotropy based on (85). For example, Hajian and Souradeep [68] have constructed an unbiased estimator 𝜅_ℓ for this function based solely on the requirement that this estimator should be rotationally invariant. Although it is true that any statistically significant 𝜅_ℓ>0 will point towards anisotropy, it is not clear what type of anisotropy is being detected by this estimator.

We can still use (85) to search for deviations of SI if we restrict its domain to a smaller and nontrivial subdomain. For example, we can take the vectors 𝑛₁̂ and 𝑛₂̂ to be the same and expand a function of the form𝐶=𝐶𝜃₁,𝜑₁,(86) which is equivalent to 𝐶 : 𝑆²→ℝ. This form of the 2pcf makes it ideal for searching for power multipole moments in CMB, once a suitable estimator is defined [37]. Unfortunately, when we take 𝑛₁̂=𝑛₂̂, we are in fact considering a one-point correlation function, which by construction does not allow us to measure correlations between different points in the sky.

It seems in principle that the functions (86) and (85) are the only possibilities besides the isotropic 2pcf (24). If not, what other combinations of the vectors 𝑛₁̂ and 𝑛₂̂ can we consider? As a matter of fact, these two vectors are geometrical quantities intuitively bound to our notion of two-point correlation functions on the sphere. From this perspective they are not fundamental quantities. In fact, we can equally well represent the 2pcf by a disc living inside the unit sphere, as shown in Figure 4.

Figure 4 

Geometrical representation of the 2pcf in terms of a unit disc (or plane).

Figure 5 

Probability density function for the theoretical variable 𝑥. The shaded area is the probability (98).

In this representation, 𝜃 is the angle between the vectors 𝑛₁̂ and 𝑛₂̂, as usual. The normal to the plane, 𝑛̂, is represented by two spherical angles (Θ,Φ). Finally, there is an overall orientation 𝜑 of the disc around its unit vector which completes the four degrees of freedom contained in (85). We have found therefore another valid geometrical representation of the most general 2pcf:𝐶=𝐶(Θ,Φ,𝜃,𝜑).(87)

The main advantage of the above representation when compared to (85) is its straightforward geometrical interpretation. First, note that the angular separation 𝜃 of the isotropic 2pcf is trivially included in this definition and does not need to be obtained as a consequence of rotational invariance (see the discussion in Section 3.2.1.) Second, by characterizing the correlation function in terms of the geometrical components of a disc, we know exactly what are the degrees of freedom involved. This makes it easier to construct estimators of statistical anisotropy, alleviating the drawbacks of model-independent approaches mentioned above.

5.1. Anisotropy through Planarity

An immediate application of the representation (87) is its use in the search for planar deviations of isotropy in CMB [28, 71]. Planar modulations of astrophysical origin may play an important role to the CMB morphology. One example is the role played by the galactic and ecliptic plane in the quadrupole-octupole/north-south anomalies [26]. Also, it is well known that our galactic plane is sensible source of foreground contamination in the construction of cleaned CMB maps. These hints indicate that CMB modulations induced by the disc in Figure 4 are not only a mathematical possibility but perhaps also a symmetry of cosmological relevance.

Since we are primarily interested in measuring planar modulations of CMB, but including the usual angular modulation as the isotropic limit of the 2pcf, we can consider only the azimuthal average of (87):𝐶⟶12𝜋^02𝜋𝐶(Θ,Φ,𝜃,𝜑)𝑑𝜑.(88) The resulting function can be easily expand in terms of simple special functions as𝐶(Θ,Φ,𝜃)=_ℓ_𝑙,𝑚2ℓ+1√4𝜋𝒞^ℓ𝑙𝑚𝑃_ℓ(cos𝜃)𝑌_𝑙𝑚(Θ,Φ),𝑙∈2ℕ,(89) where the restriction on the 𝑙-mode results from the symmetry 𝑛₁̂↔𝑛₂̂.(See [71] for more details. Incidentally, the functional dependence of (89) differs from that originaly shown in [71], which does not correctly take the planar dependence into account. We thank Yuri Shtanov for pointing this out.) The multipolar coefficients 𝒞^ℓ𝑙𝑚 correspond to a generalization of the usual angular power spectrum 𝐶_ℓ's. In fact, they can be seen as the coefficients of a spherical harmonic decomposition of the function 𝐶_ℓ(𝑛)̂, provided that this function suffers modulations as we sweep planes on the sphere.

5.1.1. Angular-Planar Power Spectrum

Since we are restricting our analysis to the Gaussian framework, the set of coefficients 𝒞^ℓ𝑙𝑚 is all we need to characterize the two-point correlation function. However, the final product of CMB observations is temperature maps, and not correlation maps. What we need then is an algebraic relation between the multipolar coefficients 𝒞^ℓ𝑙𝑚 and the temperature coefficients 𝑎_ℓ𝑚 defined in (3). At first sight, this relation could be obtained by equating expression (89) to its standard definition in (22) and then using the orthogonality of the special functions to isolate the 𝒞^ℓ𝑙𝑚's in terms of the 𝑎_ℓ𝑚's. But for that to work we need the relation between the set of angles (Θ,Φ,𝜃) and (𝜃₁,𝜑₁,𝜃₂,𝜑₂) which, depending on the reference frame we choose, is extremely complicated. Fortunately, all we need is the following relation:𝑛₁̂⋅𝑛₂̂=cos𝜃=cos𝜃₁cos𝜃₂+sin𝜃₁sin𝜃₂cos𝜑₁−𝜑₂,(90) together with a suitable choice of our coordinate system. For example, we can use the invariance of the scalar product 𝑛₁̂⋅𝑛₂̂ and choose our coordinate system such that the discs of Figure 4 lies in the 𝑥𝑦 plane. With this choice we will have(Θ,Φ)=(0,0),cos𝜃=cos𝜑₁−𝜑₂,(91) and the integration over 𝜃 becomes simple. Once this is done, we make a passive rotation of the coordinate system and then we integrate over the remaining angles Θ and Φ, which will then be given precisely by the Euler angles used in the rotation. The details are rather technical and can be found in the Appendix. The final expression is(−1)^𝑚𝒞^ℓ𝑙𝑚√2𝑙+1=2𝜋_ℓ1,𝑚1_ℓ2,𝑚2𝑎_ℓ1𝑚2𝑎_ℓ2𝑚2𝑙ℓ¹ℓ₂−𝑚𝑚₁𝑚₂𝐼^{ℓ𝑙,ℓ1ℓ2},(92) where𝐼^ℓ𝑙,ℓ_1ℓ2≡_𝑚(−1)^𝑚𝜆_ℓ1𝑚𝜆_ℓ2𝑚𝑙ℓ¹ℓ₂0𝑚−𝑚^𝜋0𝑑(−cos𝜃)𝑃_ℓ(cos𝜃)𝑒^𝑖𝑚𝜃,(93) and where the 𝜆_ℓ𝑖𝑚's form a set of coefficients resulting from the 𝜃 integration, which are zero unless ℓ_𝑖+𝑚= even (see the Appendix for more details).

Expression (92) is what we were looking for. With this relation, the angular-planar power spectrum 𝒞^ℓ𝑙𝑚 can be calculated from first principles for any model predicting a specific covariance matrix. Moreover, since the angular-planar function (89) is, after all, a correlation function, it should be possible to relate the angular-planar power spectrum 𝒞^ℓ𝑙𝑚 to the bipolar power spectrum 𝒜^ℓ𝐿𝑀_1ℓ2 of Hajian and Souradeep. In fact, by inverting expression (30) and plugging the result in (92), we find a linear relation between these two set of coefficients:𝒞^ℓ𝑙𝑚=2𝜋_ℓ1,ℓ2𝐴^ℓ𝑙𝑚_1ℓ2(−1)^ℓ1−ℓ2(2𝑙+1)𝐼^ℓ𝑙,ℓ_1ℓ2.(94) Here, the set of geometrical coefficients 𝐼^ℓ𝑙,ℓ_1ℓ2 plays a similar role to the 3-j symbols in expression (30). Note also that the angular-planar power spectrum has only three free indices while the bipolar power spectrum has four. This is a consequence of the azimuthal average we took in (88), which further constrains the degrees of freedom of the correlation function.

5.1.2. Statistical Estimators and 𝜒² Analysis

We have shown that the angular-planar power spectrum (92) is given in terms of an ensemble average of temperature maps. Evidently, we cannot calculate it directly from data, for we have only one CMB map (the one taken from our own Universe.) The best we can do is to estimate the statistical properties of (92), like its mean and variance, and see whether these quantities agree, in the statistical sense, with what we would expect to obtain from a particular model of the Universe.

The reader should note that this procedure is not new—its limitation is due to the same cosmic variance which lead us to construct an estimator for the angular power spectrum 𝐶_ℓ (see the discussion of Section 2.4). For the same reason, we will need to construct an estimator for the angular-planar power spectrum. An obvious choice is𝒞^ℓ𝑙𝑚≡2𝜋√2𝑙+1_{ℓ1,𝑚1ℓ}_2,𝑚2𝑎_ℓ1𝑚1𝑎_ℓ2𝑚2𝑙ℓ¹ℓ₂𝑚𝑚₁𝑚₂𝐼^{ℓ𝑙,ℓ1ℓ2},(95) for in this case we have an unbiased estimator:𝒞^ℓ𝑙𝑚=𝒞^ℓ𝑙𝑚,(96)

regardless of the underlying model.

If we now have a model predicting the angular-planar power spectrum, we can ask how good this model fits the observational data once it is calculated using (95). All we need is a simple chi-square (𝜒²) goodness-of-fit test, which in our case can be written in the following generalized form:𝜒^2𝜈^𝑙ℓ≡2𝑙1+1^𝑙_𝑚=−𝑙|||𝒞^ℓ𝑙𝑚−𝒞^ℓ𝑙𝑚|||²𝜎^ℓ𝑙𝑚²,(97) in which ℓ and 𝑙 are the angular and planar degrees of freedom, respectively, and where 𝜎^ℓ𝑙𝑚 is just the standard deviation of the estimator 𝒞^ℓ𝑙𝑚. The (2𝑙+1)⁻¹ factor accounts for the 2𝑙+1 planar degrees of freedom and was introduced for latter convenience.

In Section 7 we will apply this test to the 5-year WMAP temperature maps in order to check the robustness of the ΛCDM model against the hypothesis of SI. Before that, we will stop and digress a little about how observational uncertainties should be included in our analysis.

6. Theory v. Observations: Cosmic and Statistical Variances

Until now we have been concerned with the formal aspects of non-Gaussian and statistically anisotropic universes, and how model-independent tests might be designed to detect deviations of either gaussianity or statistical isotropy. We will now discuss how such tests can be carried out and interpreted once we possess cosmological data.

In a great variety of tests, statistical tools are designed to detect particular deviations of gaussianity or SI from cosmological data like, for example, the CMB. The final outcome of these tests is usually a probability (a pure number), which should be interpreted as the chance that a Universe like ours might result from an ensemble of “equally prepared” Gaussian and SI universes. An anomaly in CMB is usually understood as a measure of how unlikely a particular feature of our Universe is according to this specific test. The multipole vector statistics, for example, when applied to a large number of simulated (Gaussian and SI) CMB maps, show that only ∼0.01% of these maps have a quadrupole-octupole alignment as strong as WMAP maps [23, 26].

There are two points to keep in mind when carrying this type of analysis. The first is that a particular “detection” may always turn out to be a statistical fluctuation revealed by one specific tool. The robustness of an anomaly then depends on the number of independent tests pointing to the same result. The second is that the implementation of statistical tools is sensitive to the way we extract information from data, requiring an accurate separation between cosmological signal and astrophysical/instrumental noise.

In this Section we present a critical review of the standard procedure used to implement cosmological tests. We show that it does not account for the intrinsic uncertainties of cosmological observations, which may possibly lead to an under/overestimation of anomalies. We then present a generalization of this process which naturally accounts for these uncertainties.

6.1. Standard Calculations

Suppose 𝑥 is a random variable predicted by a particular model and that cosmological observations of this quantity returned the value 𝑥₀. (Rigorously, 𝑥 is only one realization of a random variable 𝑋, which is a real-valued function defined on a sample space. By the same reason we should not call the 𝑎_ℓ𝑚's in (3) random variables, though we will stick to this nomenclature throughout this text.)We would like to calculate the probability, according to this model, that in a random Universe we would have 𝑥≤𝑥₀. Assuming that 𝒫_th is the (normalized) probability density function (pdf) of 𝑥, this probability is commonly defined to be𝑃_≤≡^𝑥0_−∞𝒫_th(𝑥)𝑑𝑥.(98) The probability of having 𝑥>𝑥₀ is then simply given by 𝑃_>=1−𝑃_≤. See Figure 5.

If 𝑃_≤ is found to be too small or, equivalently, too high, we might be tempted to interpret 𝑥₀ as “anomalous” according to this model. However, this definition of probability assumes that 𝑥₀ was measured with infinite precision, and so it says nothing about an important question we must deal with: typically, the measurement of 𝑥₀ has an uncertainty which needs to be folded into the final probabilities that the observations match the theoretical expectations.

Moreover, since no two equally prepared experiments will ever return the same value 𝑥₀, our measurements should also be regarded as random events. In a more rigorous approach, we would have to consider 𝑥₀ itself as one realization of a random variable, conditioned to the distribution of the signal. In the case of CMB, however, this would be only part of the whole picture, since the randomness of the measurements of 𝑥₀ should also be related to the way this data is reduced to its final form. This happens because different map cleaning procedures will lead to slightly different values for 𝑥₀. This difference induces a variance in the data which reflects the remaining foreground contamination of the temperature map. We will elaborate more on this point through a concrete example, after we show how (98) may be changed in order to include the indeterminacy of cosmological measurements.

6.2. Convolving Probabilities

The question we want to answer is how to calculate the probability of 𝑥 being smaller than our measurements when the latter are also random events. Let us suppose that our measurement is described by the random variable 𝑦 and that 𝑥₀ is its most probable value. The probability of having 𝑥≤𝑦 is simply the probability of having𝑧≡𝑥−𝑦≤0(99) for some particular realization of the variable 𝑦. It should be clear by now that if we know the pdf of 𝑧, the probability we are looking for is simply the area under this distribution for −∞≤𝑧≤0. The probability 𝑃_≤ of 𝑧 being smaller or equal to zero can be calculated as𝑃_≤≡𝑃{(𝑥,𝑦)∣𝑥−𝑦≤0}=_{𝑥−𝑦≤0}𝒫(𝑥,𝑦)𝑑𝑥𝑑𝑦.(100) Now, under the hypothesis of independence of 𝑥 and 𝑦 we have 𝒫(𝑥,𝑦)=𝒫_th(𝑥)𝒫_obs(𝑦), where 𝒫_obs is the (normalized) probability distribution function of the variable 𝑦. We can therefore rewrite the last expression as𝑃_≤=^∞−∞^𝑦−∞𝒫_th(𝑥)𝑑𝑥𝒫_obs(𝑦)𝑑𝑦.(101) If we now hold 𝑦 fixed and do the change of variables 𝑥=𝑢+𝑦, we get𝑃_≤=^0−∞^∞−∞𝒫_th(𝑢+𝑦)𝒫_obs(𝑦)𝑑𝑦𝑑𝑢,(102) where we have changed the position of the integrals. The reader will now notice that term inside brackets is precisely the pdf we were looking for. Since there is nothing special about the variable 𝑦, we can equally well hold 𝑥 fixed and repeat the calculus, obtaining the symmetric version of this result. In fact, the final pdf for 𝑧 is nothing else than the convolution of the pdfs of each variable [76]:𝒫(𝑧)≡𝒫_obs∗𝒫_th(𝑧)==^∞−∞𝒫_obs(𝑦)𝒫_th(𝑧±𝑦)𝑑𝑦_∞−∞𝒫_obs(𝑥∓𝑧)𝒫_th(𝑥)𝑑𝑥,(103) where the plus (minus) sign refers to the difference (sum) of 𝑥 and 𝑦. Integrating this pdf from (−∞,0] we get our answer𝑃_≤≡^0−∞𝒫(𝑧)𝑑𝑧.(104) As a consistency check, notice that in the limit where observations are made with infinite precision, 𝒫_obs(𝑦) becomes a delta function and we have𝑃_≤=^0−∞^+∞−∞𝛿𝑦−𝑥₀𝒫_th(𝑧±𝑦)𝑑𝑦𝑑𝑧=^±𝑥0_−∞𝒫_th(𝑥)𝑑𝑥(105) which agrees with our previous definition. The reader must be careful, though, not to think of (98) as some lower bound to (104). Since none of the pdfs appearing in (103) are necessarily symmetric, a large distance from 𝑥₀ to the most probable (not the mean!) value 𝑥 would not, by itself, constitute sufficient grounds to claim that the measured value of this observable is “unusual” in any sense, simply because a large overlap between the two pdfs can render the result usual according to (104).

To illustrate this point, let us calculate 𝑃_≤ using (98) and (104) for the pdf's which appear in Figure 6. For pedagogical reasons, we have chosen 𝒫_th and 𝒫_obs as positively (Maxwell distribution) and negatively (Gumbel distribution) skewed pdf's, respectively. The convolved distribution appears as the solid (black) line. For these pdf's, (98) gives 99.2%, while (104) gives 93.6% of chance of 𝑥 being smaller than the observed 𝑥₀; all pdf's were normalized to one.

Figure 6 

Probability density functions for 𝑥 (blue, dashed line), 𝑦 (red, dot-dashed line), and 𝑧 (solid line). The shaded area gives the probability of 𝑥 being smaller than 𝑥0 (dashed vertical line). See the text for more details.

7. A 𝜒² Test of Statistical Isotropy

Although the last example was constructed to emphasize an important feature of the formalism developed in Section 6, cosmological observables designed to measure deviations of either gaussianity or statistical isotropy will often follow asymmetric distributions. The intrinsic uncertainties of cosmological measurements, specially the ones originating from map cleaning procedures, may be crucial when searching for any map's anomalies. We will now make a concrete application of this formalism using the angular-planar chi-square test developed in Section 5.1.

7.1. ΛCDM Model

In the simplest realization of the ΛCDM model, the covariance matrix of temperature maps is determined by ⟨𝑎_ℓ1𝑚1𝑎_ℓ2𝑚2⟩=(−1)^𝑚2𝐶_ℓ1𝛿_ℓ1ℓ2𝛿_{𝑚1,−𝑚2}. Using this expression in (92), we find that𝒞=^{ℓ(𝑙𝑚SI)}𝐶_ℓ𝛿_𝑙0𝛿_𝑚0.(106) On the other hand, if the only nonzero 𝒞^ℓ𝑙𝑚's are given by 𝑙=𝑚=0, then 𝒞^ℓ00=𝐶_ℓ. Therefore, statistical isotropy is achieved if and only if the angular-planar power spectrum is of the form (106). Since we are only interested in nontrivial planar modulations, we will restrict our analysis to the cases where 𝑙≠0, that is,𝒞^ℓ𝑙𝑚=0,(𝑙≥2),(107) where, we remind the reader, the parity of ℓ comes from the symmetry (27).

For this particular model, we can also calculate the covariance matrix of the estimator (95) explicitly. Using the null hypothesis above, we find after some algebra that𝒞^ℓ𝑙𝑚^∗𝒞^{𝑙′𝑚′}_ℓ=8𝜋²_ℓ1,ℓ2𝐶_ℓ1𝐶_ℓ2𝐼^ℓ𝑙,ℓ₁^ℓ₂²𝛿_𝑙𝑙^′𝛿_𝑚𝑚′.(108) This matrix has some interesting properties. First, note that the planar degrees of freedom are independent in this case (which justifies the (2𝑙+1)⁻¹) factor introduced in (97)). Second, its diagonal elements are given by the variance (𝜎^ℓ𝑙𝑚)²=⟨(𝒞^ℓ𝑙𝑚)^∗𝒞^ℓ𝑙𝑚⟩, which now becomes 𝑚-independent:𝜎^ℓ𝑙𝑚²⟶𝜎^𝑙ℓ²=8𝜋²_ℓ1,ℓ2𝐶_ℓ1𝐶_ℓ2𝐼^ℓ𝑙,ℓ₁^ℓ₂².(109)

Therefore, for the particular case of the ΛCDM model, the chi-square test (97) gets even simpler. Using (107) and (109), we find𝜒^2𝜈^𝑙ℓ=2𝑙1+1^𝑙_𝑚=−𝑙|||𝒞^ℓ𝑙𝑚|||²𝜎^𝑙ℓ².(110) It is now clear that if the data under analysis are really described by this model, then it must be true that𝜒^2𝜈^𝑙ℓ=2𝑙1+1^𝑙_𝑚=−𝑙𝒞^ℓ𝑙𝑚^∗𝒞^ℓ𝑙𝑚𝜎^𝑙ℓ²=1,(111) where we have used (109). This shows that any large deviation of this test from unity will be an indication of planar modulation in temperature maps, up, of course, to error bars. For convenience, let us define a new quantity as𝜒^𝑙ℓ≡𝜒^2𝜈^𝑙ℓ−1(112) which will quantify anisotropies whenever 𝜒^𝑙ℓ is significantly positive or negative.

This generalized chi-square test furnishes a complete prescription when searching for planar modulations of temperature in CMB maps. We emphasize, though, that for a given CMB map, the chi-square analysis must be done entirely in terms of that map's data. Since we are performing a model-independent test, we are not allowed to introduce fiducial biases in the analysis (e.g., by calculating 𝜎^𝑙ℓ using 𝐶^model_ℓ), which would only include our a priori prejudices about what the map's anisotropies should look like. Since the 𝐶_ℓ's are, by construction, a measure of statistical isotropy, consequently, an “anomalous” detection of 𝐶_ℓ's is by no means a measure of statistical anisotropy, and it is this particular value that should be used in (112) if we want to find deviations of isotropy, regardless of how high/low it is.

7.2. Searching for Planar Signatures in WMAP

In order to apply the test (112) to the 5-year WMAP data [11, 14], we will define two new variables:𝑥≡𝜒²^𝑙ℓ(_th),𝑦≡𝜒²^𝑙ℓ(_obs),(113) which will be jointly analyzed using the formalism of Section 6. Still, there remains the question of how to obtain their pdfs. These functions can be obtained numerically provided that the number of realizations of each variable is large enough, since in this case their histograms can be considered as piecewise constant functions which approximate the real pdfs. For the case of the (theoretical) variable 𝑥 defined above we have run 2×10⁴ Monte Carlo simulations of Gaussian and statistically isotropic CMB maps using the ΛCDM best-fit 𝐶_ℓ's provided by the WMAP team [77]. With these maps we have then constructed 2×10⁴ realizations of the variable 𝑥.

The simulation of the (observational) variable 𝑦 is more difficult and depends on the way we estimate contamination from residual foregrounds in CMB maps. As is well known, not only instrumental noise, but also systematic errors (e.g., in the map-making process), the inhomogeneous scanning of the sky (i.e., the exposure function of the probe), or unremoved foreground emissions (even after applying a cut-sky mask) could corrupt—at distinct levels—the CMB data.

Foreground contamination, on the other hand, may have several different sources, many of which are far beyond our present scopes. However, since different teams apply distinct procedures on the raw data in order to produce a final map, we will make the hypothesis that maps cleaned by different teams represent—to a good extent—“independent” CMB maps. Therefore, we can estimate the residual foreground contaminations by comparing these different foreground-cleaned maps.

In fact, the WMAP science team has made substantial efforts to improve the data products by minimizing the contaminating effects caused by diffuse galactic foregrounds, astrophysical point sources, artifacts from the instruments and measurement process, and systematic errors [83, 84]. As a result, multifrequency foreground-cleaned full-sky CMB maps were produced, named Internal Linear Combination maps, corresponding to three- and five-year WMAP data [14, 79]. To account for the mentioned randomness, systematic, and contaminating effects of the CMB data, we will use in our analyses several full-sky foreground-cleaned CMB maps, listed in Table 1, which were produced using both the three- and five-year WMAP data.

Table 1 

Full-sky foreground-cleaned CMB maps from WMAP data used in our analysis to estimate the variable 𝑦 (see the text for more details). Note that the reference [78] includes the analysis of maps from the three- and five-year WMAP releases.

The prescription we adopt to determine the distribution of the observational variable 𝑦 is as follows: we simulate Gaussian random 𝑎_ℓ𝑚's in such a way that their central values are given by the five-year ILC5 data [14, 79], and with a variance which is estimated from the sample standard deviation of all the maps listed in Table 1. So, for example, suppose we have 𝑛 different full-sky temperature maps at hand and we want to estimate the randomness inherent in the determination of, let's say, 𝑎₃₂; therefore, we take𝒩𝑎^ILC532,𝜎₃₂⟶𝑎₃₂,(114) with𝜎₃₂=⎷1𝑛−1^𝑛_𝑖=1𝑎^𝑖32−𝑎₃₂²,𝑎₃₂=1𝑛^𝑛_𝑖=1𝑎^𝑖32,(115) where 𝒩(𝜇,𝜎) represents a Gaussian distribution with mean 𝜇 and standard deviation 𝜎. Note that if the residual contamination is indeed weak, then the sample variance above will be small, and our procedure will reduce to the standard way of calculating probabilities. As for the use of a Gaussian in (114), this choice was dictated not only by simplicity, but rather by the fact that the propagation of uncertainties in physical experiments is usually assumed to follow a normal distribution. However; note that there are some instrumental uncertainties, such as beam or gain uncertainties, which will not in general follow normal distributions. In fact, some of them may even fail to be additive, meaning that our convolution formula will be inapplicable in these cases. In our analysis, we have focused only on foreground residuals, where the normality hypothesis is reasonable. (We thank the referee for pointing out this important aspect of the formalism.)

7.3. Full-Sky Maps

Following this procedure, we have used the full-sky maps shown in Table 1 to construct 10⁴ Gaussian random 𝑎_ℓ𝑚's, which were then used to calculate 10⁴ realizations of 𝑦=(𝜒²)^𝑙ℓ(_obs). With those variables we constructed histograms which, together with the histograms for the (full-sky) variable 𝑥, were used to calculate the final probability (104). We have restricted our analysis to the range of values (ℓ,𝑙)∈[2,10], since the low multipolar sector (i.e., large angular scales) is where most of the anomalies were reported. The resulting histograms and pdf's are shown in Figure 7, and the final probabilities we obtained are shown in Table 2.

Table 2 

Final probabilities of obtaining, in a random ΛCDM Universe, a chi-square value smaller or equal to (𝜒²)^{𝑙ℓ(𝑜𝑏𝑠)}, as given by full-sky temperature maps.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 7 

Histograms for (𝜒2)𝑙ℓ(th) (blue), (𝜒2)𝑙ℓ(obs) (purple) and for the difference (𝜒2)𝑙ℓ(th)−(𝜒2)𝑙ℓ(obs) (solid, red line). We show only a few representative figures. The final probabilities are shown in Table 2 and correspond to the area under the solid curve from −∞ to 0. All pdfs are normalized to 1.

Overall, our results show no significant planar deviations of anisotropy in WMAP data. The most unlikely individual values in Table 2 are in the sectors (𝑙,ℓ) given by (2,5), (10,5), (4,7), and (6,8) and are all above a relative chance of 5% of either being too negative [(2,5), (10,5)] or too positive [(4,7), (6,8)]. However, it is perhaps worth mentioning that not only the individual values of (𝜒²)^𝑙ℓ are relevant, but also their coherence over a range of angular or planar momenta carries interesting information. So, for example, a set of (𝜒²)^𝑙ℓ's which are all individually within the cosmic variance bounds, but which are all positive (or negative), can be an indication of an excess (or lack) of planar modulation. This type of coherent behavior appears in the following cases: (𝜒²)^𝑙2, (𝜒²)^𝑙3 and, to a lesser extent, (𝜒²)^4ℓ (see Table 2). The angular quadrupole ℓ=2, as well as the angular octupole ℓ=3, have all positive planar spectra (for all values of 𝑙 which we were able to compute), indicated by probabilities larger than 50%. The planar hexadecupole 𝑙=4 also has 8 out of 9 angular spectra assuming positive values (only ℓ=5 is negative).

The data analyzed in this Section relates to the full-sky maps, which are certainly still affected by residual galactic foregrounds. The reader interested in the complete analysis, including data from masked CMB maps, can check reference [28].

8. Conclusions

We know our Universe is not perfectly Gaussian, homogeneous, or isotropic. The deviations from an idealized picture (or the lack thereof), whether predictably small or surprisingly large, can tell us a great deal about the Universe we live in. Since the types of physical mechanisms behind deviations from perfect gaussianity, homogeneity, or isotropy are typically very different, we should try to measure these individual features separately—whenever possible or practical.

Recently, it has been suggested that some of the most discussed anomalies in the CMB can be explained away [85], or that the evidence for them is statistically weak [86]. But even if it turns out that our Universe is a plain vanilla kind of place, where everything goes according to the inflationary theorist's dreams, we would still need to analyze it with tools that allow us to check the standard picture against the data. In addition, local physics (related to the solar system, or our galaxy), as well as instrumental quirks, tend to leave imprints on the CMB which are clearly anisotropic, but have a certain coherence which can be detected, and possibly corrected for, with the help of these checks.

However, in an era where at least the large-scale maps of the CMB are likely to remain basically unchanged, we should be careful not to over analize the data with the benefit of an ever greater hindsight (put another way, a posteriori conundrums only get worse with time.) This can only be achieved if we find natural and generally agreed-upon classifications of the types of deviations that may occur, without too much guidance from what the data is telling us. We believe that focusing on the possible underlying symmetries, with perhaps some guidance from group-theoretic arguments, is one way to settle these issues. We have presented a few methods along these lines, one using multipole vectors, the other using a natural generalization of the two-point correlation function (also other methods have been presented in this paper.

Perhaps the best indication that we are on the right track is the fact that most of these methods are applicable in other areas of physics and astronomy, and that in some cases we have adapted tests of anisotropy from other areas, such as scattering theory and the theory of angular momentum in quantum mechanics. So, even if these anomalies eventually perish, they will be survived by the powerful methods that have been devised to test them.

Appendix

A. Geometrical Identities and Derivations

A.1. Wigner 𝐷-Functions

From the unitarity of the rotation operators 𝐷(𝑅), we have_𝑚𝐷^ℓ𝑚_′𝑚𝜔₁𝐷^ℓ𝑚_′′𝑚𝜔₂=𝐷^ℓ𝑚_{′𝑚′′}𝜔₁𝜔₂.(A.1) If 𝜔₂=𝜔¹⁻¹, then using the identity 𝐷^ℓ𝑚_′′𝑚(𝜔⁻¹)=𝐷^𝑚∗ℓ_′′𝑚(𝜔), we find_𝑚𝐷^ℓ𝑚_′𝑚𝜔₁𝐷^𝑚∗ℓ_′′𝑚𝜔₁=𝛿_𝑚^′_𝑚^′′.(A.2)

A.2. Gaunt Integral

The definition of the Gaunt integral used in this paper is𝒢^ℓ1ℓ2ℓ3_{𝑚1𝑚2𝑚3}=(−1)^𝑚1𝑑²𝑛̂𝑌_{ℓ1,−𝑚1}(𝑛)̂𝑌_ℓ2𝑚2(𝑛)̂𝑌_ℓ2𝑚2(𝑛)̂.(A.3)

A.3. 3-j Symbols

We present here some useful identities related to the 3-j symbols (i) Isotropic limit𝑙¹𝑚₁𝑙₂0𝑚₂0=√(−21𝑙)₁+1^{𝑙2−𝑚1}𝛿^𝑙1𝑙2𝛿_{𝑚1,−𝑚2}.(A.4)(ii) Parity and permutations 𝑙¹𝑚₁𝑙₂𝑙𝑚₂𝑚=𝑙𝑙¹𝑙₂=(−𝑚1)𝑚₁𝑚_{2𝑙1+𝑙2+𝑙}𝑙²𝑚₂𝑙₁𝑙=(−1)_{𝑙1+𝑙2+𝑙}𝑚₁𝑚𝑙¹−𝑚₁𝑙₂𝑙−𝑚₂−𝑚.(A.5)(iii) Orthogonality ^𝑙1_{𝑚1=−𝑙1}^𝑙2_{𝑚2=−𝑙2}𝑙¹𝑙₂𝑙₃𝑚₁𝑚₂𝑚₃𝑙¹𝑙₂𝑙^3=𝛿_{𝑙3𝑙′3}𝑚₁𝑚₂𝑚^3𝛿_𝑚3𝑚^′32𝑙₃+1_𝑙2+𝑙3,_{𝑙1=||𝑙2−𝑙3||𝑙}¹_{𝑚1=−𝑙1}(2𝑙+1)𝑙¹𝑙₂𝑙₃𝑚₁𝑚₂𝑚₃𝑙¹𝑙₂𝑙₄=𝛿_{𝑚2𝑚′2}𝑚₁𝑚^2𝑚^3𝛿_𝑚3𝑚^′3,_𝑙𝑚_=−𝑙(−1)^𝑙−𝑚=√𝑚𝑙−𝑙𝑚ℓ02𝑙+1𝛿_ℓ,0.(A.6)

The last expression is particularly useful in the derivation of (106).

A.4. Derivation of (92)

We start by equating expressions (89) and(22)_ℓ0𝑥0200𝑑_𝑙,𝑚2ℓ+1√4𝜋𝒞^ℓ𝑙𝑚=𝑃_ℓ(cos𝜗)𝑌_𝑙𝑚(𝑛)_ℓ1,𝑚1̂_ℓ_2,𝑚2𝑎^ℓ₁^𝑚₁𝑎^∗ℓ_2𝑚2𝑌^ℓ₁^𝑚₁𝑛₁̂𝑌_ℓ2𝑚2𝑛₂̂.(A.7) As mentioned in the main text, the inversion of 𝒞^ℓ𝑙𝑚 as a function of the 𝑎_ℓ𝑚's is not a trivial task, since the angles (Θ,Φ,𝜃) depend nonlinearly on the angles (𝜃₁,𝜑₁,𝜃₂,𝜑₂). The easiest way to achieve this goal is to pick up a coordinate system where only the 𝜃 dependence is present. After integrating it out, we rotate our coordinate system using three Euler angles to recover back the (Θ,Φ) dependence, which can then be integrated with the help of some Wigner matrices identities. We start by positioning the vectors 𝑛₁̂ and 𝑛₂̂ in the 𝑥𝑦 plane, that is; we chose 𝑛₁̂=(𝜋/2,𝜙₁), 𝑛₂̂=(𝜋/2,𝜙₂). By (90), we then have cos𝜗=cos(𝜙₁−𝜙₂). Using the relation [41]𝑌_ℓ𝑚𝜋2,𝜙=𝜆^ℓ𝑚𝜆_ℓ𝑚𝑒^𝑖𝑚𝜙,=⎧⎪⎨⎪⎩(−1)^(ℓ+𝑚)/22ℓ+10,4𝜋(ℓ+𝔄!𝑚!)!!(ℓ−𝔅!𝑚!)!!,iofthℓer+w𝑚ise∈2ℕ,(A.8) where 𝔄 denotes (ℓ+𝑚−1) and 𝔅 denotes (ℓ−𝑚−1), we can integrate the 𝜃 dependence on both sides of (A.7). This gives us1√𝜋_𝑙,𝑚𝒞^ℓ𝑙𝑚𝑌_𝑙𝑚(0,0)=_ℓ1,𝑚10𝑥0200𝑑_ℓ2,𝑚2𝑎^ℓ₁^𝑚₁𝑎^∗ℓ_2𝑚2𝐼^ℓℓ₁^𝑚₁^ℓ₂^𝑚₂,(A.9) where we have introduced the following definition:𝐼^ℓℓ_{1𝑚1ℓ2𝑚2}≡−𝜆_ℓ1𝑚1𝜆_ℓ2𝑚2×^𝜋0𝑃_ℓcos𝜑₁−𝜑₂𝑒^{𝑖(𝑚1𝜑1−𝑚2𝜑2)}𝑑cos𝜑₁−𝜑₂.(A.10)

We need now to integrate out the Θ and Φ dependence in the right-hand side of (A.9) which was hidden due to our choice of a particular coordinate system. In order to do that, we keep the vectors 𝑛₁̂ and 𝑛₂̂ fixed and make a rotation of our coordinate system using three Euler angles 𝜔={𝛼,𝛽,𝛾}. This rotation changes the coefficients 𝒞^ℓ𝑙𝑚's and 𝑎_ℓ𝑚's according to𝑎_ℓ𝑚=_𝑚′𝐷^ℓ𝑚𝑚_′(𝜔)𝑎̃_ℓ𝑚′,𝒞^ℓ𝑙𝑚=_𝑚′𝐷^𝑙𝑚𝑚_′(𝜔)𝒞^𝑙𝑚′_ℓ,(A.11) where 𝒞^ℓ𝑙𝑚 and 𝑎̃_𝑙𝑚 are the multipolar coefficients in the new coordinate system, and where 𝐷^𝑙𝑚𝑚_′(𝜔) are the elements of the Wigner rotation matrix. The advantage of positioning the vectors 𝑛₁̂ and 𝑛₂̂ in the plane 𝑥𝑦 are that now the angles Θ and Φ are given precisely by the Euler angles 𝛽 and 𝛾, regardless of the value of 𝛼_𝑙,𝑚𝒞^ℓ𝑙𝑚𝑌_𝑙𝑚(0,0)=_𝑙,𝑚′𝒞^𝑙𝑚′_ℓ_𝑚=𝐷^𝑙𝑚𝑚_′(𝛼,𝛽,𝛾)𝑌_𝑙𝑚(0,0)_𝑙,𝑚′𝒞^𝑙𝑚′_ℓ𝑌_𝑙,−𝑚(𝛽,𝛾),(A.12) where in the last step we have used 𝑌_𝑙𝑚(0,0)=√(2𝑙+1)/4𝜋𝛿_𝑚0. Therefore, in our new coordinate system we have (dropping the “”in our notation)12𝜋_𝑙,𝑚=𝒞^ℓ𝑙𝑚𝐷^𝑙0𝑚(𝜔)√2𝑙+1_{ℓ1,𝑚1ℓ}_2,𝑚2𝑎^ℓ₁^𝑚₁𝑎^∗ℓ_2𝑚2_{𝑚′1𝑚′2}𝐼^ℓℓ_{1𝑚′1ℓ2𝑚′2}𝐷^ℓ1_{𝑚′1𝑚1}(𝜔)𝐷^ℓ2∗𝑚_′2𝑚2(𝜔).(A.13) We may now isolate 𝒞^ℓ𝑙𝑚 using the identities [41] 𝑑𝜔𝐷^𝑙1∗𝑚₁^𝑚′1=8𝜋₂(𝜔)𝐷^𝑙2_{𝑚2𝑚′2}(𝜔)2𝑙₁+1𝛿^𝑙1𝑙2𝑑𝜔𝐷_{𝑙1𝑚′1𝑚1}𝛿_𝑚1𝑚2𝛿_𝑚^′1_𝑚^′2,=8𝜋₂(𝜔)𝐷^𝑙2_{𝑚′2𝑚2}(𝜔)𝐷^𝑙3_{𝑚′3𝑚3}(𝜔)𝑙¹𝑙₂𝑙₃𝑚^1𝑚^2𝑚^3𝑙¹𝑙₂𝑙₃𝑚₁𝑚₂𝑚₃,(A.14) where 𝑑𝜔=sin𝛽𝑑𝛽𝑑𝛼𝑑𝛾, to obtain 1√2𝑙+1𝒞^ℓ𝑙𝑚=2𝜋_{ℓ1,𝑚1ℓ}_2,𝑚2𝑎^ℓ₁^𝑚₁𝑎^∗ℓ_2𝑚2_{𝑚′1𝑚′2}𝐼^ℓℓ_{1𝑚′1ℓ2𝑚′2}(−1)^{𝑚2+𝑚′2+𝑚}×ℓ¹𝑚_1ℓ₂𝑙−𝑚^20ℓ¹𝑚₁ℓ₂𝑙−𝑚₂−𝑚.(A.15) If we now redefine −𝑚₂→𝑚₂, and note that the first 3-j symbol above is identically zero unless 𝑚^1=𝑚^2, we finally obtain (92).

A.5. Some Properties of the Integral (93)

The geometrical coefficients 𝐼^ℓ𝑙,ℓ_1ℓ2 defined in (93) have many interesting properties which can be explored in order to speed up numerical computation of (92). First, we note that it is symmetric under permutation of ℓ₁ and ℓ₂𝐼^ℓ𝑙,ℓ_1ℓ2=_𝑚𝐼^ℓℓ_{1𝑚ℓ2𝑚}(−1)^𝑚ℓ¹=𝑚−ℓ𝑚₂0_𝑚𝑙𝐼^ℓℓ_{2𝑚ℓ1𝑚}(−1)^{𝑚+ℓ1+ℓ2+𝑙}ℓ²=−𝑚ℓ𝑚₁0_𝑚𝑙𝐼^ℓℓ_{2𝑚ℓ1𝑚}(−1)^{𝑚+2(ℓ1+ℓ2+𝑙)}ℓ²=𝐼_{ℓ𝑙,ℓ2ℓ1}𝑚−ℓ𝑚₁0𝑙.(A.16) Some of the other properties are a consequence of the integral 𝐼^ℓℓ_{1𝑚ℓ2𝑚} defined in (A.10). We may note for example that, due to the symmetry of the 𝜆_ℓ𝑚 coefficient defined in (A.8), we will have𝐼^ℓ𝑙,ℓ_1ℓ2=0,foranyℓ₁,ℓ₂∈ℕ∣ℓ₁+ℓ₂=odd.(A.17) Furthermore, the 𝜆_ℓ𝑚 coefficients restrict the 𝑚 summation above to their values which obey 𝑚+ℓ₁+ℓ₂=even. If we further notice that (A.10) is proportional to the integral of an integral of the form ∫¹⁻¹𝑃_ℓ(cos𝜃)cos𝑚𝜃𝑑𝜃, and that this integral is zero unless ℓ+𝑚=even, we conclude that𝐼^ℓ𝑙,ℓ_1ℓ2=0,foranyℓ₁,ℓ₂,ℓ∈ℕ∣ℓ₁+ℓ₂+ℓ=odd.(A.18) Besides, using the fact that the integral ∫¹⁻¹𝑃_ℓ(cos𝜃)⋅cos𝑚𝜃𝑑𝜃 is zero for any 𝑚<𝑙, we find𝐼^ℓ𝑙,ℓ_1ℓ2=0,foranyℓ₁,ℓ₂,ℓ∈ℕ∣ℓ₁<𝑙,ℓ₂<𝑙.(A.19) We finally comment on the special case where 𝑙=0, for which we have𝐼^ℓ0,ℓ_1ℓ2=(−1)^ℓ1√2ℓ₁+1_𝑚𝐼^ℓℓ_{1𝑚ℓ1𝑚}𝛿^ℓ1ℓ2.(A.20) However, ^ℓ^′_{𝑚=−ℓ′}𝐼^ℓℓ_{′𝑚ℓ′𝑚}=×^𝜋0⎛⎜⎜⎝𝑃_ℓ(cos𝜗)_ℓ′𝑚_=−ℓ′2ℓ^+14𝜋ℓ^𝔇+!𝑚!!!ℓ^𝔊−!𝑚!!!𝑒^𝑖𝑚𝜗=2ℓ_⎞⎟⎟⎠𝑑(−cos𝜗)+14𝜋¹⁻¹=1𝑃_ℓ(𝑥)𝑃_ℓ^′(𝑥)𝑑𝑥2𝜋𝛿^ℓℓ′,(A.21) where 𝔇 denotes (ℓ^+𝑚−1), and 𝔊 denotes (ℓ^−𝑚−1), and in the derivation above we have made use of the Fourier series expansion of the Legendre polynomial. So we conclude that𝐼^ℓ0,ℓ_1ℓ2=(−1)^ℓ12𝜋√2ℓ₁+1𝛿^ℓℓ1𝛿_ℓ1ℓ2,(A.22) which is needed in the derivation of (106).

Acknowledgments

The authors would like to thank Glenn Starkman and Yuri Shtanov for their enlightening discussions during the development of this paper. This paper was supported by FAPESP and CNPq.