Reconstruction of Missing Data in the Sky using Iterative Harmonic Expansion
Atsushi J. Nishizawa∗
Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU,
WPI), The University of Tokyo, Chiba 277-8583, Japan

Kaiki Taro Inoue
Department of Science and Engineering, Kinki University, Higashi-Osaka, 577-8502, Japan(Dated: May 1, 2013)
Missing fluctuations in masked regions in the sky can be reconstructed from fluctuations in the surrounding unmasked regions if they are sufficiently smooth. We propose to reconstruct such missing fluctuations by
iteratively applying a spherical harmonic expansion to fluctuations in the unmasked region. The accuracy ofreconstruction depends on the mask geometries, the spectrum of underlying density fluctuations, and the number of iterations. For Gaussian fluctuations with the Harrison-Zel’dovich spectrum, our method provides moreaccurate restoration than naive methods using the brute-forth matrix inversion or the singular value decomposition.
I.

INTRODUCTION

After the first data release of the cosmic microwave background (CMB) temperature fluctuations observed by the
Wilkinson Microwave Anisotropy Probe (WMAP), anomalous signatures in the CMB on large angular scales, so called“large-angle anomalies” have been claimed by many authors[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] and also confirmed recently by
the Planck mission [12]. So far, the origin has been veiled inmystery. They may be due to (a) difference in a priori significance and a posteriori significance [13, 14, 15, 16] (b)incomplete subtraction of foreground emissions [17, 18, 19](c) contribution from large scale structures via the integrated
Sachs-Wolfe effect [12, 20, 21, 22, 23, 24, 25, 26, 27, 28], or
kinetic Sunyaev-Zel’dovich effect [29], (d) possible systematics from instruments [30] (e) incomplete treatment of masking[31], or (f) extensions of inflationary models [32, 33, 34, 35,
36, 37, 38, 39].
In real observation, we cannot observe the whole sky witha sufficiently high signal-to-noise (S/N) ratio. A conservative
approach is to mask out the region (e.g., the Zone of Avoidance), where the S/N ratio is low and just ignore the data onit.
For estimating the power spectrum of fluctuations in thesky, we can use deconvolution techniques [e.g. 40] assuming
an isotropy prior on the power spectrum. However, for estimating the density field itself, we need to develop methods
that can reconstruct the phases (if expressed in complex numbers) as well as the amplitudes of missing fluctuations.
In order to reconstruct missing data on the masked region,
we need to find the inverse of the masking operator. However,in general, the masking matrix is singular thus not invertible.
We need to assume a certain prior information about the underlying data such as isotropy and smoothness [31, 41, 42, 43]
or their derivatives [44] in order to regularize the inverse operator. Since the result depends on the choice of the prior, therobustness of each reconstruction method should be mutually∗ Electronic address: atsushi.nishizawa@ipmu.jp

checked.
In this paper, we propose a new method for regularizing
the inverse of masking operators using the iterative sphericalharmonic expansion (IHE). As our aim is to probe the origin of
the large-angle CMB anomalies, we ignore the effect of noisein the signal in the following, which is a good approximation
for large-angle fluctuations.
This paper is organized as follows. In section II, we describe the formulation of our reconstruction method using theiterative spherical harmonic expansion and our simulations for
random Gaussian fluctuations. In section III, we show the accuracy of the reconstruction, which depends on the pixelization of the sky, the underlying power spectrum, and the number of iteration. And then we compare the accuracy of our
method to that of the brute-force and the singular value decomposition (SVD) method. In section IV, we give our conclusion.

II. ITERATIVE HARMONIC EXPANSION
A. Theory

Suppose a density field δ(γ̂) on a unit sphere where γ̂ represents a unit vector that specifies the angular position in thesky. The observed field is given byδobs (γ̂) = W(γ̂)δ(γ̂),(1)

where W = 0 in the masked region and W = 1 elsewhere.
Then we can obtain the pseudo harmonic coefficient ãi by simply integrating over the whole sky,
Z
(2)
ãi =
dΩ δ(γ̂)W(γ̂)Yi∗ (γ̂)
R
where Wi j = dΩ WY j Yi∗ is the (i, j) component of the mode
coupling matrix due to incomplete sky coverage, Yi (γ̂) is aspherical harmonic, and dΩ is the surface element on a unit
sphere. For simplicity, we use the single subscript index i =`2 +` +m+1 where ` and m represent a multipole moment and2(a) Original map with mask(c) Pseudo-a`m reconstruction(e) Pseudo-a`m residual(b) Original smoothed map without mask(d) IHE reconstruction(f) IHE residual

FIG. 1: Plots of an original and reconstructed maps. (a) Original simulated map with `cut = 30 and a sky cut |b| < 20◦ , which is used as an
input. (b) Original simulated map smoothed with `cut = 7 on a full sky (no mask). (c) Reconstructed map with `max = 7 using pseudo-ai .
(d) Reconstructed map with `max = 7 using IHE with Nite = 10 iterations. (e) Residual map after subtraction:(b)-(c). (f) Residual map after
subtraction:(d)-(b). Color scale (-1.87 to 1.55) is the same for all the figures.

an azimuthal number, respectively. In terms of true harmoniccoefficients to be reconstructed atrue
j , Eq. (2) can be written as
ãi =

∞
X

atrue
j Wi j .

(3)

j

The mode coupling matrix W is explicitly given as [40]Z
X
Wi1 i2 =w i3
dΩ Yi∗3 Yi2 Yi1
i3#1/2(2`1 + 1)(2`2 + 1)(2`3 + 1)
4π
i3
!
!
` ` `
`1 `2 `3
× 1 2 3,
(4)0 0 0
m1 −m2 m3
R
where wi =dΩ WYi∗ and (:::) denotes the Wigner 3 jsymbols. Note that the subscript i j depends only on ` j andm j , i.e. i j = `2j + ` j + m j + 1.
As shown in Eq. (3), atrue
i ’s can be obtained by inverting
the mode coupling matrix W. However, in general, it can be
singular and not invertible. In order to regularize the inversion, we propose the following iterative scheme. The iteration
process starts from
Z
ã(1)
=
dΩ δ(γ̂)W(γ̂)Yi∗ (γ̂),
(5)
i
=

X

"

wi3 (−1)m2

which we call “the pseudo-ai ’s”. For m > 1, the m-th set ofãi ’s is constructed iteratively from two maps; the original observed map at the unmasked region and the map reconstructedfrom the inverse transform of the (m − 1)-th ãi ’s at the masked
region,
Z
h
i
(m)
ãi =
dΩ δ(γ̂)W(γ̂) + δ̃(m−1) (γ̂)S (γ̂) Yi∗ (γ̂), (6)

where
X
∞




ãi(m−1) Yi (γ̂) m ≥ 2,

(m−1)
δ̃
(γ̂) = i0m = 1,(7)and S = 1 − W. In what follows, we assume that fluctuations whose angular scales are larger than that of the masked
region do not significantly correlate with fluctuations smaller
than the masked region. In that case, the summation in Eq. (3)and Eq. (7) can be truncated at a certain multipole imax as longas we concern large-angle fluctuations corresponding to multipoles i  imax . In what follows, we sum up ai ’s up to themultipole iimax , and we omit the summation symbol when no
confusion arises.
Substituting Eq. (7) into (6) recursively, we obtain the general formula of the m-th iterated harmonic coefficients as aseries of S i j ,K
m−1ãi(m) = ã(1)]i j ),j (δi j + S i j + S ik S k j + · · · + [S(8)where δiKj is the Kronecker delta, and S i j = δiKj − Wi j . The map

between m-th iterated harmonic coefficient ãi(m) is equivalent
to the product of ã(1)i ’s and the Taylor expansion of the inverseof the mode coupling matrix W truncated at the (m − 1)-thorder,(1)K−1−1ãi(m) ' ã(1)j [(δ − S ) ]i j = ã j [W ]i j ,(9)Thus ã(m)i ’s approximately represent fluctuations in which themode coupling matrix W is deconvolved. The number of iteration Nite should be evaluated using Monte-Carlo simulations,which is discussed later in Sec.III.

3
We can regard the above process as a mapping from a observable to the estimator,
FIHE : ãi → ã(m) .

(10)

As we will see in Sec. III, the accuracy of the IHE method
depends on the maximum multipole `max to be reconstructed,
the spectrum of the underlying field δ, the mask geometry,
and the number of iteration Nite . Thus we can write FIHE =FIHE (`max , δ, W, Nite ).
Note that the computation time for estimating ã(Nite ) can besignificantly reduced if we use Eq. (6) instead of Eq. (8). The
reason is as follows. As the rank of the matrix S i j is of the2
4order `max
, it costs Npix `maxby Eq. (4), and matrix algebra of
6
Eq. (8) costs ∼ Nite `max computations. If one uses Eq. (6), it2, where Nite is the number
would be of the order Nite Npix `maxof iteration, and Npix is the number of pixels which tile thesky, and `max is the maximum multipole to be reconstructed.
For example, for given `max = 10, Npix = 12 × 10242 pixels
and 5 time iteration, the computation time will be reduced bya factor of ∼ 25.

B. Simulation

In order to assess the accuracy of reconstruction, we generate 1000 random isotropic Gaussian density fields in thesky. We use the code synfast, which is publicly available asa package in Healpix [50] [45] for generating random Gaussian maps. First, we use the Harrison-Zel’dovich spectrum as
the input power spectrum. It gives an angular power spectrumC` ∝ ln where n = −2 on large angular scales, in the Einstein de-Sitter universe which correspond to the Sachs-Wolfe
plateau of the CMB power spectrum. As a simple model ofthe zone of avoidance, we consider an azimuthally symmetric mask W(γ̂) = 0 at the Galactic latitude b smaller than 20
degrees, i.e. |b| < 20◦ and W = 1 otherwise. We use pixels
that are sufficiently smaller than the size of the mask and thefluctuation scales to be reconstructed in order to reduce theerrors due to the pixelization effect (see Sec. III A). We adopt
the Healpix resolution Nside = 512 (the total number of pixels
2on the whole sky is Npix = 12 × Nside= 3145728). The input power spectrum should be truncated at sufficiently small
scales `cut compared to the scale to be reconstructed `max , i.e.
`max  `cut . The cut-off scale `cut is inferred from the detector resolution scale in CMB observations or the shot noise
dominant scales in galaxy surveys. In our simulated maps, we
set `cut = 30, which is sufficiently smaller than the claimed
anomalous multipoles in the CMB. We find that the choice of`cut is not sensitive to our result.
III. RESULTIn this section, we first explore the pixelization effect,
which may affect the accuracy of reconstruction in Sec. III A
and then we compare our reconstruction method with otherinversion methods using the brute-force or the singular value

FIG. 2: The accuracy of map reconstruction averaged over 1000
realizations of Gaussian maps as a function of the iteration number
for a mask |b| < 20◦ with `max = 7 for various pixel resolution Nside .
Both solid and dashed lines are obtained by using Eqs. (6) and (7).
Dashed lines represent values for which approximated S i j ’s for pixelbased masks are used. Solid lines show values for which the exact
expression of Eq. (11) is used. Dotted lines indicate values obtainedby using the direct inversion of the mode coupling matrix W, thesecond equation in Eq. (9). In the bottom panel, we show the relative
differences of the accuracy in percentage. Nside = 1024 provides
accuracy better than a few percent. Nite ∼ 10 with Nside = 512 canachieve sub-percent accuracy.

decomposition (SVD) in Sec. III B. Finally, we discuss thecondition of the original maps that can be reconstructed byusing the IHE in Sec.III C.
A. Pixelization effectWe need to assess the accuracy of the mode coupling matrix W since it may be close to being singular. In other words,
the difference between W(γ̂) and W(γ̂i ) should be carefully
checked where γ̂ and γ̂i are unpixelized and pixelized positions on a unit sphere, respectively. In order to do so, we needto estimate uncertainties due to pixelization of the sky as it
could be potential sources of errors. For an azimuthally symmetric mask W = 0 for |b| ≤ b0 , we do not need pixelization,
since Wi j can be analytically calculated as
Z π/2+b0
Wi j = 2π 1 −

π/2−b0!
Yi (θ, 0)Y j (θ, 0) sin θdθ ,(11)where θ is the polar angle, and θ = π/2 − b. We compare
the eigenvalues of W obtained by using Eq. (4) and Eq. (11)to study the accuracy of the calculated matrix. The differ-

4
ence depends on the resolution of pixels Nside . We adopt
Nside = 16, 128, 512, 1024 and find that the differences in the
eigenvalues are less than 0.01% for Nside ≥ 512. In order toevaluate the accuracy of map reconstruction, we use the L2
norm of the fractional difference between the reconstructed
map and the original map,
PNpix
D2` ≡

i

[δrec (γ̂i ; `) − δtrue (γ̂i ; `)]2
PNpix 2
δtrue (γ̂i ; `)i(a) C` ∝ `0(12)where δrec and δtrue describe the reconstructed and the originaldensity fluctuations contributed from the `-th multipole mode,δ(γ̂; `) ≡X̀λk(b) C` ∝ `a`m Y`m (γ̂),

−2

(13)

m=−`

respectively. An overall fractional difference D2 is defined bysubstituting δ(γ̂) for δ(γ̂, `) in Eq. (12) whereP δ(γ̂) is a sum of
all the multipole contributions, i.e. δ(γ̂) = ` δ(γ̂; `).
In Fig. 2, we can see the difference between a pixelizedmask and an unpixelized “smooth” mask given by Eq. (11).
Using Healpix, we pixelize the mask with different resolution
Nside where the total number of pixel on the sky is defined as
2
Npix = 12 × Nside. Before looking into the pixelization effect,
we first study the overall behaviour of reconstruction accuracy
as a function of the number of iteration. As we increase thenumber of iteration, the accuracy of reconstruction gradually
improves at Nite < 10 then it begins to degrade. Finally, theD2 accuracy converges to the one obtained by using the bruteforce inversion (dotted lines in the Fig. 2 ). As we can see inthe bottom panel of Fig. 2, the differences in D2 between the
pixelized and unpixelized Wi j are significant for Nside = 16
and 128, while they are less than 1% at the D2 minimum for
Nside ≥ 512. We then conclude that the pixelization effect canbe ignored when we use sufficiently fine pixels with Nside ≥512.

λk(c) C` ∝ `

−4

λk
B.

Comparison with W −1 and SVD

Using Eq. (9), we can compare our method to the directbrute-force inversion of the matrix W and the singular value
decomposition[46].
If the mode coupling matrix W is invertible, we obtain theunique solution of the underlying density fluctuation within
the mask. However, if the matrix W contains some eigenvalues which are close to zero, where the matrix is close to being
singular, the inversion causes big errors. In such a case, we
can remove the singularity by replacing the small eigenvalues
with zero, so called the SVD method. More specifically, thematrix W can be decomposed into three matrices as,
W = U Σ V†,(14)

where V and U are imax × imax unitary matrices and a superscript † denotes the Hermitian conjugate. Σ is the diagonal
matrix which consists of the eigenvalues of W, or called as the
singular values. The order of the eigenvalues in Σ is arbitrary

FIG. 3: In the right panels, we plot D2 accuracy as a function of
iteration number. Three different rows show the different underlying
density fluctuations that have power-law spectrum of C` ∝ `0 , `−2and `−4 from top to bottom. For a red (n = 0) spectrum, we donot see any improvement. However, for a flat (n = −2) spectrum,
there exists an optimal number of iteration around Nite = 10 thatminimizes D2 . The optimal Nite depends on the `max . For a blue(n = −4) spectrum, D2 converges to a few percent accuracy aftersufficient number of iteration. For all cases, after a sufficient numberof iteration, D2 converges to the value which is obtained by the direct
inversion. In the left panels, we show the result for the SVD method
for the comparison as a function of eigenvalue threshold λk .

but they are arranged in a descending order so that the decomposition is determined uniquely. Let the k−th eigenvalue be
λk , and the eigenvalues smaller than λk are set to zero. Thenthe pseudo inversion of W is written in terms of Σ+ , the rank-k
diagonal matrix that consists of the reciprocal of the non-zero

5(a) C` ∝ `0(b) C` ∝ `−2(c) C` ∝ `−4

FIG. 4: The IHE reconstruction accuracy D2` for each multipole component. From right to left, the indices of underlying power spectrum aren = 0, −2 and −4, respectively. From top to bottom, the components of multipoles ` = 1 to 4 are shown. We reconstruct the masked map withmaximum multipoles of `max = 3, 5, 7 and 9. The error bars represent 1 σ standard deviation calculated from 1000 random realizations.

eigenvalues that are equal to or larger than λk , as
W + = V Σ+ U † ,

(15)

+
aest
i = ã j Wi j .

(16)which gives

The choice of threshold λk is not trivial, and should be carefully determined a priori because the mapping FSVD depends
on various factors including the threshold λk : i.e. FSVD =FSVD (`max , δ, W, λk ). In this context, the singularity can only
be defined in terms of λk . More specifically, the mapping
FSVD is contaminated by the eigenvalues which are close tobeing singular if ∂D2 /∂λk < 0. In such cases, the thresholdof eigenvalues should be increased to implement more accurate reconstruction. The brute-force inversion corresponds to
Finv = FSVD (`max , δ, W, λk = 0).
In the left panels of Fig. 3, D2 is shown as a function of
the minimum non-zero eigenvalue λk . The solid lines show
the different maximum multipole to be reconstructed, `max =3, 5, 7, 9 from dark to light color. The different three rows
show the different power of the input power spectra which
will be discussed in detail in the next section. Let us focus on
the middle panel. For `max = 3 and 5, D2 increases monotonically as we increase the threshold. It means that the FSVD is
not contaminated by singular eigenvalues as ∂D2 /∂λk > 0. On
the other hand, for `max = 7 and 9, D2 has a minimum aroundλk ∼ 0.1. Thus the eigenvalues smaller than ∼ 0.1 may be thecontaminant of the mapping F and we can better estimate theoriginal density fluctuations when we limit the eigenvalues as
λ < 0.1. In the right panel of Fig. 3, we show D2 for our IHE
method as a function of the number of iteration. As in Fig. 2,
D2 has the minimum around Nite ' 10 depending on the maximum multipole. After a sufficient number of iteration, D2
values converge to those obtained by using the brute-force inversion, which are shown as the horizontal dashed lines in

Fig.2. Our IHE method with a certain finite number of iteration is always better than the SVD method. However, in practice, we should know a priori the optimal number of iteration.
The number depends on the mask geometry, the maximum
multipole to be reconstructed, and the underlying spectrum of
the density fluctuation. We can estimate the optimal iteration
number by carrying out Monte-Carlo simulations. In the next
section we will see how the result changes for different typesof underlying power spectrum.
C. Dependence on underlying power spectrum

Reconstruction accuracy depends on the underlying power
spectrum C` ∝ `n of the fluctuation as each harmonic mode
is not independent in the masked incomplete sky even the underlying fluctuation is Gaussian. We consider three power law
indices, n = 0, −2 and −4, which correspond to the following
three cases. The projected two dimensional galaxy or dark
matter distribution is approximated as C`g ∝ `0 [e.g. 47] onlarge scales, and the ordinary Sachs-Wolfe spectrum is approximated as C`SW ∝ `−2 [48]. On very large-angular scales,
the integrated Sachs-Wolfe effect, which gives C`iSW ∝ `−4[e.g. 49] dominates the CMB in the standard ΛCDM scenario.
Given the typical scale of mask θM , it is impossible to reconstruct the fluctuation whose scale is smaller than `M ≥180/θM . Because of the mode coupling, fluctuations with angular scales corresponding to `M are strongly affected by fluctuations with smaller angular sizes. If the spectral index isnegative, the amplitude of a smaller scale fluctuation is weakand does not strongly disturb the large scale fluctuations. Thus
the deconvolution mapping F is less affected by singularities.
However, if the slope of the spectrum is flat or positive, large
scale modes are highly contaminated by the noisy small scale
fluctuations.
In Fig. 3, the top and bottom panels show D2 accuracy for
n = 0 and n = −4. For n = 0 case, the SVD has a min-

6
imum at λk ∼ 0.5 that is larger than the case of n = −2.
It means that FSVD (n = 0) is more affected by singularities
than FSVD (n = −2). The IHE result in the right panels shows
Nite = 1 gives the best accuracy and it gradually degrades tothe value given by the brute-force method. On the other hand,
FSVD (n = −4) is not affected by singularities because it always shows ∂D2 /∂λk < 0 for the SVD method. In Fig. 4,
we also show the reconstruction accuracy for each multipole
component, D2` . The behaviour is similar to the one describedin Fig. 3. The reconstruction accuracy is always better for low
multipoles and that the optimal number of iteration depends
also on the multipole component we are going to reconstruct.
IV. SUMMARY

the missing data on an azimuthally symmetric mask. We
have considered three types of gaussian fluctuations with thepower-law indices n = 0, −2 and −4, which correspond to thegalaxy power spectrum, the ordinary Sachs-Wolfe spectrum,and the integrated Sachs-Wolfe spectrum, respectively. For
n = −2 case, we have found that there exists an optimal finite number of iteration that makes the reconstruction moreaccurate than the SVD method or the brute-force method. For
n = 0 case, the pseudo-a`m is the best estimator for the projected density fluctuations. For n = −4 case, the brute-force
inversion method gives the best accuracy. In that case, the IHE
method can help to reduce the computation time of inversion.
It would be interesting to see how the significance of thelarge-angle CMB anomaly changes when we use different
methods of map reconstruction. The IHE method may help
to unveil the origin of the CMB anomaly. We will explore thisproblem in our future work.

We have developed a new method based on the iterative harmonic expansions (IHE) for reconstructing the missing fluctuations in the masked region in the sky. Reconstructing thedata at a masked region is known as an inverse problem. Our
method is equivalent to the brute-force inversion in the limitthat the iteration number goes to infinity. However, in somecases, the finite truncation of the iteration gives a better estimate of the underlying fluctuation. The reconstruction accuracy depends on the geometrical shape of the mask, the maximum multipole mode to be included in the IHE analysis, thetilt of the underlying power spectrum, and the multipole component of the map to be reconstructed.
As an example, we have applied the method to reconstruct

We thank Masahiro Takada, Eiichiro Komatsu, Issha Kayo,
Takahiro Nishimichi for useful discussions. AN is supported
in part by the Grant-in-Aid for the Scientific Research Fund
(No. 23340061), by JSPS Core-to-Core Program “ International Research Network for Dark Energy”, by World Premier International Research Center Initiative (WPI Initiative),
MEXT, Japan and by the FIRST program “Subaru Measurements of Images and Redshifts (SuMIRe)”, CSTP, Japan.

[1] J. P. Ralston and P. Jain, International Journal of Modern
Physics D 13, 1857 (2004), arXiv:astro-ph/0311430.
[2] A. de Oliveira-Costa, M. Tegmark, M. Zaldarriaga, and
A. Hamilton, Phys. Rev. D 69, 063516 (2004), arXiv:astroph/0307282.
[3] F. K. Hansen, A. J. Banday, and K. M. Górski, MNRAS 354,
641 (2004), arXiv:astro-ph/0404206.
[4] A. Hajian, T. Souradeep, and N. Cornish, ApJL 618, L63(2005), arXiv:astro-ph/0406354.
[5] J. W. Moffat, JCAP 10, 012 (2005), arXiv:astro-ph/0502110.
[6] K. Land and J. Magueijo, MNRAS 367, 1714 (2006),
arXiv:astro-ph/0509752.
[7] A. Bernui, T. Villela, C. A. Wuensche, R. Leonardi, and I. Ferreira, Astronomy & Astrophysics 454, 409 (2006), arXiv:astroph/0601593.
[8] C. J. Copi, D. Huterer, D. J. Schwarz, and G. D. Starkman,
Phys. Rev. D 75, 023507 (2007), arXiv:astro-ph/0605135.
[9] H. K. Eriksen, A. J. Banday, K. M. Górski, F. K. Hansen, andP. B. Lilje, ApJL 660, L81 (2007), arXiv:astro-ph/0701089.
[10] C. Monteserı́n, R. B. Barreiro, P. Vielva, E. Martı́nez-González,
M. P. Hobson, and A. N. Lasenby, MNRAS 387, 209 (2008),
0706.4289.
[11] P. K. Samal, R. Saha, P. Jain, and J. P. Ralston, MNRAS 396,
511 (2009), 0811.1639.
[12] Planck Collaboration, P. A. R. Ade, N. Aghanim, C. ArmitageCaplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont, C. Baccigalupi, A. J. Banday, et al., ArXiv e-prints(2013), 1303.5083.

[13] R. Aurich, S. Lustig, and F. Steiner, Classical and Quantum
Gravity 27, 095009 (2010), 0903.3133.
[14] A. Pontzen and H. V. Peiris, Phys. Rev. D 81, 103008 (2010),
1004.2706.
[15] G. Efstathiou, Y.-Z. Ma, and D. Hanson, MNRAS 407, 2530
(2010), 0911.5399.
[16] C. L. Bennett, R. S. Hill, G. Hinshaw, D. Larson, K. M. Smith,
J. Dunkley, B. Gold, M. Halpern, N. Jarosik, A. Kogut, et al.,ApJS 192, 17 (2011), 1001.4758.
[17] L. R. Abramo, L. Sodré, Jr., and C. A. Wuensche, Phys. Rev. D
74, 083515 (2006), arXiv:astro-ph/0605269.
[18] M. Cruz, P. Vielva, E. Martı́nez-González, and R. B. Barreiro,
MNRAS 412, 2383 (2011), 1005.1264.
[19] M. Hansen, J. Kim, A. M. Frejsel, S. Ramazanov, P. Naselsky,
W. Zhao, and C. Burigana, JCAP 10, 059 (2012), 1206.6981.
[20] K. T. Inoue and J. Silk, Astrophys. J. 648, 23 (2006),
arXiv:astro-ph/0602478.
[21] K. T. Inoue and J. Silk, Astrophys. J. 664, 650 (2007),
arXiv:astro-ph/0612347.
[22] A. Rassat, K. Land, O. Lahav, and F. B. Abdalla, MNRAS 377,
1085 (2007), arXiv:astro-ph/0610911.
[23] N. Afshordi, G. Geshnizjani, and J. Khoury, JCAP 8, 030(2009), 0812.2244.
[24] C. L. Francis and J. A. Peacock, MNRAS 406, 14 (2010),
0909.2495.
[25] A. Rassat, J.-L. Starck, and F.-X. Dupe, ArXiv e-prints (2013),
1303.4727.
[26] K. Tomita and K. T. Inoue, Phys. Rev. D 77, 103522 (2008),

Acknowledgments

70712.1291.
[27] N. Sakai and K. T. Inoue, Phys. Rev. D 78, 063510 (2008),
0805.3446.
[28] K. T. Inoue, MNRAS 421, 2731 (2012), 1109.4527.
[29] H. V. Peiris and T. L. Smith, Phys. Rev. D 81, 123517 (2010),
1002.0836.
[30] D. Hanson, A. Lewis, and A. Challinor, Phys. Rev. D 81,
103003 (2010), 1003.0198.
[31] J. Kim, P. Naselsky, and N. Mandolesi, ApJL 750, L9 (2012),
1202.0188.
[32] R. Aurich, S. Lustig, F. Steiner, and H. Then, Classical andQuantum Gravity 24, 1879 (2007), arXiv:astro-ph/0612308.
[33] A. Emir Gümrükçüoglu, C. R. Contaldi, and M. Peloso, JCAP
11, 005 (2007), 0707.4179.
[34] D. C. Rodrigues, Phys. Rev. D 77, 023534 (2008), 0708.1168.
[35] A. Bernui and W. S. Hipólito-Ricaldi, MNRAS 389, 1453
(2008), 0807.1076.
[36] M. Cruz, E. Martı́nez-González, P. Vielva, J. M. Diego,
M. Hobson, and N. Turok, MNRAS 390, 913 (2008),
0804.2904.
[37] A. Fialkov, N. Itzhaki, and E. D. Kovetz, JCAP 2, 004 (2010),
0911.2100.
[38] H. Zheng and E. F. Bunn, Phys. Rev. D 82, 063533 (2010),
1003.5548.
[39] H. Liu, A. M. Frejsel, and P. Naselsky, ArXiv e-prints (2013),

1302.6080.
[40] E. Hivon, K. M. Górski, C. B. Netterfield, B. P. Crill, S. Prunet,and F. Hansen, Astrophys. J. 567, 2 (2002).
[41] P. Abrial, Y. Moudden, J.-L. Starck, J. Fadili, J. Delabrouille,
and M. K. Nguyen, Statistical Methodology 5, 289 (2008),
0804.1295.
[42] M. Bucher and T. Louis, MNRAS 424, 1694 (2012),
1109.0286.
[43] J.-L. Starck, M. J. Fadili, and A. Rassat, Astronomy & Astrophysics 550, A15 (2013), 1210.6587.
[44] K. T. Inoue, P. Cabella, and E. Komatsu, Phys. Rev. D 77,
123539 (2008), 0804.0527.
[45] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K.
Hansen, M. Reinecke, and M. Bartelmann, Astrophys. J. 622,
759 (2005).
[46] G. Efstathiou, MNRAS 348, 885 (2004), arXiv:astroph/0310207.
[47] W. J. Frith, P. J. Outram, and T. Shanks, MNRAS 364, 593(2005), arXiv:astro-ph/0507215.
[48] R. K. Sachs and A. M. Wolfe, Astrophys. J. 147, 73 (1967).
[49] A. Cooray, Phys. Rev. D 65, 083518 (2002), arXiv:astroph/0109162.
[50] http://healpix.jpl.nasa.gov/