Issue 
A&A
Volume 571, November 2014
Planck 2013 results



Article Number  A23  
Number of page(s)  48  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201321534  
Published online  29 October 2014 
Planck 2013 results. XXIII. Isotropy and statistics of the CMB
^{1} APC, AstroParticule et Cosmologie,
Université Paris Diderot, CNRS/IN2P3, CEA/lrfu, Observatoire de Paris,
Sorbonne Paris Cité, 10 rue Alice Domon et Léonie
Duquet, 75205
Paris Cedex 13,
France
^{2} Aalto University Metsähovi Radio
Observatory and Dept of Radio Science and Engineering, PO Box 13000, 00076
Aalto,
Finland
^{3} African Institute for Mathematical
Sciences, 68 Melrose Road, 7945
Muizenberg, Cape
Town, South Africa
^{4} Agenzia Spaziale Italiana Science
Data Center, via del Politecnico snc, 00133
Roma,
Italy
^{5} Agenzia Spaziale Italiana, viale
Liegi 26, 00198
Roma,
Italy
^{6} Astrophysics Group, Cavendish
Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge
CB3 0HE,
UK
^{7} Astrophysics and Cosmology Research
Unit, School of Mathematics, Statistics and Computer Science, University of
KwaZuluNatal, Westville Campus,
Private Bag X54001, 4000
Durban, South
Africa
^{8} CITA, University of
Toronto, 60 St. George St.,
Toronto, ON
M5S 3H8,
Canada
^{9} CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028
Toulouse Cedex 4,
France
^{10} California Institute of
Technology, Pasadena,
California,
USA
^{11} Centre for Theoretical Cosmology,
DAMTP, University of Cambridge, Wilberforce Road, Cambridge
CB3 0WA,
UK
^{12} Centro de Estudios de Física del
Cosmos de Aragón (CEFCA), Plaza San Juan, 1, planta 2, 44001
Teruel,
Spain
^{13} Computational Cosmology Center,
Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{14} Consejo Superior de Investigaciones
Científicas (CSIC), Madrid, Spain
^{15} DSM/Irfu/SPP,
CEASaclay, 91191
GifsurYvette Cedex,
France
^{16} DTU Space, National Space
Institute, Technical University of Denmark, Elektrovej 327, 2800
Kgs. Lyngby,
Denmark
^{17} Département de Physique Théorique,
Université de Genève, 24 quai E.
Ansermet, 1211
Genève 4,
Switzerland
^{18} Departamento de Física Fundamental,
Facultad de Ciencias, Universidad de Salamanca, 37008
Salamanca,
Spain
^{19} Departamento de Física, Universidad
de Oviedo, Avda. Calvo Sotelo
s/n, 33007
Oviedo,
Spain
^{20} Departamento de Matemáticas,
Estadística y Computación, Universidad de Cantabria, Avda. de los Castros s/n, 39005
Santander,
Spain
^{21} Department of Astronomy and
Astrophysics, University of Toronto, 50 Saint George Street, Toronto, Ontario,
Canada
^{22} Department of Astrophysics/IMAPP,
Radboud University Nijmegen, PO Box
9010, 6500 GL
Nijmegen, The
Netherlands
^{23} Department of Electrical
Engineering and Computer Sciences, University of California,
Berkeley, California,
USA
^{24} Department of Physics &
Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver,
British Columbia,
Canada
^{25} Department of Physics and
Astronomy, Dana and David Dornsife College of Letter, Arts and Sciences, University of
Southern California, Los
Angeles, CA
90089,
USA
^{26} Department of Physics and
Astronomy, University College London, London
WC1E 6BT,
UK
^{27} Department of Physics, Florida
State University, Keen Physics
Building, 77 Chieftan Way, Tallahassee, Florida, USA
^{28} Department of Physics, Gustaf
Hällströmin katu 2a, University of Helsinki, 00014
Helsinki,
Finland
^{29} Department of Physics, Princeton
University, Princeton, New
Jersey, USA
^{30} Department of Physics, University
of Alberta, 1132289
Avenue, Edmonton,
Alberta, T6G 2G7,
Canada
^{31} Department of Physics, University
of California, Berkeley, California, USA
^{32} Department of Physics, University
of California, One Shields
Avenue, Davis,
California,
USA
^{33} Department of Physics, University
of California, Santa
Barbara, California, USA
^{34} Department of Physics, University
of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, Illinois,
USA
^{35} Dipartimento di Fisica e Astronomia
G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131
Padova,
Italy
^{36} Dipartimento di Fisica e
Astronomia, Università degli Studi di Bologna, viale Berti Pichat 6/2, 40127
Bologna,
Italy
^{37} Dipartimento di Fisica e Scienze
della Terra, Università di Ferrara, via Saragat 1, 44122
Ferrara,
Italy
^{38} Dipartimento di Fisica, Università
La Sapienza, P. le A. Moro
2, 00185
Roma,
Italy
^{39} Dipartimento di Fisica, Università
degli Studi di Milano, via Celoria,
16, 20133
Milano,
Italy
^{40} Dipartimento di Fisica, Università
degli Studi di Trieste, via A.
Valerio 2, 34127
Trieste,
Italy
^{41} Dipartimento di Fisica, Università
di Roma Tor Vergata, via della
Ricerca Scientifica, 1, Roma, Italy
^{42} Dipartimento di Matematica,
Università di Roma Tor Vergata, via
della Ricerca Scientifica, 1, 00133
Roma,
Italy
^{43} Discovery Center, Niels Bohr
Institute, Blegdamsvej
17, 2100
Copenhagen,
Denmark
^{44} Dpto. Astrofísica, Universidad de
La Laguna (ULL), 38206 La
Laguna, Tenerife,
Spain
^{45} European Space Agency, ESAC, Planck
Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo,
Villanueva de la Cañada, Madrid, Spain
^{46} European Space Agency, ESTEC,
Keplerlaan 1, 2201
AZ
Noordwijk, The
Netherlands
^{47} Helsinki Institute of Physics,
Gustaf Hällströmin katu 2, University of Helsinki, Helsinki,
Finland
^{48} INAF – Osservatorio Astronomico di
Padova, Vicolo dell’Osservatorio 5, Padova, Italy
^{49} INAF – Osservatorio Astronomico di
Roma, via di Frascati 33, Monte
Porzio Catone, Italy
^{50} INAF – Osservatorio Astronomico di
Trieste, via G.B. Tiepolo 11, Trieste, Italy
^{51} INAF Istituto di Radioastronomia,
via P. Gobetti 101, 40129
Bologna,
Italy
^{52} INAF/IASF Bologna, via Gobetti
101, 40129
Bologna,
Italy
^{53} INAF/IASF Milano, via E. Bassini
15, 20133
Milano,
Italy
^{54} INFN, Sezione di Bologna, via
Irnerio 46, 40126
Bologna,
Italy
^{55} INFN, Sezione di Roma 1, Università
di Roma Sapienza, Piazzale Aldo
Moro 2, 00185
Roma,
Italy
^{56} IPAG: Institut de Planétologie et
d’Astrophysique de Grenoble, Université Joseph Fourier, Grenoble 1/CNRSINSU, UMR
5274, 38041
Grenoble,
France
^{57} ISDC Data Centre for Astrophysics,
University of Geneva, Ch. d’Ecogia
16, 1290
Versoix,
Switzerland
^{58} IUCAA, Post Bag 4, Ganeshkhind,
Pune University Campus, 411
007
Pune,
India
^{59} Imperial College London,
Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7
2AZ, UK
^{60} Infrared Processing and Analysis
Center, California Institute of Technology, Pasadena, CA
91125,
USA
^{61} Institut Néel, CNRS, Université
Joseph Fourier Grenoble I, 25 rue
des Martyrs, 38054
Grenoble,
France
^{62} Institut Universitaire de
France, 103 bd
SaintMichel, 75005
Paris,
France
^{63} Institut d’Astrophysique Spatiale,
CNRS (UMR 8617) Université ParisSud 11, Bâtiment 121, 91405
Orsay,
France
^{64} Institut d’Astrophysique de Paris,
CNRS (UMR 7095), 98bis boulevard
Arago, 75014
Paris,
France
^{65} Institute for Space
Sciences, 077125
BucharestMagurale,
Romania
^{66} Institute of Astronomy and
Astrophysics, Academia Sinica, 106
Taipei,
Taiwan
^{67} Institute of Astronomy, University
of Cambridge, Madingley
Road, Cambridge
CB3 0HA,
UK
^{68} Institute of Theoretical
Astrophysics, University of Oslo, Blindern, 0315
Oslo,
Norway
^{69} Instituto de Astrofísica de
Canarias, C/Vía Láctea s/n, La
Laguna, 38200
Tenerife,
Spain
^{70} Instituto de Física de Cantabria
(CSICUniversidad de Cantabria), Avda. de los Castros s/n, 39005
Santander,
Spain
^{71} Jet Propulsion Laboratory,
California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, USA
^{72} Jodrell Bank Centre for
Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of
Manchester, Oxford
Road, Manchester,
M13 9PL,
UK
^{73} Kavli Institute for Cosmology
Cambridge, Madingley
Road, Cambridge,
CB3 0HA,
UK
^{74} LAL, Université ParisSud,
CNRS/IN2P3, 91898
Orsay,
France
^{75} LERMA, CNRS, Observatoire de Paris,
61 avenue de l’Observatoire, 75014
Paris,
France
^{76} Laboratoire AIM, IRFU/Service
d’Astrophysique – CEA/DSM – CNRS – Université Paris Diderot, Bât. 709,
CEASaclay, 91191
GifsurYvette Cedex,
France
^{77} Laboratoire Traitement et
Communication de l’Information, CNRS (UMR 5141) and Télécom ParisTech,
46 rue Barrault, 75634
Paris Cedex 13,
France
^{78} Laboratoire de Physique Subatomique
et de Cosmologie, Université Joseph Fourier Grenoble I, CNRS/IN2P3, Institut National
Polytechnique de Grenoble, 53 rue
des Martyrs, 38026
Grenoble Cedex,
France
^{79} Laboratoire de Physique Théorique,
Université ParisSud 11 & CNRS, Bâtiment 210, 91405
Orsay,
France
^{80} Lawrence Berkeley National
Laboratory, Berkeley,
California,
USA
^{81} MaxPlanckInstitut für Astrophysik,
KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{82} MaxPlanckInstitut für
Extraterrestrische Physik, Giessenbachstraße, 85748
Garching,
Germany
^{83} McGill Physics, Ernest Rutherford
Physics Building, McGill University, 3600 rue University, Montréal, QC, H3A 2T8,
Canada
^{84} MilliLab, VTT Technical Research
Centre of Finland, Tietotie 3, Espoo, Finland
^{85} Mullard Space Science Laboratory,
University College London, Surrey
RH5 6NT,
UK
^{86} National University of Ireland,
Department of Experimental Physics, Maynooth, Co.
Kildare, Ireland
^{87} Niels Bohr Institute,
Blegdamsvej 17, Copenhagen,
Denmark
^{88} Observational Cosmology, Mail Stop
36717, California Institute of Technology, Pasadena, CA
91125,
USA
^{89} Optical Science Laboratory,
University College London, Gower
Street, London,
UK
^{90} SBITPLPPC, EPFL,
1015, Lausanne,
Switzerland
^{91} SISSA, Astrophysics Sector, via
Bonomea 265, 34136
Trieste,
Italy
^{92} School of Physics and Astronomy,
Cardiff University, Queens
Buildings, The Parade, Cardiff, CF24
3AA, UK
^{93} School of Physics and Astronomy,
University of Nottingham, Nottingham
NG7 2RD,
UK
^{94} Space Sciences Laboratory,
University of California, Berkeley, California, USA
^{95} Special Astrophysical Observatory,
Russian Academy of Sciences, Nizhnij Arkhyz, Zelenchukskiy region, 369167
KarachaiCherkessian Republic,
Russia
^{96} Stanford University,
Dept of Physics, Varian Physics Bldg, 382 via
Pueblo Mall, Stanford, California, USA
^{97} SubDepartment of Astrophysics,
University of Oxford, Keble
Road, Oxford
OX1 3RH,
UK
^{98} Theory Division, PHTH, CERN, 1211, Geneva 23, Switzerland
^{99} UPMC Univ Paris 06,
UMR7095, 98bis boulevard
Arago, 75014
Paris,
France
^{100} Université de
Toulouse, UPSOMP,
IRAP, 31028
Toulouse Cedex 4,
France
^{101} University of Granada,
Departamento de Física Teórica y del Cosmos, Facultad de Ciencias,
Granada,
Spain
^{102} Warsaw University
Observatory, Aleje Ujazdowskie
4, 00478
Warszawa,
Poland
Received:
21
March
2013
Accepted:
11
January
2014
The two fundamental assumptions of the standard cosmological model – that the initial fluctuations are statistically isotropic and Gaussian – are rigorously tested using maps of the cosmic microwave background (CMB) anisotropy from the Planck satellite. The detailed results are based on studies of four independent estimates of the CMB that are compared to simulations using a fiducial ΛCDM model and incorporating essential aspects of the Planck measurement process. Deviations from isotropy have been found and demonstrated to be robust against component separation algorithm, mask choice, and frequency dependence. Many of these anomalies were previously observed in the WMAP data, and are now confirmed at similar levels of significance (about 3σ). However, we find little evidence of nonGaussianity, with the exception of a few statistical signatures that seem to be associated with specific anomalies. In particular, we find that the quadrupoleoctopole alignment is also connected to a low observed variance in the CMB signal. A power asymmetry is now found to persist on scales corresponding to about ℓ = 600 and can be described in the lowℓ regime by a phenomenological dipole modulation model. However, any primordial power asymmetry is strongly scaledependent and does not extend toarbitrarily small angular scales. Finally, it is plausible that some of these features may be reflected in the angular power spectrum of the data, which shows a deficit of power on similar scales. Indeed, when the power spectra of two hemispheres defined by a preferred direction are considered separately, one shows evidence of a deficit in power, while its opposite contains oscillations between odd and even modes that may be related to the parity violation and phase correlations also detected in the data. Although these analyses represent a step forward in building an understanding of the anomalies, a satisfactory explanation based on physically motivated models is still lacking.
Key words: cosmic background radiation / cosmology: observations / cosmology: miscellaneous
© ESO, 2014
1. Introduction
This paper, one of a set associated with the 2013 release of data from the Planck^{1} mission (Planck Collaboration I 2014) describes a set of studies undertaken to determine the statistical properties of the cosmic microwave background (CMB).
The standard cosmological model is described well by the FriedmannLemaîtreRobertsonWalker (FLRW) solution of the Einstein field equations. This model is characterized by a homogeneous and isotropic metric and a scale factor of the expanding Universe. At very early times it is hypothesized that the Universe went through a period of accelerated expansion, the socalled cosmological inflation, driven by a hypothetical scalar field, the inflaton. During inflation the universe behaves as a de Sitter space, providing the conditions in which some of the present properties of the universe can be realized and specifically relaxing the problem of initial conditions. In particular, the seeds that gave rise to the present largescale matter distribution via gravitational instability originated as quantum fluctuations of the inflaton about its vacuum state. These fluctuations in the inflaton produce energy perturbations that are distributed as a homogeneous and isotropic Gaussian random field. Linear theory relates those energy fluctuations to the CMB anisotropies, implying a distribution for the anisotropies very close to that of an isotropic Gaussian random field.
The scope of this work is to use Planck data to test the Gaussianity and near isotropy of the CMB in intensity, as expected in the standard cosmology paradigm. Testing these fundamental properties is crucial for the validation of the standard cosmological scenario, and has profound implications for our understanding of the physical nature of the Universe and the initial conditions of structure formation. Moreover, the confirmation of the isotropic and Gaussian nature of the CMB is essential for justifying the corresponding assumptions usually made when estimating the CMB power spectra and other quantities to be obtained from the Planck data. Conversely, the detection of significant deviations from these assumptions that are not consistent with known systematic effects or foreground residuals would necessitate major revision of current methodological approaches for derivating the mission’s many science results.
Significant deviations from Gaussianity are expected from nonlinear processes that lead to secondary anisotropies, e.g., the integrated SachsWolfe (ISW) effect and lensing. Indeed, these effects are the subject of two companion Planck papers (Planck Collaboration XVII 2014; Planck Collaboration XIX 2014, respectively). However, remarkably, a number of anomalies, by which we mean features of the observed sky that are not statistically consistent with the bestfit ΛCDM model, have been found in the WMAP data. Indeed, the WMAP team (Spergel et al. 2003) themselves initially proposed some intriguing discrepancies in the form of a lack of power on large angular scales, which was subsequently reconfirmed for the threeyear data in Copi et al. (2007). Further examples include an alignment of the loworder multipoles (Tegmark et al. 2003; Schwarz et al. 2004; Bielewicz et al. 2005; Land & Magueijo 2005a), a northsouth asymmetry in both power spectra (Eriksen et al. 2004a; Hansen et al. 2009) and various measures of nonGaussianity (Eriksen et al. 2004c, 2005; Räth et al. 2007a), parity asymmetry in the power spectrum corresponding to large angular scales (Kim & Naselsky 2010a), and a region of significant temperature decrement – the socalled cold spot (Vielva et al. 2004; Cruz et al. 2005).
Whilst WMAP have presented refutations of these anomalies, either by criticizing the robustness of the statistical methods employed (Bennett et al. 2011) or by associating them with systematic artefacts of the data processing that have been corrected in the nineyear data release (Bennett et al. 2013), Planck represents a unique opportunity to independently assess their existence. Its higher angular resolution and sensitivity and wider frequency range enable a better understanding and removal of the Galactic and extragalactic foregrounds thus allowing a larger fraction of the sky to be useful for performing isotropy and Gaussianity analysis and to confirm and interpret those anomalies.
Throughout this paper, we quantify the significance of the test statistic in terms of the pvalue. This is the probability of obtaining a test statistic at least as extreme as the observed one, under the assumption that the null hypothesis (i.e., Gaussianity and isotropy of the CMB) is true. In some tests, where it is very justified to only use a onetailed probability, the pvalue is replaced by the corresponding upper or lowertail probability. A low pvalue is indicative of a tension between the data and the assumed statistical model (the null hypothesis). This can arise either when the assumed cosmological model is incorrect, if unknown or unmodelled aspects of the foreground emission or the measurement process exist, or as a result of a natural statistical fluctuation. The most interesting possibility, of course, is that a low pvalue is an indication of new physics.
From the theoretical point of view, there are many variants of inflation that predict high levels of nonGaussianity and new scenarios motivated by string and Mtheory. In addition, there are many physical effects that might give rise to a deviation from isotropy or the presence of nonGaussianity. Those deviations may be classified according to their physical nature and origin as follows: nonstandard inflationary models, geometry and topology of the Universe, and topological defects. The main results from these areas, as well as the detailed descriptions of methodologies and of specific theoreticallymotivated model constraints, are provided in the companion papers (Planck Collaboration XXIV 2014; Planck Collaboration XXVI 2014; Planck Collaboration XXV 2014).
This paper covers all relevant aspects related to the phenomenological study of the statistical isotropy and Gaussian nature of the CMB measured by the Planck satellite. It is organized as follows. Section 2 describes the Planck data used for the analyses. Section 3 explains the main characteristics of the simulations that constitute our reference set of Gaussian sky maps representative of the null hypothesis. In Sect. 4 the null hypothesis is tested with a number of standard tests that probe different aspects of nonGaussianity. The WMAP anomalies are revisited in the light of the Planck data in Sect. 5. In Sect. 6 the implications of the found deviations of the null hypothesis on C_{ℓ} and cosmological parameters estimations are discussed. Finally, Sect. 7 provides the main conclusions of the paper.
2. Data description
In this paper, we utilize data from the Planck2013 data release corresponding to the nominal period of the Planck mission. In part, this comprises sky maps at nine frequencies, with corresponding “halfring” maps that are generated by separating the data for a given pointing period into two halves, plus maps generated from data within the first and second survey periods. The maps are provided in HEALPix^{2} (Górski et al. 2005) format, with a pixel size defined by the N_{side} parameter. This set of maps allows a variety of consistency checks to be made, together with estimates of the instrumental noise contributions to analyses and limits on timevarying systematic artefacts. Full details are provided in a series of companion papers (Planck Collaboration II 2014; Planck Collaboration VI 2014).
Our main results are based on the CMB maps resulting from sophisticated component separation algorithms applied to the frequency maps, as detailed in Planck Collaboration XII (2014). The four methods – CommanderRuler, NILC, SEVEM and SMICA – are used to generate estimates of the CMB sky with an effective angular resolution of around 7′ or better, with accompanying symmetrized beam profiles, analysis masks, and halfring maps.
In Table 1 we list the different masks that have been used in the analyses described in this paper. These masks have been produced at full resolution (N_{side} = 2048) and are described in papers (Planck Collaboration XII 2014; Planck Collaboration XV 2014). In general, we make use of a standardized common mask that merges those associated with the individual methods (listed in Table 1 as U73). Where appropriate, we manipulate the masks for use at lower resolution, and, in particular, they have been degraded to lower resolutions (N_{side} = 1024, 512, 256, 128, 64, 32 and 16). The masks with resolutions N_{side} = 128–1024 have been degraded using the following procedure: first, the masks are degraded to their final resolution using the ud_degrade HEALPix routine, and then, a conservative approach is followed setting to zero any pixels with a value lower than 0.8. If the masks have to be degraded to even lower resolutions, N_{side} = 16–64, the procedure that has been used is different. First, the fullresolution mask is smoothed with a Gaussian kernel with a FWHM equal to three times the pixel size of the low resolution mask that we want to produce. Then the mask is degraded using ud_degrade to their final resolution. Finally, those pixels with a value lower or equal than 0.5 have been set to zero and the rest have been set to 1. This criterion, less conservative than the one used for the higher resolution masks, is a compromise between minimizing the amount of sky that is being masked and the level of contamination left unmasked. Note that in some cases the more conservative criterion of a 0.8 threshold has been used for the lower resolutions, as stated in the corresponding analyses.
List of the masks that have been used for the analyses described in this paper.
3. Simulations
The derivation of results to be presented in this paper requires extensive simulations, essential aspects of which include: 1) modelling the Planck instrumental effects that affect the quality of the data, including instrumental noise and identified systematic effects; 2) replicating the foreground removal approach and estimating the extent of foreground residuals; and 3) modelling the intrinsic statistical properties, Gaussian or otherwise, of the CMB signals expected from specific models of the Universe.
The full focal plane (FFP6) simulations described in Planck Collaboration (2013) provide a complete realization of the Planck mission capturing all characteristics of the scanning strategy, telescope, detector responses, and data reduction pipeline over the nominal mission period of 15.5 months. The Planck Sky Model (PSM) is used as input, encompassing the best current estimate of the microwave sky at Planck wavelengths including Galactic and extragalactic astrophysical foreground emission. The outputs include a complete set of maps for all detectors with accompanying halfring and survey splits generated for a reference CMB sky. These have been used to test and validate various analysis tools, employed in turn to evaluate the CMB component separation algorithms as applied to the data set. This also allows an FFP6based estimate of the foreground residuals remaining in the CMB sky after component separation to be evaluated, and their impact on various statistical estimators quantified.
An accompanying set of Monte Carlo simulations, that are allowed to vary within the cosmic variance, provides us with the reference set of Gaussian sky maps used for the null tests we employ. These simulations include FEBeCoP (Mitra et al. 2011) beam convolution at each of the Planck frequencies, which are then propagated through the various component separation pipelines using the same weights as derived from the Planck nominal mission data analysis. A fiducial CMB power spectrum has been adopted (and varied in the simulations within the cosmic variance) based on an analysis of the Planck data at an advanced, but not final stage of processing. Only small changes relative to the final Planck power spectrum presented in (Planck Collaboration XV 2014; Planck Collaboration XVI 2014) are observed.
4. Are the primordial fluctuations Gaussian?
As has been previously established, it is of major interest to determine whether the statistical properties of the primordial CMB anisotropies correspond to an isotropic Gaussian random field. Recent attempts to test this hypothesis have mainly relied on the WMAP data that have less sensitivity and cover a narrower frequency interval. Planck represents a unique opportunity to probe fundamental statistical properties of the Universe with cosmic variance limited sensitivity up to ℓ ≈ 2000 and minimum foreground contamination.
There is no unique signature of nonGaussianity, however, different tests can allow us to probe different types of nonGaussianity. As a consequence, it is important to subject the data to a variety of tests, and we do so in this section using a number of nonparametric tools. Specific signatures of nonGaussianity are sought in three companion papers (Planck Collaboration XXIV 2014; Planck Collaboration XXV 2014; Planck Collaboration XXVI 2014).
Any signal on the sphere, T(x), can be written in terms of the following spectral representation: (1)where x is a unit direction vector, Y_{ℓm} the spherical harmonics and (2)For an isotropic random field the a_{ℓm} coefficients are uncorrelated: (3)where δ is the Kronecker delta and C_{ℓ} is the angular power spectrum. Note that for a Gaussian and isotropic random field, the a_{ℓm} coefficients are independent and, if the mean is zero, the angular power spectrum provides a complete characterization of its statistical distribution.
In this paper, we examine the goodnessoffit of the data to the Planck bestfit fiducial CMB model, which constitutes our null hypothesis. The methods adopted constitute a broad range of statistical tools that allow the study of complementary statistical properties of the null hypothesis in both the real and harmonic space data representations. Claims of either consistency with the fiducial Planck cosmological model or of evidence for nonGaussianity must be demonstrably robust to data selection and specifics of the data analysis. Residuals from the diffuse Galactic foreground are likely to be nonGaussian in nature, and pointsources can be a source of nonGaussianity on small angular scales. In addition, the analysis of multifrequency data must be considered in order to confirm that any claimed nonGaussianity has a thermal (cosmological) spectrum. Moreover, the combined ISWlensing effect produces secondary anisotropies that significantly deviate from Gaussianity and whose effect has been detected in the Planck data (Planck Collaboration XIX 2014). This nonGaussian effect has to be considered when testing the null hypothesis.
We address these issues by analysing the cosmologically interesting subset of Planck frequency channels. Specifically, we analyse the uncorrected sky maps at 70, 100, 143 and 217 GHz, as a function of Galactic mask, to assess the likely contamination due to Galactic foregrounds. These tests have direct relevance for the Planck likelihood approach described in Planck Collaboration XV (2014) and the parameter estimation results presented in Planck Collaboration XVI (2014). We then consider the foreground cleaned versions of these maps generated by the SEVEM algorithm (Planck Collaboration XII 2014). Such a comparison also allows a semiindependent crosscheck of the cosmological signal seen by Planck LFI (70 GHz) and HFI (100, 143, 217 GHz). Although the cosmological content of the cleaned LFI and HFI data sets are independent, the cleaning makes use of difference maps generated from the remaining Planck frequency bands. Nevertheless, since the calibration and beam responses of the data are well understood over the ful range of frequencies, there will be no leakage of cosmological signal between the instrument specific frequencies.
We then continue with analyses of the CMB sky estimates provided by four component separation approaches (CommanderRuler, NILC, SEVEM, and SMICA) described in Planck Collaboration XII (2014), together with the corresponding mask appropriate for these methods. The largest sky area possible should be used for definitive statements about Gaussianity, since, in the absence of foreground residuals or systematic artefacts, it represents a superior sample of the Universe. Conversely, overly conservative sky cuts suffer from a loss of information.
4.1. One dimensional moments
In this section we perform some of the simplest Gaussianity tests, such as comparing the sample skewness and kurtosis of the data with simulations. The skewness, γ, and kurtosis, κ, of a random variable, X, are defined as where σ^{2}, the variance is: (6)Here X represents the pixel temperature at some N_{side} and the averages are performed over the sky. The skewness is a measure of the asymmetry of the probability distribution of a realvalued random variable. Qualitatively, a positive (negative) skew indicates that the tail on the right (left) side of the probability density function is longer than the left (right) side. A zero value indicates that the values are relatively evenly distributed on both sides of the mean, typically but not necessarily, implying a symmetric distribution. The kurtosis is a measure of the peakedness of the distribution and the heaviness of its tail. A distribution with positive (negative) excess kurtosis indicates that the distribution has a more acute (wider) peak around the mean and fatter (thinner) tails. Normal random variables have zero skewness and kurtosis.
The sample variance, σ^{2} is also considered in this section as a further consistency test, although it is not a normality test statistic. We begin by analysing the full resolution combined maps, applying the U73 mask for the four different component separation methods. The results for the variance, skewness, and kurtosis estimators are shown in Table 2. All four methods show similar results. The data are consistent with simulations for the skewness and kurtosis estimators, whereas the variance is anomalously low. This inconsistency was already reported for the WMAP data in Monteserín et al. (2008) and Cruz et al. (2011) at resolution N_{side} = 256 for a mask allowing slightly less sky coverage.
Lower tail probablity for the variance, skewness, and kurtosis estimators at N_{side} = 2048, using the U73 mask and four different component separation methods.
The mask dependence of our results is studied by recalculating the estimators using the CL58 and CL37 masks, which allow sky fractions of f_{sky} = 58% and f_{sky} = 37%, respectively. The SMICA cleaned maps at full resolution are considered. The most significant lower tail probability is obtained for the CL58, mask as can be seen in Table 3. The lower tail probabilities show a small dependence on the mask used, which could indicate either the presence of Galactic foreground residuals with larger sky coverage, or the increase of the sampling variance, and consequently a less significant probability, when a smaller fraction of the sky is considered.
Lower tail probablity for the variance, skewness, and kurtosis estimators at N_{side} = 2048, for the SMICA method, using different masks.
In order to identify any foreground contamination, the frequency dependence of our estimators is analysed. We use the SEVEM cleaned maps and the U73 mask. Note that as the 70 GHz full resolution noise is high we do not consider 70 GHz in this comparison. As the 100 GHz noise is not negligible we estimate the variance taking into account the noise dispersion as described in Cruz et al. (2011). The results are similar to those found for the combined map, as can be seen in Table 4. There is a small frequency dependence since the 100 GHz and 217 GHz maps show slightly higher variance and kurtosis than the 143 GHz map. However the 143 GHz map has a dominant contribution to the combined map, hence the foreground residuals in the combined map are likely to be small. The lower tail probabilities for the variance at 100 GHz, 143 GHz, and 217 GHz are respectively 0.023, 0.014, 0.025, whereas the skewness and kurtosis are compatible with simulations.
Lower tail probablity for the variance, skewness, and kurtosis estimators at N_{side} = 2048, for the SEVEM cleaned maps at different frequencies.
In order to study the hemispherical asymmetry reported in Cruz et al. (2011) we reanalyse the SMICA data and simulations considering independently the northern and southern ecliptic hemispheres outside the U73 mask. A clear asymmetry is found in the variance, with an anomalously low value found in the northern hemisphere, as seen in Table 3. Hemispherical asymmetries are further analysed in Sect. 5.3.
The results for different resolutions using the U73 mask are shown in Table 5. Note that the N_{side} = 2048 and 512 U73 masks have f_{sky} = 73%, while the low resolution masks at N_{side} = 64, 32, and 16 have f_{sky} = 78%. The variance is anomalously low at all the considered resolutions, whereas at low resolutions, the skewness is anomalously low and the kurtosis anomalously high. These results will be further analysed in Sect. 5.2. However, it is clear that, except on the largest angular scales, there is no evidence for nonGaussian behaviour in the data using these simple statistical measures.
Lower tail probablity for the variance, skewness, and kurtosis estimators at different resolutions, for the four component separation methods, using the U73 mask.
Fig. 1 Variance, skewness, and kurtosis for the combined map of the four different component separation methods. From top row to bottom row CR, NILC, SEVEM, SMICA. 
4.2. Npdf analysis
Under the assumption of Gaussianity, the Nprobability density function (Npdf) is given by a multivariate Gaussian function: (7)where T is a vector formed from the measured temperatures T(x) over all positions allowed by the applied mask, N_{pix} is the number of pixels in the vector, C is the covariance of the Gaussian field (of size N_{pix} × N_{pix}).
Unfortunately, the calculation of TC^{1}T^{T} is computationally unfeasible for the full Planck resolution at HEALPix N_{side} = 2048. At a lower resolution, the problem is tractable, and the noise level can also be considered negligible compared to the CMB signal. That implies that under the assumption of isotropy the covariance matrix C is fully defined by the Planck angular power spectrum (C_{ℓ}): (8)where C_{ij} is the covariance between pixels i and j, and θ_{ij} is angle between them, P_{ℓ} are the Legendre polynomials, b_{ℓ} is an effective window function associated with the N_{side} resolution, and ℓ_{max} is the maximum multipole probed.
Under the multivariate Gaussian hypothesis, the argument on the exponential in Eq. (7) should follow a χ^{2} distribution with N_{pix} degrees of freedom, or, equivalently (for N_{pix} ≫ 1) a normal distribution .
We begin by analysing the χ^{2} quantity for low resolution maps at N_{side} = 32 and filtering with a 5° FWHM Gaussian. 1 μK uncorrelated regularization noise is added to the covariance matrix before inverting it. Regularization noise realizations are added to the data and simulations for consistency (see Eriksen et al. 2007b, for more details).
Fig. 2 Npdf χ^{2} for the U73 mask, CL58, CL37, ecliptic North and ecliptic South. The different colours represent the four component separation methods, namely CR (green), NILC (blue), SEVEM (red), and SMICA (orange). 
We analyse the four cleaned data maps, applying the common, CL58 and CL37 masks. The admitted fraction of the sky is respectively 78%, 58% and 37%. The northern and southern ecliptic hemispheres outside the U73 mask are also considered independently in order to check for possible hemispherical asymmetries as those analysed in Sect. 5.3. The results are shown in Fig. 2 and Table 6. In the U73 mask case, the lower tail probabilities are low. Applying the two CL58 and CL37 masks that permit less sky coverage, gives data that are consistent with simulations. The low χ^{2} value appears to be localized in the northern ecliptic hemisphere. These results are directly comparable to the anomalous variance mentioned in Sect. 4.1. Note that the four maps show similar values, but the differences are larger when using the U73 mask. This could indicate the presence of some residual foreground contamination near the Galactic plane.
Lower tail probablity for the Npdf, using different masks.
Therefore, the frequency dependence of our estimator is analysed in order to identify any possible foreground contamination. The results are shown in Fig. 3 and Table 7. A moderate frequency dependence is found when using the U73 mask, which could indicate the presence of some foreground residuals near the Galactic plane. The frequency dependence of the results vanishes when using the CL58 and CL37 masks that exclude more of the sky from analysis.
4.3. Npoint correlation functions
In this section we present tests of the nonGaussianity of the Planck CMB maps using realspace Npoint correlation functions. While harmonicspace methods are often preferred over realspace methods for studying primordial fluctuations, realspace methods may have an advantage with respect to systematics and foregrounds, since such effects are usually localized in real space. It is therefore important to analyse the data in both spaces in order to highlight different features.
An Npoint correlation function is by definition the average product of N temperatures, measured in a fixed relative orientation on the sky, (9)where the unit vectors span an Npoint polygon on the sky. By assuming statistical isotropy, the Npoint functions are only functions of the shape and size of the Npoint polygon, and not on its particular position or orientation on the sky. Hence, the smallest number of parameters that uniquely determines the shape and size of the Npoint polygon is 2N − 3. In practice, the functions are estimated by simple product averages over all sets of N pixels fulfilling the geometric requirements set by θ_{1},...,θ_{2N − 3} characterizing the shape and size of the polygon, (10)Pixel weights can be introduced in order to reduce noise or mask boundary effects. Here they represent masking by being set to 1 for included pixels and to 0 for excluded pixels.
Fig. 3 Frequency dependence for 70 GHz (green), 100 GHz (blue), 143 GHz (red), and 217 GHz (orange), for different masks. 
Frequency dependence of the lower tail probablity for the Npdf, using different masks.
The main difficulty with computing Npoint functions is their computational scaling. The number of independent pixel combinations scales as , and for each combination of N pixels, 2N − 3 angular distances must be computed to uniquely determine the properties of the corresponding polygon. Computing the full Npoint function for N> 2 and N_{pix} ≳ 10^{5} is therefore computationally challenging. However, it is not necessary to include all possible Npoint configurations in order to produce interesting results. For instance, one may focus only on small angular scales, or on configurations with some special symmetry properties. By using the methods described in Eriksen et al. (2004b), the computational expense then becomes tractable, since no CPU time is spent on excluded configurations. In this paper several such subsets are computed, covering three distinct ranges of scales, namely small (up to 3°), intermediate (up to 10°) and large angular scales (the full range between 0° and 180°).
The shapes of the polygons selected for the analysis are the pseudocollapsed and equilateral configurations for the 3point function, and the rhombic configuration for the 4point function, composed of two equilateral triangles that share a common side. We use the same definition of pseudocollapsed as in Eriksen et al. (2005) i.e., as isosceles triangle where the length of the baseline falls within the second bin of the separation angles. The length of the longer edge of the triangle, θ, parametrizes its size. Analogously, in the case of the equilateral triangle and rhombus, the size of the polygon is parametrized by the length of the edge, θ. Note that these functions are chosen because of ease of implementation, not because they are better suited for testing Gaussianity than other configurations. In the following, all results refer to the reduced 4point function, i.e., corrected for the Gaussian contribution due to Wick’s theorem.
We analyse the CMB estimates downgraded to N_{side} = 64 and N_{side} = 512 as well as at the original resolution of N_{side} = 2048. In the case of the analysis at N_{side} = 64 the maps were additionally smoothed with FWHM of 165′ (three times the pixel size for the downgraded map). Due to computational limitations, an analysis is possible on the full sky only in the case of resolution N_{side} = 64. For the higher resolution maps, we perform the analysis on a set of nonoverlapping discs. For N_{side} = 512 we uniformly retain, after masking, part of the sky with approximately 100 discs of radius 10°. Analogously to the analysis by Eriksen et al. (2005), we consider two disc sets, A and B, with a relative offset between their grids such that the centres of the discs of set B are located in parts of the sky not covered by disc set A (see Fig. 4). For studies at the original resolution N_{side} = 2048, we restrict the analysis to 20 discs with a radius of 3° located randomly on an unmasked part of the sky (Fig. 4).
Fig. 4 Two sets of discs of radius 10° A (top) and B (middle), for the N_{side} = 512 CMB estimates; and a set of 20 randomly placed discs of radius 3° superimposed on the U73 mask (blue region) for the CMB estimates at N_{side} = 2048 (bottom). 
As in Eriksen et al. (2005), we consider the Npoint correlation functions averaged over the disc sets. In order to minimize correlations between the discs, we subtract from the maps at resolutions N_{side} = 512 and N_{side} = 2048 the bestfit multipoles computed for the ranges ℓ ≤ 18 and ℓ ≤ 60, respectively. This procedure corresponds in practice to a highpass filtering of the maps. One consequence of this procedure is the reduced sensitivity of the method to diffuse Galactic foreground residuals, which contribute most significantly to the low order multipoles. However, such a contribution is tested for in the analysis at resolution N_{side} = 64, in which highpass filtering is not employed beyond the removal of the monopole and dipole.
In the case of the analysis of the maps at resolution N_{side} = 512 the Npoint correlation functions were averaged over all discs for a given set, as well as separately for discs as a function of Galactic latitude. Such analysis allows the study of possible contamination of the maps by diffuse Galactic foreground residuals distributed symmetrically around the Galactic plane.
The low resolution versions of the U73 mask described earlier were used as required. Residual monopole and dipole contributions were then removed from the maps.
A simple χ^{2} test is chosen to quantify the degree of agreement between the simulations and the observations, where χ^{2} as usual is defined by (11)Here C_{N}(θ_{i}) is the Npoint correlation function for the bin of separation angle, θ_{i} (corresponding to the length of the longest side for the pseudocollapsed configuration and the length of the polygon’s side for the equilateral triangle and rhombus), ⟨C_{N}(θ_{i})⟩ is the corresponding average from the Monte Carlo (MC) ensemble, and (12)is the covariance matrix. Although the inverse of the covariance matrix constructed from MC simulations can be biased, it is relatively small for 1000 simulations and has a negligible impact on the significance levels estimated from the simulations, as described below.
Fig. 5 Pseudocollapsed 3point function averaged over disc set A for the raw (upper figure) and SEVEM foreground corrected (lower figure) 143 GHz map at N_{side} = 512 for different values of sky coverage (f_{sky}). Estimates of the multipoles for ℓ ≤ 18 are removed from the sky maps. The black solid line indicates the mean for 1000 MC simulations and the shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively, for the CG90 (f_{sky} = 0.9) mask. See Sect. 4.3 for the definition of the separation angle θ. 
Probabilities of obtaining values of the χ^{2} statistic for the concordance ΛCDM model at least as large as the observed values of the statistic for the raw 143 GHz (first row) and foreground corrected 143 GHz SEVEM CMB maps (second row).
This statistic is optimized for studying Gaussian distributed data. However, usually it also works quite well for mildly nonGaussian distributions, and in particular symmetric ones. Nevertheless, as for any statistic constructed from MC simulations, it can also be used for nonGaussian and asymmetrically distributed data. Below, we quote the significance level in terms of the fraction of simulations with a larger χ^{2} value than the observed map.
We analyse the mask dependence of the nonGaussianity of the maps using the pseudocollapsed 3point correlation function. The low order multipoles were estimated and subtracted from the maps for each mask separately. These estimates thus depend on the specific sky coverage employed, so that the comparison of the correlation functions for different masks is nontrivial. The function averaged over disc set A is shown in Fig. 5. The corresponding probabilities of obtaining values of the χ^{2} statistic for the concordance ΛCDM model at least as large as the observed values are given in Table 8.
Fig. 6 The 2point (upper left), pseudocollapsed (upper right), equilateral (lower left) 3point and reduced rhombic 4point (lower right) functions for the N_{side} = 64 CMB estimates. The black solid line indicates the mean for 1000 Monte Carlo simulations and the shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 
Fig. 7 The 2point (upper left), pseudocollapsed (upper right), equilateral (lower left) 3point and reduced rhombic 4point (lower right) functions averaged over disc set A for the N_{side} = 512 CMB estimates. Estimates of the multipoles for ℓ ≤ 18 are removed from the sky maps. The black solid line indicates the mean for 1000 Monte Carlo simulations and the shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 
Fig. 9 The 2point (upper left), pseudocollapsed (upper right), equilateral (lower left) 3point and reduced rhombic 4point (lower right) functions averaged over the disc set for the N_{side} = 2048 CMB estimates. Estimates of the multipoles for ℓ ≤ 60 are removed from the sky maps. The black solid line indicates the mean for 1000 Monte Carlo simulations and the shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 
In summary, the pseudocollapsed 3point function does not show any significant deviation from Gaussianity for the raw 143 GHz map masked with the CG60 (f_{sky} = 0.6) and CG70 (f_{sky} = 0.7) masks. To a lesser extent, this is true also for the CG80 (f_{sky} = 0.8) mask. We do not see any significant deviation for any of the analysed masks after cleaning the 143 GHz map using the SEVEM method.
The correlation functions for the four component separation methods and resolution parameters N_{side} = 64, N_{side} = 512 and N_{side} = 2048 are shown in Fig. 6, Fig. 7 (disc set A), Fig. 8 (disc set B). and Fig. 9, respectively. The probabilities of obtaining values of the χ^{2} statistic for the Planck fiducial ΛCDM model at least as large as the observed values are given in the Tables 9, 10, and 12, respectively. In addition, Table 11 reports the probabilities for the Npoint functions averaged over sets of discs as a function of Galactic latitude for the SMICA map at N_{side} = 512.
The results show consistency between the CMB maps estimated using the different component separation methods. No statistically significant deviations of the CMB maps from Gaussianity are found on any scale, with the exception of some analyses performed on sets of discs centred at low Galactic latitudes. However, for disc set B the fraction of such discs available for analysis after masking is very small, reducing the robustness of the results. A more intriguing case with very small probability arises for the disc set A 4point function averaged over discs centred at Galactic latitude ±10°. The fraction of discs available for the analysis is much larger in this case than for disc set B, and may indicate the presence of Galactic emission residuals.
It is also worth noticing that the CMB map smoothed and downgraded to N_{side} = 64 show the largest deviation, with a tendency towards superGaussianity, especially for the 4point correlation function, in comparison to the intermediate and small angular scale analyses. Since the CMB maps at this resolution were not highpass filtered, they may be more susceptible to Galactic foreground residual artefacts that survive the component separation processes.
4.4. Minkowski functionals
Minkowski functionals (MFs) describe the morphology of fields in any dimension and have long been used as estimators of nonGaussianity and anisotropy in the CMB (see e.g., Mecke et al. 1994; Schmalzing & Buchert 1997; Schmalzing & Gorski 1998; Komatsu et al. 2003; Eriksen et al. 2004c; Curto et al. 2007, 2008; De Troia et al. 2007; Spergel et al. 2007; Curto et al. 2008; Hikage et al. 2008; Komatsu et al. 2009). They are additive for disjoint regions of the sky and invariant under rotations and translations. Traditionally in the literature, the contours are defined by a threshold ν, usually given in units of the sky standard deviation (σ_{0}). We compute MFs for the regions colder and hotter than a given threshold ν. Thus, the three MFs, the area V_{0}(ν) = A(ν), the perimeter V_{1}(ν) = C(ν), and the genus V_{2}(ν) = G(ν), are defined respectively as: Here N_{ν} is the number of pixels where ΔT/σ_{0}>ν, N_{pix} is the total number of available pixels, A_{tot} is the total area of the available sky, N_{hot} is the number of compact hot spots, N_{cold} is the number of compact cold spots, and S_{i} is the contour length of each hot spot. We construct a fourth functional V_{3}(ν) = N_{cluster}(ν) which corresponds to N_{cold} for negative ν and N_{hot} for positive ν (Ducout et al. 2012). Analytical expressions for a Gaussian random field can be derived in terms of ν (see e.g. Vanmarcke 1983; Matsubara 2010) and give (16)with and (19)The amplitude A_{k} only depends on the shape of the power spectrum C_{ℓ}: where ω_{k} ≡ π^{k/ 2}/ Γ(k/ 2 + 1), which gives ω_{0} = 1, ω_{1} = 2, and ω_{2} = π with σ_{0} and σ_{1} beign respectively, the rms of the field and its first derivatives. These analytical expressions represent useful descriptions of the MFs which, for the case of a Gaussian random field, can be factorized as a function of the threshold and another of the shape of the C_{ℓ}. We will use both the unnormalized (V_{k}) and normalized (ν_{k}) MFs in the Gaussianity tests performed in this section. The unnormalized functionals are computed with a code that was used for the analysis of Archeops data (Curto et al. 2007) and has been thoroughly validated with Planck simulations, while for the normalized ones a code adapted to the high resolution Planck data and described in Ducout et al. (2012) is used.
By combining the MFs curves into a vector y of size n = n_{thresholds} × n_{functionals}, a null hypothesis test can be performed using a χ^{2} statistic given by (22)where y represents the MFs of the data, y_{G} those of the simulations and C is the covariance matrix. In order to assure convergence, in the case of the four normalized MFs C is estimated from 10^{4} Gaussian simulations, drawn from the Planck fiducial power spectrum, having the same instrumental properties of effective beam and noise as the data, the same applied mask and which have been processed in the same way to reach the corresponding resolution. For the three unnormalized MFs, C was estimated from only 10^{3} FFP6 simulations that proved to be sufficient for convergence. We compare the value obtained from the data to the χ^{2} obtained from those simulations, and report the probability of having a value of χ^{2} larger than the measured one, . We explore different resolutions represented by the parameter N_{side}, different methods of component separation (CommanderRuler, NILC, SEVEM, and SMICA) and different sky coverages.
First, the three unnormalized MFs (V_{k} as a function of ν, k = 0, 1, 2) are used to construct a test of the null hypothesis. The test assesses not only the primordial Gaussian hypothesis, but also whether the data is correctly represented by the simulations in terms of power spectrum, systematics and the lensing effect. A set of 17 thresholds in units of T/σ_{0} between −4 and + 4 in steps of 0.5 are considered. The comparison between the MFs of the data provided by the four component separation methods and those corresponding to each of the four sets of 10^{3} FFP6 simulations representing each method, for the standard U73 mask, are shown in Fig. 10. From that figure, a deviation at a level of about 2σ can be seen for the contour and genus curves at a resolution N_{side} = 512. The situation is very similar for the analyses performed at other resolutions, N_{side} = 1024, 256 and 128. Although the deviation is not particularly compelling, because of the correlations among neighbouring thresholds, it is worth mentioning that a possible explanation is the background of unresolved sources that has been detected in Planck data with the bispectrum estimators (see Planck Collaboration XXIV 2014). In order to understand the effect of unresolved sources on the MFs, we added the point source residuals derived from the FFP6 simulations as processed by the SEVEM algorithm to 100 realizations which were then analysed. We conclude that the background of unresolved sources may be responsible for at least part of the excess signal that is detected. The corresponding probabilities derived from the MF values for each of the four component separation methods and resolutions are given in Table 13. The full resolution maps have been degraded to the lower resolution ones following the procedure described in Sect. 2. All the cases considered are compatible with the null hypothesis.
Probabilities of obtaining values for the χ^{2} statistic of the Npoint functions shown in Fig. 6 for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck CMB maps with resolution parameter N_{side} = 64, estimated using the CR, NILC, SEVEM and SMICA methods.
Probabilities of obtaining values of the χ^{2} statistic of the Npoint functions shown in Figs. 7 and 8 for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for Planck CMB maps with resolution parameter N_{side} = 512. estimated using the CR, NILC, SEVEM and SMICA methods.
Probabilities of obtaining values of the χ^{2} statistic of the Npoint functions shown in Fig. 9 for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for Planck CMB maps with resolution parameter N_{side} = 2048 estimated using the CR, NILC, SEVEM, and SMICA methods.
Fig. 10 Difference of the data MFs (unnormalized) with respect to the average of the curves obtained with realistic Planck simulations for several cleaned maps. From top to bottom: Area, Contour, Genus. The errorbars represent the ±1σ (68% CL) dispersions around the mean obtained with simulations. 
In the second case, the four normalized MFs v_{k} = V_{k}/A_{k} (k = 0, 1, 2, 3) are used for the null hypothesis test. A set of 26 thresholds equally spaced between −3.5 and + 3.5, in units of T/σ_{0}, are considered. The normalization factor A_{k} is estimated directly from the maps, having previously computed the moments σ_{0} and σ_{1}. This normalization minimizes the dependence of the MFs on the power spectrum, thereby decreasing the cosmic variance and improving their sensitivity to deviations from Gaussianity. The resolutions considered in this case are N_{side} = 2048,1024,512,256, and 128. For the highest resolution N_{side} = 2048, the map is smoothed with a Gaussian smoothing kernel with a width θ_{FWHM} = 5′, in order to decrease the noise level. We use the standard U73 mask, inpainting the smallest point sources. The maps at lower resolution are constructed by the standard simple degrading process applied to the original map at N_{side} = 2048, and the corresponding masks are degraded following a conservative procedure such that any degraded pixel with a value <0.8 is set to zero (as explained in Sect. 2). The results of the analysis performed on the SMICA map at different resolutions are presented in Table 14.
Nondirectional Gaussianity tests using unnormalized MFs: as a function of sky resolution for different component separation methods.
Nondirectional Gaussianity tests using normalized MFs: dependence of on resolution.
The results of the analysis performed on the different component separation methods at the highest resolution (N_{side} = 2048) are presented in Table 15. The difference of the normalized MFs with respect to the expected values of the null hypothesis as a function of the threshold ν are shown in Fig. 11. A slight deviation in N_{cluster}(ν) is noticeable at thresholds ν ≈ 0, however it is not very compelling, since the values at neighboring thresholds are very correlated and this correlation is taken into account in the χ^{2} statistics.
Nondirectional Gaussianity tests using normalized MFs: dependence of on component separation method.
Finally, we analyse the dependence of the normalized MFs on the sky coverage. We use the standard U73 mask and then decrease the sky coverage by using CL65, CL48 and CL25 masks in combination with a special point source mask that is based on the U73 mask. The fraction of sky left unmasked in the combined masks is 62%, 46% and 23%, respectively. The point source mask is inpainted previously to the analysis. The curves obtained for the different sky coverages are presented in Fig. 12, for the SMICA method. Results of the χ^{2} analysis of the data as a function of sky coverage are compiled in Table 16. All the cases considered are compatible with the null hypothesis.
In summary, we find that the data are globally consistent with the primordial Gaussian hypothesis, and no strong deviation is found between the data and realistic simulations for both the unnormalized and normalized MFs. We would like to remark that a certain level of nonGaussianity is expected from lensing and, in particular, from the ISWlensing signal, thus it is important to compare the data to realistic lensed simulations.
4.5. Wavelet statistics
A broad range of wavelets have been used in the analysis of CMB data, but in this paper we specifically consider the spherical Mexican hat wavelet (SMHW, MartínezGonzález et al. 2002).
The SMHW is an example of a continuous, nonorthogonal wavelet. Given a signal on the sky, T(p), where p represents a given position/pixel wich is a function of the colatitude θ and longitude φ (also defined by the unit direction vector x), the SMHW coefficients at a given scale R, ω_{T}(R,p), are obtained by convolution: (23)where is the window function associated with the SMHW, ℓ_{max} is the maximum multipole allowed by the corresponding HEALPix pixelization, Y_{ℓm}(p) is the spherical harmonic basis, and t_{ℓm} are the spherical harmonic coefficients of the analysed map: (24)where dΩ = dθsinθdφ and ^{∗} denotes complex conjugation.
Several statistics can be computed from the wavelet coefficients map, in particular, the first moments: the dispersion σ_{R}, the skewness S_{R}, and the kurtosis K_{R} (as a function of scale R). It is interesting to notice that in the case of Gaussian temperature fluctuations the linear transformation involved in the determination of the wavelet coefficients (Eqs. (23) and (24)) guarantees that Gaussianity is preserved.
The study of the moments of the distribution of the CMB temperature fluctuations, as a function of the scale, is a standard approach to test the null hypothesis. We have performed a full resolution multiscale analysis of the four CMB clean maps and computed the quantities σ_{R}, S_{R}, and K_{R} from the SMHW wavelet coefficients at 18 scales, R = {2, 4, 7, 14, 25, 50, 75, 100, 150, 200, 250, 300, 400, 500, 600, 750, 900, 1050}, in arcminutes. These are compared to the standard Planck simulations.
As explained in Vielva et al. (2004), when computing the SMHW coefficients of a masked data set, artefacts are introduced close to the mask that degrade the performance of any null hypothesis tests. We therefore define a set of exclusion masks such that, at each scale, an extra region of the sky is excluded when performing any statistical test. The exclusion mask for a given scale R is defined as follows: we build an auxiliary mask by removing from the U73 mask all the features associated with compact objects, and degrade this auxiliary mask to N_{side} = 1024 (imposing a restrictive cut); a first temporary mask is obtained by extending the borders of this auxiliary mask by a distance of twice R; a second temporary mask is obtained, first, by convolving the auxiliary mask with the SMHW at that particular scale R and, second, by imposing that any pixel of that second temporary mask with an absolute value lower than 0.1 is masked, whereas the remaining ones are set to 1; the two temporary masks are multiplied to yield a single mask that is upgraded to N_{side} = 2048; and finally, the final exclusion mask is obtained by multiplying this mask by the parent U73 mask.
Fig. 11 Difference of the normalized MFs obtained from the data with respect to the expected values of the null hypothesis, for the different component separation methods. From left to right and top to bottom: area, Contour, Genus and N_{cluster}. The grey bands represent the ±1 and ±2σ dispersions around zero, based on realistic Planck simulations including lensing, for the CR method. 
The comparison of the four CMB maps with the corresponding simulations is summarized in Fig. 13. The three panels show (from left to right) the statistical significance of the standard deviation, the skewness and the kurtosis (as a function of the SMHW scale). The points represent the upper tail probability associated to a given statistic, i.e., the fraction of the simulations that present a value of a given statistic equal to or greater than the one obtained for the data. In fact, we define a modified upper tail probability: if an upper tail probability is larger than 0.5, then a new quantity is defined as 1 minus that upper tail probability. Hence, this modified definition of the upper tail probability is constrained between 10^{3} (the minimum value that can be imposed with 1000 simulations) and 0.5.
Fig. 12 Difference of the normalized MFs obtained from the data with respect to the expected value of the null hypothesis for several sky coverages. The SMICA map is considered. From left to right and top to bottom: area, Contour, Genus and N_{cluster}. The grey bands represent the 1 and 2σ dispersions around zero, based on realistic Planck simulations including lensing, for f_{sky} = 0.23. 
Overall, the agreement between the four CMB maps is quite good, showing that all of them provide a consistent estimation of the true CMB. However, several aspects need to be discussed. Let us clarify that the differences among the CMB methods for small modified upper tail probabilities are expected to be larger than for large modified upper tail probabilities. This is because a small modified upper tail probability is determined by a small number of simulations and, therefore, has a relatively large error bar. In other words, the tails of the distributions of the different statistics are quite sparsely sampled.
We will distinguish between the small (R ≲ 10′), intermediate (10′ ≲ R ≲ 500′), and large (R ≳ 500′) scale regimes. Let us focus on the three statistics independently. We will highlight the most important features and, afterwards, we will try to find an explanation for them.

On the smallest scales, the four CMB maps show a dispersion inSMHW coefficients significantly larger than seen in thesimulations. However on larger scales, the dispersion issystematically below the median of the simulations and, on scalesof R ≈ 5°, the modified upper tail probability is approximately 0.015.

Regarding the skewness, all four maps yield a value that is significantly lower (with a modified upper tail probability of around 0.004) than expected from the simulations in the small scale regime (except for the smallest one, where the deviation is around 0.07). The rest of the scales are fairly compatible with the null hypothesis.

The kurtosis is also smaller than expected in the small scale regime. Overall, the modified upper tail probability is about 0.03. At scales of around 300′, an anomalously large value (modified upper tail probability of approximately 0.01) is found.
These results are compatible with the values reported for WMAP data (Vielva et al. 2004; Cruz et al. 2005), over the scales common to both experiments (i.e., R> 10′). In particular, the large value of the kurtosis has been associated with the Cold Spot (Vielva et al. 2004). We will return to this topic specifically in Sect. 5.9. The low variance of the wavelet coefficients was previously seen in Vielva et al. (2004) and Wiaux et al. (2008). In addition, the low dispersion at scales above a few degrees is likely to be related to the low variance anomaly detected in WMAP (Monteserín et al. 2008; Cruz et al. 2011), that is also seen in the Planck data (see Sect. 4.1).
We have also studied the robustness of the results for different masking scenarios. In particular, we have investigated variations in the results when we adopt, as auxiliary masks to define the exclusion masks, the two CG70 and CG60 masks removing 30% and 40% of the sky, respectively. Note that the auxiliary masks obtained from the U73 mask already cut around 20% of the sky. The corresponding results for the SMICA map are presented in Fig. 14. The conclusions are similar for the other CMB maps. For the dispersion of the wavelet coefficients, we do not notice any significant change in the anomalously high value obtained for the SMICA map at the smallest scales However, some changes are observed at larger scales. In this regime, it seems that the most significant deviation occurs for the CG70 mask (modified upper tail probability of around 0.005), whereas similar (and slightly less significant) modified upper tail probabilities are obtained for both the U73 (modified upper tail probability of approximately 0.015) and the CG60 (modified upper tail probability of about 0.01) masks. A possible explanation for this behaviour would be that a less restrictive mask admits some residual contamination from Galactic foregrounds, thus increasing the dispersion of the wavelet coefficients, and artificially increasing their inconsistency with the null hypothesis. In principle, the larger the Galactic cut, the lower would be the dispersion of the wavelet coefficients (assuming that some residual contamination of the Galactic foregrounds is left) and, therefore, the smaller the upper tail probability. However, as we already said, the modified upper tail probability for the dispersion is higher for the CG60 mask than for the CG70 mask. This apparent contradiction could be resolved by accounting for the larger sampling variance for smaller areas, which would result in a lower significance for the anomaly.
The anomalous kurtosis at scales of R ≈ 300′ shows an overall stable modified upper tail probability of around 0.01–0.03. In the small scale regime, the differences are better defined: the smaller the mask, the more significant the deviation (characterized by the low value of the kurtosis). In particular, the modified upper tail probability associated with the CG60 mask is 0.001, around 0.009 for CG70, and approximately 0.03 for the U73 mask. A similar pattern is also observed for the skewness on these scales, although the three masks results in more similar upper tail probabilities, between around 0.001 and 0.007 (except for the smallest scale).
It is therefore clear that there is some inconsistency between the CMB data and the corresponding simulations. On intermediate scales, both the low dispersion and the high kurtosis could be related to previously known anomalies: the low variance and the Cold Spot. On the smallest scales, the three statistics report a low upper tail probability independently of the mask coverage – it is important to determine what this inconsistency is due to. Besides the possibility that it is an intrinsic cosmological signal, the nonGaussianity could be caused either by instrumental systematics or residual foreground contamination.
Nondirectional Gaussianity tests using normalized MFs: dependence of on sky coverage.
Fig. 13 Standard deviation (left), skewness (centre) and kurtosis (right) of the SMHW coefficients as a function of the wavelet scale R. Results are given for the four Planck CMB maps (green, CommanderRuler; lightblue, NILC; red, SEVEM; and yellow, SMICA). Modified upper tail probabilities (mUTP, see text for details) are obtained by comparing with 1000 simulations processed through the component separation pipelines. Squares represent modified upper tail probabilities that correspond to an actual upper tail probability above 0.5; diamonds represent upper tail probabilities below 0.5. 
Fig. 14 Standard deviation (left), skewness (centre) and kurtosis (right) of the SMHW coefficients as a function of the wavelet scale R. Results are given for the SMICA CMB map. Several masking scenarios are compared: red, CG60 mask (cutting out 40% of the sky); green, CG70 mask (cutting out 30% of the sky); and blue, U73 mask. The modified upper tail probabilities (mUTP) are defined in the text. 
In the former case, we have considered whether the origin of the signal could be related to properties of the noise that are inadequately modelled by the simulations. In particular, we have studied the statistical properties of the halfring halfdifference maps generated by the four component separation algorithms as proxies for the the noise present in the CMB maps. Although in detail there are some discrepancies between these noise estimates and the simulated ones, they are not compatible with the inconsistencies observed between the CMB map and simulations. Therefore, a systematic effect associated with the instrumental noise does not provide a satisfactory explanation for the smallscale deviations.
In the latter case, an obvious candidate is due to the contribution from residual unresolved point sources in the clean CMB maps. Although the brightest point sources are masked, and the component separation process itself can suppress the amplitude of the unresolved background of point sources, some signal will remain. Indeed, in Planck Collaboration XXIV (2014) it has been determined that the bispectrum of this contribution is clearly detected in the four CMB Planck maps, at a significance in excess of 4σ. In addition, the dispersion of the wavelet coefficients is higher than expected, which is also compatible with the presence of an additional signal. We therefore consider this as the most likely nonCMB explanation for the observed signal.
4.6. Bispectrum
The CMB bispectrum is the three point correlator of the a_{ℓm} coefficients: . In this paper, we focus on the bispectrum reconstruction as a blind test of nonGaussianity. Therefore, we assume we are seeking a nontrivial bispectrum that has arisen through a physical process which is statistically isotropic, that is, we can employ the angleaveraged bispectrum B_{ℓ1ℓ2ℓ3} and the reduced bispectrum b_{ℓ1ℓ2ℓ3}. Note that b_{ℓ1ℓ2ℓ3} is defined on a tetrahedral domain of multipole triples { ℓ_{1}ℓ_{2}ℓ_{3} } satisfying both a triangle condition and a limit given by the maximum resolution ℓ_{max} of the experiment. A much more extensive introduction to the bispectrum can be found in Planck Collaboration XXIV (2014).
Modal, wavelet and binned bispectrum estimators filter the CMB map with separable basis functions (25)to find the corresponding modal coefficients. For appropriately orthonormalized basis functions Q_{ijk}(ℓ_{1},ℓ_{2},ℓ_{3}), these coefficients can be used to reconstruct the CMB bispectrum through a signaltonoise weighted expansion (Planck Collaboration XXIV 2014). This reconstruction method has been extensively validated, showing the accurate recovery of CMB bispectra from nonGaussian simulated maps, and it has been applied to the WMAP seven year data to reconstruct the full 3D CMB bispectrum (Fergusson et al. 2010). To quantify whether or not there is a modelindependent deviation from Gaussianity, we can consider the total integrated bispectrum. By summing over all multipoles, we use the integrated nonlinearity parameter defined in Planck Collaboration XXIV (2014). For ideal Gaussian CMB maps, the quantity should obey a χ^{2}distribution with a mean given by the number of degrees of freedom (the modes) μ = n_{max} and with a variance σ^{2} = 2n_{max}. Assuming that the threepoint correlator is the leading nonGaussian contribution, then provides a blind test for the presence of any integrated CMB bispectrum (once the expected twopoint term is subtracted). We note that this is less sensitive than targeted searches for particular bispectrum shapes.
First, we discuss reconstructions from the modal estimator which has passed successfully through the full suite of nonGaussian bispectrum validation tests (for further details about bispectrum estimators, see Planck Collaboration XXIV 2014). We have applied this to the Planck temperature maps for the foregroundseparation techniques NILC SEVEM and SMICA, using two alternative sets of hybrid basis functions in order to crosscheck results and identify particular signal (Planck Collaboration XXIV 2014). As with all the other bispectrum analyses based on spherical harmonic coefficients, we used the U73 mask to which we applied inpainting. Together with the foreground separated maps, noise simulations were provided which were used to calibrate the estimator’s linear correction term and to determine the variance.
The modal coefficients extracted from the PlanckNILC, SEVEM, and SMICA maps show remarkable consistency between the different maps, with shape crosscorrelations better than 96% and the overall amplitudes to within 7% agreement (see corresponding figure in Planck Collaboration XXIV 2014). This demonstrates that the indendent foreground separation techniques do not appear to be introducing spurious nonGaussianity. The corresponding quantity shows consistency with the null hypothesis (Planck Collaboration XXIV 2014). Using the modal expansion, we have reconstructed the full 3D Planck bispectrum (Planck Collaboration XXIV 2014) for SMICA but also NILC and SEVEM; the reconstructions are visually indistinguishable.
Fig. 15 Comparison between the WMAP sevenyear bispectrum signal (left) (Fergusson et al. 2010) and the lowℓ signal of Planck (right) reconstructed from the SMICA foregroundseparated map (in both cases using polynomial modes). The same basic patterns are observed in both bispectra, including an apparent central “oscillatory” feature. 
In Fig. 15, we show a comparison of the ℓ < 500Planck bispectrum signal and that reconstructed from the WMAP sevenyear data (Fergusson et al. 2010). Here for consistency we show the Planck signal using polynomial modes. The Planck signal pattern correlates well with the WMAP bispectrum obtained previously, despite the different domains used for the modal analysis of the two different experiments.
Similarly to the modal bispectrum, a wavelet decomposition can be used to reconstruct the bispectrum. Here we use the continuous, nonorthogonal SMHW (MartínezGonzález et al. 2002). Cubic moments q_{ijk} are defined in terms of the SMHW coefficients for three different angular scales R_{i}, R_{j}, R_{k} (Curto et al. 2009b,a, 2010, 2011a, 2012). See Planck Collaboration XXIV (2014) for more details on the implementation of this wavelet estimator for Planck data analysis. Considering the covariance matrix of the q_{ijk} moments, C ≡ ⟨ qq^{T} ⟩, and its eigenvector decomposition, C = RDR^{T}, with R the eigenvector matrix and D the eigenvalue matrix, a new set of quantities y ≡ D^{− 1 / 2}R^{T}q is defined. Considering the decorrelation produced by the convolution of the SMHW on the temperature anisotropies and applying the central limit theorem to the averages defined in q_{ijk}, then these quantities are expected to have a nearly Gaussian distribution. Therefore, the y quantities are nearly multinormal and satisfy ⟨ yy^{T} ⟩ = I and ⟨ y ⟩ = 0 (Curto et al. 2011a).
We have computed this reconstruction using the Planck data and compared with the null hypothesis (Gaussian Planck simulations). The considered data map is the resulting map after foreground cleaning based on different cleaning procedures: CommanderRuler, NILC, SEVEM, and SMICA. The mask used is the U73 one (contrary to the modal reconstruction, no inpainting of the point sources is made in this case). In Fig. 16 the y statistics corresponding to the Planck data are plotted and compared with the 3σ errorbars obtained with Planck Gaussian simulations. From the list of different q_{ijk} statistics corresponding to the 16 angular scales described in Planck Collaboration XXIV (2014), there are 11, 4, 3, 3 statistics with  y_{i}  ≥ 3 (corresponding to CommanderRuler, NILC, SEVEM, and SMICA, respectively). The errorbars are obtained with Planck simulations for each type of component separation cleaned map. The errorbars of the y_{i} statistics for low indices i are associated to large scales where the q statistics have a less Gaussianlike shape. The y statistics are combined into a χ^{2} test after a principal component analysis with a threshold of 10^{12} (Curto et al. 2011a) and compared with the χ^{2} statistics obtained from Planck Gaussian simulations for each type of component separation method (see Table 17). The χ^{2} statistic corresponding to the data is compatible with the values obtained from Gaussian simulations according to the cumulative probability , as can be seen in Table 17. Therefore the wavelet bispectrum reconstruction does not detect a significant amplitude of bispectrum in the considered data maps. Details on the constraints on the amplitude of different bispectrum shapes are presented in Planck Collaboration XXIV (2014).
Fig. 16 Wavelet bispectrum reconstruction y_{i} statistics for the foreground cleaned Planck data map. Considered data map: combined map cleaned with CR, NILC, SEVEM and SMICA. The solid yellow lines represent the 3σ errorbars for SMICA (similar errorbars are obtained for CR, NILC, and SEVEM maps). 
χ^{2} statistics based on the wavelet bispectrum reconstruction y_{i} statistics for the foreground cleaned data map.
5. Intriguing inconsistencies – WMAP anomalies revisited
In the previous section, we established that the Planck data show little evidence for nonGaussianity beyond that expected due to the ISWlensing effect (which is accounted for directly by simulations) and contributions from residual unresolved point sources. The exceptions are on largeangular scales, where features consistent with various anomalies previously seen in the WMAP data have been observed. In this section, we explicitly consider several of the most important anomalies detected in the WMAP data, namely the quadrupoleoctopole alignment (Sect. 5.1), the low variance (Sect. 5.2), the hemispherical asymmetry (Sect. 5.3), phase correlations (Sect. 5.4), power asymmetry (Sect. 5.5), dipolar modulation (Sect. 5.6), generalized power modulation (Sect. 5.7), parity asymmetry (Sect. 5.8), and the Cold Spot (Sect. 5.9). Each of these anomalies may represent different violations of the fundamental properties of isotropy and/or Gaussianity in the CMB data that are assumed for the estimation of the CMB power spectrum.
There is an ongoing debate about the significance of these anomalies in the literature. A critical issue relates to the role of a posteriori choices – whether interesting features in the data bias the choice of statistical tests, or if arbitrary choices in the subsequent data analysis enhance the significance of the features. Indeed, the WMAP team (Bennett et al. 2011) contends that the anomalies are significantly overinterpreted due to such selections, whilst other authors claim highly significant and robust detections. Therefore, care must be taken to address the issue of significance, since our analyses necessarily follow up tests of the previous WMAP investigations. However, a careful and fair statistical treatment can enable us to study possible links among the anomalies and to search for a physical interpretation.
5.1. Mode alignment
Fig. 17 Upper: Wienerfiltered SMICA CMB sky (temperature range ±400 μK). Middle: derived quadrupole (temperature range ±35 μK). Lower: derived octopole (temperature range ±35 μK). The plus and star symbols indicate the axes of the quadrupole and octopole, respectively, around which the angular momentum dispersion is maximized. The diamond symbols correspond to the quadrupole axes after correction for the kinematic quadrupole. 
Tegmark et al. (2003) first reported a significant alignment between the orientation of the quadrupole and the octopole in the WMAP first year temperature data. We study this quadrupoleoctopole alignment in the Planck data using the maximization of the angular momentum dispersion as described in de OliveiraCosta et al. (2004). Specifically, we determine the orientation of the multipoles by finding the axis n around which the angular momentum dispersion (26)is maximized. Here, a_{ℓm}(n) denotes the spherical harmonic coefficients of the CMB map in a rotated coordinate system with its zaxis in the ndirection. This definition of the multipoleorientation has been devised for planar multipoles and is simply the direction perpendicular to the plane in which most of the power of the multipole lies. It is thus intuitive and easy to use. Note that the value of the statistic in Eq. (26) is the same for −n as for n, i.e., the multipole orientation is only defined up to a sign.
An alternative method, based on the multipole vector decomposition of the data (Copi et al. 2004; Schwarz et al. 2004; Bielewicz et al. 2005; Bielewicz & Riazuelo 2009), has also been used to verify the robustness of the results presented here, and excellent consistency is found.
Residual foregrounds (mostly on the Galactic plane) present in the four Planck CMB estimates could influence the reconstruction of the loworder multipoles. However, when a mask is applied, the resulting modecoupling can also affect the reconstruction of the lowℓ multipoles. Since our understanding of the relationship between the size of the mask and the extent to which residual foregrounds and modecoupling can affect the result is incomplete, we therefore utilize Wienerfiltered maps computed from the data to which the U73 mask has been applied. Specifically, we utilize the same implementation of the Wiener filter as used in Planck Collaboration XXIV (2014), i.e., a “messenger” method, as described by Elsner & Wandelt (2013). A direct inversion method for masked data (Efstathiou 2004; Bielewicz et al. 2004, 2013) is a possible alternative, but the Wienerfiltered maps result in a significantly smaller uncertainty in the reconstructed orientation of the multipoles.
We search for the preferred orientation by explicitly rotating the CMB map such that the zaxis pierces the centres of each of the low resolution pixels defined at N_{side} = 16, and then subsequently refine the search by testing the N_{side} = 2048 pixels that lie within the preferred lower resolution pixel. The angular resolution for the orientation of the multipoles is thus given by the distance between the pixel centres of the N_{side} = 2048 map, which is of order 1.94′. Figure 17 shows the Wienerfiltered SMICA CMB sky, with the corresponding reconstruction of the quadrupole and octopole moments. The reconstructed orientations are quite robust with respect to the component separation method used for reconstructing the CMB. The significance of the alignment between the quadrupole and the octopole is assessed from the scalar product of their orientations, compared to values derived from the standard set of 1000 MC simulations. The orientation, the angular distance, the scalarproduct between quadrupole and octopole, and the probability of at least as strong an alignment to occur in an isotropic universe, are summarized in Table 18 for each CMB map.
We find that, depending on the component separation method, the quadrupole and octopole orientations are misaligned by an amount between 9° and 13°. Such a range of results is reasonable given the observed differences, typically of order 5μK in amplitude at high Galactic latitudes, between the CMB estimates, as shown in Figs. 5 and 6 of Planck Collaboration XII (2014). This then implies that it may not be possible to determine the alignment with a precision of better than a few degrees, due to uncertainties related to component separation.
Note that the misalignment determined from the Planck data is larger than the 3° reported^{3} recently by Bennett et al. (2013) for the 9year WMAP ILC map. However, when considering the plausible range of errors introduced by the covariance between the CMB and foregrounds and their impact on the derived ILC map, Bennett et al. (2013) cite a median misalignment angle of 6°.
Finally, as noted by Tegmark et al. (2003), the observed quadrupole includes a contribution due to our motion relative to the CMB rest frame. This kinematic quadrupole (hereafter KQ) has a welldefined morphology and frequency dependence, as described in Kamionkowski & Knox (2003). The latter is of particular importance for the HFI channels. Schwarz et al. (2004) have demonstrated that correcting for the KQ contribution can be relevant for the evaluation of the significance of the alignment between the cosmological quadrupole and octopole. We therefore apply mapbased corrections to the NILC, SEVEM, and SMICA maps, taking into account both the weighted contributions of each Planck frequency to the final CMB maps and the corresponding KQ scaling factors. Such corrections are possible for these maps, since, unlike for CommanderRuler, the frequency weights used to generate them do not vary locally on the scales of interest. As seen in Table 18, after KQcorrection the three estimates of the quadrupoleoctopole alignment are more consistent, with the misalignment angle decreasing to approximately 8°, and the significance increasing to 99%.
5.2. Variance, skewness, and kurtosis anomalies
A low value for the variance on the CMB sky was previously observed in the WMAP data by Monteserín et al. (2008) and Cruz et al. (2011), and confirmed for Planck in Sect. 4.1. Furthermore, the effect has also been seen in the wavelet analysis of Sect. 4.5, where the variance of the SMHW coefficients is low at scales between 400 and 600 arcmin (Fig. 13). In addition, anomalous behaviour was also observed for the skewness and kurtosis in low resolution maps at N_{side} = 16. Here, we reassess these results and determine their robustness to masking and data selection. The former will allow us to determine whether a particular region is causing the anomalous behaviour, whilst the latter can establish whether foreground residuals could be responsible.
Table 19 and Fig. 18 present the results for the variance, skewness and kurtosis determined from the four CMB maps with the U73, CL58 and CL37 masks applied. Results are also computed for data within the ecliptic hemispheres surviving the U73 mask. The variance is low in all cases, with only small differences in significance observed for the different maps. Interestingly, the low variance seems to be localized in the northern ecliptic hemisphere. Conversely, anomalous values for the skewness and kurtosis only are apparent for the southern ecliptic hemisphere.
Since these results might be indicative of the presence of Galactic foreground residuals near the Galactic plane, we analyse the frequency dependence of the statistics as summarized in Table 20 and Fig. 19. The variance shows little frequency dependence for the considered masks and regions, whereas the skewness and kurtosis show a moderate frequency dependence when the U73 mask is applied, as also seen for the Npdf analysis in Sect. 4.2. Cruz et al. (2011) found that a small region of the sky localized to both the ecliptic and Galactic south and near to the Galactic plane (their socalled “gp10” region) exhibited particularly high variance. Thus, since the skewness is negative, we consider a prominent cold spot at b = − 8°, l = 32°, partially masked by the Galactic plane. However, when masking the seven coldest pixels of the spot, the significance of the skewness and kurtosis only drops slightly, with lower tail probabilities of approximately 0.03 and 0.93 respectively. If the whole “gp10” region (f_{sky} = 7%) is masked, the skewness and kurtosis drop drastically and have lower tail probabilities of approximately 0.30 and 0.50, respectively, whereas the variance is highly significant since none of the 1000 simulations has a variance below that of the data. In order to check the possible leakage of Galactic contamination due to the Gaussian smoothing applied to the low resolution data, we repeated our calculations for the Wiener filtered maps used in Sect. 5.1, but found little variation to the existing results. Therefore, it is unlikely that any leakage impacts the estimators.
Lower tail probablity for the variance, skewness, and kurtosis at N_{side} = 16, using different masks.
The incompatibility of the observed variance with simulations based on a cosmological model that has been determined from the same data set might appear puzzling at first, but can be understood as follows. The mapbased variance is dominated by contributions from large angular scales on the sky, whilst the cosmological parameter fits are relatively insensitive to these loworder ℓmodes, and are instead largely dominated by scales corresponding to ℓ > 50. Thus, the bestfit spectrum in the context of a 6parameter ΛCDM model can have a mismatch with the data on these scales, so that the corresponding simulations will not adequately capture the dearth of power at lowℓ. The results presented here do indeed imply that the largeangular scale power is low relative to the fiducial sky model. In fact, when subtracting the quadrupole and octupole from both the data and simulations outside the U73 mask, the results are more consistent. In this case, the lower tail probabilities for the variance, skewness and kurtosis are 0.192, 0.637 and 0.792, respectively. This result was already found in Cruz et al. (2011). It is then plausible that the low multipole alignment could have the same cause as the anomalies considered here. Interestingly, Sarkar et al. (2011) have demonstrated that the lowvariance and mode alignment is very unlikely to have a common origin in Gaussian random and statistically isotropic CMB anisotropies. However, when subtracting the quadrupole and octupole outside the CL58 mask, the lower tail probability for the low variance is 0.036, which remains rather low. The connection with the very low power in the ecliptic northern hemisphere also remains to be explored.
Fig. 18 Variance, skewness and kurtosis at N_{side} = 16, for the U73 mask, CL58, CL37, ecliptic North, and ecliptic South (from top to bottom). The different lines represent the four component separation methods CR (green), NILC (blue), SEVEM (red), and SMICA (orange). 
Frequency dependence of the lower tail probablity for the variance skewness and kurtosis at N_{side} = 16, using different masks.
Fig. 19 Variance, skewness and kurtosis at N_{side} = 16, for the U73 mask, CL58, CL37, ecliptic North, and ecliptic South (from top to bottom). The different lines represent the four considered frequencies, namely 70 GHz (green), 100 GHz (blue), 143 GHz (red), and 217 GHz (orange). 
5.3. Hemispherical asymmetry
In Eriksen et al. (2004a) and Hansen et al. (2004) it was discovered that the angular power spectrum of the first year WMAP data, when estimated locally at different positions on the sphere, appears not to be isotropic. In particular, the power spectrum calculated for a hemisphere centred at (l,b) = (237°,−20°) (in Galactic longitude and latitude) was larger than when calculated in the opposite hemisphere over the multipole range ℓ = 2–40. Simultaneously, Park (2004) also presented evidence for the existence of such hemispherical asymmetry – in which a particular statistical measure is considered to change discontinuously between two hemispheres on the sky – with the application of MFs to the WMAP data. Since the preferred direction of Eriksen et al. (2004a) lies close to the ecliptic plane, it was also demonstrated that the largeangular scale Npoint correlation functions showed a difference in behaviour when computed on ecliptic hemispheres. Many studies have subsequently been undertaken focusing on hemispheres in the ecliptic coordinate system, with Schwarz et al. (2004) particularly emphasizing the connection. Hemispherical asymmetry has also been seen with other measures of nonGaussianity (Eriksen et al. 2004c, 2005; Räth et al. 2007a).
Here, we repeat the analysis of Eriksen et al. (2005) on the Planck component separated data, smoothed and then downgraded to N_{side} = 64, as described in Sect. 2. However, in this section the Npoint correlation functions are not averaged over the full sky and depend on a choice of specific direction, thus, they constitute tools to study statistical isotropy rather than nonGaussianity (Ferreira & Magueijo 1997). The latter study was already presented in Sect. 4.3, where it was found that the results for the low resolution maps are the most deviant relative to the MC simulations based on the Planck fiducial model.
The Npoint correlation functions computed on the northern and southern hemispheres determined in the ecliptic coordinate frame and using the U73 mask are shown in Fig. 20. The correlation functions for the four Planck maps are very consistent, and the observed behaviour is in agreement with that seen in the WMAP data (Eriksen et al. 2004a). Specifically, the northern hemisphere correlation functions are relatively featureless (both the 3 and 4point functions lie very close to zero), whereas the southern hemisphere functions exhibit a level of structure that is in good agreement with the confidence regions computed from the Gaussian simulations.
Fig. 20 The 2point (upper left), pseudocollapsed (upper right), equilateral 3point (lower left), and rhombic 4point (lower right) correlation functions (N_{side} = 64). Correlation functions are shown for the analysis performed on northern (blue) and southern (red) hemispheres determined in the Ecliptic coordinate frame. The shaded dark and light grey bands indicate the 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 
Fig. 21 The 2point (upper left), pseudocollapsed (upper right), equilateral 3point (lower left), and rhombic 4point (lower right) correlation functions for the SMICA map (N_{side} = 64). Correlation functions are shown for the analysis performed on negative (blue) and positive (red) hemispheres determined in the Ecliptic, Galactic, Doppler boost (DB), dipole modulation (DM), and power asymmetry (PA) coordinate frames. The shaded dark and light grey bands indicate the 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 
The probabilities of obtaining a value for the χ^{2} statistic for the Planck fiducial ΛCDM model at least as large as the observed values for different CMB componentseparated maps are presented in Table 21. The probabilities for the 3point and 4point functions in the northern ecliptic hemisphere are especially large, and in the case of the pseudocollapsed configuration all simulations yielded a larger value of the χ^{2} than observed. This suggests that the tendency to superGaussianity in the maps downgraded to N_{side} = 64, as noted previously in Sect. 4.3, is associated with the northern hemisphere. Nominally, this value is even more remarkable than that found with the WMAP data (Eriksen et al. 2004a), although to interpret it correctly one has to keep in mind that the analysis presented here is an example of an a posteriori statistic. Specific choices have been made about the smoothing scale used for downgrading the data, and, in particular, for the reference direction that defines the hemispheres. This results in a tendency to overestimate the significance of the results. Nevetheless, the observed properties of the Planck data are consistent with a remarkable lack of power in a direction towards the north ecliptic pole, consistent with the simpler onepoint statistics presented in Sect. 5.2.
We have also carried out the corresponding calculations in Galactic coordinates, as well as in reference frames preferred by other analyses, such as the power asymmetry (PA, reference frame for ℓ_{max} = 1500, see Sect. 5.5), the dipole modulation (DM, see Sect. 5.6), and the Doppler boost (DB, see Planck Collaboration XXVII 2014). The corresponding Npoint correlation functions are shown in Fig. 21. Note that the positive hemispheres for the ecliptic and Galactic reference frames correspond to the southern hemispheres in other tables/figures. Interestingly, the largest asymmetry is seen for the hemispheres determined in ecliptic coordinates.
In Table 22 we show the corresponding probabilities for the Npoint functions of the SMICA map in these reference frames. Whilst the largest asymmetry is seen for ecliptic coordinates, a substantial asymmetry is present also for Galactic coordinate hemispheres and, to a lesser extent, for the DM direction. The asymmetry for the latter can be explained by the fact that the DM direction is relatively close to the ecliptic pole. For the PA and DB directions we do not find any significant asymmetry. In the latter case, this is not surprising, since the Doppler boost effect is significant only for small angular scales, while the map used in our analysis was smoothed and downgraded to much lower resolution.
Since the largest asymmetry is seen in ecliptic coordinates, we might conjecture an association with the peculiar configuration of the quadrupole and octopole (although the Planck results from Sect. 5.1 do not favour as significant an alignment as with the WMAP data). In order to study the influence of the lowest order multipoles on the hemispherical asymmetry, we therefore reanalysed the SMICA map corrected for multipoles of order ℓ ≤ 3 and ℓ ≤ 5 and present the results in Table 23. It should be apparent that, in general, the asymmetry is less significant after removing the lowest order multipoles. However, in a few cases, such as for the 4point function computed in the ecliptic and DM frames, we do not see such behaviour. Furthermore, in all cases, with the exception of the DB frame, we observe a tendency to superGaussianity for the negative hemisphere. This presumably indicates that the hemispherical asymmetry is not only related to the properties of the very lowest order multipoles, but also extends beyond this.
Probabilities of obtaining values of the χ^{2} statistic for the Planck fiducial model at least as large as the observed values of the statistic for the Planck maps with resolution parameter N_{side} = 64, estimated using the CR, NILC, SEVEM and SMICA methods.
5.4. Phase correlations
The statistical properties of a Gaussian random field on the sphere can be fully described by its twopoint correlation function (or power spectrum) and should exhibit Fourier phases that are independent and identically distributed (i.i.d.) and follow a uniform distribution in the interval [− π,π]. Thus, demonstrating the existence of Fourier phase correlations in CMB maps could indicate the presence of nonGaussianities in the primordial density fluctuations.
Probabilities of obtaining values of the χ^{2} statistic for the Planck fiducial model at least as large as the observed values for the SMICA map at N_{side} = 64.
One could search for phase correlations using different approaches. However, without there being a specific model for the origin of any such correlations, one would like to test for scaledependent nonGaussianity and deviations from isotropy in a completely modelindependent (“blind”) way. One such approach is offered through the method of “scaling indices” and “surrogates” (which we explain below).
Previous studies using these methods, based on the WMAP three, five and sevenyear data (Räth et al. 2009, 2011; Rossmanith et al. 2012; Modest et al. 2013), showed significant evidence for intrinsic phase correlations at low ℓ values in the CMB. The signal was demonstrated to be robust with respect to the WMAP data release, to the component separation methods and to the selected test statistics. In this section we apply these methods to the Planck component separated data sets.
The scaling index method represents one way to estimate the local scaling properties of a point set in an arbitrary ddimensional embedding space. The technique provides the possibility of revealing local structural characteristics of a given point distribution. A number of analyses have used scaling indices to test the Gaussian nature and statistical isotropy of the CMB as represented by the WMAP data (Räth et al. 2007a, 2009; Rossmanith et al. 2009a; Räth et al. 2011).
Probabilities of obtaining values of the χ^{2} statistic for the Planck fiducial model at least as large as the observed values of the statistic for the PlanckSMICA map with resolution parameter N_{side} = 64 and removed multipoles of order ℓ ≤ 3 (columns denoted as “ℓ > 3”) or order ℓ ≤ 5 (columns denoted as “ℓ > 5”) .
In general, the method is a mapping that calculates, for each member p_{i},i = 1,...,N_{pix} of a point set P, a single value that depends on the spatial position of p_{i} relative to the group of other points in its neighbourhood, in which the point under consideration is embedded. A threedimensional point set P is generated for twodimensional spherical CMBdata by transforming the temperature values T(θ_{i},φ_{i}) of each pixel to a radial jitter around a sphere of radius R at the position of the pixel centre (θ_{i},φ_{i}). To obtain scaling indices the local weighted cumulative point distribution is calculated: (27)Here r describes the scaling range, while s_{r} and d denote a shaping function and a distance measure, respectively. The scaling index α(p_{i},r) is then defined as the logarithmic derivative of ρ(p_{i},r) with respect to r: (28)Using a quadratic Gaussian shaping function s_{r}(x) = exp( − (x/r)^{2}) and an isotropic Euclidian norm d(p_{i},p_{j}) = ∥ p_{i} − p_{j} ∥ as distance measure, one obtains the following analytic formula for the scaling indices: (29)where we use the abbreviation d_{ij} ≡ d(p_{i},p_{j}).
As seen from the derivation of the expression for α, the scaling index analysis can be considered as a structural decomposition of the point set under study. Points in welllocalized, clusterlike structures lead to scaling indices of α ≈ 0, filamentary structures have α ≈ 1, and sheetlike structures α ≈ 2. A uniform, threedimensional distribution of points would result in α ≈ 3. One can schematically think of the scaling indices applied to CMB data as estimating the roughness of the last scattering surface on different scales.
As should be clear from Eq. (29), the calculation of scaling indices depends on the scale parameter r. Ten scaling range parameters r_{k} = 0.05,0.1,...,0.5, k = 1,2,...,10, in the notation of Räth et al. (2007a), are used in this analysis. In order to calculate scaling indices on large scales, as in previous studies, we couple the rjitter a to r_{k} via a = 0.5r_{k}. The mean ⟨ α(r_{k}) ⟩ and standard deviation σ_{α(rk)}, derived from the full sky and from a set of 768 rotated hemispheres, are used to test for nonGaussianity and deviations from statistical isotropy.
In order to quantify the significance of the scaling index results, and focus the study on the phase properties of the observed CMB sky, we utilize the method of “surrogate maps” (Räth et al. 2009). This surrogate approach suppresses the sensitivity of the null tests to the assumed fiducial power spectrum, which is particularly pertinent, given the potential mismatch of the Planck data to the fiducial spectrum on large angular scales (Planck Collaboration XV 2014).
The basic idea of the method of surrogates is to generate for a given data set an ensemble of socalled surrogate data, which mimic some properties of the original data being defined by a null hypothesis, while all other properties are subject to randomization. Analysing both the original and surrogate data sets with test statistics, which are sensitive to the complement of the null hypothesis, then reveals whether the original data were consistent with the null hypothesis or not. In our case, the null hypothesis is that the CMB map is Gaussian and thus fully described by the power spectrum and free of any phase correlations. The test statistics being employed are those that are sensitive to higherorder correlations (i.e., correlations beyond the twopoint correlation function) in real space. Alternatively, one can think of introducing Fourier phase statistics, but those methodologies are currently less developed.
However, the Gaussianity of the temperature distribution and the randomness of the set of Fourier phases of the map to be studied are a necessary prerequisite for the application of the surrogategenerating algorithm. Therefore the following preprocessing steps are applied to generate a zero order surrogate map. First, the maps are remapped onto a Gaussian distribution in a rankordered way. Then we ensure the randomness of the set of Fourier phases by making a rankordered remapping of the phases onto a set of uniformly distributed ones. These two preprocessing steps only have marginal influence on the maps. Now the set of surrogates to be used for assessing the statistical properties of the data sets can be generated by shuffling the phases in the space of the spherical harmonics, while exactly preserving the modulus of the a_{ℓm}. Moreover, by introducing a twostep shuffling scheme for previously specified ℓranges, a scaledependent analysis is made possible. It is worth noting that while in all surrogate maps the modulus of each a_{ℓm} is exactly preserved, null tests involving a comparison to an assumed fiducial power spectrum only preserve the C_{ℓ} values, which are average values of the a_{ℓm} when summed over m. Thus, the linear properties of the surrogate maps are more tightly constrained, and specifically kept constant, than they would be in tests involving simulated maps generated on the basis of the C_{ℓ}s.
Socalled first and secondorder surrogate maps are then obtained as follows. We initially generate a firstorder surrogate map, in which any phase correlations for the scales that are not of interest are randomized by shuffling the phases φ_{ℓm} for ℓ ∉ Δl = [ℓ_{min},ℓ_{max}] ,0 <m ≤ ℓ. In a second step, N (=1000) realizations of secondorder surrogate maps are generated from the firstorder surrogate map, in which the remaining phases φ_{ℓm} with ℓ ∈ Δℓ, 0 <m ≤ ℓ are shuffled, while the previously randomized phases for the other scales are preserved. The generation of surrogates is always performed using the maps with the highest resolution, i.e., N_{side} = 2048. Given the evidence for anomalies on the largest angular scales, and to ensure consistency with the previous WMAP analyses, we perform dedicated scaledependent tests for the scales defined by Δℓ = [2,20].
The exclusion of the phases belonging to the zonal modes (m = 0) from the shuffling raises the question as to whether the results depend on the choice of coordinate system, since these modes are only defined with respect to the selected coordinate system. This issue is particularly important in the specific case of the CMB when analysing data in the Galactic coordinate system, because possible residuals in and around the Galactic plane are likely to be reflected mostly in these modes. To test the dependence of the results on the coordinate system, we repeated the surrogate analysis with N = 100 secondorder surrogates for each member of a set of 48 rotated CMB maps, where the rotations, defined by a HEALPix N_{side} = 2 map, such that the pixel centres indicate the orientation of the zaxis of the rotated frame, uniformly cover the unit sphere.
Since the methodology in its simplest form requires orthonormality of the set of basis functions Y_{ℓm}, we apply the method to the full sky foregroundcleaned maps as obtained after component separation with CommanderRuler, NILC, SEVEM and SMICA. For the selected ℓinterval Δℓ = [2,20], the generation of the firstorder surrogate map removes the phase signature of the small scale residuals in the data and this can be interpreted as a form of inpainting procedure for small masked patches in the Galactic plane. The differences between the first and secondorder surrogates are quantified by the σnormalized deviation (30)with Y = ⟨ α(r_{k}) ⟩ ,σ_{α(rk)},χ^{2}. Here, χ^{2} represents either a diagonal combination of the mean ⟨ α(r_{k}) ⟩ and standard deviation σ_{α(rk)} at a certain scale r_{k}, or for the full scaleindependent χ^{2} statistics: (31)where the test statistics to be combined are contained in the vector M and C is obtained by crosscorrelating the elements of M. With the mean and the standard deviation as input for M we obtain and statistics with M^{T} = ( ⟨ α(r_{1}) ⟩ ,..., ⟨ α(r_{10}) ⟩) and M^{T} = (σ_{α(r1)},...,σ_{α(r10)}), respectively.
Fig. 22 Deviations S( ⟨ α ⟩) of the rotated hemispheres derived from the mean of the scaling indices for the scale r_{5} determined from the SMICA map. Upper panel: result for the calculations with respect to the Galactic coordinate system. Lower panel: result obtained after averaging over 48 rotated coordinate systems. 
Figure 22 shows the deviation S( ⟨ α(r_{5} ⟩) values for the set of rotated hemispheres for the SMICA map. Each pixel of the full sky map with a HEALPix resolution of N_{side} = 8 specifies one of the 768 Svalues for a rotated hemisphere, where the pixel position indicates the orientation of the zaxis of the rotated coordinate system. We find a large number of hemispheres for which S is larger than 3, indicating a significant difference between the first and secondorder surrogate maps. Therefore, one has to conclude that the lowℓ phases are correlated. In addition, one can clearly see that the deviation S goes in opposite directions for a number of pairs of opposite hemispheres. This indicates the presence of pronounced anisotropy in the CMB. We further find that the Smap derived for the Galactic coordinate system and the Smap from the average of rotated coordinate systems are nearly indistinguishable. This shows that the zonal modes are not special and that the choice of Galactic coordinates as reference frame does not induce any bias. Therefore, all of the results that follow are only shown for Galactic coordinates. We also find very similar yet less significant patterns of asymmetry when considering the standard deviation of the scaling indices σ_{α}, which leads to a significant detection of anomalies when combining both quantities into a χ^{2}statistic (see Fig. 23). The results are consistent for CommanderRuler, NILC, SEVEM and SMICA.
Fig. 23 Deviations S(χ^{2}) of the rotated hemispheres derived from a combination of the mean and the standard deviation of the scaling indices for the scale r_{5} determined from the SMICA map. 
Fig. 24 Deviations  S(r)  for the SMICA map as a function of the scale parameter r for the full sky (black) and positive (red) and negative (blue) rotated hemisphere. The reference frame for defining the positive and negative hemisphere is chosen such that the difference ΔS = S_{pos} − S_{neg} becomes maximal for ⟨ α(r_{5}) ⟩. The plus signs denote the results for the mean ⟨ α(r_{k}) ⟩, the starsigns represent the standard deviation σ_{α(rk)}. The dashed (dotted) line indicates the 1 (3) σ significance level. 
Deviations S and empirical probabilities p for the scaleindependent χ^{2}statistics derived from the CR, NILC, SEVEM, and SMICA maps.
In Fig. 24, the deviations S(Y) are displayed for the mean and standard deviation. We only show the results for the SMICA map. The other three maps yield very similar results. For all four maps the values for S( ⟨ α ⟩) extend beyond 3 for r = 0.2–0.25 when rotated hemispheres are considered separately. Since the effect in the separate hemispheres goes in opposing directions, no signal is observed for the full sky. The results for the scaleindependent χ^{2} statistics are summarized in Table 24. The results suggest a significant detection of both nonGaussianity and anisotropy in the Planck data, consistent with those obtained previously with WMAP data (for comparison see Modest et al. 2013).
Fig. 25 Probability density P( ⟨ α(r_{5}) ⟩) for the secondorder surrogates using the SMICA map for the positive (red) and negative (blue) rotated hemispheres. The vertical lines show the respective values for ⟨ α(r_{5}) ⟩ for the firstorder surrogates. 
Figure 25 presents results that allow us to attempt to elucidate the origin of the detected asymmetry and its relationship to other anomalies. Considering the firstorder surrogate maps, in which the lowℓ phases are unchanged by the shuffling procedure, there is a clear shift between the mean of the scaling indices for the positive (red) and negative (blue) hemispheres. If we now examine the probability density of the mean of the scaling indices for the positive and negative rotated hemispheres, as determined from the set of secondorder surrogates, again a systematic shift is observed. Since the phases are uncorrelated in the secondorder surrogate maps, the effect can only be attributed to properties of the distribution of a_{lm}, i.e., to power asymmetries. However, it should be apparent that the shift in mean values is in the opposite sense to the firstorder surrogates, implying that an additional pattern of asymmetry is induced solely by correlations in the phases.
Finally, it should be noted that certain choices made in the analysis, e.g., the selection of specific values for r_{k} and a parameters, together with the decision to make an analysis on hemispheres using the method of surrogates, constitute a posteriori selections and can lead to an overestimate of the significance of the results. Nevertheless, the claimed asymmetry is certainly consistent with results presented in the following sections. The analysis involving surrogate maps reveals that there are phase correlations at low ℓ.
5.5. Power asymmetry
In their analysis of the 5year WMAP data, Hansen et al. (2009) specifically searched for dipolar power asymmetry on the sky. In particular, a simple test was performed in which the power spectrum on discs was computed and binned into independent blocks of 100 multipoles from ℓ = 2 to ℓ = 600, then each block fitted for a dipolar asymmetry in the power distribution. The six ℓ ranges considered showed evidence of a consistent dipole direction. The significance of this was shown by the fact that not a single realization in a set of 10 000 simulations showed a similarly strong asymmetry. A further extension of the analysis introduced a model selection procedure taking into account the statistical penalty for introducing an asymmetric model with additional parameters (direction of asymmetry, amplitude of asymmetry and asymmetric multipole ranges). Even in this case, the asymmetry was found to be highly significant for the whole range ℓ = 2 to 600.
Unfortunately, this procedure is highly CPUintensive. Given the increased sensitivity and angular resolution of the Planck data, we have therefore chosen to focus on the simpler test, following the modified approach of Axelsson et al. (2013) as applied to the 9year WMAP data. Specifically, we estimate the power spectrum amplitude on 12 nonoverlapping patches of the sky in ℓbins of 16 multipoles each, and fit a dipole to the spatial distribution of the amplitude for each ℓbin. This allows us to probe further into a previously unexplored highℓ range and to answer at least in part any a posteriori criticisms of the study. Since the analysis is powerspectrum based, the halfring data sets for the different CMB estimators are used^{4}. The approach is as follows.

1.
We apodize the U73 Galactic, and point source mask bysmoothing it twice with a 30′ Gaussian beamfollowed by a 10′ smoothing that is repeated threetimes. The apodized mask is multiplied with the original mask between each of the smoothing operations, thus all pixels which arezero in the U73 mask remain so in order to continueto fully eliminate foreground residuals. In addition, a smoothtransition between masked and unmasked regions is created.The power spectrum estimation method has been shown to takeinto account the mode couplings introduced by sharp edges in themask. Nevertheless, we have adopted an apodized mask, sincethe error bars at high ℓ decrease slightly due to the reduced couplings. The effective sky fraction of the U73 maskfalls from 73% to 69% due to this procedure.

2.
The MASTER technique (Hivon et al. 2002) is used to estimate local cross power spectra between the two halfring maps on 12 nonoverlapping patches on the sky. A given patch corresponds to one of the twelve base pixels of an N_{side} = 1HEALPix map, masked as appropriate by the apodized U73 Galactic and point source mask. This results in usable sky coverages of approximately 7% for each of the eight patches surrounding the poles, and of order 3% for each of the four patches aligned with the Galactic plane. This approach differs from that of the earlier analysis of Hansen et al. (2009) in which 3072 highly overlapping discs were used, but correlations were ignored. The revised method results in weaker correlations between the nonoverlapping patches.

3.
Following Hansen et al. (2009) the 12 local spectra are estimated in bins of 16 multipoles. The highest multipole included in the analysis is ℓ_{max} = 1500, since both noise and highℓ foregrounds become important at higher multipoles (Planck Collaboration XII 2014; Planck Collaboration XV 2014; Planck Collaboration XVI 2014). For some comparisons, it is useful to follow Hansen et al. (2009) and further combine these 16multipole bins into blocks containing approximately 100 multipoles. Although no single block contains exactly 100 multipoles, they will be referred to as such. However, our main results are based on the 16multipole blocks, which offer improved statistics over the more coarsely binned data for our statistic of choice.

4.
Each 16multipole block now has an associated local power map at N_{side} = 1. The HEALPix routine remove_dipole is then used with inverse variance weighting (with the weight for each region determined from the variance of the local power spectra derived from simulations) to obtain the dipole direction of the power distribution for each map.

5.
The alignment of the dipole directions between the different multipole blocks is then used to construct a measure of the power spectrum asymmetry. Despite the maskinduced correlations between adjacent multipoles, the power spectra estimated in 16multipole blocks should be highly independent and the dipole directions determined for a Gaussian, isotropic CMB realization should be random, as confirmed by simulation.

6.
The same procedure is then applied to a set of 500 simulated maps of the CMB and noise to derive significance.
In order to assess the significance of the asymmetry, it is necessary to compare the clustering of the dipole directions evaluated for the different scales to that observed in simulated maps. For this purpose, we compute the mean separation angle, θ_{mean}, between all possible pairwise combinations of the 16multipole dipole directions up to a given ℓ_{max}. The expectation value for this statistic is 90° when determined from the dipole directions computed on the full sky for an isotropic, Gaussian realization of the CMB. We have shown that this remains true for an analysis when applying the apodized U73 mask. We determine θ_{mean} for the data as well as for all 500 simulated maps. The fraction of simulations with a higher value for θ_{mean} than the data defines a pvalue for the measurement. One advantage of this directional analysis is that it focuses on a central issue for tests of deviation from isotropy – whether there is a preferred direction. Moreover, as noted by Bunn & Scott (2000), the CMB may exhibit a pattern that cannot be identified from the power spectrum, but which would indicate some nontrivial largescale structure. In the context of our analysis this means that the amplitudes of the power dipoles may not be discernible, in a frequentist sense, from the distribution of values expected for an ensemble of Gaussian, isotropic CMB realizations. However, evidence for the close correlation and alignment of directions on different angular scales presents a clear signature of power asymmetry, since in the standard model, these directions should all be independent random variables. Our claims about power asymmetry in this section are based solely on this statistic.
Fig. 26 Dipole directions for 100multipole bins of the local power spectrum distribution from ℓ = 2 to 1500 in the SMICA map with the apodized U73 mask applied. We also show the total direction for ℓ_{max} = 600 determined from WMAP9, as well as the preferred dipolar modulation axis (labelled as lowℓ) derived in Sect. 5.6. 
In Fig. 26 we show the dipole directions of the 15 lowest 100multipole bins for the SMICA map. The preferred direction for WMAP9 over the range ℓ = 2 to 600 is also indicated, as is the preferred lowℓ modulation direction determined in Sect. 5.6. Similar behaviour is seen for all of the Planck foregroundcleaned maps.
Summary of the bestfit power dipole directions on the sky up to ℓ_{max} = 600/1500, for each of the four component separation methods, together with the 143 GHz map cleaned using the SEVEM method and denoted SEVEM143.
In Table 25 we present a summary of the power asymmetry results from the Planck data processed by all four Planck foreground cleaning methods – CommanderRuler, NILC, SEVEM, and SMICA – computed on the apodized U73 mask. For comparison, we also include the latest WMAP 9year result^{5} computed with their KQ85 mask (which has a usable sky fraction of 74.8%). It should be apparent that the clustering seen visually in Fig. 26 is both significant and consistent for all foreground subtraction methods.
However, Doppler boosting, due to our motion with respect to the CMB rest frame, is expected to induce both a dipolar modulation of the temperature anisotropies, and an aberration that corresponds to a change in the apparent arrival directions of the CMB photons (Challinor & van Leeuwen 2002). Both of these effects are aligned with the CMB dipole. In a companion paper (Planck Collaboration XXVII 2014), a statistically significant detection of such an effect at small angular scales is reported. From an inspection of Fig. 26 it should be apparent that whilst the asymmetry direction on large and intermediate angular scales are in general tightly clustered around the same direction as found for WMAP, weak indications of this second preferred power asymmetry direction might be seen. However, we have determined that our statistic of choice is unable to cleanly separate a Doppler dipolar modulation from a power asymmetry imprinted on the data.
In order to investigate the power asymmetry further, it is necessary to analyse data that have been corrected for the Doppler modulation effect. As shown in Planck Collaboration XXVII (2014), the modulation amplitude is frequency dependent. It is therefore difficult to correct for the effect in the component separated maps used in this paper, since they are constructed from frequency maps to which both scaledependent and spatially varying weighting schemes are applied. The SEVEM approach, however, adopts a single weight per frequency. In what follows, we utilize a version of the 143 GHz data cleaned using this method (which we will refer to as SEVEM143), with the weights summarized in Table C.1. of Planck Collaboration XII (2014). The Doppler modulation can then be removed using the prescription and weights derived in Planck Collaboration XXVII (2014). We will refer to this process as deboosting, even though a correction for aberration is not performed. Such a deaberration is not needed since the FFP6 simulations used in this analysis contain the aberration effect, but not the Doppler modulation of the CMB anisotropies. However, in order to match the noise properties of the demodulated data, the noise component of each simulation is separately demodulated.
In Fig. 27 we show the dipole directions of the 100multipole bins after the deboosting correction. As compared to Fig. 26, the dipole directions have now moved away from the Doppler dipole direction and cluster more around the WMAP9 and low multipole asymmetry directions. The mean dipole direction for ℓ = 2 to 1500 is now (l,b) = (218°,−21°), displaced by about 25° compared to the nondeboosted directions in Table 25. Note that the mean asymmetry directions for the ℓ = 2 to 600 range for the deboosted SEVEM143 and WMAP9 data differ by about 20°, most likely due to the differences in the masks used and noise properties.
Figure 28 presents the pvalues (the fraction of simulations with a larger mean separation angle than that determined from the real sky) as a function of ℓ_{max}. It is certainly the case that power asymmetry is observed to ℓ_{max} ~ 600, even for deboosted data. Cursory inspection may lead to the conclusion that asymmetry at a significance level of 98–99% persists to at least ℓ_{max} = 1300. However, this is a consequence of the cumulative nature of the statistic. Investigation of θ_{mean} for individual multipole ranges (see below) does not indicate the presence of asymmetry above ℓ_{max} ~ 600. Similar conclusions were reached by Flender & Hotchkiss (2013) and Notari et al. (2013), although they only looked at the power ratio which is not directly related to the test performed here.
Fig. 27 Dipole directions for 100multipole bins of the local power spectrum distribution from ℓ = 2 to 1500 in the deboosted SEVEM 143 GHz map (SEVEM143DB) with the apodized U73 mask applied. We also show the direction for ℓ_{max} = 600 determined from WMAP9, as well as the preferred dipolar modulation axis (labelled as lowℓ) derived in Sect. 5.6. 
We have also attempted to assess the degree to which the significance of the asymmetry depends on a specific choice for ℓ_{max} by implementing several global statistics. Values for these are determined as follows.

1.
The angle θ_{mean} is calculated as a function of ℓ_{max} both for the data and thesimulations.

2.
For a given ℓ_{max}, the fraction of simulations (out of 500) with a larger mean separation angle than for the real data is defined as the pvalue for that ℓ_{max}.

3.
This procedure is then repeated, treating each of the 500 simulations in turn as if it were the real data, and computing a corresponding pvalue from the fraction of the other 499 simulations that gives a larger θ_{mean} than the selected simulation. Examples of such pvalues as a function of ℓ_{max} for the data are shown in Figs. 28 and 29.

4.
We define three global statistics in the following manner.

“Number of p = 100%”: for certain values of ℓ_{max}, each of the 500 simulations isfound to have a larger mean separation angle than is determinedfor the real data, implying a corresponding pvalue of 100%. Wethen evaluate the number of ℓ_{max} values for which this is the case. Eachof the simulations is then treated in turn as if it were the real data,and effectively the same quantity is determined by comparison tothe other 499 simulations. The fraction of all simulations thatyield a smaller “number of p = 100%” than the real data is then used as aglobal statistic, shown in the first row in Table 26.

“Number of p> 99%”: for certain values of ℓ_{max}, more than 99% of the simulations exhibit larger mean separation angles than the data, implying a pvalue of >99%. We count the number of ℓ_{max} values for which this applies. The same quantity is then evaluated for each simulation, treated in turn as if it were the real data (as above). The fraction of all simulations that yield a smaller “number of p> 99%” than the real data is used as a second global statistic.

“Mean pvalue”: the mean pvalue for the data is computed as an average over the set of pvalues determined for each ℓ_{max} (as in 2 above). A corresponding quantity is evaluated for each simulation, treated in turn as if it were the real data. The fraction of simulations with a lower “mean pvalue” than for the data then forms the third global statistic.

We provide results for all three global statistics, since they reflect different aspects of the significance of the asymmetry. The fact that the global statistics indicate similar significance levels is evidence for the robustness of the results. We determine their values both for all ℓ_{max} available (ℓ_{max} = 2–1500) and after restricting ℓ_{max} to 2–600, since we do not claim any asymmetry beyond this value. The results are shown in Table 26. In addition, significance levels are also provided for the case when only multipoles ℓ > 100 are included in the analysis. In this case, strong evidence (>99%) is still found for asymmetry from all global measures.
We can also study the range of scales over which the asymmetry persists using these global statistics and a modified definition of θ_{mean}. Here, this is computed as the mean of the angular separations determined between all pairs of dipole directions where one direction corresponds to the range ℓ_{min}<ℓ <ℓ_{lim} and the second direction to ℓ_{lim}<ℓ <ℓ_{max}. Figure 29 shows an example of the pvalues corresponding to this definition of θ_{mean} for ℓ_{lim} = 500. The high pvalues indicated in the figure for multipoles in the range ℓ = 500–600 show that the corresponding dipole directions are strongly correlated with the dipoles in the range ℓ = 2–500. For most of the calculations, ℓ_{min} = 2 and ℓ_{max} = 1500. Table 27 presents the global significance levels for ℓ_{lim} = 100,200,300,400,500, and 600. We find that the dipole directions at high and low multipoles are strongly correlated up to ℓ = 600, although the significance is weakening for ℓ > 500. Indeed, when we consider multipoles ℓ > 600 and ℓ_{lim}> 600, we find no evidence for asymmetry. The table also provides similar significance levels for the case ℓ_{min} = 100.
Although the global significance values in Table 27 show that the asymmetry is almost insignificant for ℓ_{lim} = 500, we must take into account the fact that the global statistics look at the entire range 2–1500. If the asymmetry disappears at, for example, ℓ = 600, then the statistics would be too poor in the interval ℓ = 500–600 for this to show up in the global significance measure for ℓ_{lim} = 500. For this reason we show in Fig. 29 the plot of pvalues as a function of ℓ_{max} for ℓ_{lim} = 500, showing that there is a strong signal all the way to ℓ_{max} = 600, after which it disappears. The continued evidence for asymmetry excludes the possibility that the results in Fig. 28 are dominated by the known lowℓ signal, and supports the claim for asymmetry persisting over the range ℓ = 100 to 600.
An alternative approach to studying power asymmetry is to determine the ratio of the local power spectra computed in two opposing directions on the sky (e.g., Eriksen et al. 2004a). Here, we consider such a ratio defined for the two hemispheres centred on the positive and negative poles of the power dipole fitted over a given ℓ range. A statistic can then be defined through the fractional power ratio as follows: (32)where corresponds to the power spectrum computed for the hemisphere centred on the positive pole, and to the spectrum in the antepodal direction. This can be compared to an ensemble of isotropic and Gaussian simulations to determine whether significant excess power is observed. Figure 30 presents this quantity, binned into blocks of 100 multipoles, for the hemispheres centred on the preferred dipole direction determined for the SEVEM 143 GHz map over the ℓ range 2 to 600, both before and after deboosting. It should be apparent that, although the ratio lies systematically above zero for ℓ < 600, only a few bins lie significantly outside the range of values generated from simulations. The most significant bins are those centred on ℓ = 50 and 150. Note that the observed values are not directly comparable to the explicit dipolar modulation fits in Sects. 5.6 and 5.7. The ratio on other scales is of smaller amplitude and lower significance. We reemphasize that the claims of significant asymmetry presented in this section are based on the alignment of the power distribution as a function of angular scale, not on the corresponding amplitudes, nor on the ratio of power in the antipodal hemispheres. The ratio of power in the antipodal hemispheres is shown here for only illustrative purposes.
Fig. 28 Derived pvalues as a function of ℓ_{max}, for the SEVEM 143 GHz foregroundcleaned map, determined both before and after deboosting. The pvalues are computed using 500 FFP6 simulated maps. 
Fig. 29 Derived pvalues as a function of ℓ_{max} for the deboosted SEVEM 143 GHz foregroundcleaned map. The pvalues in this plot are based on the the mean of the angular separations determined between all pairs of dipole directions where one direction falls in the range [ℓ_{lim} = 500,ℓ_{max}] and the second direction in the range [2,ℓ_{lim} = 500]. The significance is computed using 500 FFP6 simulated maps. 
Significance of the asymmetry using several global significance measures.
Significance of the correlations of dipole directions between high and low multipoles.
Fig. 30 Fractional power ratio, ΔC_{ℓ}/C_{ℓ}, from antipodal sky regions, computed from the SEVEM 143 GHz map before and after deboosting along the mean dipole direction for ℓ = 2 to 600. All spectra are evaluated on hemispheres using an apodized mask. The grey lines show the same quantity evaluated for each of the 500 FFP6 simulations along their respective asymmetry axes. The green band shows the 68% confidence region from these simulations. 
In summary, we have presented evidence for power asymmetry in the Planck data. At high ℓ, this is expected, since a dipolar modulation of the temperature anisotropy due to Doppler effects has been predicted, and subsequently detected, as detailed in a companion paper (Planck Collaboration XXVII 2014). When this is taken into account, significant power asymmetry can be claimed up to ℓ ~ 600.
5.6. Dipole modulation
In Sect. 5.5 it was shown that power asymmetry is observed over a large range of angular scales in the Planck data, with a fairly consistent preferred axis. No explicit parametric model was assumed in the analysis. In this section, however, we only consider large angular scales and revisit the phenomenological model due to Gordon et al. (2005), who proposed that the power asymmetry could be described in terms of a multiplicative dipole modulation model of the form d = (1 + Ap·n)s_{iso} + n ≡ Ms_{iso} + n, where A is the dipole amplitude, p is the dipole direction, n denotes instrumental noise, and s_{iso} is an underlying isotropic CMB field. Both s_{iso} and n are assumed to be Gaussian random fields with covariance matrices S and N, respectively. Since s_{iso} is assumed to be isotropic, its covariance may be fully specified by some angular power spectrum C_{ℓ,iso}.
In the following we present the results from a direct likelihood analysis of this model, similar to those described by Eriksen et al. (2007a) and Hoftuft et al. (2009) for the 3 and 5year WMAP data, respectively. Since this method requires matrix inversions and determinant evaluations, the computational expense scales as , and it is therefore only feasible at low resolution. Specifically, we consider maps downgraded to a HEALPix pixel resolution of N_{pix} = 32, smoothed to angular resolutions ranging from 5° to 10°, ensuring sufficient bandwidth limitation at this pixelization. All four Planck CMB map solutions are included in the analysis; however, note that the Galactic plane is handled differently in each of the four approaches. Specifically: for the Commander map the region inside the corresponding analysis mask has been replaced with a Gaussian constrained realization, eliminating the possibility that bright Galactic foreground residuals might leak outside the mask during degradation (Planck Collaboration XV 2014); for SMICA and NILC a smaller region is replaced with Wiener filtered data; and for SEVEM no special precautions are taken.
Summary of dipole modulation likelihood results as a function of scale for all four Planck CMB solutions.
After degrading each map to the appropriate resolution, we add random uniform Gaussian noise of 1 μK rms to each pixel to regularize the covariance matrix. All pixels inside the U73 mask are excluded, and we adopt the difference maps between the raw Planck LFI 30 GHz and HFI 353 GHz maps and the SMICA CMB solution as two foreground templates, tracing low and highfrequency foregrounds, respecively. We marginalize over these Galactic foreground templates, f, as well as four monopole and dipole templates, by adding corresponding terms of the form αff^{T} to the total data covariance matrix, where α is set to a numerically large value.
Before writing down the likelihood for A and p, a choice has to be made for the power spectrum,C_{ℓ,iso}. We follow Eriksen et al. (2007a), and adopt a simple twoparameter model of the form C_{ℓ,iso} = q^{(}ℓ/ℓ_{pivot}^{)}^{n}C_{ℓ,fid}, where the fiducial spectrum, C_{ℓ,fid}, is the bestfit Planck spectrum, and q and n describe an amplitude scaling and spectral tilt with respect to this. The full model therefore includes five free parameters, namely three dipole parameters and two power spectrum parameters. Introducing the two parameters q and n addresses the known issue that the lowℓ power spectrum is low by about 2–2.5σ compared with the overall bestfit ΛCDM spectrum (Planck Collaboration XV 2014; Planck Collaboration XVI 2014). Ignoring this creates a tension with the underlying isotropic model that results in the analysis measuring a combination of both asymmetry and power mismatch.
Taking advantage of the fact that both the signal and noise are assumed Gaussian, the exact likelhiood may be written down in a convenient closed form: (33)This expression is the basis of all calculations presented in the rest of this section.
Fig. 31 Marginalized dipole modulation amplitude (top), power spectrum amplitude (middle) and power spectrum tilt (bottom) probability distributions as a function of smoothing scale, shown for the Commander CMB solution. 
Fig. 32 Consistency between component separation algorithms as measured by the dipole modulation likelihood. The top panel shows the marginal power spectrum amplitude for the 5° smoothing scale, the middle panel shows the dipole modulation amplitude, and the bottom panel shows the preferred dipole directions. The coloured area indicates the 95% confidence region for the Commander solution, while the dots shows the maximumposterior directions for the other maps. 
Due to the high computational expense associated with these evaluations, we do not compute the full joint fiveparameter model in this analysis, only conditionals of it. However, we iterate once in a Gibbssampling like approach, by maximizing each conditional to obtain an approximation to the full maximumlikelihood solution. That is, we first map out the dipole likelihood for the 5° FWHM case, fixing the power spectrum at the fiducial spectrum, ℒ(A,p  q = 1,n = 0), and locate the maximumlikelihood dipole parameters. Then we map out the corresponding power spectrum conditional, ℒ(q,n  A_{mℓ},p_{mℓ}). Finally, we update the dipole likelihood with these power spectrum parameters, and evaluate the final results. Note that the power spectrum and dipole modulation parameters are only weakly correlated, and this procedure is therefore close to optimal. Furthermore, the approach is also conservative, in the sense that it will always underestimate the significance of the dipole modulation model; the derived maximumlikelihood value will always lie slightly below the true maximumlikelihood point.
The results from these calculations are summarized in Table 28, listing results for all four Planck CMB maps at angular scales between 5 and 10° FWHM. For easy reference, we also list the results from the corresponding 3 and 5year WMAP analyses (Eriksen et al. 2007a; Hoftuft et al. 2009). Note that the former was performed at a HEALPix resolution of N_{side} = 16 and the latter at an angular resolution of 4.5° FWHM.
Figure 31 shows marginals for A, q and n, as derived from the Commander CMB solution for all smoothing scales. At least two interesting points can be seen here. First, while there is clearly significant scatter in the derived dipole modulation amplitude for different smoothing scales, as originally pointed out by Hanson & Lewis (2009), all curves appear to be consistent with a single value of A ~ 0.07. No other single value fits all scales equally well. Second, it is interesting to note that the lowℓ power spectrum derived here is consistent, but not without some tension, with the fiducial spectrum, (q,n) = (1,0), around 1.5–2σ. In particular, there appears to be a slight trend toward a steeper positive spectral index as more weight is put on the larger scales, a result already noted by COBEDMR. The same conclusion is reached using the lowℓPlanck likelihood, as described in Planck Collaboration XV (2014).
In Fig. 32 we compare the results from all four CMB solutions for the 5° FWHM smoothing scale. Clearly the results are consistent, despite the use of different algorithms and different treatments of the Galactic plane, demonstrating robustness with respect to the details of the analysis methods. Further, we also note that these results are consistent with those derived from the 5year WMAP ILC map by Eriksen et al. (2007a), demonstrating robustness across experiments. On the other hand, it is notable that a higher dipole amplitude was found for the 3year WMAP ILC map using a larger mask at 9° FWHM.
Fig. 33 Loglikelihood difference between the bestfit dipole modulation model and the fiducial isotropic model as a function of smoothing scale. Horizontal dashed lines indicate 1, 2, and 3σ thresholds. 
In Fig. 33 we show the loglikelihood difference between the derived maximumlikelihood point and the isotropic model, A = 0, as a function of smoothing scale. The power spectrum parameters are kept fixed at the bestfit values for both points, leaving three additional parameters for the dipole model. The dashed horizontal lines indicate the 1, 2, and 3σ confidence regions for three degrees of freedom. As has been noted previously in the literature, these significance levels vary with smoothing scale. Taken at face value, the results presented here are suggestive but clearly not decisive, resulting in an unchanged situation with respect to earlier reports. This is of course not unexpected, given that WMAP is already strongly cosmic variance limited at these angular scales.
The critical question is whether the trend seen at smaller angular scales in Fig. 33 continues, or if the apparent likelihood peak at 5° FWHM happens to be a local maximum. Hanson & Lewis (2009), and later Bennett et al. (2011), address this question through a computationally cheaper quadratic estimator, allowing them to extend a similar analysis to small scales. In doing so, they claim that the apparent likelihood peak at 5° is indeed a local maximum, and the evidence for the modulation model falls off when more data are included.
In this respect, it should be noted that the dipole modulation model was originally proposed by Gordon et al. (2005) as a simple phenomenological characterization of the more general power asymmetry. In particular, it assumes that the modulation amplitude, A, is equally strong on all scales. From both the results presented by Hanson & Lewis (2009) and Bennett et al. (2011) and qualitatively shown in Fig. 30, this appears not to hold, as the fractional hemispherical power difference is clearly smaller at ℓ > 300 than at ℓ < 100. Furthermore, an alternative power modulation analysis in Planck Collaboration XXIV (2014) finds that the asymmetry on these scales has an amplitude of order 1%. The data, therefore, are not consistent with a simple constantamplitude dipole modulation of the power. On the other hand, the preferred direction derived from the current lowℓ analysis is remarkably consistent with the highℓ direction derived in Sect. 5.5. A proper modulation model may therefore need additional spatial structure beyond that proposed by Gordon et al. (2005), as already suggested by Hoftuft et al. (2009) and Moss et al. (2011).
Finally, Doppler boosting as presented in Planck Collaboration XXVII (2014) creates a signature similar to that observed here. However, the present effect is clearly distinct from this, both because the magnitudes of the two effects are very different – A_{boost} ~ 0.002 versus A_{asym} ~ 0.07 – and because the two preferred directions are different – (l,b)_{boost} = (264°,48°) versus (l,b)_{asym} = (227°,−15°). We have analysed the 143 and 217 GHz frequency maps, cleaned from foregrounds using the SEVEM algorithm and deboosted using the bestfit parameters from Planck Collaboration XXVII (2014), and consistently find (very slightly) higher significance levels for these data. The bestfit axis moves by approximately 5° and 2° at 143 and 217 GHz, respectively. The results presented here, therefore, are conservative with respect to boosting corrections.
5.7. Generalized modulation
In this section, we study a generalization of the dipolar modulation field analysed in Sect. 5.6 using the Bipolar Spherical Harmonic (BipoSH) formalism. For a statistically isotropy sky, the spherical harmonic space twopoint correlation matrix is diagonal, and, given by the angular power spectrum C_{ℓ}. The BipoSH representation provides a natural, mathematically complete, generalization of the angular power spectrum that captures statistical isotropy violations via coefficients that are a completely equivalent representation of the spherical harmonic correlation matrix, (34)This relationship combines the offdiagonal spherical harmonic correlations into a bipolar multipole L,M – analogous to the total angular momentum addition of states. The CMB angular power spectrum corresponds to the L = 0 BipoSH coefficients .
A simple model that results in the violation of statistical isotropy arises from the modulation of the of the CMB sky, (35)where T(n) represents the modulated CMB sky, T_{0}(n) is the underlying (statistically isotropic) random CMB sky and M(n) is a fixed, zeromean, dimensionless, modulation field. The modulation signal, if any, is expected to be weak and allows quadratic terms in M to be neglected. The BipoSH coefficients for the modulated CMB field (L> 0) are then given by (36)Here corresponds to the BipoSH coefficients of the unknown, but statistically isotropic, unmodulated CMB field, m_{LM} are the spherical harmonic coefficients of the modulating field (L> 0), and C_{ℓ} is the bestfit CMB angular power spectrum. The statistically isotropic nature of the unmodulated CMB sky implies that the expectation values of vanish for L> 0, leading to the estimator for the modulation field harmonics, (37)denoted by the overhat (Hanson & Lewis 2009). The weights for a minimum variance estimate for the modulation field correspond to (38)where N is a normalization chosen such that . The BipoSH representation further enables an estimate of the modulation field to be made over specific angular scales by windowing regions in multipole space in the sum over multipoles ℓ_{1},ℓ_{2} in Eq. (37). This additional information could be very useful in identifying the origin of the statistical isotropy violation, which could be either cosmological or due to systematic artefacts (see Hajian & Souradeep 2003; Hajian & Souradeep 2006).
First, we limit our analysis to the four low resolution N_{side} = 32 CMB maps used in Sect. 5.6 and reconstruct the modulation maps for each of them at the same low resolution. The U73 mask is applied to the reconstructed modulation maps before computing m_{LM}. The pseudopower m_{L} is corrected for the mask applied to the modulation maps. Specifically, for the case of dipole modulation, the pseudopower m_{L} is related to the dipole amplitude by .
Fig. 34 Significance of the modulation power, L(L + 1)m_{L}/ 2π, at bipolar multipoles L. The modulation spectra obtained from the four component separation maps (CR, NILC, SEVEM, and SMICA) are consistent with each other. Dipole (L = 1) modulation power is detected in all the spectra at a significance ranging from 2.9 to 3.7σ. The solid black lines denote the 3σ significance thresholds. There is no significant power detected at higher multipole of the modulation field, 1 <L ≤ 32. 
Amplitude (A) and direction of the dipole modulation in Galactic coordinates.
A dipole modulation (L = 1) signal is detected at 3σ significance in all the maps, as shown in Fig. 34. The amplitude and direction of the dipole modulation match those obtained via a likelihood analysis in Sect. 5.6. The BipoSH representation of modulation confirms the dipole modulation signal found in the lowresolution map. Since this approach enables the reconstruction of any general small amplitude modulation field, the BipoSH representation places constraints on the power in the modulation field at all higher (bipolar) multipoles allowed by the resolution of the CMB maps.
Fig. 35 Measured dipole modulation (L = 1) power in CMB multipole bins. This uses the CMB multipole dependence of the BipoSH (modulation) power L(L + 1)m_{L}/ 2π, which can be dissected into bins in ℓspace. We establish that significant power in the dipole modulation is limited to ℓ = 2–64 and does not extend to the higher CMB multipoles considered. The vertical grid lines denote the CMB multipole ℓbins. 
We then extend the analysis to higher resolution using maps at N_{side} = 256 for Commander and N_{side} = 2048 for NILC, SEVEM, and SMICA, in order to study the above effect in more detail. We repeat the analysis on these higher resolution maps using the U73 mask. Contrary to our expectations based on a scaleindependent (i.e., no ℓdependence) model, the significance of the dipole does not increase in the high resolution maps. We then subdivide the ℓrange up to ℓ_{max} = 384 into uniform bins of size Δℓ = 64. As seen in Fig. 35, we recover the dipole modulation at over 3 σ significance only for the lowest bin (ℓ = 2–64). This is consistent with the results in Sect. 5.6 and the BipoSH analysis on the corresponding low resolution maps shown in Fig. 34. However, the amplitude of the dipole is consistent with zero within 3 σ for all of the higher ℓbins considered. This suggests that the simple modulation model in Eq. (35) is inadequate and should minimally allow for the amplitude, A(ℓ), of the dipole to depend on CMB multipole, ℓ. Although this may appear to be a more complex model, it does not necessarily lack motivation. It is readily conceivable that physical mechanisms that cause a dipolar modulation of the random CMB sky would be scaledependent and possibly significant only at low wavenumbers. More importantly, such a dipole modulation has also been noted in low resolution WMAP data (Eriksen et al. 2007a; Hoftuft et al. 2009). More recently, Bennett et al. (2011) also comment (without being quantitative) that the effect is present in the WMAP maps, but limited to low ℓ and conclude that the ℓdependence rules out a simple modulation explanation. The fact that two independent experiments find this intriguing statistical isotropy violation points to a noninstrumental origin.
The search for modulation power recovered from higher multipoles of the CMB maps yields a null result, as seen in Fig. 35. However, due to our motion with respect to the CMB rest frame, it is expected that the observed CMB map is statistically anisotropic, and this has been demonstrated in Planck Collaboration XXVII (2014). To understand why a Doppler boost induced anisotropy is not detected by the modulation estimator, we first implement an equivalent description in terms of the Doppler boost induced BipoSH coefficients, (39)where β_{1M} = ^{∫}dnY_{1M}(n)β·n, β = v/c denotes the peculiar velocity of our local rest frame with respect to the CMB, and b_{ν} is the frequency dependent boost factor, as discussed in more detail in Planck Collaboration XXVII (2014). Since we perform our analysis on componentseparated maps, we use the same approximation, b_{ν} = 2.5, as advocated there.
A minimum variance estimator for β_{1M} is given by the expression in Eq. (37), with the shape function replaced by the corresponding Doppler boost term given in Eq. (39).
Fig. 36 Dipolar (L = 1) power measured from different CMB multipole bins, as denoted by the vertical grid lines. Here the Doppler estimator is used to reconstruct the Doppler field () from SMICA. Significant dipolar power is limited to ℓ = 2–64 and does not extend to the higher CMB multipoles considered. 
Fig. 37 Reconstruction noise power associated with the modulation (blue) and Doppler (red) estimators as a function of the maximum CMB multipole used. The solid lines correspond to the case of noisefree CMB maps, while the dashed lines correspond to the case where the characteristic Planck instrument noise is taken into account. The red dashdotted horizontal line denotes the expected signal level of the boost due to our local motion, while the blue dashdotted horizontal line denotes the enhanced boost amplitude as seen by the modulation estimator. The blue dotted line denotes the effective Doppler boost signal as seen by the modulation estimator after accounting for the effects of aberration induced by the boost. Also note that the reconstruction errors shown in red are derived for the total dipole power, rather than the amplitude in a particular direction, and hence are higher by a factor of as compared to the errors for individual directions reported in Planck Collaboration XXVII (2014). 
The significant power detected in the reconstructed Doppler field (), from the CMB multipole bin ℓ = 2 to 64, shown in Fig. 36, corresponds to an effective  β  = 2.6 × 10^{2}. Note that this value of  β  is an order of magnitude larger than that of our local motion. Further, we find that the dipole direction is consistent with Fig. 3 of Planck Collaboration XXVII (2014), which depicts the reconstructed direction of the Doppler boost as a function of the CMB multipole.
Fig. 38 BipoSH power spectra and obtained from the four component separation maps (CR, NILC, SEVEM, and SMICA) are consistent with each other. Note that no significant detections are found. This independently establishes the fact that the quadrupolar BipoSH detections made by WMAP were due to WMAPspecific instrument systematics. 
Figure 37 presents a comparative study of the sensitivity of the two estimators, Doppler and modulation, and depicts the reconstruction noise for the two estimators as a function of the maximum CMB multipole used. The horizontal red dashdotted line at  β  = 1.23 × 10^{3} denotes the expected value of the Doppler signal induced by our local motion. The blue dashdotted line denotes the Doppler boost signal as seen by the modulation estimator and is expected to be enhanced by an effective factor of b_{ν} for the component separated maps.
For the case of the modulation estimator, it is seen that the Doppler boost is not expected to significantly contaminate the modulation signal up to ℓ_{max} ≃ 700, establishing the fact that the boost effect corresponding to our local motion is not strong enough to affect the modulation signal seen at low CMB multipoles. However, it is expected to significantly add to the modulation signal at higher CMB multipoles.
Figure 37 also clearly shows that the Doppler estimator is expected to recover the Doppler boost signal at a high significance. The figure also establishes that at ℓ ≲ 500, the local motion contribution is not detectable in the Doppler boost search, and hence the signal measured by the Doppler estimator in the range ℓ = 2–64, depicted in Fig. 36, is expected to be linked to the modulation signal seen at low CMB multipoles.
It is, of course, also possible to extract the BipoSH coefficients , up to the maximum multipole ℓ_{max} allowed by the full resolution Planck maps at modest computational expense. This allows us to address a specific indication of statistical isotropy violation previously reported in the literature. Bennett et al. (2011) found nonzero BipoSH power spectra, and , at very high statistical significance in the WMAP maps as determined in ecliptic coordinates, corresponding to a quadrupolar power asymmetry in the CMB sky. The BipoSH spectra peaked at around ℓ = 250, and the differences in the BipoSH signal determined from two different frequency bands indicated a noncosmological origin. Furthermore, the azimuthal symmetry of this BipoSH signal in ecliptic coordinates suggested that it had its origin in some unaccountedfor systematic effect. The findings of Hanson et al. (2010) and Joshi et al. (2012) strongly suggest that the signal arises due to an incomplete treatment of beam asymmetries in the data. Bennett et al. (2013) have subsequently noted that analysis of the WMAP9 beamdeconvolved maps no longer shows the signal.
We have computed the and values in ecliptic coordinates for the full resolution Planck CMB maps, as shown in Fig. 38. The analysis yields no evidence for BipoSH coefficients that deviate significantly from zero. This provides conclusive observational evidence from independent CMB measurements that the WMAP result could only have arisen due to instrumental artefacts in that data set.
In summary, a generalized search for power in the modulation field reveals a statistically significant dipolar modulation of the CMB sky on large angular scales (ℓ ≲ 64), with an amplitude of about 0.07. Conversely, no quadrupolar modulation is observed, reconfirming the systematic origin of the corresponding anomaly seen in the WMAP7 data.
5.8. Parity asymmetry
5.8.1. Pointparity asymmetry
The CMB sky map may be considered as the sum of even and odd parity functions. Previously, an odd pointparity preference (hereafter parity asymmetry) was observed in the WMAP 3, 5, and 7year data releases (Land & Magueijo 2005b; Kim & Naselsky 2010a,b; Gruppuso et al. 2011; Aluri & Jain 2012; Naselsky et al. 2012). In this section we investigate the parity asymmetry for the Planck temperature anisotropy power spectra derived with a quadratic maximum likelihood (QML) estimator applied to the CommanderRuler, NILC, SEVEM, and SMICA maps at N_{side} = 32, and with a pseudoC_{ℓ} estimator at N_{side} = 64.
From the CMB anisotropy field defined on the sky, T(n), we may construct symmetric and antisymmetric functions using the coordinate inversion n → − n: (40)Therefore, T^{+}(n) and T^{−}(n) have even and odd parity, respectively. When combined with the parity property of spherical harmonics, Y_{ℓm}(n) = ( − 1)^{ℓ}Y_{ℓm}( − n), we obtain (41)where n is an integer, and Γ^{+}(ℓ) = cos^{2}(ℓπ/ 2), and Γ^{−}(ℓ) = sin^{2}(ℓπ/ 2).
A significant power asymmetry between even and odd multipoles may thus be interpreted as a preference for a particular parity of the anisotropy pattern, connected to the parity asymmetry of the metric perturbations at scales above 1–4 Gpc (Kim & Naselsky 2010a). For investigation of the parity asymmetry we may consider the following quantities (Kim & Naselsky 2010a): (42)Here P^{+} and P^{−} are the sum of n(n + 1)/2πC_{n} for even and odd multipoles respectively; the ratio P^{+}/P^{−} is associated with the degree of parity asymmetry, where a value of P^{+}/P^{−}< 1 indicates oddparity preference, and P^{+}/P^{−}> 1 indicates evenparity preference.
Fig. 39 Upper: the parity estimator g(ℓ) versus ℓ for CR (black diamonds), NILC (red diamonds), SEVEM (blue diamonds), and SMICA (green diamonds). Lower: pvalue for CR (black solid line), NILC (red line), SEVEM (blue line), and SMICA (green line). 
Following Kim & Naselsky (2010a), we will discuss the range of multipoles 2 ≤ ℓ ≤ 30, which belongs to the SachsWolfe plateu of the TT power spectrum, where ℓ(ℓ + 1)C_{ℓ} ≈ constant. In order to make a rigorous assessment of the statistical significance of parity asymmetry at low ℓ, we have compared g(ℓ) for the Planck power spectra with 1000 simulated CMB maps based on the fiducial Planck cosmological model. We compute power spectra using a QMLestimator (Gruppuso et al. 2009) as applied to data at N_{side} = 32 with the U73 mask applied. This yields practically identical power spectrum results for the same ℓrange determined with a pseudoC_{ℓ} estimator applied to maps at N_{side} = 64.
In Fig 39 we show the g(ℓ)parameter for the Planck power spectra and the corresponding pvalue(ℓ). The pvalue denotes the fraction of simulations in which the obtained value of P^{+}/P^{−} is as low as that observed in the data. Note that the results from the different Planck CMB maps yield consistent shapes for the g(ℓ) parameter and pvalue(ℓ). The parity asymmetry at ℓ = 22 is most anomalous, with a corresponding pvalue in the range 0.002–0.004. Finally, the statistical significance of the parity asymmetry (i.e., low pvalue) increases when we increase ℓ_{max} up to 22–25. Therefore, the odd parity preference cannot simply be attributed to the low quadrupole power. It is plausible the low quadrupole power is not an isolated anomaly, but that it shares an origin with the odd parity preference (see for details Kim & Naselsky 2010a,b; Naselsky et al. 2012).
5.8.2. Mirror parity
In this section we investigate the properties of the Planck temperature lowresolution maps under reflection with respect to a plane. This search for hidden mirror symmetries and antisymmetries complements the tests for parity asymmetry, presented in Sect. 5.8.1. Starobinsky (1993) showed how a hidden mirror symmetry might be connected to the noncompact T^{1} topology, or to a compact T^{3} topology in which one topological scale is much shorter than the others. The CMB pattern would then exhibit a mirror symmetry with respect to the plane defined by the two large dimensions. Mirror symmetry has been searched for in the COBEDMR data in de OliveiraCosta et al. (1996), resulting in a lower limit for the scale of the compact dimension as 4 Gpc (see also Gurzadyan et al. 2007; BenDavid et al. 2012 for other more recent analysis). Finelli et al. (2012) analysed hidden mirror symmetry and antisymmetry properties of the WMAP 7year ILC temperature map, finding a preferred direction that could be considered anomalous at the 93% confidence level with antisymmetry properties. This direction lies close to the one defining the hemispherical asymmetry.
Following Finelli et al. (2012), we consider the following estimators (43)where the sum is meant over the HEALPix pixels, N_{pix}, δT/T(n_{j}) is the CMB temperature anisotropy measured at the pixel pointed by the unit vector n_{j}, and n_{k} is the opposite direction of n_{j} with respect to the plane defined by n_{i}, i.e., (44)The estimator S^{+} (S^{−}) can be written explicitly in the harmonic domain and its averaged value over isotropic random realization is (), where B_{ℓ} is the Gaussian beam applied to the lowresolution maps.
We compute the quantities S^{±} for each of the 3072 directions defined by HEALPix resolution N_{side} = 16 map, by allowing the j and k indices to run over the pixels of the low resolution foreground cleaned fullsky maps. We perform the same analysis on 1000 simulated skies and store the minimum and maximum value of S^{±} for each of these to correct our statistics for the elsewhere effect.
The minimum value for the S^{+} estimator is reached for the plane defined by Galactic coordinates (l,b) = (264°, − 17°) for all four foreground cleaned maps, with a significance of 0.5% (CommanderRuler), 1.4% (NILC), 3.1% (SEVEM), and 8.9% (SMICA). The top panel of Fig. 40 shows the minimum value of S^{+} for each of the four methods and compared to the MC simulations computed for CommanderRuler, which is considered to be representative.
The minimum value for the S^{−} estimator is found for a direction which depends slightly on the foreground cleaned CMB map considered and is close to the one associated with the cosmological dipole. This minimum value is not statistically significant for any of the foreground cleaned CMB maps (see bottom panel of Fig. 40).
Fig. 40 Top panel: S^{+} statistic. The vertical lines show the minimum value for the estimator as computed on low resolution CR, NILC, SEVEM, and SMICA maps. The grey histogram shows the same quantity computed from 1000 simulated maps processed by CR. Bottom panel: as above for S^{−}. 
The anomalous antisymmetry direction found in the Planck CMB data is close to that found for the dipolar modulation in Sect. 5.6, suggesting possible connections between them. The direction which minimizes S^{+} (S^{−}) for Planck is almost exactly the same as the one found for the WMAP 7year ILC map in Finelli et al. (2012).
5.9. The Cold Spot
The Cold Spot was identified in the WMAP first year data (Vielva et al. 2004) through the estimation of the kurtosis of the SMHW (e.g., MartínezGonzález et al. 2002) coefficients, and confirmed (Cruz et al. 2005) by analysing the area of the SMHW coefficients above/below a given threshold. Since its detection, the Cold Spot has been extensively studied and verified with a large battery of statistical probes (e.g., Mukherjee & Wang 2004; Cayón et al. 2005; McEwen et al. 2005; Cruz et al. 2007a; Räth et al. 2007b; Vielva et al. 2007; Pietrobon et al. 2008; Gurzadyan et al. 2009; Rossmanith et al. 2009b). A complete review of the Cold Spot can be found in Vielva (2010), including a discussion on possible explanations of its nature.
The analysis of the kurtosis of the SMHW coefficients has already been addressed in Sect. 4.5. We have checked that the kurtosis of the coefficients corresponding to the four Planck cleaned frequency maps is larger than the expected value obtained from simulations, with a modified upper tail probability of around 0.01. This is compatible with the value obtained from WMAP.
Nevertheless, the Cold Spot is more robustly described in terms of a morphological quantity: the area of the SMHW coefficients above/below a given threshold. At a given scale R and threshold ν, the cold () and hot () areas of the SMHW coefficients are defined as: (45)Here # represents the number operator, i.e., it indicates for how many pixels p, the specific condition defined between the braces is satisfied.
Table 30 summarizes the results for the hot and cold areas determined for the four CMB maps analysed with the U73 mask (and its associated exclusions masks). The cold area is anomalous at scales between R = 200 and R = 300′, similar to the sizes already highlighted with the kurtosis analysis. We see that the higher the threshold, the smaller the upper tail probability associated with the Planck CMB map. In particular, the cold area has a upper tail probability of 0.003 at ν> 4σ_{R} and for R = 300′.
Upper tail probability (UTP, in %) associated to the cold (top) and hot (bottom) areas.
Notice that the most significant deviation comes from the cold area, although the hot area is marginally compatible. However, the cold area represents the most robust detection of an anomaly, since it is robust to the mask employed (see Tables 31 and 32).
Upper tail probability (in %) associated to the cold (top) and hot (bottom) areas.
Upper tail probabilities (in %) associated with the cold (top) and hot (bottom) areas.
The information provided in the previous tables is also represented (for the R = 300′ scale) in Fig. 41. In these panels we show the anomalous cold (in blue) and hot (in red) areas for thresholds ν> 3.0σ_{R} as determined from the SMICA map. The most extreme value (in terms of σ_{R}) of each area is indicated. The coldest area corresponds to the Cold Spot with the minimum value of the wavelet coeficient at the position (209°,−57°) in Galactic coordinates, whereas the hotest area has already been identified in the WMAP data (e.g., Vielva et al. 2007) as an anomalous hot spot. This does not depends signifincantly on the mask considered in the analysis. From these results it is clear that the Cold Spot anomaly is present in both the WMAP and Planck data.
5.10. Interpretation of anomalies
The results presented here in Sect. 5 demonstrate that many features previously observed in the WMAP data are present also in the Planck sky. This agreement between two independent experiments effectively rules out the possibility that their origin lies in systematic artefacts present in either data set. In particular, there is evidence for a violation of statistical isotropy at least on large angular scales in the context of the Planck fiducial sky model. Moreover, a power asymmetry extends to scales corresponding to ℓ ≃ 600, whilst fits at low ℓ to a model containing a dipole modulation field yield results in excess of 3σ significance. In addition, there is evidence from such fits that the lowℓ spectrum of the Planck data departs from the fiducial spectrum in both amplitude and slope. These results could have profound implications for cosmology. It is therefore pertinent to consider whether a model can be proposed to provide a common origin for the anomalies.
Fig. 41 SMHW coefficients at R = 300 arcmin, above (and below) a + 3.0σ (−3.0σ) threshold. Results for the three masks considered in the analysis are shown: U73 mask (top); CG70 (middle); and CG60 (bottom). 
The microwave sky is manifestly nonGaussian and anisotropic, with known contributions from Galactic astrophysical foregrounds, lensing of CMB anistropies by the intervening matter distribution, and the ISW. However, the excellent performance of the component separation algorithms used here in rejecting diffuse foregrounds argues strongly against known Galactic emission as the source of the anomalies.
Schwarz et al. (2004), Copi et al. (2007), Maris et al. (2011) and Hansen et al. (2012) suggested that diffuse Solar System emission could contribute to the observed structure on large angular scales, although it is not expected that the classical Zodiacal Light Emission or Kuiper Belt objects are responsible. Planck Collaboration XIV (2014) presents the current Planck contribution to the modelling of the Zodiacal cloud. Moreover, the association of the anomalies with the ecliptic reference frame seems, at best, to be confined to the very largest angular scales, since the preferred direction determined for the dipolar modulation field analysis is separated from the south ecliptic pole by more than 45°, with even larger separations found for the power asymmetry direction.
Another possibility is that the anomalies have their origin in the local Universe. According to Francis & Peacock (2009), the removal of the ISW signal originating within the volume at z< 0.3 from WMAP data reduces the significance of the apparent alignment between the CMB quadrupole and octopole and the Cold Spot. Efstathiou et al. (2010) have used the same correction to yield an increase in the structure of the twopoint correlation function for angular separations less than 60°, that had been noted as apparently anomalous since the first WMAP data release. A future possibility is that Planck itself will be able to reconstruct the ISW signal and test its impact on issues related to isotropy and nonGaussianity. Planck Collaboration XIX (2014) presents maps of the effect based on the current data release.
Of more interest to us is that the anomalies are genuinely cosmological in origin. In that context, obvious candidate models include those with simply or multiconnected topology. In a companion paper (Planck Collaboration XXVI 2014), a subset of such models are considered and the signatures of their specific correlation structures on the sky are searched for. However, no detections are found, but rather the scale of topology is limited to be of order the diameter of the lastscattering surface or greater. More interestingly, they reconsider Bianchi VII_{h} models that were previously demonstrated to show statistical correlation with the WMAP data (Jaffe et al. 2005, 2006; Bridges et al. 2007; McEwen et al. 2013), albeit with parameters inconsistent with standard cosmological parameters. In this new analysis, the Bianchi parameters are physically coupled to the cosmological ones, yielding no evidence for a Bianchi VII_{h} cosmology. However, as before, when treated simply as a template for structure contained in the CMB sky, a bestfit pattern is found to be in good agreement with the old results. Previous analyses (Jaffe et al. 2005; Cayón et al. 2006; McEwen et al. 2006) have shown that when the CMB sky is corrected for such a template, many of the largescale anomalies are no longer present at a statistically significant level. It is likely that such an effect will persist for Bianchicorrected Planck data, and we have made an explicit test as to whether the bestfit Bianchi template can also explain the presence of phase correlations. We therefore repeated the surrogate analysis from Sect. 5.4 for the appropriately corrected SMICA map. Figure 42 presents the result for the corresponding significance map. It is clear that the signature for hemispherical asymmetry is drastically reduced, thereby rendering the signal formally statistically insignificant. Thus, the bestfit Bianchi model can also account for the asymmetries induced by higher order phase correlations. It should also be noted that subtracting the bestfit Bianchi template from the data, outside the U73 mask, explains the anomalous skewness and kurtosis values but not the variance, for which the corresponding lower tail probabilities are 0.008, 0.166, and 0.306, respectively. Given the lack of consistency of the physical parameters of the model with the Planck cosmological model, the results obtained using Bianchisubtracted input maps might be considered moot, however, the morphology of the maps may provide insight into the type of underlying structures associated with the anomalies.
Although the Cold Spot is also rendered statistically insignificant by the Bianchi template, other possible explanations about its nature have been considered, including the late evolution of the largescale structure (e.g., Inoue & Silk 2006, 2007), the SunyaevZeldovich effect (e.g., Cruz et al. 2008), residual foregrounds (Cruz et al. 2006), or a cosmic texture (e.g., Cruz et al. 2007b). Further tests based on gravitational lensing (Das & Spergel 2009) or CMB polarization (Vielva et al. 2011) have been suggested to confirm or reject some of the proposed explanations.
The presence of primordial magnetic fields (PMFs) due to either pre or postrecombination mechanisms could also provide a physical basis for some of the anomalies discussed in this paper. Specifically, PMFs with coherence scales comparable to the present day horizon could result in Alfvén waves in the early Universe that generate specific signatures on the sky via the Doppler and integrated SachsWolfe effects. In particular, a preferred angular direction in the CMB anisotropy can be induced (Durrer et al. 1998; Kim & Naselsky 2009) that leads to structure in the spherical harmonic mode correlation matrix (Kahniashvili et al. 2008). Appendix A presents a search for the predicted correlations between harmonic modes separated by Δℓ = 0, ± 2, and Δm = 0, ± 1, ± 2, allowing constraints to be placed on the Alfvén wave amplitude. Further constraints on PMFs based on the power spectrum and bispectrum have been provided in Planck Collaboration XVI (2014) and Planck Collaboration XXIV (2014), respectively.
To conclude, when analysing a data set as complex and rich as that provided by Planck, some statistical outliers will be expected. However, it should be clear that the evidence for some of the largeangular scale anomalies is significant indeed, yet few physically compelling models have been proposed to account for them, and none so far that provide a common origin. The dipole modulation model that was analysed here was phenomenologically motivated, and is detected in the data at relatively high significance. Whether it can resolve the anomalous nature of other observed features remains to be evaluated.
6. Implications for C_{ℓ} and cosmological parameter estimation
The approach to C_{ℓ} estimation, the construction of the Planck likelihood and subsequent inference of cosmological parameters are described in the accompanying papers Planck Collaboration XV (2014) and Planck Collaboration XVI (2014). For these studies, specific assumptions are made about the isotropy and Gaussianity of the primordial fluctuations observed in the CMB. The latter in particular seems to be wellsupported by the comprehensive set of tests applied to the Planck data in Sect. 4. The most significant discrepancies are seen in association with the Cold Spot (Sect. 5.9), which constitutes a localized region of sky of about 10° in size. Its impact on cosmological parameters is then likely to be relatively insignificant, and masking of the region could easily test this assertion.
It is wellknown that the quadrupole and octopole have low amplitudes relative to the bestfit cosmological powerspectrum. The contribution of those multipoles to cosmological parameter estimation is very small due to the associated cosmic variance on these scales, and thus the direct impact of their alignment (as discussed in Sect. 5.1) is also likely to be small. Remarkably, however, Planck Collaboration XV (2014) presents evidence that the lowℓ multipole range from 2 to 30 is coherently low, and is not well accounted for by the standard ΛCDM model. Moreover, this conclusion is a consequence of the fact that the cosmological parameters are strongly influenced by the ℓ = 1000–1500 range, previously inaccessible to WMAP. Consistent findings have been presented here in the form of the low variance of the data in Sects. 4.1 and 5.2, although this is largely driven by the quadupole and octopole alone.
The question therefore remains as to whether there is a deeper connection with the cosmological anomalies seen in both the WMAP and Planck data sets particularly on large angular scales. Indeed, the hemispherical asymmetry and dipolar power modulation discussed in Sect. 5 could have a more important impact in that they directly address whether a broader class of cosmological models should be considered. Indeed, the lowℓ signature seen in the data has previously been associated with missing power in a universe with simply or multiplyconnected topology. However, there are specific morphological signatures of such topologies that have not been detected in the Planck data (Planck Collaboration XXVI 2014).
Fig. 43 The angular power spectrum on large angular scales computed from the SMICA map at N_{side} = 16 on the opposing hemispheres defined by the maximal asymmetry direction (l,b) = (220°,4°) determined in Sect. 5.5. 
However, the phenomenologically motivated dipole modulation model due to Gordon et al. (2005) yields a significant fit to the data, as seen in Sects. 5.6 and 5.7. The former also shows some evidence for a departure from the Planck fiducial power spectrum in both amplitude and slope. Both of these analyses are in good agreement in terms of the direction of the dipolar modulation field with the model independent power asymmetry analysis of Sect. 5.5.
A qualitative exploration as to how these features are reflected in the lowℓ power spectrum is provided in Fig. 43. Specifically, the plot presents the angular power spectra computed using a quadratic maximumlikelihood (QML) estimator (Paci et al. 2010, 2013) from the N_{side} = 16SMICA map after application of the U73 mask used in this paper. The Planck fiducial power spectrum is also shown for comparison. Clearly, there is a deficit of power as expected when no further partitioning of the sky is applied. However, further interesting properties of the data are revealed when spectra are computed for the two opposing hemispheres defined by the preferred direction from the ℓ_{max} = 1500 analysis in Sect. 5.5. In the positive direction, there is improved agreement between the derived spectrum and the Planck fiducial sky, but with an interesting oscillation between odd and even modes. For the negative direction, an overall suppression of power is again seen. It would be interesting to test the connection between these spectral features and the phase correlations detected in Sect. 5.4 or the evidence for parity violation presented in Sect. 5.8. The observations may, in part, reflect the presence of visually striking features noted by Bennett et al. (2011) – the four elongated cold fingers stretching from near the Galactic equator to the south Galactic pole and a prominent cold spot near the centre of the map. Very similar results are found when the preferred direction is specified either by the ℓ_{max} = 600 power asymmetry analysis, or the dipolar modulation direction in Sect. 5.6.
However, Fig. 30 and the corresponding analysis suggest that the asymmetry in power between hemispheres extends to much smaller angular scales, whilst Planck Collaboration XXVII (2014) demonstrates the presence of a modulation due to Doppler boosting up to ℓ ≃ 2000. Whether such properties of the data would have implications for parameter estimation may need further exploration.
7. Conclusions
In this paper, we have tested the statistical isotropy and Gaussianity of the CMB using data from the Planck satellite. We have demonstrated that little evidence is seen for nonGaussianity, although some deviations from isotropy are found.
Most of the tests performed in Sect. 4 showed an overall consistency with the null hypothesis, as represented by a set of realistic simulations based on a Planck fiducial sky model and including the secondary ISWlensing effect (as detected for the first time with the Planck data, see Planck Collaboration XIX 2014). However, two important exceptions were seen. The variance of the CMB signal was found to be significantly lower than expected, with the anomalously low signal seemingly localized in the northern ecliptic hemisphere (Sect. 4.1). This result was also confirmed with the low variance of the wavelet coefficients that was seen on scales above a few degrees (see Sect. 4.5). Moreover, a significant deviation from Gaussianity was found in the form of a positive kurtosis of the wavelet coefficients.
These results correspond to statistical features on large angular scales where numerous anomalies were previously observed in the WMAP data. In Sect. 5, we revisited these in the light of the Planck data and found results in excellent agreement with those for WMAP. In particular, the most significant anomalies, namely the quadrupoleoctopole alignment (Sect. 5.1), the low variance (Sect. 5.2), hemispherical asymmetry (Sect. 5.3), phase correlations (Sect. 5.4), power asymmetry (Sect. 5.5), dipolar modulation (Sect. 5.6), generalized power modulation (Sect. 5.7), parity asymmetry (Sect. 5.8), and the Cold Spot (Sect. 5.9) have been confirmed with the Planck data. Attempts to explain the observed features in terms of systematic artefacts, local astrophysical sources of emission, or structure in the local Universe have not been successful. It is clear that these anomalies represent real features of the CMB sky.
However, it is difficult to make a detailed interpretation of the anomalies in the absence of theoretical models, in particular with regard to the role of a posterior choices. Nevertheless, Planck does offer new possibilities to check the a posteriori claims in this context as a consequence of its superior multipole content that cannot easily been probed by WMAP.
Phenomenological models have been suggested to account for the observations. The dipolar power modulation approach due to Gordon et al. (2005) was explicitly tested in Sect. 5.6 and found to represent a good fit to the large scale asymmetry, corresponding to a detection at about 3σ significance. This result was confirmed by the more generalized modulation study in Sect. 5.7, which also ruled out the presence of modulation fields of higher order. Alternatively, a Bianchi template fit to the data performed in Planck Collaboration XXVI (2014) can provide a good fit to the hemispherical asymmetry, the Cold Spot and the phase correlations, but corresponds to values of the cosmological parameters incompatible with those derived in Planck Collaboration XVI (2014). Clearly, these do not provide complete and satisfactory explanations for the observations, and more physically motivated models should be sought.
This may also be indicated by the cosmological parameter studies presented in Planck Collaboration XV (2014) and Planck Collaboration XVI (2014). Here, it was demonstrated that while the power spectrum determined from the Planck temperature data is extremely consistent with a basic 6parameter ΛCDM model, the lowℓ multipoles (ℓ ≤ 30) deviate from the bestfit model although at a significance that does not appear to exceed 2.7σ. However, this is precisely the regime where many of the anomalies presented in this paper seem to manifest themselves, and where qualitatively interesting differences are observed in the powerspectra for two hemispheres defined by the preferred direction for the dipolar power modulation.
Finally, it is expected that the polarization data that will become available with the 2014 data release should provide valuable information on the nature of the CMB anomalies. Then, the presence, or even absence, of a specific signature in the data should help to elucidate the physical mechanism that is causing the anomaly (see Vielva et al. 2011; Frommert & Enßlin 2010; Dvorkin et al. 2008 for examples related to the Cold Spot, mode alignment, and dipolar modulation, respectively). In particular, a deviation of isotropy present at recombination should be reflected in both the temperature and polarization data, with a correlated signal. It may be that the statistical anomalies currently described in this paper are a hint of more profound physical phenomena that are yet to be revealed.
Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states (in particular the lead countries France and Italy), with contributions from NASA, (USA) and telescope reflectors provided by a collaboration between ESA and a scientific conosrtium led and funded by Denmark.
Although Bennett et al. (2013) used a variant on the estimator adopted in this paper, we also find a 3° misalignment of the quadrupole and octopole determined from the ILC map.
Note that the WMAP direction and pvalue is slightly different from the numbers found in Axelsson et al. (2013), due to small differences in the analysis. Here, we use 16multipole bins and the direction is determined from the mean dipole direction of all bins in the specified multipole range rather than from one single bin spanning the full multipole range.
Acknowledgments
The development of Planck has been supported by: ESA; CNES and CNRS/INSUIN2P3INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN, JA and RES (Spain); Tekes, AoF and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and PRACE (EU). A description of the Planck Collaboration and a list of its members, including the technical or scientific activities in which they have been involved, can be found at http://www.sciops.esa.int/index.php?project=planck&page=Planck_Collaboration. We acknowledge the use of resources from the Norewegian national super computing facilities NOTUR. The modal and KSW bispectrum estimator analysis was performed on the COSMOS supercomputer, part of the STFC DiRAC HPC Facility. We further acknowledge the computer resources and technical assistance provided by the Spanish Supercomputing Network nodes at Universidad de Cantabria and Universidad Politécnica de Madrid as well as by the Advanced Computing and eScience team at IFCA. Some of the results in this paper have been derived using the HEALPix package.
References
 Aluri, P. K., & Jain, P. 2012, MNRAS, 419, 3378 [NASA ADS] [CrossRef] [Google Scholar]
 Axelsson, M., Fantaye, Y., Hansen, F. K., et al. 2013, ApJ, 773, L3 [NASA ADS] [CrossRef] [Google Scholar]
 BenDavid, A., Kovetz, E. D., & Itzhaki, N. 2012, ApJ, 748, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, C. L., Hill, R. S., Hinshaw, G., et al. 2011, ApJS, 192, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20 [Google Scholar]
 Bielewicz, P., & Riazuelo, A. 2009, MNRAS, 396, 609 [NASA ADS] [CrossRef] [Google Scholar]
 Bielewicz, P., Górski, K. M., & Banday, A. J. 2004, MNRAS, 355, 1283 [NASA ADS] [CrossRef] [Google Scholar]
 Bielewicz, P., Eriksen, H. K., Banday, A. J., Górski, K. M., & Lilje, P. B. 2005, ApJ, 635, 750 [NASA ADS] [CrossRef] [Google Scholar]
 Bielewicz, P., Wandelt, B. D., & Banday, A. J. 2013, MNRAS, 429, 1376 [NASA ADS] [CrossRef] [Google Scholar]
 Bridges, M., McEwen, J. D., Lasenby, A. N., & Hobson, M. P. 2007, MNRAS, 377, 1473 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bunn, E. F., & Scott, D. 2000, MNRAS, 313, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Cayón, L., Jin, J., & Treaster, A. 2005, MNRAS, 362, 826 [NASA ADS] [CrossRef] [Google Scholar]
 Cayón, L., Banday, A. J., Jaffe, T., et al. 2006, MNRAS, 369, 598 [NASA ADS] [CrossRef] [Google Scholar]
 Challinor, A., & van Leeuwen, F. 2002, Phys. Rev. D, 65, 103001 [NASA ADS] [CrossRef] [Google Scholar]
 Copi, C. J., Huterer, D., & Starkman, G. D. 2004, Phys. Rev. D, 70, 043515 [NASA ADS] [CrossRef] [Google Scholar]
 Copi, C. J., Huterer, D., Schwarz, D. J., & Starkman, G. D. 2007, Phys. Rev. D, 75, 023507 [Google Scholar]
 Cruz, M., MartínezGonzález, E., Vielva, P., & Cayón, L. 2005, MNRAS, 356, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Cruz, M., Tucci, M., MartínezGonzález, E., & Vielva, P. 2006, MNRAS, 369, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Cruz, M., Cayón, L., MartínezGonzález, E., Vielva, P., & Jin, J. 2007a, ApJ, 655, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Cruz, M., Turok, N., Vielva, P., MartínezGonzález, E., & Hobson, M. 2007b, Science, 318, 1612 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Cruz, M., MartínezGonzález, E., Vielva, P., et al. 2008, MNRAS, 390, 913 [NASA ADS] [CrossRef] [Google Scholar]
 Cruz, M., Vielva, P., MartínezGonzález, E., & Barreiro, R. B. 2011a, MNRAS, 412, 2383 [NASA ADS] [CrossRef] [Google Scholar]
 Curto, A., Aumont, J., MacíasPérez, J. F., et al. 2007, A&A, 474, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Curto, A., MacíasPérez, J. F., MartínezGonzález, E., et al. 2008, A&A, 486, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Curto, A., MartínezGonzález, E., & Barreiro, R. B. 2009a, ApJ, 706, 399 [NASA ADS] [CrossRef] [Google Scholar]
 Curto, A., MartínezGonzález, E., Mukherjee, P., et al. 2009b, MNRAS, 393, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Curto, A., MartínezGonzález, E., & Barreiro, R. B. 2010, in Highlights of Spanish Astrophysics V, eds. J. M. Diego, L. J. Goicoechea, J. I. GonzálezSerrano, & J. Gorgas (Berlin, Heidelberg: Springer Verlag), 277 [Google Scholar]
 Curto, A., MartínezGonzález, E., & Barreiro, R. B. 2011a, MNRAS, 412, 1038 [NASA ADS] [Google Scholar]
 Curto, A., MartínezGonzález, E., Barreiro, R. B., & Hobson, M. P. 2011b, MNRAS, 417, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Curto, A., MartinezGonzalez, E., & Barreiro, R. B. 2012, MNRAS, 426, 1361 [NASA ADS] [CrossRef] [Google Scholar]
 Das, S., & Spergel, D. N. 2009, Phys. Rev. D, 79, 043007 [NASA ADS] [CrossRef] [Google Scholar]
 de OliveiraCosta, A., Smoot, G. F., & Starobinsky, A. A. 1996, ApJ, 468, 457 [NASA ADS] [CrossRef] [Google Scholar]
 de OliveiraCosta, A., Tegmark, M., Zaldarriaga, M., & Hamilton, A. 2004, Phys. Rev. D, 69, 063516 [NASA ADS] [CrossRef] [Google Scholar]
 De Troia, G., Ade, P. A. R., Bock, J. J., et al. 2007, ApJ, 670, L73 [NASA ADS] [CrossRef] [Google Scholar]
 Ducout, A., Bouchet, F. R., Colombi, S., Pogosyan, D., & Prunet, S. 2012, MNRAS, 423 [Google Scholar]
 Durrer, R., Kahniashvili, T., & Yates, A. 1998, Phys. Rev. D, 58, 123004 [NASA ADS] [CrossRef] [Google Scholar]
 Dvorkin, C., Peiris, H. V., & Hu, W. 2008, Phys. Rev. D, 77, 063008 [NASA ADS] [CrossRef] [Google Scholar]
 Efstathiou, G. 2004, MNRAS, 348, 885 [NASA ADS] [CrossRef] [Google Scholar]
 Efstathiou, G., Ma, Y.Z., & Hanson, D. 2010, MNRAS, 407, 2530 [NASA ADS] [CrossRef] [Google Scholar]
 Elsner, F., & Wandelt, B. D. 2013, A&A, 549, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eriksen, H. K., Hansen, F. K., Banday, A. J., Górski, K. M., & Lilje, P. B. 2004a, ApJ, 605, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Lilje, P. B., Banday, A. J., & Górski, K. M. 2004b, ApJS, 151, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Novikov, D. I., Lilje, P. B., Banday, A. J., & Górski, K. M. 2004c, ApJ, 612, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Banday, A. J., Górski, K. M., & Lilje, P. B. 2005, ApJ, 622, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Banday, A. J., Górski, K. M., Hansen, F. K., & Lilje, P. B. 2007a, ApJ, 660, L81 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Huey, G., Saha, R., et al. 2007b, ApJ, 656, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Fergusson, J., Liguori, M., & Shellard, E. 2010 [Google Scholar]
 Ferreira, P. G., & Magueijo, J. 1997, Phys. Rev. D, 56, 4578 [NASA ADS] [CrossRef] [Google Scholar]
 Finelli, F., Gruppuso, A., Paci, F., & Starobinsky, A. 2012, JCAP, 1207, 049 [NASA ADS] [CrossRef] [Google Scholar]
 Fixsen, D. J. 2009, ApJ, 707, 916 [Google Scholar]
 Flender, S., & Hotchkiss, S. 2013, J. Cosmol. Astropart. Phys., 9, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Francis, C. L., & Peacock, J. A. 2009, MNRAS, 406, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Frommert, M., & Enßlin, T. A. 2010, MNRAS, 403, 1739 [NASA ADS] [CrossRef] [Google Scholar]
 Gordon, C., Hu, W., Huterer, D., & Crawford, T. 2005, Phys. Rev. D, 72, 103002 [NASA ADS] [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Gruppuso, A., de Rosa, A., Cabella, P., et al. 2009, MNRAS, 400, 463 [NASA ADS] [CrossRef] [Google Scholar]
 Gruppuso, A., Finelli, F., Natoli, P., et al. 2011, MNRAS, 411, 1445 [NASA ADS] [CrossRef] [Google Scholar]
 Gurzadyan, V., Starobinsky, A., Kashin, A., Khachatryan, H., & Yegorian, G. 2007, Mod. Phys. Lett., A22, 2955 [Google Scholar]
 Gurzadyan, V. G., Allahverdyan, A. E., Ghahramanyan, T., et al. 2009, A&A, 497, 343 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hajian, A., & Souradeep, T. 2003, ApJ, 597, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Hajian, A., & Souradeep, T. 2006, Phys. Rev. D, 74, 123521 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, F. K., Banday, A. J., & Górski, K. M. 2004, MNRAS, 354, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, F. K., Banday, A. J., Górski, K. M., Eriksen, H. K., & Lilje, P. B. 2009, ApJ, 704, 1448 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, M., Kim, J., Frejsel, A. M., et al. 2012, J. Cosmol. Astropart. Phys., 10, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Hanson, D., & Lewis, A. 2009, Phys. Rev. D, 80, 063004 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hanson, D., Lewis, A., & Challinor, A. 2010, Phys. Rev. D, 81, 103003 [NASA ADS] [CrossRef] [Google Scholar]
 Hikage, C., Matsubara, T., Coles, P., et al. 2008, MNRAS, 389, 1439 [NASA ADS] [CrossRef] [Google Scholar]
 Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Hoftuft, J., Eriksen, H. K., Banday, A. J., et al. 2009, ApJ, 699, 985 [NASA ADS] [CrossRef] [Google Scholar]
 Inoue, K. T., & Silk, J. 2006, ApJ, 648, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Inoue, K. T., & Silk, J. 2007, ApJ, 664, 650 [NASA ADS] [CrossRef] [Google Scholar]
 Jaffe, T. R., Banday, A. J., Eriksen, H. K., Górski, K. M., & Hansen, F. K. 2005, ApJ, 629, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Jaffe, T. R., Hervik, S., Banday, A. J., & Górski, K. M. 2006, ApJ, 644, 701 [NASA ADS] [CrossRef] [Google Scholar]
 Joshi, N., Das, S., Rotti, A., Mitra, S., & Souradeep, T. 2012 [arXiv:1210.7318] [Google Scholar]
 Kahniashvili, T., Lavrelashvili, G., & Ratra, B. 2008, Phys. Rev. D, 78, 063012 [NASA ADS] [CrossRef] [Google Scholar]
 Kamionkowski, M., & Knox, L. 2003, Phys. Rev. D, 67, 063001 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J., & Naselsky, P. 2009, J. Cosmol. Astropart. Phys., 7, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J., & Naselsky, P. 2010a, ApJ, 714, L265 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J., & Naselsky, P. 2010b, Phys. Rev. D, 82, 063002 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., Kogut, A., Nolta, M. R., et al. 2003, ApJS, 148, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., Dunkley, J., Nolta, M. R., et al. 2009, ApJS, 180, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Kronberg, P. P. 2009, in IAU Symp., 259, 499 [Google Scholar]
 Land, K., & Magueijo, J. 2005a, Phys. Rev. Lett., 95, 071301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Land, K., & Magueijo, J. 2005b, Phys. Rev. D, 72, 101302 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511 [NASA ADS] [CrossRef] [Google Scholar]
 Maris, M., Burigana, C., Gruppuso, A., Finelli, F., & Diego, J. M. 2011, MNRAS, 415, 2546 [NASA ADS] [CrossRef] [Google Scholar]
 MartínezGonzález, E., Gallegos, J. E., Argüeso, F., Cayón, L., & Sanz, J. L. 2002, MNRAS, 336, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Matsubara, T. 2010, Phys. Rev. D, 81, 083505 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, J. D., Hobson, M. P., Lasenby, A. N., & Mortlock, D. J. 2005, MNRAS, 359, 1583 [Google Scholar]
 McEwen, J. D., Hobson, M. P., Lasenby, A. N., & Mortlock, D. J. 2006, MNRAS, 369, 1858 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, J. D., Josset, T., Feeney, S. M., Peiris, H. V., & Lasenby, A. N. 2013, MNRAS, 436, 3680 [NASA ADS] [CrossRef] [Google Scholar]
 Mecke, K. R., Buchert, T., & Wagner, H. 1994, A&A, 288, 697 [NASA ADS] [Google Scholar]
 Mitra, S., Rocha, G., Górski, K. M., et al. 2011, ApJS, 193, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Modest, H. I., Räth, C., Banday, A. J., et al. 2013, MNRAS, 428, 551 [NASA ADS] [CrossRef] [Google Scholar]
 Monteserín, C., Barreiro, R. B., Vielva, P., et al. 2008, MNRAS, 387, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Moss, A., Scott, D., Zibin, J. P., & Battye, R. 2011, Phys. Rev. D, 84, 023014 [NASA ADS] [CrossRef] [Google Scholar]
 Mukherjee, P., & Wang, Y. 2004, ApJ, 613, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Naselsky, P., Zhao, W., Kim, J., & Chen, S. 2012, ApJ, 749, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Notari, A., Quartin, M., & Catena, R. 2013, [arXiv:1304.3506] [Google Scholar]
 Paci, F., Gruppuso, A., Finelli, F., et al. 2010, MNRAS, 407, 399 [NASA ADS] [CrossRef] [Google Scholar]
 Paci, F., Gruppuso, A., Finelli, F., et al. 2013, MNRAS, 434, 3071 [NASA ADS] [CrossRef] [Google Scholar]
 Park, C.G. 2004, MNRAS, 349, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Pietrobon, D., Amblard, A., Balbi, A., et al. 2008, Phys. Rev. D, 78, 103504 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration 2013, The Explanatory Supplement to the Planck 2013 results, http://www.sciops.esa.int/wikiSI/planckpla/index.php?title=Main_Page (ESA) [Google Scholar]
 Planck Collaboration I. 2014, A&A, 571, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2014, A&A, 571, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2014, A&A, 571, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2014, A&A, Vol. A4 [Google Scholar]
 Planck Collaboration V. 2014, A&A, 571, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2014, A&A, 571, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2014, A&A, 571, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2014, A&A, 571, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2014, A&A, 571, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2014, A&A, 571, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2014, A&A, 571, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2014, A&A, 571, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIV. 2014, A&A, 571, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XV. 2014, A&A, 571, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVII. 2014, A&A, 571, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2014, A&A, 571, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2014, A&A, 571, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XX. 2014, A&A, 571, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXI. 2014, A&A, 571, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2014, A&A, 571, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIII. 2014, A&A, 571, A23 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Planck Collaboration XXIV. 2014, A&A, 571, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXV. 2014, A&A, 571, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2014, A&A, 571, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2014, A&A, 571, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVIII. 2014, A&A, 571, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIX. 2014, A&A, 571, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXXI. 2014, A&A, 571, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Räth, C., Schuecker, P., & Banday, A. J. 2007a, MNRAS, 380, 466 [NASA ADS] [CrossRef] [Google Scholar]
 Räth, C., Schuecker, P., & Banday, A. J. 2007b, MNRAS, 380, 466 [NASA ADS] [CrossRef] [Google Scholar]
 Räth, C., Morfill, G. E., Rossmanith, G., Banday, A. J., & Górski, K. M. 2009, Phys. Rev. Lett., 102, 131301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Räth, C., Banday, A. J., Rossmanith, G., et al. 2011, MNRAS, 415, 2205 [NASA ADS] [CrossRef] [Google Scholar]
 Rossmanith, G., Räth, C., Banday, A. J., & Morfill, G. 2009a, MNRAS, 399, 1921 [NASA ADS] [CrossRef] [Google Scholar]
 Rossmanith, G., Räth, C., Banday, A. J., & Morfill, G. 2009b, MNRAS, 399, 1921 [NASA ADS] [CrossRef] [Google Scholar]
 Rossmanith, G., Modest, H., Räth, C., et al. 2012, Phys. Rev. D, 86, 083005 [NASA ADS] [CrossRef] [Google Scholar]
 Sarkar, D., Huterer, D., Copi, C. J., Starkman, G. D., & Schwarz, D. J. 2011, Astropart. Phys., 34, 591 [NASA ADS] [CrossRef] [Google Scholar]
 Schmalzing, J., & Buchert, T. 1997, ApJ, 482, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Schmalzing, J., & Gorski, K. M. 1998, MNRAS, 297, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Schwarz, D. J., Starkman, G. D., Huterer, D., & Copi, C. J. 2004, Phys. Rev. Lett., 93, 221301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Spergel, D. N., Bean, R., Doré, O., et al. 2007, ApJS, 170, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Starobinsky, A. A. 1993, JETP Lett., 57, 622 [NASA ADS] [Google Scholar]
 Tegmark, M., de OliveiraCosta, A., & Hamilton, A. 2003, Phys. Rev. D, 68, 123523 [NASA ADS] [CrossRef] [Google Scholar]
 Vanmarcke, E. 1983, in Random Fields (Cambridge, Massachusetts, USA: The MIT Press), 372 [Google Scholar]
 Vielva, P. 2010, Adv. Astron., 2010, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Vielva, P., MartínezGonzález, E., Barreiro, R. B., Sanz, J. L., & Cayón, L. 2004, ApJ, 609, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Vielva, P., Wiaux, Y., MartínezGonzález, E., & Vandergheynst, P. 2007, MNRAS, 381, 932 [NASA ADS] [CrossRef] [Google Scholar]
 Vielva, P., MartínezGonzález, E., Cruz, M., Barreiro, R. B., & Tucci, M. 2011, MNRAS, 410, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Wiaux, Y., Vielva, P., Barreiro, R. B., MartínezGonzález, E., & Vandergheynst, P. 2008, MNRAS, 385, 939 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Constraints on Alfvén waves
Observations of synchrotron emission and Faraday rotation provide increasing evidence that largescale astrophysical systems in the Universe are pervaded by magnetic fields. These huge systems include the Lyα forest and intercluster regions (see Kronberg 2009, for a review). Both pre and postrecombination mechanisms could result in a background of nanogauss fields that might be detectable in largescale structures or the CMB, although at present no imprints of these Primordial Magnetic Fields (PMFs) have been detected therein.
Here, we report our findings based on an analysis of the Planck data to search for the the predicted signature of statistical anisotropy due to PMFs. Specifically, PMFs with coherence scales comparable to the present day horizon may induce and sustain Alfvén waves in the early Universe that can leave observable imprints on the CMB via the Doppler and integrated SachsWolfe effects. In particular, this results in a preferred angular direction in the CMB anisotropy, therefore breaking statistical isotropy (Durrer et al. 1998; Kim & Naselsky 2009).
Durrer et al. (1998) showed that cosmological Alfvén waves generate a fractional CMB anisotropy for a Fourier mode k:(A.1)where n denotes sky direction, B is a unit vector in the direction of the coherent PMF, Ω(k,η_{last}) is the Gauge invariant linear combination associated with vector perturbations, η_{last} denotes the conformal time at the moment of baryonphoton decoupling, and T_{0} is the CMB monopole temperature 2.7255 K(Fixsen 2009). Durrer et al. (1998) assumed that the vector perturbations are initially created by some random stochastic PMF and have the following statistical properties over an ensemble of universes: (A.2)Here, P(k) is the power spectrum assumed to follow a simple power law (Durrer et al. 1998): (A.3)where k_{0} is a pivot wavenumber set to 0.05 / Mpc in this analysis. The Alfvén wave velocity is given by (Durrer et al. 1998): (A.4)where ρ_{r} and p_{r} are the density and the pressure of photons, and the speed of light c is set to 1.
Kahniashvili et al. (2008) showed that the presence of Alfvén waves in the early Universe leads to specific correlations of the CMB in harmonic space:
Here C_{ℓ} is the power spectrum in the absence of Alfvén waves, θ_{B} and φ_{B} are the spherical angles of a PMF direction B, and is given by where η_{0} is the present conformal time, and k_{D} denotes the comoving wavenumber of the dissipation scale, due to photon viscosity and given by approximately 10 /η_{last} (Durrer et al. 1998). The damping effect becomes significant on multipoles ℓ ≳ 500 (Durrer et al. 1998). As shown above, Alfvén waves in the early Universe produce correlations between harmonic modes separated by Δℓ = 0, ± 2, and Δm = 0, ± 1, ± 2. Investigating these imprints, we may impose a constraint on the Alfvén waves. In the weak Alfvén wave limit, the CMB data loglikelihood can be expanded as follows: (A.5)Since all correlations produced by Alfvén waves are proportional to , the first term in Eq. (A.5) is simply equal to the likelihood of the standard cosmological model. The first and second derivative of the likelihood are obtained by (A.6)where (A.7)Here a is the data vector, consisting of the spherical harmonic coefficients, a_{ℓm}, of the CMB anisotropy data, and C is their covariance matrix.
Planck constraints on the Alfvén wave amplitude .
In our analysis, we consider the four foreground cleaned CMB maps CommanderRuler, NILC, SEVEM, and SMICA and apply the common mask. We assume the fiducial Planck cosmological model, and use MC simulations to estimate the ensemble average values for signal and noise, as required in Eq. (A.6). The quantity C^{1}a from Eq. (A.7) was then determined by the messenger field method (Elsner & Wandelt 2013). The CosmoMC package (Lewis & Bridle 2002) is then used as a generic sampler in order to obtain the posterior probability for the Alfvén wave parameters .
In Table A.1, we show the upper bounds on the Alfvén wave amplitude at various confidence levels, after marginalizing over the spectral index n_{v} and the direction θ_{B},φ_{B}. From the analysis of the Planck data, we impose an upper bound on the Alfvén wave amplitude that is tighter than that from the WMAP data by more than one order of magnitude.
All Tables
Lower tail probablity for the variance, skewness, and kurtosis estimators at N_{side} = 2048, using the U73 mask and four different component separation methods.
Lower tail probablity for the variance, skewness, and kurtosis estimators at N_{side} = 2048, for the SMICA method, using different masks.
Lower tail probablity for the variance, skewness, and kurtosis estimators at N_{side} = 2048, for the SEVEM cleaned maps at different frequencies.
Lower tail probablity for the variance, skewness, and kurtosis estimators at different resolutions, for the four component separation methods, using the U73 mask.
Frequency dependence of the lower tail probablity for the Npdf, using different masks.
Probabilities of obtaining values of the χ^{2} statistic for the concordance ΛCDM model at least as large as the observed values of the statistic for the raw 143 GHz (first row) and foreground corrected 143 GHz SEVEM CMB maps (second row).
Probabilities of obtaining values for the χ^{2} statistic of the Npoint functions shown in Fig. 6 for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck CMB maps with resolution parameter N_{side} = 64, estimated using the CR, NILC, SEVEM and SMICA methods.
Probabilities of obtaining values of the χ^{2} statistic of the Npoint functions shown in Figs. 7 and 8 for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for Planck CMB maps with resolution parameter N_{side} = 512. estimated using the CR, NILC, SEVEM and SMICA methods.
Probabilities of obtaining values of the χ^{2} statistic of the Npoint functions shown in Fig. 9 for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for Planck CMB maps with resolution parameter N_{side} = 2048 estimated using the CR, NILC, SEVEM, and SMICA methods.
Nondirectional Gaussianity tests using unnormalized MFs: as a function of sky resolution for different component separation methods.
Nondirectional Gaussianity tests using normalized MFs: dependence of on resolution.
Nondirectional Gaussianity tests using normalized MFs: dependence of on component separation method.
Nondirectional Gaussianity tests using normalized MFs: dependence of on sky coverage.
χ^{2} statistics based on the wavelet bispectrum reconstruction y_{i} statistics for the foreground cleaned data map.
Lower tail probablity for the variance, skewness, and kurtosis at N_{side} = 16, using different masks.
Frequency dependence of the lower tail probablity for the variance skewness and kurtosis at N_{side} = 16, using different masks.
Probabilities of obtaining values of the χ^{2} statistic for the Planck fiducial model at least as large as the observed values of the statistic for the Planck maps with resolution parameter N_{side} = 64, estimated using the CR, NILC, SEVEM and SMICA methods.
Probabilities of obtaining values of the χ^{2} statistic for the Planck fiducial model at least as large as the observed values for the SMICA map at N_{side} = 64.
Probabilities of obtaining values of the χ^{2} statistic for the Planck fiducial model at least as large as the observed values of the statistic for the PlanckSMICA map with resolution parameter N_{side} = 64 and removed multipoles of order ℓ ≤ 3 (columns denoted as “ℓ > 3”) or order ℓ ≤ 5 (columns denoted as “ℓ > 5”) .
Deviations S and empirical probabilities p for the scaleindependent χ^{2}statistics derived from the CR, NILC, SEVEM, and SMICA maps.
Summary of the bestfit power dipole directions on the sky up to ℓ_{max} = 600/1500, for each of the four component separation methods, together with the 143 GHz map cleaned using the SEVEM method and denoted SEVEM143.
Significance of the correlations of dipole directions between high and low multipoles.
Summary of dipole modulation likelihood results as a function of scale for all four Planck CMB solutions.
Upper tail probability (UTP, in %) associated to the cold (top) and hot (bottom) areas.
Upper tail probability (in %) associated to the cold (top) and hot (bottom) areas.
Upper tail probabilities (in %) associated with the cold (top) and hot (bottom) areas.
All Figures
Fig. 1 Variance, skewness, and kurtosis for the combined map of the four different component separation methods. From top row to bottom row CR, NILC, SEVEM, SMICA. 

In the text 
Fig. 2 Npdf χ^{2} for the U73 mask, CL58, CL37, ecliptic North and ecliptic South. The different colours represent the four component separation methods, namely CR (green), NILC (blue), SEVEM (red), and SMICA (orange). 

In the text 
Fig. 3 Frequency dependence for 70 GHz (green), 100 GHz (blue), 143 GHz (red), and 217 GHz (orange), for different masks. 

In the text 
Fig. 4 Two sets of discs of radius 10° A (top) and B (middle), for the N_{side} = 512 CMB estimates; and a set of 20 randomly placed discs of radius 3° superimposed on the U73 mask (blue region) for the CMB estimates at N_{side} = 2048 (bottom). 

In the text 
Fig. 5 Pseudocollapsed 3point function averaged over disc set A for the raw (upper figure) and SEVEM foreground corrected (lower figure) 143 GHz map at N_{side} = 512 for different values of sky coverage (f_{sky}). Estimates of the multipoles for ℓ ≤ 18 are removed from the sky maps. The black solid line indicates the mean for 1000 MC simulations and the shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively, for the CG90 (f_{sky} = 0.9) mask. See Sect. 4.3 for the definition of the separation angle θ. 

In the text 
Fig. 6 The 2point (upper left), pseudocollapsed (upper right), equilateral (lower left) 3point and reduced rhombic 4point (lower right) functions for the N_{side} = 64 CMB estimates. The black solid line indicates the mean for 1000 Monte Carlo simulations and the shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 

In the text 
Fig. 7 The 2point (upper left), pseudocollapsed (upper right), equilateral (lower left) 3point and reduced rhombic 4point (lower right) functions averaged over disc set A for the N_{side} = 512 CMB estimates. Estimates of the multipoles for ℓ ≤ 18 are removed from the sky maps. The black solid line indicates the mean for 1000 Monte Carlo simulations and the shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 

In the text 
Fig. 8 As Fig. 7 for disc set B. 

In the text 
Fig. 9 The 2point (upper left), pseudocollapsed (upper right), equilateral (lower left) 3point and reduced rhombic 4point (lower right) functions averaged over the disc set for the N_{side} = 2048 CMB estimates. Estimates of the multipoles for ℓ ≤ 60 are removed from the sky maps. The black solid line indicates the mean for 1000 Monte Carlo simulations and the shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 

In the text 
Fig. 10 Difference of the data MFs (unnormalized) with respect to the average of the curves obtained with realistic Planck simulations for several cleaned maps. From top to bottom: Area, Contour, Genus. The errorbars represent the ±1σ (68% CL) dispersions around the mean obtained with simulations. 

In the text 
Fig. 11 Difference of the normalized MFs obtained from the data with respect to the expected values of the null hypothesis, for the different component separation methods. From left to right and top to bottom: area, Contour, Genus and N_{cluster}. The grey bands represent the ±1 and ±2σ dispersions around zero, based on realistic Planck simulations including lensing, for the CR method. 

In the text 
Fig. 12 Difference of the normalized MFs obtained from the data with respect to the expected value of the null hypothesis for several sky coverages. The SMICA map is considered. From left to right and top to bottom: area, Contour, Genus and N_{cluster}. The grey bands represent the 1 and 2σ dispersions around zero, based on realistic Planck simulations including lensing, for f_{sky} = 0.23. 

In the text 
Fig. 13 Standard deviation (left), skewness (centre) and kurtosis (right) of the SMHW coefficients as a function of the wavelet scale R. Results are given for the four Planck CMB maps (green, CommanderRuler; lightblue, NILC; red, SEVEM; and yellow, SMICA). Modified upper tail probabilities (mUTP, see text for details) are obtained by comparing with 1000 simulations processed through the component separation pipelines. Squares represent modified upper tail probabilities that correspond to an actual upper tail probability above 0.5; diamonds represent upper tail probabilities below 0.5. 

In the text 
Fig. 14 Standard deviation (left), skewness (centre) and kurtosis (right) of the SMHW coefficients as a function of the wavelet scale R. Results are given for the SMICA CMB map. Several masking scenarios are compared: red, CG60 mask (cutting out 40% of the sky); green, CG70 mask (cutting out 30% of the sky); and blue, U73 mask. The modified upper tail probabilities (mUTP) are defined in the text. 

In the text 
Fig. 15 Comparison between the WMAP sevenyear bispectrum signal (left) (Fergusson et al. 2010) and the lowℓ signal of Planck (right) reconstructed from the SMICA foregroundseparated map (in both cases using polynomial modes). The same basic patterns are observed in both bispectra, including an apparent central “oscillatory” feature. 

In the text 
Fig. 16 Wavelet bispectrum reconstruction y_{i} statistics for the foreground cleaned Planck data map. Considered data map: combined map cleaned with CR, NILC, SEVEM and SMICA. The solid yellow lines represent the 3σ errorbars for SMICA (similar errorbars are obtained for CR, NILC, and SEVEM maps). 

In the text 
Fig. 17 Upper: Wienerfiltered SMICA CMB sky (temperature range ±400 μK). Middle: derived quadrupole (temperature range ±35 μK). Lower: derived octopole (temperature range ±35 μK). The plus and star symbols indicate the axes of the quadrupole and octopole, respectively, around which the angular momentum dispersion is maximized. The diamond symbols correspond to the quadrupole axes after correction for the kinematic quadrupole. 

In the text 
Fig. 18 Variance, skewness and kurtosis at N_{side} = 16, for the U73 mask, CL58, CL37, ecliptic North, and ecliptic South (from top to bottom). The different lines represent the four component separation methods CR (green), NILC (blue), SEVEM (red), and SMICA (orange). 

In the text 
Fig. 19 Variance, skewness and kurtosis at N_{side} = 16, for the U73 mask, CL58, CL37, ecliptic North, and ecliptic South (from top to bottom). The different lines represent the four considered frequencies, namely 70 GHz (green), 100 GHz (blue), 143 GHz (red), and 217 GHz (orange). 

In the text 
Fig. 20 The 2point (upper left), pseudocollapsed (upper right), equilateral 3point (lower left), and rhombic 4point (lower right) correlation functions (N_{side} = 64). Correlation functions are shown for the analysis performed on northern (blue) and southern (red) hemispheres determined in the Ecliptic coordinate frame. The shaded dark and light grey bands indicate the 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 

In the text 
Fig. 21 The 2point (upper left), pseudocollapsed (upper right), equilateral 3point (lower left), and rhombic 4point (lower right) correlation functions for the SMICA map (N_{side} = 64). Correlation functions are shown for the analysis performed on negative (blue) and positive (red) hemispheres determined in the Ecliptic, Galactic, Doppler boost (DB), dipole modulation (DM), and power asymmetry (PA) coordinate frames. The shaded dark and light grey bands indicate the 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 

In the text 
Fig. 22 Deviations S( ⟨ α ⟩) of the rotated hemispheres derived from the mean of the scaling indices for the scale r_{5} determined from the SMICA map. Upper panel: result for the calculations with respect to the Galactic coordinate system. Lower panel: result obtained after averaging over 48 rotated coordinate systems. 

In the text 
Fig. 23 Deviations S(χ^{2}) of the rotated hemispheres derived from a combination of the mean and the standard deviation of the scaling indices for the scale r_{5} determined from the SMICA map. 

In the text 
Fig. 24 Deviations  S(r)  for the SMICA map as a function of the scale parameter r for the full sky (black) and positive (red) and negative (blue) rotated hemisphere. The reference frame for defining the positive and negative hemisphere is chosen such that the difference ΔS = S_{pos} − S_{neg} becomes maximal for ⟨ α(r_{5}) ⟩. The plus signs denote the results for the mean ⟨ α(r_{k}) ⟩, the starsigns represent the standard deviation σ_{α(rk)}. The dashed (dotted) line indicates the 1 (3) σ significance level. 

In the text 
Fig. 25 Probability density P( ⟨ α(r_{5}) ⟩) for the secondorder surrogates using the SMICA map for the positive (red) and negative (blue) rotated hemispheres. The vertical lines show the respective values for ⟨ α(r_{5}) ⟩ for the firstorder surrogates. 

In the text 
Fig. 26 Dipole directions for 100multipole bins of the local power spectrum distribution from ℓ = 2 to 1500 in the SMICA map with the apodized U73 mask applied. We also show the total direction for ℓ_{max} = 600 determined from WMAP9, as well as the preferred dipolar modulation axis (labelled as lowℓ) derived in Sect. 5.6. 

In the text 
Fig. 27 Dipole directions for 100multipole bins of the local power spectrum distribution from ℓ = 2 to 1500 in the deboosted SEVEM 143 GHz map (SEVEM143DB) with the apodized U73 mask applied. We also show the direction for ℓ_{max} = 600 determined from WMAP9, as well as the preferred dipolar modulation axis (labelled as lowℓ) derived in Sect. 5.6. 

In the text 
Fig. 28 Derived pvalues as a function of ℓ_{max}, for the SEVEM 143 GHz foregroundcleaned map, determined both before and after deboosting. The pvalues are computed using 500 FFP6 simulated maps. 

In the text 
Fig. 29 Derived pvalues as a function of ℓ_{max} for the deboosted SEVEM 143 GHz foregroundcleaned map. The pvalues in this plot are based on the the mean of the angular separations determined between all pairs of dipole directions where one direction falls in the range [ℓ_{lim} = 500,ℓ_{max}] and the second direction in the range [2,ℓ_{lim} = 500]. The significance is computed using 500 FFP6 simulated maps. 

In the text 
Fig. 30 Fractional power ratio, ΔC_{ℓ}/C_{ℓ}, from antipodal sky regions, computed from the SEVEM 143 GHz map before and after deboosting along the mean dipole direction for ℓ = 2 to 600. All spectra are evaluated on hemispheres using an apodized mask. The grey lines show the same quantity evaluated for each of the 500 FFP6 simulations along their respective asymmetry axes. The green band shows the 68% confidence region from these simulations. 

In the text 
Fig. 31 Marginalized dipole modulation amplitude (top), power spectrum amplitude (middle) and power spectrum tilt (bottom) probability distributions as a function of smoothing scale, shown for the Commander CMB solution. 

In the text 
Fig. 32 Consistency between component separation algorithms as measured by the dipole modulation likelihood. The top panel shows the marginal power spectrum amplitude for the 5° smoothing scale, the middle panel shows the dipole modulation amplitude, and the bottom panel shows the preferred dipole directions. The coloured area indicates the 95% confidence region for the Commander solution, while the dots shows the maximumposterior directions for the other maps. 

In the text 
Fig. 33 Loglikelihood difference between the bestfit dipole modulation model and the fiducial isotropic model as a function of smoothing scale. Horizontal dashed lines indicate 1, 2, and 3σ thresholds. 

In the text 
Fig. 34 Significance of the modulation power, L(L + 1)m_{L}/ 2π, at bipolar multipoles L. The modulation spectra obtained from the four component separation maps (CR, NILC, SEVEM, and SMICA) are consistent with each other. Dipole (L = 1) modulation power is detected in all the spectra at a significance ranging from 2.9 to 3.7σ. The solid black lines denote the 3σ significance thresholds. There is no significant power detected at higher multipole of the modulation field, 1 <L ≤ 32. 

In the text 
Fig. 35 Measured dipole modulation (L = 1) power in CMB multipole bins. This uses the CMB multipole dependence of the BipoSH (modulation) power L(L + 1)m_{L}/ 2π, which can be dissected into bins in ℓspace. We establish that significant power in the dipole modulation is limited to ℓ = 2–64 and does not extend to the higher CMB multipoles considered. The vertical grid lines denote the CMB multipole ℓbins. 

In the text 
Fig. 36 Dipolar (L = 1) power measured from different CMB multipole bins, as denoted by the vertical grid lines. Here the Doppler estimator is used to reconstruct the Doppler field () from SMICA. Significant dipolar power is limited to ℓ = 2–64 and does not extend to the higher CMB multipoles considered. 

In the text 
Fig. 37 Reconstruction noise power associated with the modulation (blue) and Doppler (red) estimators as a function of the maximum CMB multipole used. The solid lines correspond to the case of noisefree CMB maps, while the dashed lines correspond to the case where the characteristic Planck instrument noise is taken into account. The red dashdotted horizontal line denotes the expected signal level of the boost due to our local motion, while the blue dashdotted horizontal line denotes the enhanced boost amplitude as seen by the modulation estimator. The blue dotted line denotes the effective Doppler boost signal as seen by the modulation estimator after accounting for the effects of aberration induced by the boost. Also note that the reconstruction errors shown in red are derived for the total dipole power, rather than the amplitude in a particular direction, and hence are higher by a factor of as compared to the errors for individual directions reported in Planck Collaboration XXVII (2014). 

In the text 
Fig. 38 BipoSH power spectra and obtained from the four component separation maps (CR, NILC, SEVEM, and SMICA) are consistent with each other. Note that no significant detections are found. This independently establishes the fact that the quadrupolar BipoSH detections made by WMAP were due to WMAPspecific instrument systematics. 

In the text 
Fig. 39 Upper: the parity estimator g(ℓ) versus ℓ for CR (black diamonds), NILC (red diamonds), SEVEM (blue diamonds), and SMICA (green diamonds). Lower: pvalue for CR (black solid line), NILC (red line), SEVEM (blue line), and SMICA (green line). 

In the text 
Fig. 40 Top panel: S^{+} statistic. The vertical lines show the minimum value for the estimator as computed on low resolution CR, NILC, SEVEM, and SMICA maps. The grey histogram shows the same quantity computed from 1000 simulated maps processed by CR. Bottom panel: as above for S^{−}. 

In the text 
Fig. 41 SMHW coefficients at R = 300 arcmin, above (and below) a + 3.0σ (−3.0σ) threshold. Results for the three masks considered in the analysis are shown: U73 mask (top); CG70 (middle); and CG60 (bottom). 

In the text 
Fig. 42 Same as Fig. 23 but with the best fit Bianchi template subtracted from the SMICA map. 

In the text 
Fig. 43 The angular power spectrum on large angular scales computed from the SMICA map at N_{side} = 16 on the opposing hemispheres defined by the maximal asymmetry direction (l,b) = (220°,4°) determined in Sect. 5.5. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.