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

Authors: L. Raul Abramo; Thiago S. Pereira
Journal: Advances in Astronomy
 (2010)

---

## 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”.


---

## Body

## 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 some10-35 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 (z≈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:(1)Θ(n̂)≡ΔT(n̂;η0)T(η0)=∫d3x′∫0η0dη′∑iKi(x⃗′,η′;n̂)Si(x⃗′,η′),
where η′ is conformal time, and 𝒮i denote the initial conditions of all matter and metric fields (as well as their time derivatives, if the initial conditions are nonadiabatic). Here Ki 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 wavenumberk=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 GH2≃10-12, where G is Newton's constant and H 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 ofad 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 ofaprioristic 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 (Section2). 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 directionn̂ in the CMB sky. So let us consider the dimensionless function on S2 as follows:(2)Θ(n̂)≡T(n̂)-T0T0,
where T0=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 vectorn⃗, 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 S2, for which the following decomposition holds:(3)Θ(n̂)=∑ℓ,maℓmYℓm(n̂).
Since the spherical harmonics Yℓm(n̂) obey the symmetry Yℓm*(n̂)=(-1)mYℓ,-m(n̂), the fact that the temperature field is a real function implies the identity aℓm*=(-1)maℓ,-m. 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-sightn̂ would reflect simply the density and gravitational potentials (the “sources”) at the position Rn̂, where R 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ηR, the linear kernels of (1) reduce to Ki(x⃗′,η′;n̂)→βiδ(η′-ηR)δ(x⃗′-n̂R), where R=η0-ηR and βi are constant coefficients. The photon distribution that we measure on Earth would therefore be given by(4)Θ(n̂)≈∑i‍βiSi(x⃗′=n̂R,η′=ηR).
We can also express this result in terms of the Fourier spectrum of the sources as follows:(5)Θ(n̂)≈∑i‍βi∫d3k(2π)3eik⃗·n̂RSi(k⃗,ηR).
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:(6)eik⃗·x⃗=4π∑ℓmiℓjℓ(kx)Yℓm*(k̂)Yℓm(x̂),
where jℓ(z) are the spherical Bessel functions. Substituting (6) into (5) we obtain that(7)aℓm=∫d2n̂Yℓm*(n̂)Θ(n̂)≈∫d3k(2π)3Θ(k⃗)×4πiℓjℓ(kR)Yℓm*(k̂),
where we have loosely collected the sources into the term Θ(k⃗)≡∑iβiSi(k⃗,ηR). 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 aℓm are inherited from those of the Fourier modes Θ(k⃗) of the underlying matter and metric fields. Notice that the properties of the aℓms 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 Yℓm.


### 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)C(k⃗,k′⃗)≡〈Θ(k⃗)Θ(k′⃗)〉.
In the absence of any symmetries, this would be a generic function of the arguments k⃗ and k′⃗, with only two constraints: first, because Θ(x⃗) is a real function, Θ*(k⃗)=Θ(-k⃗), hence, in our definition C*(k⃗,k′⃗)=C(-k⃗,-k′⃗); second, due to the associative nature of the expectation value, C(k⃗,k′⃗)=C(k′⃗,k⃗). It is obvious how to generalize this definition to 3, 4, or an arbitrary number of fields at different k⃗s (or “points”).Let us first discuss the issue of gaussianity. If we say that the variablesΘ(k⃗) are Gaussian random numbers, then all the information that characterizes their distribution is contained in their two-point function C(k⃗,k′⃗). The probability distribution function (pdf) is then formally given by (9)P[Θ(k⃗),Θ(k′⃗)]~exp[-Θ(k⃗)Θ(k′⃗)2C(k⃗,k′⃗)].
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:  (10)〈Θ(k⃗1)Θ(k⃗2)⋯Θ(k⃗2N)〉G=∑i,j∏α=1NBi,jα〈Θ(k⃗i)Θ(k⃗j)〉,
where the sum runs over all permutations of the pairs of wave vectors and Bi,j 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 theN-point functions in real space, we should have the following  (11)〈Θ(x⃗1)Θ(x⃗2)⋯Θ(x⃗N)〉→Homog.CN(x⃗1-x⃗2,…,x⃗N-1-x⃗N).
Writing this expression in terms of the Fourier modes, we get the following(12)〈Θ(x⃗1)Θ(x⃗2)⋯Θ(x⃗N)〉=∫d3k1d3k1⋯d3kN(2π)3Ne-ik⃗1·x⃗1e-ik⃗2·x⃗2⋯e-ik⃗N·x⃗N×〈Θ(k⃗1)Θ(k⃗2)⋯Θ(k⃗N)〉.
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 δ(k⃗1+k⃗2+⋯+k⃗N). In other words, the hypothesis of homogeneity constrains the N-point function in Fourier space to be of the following form:  (13)〈Θ(k⃗1)Θ(k⃗2)⋯Θ(k⃗N)〉H=(2π)3Ñ(k⃗1,k⃗2,…,k⃗N)×δ(k⃗1+k⃗2+⋯+k⃗N).
Notice that the “(N-1)-spectrum” in Fourier space, Ñ, can still be a function of the directions of the wavenumbers k⃗i (it will be, in fact, a function of N-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 theN-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 k⃗i in (12). Let us instead rewrite the δ-functions in terms of triangles, so for the 4-point function we have  (14)δ(k⃗1+k⃗2+k⃗3+k⃗4)=∫d3qδ(k⃗1+k⃗2-q⃗)δ(k⃗3+k⃗4+q⃗),
whereas for the 5-point function we have  (15)δ(k⃗1+k⃗2+k⃗3+k⃗4+k⃗5)=∫d3qd3q′δ(k⃗1+k⃗2-q⃗)δ(q⃗+k⃗3-q⃗′)×δ(q⃗′+k⃗4+k⃗5),
and so on, so that the N-point δ-function is reduced to N-2 triangles with N-3 “internal momenta” (the idea is nicely illustrated in Figure 1.) Substituting the expression for the N-point δ-function into (12) and integrating out all external momenta but the first (k⃗1) and last (k⃗N), the result is as follows(16)〈Θ(x⃗1)Θ(x⃗2)⋯Θ(x⃗N)〉=1(2π)3N∫d3k1d3q1⋯d3qN-3d3kN×eik⃗1·(x⃗1-x⃗2)eiq⃗1·(x⃗2-x⃗3)⋯eiq⃗N-3·(x⃗N-2-x⃗N-1)×eik⃗N·(x⃗N-1-x⃗N)〈Θ(k⃗1)Θ(q⃗1-k⃗1)⋯Θ(k⃗N)〉.
This expression shows explicitly that the real-space N-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 theN-point functions that come from isotropy alone? Clearly, no dependence on the directions defined by the points, x⃗i-x⃗j, can arise in the final expression for the N-point functions in real space, so from (12) we see that the N-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 theN-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,Θ(n̂), which are taken over the unit sphere S2. We will discuss mostly the two-point function here, and we defer a fuller discussion of N-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  (17)〈aℓmaℓ′m′*〉=∫d3kd3k′(2π)6(4π)2iℓ(-i)ℓ′jℓ(kR)jℓ′(k′R)Yℓm(k̂)×Yℓ′m′*(k̂′)〈Θ(k⃗)Θ*(k⃗′)〉.
Under the hypothesis of homogeneity, this expression simplifies considerably, leading to  (18)〈aℓmaℓ′m′*〉H=∫d3k2πiℓ(-i)ℓ′jℓ(kR)jℓ′(kR)Yℓm(k̂)Yℓ′m′*(k̂)×Ñ2(k⃗).
If, in addition to homogeneity, we also assume isotropy, then Ñ2→P(k), and the integration over angles factors out, leading to the orthogonality condition for spherical harmonics as follows (19)∫d2k̂Yℓm(k̂)Yℓ′m′*(k̂)=δℓℓ′δmm′,
and as a result the covariance of the aℓms becomes diagonal as follows(20)〈aℓmaℓ′m′*〉H,I=δℓℓ′δmm′∫dkkjℓ2(kR)2πk3P(k)=4πδℓℓ′δmm′∫dlogkjℓ2(kR)ΔT2(k)≡Cℓδℓℓ′δmm′,
where we have defined the usual temperature power spectrumΔT(k)=k3P(k)/2π2 in the middle line, and the angular power spectrum Cℓ 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 logk as follows:  (21)〈Θ(x⃗)Θ(x⃗′)〉H,I=∫dlogksin(k|x⃗-x⃗′|)k|x⃗-x⃗′|ΔT2(k).In an analogous manner to what was done above, we can also construct theangular two-point correlation function in harmonic space as follows:(22)〈Θ(n̂)Θ(n̂′)〉=∑ℓm∑ℓ′m′〈aℓmaℓ′m′*〉Yℓm(n̂)Yℓm*(n̂′).
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:(23)〈Θ(n̂)Θ(n̂′)〉H,I=∑ℓm‍∑ℓ′m′Cℓδℓℓ′δmm′Yℓm(n̂)Yℓm*(n̂′)=∑ℓCℓ2ℓ+14πPℓ(n̂·n̂′).
Clearly, not only is this expression the analogous in S2 of (21), but in fact the Fourier power spectrum ΔT2(k) and the angular power spectrum Cℓ are defined in terms of each other as indicated in (23) as follows:  (24)Cℓ=4π∫dlogkjℓ2(kR)ΔT2(k).
Now, using the facts that the spherical Bessel function of order ℓ peaks when its argument is approximately given by ℓ, and that ∫dlogzjℓ2(z)=1/(2ℓ(ℓ+1)), we obtain the following (this is one type of  what has become known in the literature as Limber's approximations)(25)Cℓ≈2πℓ(ℓ+1)ΔT2(k=ℓR).
Incidentally, from this expression it is clear why it is customary to define (26)𝒞ℓ≡ℓ(ℓ+1)2πCℓ≈ΔT2(k=ℓR).Using (12), we can easily generalize the results of this subsection  to N-point functions in S2 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 N-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)C(n̂1,n̂2)=C(n̂2,n̂1).
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, (ℓ1,m1;ℓ2,m2), we could equally well have used the basis of total angular momentum (L,M;ℓ1,ℓ2) and decomposed that expression as  (28)C(n̂1,n̂2)=∑L,M∑ℓ1,ℓ2𝒜ℓ1ℓ2LM𝒴ℓ1ℓ2LM(n̂1,n̂2),
where 𝒴ℓ1ℓ2LM are known as the bipolar spherical harmonics, defined by [41](29)𝒴ℓ1ℓ2LM(n̂1,n̂2)=[Yℓ1(n̂1)⊗Yℓ2(n̂2)]LM,
where L and M=m1+m2 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(30)𝒜ℓ1ℓ2LM=∑m1m2〈aℓ1m1aℓ2m2〉(-1)M+ℓ1-ℓ22L+1(ℓ1ℓ2Lm1m2-M),
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 aℓms, 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 aℓms. 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 aℓm as statistically independent Gaussian random variables, and define a statistical estimator for their variances as the average(31)Ĉℓ≡12ℓ+1∑m=-ℓℓ|aℓm|2.
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 have(32)limℓ→∞Ĉℓ=Cℓ.
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 Cℓ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 spectrumCℓ 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 aℓms are Gaussian, then we can rewrite expression (31) as(33)(2ℓ+1)CℓĈℓ=Xℓ,Xℓ=∑m=-ℓℓ|aℓm|2Cℓ,
where Xℓ is a chi-square random variable with 2ℓ+1 degrees of freedom. According to the central limit theorem, when ℓ→∞, Xℓ approaches a standard normal variable (A standard normal variable is a Gaussian variable X with zero mean and unit variance. Any other Gaussian variable Y with mean μ and variance σ can be obtained from X through Y=σX+μ.) which implies that Ĉℓ will itself follow a Gaussian distribution. Its mean can be easily calculated using (20) and (31) and is of course given by(34)〈Ĉℓ〉=Cℓ,
which shows that the Ĉℓs are unbiased estimators of the Cℓs. It is also straightforward to calculate its variance (valid for any ℓ)(35)〈(Ĉℓ-Cℓ)(Ĉℓ′-Cℓ′)〉=22ℓ+1Cℓ2δℓℓ′.
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. Ĉℓ is therefore the best estimator we can build to measure the statistical properties of the multipolar coefficients aℓm when both statistical isotropy and gaussianity hold.In later Sections, we will explore angular or harmonicN-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 aℓms. The estimator for the angular power spectrum, Ĉℓ, takes into account all the aℓms by dividing them into the different ℓs and summing over all m∈(-ℓ,ℓ). Clearly, it will inherit a sample variance for small ℓs, when the aℓms can only be divided into a few “independent parts”. As we try to estimate higher-order objects such as the N-point functions, we will have to subdivide the aℓm'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 Cℓ's, are m-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 m) 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 anyN>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 〈a2m1a3m2*〉, 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(36)P(α)={1,α∈[0,1],0,otherwise.
Let us now define two other variables x=cos(2πα) and y=sin(2πα). From these definitions, it follows that x and y are statistically dependent variables, since knowledge of the mean/variance of x automatically gives the mean/variance of y. However, these variables are clearly uncorrelated(37)〈xy〉=12π∫02πcosηsinηdη=0.
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.


## 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-sightn̂ would reflect simply the density and gravitational potentials (the “sources”) at the position Rn̂, where R 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ηR, the linear kernels of (1) reduce to Ki(x⃗′,η′;n̂)→βiδ(η′-ηR)δ(x⃗′-n̂R), where R=η0-ηR and βi are constant coefficients. The photon distribution that we measure on Earth would therefore be given by(4)Θ(n̂)≈∑i‍βiSi(x⃗′=n̂R,η′=ηR).
We can also express this result in terms of the Fourier spectrum of the sources as follows:(5)Θ(n̂)≈∑i‍βi∫d3k(2π)3eik⃗·n̂RSi(k⃗,ηR).
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:(6)eik⃗·x⃗=4π∑ℓmiℓjℓ(kx)Yℓm*(k̂)Yℓm(x̂),
where jℓ(z) are the spherical Bessel functions. Substituting (6) into (5) we obtain that(7)aℓm=∫d2n̂Yℓm*(n̂)Θ(n̂)≈∫d3k(2π)3Θ(k⃗)×4πiℓjℓ(kR)Yℓm*(k̂),
where we have loosely collected the sources into the term Θ(k⃗)≡∑iβiSi(k⃗,ηR). 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 aℓm are inherited from those of the Fourier modes Θ(k⃗) of the underlying matter and metric fields. Notice that the properties of the aℓms 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 Yℓm.


## 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)C(k⃗,k′⃗)≡〈Θ(k⃗)Θ(k′⃗)〉.
In the absence of any symmetries, this would be a generic function of the arguments k⃗ and k′⃗, with only two constraints: first, because Θ(x⃗) is a real function, Θ*(k⃗)=Θ(-k⃗), hence, in our definition C*(k⃗,k′⃗)=C(-k⃗,-k′⃗); second, due to the associative nature of the expectation value, C(k⃗,k′⃗)=C(k′⃗,k⃗). It is obvious how to generalize this definition to 3, 4, or an arbitrary number of fields at different k⃗s (or “points”).Let us first discuss the issue of gaussianity. If we say that the variablesΘ(k⃗) are Gaussian random numbers, then all the information that characterizes their distribution is contained in their two-point function C(k⃗,k′⃗). The probability distribution function (pdf) is then formally given by (9)P[Θ(k⃗),Θ(k′⃗)]~exp[-Θ(k⃗)Θ(k′⃗)2C(k⃗,k′⃗)].
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:  (10)〈Θ(k⃗1)Θ(k⃗2)⋯Θ(k⃗2N)〉G=∑i,j∏α=1NBi,jα〈Θ(k⃗i)Θ(k⃗j)〉,
where the sum runs over all permutations of the pairs of wave vectors and Bi,j 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 theN-point functions in real space, we should have the following  (11)〈Θ(x⃗1)Θ(x⃗2)⋯Θ(x⃗N)〉→Homog.CN(x⃗1-x⃗2,…,x⃗N-1-x⃗N).
Writing this expression in terms of the Fourier modes, we get the following(12)〈Θ(x⃗1)Θ(x⃗2)⋯Θ(x⃗N)〉=∫d3k1d3k1⋯d3kN(2π)3Ne-ik⃗1·x⃗1e-ik⃗2·x⃗2⋯e-ik⃗N·x⃗N×〈Θ(k⃗1)Θ(k⃗2)⋯Θ(k⃗N)〉.
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 δ(k⃗1+k⃗2+⋯+k⃗N). In other words, the hypothesis of homogeneity constrains the N-point function in Fourier space to be of the following form:  (13)〈Θ(k⃗1)Θ(k⃗2)⋯Θ(k⃗N)〉H=(2π)3Ñ(k⃗1,k⃗2,…,k⃗N)×δ(k⃗1+k⃗2+⋯+k⃗N).
Notice that the “(N-1)-spectrum” in Fourier space, Ñ, can still be a function of the directions of the wavenumbers k⃗i (it will be, in fact, a function of N-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 theN-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 k⃗i in (12). Let us instead rewrite the δ-functions in terms of triangles, so for the 4-point function we have  (14)δ(k⃗1+k⃗2+k⃗3+k⃗4)=∫d3qδ(k⃗1+k⃗2-q⃗)δ(k⃗3+k⃗4+q⃗),
whereas for the 5-point function we have  (15)δ(k⃗1+k⃗2+k⃗3+k⃗4+k⃗5)=∫d3qd3q′δ(k⃗1+k⃗2-q⃗)δ(q⃗+k⃗3-q⃗′)×δ(q⃗′+k⃗4+k⃗5),
and so on, so that the N-point δ-function is reduced to N-2 triangles with N-3 “internal momenta” (the idea is nicely illustrated in Figure 1.) Substituting the expression for the N-point δ-function into (12) and integrating out all external momenta but the first (k⃗1) and last (k⃗N), the result is as follows(16)〈Θ(x⃗1)Θ(x⃗2)⋯Θ(x⃗N)〉=1(2π)3N∫d3k1d3q1⋯d3qN-3d3kN×eik⃗1·(x⃗1-x⃗2)eiq⃗1·(x⃗2-x⃗3)⋯eiq⃗N-3·(x⃗N-2-x⃗N-1)×eik⃗N·(x⃗N-1-x⃗N)〈Θ(k⃗1)Θ(q⃗1-k⃗1)⋯Θ(k⃗N)〉.
This expression shows explicitly that the real-space N-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 theN-point functions that come from isotropy alone? Clearly, no dependence on the directions defined by the points, x⃗i-x⃗j, can arise in the final expression for the N-point functions in real space, so from (12) we see that the N-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 theN-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,Θ(n̂), which are taken over the unit sphere S2. We will discuss mostly the two-point function here, and we defer a fuller discussion of N-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  (17)〈aℓmaℓ′m′*〉=∫d3kd3k′(2π)6(4π)2iℓ(-i)ℓ′jℓ(kR)jℓ′(k′R)Yℓm(k̂)×Yℓ′m′*(k̂′)〈Θ(k⃗)Θ*(k⃗′)〉.
Under the hypothesis of homogeneity, this expression simplifies considerably, leading to  (18)〈aℓmaℓ′m′*〉H=∫d3k2πiℓ(-i)ℓ′jℓ(kR)jℓ′(kR)Yℓm(k̂)Yℓ′m′*(k̂)×Ñ2(k⃗).
If, in addition to homogeneity, we also assume isotropy, then Ñ2→P(k), and the integration over angles factors out, leading to the orthogonality condition for spherical harmonics as follows (19)∫d2k̂Yℓm(k̂)Yℓ′m′*(k̂)=δℓℓ′δmm′,
and as a result the covariance of the aℓms becomes diagonal as follows(20)〈aℓmaℓ′m′*〉H,I=δℓℓ′δmm′∫dkkjℓ2(kR)2πk3P(k)=4πδℓℓ′δmm′∫dlogkjℓ2(kR)ΔT2(k)≡Cℓδℓℓ′δmm′,
where we have defined the usual temperature power spectrumΔT(k)=k3P(k)/2π2 in the middle line, and the angular power spectrum Cℓ 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 logk as follows:  (21)〈Θ(x⃗)Θ(x⃗′)〉H,I=∫dlogksin(k|x⃗-x⃗′|)k|x⃗-x⃗′|ΔT2(k).In an analogous manner to what was done above, we can also construct theangular two-point correlation function in harmonic space as follows:(22)〈Θ(n̂)Θ(n̂′)〉=∑ℓm∑ℓ′m′〈aℓmaℓ′m′*〉Yℓm(n̂)Yℓm*(n̂′).
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:(23)〈Θ(n̂)Θ(n̂′)〉H,I=∑ℓm‍∑ℓ′m′Cℓδℓℓ′δmm′Yℓm(n̂)Yℓm*(n̂′)=∑ℓCℓ2ℓ+14πPℓ(n̂·n̂′).
Clearly, not only is this expression the analogous in S2 of (21), but in fact the Fourier power spectrum ΔT2(k) and the angular power spectrum Cℓ are defined in terms of each other as indicated in (23) as follows:  (24)Cℓ=4π∫dlogkjℓ2(kR)ΔT2(k).
Now, using the facts that the spherical Bessel function of order ℓ peaks when its argument is approximately given by ℓ, and that ∫dlogzjℓ2(z)=1/(2ℓ(ℓ+1)), we obtain the following (this is one type of  what has become known in the literature as Limber's approximations)(25)Cℓ≈2πℓ(ℓ+1)ΔT2(k=ℓR).
Incidentally, from this expression it is clear why it is customary to define (26)𝒞ℓ≡ℓ(ℓ+1)2πCℓ≈ΔT2(k=ℓR).Using (12), we can easily generalize the results of this subsection  to N-point functions in S2 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 N-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)C(n̂1,n̂2)=C(n̂2,n̂1).
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, (ℓ1,m1;ℓ2,m2), we could equally well have used the basis of total angular momentum (L,M;ℓ1,ℓ2) and decomposed that expression as  (28)C(n̂1,n̂2)=∑L,M∑ℓ1,ℓ2𝒜ℓ1ℓ2LM𝒴ℓ1ℓ2LM(n̂1,n̂2),
where 𝒴ℓ1ℓ2LM are known as the bipolar spherical harmonics, defined by [41](29)𝒴ℓ1ℓ2LM(n̂1,n̂2)=[Yℓ1(n̂1)⊗Yℓ2(n̂2)]LM,
where L and M=m1+m2 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(30)𝒜ℓ1ℓ2LM=∑m1m2〈aℓ1m1aℓ2m2〉(-1)M+ℓ1-ℓ22L+1(ℓ1ℓ2Lm1m2-M),
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 aℓms, 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 aℓms. 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 aℓm as statistically independent Gaussian random variables, and define a statistical estimator for their variances as the average(31)Ĉℓ≡12ℓ+1∑m=-ℓℓ|aℓm|2.
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 have(32)limℓ→∞Ĉℓ=Cℓ.
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 Cℓ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 spectrumCℓ 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 aℓms are Gaussian, then we can rewrite expression (31) as(33)(2ℓ+1)CℓĈℓ=Xℓ,Xℓ=∑m=-ℓℓ|aℓm|2Cℓ,
where Xℓ is a chi-square random variable with 2ℓ+1 degrees of freedom. According to the central limit theorem, when ℓ→∞, Xℓ approaches a standard normal variable (A standard normal variable is a Gaussian variable X with zero mean and unit variance. Any other Gaussian variable Y with mean μ and variance σ can be obtained from X through Y=σX+μ.) which implies that Ĉℓ will itself follow a Gaussian distribution. Its mean can be easily calculated using (20) and (31) and is of course given by(34)〈Ĉℓ〉=Cℓ,
which shows that the Ĉℓs are unbiased estimators of the Cℓs. It is also straightforward to calculate its variance (valid for any ℓ)(35)〈(Ĉℓ-Cℓ)(Ĉℓ′-Cℓ′)〉=22ℓ+1Cℓ2δℓℓ′.
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. Ĉℓ is therefore the best estimator we can build to measure the statistical properties of the multipolar coefficients aℓm when both statistical isotropy and gaussianity hold.In later Sections, we will explore angular or harmonicN-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 aℓms. The estimator for the angular power spectrum, Ĉℓ, takes into account all the aℓms by dividing them into the different ℓs and summing over all m∈(-ℓ,ℓ). Clearly, it will inherit a sample variance for small ℓs, when the aℓms can only be divided into a few “independent parts”. As we try to estimate higher-order objects such as the N-point functions, we will have to subdivide the aℓm'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 Cℓ's, are m-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 m) 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 anyN>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 〈a2m1a3m2*〉, 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(36)P(α)={1,α∈[0,1],0,otherwise.
Let us now define two other variables x=cos(2πα) and y=sin(2πα). From these definitions, it follows that x and y are statistically dependent variables, since knowledge of the mean/variance of x automatically gives the mean/variance of y. However, these variables are clearly uncorrelated(37)〈xy〉=12π∫02πcosηsinηdη=0.
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 Section3.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 ofN-Point Correlation Functions

We turn now to the question of non-Gaussian but statistically isotropic probabilities distributions. We will keep working with theN-point correlation function defined in harmonic space,(38)〈aℓ1m1aℓ2m2…aℓNmN〉,
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 N-point correlation function which is invariant under arbitrary 3-dimensional spatial rotations. When rotated to a new (primed) coordinate system, the N-point correlation function transforms as(39)〈aℓ1m1′aℓ2m2′…aℓNmN′〉=∑allm〈aℓ1m1aℓ2m2…aℓNmN〉Dm1′m1ℓ1Dm2′m2ℓ2⋯DmN′mNℓn,
where the Dmimi′ℓ(α,β,γ)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 ms. For the 2-point correlation function, we have already seen that the well-known expression(40)〈aℓ1m1aℓ2m2〉=(-1)m2Cℓ1δℓ1ℓ2δm1,-m2
does the job (41)〈aℓ1m1′aℓ2m2′〉=Cℓ1(∑m1(-1)m1Dm1′m1ℓ1Dm2′-m1ℓ1)δℓ1ℓ2=(-1)m2′Cℓ1δℓ1ℓ2δm1′,-m2′.
Note the importance of the angular spectrum, Cℓ, being an m-independent function.What about the 3-point function? In this case, the invariant combination is found to be(42)〈aℓ1m1aℓ2m2aℓ3m3〉=Bℓ1ℓ2ℓ3(ℓ1ℓ2ℓ3m1m2m3)
which can be verified by straightforward calculations. Again, the nontrivial physical content of this statistical moment is contained in an arbitrary but otherwise m-independent function: the bispectrumBℓ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 Bℓ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 theN-point correlation function. For the 4-point correlation function, for example, Hu [48] has found the following rotationally invariant combination(43)〈aℓ1m1…aℓ4m4〉=∑LMQℓ3ℓ4ℓ1ℓ2(L)(-1)M(ℓ1ℓ2Lm1m2-M)(ℓ3ℓ4Lm3m4M),
where the Qℓ3ℓ4ℓ1ℓ2(L) function is known as the trispectrum, and L is an internal angular momentum needed to ensure parity invariance. In a likewise manner, it can be verified that the following expression (44)〈aℓ1m1⋯aℓ5m5〉=∑LM‍∑L′M′Pℓ3ℓ4ℓ5ℓ1ℓ2(L,L′)(-1)M+M′(ℓ1ℓ2Lm1m2-M)×(ℓ3ℓ4L′m3m4-M′)(ℓ5LL′m5MM′)
gives the rotationally invariant quadrispectrum Pℓ3ℓ4ℓ5ℓ1ℓ2(L,L′).The examples above should be enough to show how the general structure of these functions emerges under SI; apart from anm-independent function, every pair of momenta ℓi 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 (ℓ,m).Figure 3
Geometrical representation of the 2pcf in terms of two unit vectors.Although we have always shownN-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 Ns. Here we will present a recipe for doing that, which also guarantees the uniqueness of the solution.The general recipe for obtaining the rotationally invariantN-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 (45)Dm1′m1ℓ1(ω)Dm2′m2ℓ2(ω)=∑L,M,M′(ℓ1ℓ2Lm1′m2′-M)(ℓ1ℓ2Lm1m2-M′)×(2L+1)(-1)M+M′DM′ML(ω),
and where ω={α,β,γ} is a shortcut notation for the three Euler angles. Once this contraction is done, there will remain N/2D-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:(46)〈aℓ1m1⋯aℓNmN〉=∑allm′〈aℓ1m1′⋯aℓNmN′〉×∑geometricalfactors×DMM′L(ω).
Now, we see that the only way for this combination to be rotationally invariant is when the remaining DMM′L function above does not depend on ω, that is, DMM′L(ω)=δL0δM0δM′0. Once this identity is applied to the geometrical factors, we are done, and the remaining terms inside the primed m-summation will give the rotationally invariant (N-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(47)〈aℓ1m1aℓ2m2〉=[∑m1′〈|aℓ1m1′|2〉2ℓ1+1](-1)m2δℓ1ℓ2δm1,-m2,
where, of course(48)Cℓ≡12ℓ+1∑m〈|aℓm|2〉
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(49)〈aℓ1m1aℓ2m2aℓ3m3〉=[∑m1′,m2′,m3′〈aℓ1m1′aℓ2m2′aℓ3m3′〉(ℓ1ℓ2ℓ3m1′m2′m3′)](ℓ1ℓ2ℓ3m1m2m3).
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:(50)Bℓ1ℓ2ℓ3≡∑m1,m2,m3〈aℓ1m1aℓ2m2aℓ3m3〉(ℓ1ℓ2ℓ3m1m2m3).It should be mentioned that this recipe not only enables us to establish the rotational invariance of anyN-point correlation function, but it also furnishes a straightforward definition of unbiased estimators for the N-point functions. All we have to do is to drop the ensemble average of the primed aℓm's. So, for example, for the 2- and 3-point functions above, the unbiased estimators are given, respectively, by (51)Ĉℓ=12ℓ+1∑maℓmaℓm*,B̂ℓ1ℓ2ℓ3=∑m1,m2,m3aℓ1m1aℓ2m2aℓ3m3(ℓ1ℓ2ℓ3m1m2m3).Notice that isotropy plays the same role, inS2, that homogeneity plays in ℝ3. What enforces homogeneity in ℝ3 is the Fourier-space δ-functions, as in the discussion around Figure 1. However, in S2, 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 N-point angular correlation function. Hence, the diagrammatic representations of the constituents of the N-point functions in Fourier (Figure 1) and in harmonic space (Figure 2) really do convey the same physical idea—one in ℝ3, the other in S2.


### 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 anyN-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)P(k⃗)=∑ℓ,mrℓm(k)Yℓm(k̂),
where the reality of P(k⃗) requires that rℓm(k)=(-1)mrℓ,-m*(k). Given that temperature perturbations Θ(x⃗) are real, their Fourier components must satisfy the relation Θ(k⃗)=Θ*(-k⃗). This property taken together with the definition (13) implies that(53)P(k⃗)=P(-k⃗),
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 z→-z, we find that (-1)ℓ-m=1. Similarly, invariance under the transformations x→-x and y→-y implies the conditions rℓm=(-1)mrℓ,-m and rℓm=rℓ,-m, respectively. Gathering all these constraints with the parity of ℓ, we conclude that(54)rℓm∈ℝ,ℓ,m∈2ℕ.
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(55)〈aℓ1m1aℓ2m2*〉=∑ℓ3,m3𝒢m1m2m3ℓ1ℓ2ℓ3Hℓ1ℓ2ℓ3m3,
where(56)𝒢m1m2m3ℓ1ℓ2ℓ3=(-1)m1(2ℓ1+1)(2ℓ2+1)(2ℓ3+1)4π×(ℓ1ℓ2ℓ3000)(ℓ1ℓ2ℓ3-m1m2m3)
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:(57)ℓ1+ℓ2+ℓ3∈2ℕ,m1+m2+m3=0,|ℓi-ℓj|≤ℓk≤ℓi+ℓj,∀i,j,k∈{1,2,3}.
The remaining coefficients in (55) are given by(58)Hℓ1ℓ2ℓ3m3=4πiℓ1-ℓ2∫0∞dlogkrℓ3m3(k)jℓ1(kR)jℓ2(kR)
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 ℓ3 is even, the quantity ℓ1±ℓ2 must also be even, that is, multipoles with different parity do not couple to each other in this class of models(59)〈aℓ1m1aℓ2m2*〉=0,ℓ1+ℓ2=odd.
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 Cℓmδℓℓ′δmm′. If the Cℓ,0 happens to be zero, for example, then all multipoles will present a preferred direction (in this case, the z-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 theCℓ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(60)〈aℓmaℓm*〉=Cℓ+ϵ∑ℓ′>0𝒢m,-m,0ℓℓℓ′Hℓℓℓ′0,
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〈aℓmaℓm*〉. Since we are looking for the diagonal entries of the matrix (55), a first guess would be(61)Ĉℓeff=12ℓ+1∑m=-ℓℓaℓmaℓm*.
To check whether this is an unbiased estimator of the effective angular spectrum, we apply it to (55) and take its average(62)〈Cℓeff〉=Cℓ+ϵ∑ℓ′>0,m𝒢m,-m,0ℓℓℓ′Hℓℓℓ′0.
Using the definition (56) of the Gaunt coefficients, the m-summation in the expression above becomes [41](63)∑m(-1)ℓ-m(ℓℓℓ′m-m0)=2ℓ+1δℓ′,0=0,
where the last equality follows because ℓ'>0. Consequently, we conclude that(64)〈Ĉℓeff〉=Cℓ≠〈aℓmaℓm*〉.
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,(65)〈(Ĉℓeff-Cℓ)(Ĉℓ′eff-Cℓ′)〉=22ℓ+1Cℓ2δℓℓ′+𝒪(ϵ2).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.


## 3.1. Non-Gaussian and SI Models

### 3.1.1. Rotational Invariance ofN-Point Correlation Functions

We turn now to the question of non-Gaussian but statistically isotropic probabilities distributions. We will keep working with theN-point correlation function defined in harmonic space,(38)〈aℓ1m1aℓ2m2…aℓNmN〉,
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 N-point correlation function which is invariant under arbitrary 3-dimensional spatial rotations. When rotated to a new (primed) coordinate system, the N-point correlation function transforms as(39)〈aℓ1m1′aℓ2m2′…aℓNmN′〉=∑allm〈aℓ1m1aℓ2m2…aℓNmN〉Dm1′m1ℓ1Dm2′m2ℓ2⋯DmN′mNℓn,
where the Dmimi′ℓ(α,β,γ)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 ms. For the 2-point correlation function, we have already seen that the well-known expression(40)〈aℓ1m1aℓ2m2〉=(-1)m2Cℓ1δℓ1ℓ2δm1,-m2
does the job (41)〈aℓ1m1′aℓ2m2′〉=Cℓ1(∑m1(-1)m1Dm1′m1ℓ1Dm2′-m1ℓ1)δℓ1ℓ2=(-1)m2′Cℓ1δℓ1ℓ2δm1′,-m2′.
Note the importance of the angular spectrum, Cℓ, being an m-independent function.What about the 3-point function? In this case, the invariant combination is found to be(42)〈aℓ1m1aℓ2m2aℓ3m3〉=Bℓ1ℓ2ℓ3(ℓ1ℓ2ℓ3m1m2m3)
which can be verified by straightforward calculations. Again, the nontrivial physical content of this statistical moment is contained in an arbitrary but otherwise m-independent function: the bispectrumBℓ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 Bℓ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 theN-point correlation function. For the 4-point correlation function, for example, Hu [48] has found the following rotationally invariant combination(43)〈aℓ1m1…aℓ4m4〉=∑LMQℓ3ℓ4ℓ1ℓ2(L)(-1)M(ℓ1ℓ2Lm1m2-M)(ℓ3ℓ4Lm3m4M),
where the Qℓ3ℓ4ℓ1ℓ2(L) function is known as the trispectrum, and L is an internal angular momentum needed to ensure parity invariance. In a likewise manner, it can be verified that the following expression (44)〈aℓ1m1⋯aℓ5m5〉=∑LM‍∑L′M′Pℓ3ℓ4ℓ5ℓ1ℓ2(L,L′)(-1)M+M′(ℓ1ℓ2Lm1m2-M)×(ℓ3ℓ4L′m3m4-M′)(ℓ5LL′m5MM′)
gives the rotationally invariant quadrispectrum Pℓ3ℓ4ℓ5ℓ1ℓ2(L,L′).The examples above should be enough to show how the general structure of these functions emerges under SI; apart from anm-independent function, every pair of momenta ℓi 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 (ℓ,m).Figure 3
Geometrical representation of the 2pcf in terms of two unit vectors.Although we have always shownN-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 Ns. Here we will present a recipe for doing that, which also guarantees the uniqueness of the solution.The general recipe for obtaining the rotationally invariantN-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 (45)Dm1′m1ℓ1(ω)Dm2′m2ℓ2(ω)=∑L,M,M′(ℓ1ℓ2Lm1′m2′-M)(ℓ1ℓ2Lm1m2-M′)×(2L+1)(-1)M+M′DM′ML(ω),
and where ω={α,β,γ} is a shortcut notation for the three Euler angles. Once this contraction is done, there will remain N/2D-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:(46)〈aℓ1m1⋯aℓNmN〉=∑allm′〈aℓ1m1′⋯aℓNmN′〉×∑geometricalfactors×DMM′L(ω).
Now, we see that the only way for this combination to be rotationally invariant is when the remaining DMM′L function above does not depend on ω, that is, DMM′L(ω)=δL0δM0δM′0. Once this identity is applied to the geometrical factors, we are done, and the remaining terms inside the primed m-summation will give the rotationally invariant (N-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(47)〈aℓ1m1aℓ2m2〉=[∑m1′〈|aℓ1m1′|2〉2ℓ1+1](-1)m2δℓ1ℓ2δm1,-m2,
where, of course(48)Cℓ≡12ℓ+1∑m〈|aℓm|2〉
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(49)〈aℓ1m1aℓ2m2aℓ3m3〉=[∑m1′,m2′,m3′〈aℓ1m1′aℓ2m2′aℓ3m3′〉(ℓ1ℓ2ℓ3m1′m2′m3′)](ℓ1ℓ2ℓ3m1m2m3).
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:(50)Bℓ1ℓ2ℓ3≡∑m1,m2,m3〈aℓ1m1aℓ2m2aℓ3m3〉(ℓ1ℓ2ℓ3m1m2m3).It should be mentioned that this recipe not only enables us to establish the rotational invariance of anyN-point correlation function, but it also furnishes a straightforward definition of unbiased estimators for the N-point functions. All we have to do is to drop the ensemble average of the primed aℓm's. So, for example, for the 2- and 3-point functions above, the unbiased estimators are given, respectively, by (51)Ĉℓ=12ℓ+1∑maℓmaℓm*,B̂ℓ1ℓ2ℓ3=∑m1,m2,m3aℓ1m1aℓ2m2aℓ3m3(ℓ1ℓ2ℓ3m1m2m3).Notice that isotropy plays the same role, inS2, that homogeneity plays in ℝ3. What enforces homogeneity in ℝ3 is the Fourier-space δ-functions, as in the discussion around Figure 1. However, in S2, 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 N-point angular correlation function. Hence, the diagrammatic representations of the constituents of the N-point functions in Fourier (Figure 1) and in harmonic space (Figure 2) really do convey the same physical idea—one in ℝ3, the other in S2.


## 3.1.1. Rotational Invariance ofN-Point Correlation Functions

We turn now to the question of non-Gaussian but statistically isotropic probabilities distributions. We will keep working with theN-point correlation function defined in harmonic space,(38)〈aℓ1m1aℓ2m2…aℓNmN〉,
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 N-point correlation function which is invariant under arbitrary 3-dimensional spatial rotations. When rotated to a new (primed) coordinate system, the N-point correlation function transforms as(39)〈aℓ1m1′aℓ2m2′…aℓNmN′〉=∑allm〈aℓ1m1aℓ2m2…aℓNmN〉Dm1′m1ℓ1Dm2′m2ℓ2⋯DmN′mNℓn,
where the Dmimi′ℓ(α,β,γ)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 ms. For the 2-point correlation function, we have already seen that the well-known expression(40)〈aℓ1m1aℓ2m2〉=(-1)m2Cℓ1δℓ1ℓ2δm1,-m2
does the job (41)〈aℓ1m1′aℓ2m2′〉=Cℓ1(∑m1(-1)m1Dm1′m1ℓ1Dm2′-m1ℓ1)δℓ1ℓ2=(-1)m2′Cℓ1δℓ1ℓ2δm1′,-m2′.
Note the importance of the angular spectrum, Cℓ, being an m-independent function.What about the 3-point function? In this case, the invariant combination is found to be(42)〈aℓ1m1aℓ2m2aℓ3m3〉=Bℓ1ℓ2ℓ3(ℓ1ℓ2ℓ3m1m2m3)
which can be verified by straightforward calculations. Again, the nontrivial physical content of this statistical moment is contained in an arbitrary but otherwise m-independent function: the bispectrumBℓ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 Bℓ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 theN-point correlation function. For the 4-point correlation function, for example, Hu [48] has found the following rotationally invariant combination(43)〈aℓ1m1…aℓ4m4〉=∑LMQℓ3ℓ4ℓ1ℓ2(L)(-1)M(ℓ1ℓ2Lm1m2-M)(ℓ3ℓ4Lm3m4M),
where the Qℓ3ℓ4ℓ1ℓ2(L) function is known as the trispectrum, and L is an internal angular momentum needed to ensure parity invariance. In a likewise manner, it can be verified that the following expression (44)〈aℓ1m1⋯aℓ5m5〉=∑LM‍∑L′M′Pℓ3ℓ4ℓ5ℓ1ℓ2(L,L′)(-1)M+M′(ℓ1ℓ2Lm1m2-M)×(ℓ3ℓ4L′m3m4-M′)(ℓ5LL′m5MM′)
gives the rotationally invariant quadrispectrum Pℓ3ℓ4ℓ5ℓ1ℓ2(L,L′).The examples above should be enough to show how the general structure of these functions emerges under SI; apart from anm-independent function, every pair of momenta ℓi 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 (ℓ,m).Figure 3
Geometrical representation of the 2pcf in terms of two unit vectors.Although we have always shownN-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 Ns. Here we will present a recipe for doing that, which also guarantees the uniqueness of the solution.The general recipe for obtaining the rotationally invariantN-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 (45)Dm1′m1ℓ1(ω)Dm2′m2ℓ2(ω)=∑L,M,M′(ℓ1ℓ2Lm1′m2′-M)(ℓ1ℓ2Lm1m2-M′)×(2L+1)(-1)M+M′DM′ML(ω),
and where ω={α,β,γ} is a shortcut notation for the three Euler angles. Once this contraction is done, there will remain N/2D-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:(46)〈aℓ1m1⋯aℓNmN〉=∑allm′〈aℓ1m1′⋯aℓNmN′〉×∑geometricalfactors×DMM′L(ω).
Now, we see that the only way for this combination to be rotationally invariant is when the remaining DMM′L function above does not depend on ω, that is, DMM′L(ω)=δL0δM0δM′0. Once this identity is applied to the geometrical factors, we are done, and the remaining terms inside the primed m-summation will give the rotationally invariant (N-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(47)〈aℓ1m1aℓ2m2〉=[∑m1′〈|aℓ1m1′|2〉2ℓ1+1](-1)m2δℓ1ℓ2δm1,-m2,
where, of course(48)Cℓ≡12ℓ+1∑m〈|aℓm|2〉
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(49)〈aℓ1m1aℓ2m2aℓ3m3〉=[∑m1′,m2′,m3′〈aℓ1m1′aℓ2m2′aℓ3m3′〉(ℓ1ℓ2ℓ3m1′m2′m3′)](ℓ1ℓ2ℓ3m1m2m3).
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:(50)Bℓ1ℓ2ℓ3≡∑m1,m2,m3〈aℓ1m1aℓ2m2aℓ3m3〉(ℓ1ℓ2ℓ3m1m2m3).It should be mentioned that this recipe not only enables us to establish the rotational invariance of anyN-point correlation function, but it also furnishes a straightforward definition of unbiased estimators for the N-point functions. All we have to do is to drop the ensemble average of the primed aℓm's. So, for example, for the 2- and 3-point functions above, the unbiased estimators are given, respectively, by (51)Ĉℓ=12ℓ+1∑maℓmaℓm*,B̂ℓ1ℓ2ℓ3=∑m1,m2,m3aℓ1m1aℓ2m2aℓ3m3(ℓ1ℓ2ℓ3m1m2m3).Notice that isotropy plays the same role, inS2, that homogeneity plays in ℝ3. What enforces homogeneity in ℝ3 is the Fourier-space δ-functions, as in the discussion around Figure 1. However, in S2, 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 N-point angular correlation function. Hence, the diagrammatic representations of the constituents of the N-point functions in Fourier (Figure 1) and in harmonic space (Figure 2) really do convey the same physical idea—one in ℝ3, the other in S2.


## 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 anyN-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)P(k⃗)=∑ℓ,mrℓm(k)Yℓm(k̂),
where the reality of P(k⃗) requires that rℓm(k)=(-1)mrℓ,-m*(k). Given that temperature perturbations Θ(x⃗) are real, their Fourier components must satisfy the relation Θ(k⃗)=Θ*(-k⃗). This property taken together with the definition (13) implies that(53)P(k⃗)=P(-k⃗),
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 z→-z, we find that (-1)ℓ-m=1. Similarly, invariance under the transformations x→-x and y→-y implies the conditions rℓm=(-1)mrℓ,-m and rℓm=rℓ,-m, respectively. Gathering all these constraints with the parity of ℓ, we conclude that(54)rℓm∈ℝ,ℓ,m∈2ℕ.
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(55)〈aℓ1m1aℓ2m2*〉=∑ℓ3,m3𝒢m1m2m3ℓ1ℓ2ℓ3Hℓ1ℓ2ℓ3m3,
where(56)𝒢m1m2m3ℓ1ℓ2ℓ3=(-1)m1(2ℓ1+1)(2ℓ2+1)(2ℓ3+1)4π×(ℓ1ℓ2ℓ3000)(ℓ1ℓ2ℓ3-m1m2m3)
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:(57)ℓ1+ℓ2+ℓ3∈2ℕ,m1+m2+m3=0,|ℓi-ℓj|≤ℓk≤ℓi+ℓj,∀i,j,k∈{1,2,3}.
The remaining coefficients in (55) are given by(58)Hℓ1ℓ2ℓ3m3=4πiℓ1-ℓ2∫0∞dlogkrℓ3m3(k)jℓ1(kR)jℓ2(kR)
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 ℓ3 is even, the quantity ℓ1±ℓ2 must also be even, that is, multipoles with different parity do not couple to each other in this class of models(59)〈aℓ1m1aℓ2m2*〉=0,ℓ1+ℓ2=odd.
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 Cℓmδℓℓ′δmm′. If the Cℓ,0 happens to be zero, for example, then all multipoles will present a preferred direction (in this case, the z-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 theCℓ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(60)〈aℓmaℓm*〉=Cℓ+ϵ∑ℓ′>0𝒢m,-m,0ℓℓℓ′Hℓℓℓ′0,
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〈aℓmaℓm*〉. Since we are looking for the diagonal entries of the matrix (55), a first guess would be(61)Ĉℓeff=12ℓ+1∑m=-ℓℓaℓmaℓm*.
To check whether this is an unbiased estimator of the effective angular spectrum, we apply it to (55) and take its average(62)〈Cℓeff〉=Cℓ+ϵ∑ℓ′>0,m𝒢m,-m,0ℓℓℓ′Hℓℓℓ′0.
Using the definition (56) of the Gaunt coefficients, the m-summation in the expression above becomes [41](63)∑m(-1)ℓ-m(ℓℓℓ′m-m0)=2ℓ+1δℓ′,0=0,
where the last equality follows because ℓ'>0. Consequently, we conclude that(64)〈Ĉℓeff〉=Cℓ≠〈aℓmaℓm*〉.
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,(65)〈(Ĉℓeff-Cℓ)(Ĉℓ′eff-Cℓ′)〉=22ℓ+1Cℓ2δℓℓ′+𝒪(ϵ2).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.


## 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)P(k⃗)=∑ℓ,mrℓm(k)Yℓm(k̂),
where the reality of P(k⃗) requires that rℓm(k)=(-1)mrℓ,-m*(k). Given that temperature perturbations Θ(x⃗) are real, their Fourier components must satisfy the relation Θ(k⃗)=Θ*(-k⃗). This property taken together with the definition (13) implies that(53)P(k⃗)=P(-k⃗),
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 z→-z, we find that (-1)ℓ-m=1. Similarly, invariance under the transformations x→-x and y→-y implies the conditions rℓm=(-1)mrℓ,-m and rℓm=rℓ,-m, respectively. Gathering all these constraints with the parity of ℓ, we conclude that(54)rℓm∈ℝ,ℓ,m∈2ℕ.
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(55)〈aℓ1m1aℓ2m2*〉=∑ℓ3,m3𝒢m1m2m3ℓ1ℓ2ℓ3Hℓ1ℓ2ℓ3m3,
where(56)𝒢m1m2m3ℓ1ℓ2ℓ3=(-1)m1(2ℓ1+1)(2ℓ2+1)(2ℓ3+1)4π×(ℓ1ℓ2ℓ3000)(ℓ1ℓ2ℓ3-m1m2m3)
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:(57)ℓ1+ℓ2+ℓ3∈2ℕ,m1+m2+m3=0,|ℓi-ℓj|≤ℓk≤ℓi+ℓj,∀i,j,k∈{1,2,3}.
The remaining coefficients in (55) are given by(58)Hℓ1ℓ2ℓ3m3=4πiℓ1-ℓ2∫0∞dlogkrℓ3m3(k)jℓ1(kR)jℓ2(kR)
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 ℓ3 is even, the quantity ℓ1±ℓ2 must also be even, that is, multipoles with different parity do not couple to each other in this class of models(59)〈aℓ1m1aℓ2m2*〉=0,ℓ1+ℓ2=odd.
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 Cℓmδℓℓ′δmm′. If the Cℓ,0 happens to be zero, for example, then all multipoles will present a preferred direction (in this case, the z-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 theCℓ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(60)〈aℓmaℓm*〉=Cℓ+ϵ∑ℓ′>0𝒢m,-m,0ℓℓℓ′Hℓℓℓ′0,
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〈aℓmaℓm*〉. Since we are looking for the diagonal entries of the matrix (55), a first guess would be(61)Ĉℓeff=12ℓ+1∑m=-ℓℓaℓmaℓm*.
To check whether this is an unbiased estimator of the effective angular spectrum, we apply it to (55) and take its average(62)〈Cℓeff〉=Cℓ+ϵ∑ℓ′>0,m𝒢m,-m,0ℓℓℓ′Hℓℓℓ′0.
Using the definition (56) of the Gaunt coefficients, the m-summation in the expression above becomes [41](63)∑m(-1)ℓ-m(ℓℓℓ′m-m0)=2ℓ+1δℓ′,0=0,
where the last equality follows because ℓ'>0. Consequently, we conclude that(64)〈Ĉℓeff〉=Cℓ≠〈aℓmaℓm*〉.
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,(65)〈(Ĉℓeff-Cℓ)(Ĉℓ′eff-Cℓ′)〉=22ℓ+1Cℓ2δℓℓ′+𝒪(ϵ2).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 (E- and B-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 thebottom-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 thetop-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) functionh that satisfies Laplace's equation(66)∇2h=0,
where ∇2 is the Laplace operator. In spherical coordinates, the formal solution to Laplace's equation which is regular at the origin (r=0) is(67)h=∑ℓ=0∞hℓ(r,θ,φ),hℓ=∑m=-ℓℓaℓmrℓYℓm(θ,φ).
The functions rℓYℓm 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 ℓ:(68)h=∑ℓ=0∞hℓ(x,y,z),hℓ=∑abcℓλabcxaybzc,(a+b+c=ℓ).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 formv⃗·∇≡∇v⃗ over the monopole potential 1/r, where r=x2+y2+z2 and v⃗ is a unit vector, he could construct solutions of the form (68). That is,(69)fℓ(x,y,z)=λℓ∇v⃗ℓ⋯∇v⃗2∇v⃗11r|r=1,
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(70)f0=1r,f1=∇v⃗1f0=-v⃗1·r⃗r3,f2=∇v⃗2f1=3(v⃗2·r⃗)(v⃗1·r⃗)-r2(v⃗1·v⃗2)r3,
and so on. By induction, one can show that the general expression will be given by(71)fℓ=(-1)ℓ(2ℓ-1)!!∏i=1ℓv⃗i·r⃗+r2Qℓ-2r2ℓ+1,
where Qℓ-2 is a homogeneous polynomial of order ℓ-2 which only involves combinations of the vectors v⃗i and r⃗.The numerator of the functionfℓ given by (71) is clearly a homogeneous polynomial of order ℓ (as one can easily check for some ℓs and also prove by mathematical induction). A not so obvious result is that this polynomial is also harmonic. To prove that, let us define(72)gℓ=(-1)ℓ(2ℓ-1)!!∏i=1ℓv⃗i·r⃗+r2Qℓ-2
and consider the application of the operator ∇2 over the combination rαgℓ. For the x-component we get(73)∂x2(rαgℓ)=rα∂x2gℓ+2αrα-2x∂xgℓ+[αrα-2+α(α-2)x2rα-4]gℓ.
Repeating this process for the y- and z-components and then adding the results, we find (74)∇2(rαgℓ)=rα∇2gℓ+2αrα-2(x∂xgℓ+y∂ygℓ+z∂zgℓ)+α(α+1)rα-2gℓ=rα∇2gℓ+α(α+2ℓ+1)(rα-2gℓ),
where in the last step we have used Euler's theorem on homogeneous functions, that is, r⃗·∇gℓ=ℓgℓ. If we now choose α=-(2ℓ+1), we find immediately that(75)∇2(gℓr2ℓ+1)=∇2gℓr2ℓ+1.By construction, the left-hand side of the above expression is equal to∇2fℓ. But according to the definition (69), this quantity is also zero, since Laplace's operator commutes with directional derivatives and ∇2(1/r)=0 for r>0. Therefore, ∇2gℓ=0 and gℓ 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 aℓm's) and ℓ unit vectors v⃗i. 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 definedirections, 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 v⃗ℓ,p (p=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)w⃗ℓ,pp′≡v⃗ℓ,p⊗v⃗ℓ,p′.
Mathematically, the sign ambiguity implies that the multipole vectors belong to the quotient space S2/ℤ2 (also known as ℝP2) while the normal vectors belong to ℝ3/ℤ2.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. TheR Statistic

The first example of a test involving the multipole vectors would be a scalar product, such asv⃗·v⃗′. The most natural test would involve asking whether the ℓ multipole vectors of the given multipole ℓ are especially aligned or not. This means computing(77)Rℓℓ=2ℓ(ℓ-1)∑p,p′>pℓ|v⃗ℓ,p·v⃗ℓ,p′|,
where the normalization was introduced to make 0≤Rℓℓ≤1.This idea could be generalized to test alignments between multipole vectors at different multipoles:(78)Rℓℓ′=1ℓℓ′∑p,p′|v⃗ℓ,p·v⃗ℓ′,p′|.
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 R23 is of the order of 90%–95% CL [35].


#### 4.2.2. TheS 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)w⃗ℓ,pp′=v⃗ℓ,p⊗v⃗ℓ′,p′.
Notice that the number of normal vectors for a given multipole ℓ is l=ℓ(ℓ-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(80)Sℓℓ=2l(l-1)∑p,p′>pl|w⃗ℓ,p·w⃗ℓ,p′|,
where the normalization was introduced again to make 0≤Sℓℓ≤1.And yet again this idea is easily generalized to test alignments between normal vectors at different multipoles:(81)Sℓℓ′=1ll′∑p,p′|w⃗ℓ,p·w⃗ℓ′,p′|.
With this test, the quadrupole-octupole alignment is much more significant: for essentially all CMB maps the test S23 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:(82)Dℓℓ′=1ll′∑p,p′|ŵℓ,p·ŵℓ′,p′|.
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:(83)Bℓℓ′=1ℓl′∑p,p′|v⃗ℓ,p·w⃗ℓ′,p′|.
Of course, this test cannot be easily generalized to a self-alignment. With the B23 test, the quadrupole-octupole alignment is significant to about 95%–98% CLWe could go on here, but it should be clear that all the information (the2ℓ+1 real degrees of freedom for each multipole ℓ) has already been exhausted in the tests above.


## 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) functionh that satisfies Laplace's equation(66)∇2h=0,
where ∇2 is the Laplace operator. In spherical coordinates, the formal solution to Laplace's equation which is regular at the origin (r=0) is(67)h=∑ℓ=0∞hℓ(r,θ,φ),hℓ=∑m=-ℓℓaℓmrℓYℓm(θ,φ).
The functions rℓYℓm 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 ℓ:(68)h=∑ℓ=0∞hℓ(x,y,z),hℓ=∑abcℓλabcxaybzc,(a+b+c=ℓ).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 formv⃗·∇≡∇v⃗ over the monopole potential 1/r, where r=x2+y2+z2 and v⃗ is a unit vector, he could construct solutions of the form (68). That is,(69)fℓ(x,y,z)=λℓ∇v⃗ℓ⋯∇v⃗2∇v⃗11r|r=1,
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(70)f0=1r,f1=∇v⃗1f0=-v⃗1·r⃗r3,f2=∇v⃗2f1=3(v⃗2·r⃗)(v⃗1·r⃗)-r2(v⃗1·v⃗2)r3,
and so on. By induction, one can show that the general expression will be given by(71)fℓ=(-1)ℓ(2ℓ-1)!!∏i=1ℓv⃗i·r⃗+r2Qℓ-2r2ℓ+1,
where Qℓ-2 is a homogeneous polynomial of order ℓ-2 which only involves combinations of the vectors v⃗i and r⃗.The numerator of the functionfℓ given by (71) is clearly a homogeneous polynomial of order ℓ (as one can easily check for some ℓs and also prove by mathematical induction). A not so obvious result is that this polynomial is also harmonic. To prove that, let us define(72)gℓ=(-1)ℓ(2ℓ-1)!!∏i=1ℓv⃗i·r⃗+r2Qℓ-2
and consider the application of the operator ∇2 over the combination rαgℓ. For the x-component we get(73)∂x2(rαgℓ)=rα∂x2gℓ+2αrα-2x∂xgℓ+[αrα-2+α(α-2)x2rα-4]gℓ.
Repeating this process for the y- and z-components and then adding the results, we find (74)∇2(rαgℓ)=rα∇2gℓ+2αrα-2(x∂xgℓ+y∂ygℓ+z∂zgℓ)+α(α+1)rα-2gℓ=rα∇2gℓ+α(α+2ℓ+1)(rα-2gℓ),
where in the last step we have used Euler's theorem on homogeneous functions, that is, r⃗·∇gℓ=ℓgℓ. If we now choose α=-(2ℓ+1), we find immediately that(75)∇2(gℓr2ℓ+1)=∇2gℓr2ℓ+1.By construction, the left-hand side of the above expression is equal to∇2fℓ. But according to the definition (69), this quantity is also zero, since Laplace's operator commutes with directional derivatives and ∇2(1/r)=0 for r>0. Therefore, ∇2gℓ=0 and gℓ 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 aℓm's) and ℓ unit vectors v⃗i. 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 definedirections, 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 v⃗ℓ,p (p=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)w⃗ℓ,pp′≡v⃗ℓ,p⊗v⃗ℓ,p′.
Mathematically, the sign ambiguity implies that the multipole vectors belong to the quotient space S2/ℤ2 (also known as ℝP2) while the normal vectors belong to ℝ3/ℤ2.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. TheR Statistic

The first example of a test involving the multipole vectors would be a scalar product, such asv⃗·v⃗′. The most natural test would involve asking whether the ℓ multipole vectors of the given multipole ℓ are especially aligned or not. This means computing(77)Rℓℓ=2ℓ(ℓ-1)∑p,p′>pℓ|v⃗ℓ,p·v⃗ℓ,p′|,
where the normalization was introduced to make 0≤Rℓℓ≤1.This idea could be generalized to test alignments between multipole vectors at different multipoles:(78)Rℓℓ′=1ℓℓ′∑p,p′|v⃗ℓ,p·v⃗ℓ′,p′|.
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 R23 is of the order of 90%–95% CL [35].


### 4.2.2. TheS 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)w⃗ℓ,pp′=v⃗ℓ,p⊗v⃗ℓ′,p′.
Notice that the number of normal vectors for a given multipole ℓ is l=ℓ(ℓ-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(80)Sℓℓ=2l(l-1)∑p,p′>pl|w⃗ℓ,p·w⃗ℓ,p′|,
where the normalization was introduced again to make 0≤Sℓℓ≤1.And yet again this idea is easily generalized to test alignments between normal vectors at different multipoles:(81)Sℓℓ′=1ll′∑p,p′|w⃗ℓ,p·w⃗ℓ′,p′|.
With this test, the quadrupole-octupole alignment is much more significant: for essentially all CMB maps the test S23 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:(82)Dℓℓ′=1ll′∑p,p′|ŵℓ,p·ŵℓ′,p′|.
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:(83)Bℓℓ′=1ℓl′∑p,p′|v⃗ℓ,p·w⃗ℓ′,p′|.
Of course, this test cannot be easily generalized to a self-alignment. With the B23 test, the quadrupole-octupole alignment is significant to about 95%–98% CLWe could go on here, but it should be clear that all the information (the2ℓ+1 real degrees of freedom for each multipole ℓ) has already been exhausted in the tests above.


## 4.2.1. TheR Statistic

The first example of a test involving the multipole vectors would be a scalar product, such asv⃗·v⃗′. The most natural test would involve asking whether the ℓ multipole vectors of the given multipole ℓ are especially aligned or not. This means computing(77)Rℓℓ=2ℓ(ℓ-1)∑p,p′>pℓ|v⃗ℓ,p·v⃗ℓ,p′|,
where the normalization was introduced to make 0≤Rℓℓ≤1.This idea could be generalized to test alignments between multipole vectors at different multipoles:(78)Rℓℓ′=1ℓℓ′∑p,p′|v⃗ℓ,p·v⃗ℓ′,p′|.
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 R23 is of the order of 90%–95% CL [35].


## 4.2.2. TheS 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)w⃗ℓ,pp′=v⃗ℓ,p⊗v⃗ℓ′,p′.
Notice that the number of normal vectors for a given multipole ℓ is l=ℓ(ℓ-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(80)Sℓℓ=2l(l-1)∑p,p′>pl|w⃗ℓ,p·w⃗ℓ,p′|,
where the normalization was introduced again to make 0≤Sℓℓ≤1.And yet again this idea is easily generalized to test alignments between normal vectors at different multipoles:(81)Sℓℓ′=1ll′∑p,p′|w⃗ℓ,p·w⃗ℓ′,p′|.
With this test, the quadrupole-octupole alignment is much more significant: for essentially all CMB maps the test S23 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:(82)Dℓℓ′=1ll′∑p,p′|ŵℓ,p·ŵℓ′,p′|.
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:(83)Bℓℓ′=1ℓl′∑p,p′|v⃗ℓ,p·w⃗ℓ′,p′|.
Of course, this test cannot be easily generalized to a self-alignment. With the B23 test, the quadrupole-octupole alignment is significant to about 95%–98% CLWe could go on here, but it should be clear that all the information (the2ℓ+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 functionC of the form(84)C:S2×S2→ℝ.
We have seen in Section 2.3 that if we choose spherical coordinates (θi,φi) to describe each vector n̂i, the function above can be decomposed either in terms of two spherical harmonics Yℓimi or in terms of the bipolar spherical harmonic 𝒴ℓ1ℓ2LM. In any case, the 2pcf will have the following functional dependence:(85)C=C(θ1,φ1,θ2,φ2).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 Figure5).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 n̂1 and n̂2 to be the same and expand a function of the form(86)C=C(θ1,φ1),
which is equivalent to C : S2→ℝ. 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 n̂1=n̂2, 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 n̂1 and n̂2 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 variablex. The shaded area is the probability (98).In this representation,θ is the angle between the vectors n̂1 and n̂2, as usual. The normal to the plane, n̂, 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)C=C(Θ,Φ,θ,φ).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):(88)C→12π∫02πC(Θ,Φ,θ,φ)dφ.
The resulting function can be easily expand in terms of simple special functions as(89)C(Θ,Φ,θ)=∑ℓ∑l,m2ℓ+14π𝒞ℓlmPℓ(cosθ)Ylm(Θ,Φ),l∈2ℕ,
where the restriction on the l-mode results from the symmetry n̂1↔n̂2.(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 𝒞ℓlm correspond to a generalization of the usual angular power spectrum Cℓ's. In fact, they can be seen as the coefficients of a spherical harmonic decomposition of the function Cℓ(n̂), 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𝒞ℓlm 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 𝒞ℓlm and the temperature coefficients aℓm 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 𝒞ℓlm's in terms of the aℓm's. But for that to work we need the relation between the set of angles (Θ,Φ,θ) and (θ1,φ1,θ2,φ2) which, depending on the reference frame we choose, is extremely complicated. Fortunately, all we need is the following relation:(90)n̂1·n̂2=cosθ=cosθ1cosθ2+sinθ1sinθ2cos(φ1-φ2),
together with a suitable choice of our coordinate system. For example, we can use the invariance of the scalar product n̂1·n̂2 and choose our coordinate system such that the discs of Figure 4 lies in the xy plane. With this choice we will have(91)(Θ,Φ)=(0,0),cosθ=cos(φ1-φ2),
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(92)(-1)m𝒞ℓlm2l+1=2π∑ℓ1,m1∑ℓ2,m2〈aℓ1m2aℓ2m2〉(lℓ1ℓ2-mm1m2)Iℓ1ℓ2l,ℓ,
where(93)Iℓ1ℓ2l,ℓ≡∑m(-1)mλℓ1mλℓ2m(lℓ1ℓ20m-m)∫0πd(-cosθ)Pℓ(cosθ)eimθ,
and where the λℓim's form a set of coefficients resulting from the θ integration, which are zero unless ℓi+m= even (see the Appendix for more details).Expression (92) is what we were looking for. With this relation, the angular-planar power spectrum 𝒞ℓlm 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 𝒞ℓlm to the bipolar power spectrum 𝒜ℓ1ℓ2LM 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:(94)𝒞ℓlm=2π∑ℓ1,ℓ2Aℓ1ℓ2lm(-1)ℓ1-ℓ2(2l+1)Iℓ1ℓ2l,ℓ.
Here, the set of geometrical coefficients Iℓ1ℓ2l,ℓ 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χ2 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 spectrumCℓ (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(95)𝒞̂ℓlm≡2π2l+1∑ℓ1,m1∑ℓ2,m2aℓ1m1aℓ2m2(lℓ1ℓ2mm1m2)Iℓ1ℓ2l,ℓ,
for in this case we have an unbiased estimator:(96)〈𝒞̂ℓlm〉=𝒞ℓlm,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 (χ2) goodness-of-fit test, which in our case can be written in the following generalized form:(97)(χν2)ℓl≡12l+1∑m=-ll|𝒞̂ℓlm-𝒞ℓlm|2(σℓlm)2,
in which ℓ and l are the angular and planar degrees of freedom, respectively, and where σℓlm is just the standard deviation of the estimator 𝒞̂ℓlm. The (2l+1)-1 factor accounts for the 2l+1 planar degrees of freedom and was introduced for latter convenience.In Section7 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.


## 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):(88)C→12π∫02πC(Θ,Φ,θ,φ)dφ.
The resulting function can be easily expand in terms of simple special functions as(89)C(Θ,Φ,θ)=∑ℓ∑l,m2ℓ+14π𝒞ℓlmPℓ(cosθ)Ylm(Θ,Φ),l∈2ℕ,
where the restriction on the l-mode results from the symmetry n̂1↔n̂2.(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 𝒞ℓlm correspond to a generalization of the usual angular power spectrum Cℓ's. In fact, they can be seen as the coefficients of a spherical harmonic decomposition of the function Cℓ(n̂), 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𝒞ℓlm 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 𝒞ℓlm and the temperature coefficients aℓm 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 𝒞ℓlm's in terms of the aℓm's. But for that to work we need the relation between the set of angles (Θ,Φ,θ) and (θ1,φ1,θ2,φ2) which, depending on the reference frame we choose, is extremely complicated. Fortunately, all we need is the following relation:(90)n̂1·n̂2=cosθ=cosθ1cosθ2+sinθ1sinθ2cos(φ1-φ2),
together with a suitable choice of our coordinate system. For example, we can use the invariance of the scalar product n̂1·n̂2 and choose our coordinate system such that the discs of Figure 4 lies in the xy plane. With this choice we will have(91)(Θ,Φ)=(0,0),cosθ=cos(φ1-φ2),
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(92)(-1)m𝒞ℓlm2l+1=2π∑ℓ1,m1∑ℓ2,m2〈aℓ1m2aℓ2m2〉(lℓ1ℓ2-mm1m2)Iℓ1ℓ2l,ℓ,
where(93)Iℓ1ℓ2l,ℓ≡∑m(-1)mλℓ1mλℓ2m(lℓ1ℓ20m-m)∫0πd(-cosθ)Pℓ(cosθ)eimθ,
and where the λℓim's form a set of coefficients resulting from the θ integration, which are zero unless ℓi+m= even (see the Appendix for more details).Expression (92) is what we were looking for. With this relation, the angular-planar power spectrum 𝒞ℓlm 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 𝒞ℓlm to the bipolar power spectrum 𝒜ℓ1ℓ2LM 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:(94)𝒞ℓlm=2π∑ℓ1,ℓ2Aℓ1ℓ2lm(-1)ℓ1-ℓ2(2l+1)Iℓ1ℓ2l,ℓ.
Here, the set of geometrical coefficients Iℓ1ℓ2l,ℓ 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χ2 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 spectrumCℓ (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(95)𝒞̂ℓlm≡2π2l+1∑ℓ1,m1∑ℓ2,m2aℓ1m1aℓ2m2(lℓ1ℓ2mm1m2)Iℓ1ℓ2l,ℓ,
for in this case we have an unbiased estimator:(96)〈𝒞̂ℓlm〉=𝒞ℓlm,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 (χ2) goodness-of-fit test, which in our case can be written in the following generalized form:(97)(χν2)ℓl≡12l+1∑m=-ll|𝒞̂ℓlm-𝒞ℓlm|2(σℓlm)2,
in which ℓ and l are the angular and planar degrees of freedom, respectively, and where σℓlm is just the standard deviation of the estimator 𝒞̂ℓlm. The (2l+1)-1 factor accounts for the 2l+1 planar degrees of freedom and was introduced for latter convenience.In Section7 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.


## 5.1.1. Angular-Planar Power Spectrum

Since we are restricting our analysis to the Gaussian framework, the set of coefficients𝒞ℓlm 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 𝒞ℓlm and the temperature coefficients aℓm 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 𝒞ℓlm's in terms of the aℓm's. But for that to work we need the relation between the set of angles (Θ,Φ,θ) and (θ1,φ1,θ2,φ2) which, depending on the reference frame we choose, is extremely complicated. Fortunately, all we need is the following relation:(90)n̂1·n̂2=cosθ=cosθ1cosθ2+sinθ1sinθ2cos(φ1-φ2),
together with a suitable choice of our coordinate system. For example, we can use the invariance of the scalar product n̂1·n̂2 and choose our coordinate system such that the discs of Figure 4 lies in the xy plane. With this choice we will have(91)(Θ,Φ)=(0,0),cosθ=cos(φ1-φ2),
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(92)(-1)m𝒞ℓlm2l+1=2π∑ℓ1,m1∑ℓ2,m2〈aℓ1m2aℓ2m2〉(lℓ1ℓ2-mm1m2)Iℓ1ℓ2l,ℓ,
where(93)Iℓ1ℓ2l,ℓ≡∑m(-1)mλℓ1mλℓ2m(lℓ1ℓ20m-m)∫0πd(-cosθ)Pℓ(cosθ)eimθ,
and where the λℓim's form a set of coefficients resulting from the θ integration, which are zero unless ℓi+m= even (see the Appendix for more details).Expression (92) is what we were looking for. With this relation, the angular-planar power spectrum 𝒞ℓlm 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 𝒞ℓlm to the bipolar power spectrum 𝒜ℓ1ℓ2LM 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:(94)𝒞ℓlm=2π∑ℓ1,ℓ2Aℓ1ℓ2lm(-1)ℓ1-ℓ2(2l+1)Iℓ1ℓ2l,ℓ.
Here, the set of geometrical coefficients Iℓ1ℓ2l,ℓ 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χ2 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 spectrumCℓ (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(95)𝒞̂ℓlm≡2π2l+1∑ℓ1,m1∑ℓ2,m2aℓ1m1aℓ2m2(lℓ1ℓ2mm1m2)Iℓ1ℓ2l,ℓ,
for in this case we have an unbiased estimator:(96)〈𝒞̂ℓlm〉=𝒞ℓlm,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 (χ2) goodness-of-fit test, which in our case can be written in the following generalized form:(97)(χν2)ℓl≡12l+1∑m=-ll|𝒞̂ℓlm-𝒞ℓlm|2(σℓlm)2,
in which ℓ and l are the angular and planar degrees of freedom, respectively, and where σℓlm is just the standard deviation of the estimator 𝒞̂ℓlm. The (2l+1)-1 factor accounts for the 2l+1 planar degrees of freedom and was introduced for latter convenience.In Section7 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

Supposex is a random variable predicted by a particular model and that cosmological observations of this quantity returned the value x0. (Rigorously, x is only one realization of a random variable X, which is a real-valued function defined on a sample space. By the same reason we should not call the aℓm'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 x≤x0. Assuming that 𝒫th is the (normalized) probability density function (pdf) of x, this probability is commonly defined to be(98)P≤≡∫-∞x0𝒫th(x)dx.
The probability of having x>x0 is then simply given by P>=1-P≤. See Figure 5.IfP≤ is found to be too small or, equivalently, too high, we might be tempted to interpret x0 as “anomalous” according to this model. However, this definition of probability assumes that x0 was measured with infinite precision, and so it says nothing about an important question we must deal with: typically, the measurement of x0 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 valuex0, our measurements should also be regarded as random events. In a more rigorous approach, we would have to consider x0 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 x0 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 x0. 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 ofx being smaller than our measurements when the latter are also random events. Let us suppose that our measurement is described by the random variable y and that x0 is its most probable value. The probability of having x≤y is simply the probability of having(99)z≡x-y≤0
for some particular realization of the variable y. It should be clear by now that if we know the pdf of z, the probability we are looking for is simply the area under this distribution for -∞≤z≤0. The probability P≤ of z being smaller or equal to zero can be calculated as(100)P≤≡P{(x,y)∣x-y≤0}=∬x-y≤0𝒫(x,y)dxdy.
Now, under the hypothesis of independence of x and y we have 𝒫(x,y)=𝒫th(x)𝒫obs(y), where 𝒫obs is the (normalized) probability distribution function of the variable y. We can therefore rewrite the last expression as(101)P≤=∫-∞∞[∫-∞y𝒫th(x)dx]𝒫obs(y)dy.
If we now hold y fixed and do the change of variables x=u+y, we get(102)P≤=∫-∞0[∫-∞∞𝒫th(u+y)𝒫obs(y)dy]du,
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 y, we can equally well hold x fixed and repeat the calculus, obtaining the symmetric version of this result. In fact, the final pdf for z is nothing else than the convolution of the pdfs of each variable [76]:(103)𝒫(z)≡(𝒫obs*𝒫th)(z)=∫-∞∞𝒫obs(y)𝒫th(z±y)dy=∫-∞∞𝒫obs(x∓z)𝒫th(x)dx,
where the plus (minus) sign refers to the difference (sum) of x and y. Integrating this pdf from (-∞,0] we get our answer(104)P≤≡∫-∞0𝒫(z)dz.
As a consistency check, notice that in the limit where observations are made with infinite precision, 𝒫obs(y) becomes a delta function and we have(105)P≤=∫-∞0∫-∞+∞δ(y-x0)𝒫th(z±y)dydz=∫-∞±x0𝒫th(x)dx
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 x0 to the most probable (not the mean!) value x 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 calculateP≤ 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 x being smaller than the observed x0; all pdf's were normalized to one.Figure 6
Probability density functions forx (blue, dashed line), y (red, dot-dashed line), and z (solid line). The shaded area gives the probability of x being smaller than x0 (dashed vertical line). See the text for more details.


## 6.1. Standard Calculations

Supposex is a random variable predicted by a particular model and that cosmological observations of this quantity returned the value x0. (Rigorously, x is only one realization of a random variable X, which is a real-valued function defined on a sample space. By the same reason we should not call the aℓm'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 x≤x0. Assuming that 𝒫th is the (normalized) probability density function (pdf) of x, this probability is commonly defined to be(98)P≤≡∫-∞x0𝒫th(x)dx.
The probability of having x>x0 is then simply given by P>=1-P≤. See Figure 5.IfP≤ is found to be too small or, equivalently, too high, we might be tempted to interpret x0 as “anomalous” according to this model. However, this definition of probability assumes that x0 was measured with infinite precision, and so it says nothing about an important question we must deal with: typically, the measurement of x0 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 valuex0, our measurements should also be regarded as random events. In a more rigorous approach, we would have to consider x0 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 x0 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 x0. 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 ofx being smaller than our measurements when the latter are also random events. Let us suppose that our measurement is described by the random variable y and that x0 is its most probable value. The probability of having x≤y is simply the probability of having(99)z≡x-y≤0
for some particular realization of the variable y. It should be clear by now that if we know the pdf of z, the probability we are looking for is simply the area under this distribution for -∞≤z≤0. The probability P≤ of z being smaller or equal to zero can be calculated as(100)P≤≡P{(x,y)∣x-y≤0}=∬x-y≤0𝒫(x,y)dxdy.
Now, under the hypothesis of independence of x and y we have 𝒫(x,y)=𝒫th(x)𝒫obs(y), where 𝒫obs is the (normalized) probability distribution function of the variable y. We can therefore rewrite the last expression as(101)P≤=∫-∞∞[∫-∞y𝒫th(x)dx]𝒫obs(y)dy.
If we now hold y fixed and do the change of variables x=u+y, we get(102)P≤=∫-∞0[∫-∞∞𝒫th(u+y)𝒫obs(y)dy]du,
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 y, we can equally well hold x fixed and repeat the calculus, obtaining the symmetric version of this result. In fact, the final pdf for z is nothing else than the convolution of the pdfs of each variable [76]:(103)𝒫(z)≡(𝒫obs*𝒫th)(z)=∫-∞∞𝒫obs(y)𝒫th(z±y)dy=∫-∞∞𝒫obs(x∓z)𝒫th(x)dx,
where the plus (minus) sign refers to the difference (sum) of x and y. Integrating this pdf from (-∞,0] we get our answer(104)P≤≡∫-∞0𝒫(z)dz.
As a consistency check, notice that in the limit where observations are made with infinite precision, 𝒫obs(y) becomes a delta function and we have(105)P≤=∫-∞0∫-∞+∞δ(y-x0)𝒫th(z±y)dydz=∫-∞±x0𝒫th(x)dx
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 x0 to the most probable (not the mean!) value x 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 calculateP≤ 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 x being smaller than the observed x0; all pdf's were normalized to one.Figure 6
Probability density functions forx (blue, dashed line), y (red, dot-dashed line), and z (solid line). The shaded area gives the probability of x being smaller than x0 (dashed vertical line). See the text for more details.


## 7. Aχ2 Test of Statistical Isotropy

Although the last example was constructed to emphasize an important feature of the formalism developed in Section6, 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 〈aℓ1m1aℓ2m2〉=(-1)m2Cℓ1δℓ1ℓ2δm1,-m2. Using this expression in (92), we find that(106)𝒞ℓlm=(SI)Cℓδl0δm0.
On the other hand, if the only nonzero 𝒞ℓlm's are given by l=m=0, then 𝒞ℓ00=Cℓ. 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 l≠0, that is,(107)𝒞ℓlm=0,(l≥2),
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(108)〈(𝒞̂ℓlm)*𝒞̂ℓl′m′〉=8π2∑ℓ1,ℓ2Cℓ1Cℓ2(Iℓ1ℓ2l,ℓ)2δll′δmm′.
This matrix has some interesting properties. First, note that the planar degrees of freedom are independent in this case (which justifies the (2l+1)-1) factor introduced in (97)). Second, its diagonal elements are given by the variance (σℓlm)2=〈(𝒞̂ℓlm)*𝒞̂ℓlm〉, which now becomes m-independent:(109)(σℓlm)2→(σℓl)2=8π2∑ℓ1,ℓ2Cℓ1Cℓ2(Iℓ1ℓ2l,ℓ)2.Therefore, for the particular case of theΛCDM model, the chi-square test (97) gets even simpler. Using (107) and (109), we find(110)(χν2)ℓl=12l+1∑m=-ll|𝒞̂ℓlm|2(σℓl)2.
It is now clear that if the data under analysis are really described by this model, then it must be true that(111)〈(χν2)ℓl〉=12l+1∑m=-ll〈(𝒞̂ℓlm)*𝒞̂ℓlm〉(σℓl)2=1,
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(112)χ¯ℓl≡(χν2)ℓl-1
which will quantify anisotropies whenever χ¯ℓl 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σℓl using Cℓmodel), which would only include our a priori prejudices about what the map's anisotropies should look like. Since the Cℓ's are, by construction, a measure of statistical isotropy, consequently, an “anomalous” detection of Cℓ'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:(113)x≡(χ¯2)ℓ(th)l,y≡(χ¯2)ℓ(obs)l,
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 x defined above we have run 2×104 Monte Carlo simulations of Gaussian and statistically isotropic CMB maps using the ΛCDM best-fit Cℓ's provided by the WMAP team [77]. With these maps we have then constructed 2×104 realizations of the variable x.The simulation of the (observational) variabley 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 variabley (see the text for more details). Note that the reference [78] includes the analysis of maps from the three- and five-year WMAP releases.
Full sky mapsReferencesHinshaw et al.[14, 79]de Oliveira-Costa et al.[80]Kim et al.[78]Park et al.[81]Delabrouille et al.[82]The prescription we adopt to determine the distribution of the observational variabley is as follows: we simulate Gaussian random aℓm'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 n different full-sky temperature maps at hand and we want to estimate the randomness inherent in the determination of, let's say, a32; therefore, we take(114)𝒩(a32ILC5,σ32)→a32,
with(115)σ32=1n-1∑i=1n(a32i-a¯32)2,a¯32=1n∑i=1na32i,
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 Table1 to construct 104 Gaussian random aℓm's, which were then used to calculate 104 realizations of y=(χ¯2)ℓ(obs)l. With those variables we constructed histograms which, together with the histograms for the (full-sky) variable x, were used to calculate the final probability (104). We have restricted our analysis to the range of values (ℓ,l)∈[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 (χ¯2)ℓ(obs)l, as given by full-sky temperature maps.
l∖ℓ2345678910281.1%73.7%54.1%6.1%80.7%46.4%36.9%47.5%81.8%474.0%72.6%55.0%39.6%74.2%93.2%51.6%55.4%56.3%678.1%80.7%69.3%52.3%33.6%80.0%95.0%50.3%82.2%863.5%87.7%18.8%51.5%21.4%66.4%31.6%27.5%82.3%1067.9%50.0%61.0%8.7%37.7%59.5%36.6%29.2%35.7%Histograms for(χ¯2)ℓ(th)l (blue), (χ¯2)ℓ(obs)l (purple) and for the difference (χ¯2)ℓ(th)l-(χ¯2)ℓ(obs)l (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.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Overall, our results show no significant planar deviations of anisotropy in WMAP data. The most unlikely individual values in Table2 are in the sectors (l,ℓ) 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 (χ¯2)ℓl are relevant, but also their coherence over a range of angular or planar momenta carries interesting information. So, for example, a set of (χ¯2)ℓl'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)2l, (χ¯2)3l and, to a lesser extent, (χ¯2)ℓ4 (see Table 2). The angular quadrupole ℓ=2, as well as the angular octupole ℓ=3, have all positive planar spectra (for all values of l which we were able to compute), indicated by probabilities larger than 50%. The planar hexadecupole l=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].


## 7.1.ΛCDM Model

In the simplest realization of theΛCDM model, the covariance matrix of temperature maps is determined by 〈aℓ1m1aℓ2m2〉=(-1)m2Cℓ1δℓ1ℓ2δm1,-m2. Using this expression in (92), we find that(106)𝒞ℓlm=(SI)Cℓδl0δm0.
On the other hand, if the only nonzero 𝒞ℓlm's are given by l=m=0, then 𝒞ℓ00=Cℓ. 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 l≠0, that is,(107)𝒞ℓlm=0,(l≥2),
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(108)〈(𝒞̂ℓlm)*𝒞̂ℓl′m′〉=8π2∑ℓ1,ℓ2Cℓ1Cℓ2(Iℓ1ℓ2l,ℓ)2δll′δmm′.
This matrix has some interesting properties. First, note that the planar degrees of freedom are independent in this case (which justifies the (2l+1)-1) factor introduced in (97)). Second, its diagonal elements are given by the variance (σℓlm)2=〈(𝒞̂ℓlm)*𝒞̂ℓlm〉, which now becomes m-independent:(109)(σℓlm)2→(σℓl)2=8π2∑ℓ1,ℓ2Cℓ1Cℓ2(Iℓ1ℓ2l,ℓ)2.Therefore, for the particular case of theΛCDM model, the chi-square test (97) gets even simpler. Using (107) and (109), we find(110)(χν2)ℓl=12l+1∑m=-ll|𝒞̂ℓlm|2(σℓl)2.
It is now clear that if the data under analysis are really described by this model, then it must be true that(111)〈(χν2)ℓl〉=12l+1∑m=-ll〈(𝒞̂ℓlm)*𝒞̂ℓlm〉(σℓl)2=1,
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(112)χ¯ℓl≡(χν2)ℓl-1
which will quantify anisotropies whenever χ¯ℓl 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σℓl using Cℓmodel), which would only include our a priori prejudices about what the map's anisotropies should look like. Since the Cℓ's are, by construction, a measure of statistical isotropy, consequently, an “anomalous” detection of Cℓ'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:(113)x≡(χ¯2)ℓ(th)l,y≡(χ¯2)ℓ(obs)l,
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 x defined above we have run 2×104 Monte Carlo simulations of Gaussian and statistically isotropic CMB maps using the ΛCDM best-fit Cℓ's provided by the WMAP team [77]. With these maps we have then constructed 2×104 realizations of the variable x.The simulation of the (observational) variabley 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 variabley (see the text for more details). Note that the reference [78] includes the analysis of maps from the three- and five-year WMAP releases.
Full sky mapsReferencesHinshaw et al.[14, 79]de Oliveira-Costa et al.[80]Kim et al.[78]Park et al.[81]Delabrouille et al.[82]The prescription we adopt to determine the distribution of the observational variabley is as follows: we simulate Gaussian random aℓm'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 n different full-sky temperature maps at hand and we want to estimate the randomness inherent in the determination of, let's say, a32; therefore, we take(114)𝒩(a32ILC5,σ32)→a32,
with(115)σ32=1n-1∑i=1n(a32i-a¯32)2,a¯32=1n∑i=1na32i,
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 Table1 to construct 104 Gaussian random aℓm's, which were then used to calculate 104 realizations of y=(χ¯2)ℓ(obs)l. With those variables we constructed histograms which, together with the histograms for the (full-sky) variable x, were used to calculate the final probability (104). We have restricted our analysis to the range of values (ℓ,l)∈[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 (χ¯2)ℓ(obs)l, as given by full-sky temperature maps.
l∖ℓ2345678910281.1%73.7%54.1%6.1%80.7%46.4%36.9%47.5%81.8%474.0%72.6%55.0%39.6%74.2%93.2%51.6%55.4%56.3%678.1%80.7%69.3%52.3%33.6%80.0%95.0%50.3%82.2%863.5%87.7%18.8%51.5%21.4%66.4%31.6%27.5%82.3%1067.9%50.0%61.0%8.7%37.7%59.5%36.6%29.2%35.7%Histograms for(χ¯2)ℓ(th)l (blue), (χ¯2)ℓ(obs)l (purple) and for the difference (χ¯2)ℓ(th)l-(χ¯2)ℓ(obs)l (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.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Overall, our results show no significant planar deviations of anisotropy in WMAP data. The most unlikely individual values in Table2 are in the sectors (l,ℓ) 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 (χ¯2)ℓl are relevant, but also their coherence over a range of angular or planar momenta carries interesting information. So, for example, a set of (χ¯2)ℓl'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)2l, (χ¯2)3l and, to a lesser extent, (χ¯2)ℓ4 (see Table 2). The angular quadrupole ℓ=2, as well as the angular octupole ℓ=3, have all positive planar spectra (for all values of l which we were able to compute), indicated by probabilities larger than 50%. The planar hexadecupole l=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.



---

*Source: 378203-2010-08-24.xml*