Issue 
A&A
Volume 594, October 2016
Planck 2015 results



Article Number  A16  
Number of page(s)  62  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201526681  
Published online  20 September 2016 
Planck 2015 results
XVI. 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, 6–8 Melrose Road, Muizenberg, Cape Town, South Africa
^{4} Agenzia Spaziale Italiana Science Data Center, via del Politecnico snc, 00133 Roma, Italy
^{5} Aix Marseille Université, CNRS, LAM (Laboratoire d’Astrophysique de Marseille) UMR 7326, 13388 Marseille, France
^{6} Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{7} Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & 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} CRANN, Trinity College, Dublin, Ireland
^{11} California Institute of Technology, Pasadena, California CA 91125, USA
^{12} Centre for Theoretical Cosmology, DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
^{13} Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan, 1, planta 2, 44001 Teruel, Spain
^{14} Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{15} Consejo Superior de Investigaciones Científicas (CSIC), Madrid, Spain
^{16} DSM/Irfu/SPP, CEASaclay, 91191 GifsurYvette Cedex, France
^{17} DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
^{18} Département de Physique Théorique, Université de Genève, 24 Quai E. Ansermet, 1211 Genève 4, Switzerland
^{19} Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{20} Departamento de Física, Universidad de Oviedo, Avda. Calvo Sotelo s/n, 33003 Oviedo, Spain
^{21} Departamento de Matemáticas, Estadística y Computación, Universidad de Cantabria, Avda. de los Castros s/n, 39005 Santander, Spain
^{22} Department of Astronomy and Astrophysics, University of Toronto, 50 Saint George Street, Toronto, Ontario, Canada
^{23} Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{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, 00100 Helsinki, Finland
^{29} Department of Physics, Princeton University, Princeton, New Jersey NJ 08544, USA
^{30} Department of Physics, University of California, Santa Barbara, California CA 93106, USA
^{31} Department of Physics, University of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, Illinois, USA
^{32} Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131 Padova, Italy
^{33} Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{34} Dipartimento di Fisica, Università La Sapienza, P. le A. Moro 2, Roma, Italy
^{35} Dipartimento di Fisica, Università degli Studi di Milano, via Celoria, 16 Milano, Italy
^{36} Dipartimento di Fisica, Università degli Studi di Trieste, via A. Valerio 2, Trieste, Italy
^{37} Dipartimento di Matematica, Università di Roma Tor Vergata, via della Ricerca Scientifica, 1 Roma, Italy
^{38} Discovery Center, Niels Bohr Institute, 17 Blegdamsvej, Copenhagen, Denmark
^{39} European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain
^{40} European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{41} Facoltà di Ingegneria, Università degli Studi eCampus, via Isimbardi 10, 22060 Novedrate (CO), Italy
^{42} Gran Sasso Science Institute, INFN, viale F. Crispi 7, 67100 L’ Aquila, Italy
^{43} HGSFP and University of Heidelberg, Theoretical Physics Department, Philosophenweg 16, 69120 Heidelberg, Germany
^{44} Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland
^{45} INAF–Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, Padova, Italy
^{46} INAF–Osservatorio Astronomicodi Roma, via di Frascati 33, Monte Porzio Catone, Italy
^{47} INAF–Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, Trieste, Italy
^{48} INAF/IASF Bologna, via Gobetti 101, Bologna, Italy
^{49} INAF/IASF Milano, via E. Bassini 15, Milano, Italy
^{50} INFN, Sezione di Bologna, via Irnerio 46, 40126 Bologna, Italy
^{51} INFN, Sezione di Roma 1, Università di Roma Sapienza, Piazzale Aldo Moro 2, 00185 Roma, Italy
^{52} INFN, Sezione di Roma 2, Università di Roma Tor Vergata, via della Ricerca Scientifica, 1 Roma, Italy
^{53} INFN/National Institute for Nuclear Physics, via Valerio 2, 34127 Trieste, Italy
^{54} IPAG: Institut de Planétologie et d’Astrophysique de Grenoble, Université Grenoble Alpes, IPAG; CNRS, IPAG, 38000 Grenoble, France
^{55} IUCAA, Post Bag 4, Ganeshkhind, Pune University Campus, 411 007 Pune, India
^{56} Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK
^{57} Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA
^{58} Institut Néel, CNRS, Université Joseph Fourier Grenoble I, 25 rue des Martyrs, Grenoble, France
^{59} Institut Universitaire de France, 103 bd SaintMichel, 75005 Paris, France
^{60} Institut d’Astrophysique Spatiale, CNRS (UMR 8617) Université ParisSud 11, Bâtiment 121, 91440 Orsay, France
^{61} Institut d’Astrophysique de Paris, CNRS (UMR 7095), 98 bis Boulevard Arago, 75014 Paris, France
^{62} Institut für Theoretische Teilchenphysik und Kosmologie, RWTH Aachen University, 52056 Aachen, Germany
^{63} Institute for Space Sciences, BucharestMagurale, 077125 Romania
^{64} Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{65} Institute of Theoretical Astrophysics, University of Oslo, Blindern, 0371 Oslo, Norway
^{66} Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, La Laguna, 38205 Tenerife, Spain
^{67} Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, 39005 Santander, Spain
^{68} Istituto Nazionale di Fisica Nucleare, Sezione di Padova, via Marzolo 8, 35131 Padova, Italy
^{69} Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, USA
^{70} Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13 9PL, UK
^{71} Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA
^{72} Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
^{73} Kazan Federal University, 18 Kremlyovskaya St., 420008 Kazan, Russia
^{74} LAL, Université ParisSud, CNRS/IN2P3, 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 Cosmologie, Université GrenobleAlpes, CNRS/IN2P3, 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} Lebedev Physical Institute of the Russian Academy of Sciences, Astro Space Centre, 84/32 Profsoyuznaya st., GSP7, 117997 Moscow, Russia
^{82} Leung Center for Cosmology and Particle Astrophysics, National Taiwan University, 10617 Taipei, Taiwan
^{83} MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{84} McGill Physics, Ernest Rutherford Physics Building, McGill University, 3600 rue University, Montréal, QC, H3A 2T8, Canada
^{85} National University of Ireland, Department of Experimental Physics, Maynooth, Co. Kildare, Ireland
^{86} Nicolaus Copernicus Astronomical Center, Bartycka 18, 00716 Warsaw, Poland
^{87} Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark
^{88} Optical Science Laboratory, University College London, Gower Street, London, UK
^{89} SISSA, Astrophysics Sector, via Bonomea 265, 34136 Trieste, Italy
^{90} School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff, CF24 3AA, UK
^{91} School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{92} Simon Fraser University, Department of Physics, 8888 University Drive, Burnaby BC, Canada
^{93} Sorbonne UniversitéUPMC, UMR 7095, Institut d’Astrophysique de Paris, 98 bis Boulevard Arago, 75014 Paris, France
^{94} Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, 117997 Moscow, Russia
^{95} Space Sciences Laboratory, University of California, Berkeley, California, USA
^{96} Special Astrophysical Observatory, Russian Academy of Sciences, Nizhnij Arkhyz, Zelenchukskiy region, 369167 KarachaiCherkessian Republic, Russia
^{97} Stanford University, Dept of Physics, Varian Physics Bldg, 382 via Pueblo Mall, Stanford, California, USA
^{98} SubDepartment of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
^{99} Theory Division, PHTH, CERN, 1211, Geneva 23, Switzerland
^{100} UPMC Univ Paris 06, UMR 7095, 98 bis Boulevard Arago, 75014 Paris, France
^{101} Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex 4, France
^{102} University of Granada, Departamento de Física Teórica y del Cosmos, Facultad de Ciencias, 18071 Granada, Spain
^{103} University of Granada, Instituto Carlos I de Física Teórica y Computacional, 18071 Granada, Spain
^{104} University of Heidelberg, Institute for Theoretical Physics, Philosophenweg 16, 69120 Heidelberg, Germany
^{105} Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
^{⋆}
Corresponding author: A. J. Banday, email: anthony.banday@irap.omp.eu
Received: 5 June 2015
Accepted: 7 September 2015
We test the statistical isotropy and Gaussianity of the cosmic microwave background (CMB) anisotropies using observations made by the Planck satellite. Our results are based mainly on the full Planck mission for temperature, but also include some polarization measurements. In particular, we consider the CMB anisotropy maps derived from the multifrequency Planck data by several componentseparation methods. For the temperature anisotropies, we find excellent agreement between results based on these sky maps over both a very large fraction of the sky and a broad range of angular scales, establishing that potential foreground residuals do not affect our studies. Tests of skewness, kurtosis, multinormality, Npoint functions, and Minkowski functionals indicate consistency with Gaussianity, while a power deficit at large angular scales is manifested in several ways, for example low map variance. The results of a peak statistics analysis are consistent with the expectations of a Gaussian random field. The “Cold Spot” is detected with several methods, including map kurtosis, peak statistics, and mean temperature profile. We thoroughly probe the largescale dipolar power asymmetry, detecting it with several independent tests, and address the subject of a posteriori correction. Tests of directionality suggest the presence of angular clustering from large to small scales, but at a significance that is dependent on the details of the approach. We perform the first examination of polarization data, finding the morphology of stacked peaks to be consistent with the expectations of statistically isotropic simulations. Where they overlap, these results are consistent with the Planck 2013 analysis based on the nominal mission data and provide our most thorough view of the statistics of the CMB fluctuations to date.
Key words: cosmology: observations / cosmic background radiation / polarization / methods: data analysis / methods: statistical
© ESO, 2016
1. Introduction
This paper, one of a set associated with the 2015 release of data from the Planck^{1} mission (Planck Collaboration I 2016), describes a set of studies undertaken to determine the statistical properties of both the temperature and polarization anisotropies of the cosmic microwave background (CMB).
The standard cosmological model is described well by the FriedmannLemaîtreRobertsonWalker solution of the Einstein field equations. This model is characterized by a homogeneous and isotropic background metric and a scale factor of the expanding Universe. It is hypothesized that at very early times 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 approximately as a de Sitter space, providing the conditions by which some of its present properties 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 density perturbations that are distributed as a statistically homogeneous and isotropic Gaussian random field. Linear theory relates those perturbations to the temperature and polarization anisotropies of the CMB, implying a distribution for the anisotropies very close to that of a statistically isotropic Gaussian random field.
The aim of this paper is to use the full mission Planck data to test the Gaussianity and isotropy of the CMB as measured in both intensity and, in a more limited capacity, polarization. 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 statistically 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. Indeed, the isotropy and Gaussianity of the CMB anisotropies are implicitly assumed in critical science papers from the 2015 release, in particular those describing the likelihood and the derivation of cosmological parameter constraints (Planck Collaboration XI 2016; Planck Collaboration XIII 2016). Conversely, if the detection of significant deviations from these assumptions cannot be traced to known systematic effects or foreground residuals, the presence of which should be diagnosed by the statistical tests set forth in this paper, this would necessitate a major revision of the current methodological approaches adopted in deriving the mission’s many science results.
Wellunderstood physical processes due to the integrated SachsWolfe (ISW) effect (Planck Collaboration XVII 2014; Planck Collaboration XXI 2016) and gravitational lensing (Planck Collaboration XIX 2014; Planck Collaboration XV 2016) lead to secondary anisotropies that exhibit marked deviation from Gaussianity. In addition, Doppler boosting, due to our motion with respect to the CMB rest frame, induces 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, and were detected at a statistically significant level on small angular scales in Planck Collaboration XXVII (2014). Beyond these, Planck Collaboration XXIII (2014, hereafter PCIS13 established that the Planck 2013 data set showed little evidence for nonGaussianity, with the exception of a number of CMB temperature anisotropy anomalies on large angular scales that confirmed earlier claims based on WMAP data. Moreover, given that the broader frequency coverage of the Planck instruments allowed improved component separation methods to be applied in the derivation of foregroundcleaned CMB maps, it was generally considered that the case for anomalous features in the CMB had been strengthened. Hence, such anomalies have attracted considerable attention in the community, since they could be the visible traces of fundamental physical processes occurring in the early Universe.
However, the literature also supports an ongoing debate about the significance of these anomalies. The central issue in this discussion is connected with 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) base their rejection of the presence of anomalies in the CMB on such arguments. Of course, one should attempt to correct for any choices that were made in the process of detecting an anomaly. However, in the absence of an alternative model for comparison to the standard Gaussian, statistically isotropic one adopted to quantify significance, this is often simply not possible. In this work, whilst it is recognized that care must be taken in the assessment of significance, we proceed on the basis that allowing a posteriori reasoning permits us to challenge the limits of our existing knowledge (Pontzen & Peiris 2010). That is, by focusing on specific properties of the observed data that are shown to be empirically interesting, we may open up new paths to a better theoretical understanding of the Universe. We will clearly describe the methodology applied to the data, and attempt to study possible links among the anomalies in order to search for a physical interpretation.
The analysis of polarization data introduces a new opportunity to explore the statistical properties of the CMB sky, including the possibility of improvement of the significance of detection of largescale anomalies. However, this cannot be fully included in the current data assessment, since the componentseparation products in polarization are highpass filtered to remove large angular scales (Planck Collaboration IX 2016), owing to the persistence of significant systematic artefacts originating in the High Frequency Instrument (HFI) data (Planck Collaboration VII 2016; Planck Collaboration VIII 2016). In addition, limitations of the simulations with which the data are to be compared (Planck Collaboration XII 2016), in particular a significant mismatch in noise properties, limit the extent to which any polarization results can be included. Therefore, we only present a stacking analysis of the polarized data, although this is a significant extension of previous approaches found in the literature.
With future Planck data releases, it will be important to determine in more detail whether there are any pecularities in the CMB polarization, and if so, whether they are related to existing features in the CMB temperature field. Conversely, the absence of any corresponding features in polarization might imply that the temperature anomalies (if they are not simply flukes) could be due to a secondary effect such as the ISW effect, or alternative scenarios in which the anomalies arise from physical processes that do not correlate with the temperature, e.g., texture or defect models. Either one of these possible outcomes could yield a breakthrough in understanding the nature of the CMB anomalies. Of course, there also remains the possibility that anomalies may be found in the polarization data that are unrelated to existing features in the temperature measurements.
Following the approach established in Planck Collaboration XXIII (2014), throughout this paper we quantify the significance of a 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., primordial Gaussianity and isotropy of the CMB) is true. In some tests, where it is clearly justified to only use a onetailed probability, the pvalue is replaced by the corresponding upper or lowertail probability.
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. Specific theoreticallymotivated model constraints on isotropy or nonGaussianity, as might arise from nonstandard inflationary models, the geometry and topology of the Universe, and primordial magnetic fields are provided in the companion papers (Planck Collaboration XVII 2016; Planck Collaboration XX 2016; Planck Collaboration XVIII 2016; Planck Collaboration XIX 2016). The paper is organized as follows. Section 2 summarizes the Planck full mission data used for the analyses, and important limitations of the polarization maps that are studied. Section 3 describes the 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. Several important anomalous features of the CMB sky, originally detected with the WMAP data and subsequently confirmed in PCIS13, are reassessed in Sect. 5. Aspects of the CMB fluctuations specifically related to dipolar asymmetry are examined in Sect. 6. The sensitivity of the results for a number of statistical tests to the sky fraction is examined in Sect. 7. Section 8 presents tests of the statistical nature of the polarization signal observed by Planck using a local analysis of stacked patches of the sky. Finally, Sect. 9 provides the main conclusions of the paper.
2. Data description
In this paper, we use data from the Planck2015 full mission data release. This contains approximately 29 months of data for the HFI and 50 months for the Low Frequency Instrument (LFI). The release includes sky maps at nine frequencies in intensity (seven in polarization), with corresponding “halfmission” maps that are generated by splitting the fullmission data sets in various ways. The maps are provided in HEALPix format (Górski et al. 2005)^{2}, 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 our analyses and limits on timevarying systematic artefacts. Full details are provided in a series of companion papers (Planck Collaboration II 2016; Planck Collaboration III 2016; Planck Collaboration IV 2016; Planck Collaboration V 2016; Planck Collaboration VI 2016; Planck Collaboration VII 2016; Planck Collaboration VIII 2016).
Our main results are based on estimates of the CMB generated by four distinct componentseparation algorithms – Commander, NILC, SEVEM, and SMICA – as described in Planck Collaboration IX (2016). These effectively combine the raw Planck frequency maps in such a way as to minimize foreground residuals from diffuse Galactic emission. Note that the additional information in the full mission data set allows us to improve the reconstruction noise levels by roughly a factor of 2 (in temperature) as compared to the Planck2013 nominal mission data release. The CMB intensity maps were derived using all channels, from 30 to 857 GHz, and provided at a common angular resolution of 5′ FWHM and N_{side} = 2048. The intensity maps are only partially corrected for the second order temperature quadrupole (Kamionkowski & Knox 2003). Therefore, where appropriate, the componentseparated maps should be corrected for the residual contribution (Notari & Quartin 2015), specifically as described in Planck Collaboration IX (2016). The polarization solutions include all channels sensitive to polarization, from 30 to 353 GHz, at a resolution of 10′ FWHM and N_{side} = 1024. Possible residual emission is then mitigated in our analyses by the use of skycoverage masks, provided for both intensity and polarization.
Since in some cases it is important to study the frequency dependence of the cosmological signal, either to establish its primordial origin or to test for the frequency dependence associated with specific effects such as Doppler boosting (see Sect. 6.4), we also consider the foregroundcleaned versions of the 100, 143, and 217 GHz sky maps generated by the SEVEM algorithm (Planck Collaboration IX 2016), hereafter referred to as SEVEM100, SEVEM143, and SEVEM217, respectively.
For the present release, a postprocessing highpassfiltering has been applied to the CMB polarization maps in order to mitigate residual largescale systematic errors in the HFI channels (Planck Collaboration VII 2016). The filter results in the elimination of structure in the maps on angular scales larger than about 10°, and a weighted suppression of power down to scales of 5°, below which the maps remain unprocessed.
Lowerresolution versions of these data sets are also used in the analyses presented in this paper. The downgrading procedure for maps is to decompose them into spherical harmonics on the full sky at the input HEALPix resolution. The spherical harmonic coefficients, a_{ℓm}, are then convolved to the new resolution using (1)where b_{ℓ} is the beam transfer function, p_{ℓ} is the HEALPix pixel window function, and the “in” and “out” superscripts denote the input and output resolutions. They are then synthesized into a map directly at the output HEALPix resolution. Masks are downgraded in a similar way. The binary mask at the starting resolution is first downgraded like a temperature map. The smooth downgraded mask is then thresholded by setting pixels where the value is less than 0.9 to zero and all others to unity in order to make a binary mask. Table 1 lists the N_{side} and FWHM values defining the resolution of these maps, together with the different masks and their sky coverages that accompany the signal maps. In general, we make use of standardized masks that are the union of those associated with the individual componentseparation methods.
As recommended in Planck Collaboration IX (2016), the mask UT78 is adopted for all highresolution analyses of temperature data. UTA76 is an extended version of this mask more suitable for some nonGaussianity studies. The mask preferred for polarization studies, UPB77, is again the union of those determined for each component separation method, but in addition the polarized point sources detected at each frequency channel are excluded. These masks are then downgraded for lowerresolution studies. As a consequence of the common scheme applied in order to generate such lowresolution masks, they are generally more conservative than the corresponding ones used in the 2013 analyses.
In what follows, we will undertake analyses of the data at a given resolution denoted by a specific N_{side} value. Unless otherwise stated, this implies that the data have been smoothed to a corresponding FWHM as described above, and a standardized mask employed. Irrespective of the resolution in question, we will then often simply refer to the latter as the “common mask”.
Standardized data sets used in this paper.
3. Simulations
The results presented in this paper are derived using the extensive full focal plane (FFP8) simulations described in Planck Collaboration XII (2016). Of most importance are the Monte Carlo (MC) simulations that provide the reference set of Gaussian sky maps used for the null tests employed here. They also form the basis of any debiasing in the analysis of the real data as required by certain statistical methods.
The simulations include both CMB signal and instrumental noise realizations that capture important characteristics of the Planck scanning strategy, telescope, detector responses, and data reduction pipeline over the full mission period. In particular, the signal realizations include FEBeCoP (Mitra et al. 2011) beam convolution at each of the Planck frequencies, and are propagated through the various componentseparation pipelines using the same weights as derived from the Planck full mission data analysis.
The FFP8 fiducial CMB power spectrum has been adopted from our best estimate of the cosmological parameters from the first Planck data release (Planck Collaboration I 2014). This corresponds to a cosmology with baryon density given by ω_{b} = Ω_{b}h^{2} = 0.0222, cold dark matter (CDM) density ω_{c} = Ω_{c}h^{2} = 0.1203, neutrino energy density ω_{ν} = Ω_{ν}h^{2} = 0.00064, density parameter for the cosmological constant Ω_{Λ} = 0.6823, Hubble parameter H_{0} = 100h km s^{1} Mpc^{1} with h = 0.6712, spectral index of the power spectrum of the primordial curvature perturbation n_{s} = 0.96, and amplitude of the primordial power spectrum (at k = 0.05 Mpc^{1}) A_{s} = 2.09 × 10^{9}, and with the Thomson optical depth through reionization defined to be τ = 0.065. Each realization of the CMB sky is generated including lensing, Rayleigh scattering, and Doppler boosting effects, the latter two of which are frequencydependent. Unfortunately, the aberration contribution to the Doppler boost was erroneously omitted from the simulations, but, with possible exceptions described in Sect. 6, this does not lead to any significant impact on the results in this paper. A second order temperature quadrupole (Kamionkowski & Knox 2003) is added to each simulation with an amplitude corresponding to the residual uncorrected contribution present in the observed data, as described in Planck Collaboration XII (2016).
However, the Planck maps were effectively renormalized by approximately 2% to 3% in power in the time between the generation of the FFP8 simulations and the final maps. As discussed in Planck Collaboration XII (2016), correction for this calibration effect should have no significant impact on cosmological parameters. As recommended, in this paper the CMB component of the simulations is simply rescaled by a factor of 1.0134 before analysis.
Of somewhat more importance is an observed noise mismatch between the simulations and the data. Whilst this has essentially no impact on studies of temperature anisotropy, it imposes important limitations on the statistical studies of polarization sky maps that can be included here. Conversely, analyses based on 1point statistics, such as the variance, and the Npoint correlation functions have played important roles in establishing the nature of this mismatch, which seems to be scaledependent with an amplitude around 20% at lower resolutions but falling to a few per cent at higher resolution. As a consequence, this paper only includes results from a stacking analysis of the polarized data, in which the stacking of the data themselves necessarily acts to lower the effect of the noise mismatch. Polarization studies that do not rely on autostatistics can still yield interesting new results, as found in Planck Collaboration XIII (2016); Planck Collaboration XVII (2016); Planck Collaboration XVIII (2016).
4. Tests of nonGaussianity
There is no unique signature of nonGaussianity, but the application of a variety of tests over a range of angular scales allows us to probe the data for departures from theoretically motivated Gaussian statistics. One of the more important tests in the context of inflationary cosmology is related to the analysis of the bispectrum. This is discussed thoroughly in Planck Collaboration XVII (2016), and is therefore not discussed further in this paper. In this section, we present the results from a variety of statistical tools. Unless otherwise specified, the analyses are applied to all four component separation products (Commander, NILC, SEVEM, and SMICA) at a given resolution with the accompanying common mask, and significance levels are determined by comparison with the corresponding results derived from the FFP8 simulations, with typically 1000 being used for this purpose. Establishing the consistency of the results derived from the different componentseparation techniques is essential in order to be able to make robust claims about the statistical nature of the observed temperature fluctuations, and potential deviations from Gaussianity.
4.1. Onedimensional moments
In this section we consider simple tests of Gaussianity based on the variance, skewness, and kurtosis of the CMB temperature maps. Previous analyses found an anomalously low variance in the WMAP sky maps (Monteserín et al. 2008; Cruz et al. 2011), which was subsequently confirmed in an analysis of the Planck 2013 data (PCIS13).
Cruz et al. (2011) developed the unit variance estimator to determine the variance, , of the CMB signal on the sky in the presence of noise. The normalized CMB map, u^{X}, is given by (2)where X_{i} is the observed temperature at pixel i and is the noise variance for that pixel. Although this estimator is not optimal, Cruz et al. (2011) and Monteserín et al. (2008) have demonstrated that it is unbiased and sufficiently accurate for our purposes. The noise variance is estimated from the noise simulations for each componentseparation algorithm. The CMB variance is then estimated by requiring that the variance of the normalized map u^{X} is unity. The skewness and kurtosis can then be obtained from the appropriately normalized map.
Figure 1 presents results for the variance, skewness, and kurtosis determined from the data at a resolution of 5′, N_{side} = 2048. Good agreement between the component separation techniques is found, with small discrepancies likely due to sensitivity to the noise properties and their variation between methods.
Fig. 1 Variance, skewness, and kurtosis for the four different componentseparation methods – Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) – compared to the distributions derived from 1000 Monte Carlo simulations. 
Lowertail probabilities for the variance, skewness, and kurtosis of the componentseparated maps.
Table 2 summarizes the lowertail probabilities, defined as the percentage of MC simulations that show a lower variance, skewness, or kurtosis than the observed map, for these analyses. The results are in good agreement with PCIS13; the skewness and kurtosis are compatible with simulations, but the variance is marginally lower than in the simulations.
Although the variance is observed to be low, the results could still be affected by the presence of residual foregrounds at small scales in these maps, so that the true variance would be lower still. We assess this by application of the estimator to the cleaned frequency maps SEVEM100, SEVEM143, and SEVEM217. The results, also presented in Table 2, are similar to those found for the combined map, although slightly less significant, which is most likely attributable to higher noise in the cleaned frequency maps.
In conclusion, a simple statistical assessment of the Planck 2015 data using skewness and kurtosis shows no evidence for nonGaussianity, although a low variance is found, which we will readdress in Sect. 5.1.
4.2. Testing the multinormality of the CMB
Under the assumption of Gaussianity, the probability density function (PDF) of the Ndimensional pixelized temperature map is given by a multivariate Gaussian function: (3)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, and C is the covariance of the Gaussian field (of size N_{pix} × N_{pix}).
Although the calculation of TC^{1}T^{T} can be achieved by conjugate gradient methods, the evaluation of detC remains computationally difficult 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_{ℓ}): (4)where C_{ij} is the covariance between pixels i and j, θ_{ij} is the angle between them, P_{ℓ} are the Legendre polynomials, b_{ℓ} is an effective window function describing the combined effects of the instrumental beam and pixel window at resolution N_{side}, and ℓ_{max} is the maximum multipole probed.
Under the multivariate Gaussian hypothesis, the argument of the exponential in Eq. (3) should follow a χ^{2} distribution with N_{pix} degrees of freedom, or, equivalently (for N_{pix} ≫ 1) a normal distribution .
These χ^{2} statistics are computed for the Planck 2015 componentseparated CMB maps at N_{side} = 16 and 32, then compared with the equivalent quantities derived from the corresponding FFP8 simulations. For those cases in which the covariance matrix is illconditioned, we use a principal component analysis approach to remove the lowest degenerate eigenvalues of the covariance matrix (see, e.g., Curto et al. 2011). This process is equivalent to adding uncorrelated regularization noise of amplitude ≈ 1 μK to the data before inversion. The results of the analysis are presented in Table 3 and indicate that the data are consistent with Gaussianity. We note that the lowertail probabilities for the Npdf decrease when the resolution of the data is increased from N_{side} = 16 to 32. However, this behaviour is consistent with that seen for simulations, and should not be considered to be significant.
4.3. Npoint correlation functions
In this section, we present tests of the nonGaussianity of the Planck 2015 temperature CMB maps using realspace Npoint correlation functions. While harmonicspace methods are often preferred over realspace methods for studying primordial fluctuations, realspace methods have an advantage with respect to systematic errors 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 defined as the average product of N temperatures, measured in a fixed relative orientation on the sky, (5)where the unit vectors span an Npoint polygon. Under the assumption of statistical isotropy, these functions depend only on 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.
Lowertail probabilities for the Npdf χ^{2} statistics derived from the Planck 2015 componentseparated maps at N_{side} = 16 and 32.
The correlation 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, (6)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.
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., an 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, θ, parameterizes its size. Analogously, in the case of the equilateral triangle and rhombus, the size of the polygon is parameterized by the length of the edge, θ. Note that these functions are chosen for ease of implementation, not because they are better suited for testing Gaussianity than other configurations. For a Gaussian field, Wick’s theorem (Wick 1950) means that the ensemble average of the 4point function may be written in terms of the 2point function. In the following, all results refer to the connected 4point function, i.e., are corrected for this Gaussian contribution.
We use a simple χ^{2} statistic to quantify the agreement between the observed data and simulations, defined by (7)Here, Ĉ_{N}(θ_{i}) is the Npoint correlation function for the bin with separation angle θ_{i}, ⟨C_{N}(θ_{i})⟩ is the corresponding average from the MC simulation ensemble, and N_{bin} is the number of bins used for the analysis. If is the kth simulated Npoint correlation function and N_{sim} is the number of simulations, then the covariance matrix M_{ij} is given by (8)where . Following Hartlap et al. (2007), we then correct for bias in the inverse covariance matrix by multiplying it by the factor . 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 CMB estimates at a resolution of N_{side} = 64, this being constrained by computational limitations. The results are presented in Fig. 2, where we compare the Npoint functions for the data and the mean values estimated from 1000 MC simulations. 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 Table 4.
Fig. 2 Npoint correlation functions determined from the N_{side} = 64Planck CMB 2015 temperature maps. Results are shown for the 2point, pseudocollapsed 3point (upper left and right panels, respectively), equilateral 3point, and connected rhombic 4point functions (lower left and right panels, respectively). The red dotdotdotdashed, orange dashed, green dotdashed, and blue long dashed lines correspond to the Commander, NILC, SEVEM, and SMICA maps, respectively. Note that the lines lie on top of each other. The black solid line indicates the mean determined from 1000 SMICA simulations. The shaded dark and light grey regions indicate the corresponding 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 
Probabilities of obtaining values for the χ^{2} statistic of the Npoint functions for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2015 temperature CMB maps with resolution parameter N_{side} = 64, estimated using the Commander, NILC, SEVEM, and SMICA methods.
It is worth noting that the values of the Npoint functions for different angular separations are strongly correlated, and for this reason the figures show only one profile of each function in multidimensional space. Since the estimated probabilities take into account the correlations, they provide more reliable information on the goodness of fit between the data and a given model than simple inspection of the figures.
The results show excellent consistency between the CMB maps estimated using the different componentseparation methods. No statistically significant deviations of the CMB maps from Gaussianity are found. Indeed, the slight preference for superGaussianity of the equilateral 3point and 4point functions observed for the 2013 data is now less pronounced. That may be caused by differences between the masks used for the analysis. Interestingly, the 2point function shows clear evidence of a lack of structure for large separation angles. Such behaviour was originally noted for the WMAP firstyear data by Bennett et al. (2003), and has subsequently been discussed at length in the literature (Efstathiou 2004; Copi et al. 2007, 2015). We will return to this issue in Sect. 5.2.
4.4. Minkowski functionals
The Minkowski functionals (hereafter 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. 2004b; Curto et al. 2007; De Troia et al. 2007; Spergel et al. 2007; Curto et al. 2008; Hikage et al. 2008; Komatsu et al. 2009; Planck Collaboration XXIII 2014). They are additive for disjoint regions of the sky and invariant under rotations and translations. In the literature, the contours are traditionally 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, namely the area V_{0}(ν) = A(ν), the perimeter V_{1}(ν) = C(ν), and the genus V_{2}(ν) = G(ν), are defined respectively as (9)
where 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.
For a Gaussian random field in pixel space, the MFs can be written in terms of two functions: A_{k}, which depends only on the power spectrum, and v_{k}, which is a function only of the threshold ν (see, e.g., Vanmarcke 1983; Pogosyan et al. 2009; Gay et al. 2012; Matsubara 2010; Fantaye et al. 2015). The analytical expressions are (12)with and (15)The amplitude A_{k} depends only on the shape of the power spectrum C_{ℓ} through the rms of the field σ_{0} and its first derivative σ_{1}: where ω_{k} ≡ π^{k/ 2}/ Γ(k/ 2 + 1).
Since this factorization is still valid in the weakly nonGaussian case, we can use the normalized MFs, v_{k}, to focus on deviations from Gaussianity, with a reduced sensitivity to cosmic variance.
Apart from the characterization of the MFs using fullresolution temperature sky maps, we also consider results at different angular scales. In this paper, two different approaches are considered to study these degrees of freedom: in real space via a standard Gaussian smoothing and degradation of the maps; and, for the first time, in harmonic space by using needlets. Such a complete investigation provides an insight regarding the harmonic and spatial nature of possible nonGaussian features detected with the MFs.
Probability as a function of resolution for the unnormalized, classical Minkowski functionals.
First, we apply scaledependent analyses in real space by considering the sky maps at different resolutions. The three classical MFs – area, contour length, and genus – are evaluated over the threshold range −3 ≤ ν ≤ 3 in σ units, with a step of 0.5. This provides a total of 39 different statistics. The values of these statistics for the Planck data are all within the 95% confidence region when compared with Gaussian simulations for all of the resolutions considered. A χ^{2} value is computed for each componentseparation method by combining the 39 statistics and taking into account their correlations (see e.g., Curto et al. 2007, 2008). The corresponding covariance matrix is computed using 1000 simulations. The pvalue of this χ^{2} test is presented in Table 5 for each component separation technique and for map resolutions between N_{side} = 1024 and N_{side} = 16. We find no significant deviations from Gaussianity for any of the resolutions considered.
Then we consider the four normalized functionals described above. For every scale we used 26 thresholds ranging between −3.5 and 3.5 in σ units, except for θ = 640′ where 13 thresholds between −3.0 and 3.0 in σ units were more appropriate. Table 6 indicates that no significant deviation from Gaussianity is found.
Probability as a function of resolution determined using normalized MFs.
Fig. 3 Needlet space MFs for Planck 2015 data using the four componentseparated maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue); the grey regions, from dark to light, correspond, respectively, to 1, 2, and 3σ confidence regions estimated from the 1000 FFP8 simulations processed by the Commander method. The columns from left to right correspond to the needlet parameters j = 4,6, and 8, respectively; the jth needlet parameter has compact support over multipole ranges [2^{j−1},2^{j + 1}]. The ℓ_{c} = 2^{j} value indicates the central multipole of the corresponding needlet map. Note that to have the same range at all the needlet scales, the vertical axis has been multiplied by a factor that takes into account the steady decrease of the variance of the MFs as a function of scale. 
Third, we tested MFs on needlet components. The needlet components of the CMB field as defined by Marinucci et al. (2008) and Baldi et al. (2009) are given by: (18)Here, denotes the component at multipole ℓ of the CMB map , i.e., (19)where denotes the pointing direction, B is a fixed parameter (usually taken to be between 1 and 2) and b(.) is a smooth function such that ∑ _{j}b^{2}(ℓ/B^{j}) = 1 for all ℓ. Fantaye et al. (2015) show in a rigorous way that a general analytical expression for MFs at a given needlet scale j, which deals with an arbitrary mask and takes into account the spherical geometry of the sky, can be written as (20)where t_{0} = 2, t_{1} = 0, and t_{2} = 4π are respectively the EulerPoincaré characteristic, boundary length, and area of the full sphere. The quantities v_{k} are the normalized MFs given in Eq. (13), while the needlet scale amplitudes have a similar form as A_{k} but with the variances of the map and its first derivative given by Implementing the MFs in needlet space has several advantages: the needlet filter is localized in pixel space, hence the needlet component maps are minimally affected by masked regions, especially at highfrequency j; and the doublelocalization properties of needlets (in real and harmonic space) allow a much more precise, scalebyscale, interpetation of any possible anomalies. While the behaviour of standard allscale MFs is contaminated by the large cosmic variance of the low multipoles, this is no longer the case for MFs evaluated at the highest needlet scales; in such circumstances, the variance of normalized components may be shown to decrease steadily, entailing a much greater detection power in the presence of anomalies. Finally, and most importantly, the needlet MFs are more sensitive to the shape of the power spectrum than the corresponding allscale MFs.
Fig. 4 Histograms of χ^{2} for the Planck 2015 Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) foregroundcleaned maps analysed with the common mask. The χ^{2} is obtained by combining the three MFs in needlet space with an appropiate covariance matrix. The histograms are for the FFP8 simulations, while the vertical lines are for the data. The figures from left to right are for the needlet scales j = 4,6, and 8, with the central multipoles ℓ_{c} = 2^{j} shown in each panel. 
The needlet parameters we use in this analysis are B = 2, j = 3,4,5,6,7,8. Since the masks in pixel space are mapresolution dependent, we also use different masks for each needlet scale. These new masks are constructed by multiplying the highresolution common mask with the upgraded version of the appropriate lowresolution common mask. For needlet scales j = 2 and j = 3, we use the common mask defined at N_{side} = 16, and upgraded to N_{side} = 2048. Similarly, for the higher needlet scales, j = 2^{n}, where n = 4,5,6,7,8, we use upgraded versions of the common masks defined at N_{side} = 2^{n}.
The results concerning needlet MFs from the Commander, NILC, SEVEM, and SMICA foregroundcleaned temperature maps for needlet scales B = 2, j = 4,6,8 are shown in Fig. 3. All cases are computed using 26 thresholds ranging between −3.5 and 3.5 in σ units.
Probability as a function of needlet scale.
The figure shows the fractional difference between the Planck data and the FFP8 simulations in area (top panels), boundary length (middle panels), and genus (bottom panels) for different needlet scales. The jth needlet scale has compact support over the multipole ranges [2^{j−1},2^{j + 1}]. All the scales we considered are consistent with the Gaussian FFP8 simulations. This can be seen in Fig. 4, where we compare the data and simulation χ^{2} values, which are computed by combining the three MFs with an appropiate covariance matrix. The vertical lines in these figures represent the data, while the histogram shows the results for the 1000 FFP8 simulations. We also show in Table 7 the pvalues for the four componentseparation methods, as well as all needlet scales we considered. Despite the relatively small pvalues for some scales, the Planck temperature maps show no significant deviation from the Gaussian simulations up to ℓ_{max} = 512, which corresponds to the maximum multipole of our highestfrequency needlet map.
4.5. Multiscale analyses
Multiscale data analysis is a powerful approach for probing the fundamental hypotheses of the isotropy and Gaussianity of the CMB. The exploration of different scales (in an almost independent manner) not only helps to test the specific predictions of a given scenario for the origin and evolution of the fluctuations, but also is an important check on the impact of systematic errors or other contaminants on the cosmological signal.
There are several ways of performing a multiscale analysis, the simplest being to smooth/degrade the CMB map to different resolutions. However, in this section, we will focus on image processing techniques related to the application of wavelets and more general bandpass filtering kernels to the original CMB fluctuations. The advantage of waveletlike analyses over scale degradation is clear: they allow the exploration of characteristics of the data that are related to specific angular scales. Wavelets have already been extensively used in the study of the Gaussianity and isotropy of the CMB (e.g., McEwen et al. 2007; Vielva 2007). Indeed, a waveletbased (needlet) analysis of the Planck 2015 data has already been presented in Sect. 4.4.
We recall that in the 2013 analysis, some of the applied estimators deviated from the null hypothesis. In particular, it was determined that the cold area of the spherical Mexican hat wavelet (SMHW, MartinezGonzález et al. 2002) coefficients at scales of around 5° yielded a pvalue of 0.3%. In addition, we also found an excess in the kurtosis of the wavelet coefficients on the same scales. Previous analyses (for a review, see Vielva 2010) have suggested that the “Cold Spot” (see Sect. 5.7) was the major contributor to these statistical outliers.
In what follows, we will consider the application of the SMHW, together with matched filters for a 2DGaussian profile (GAUSS), and for generalized spherical SavitzkyGolay kernels (SSG, Savitzky & Golay 1964, see Appendix A).
The application of a filter ψ(R,p) to a signal on the sky S(p) can be written as (23)where p represents a given position/pixel, R parameterizes a characteristic scale for the filter (e.g., a wavelet scale), is the window function associated with the filter ψ(R,p), ℓ_{max} is the maximum multipole allowed by the corresponding HEALPix pixelization, and Y_{ℓm}(p) is the spherical harmonic basis. Here, s_{ℓm}, the spherical harmonic coefficients of the analysed map, are given by (24)where dΩ = dθsinθdφ and the asterisk denotes complex conjugation. Note that the filtered map (or the wavelet coefficient map, if ψ(R,p) is a continous wavelet) conserves the statistical properties of the original map, since the convolution is a linear operation. In particular, if S(p) is a Gaussian and statistically isotropic random signal, ω_{S}(R,p) is also Gaussian and statistically isotropic.
In the present work, the signal S(p) corresponds to a temperature map T(p). Several statistics can then be computed from the derived filtered map as a function of the filter scale, in particular, the first moments (the dispersion σ_{R}, the skewness S_{R}, and the kurtosis K_{R}), the total area above/below a given threshold, and the peak distribution. These statistics are compared to the corresponding results determined from the FFP8 simulations to establish the degree of compatibility with the null hypothesis.
4.5.1. Firstorder moments of the multiscale maps
For the three filters considered (SMHW, GAUSS, and SSG84^{3}) the variance, skewness, and kurtosis are computed at 18 scales, R(arcmin) = {2, 4, 7, 14, 25, 50, 75, 100, 150, 200, 250, 300, 400, 500, 600, 750, 900, 1050}. These scales are chosen to be consistent with previous analyses. They cover a wide angular range, and are selected so that the intervals between them increase with scale. Notice that, for a given scale, the three filters do not cover exactly the same multipole range, since that depends on the specific filter definition. This can be seen in Fig. 5: the SMHW is the narrowest filter, followed by SSG84, then GAUSS. The three filters have an equivalent effective ℓ_{max}, but differ in the effective ℓ_{min}. Overall, the differences between the filters become smaller with increasing effective scale. In this paper, we refer to both the scale, R, and FWHM as parameters defining the size of the filters. For the SMHW, these are related by , whereas for the GAUSS and SSG84 filters, the scale is defined to be half the FWHM. The latter definition is appropriate for filters that include prewhitening since it is simple yet matches the ℓspace bandwidth reasonably well.
Fig. 5 Comparison of the window functions (normalized to have equal area) for the SMHW (blue), GAUSS (yellow), and SSG84 (magenta) filters. The scales shown are 25′ (top) and 250′ (bottom). 
Following the procedure explained in PCIS13, after convolution with a given filter, the common mask is extended to omit pixels from the analysis that could be contaminated by the mask. These pixels introduce an extra correlation between the data and the simulations, degrading the statistical power of the comparison with the null hypothesis (see, e.g., Vielva et al. 2004). For a given scale R, the exclusion mask is defined by extending an auxiliary mask by a distance 2R from its border, where the auxiliary mask is that part of the common mask related to residual diffuse Galacic emission (i.e., the auxiliary mask does not mask point sources).
The following figures represent the uppertail probability (UTP) for a given statistic, i.e., the fraction of simulations that yield a value equal to or greater than that obtained for the data. In fact, as explained in PCIS13, if a given UTP is larger than 0.5, a new quantity is defined as mUTP = 1−UTP. Therefore, mUTP is constrained to lie between 1 /N and 0.5, where N is the number of simulations used for each statistic.
Figure 6 presents the comparison of the CMB temperature maps with the corresponding simulations for the SMHW, GAUSS, and SSG84 filters. The full mission Planck data confirm the results already obtained with the 2013 release for temperature. In particular, for the SMHW, we find (i) an excess of kurtosis (≈ 0.8%) at scales of around 300′; (ii) that the dispersion of the wavelet coefficients at these scales and at around 700′ is relatively low (≈ 1%); and (iii) that the dispersion of the wavelet coefficients at scales below 5′ is significantly high (≲ 0.1%).
The excess of kurtosis has been previously associated with the Cold Spot (e.g., Vielva et al. 2004), and the low value of the standard deviation of the coefficients on large scales could be related to the low variance discussed in Sect. 5.1. Regarding the large dispersion of the coefficients on the smallest scales, this can be understood either by the presence of residual foreground contributions (extragalactic point sources) or by incomplete characterization of the true instrumental noise properties by the FFP8 simulations. We explore these possibilities with two additional tests undertaken with the SMHW.
Figure 7 shows the significance of the statistics derived from the SEVEM100, SEVEM143, and SEVEM217 maps. The three cleaned maps yield very consistent values of the mUTP for the standard deviation, skewness, and kurtosis of the wavelet coefficients, with only small differences seen at small scales. This frequencyindependence of the results argues against the foreground residuals hypothesis. Figure 8 presents the same statistics as applied to an estimator of the noise properties of the CMB maps. This is derived from the halfdifference of the halfring data sets, which provides the best estimate of the noise properties of the full mission data set. However, since there is still a known mismatch in noise properties, any conclusions will be more qualitative than quantitative. Nevertheless, the noise study reveals that, at the smallest scales, there are some discrepancies with the FFP8 simulations, and in particular the estimated dispersion of the SMHW noise coefficients is higher than predicted.
Fig. 6 Modified upper tail probabilities (mUTP) obtained from the analyses of the filter coefficients as a function of the filter scale R for the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) sky maps. From left to right, the panels correspond to standard deviation, skewness, and kurtosis results, when determined using the SMHW (top), GAUSS (middle), and SSG84 (bottom) filters. The squares represent UTP values above 0.5, whereas circles represent UTP values below 0.5. 
Fig. 7 Modified upper tail probabilities (mUTP) obtained from the analyses of the SMHW coefficients as a function of the wavelet scale R for the SEVEM100 (blue), SEVEM143 (yellow), SEVEM217 (magenta), and SEVEM (green) maps. From left to right, the panels correspond to the standard deviation, skewness, and kurtosis. 
Fig. 8 Modified upper tail probabilities (mUTP) obtained from the analyses of the SMHW coefficients as a function of the wavelet scale R for the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) halfring halfdifference noise estimates. From left to right, the panels correspond to the standard deviation, skewness, and kurtosis. 
4.5.2. The area above/below a threshold
In the context of the study of the Cold Spot, the area above/below a given threshold, as a function of the SMHW wavelet scale, has been demonstrated to provide a useful and robust statistic (e.g., Cruz et al. 2005), since it is rather independent of any masking required. Our previous analysis (PCIS13) confirmed that the CMB temperature fluctuations exhibit an anomalously large cold area on scales of around 10°, which can be mostly associated with the Cold Spot. Here, we extend the analysis by including results derived using the GAUSS and SSG84 filters.
At a given scale R and threshold ν, the cold () and hot () areas of a filtered map are defined as where the operator # represents the number of pixels p in which the condition defined between the braces is satisfied.
Table 8 summarizes the results for the hot and cold areas determined for the four CMB temperature maps analysed with the common mask (and its associated exclusion masks). The results are similar to those obtained in 2013, with some small differences on those scales related to the Cold Spot (between 200′ and 400′). Specifically, the cold area is slightly less significant for smaller values of R, whereas the anomalous behaviour remains for larger filter scales. The three filters are in reasonable agreement, but, as expected from Fig. 6, the SMHW yields higher significance levels than the SSG84 and GAUSS filters. However, it is worth recalling that, for a given scale, the three filters are not probing exactly the same multipole range and therefore some differences should be expected.
In Fig. 9 we plot the areas for thresholds ν > 3.0, where the threshold is defined in units of σ_{R}, as determined from the SEVEM temperature map. The results for Commander, NILC, and SMICA are in good agreement with these. The panels refer to SMHW scales of R = 200′, 250′, 300′, and 400′. The most extreme value (in terms of σ_{R}) for each area is indicated.
Fig. 9 Cold and hot areas for thresholds ν> 3.0 as determined from the SEVEM temperature map. From top to bottom, the maps are for SMHW scales of R = 200′, R = 250′, R = 300′, and R = 400′. 
The coldest area corresponds to the Cold Spot with the minimum value of the wavelet coefficient at the position (209°,−57°) in Galactic coordinates. The hottest area has already been identified in the WMAP data (e.g., Vielva et al. 2007). The results are insensitive to the choice of CMB temperature map that is adopted. It is clear that the southern Galactic hemisphere yields more anomalous signatures than the northern one. These results confirm the importance of the Cold Spot as the most extreme feature in the analysed sky. More insights about its nature are provided in Sect. 5.7.
Modified upper tail probability (mUTP ) for the cold (top) and hot (bottom) areas.
4.5.3. Peak statistics
The statistical properties of local extrema (both minima and maxima, which we refer to collectively as “peaks”) provide an alternative approach to search for evidence of nonGaussianity in the data. Such peaks, defined as pixels whose amplitudes are higher or lower than the corresponding values for all of their nearest neighbours, trace topological properties of the data. Peak locations and amplitudes, and various derived quantities, such as their correlation functions, have previously been used to characterize the WMAP sky maps by Larson & Wandelt (2004, 2005) and Hou et al. (2009).
The statistical properties of peaks for a statistically isotropic Gaussian random field were derived by Bond & Efstathiou (1987). The integrated number density of peaks, n_{pk} (composed of maxima and minima with corresponding densities n_{max} and n_{min}), with amplitudes x above a certain threshold ν = x/σ is given by (27)where σ is the rms fluctuation amplitude measured on the sky, and γ is the spectral shape parameter of the underlying field. Uncharacteristically cold and hot spots are then manifested as extreme outliers in the peak values, and can constitute evidence for nonGaussianity or deviation from isotropy.
Here, we consider the peak statistics of the Planck componentseparated temperature maps at N_{side} = 2048. The maps are prewhitened as described in Appendix A. This step allows the construction of an estimator that is nearly optimal with respect to the fiducial CMB properties. After application of the common mask, weighted convolutions of the data are performed with either SSG or GAUSS kernels of variable scale. In order to avoid potential contamination by boundary effects, the mask is extended by rejecting pixels with an effective convolution weight that differs from unity by more than 12%. Peaks are extracted from the filtered map (removing any that are adjacent to masked pixels), their positions and values are recorded for further analysis, and their cumulative density function (CDF) is constructed by sorting peak values. Table 9 presents peak counts for the componentseparated sky maps for several different kernels and representative filtering scales, together with the number of peaks that are common to all maps. There is excellent agreement between the various CMB estimates. All statistical inference is then performed by comparison of the peak distributions derived from the data with equivalently processed simulations. As an internal consistency check, the properties of the FFP8 simulations are found to be in agreement with the predictions of Eq. (27).
Peak counts in maps filtered to different scales.
Figure 10 presents the distributions of peaks for the SMICA CMB map filtered with two representative kernels on scales of 40′ and 800′ FWHM. The lower panels show the empirical peak CDFs as a function of peak value x, defined for a set of n peaks at values { X_{i} } as (28)For plotting purposes alone, the horizontal axis is scaled in units of σ defined by Eq. (27) and derived from the underlying median CDF, , of the simulations. The upper panels show the difference between the observed and median simulated CDF values, , with the grey bands representing the 68.3%, 95.4%, and 99.7% regions of the simulated CDF distributions. The maximal value of this difference defines a KolmogorovSmirnov (KS) deviation estimator: (29)This forms the basis of a standard KS test of consistency between the two distributions. Although the KS deviation has a known limiting distribution, we also derive its CDF directly from the simulations.
Fig. 10 Cumulative density function of the peak distribution for the SMICA CMB temperature map. The top row shows the peak CDF for maps filtered with a GAUSS kernel of 40′ FWHM. The bottom row shows the corresponding peak CDF for an SSG84 kernel of 800′ FWHM. The spectral shape parameter γ (see Eq. (27)) is the bestfit value for the simulated ensemble, as indicated by the cyan circle in Fig. 11. Similar results are obtained for the other componentseparation methods. 
Fig. 11 Distribution of bestfit Gaussian peak CDF spectral shape parameters, σ and γ (as defined in Eq. (27)), recovered from 1000 simulations, as indicated by the black dots and the smoothed density map and compared to those derived for the observed sky (shown by the red star). The blue contours enclose 68% and 95% of the parameter distribution, and the cyan circle represents the bestfit parameters for the median peak CDF determined from simulations. The upper panel shows the peak CDF parameters for the SMICA map filtered with a GAUSS kernel of 40′ FWHM. The lower panel shows the corresponding peak CDF for an SSG84 kernel of 800′ FWHM. Similar results are obtained for the other componentseparation methods. 
Fig. 12 Cumulative density function of the mean amplitude of all extrema, maxima (red) and minima (blue), derived from simulations, compared to the equivalent values observed for the SMICA CMB temperature map. The upper panel shows the peak mean amplitudes for maps filtered with a GAUSS kernel of 40′ FWHM. The lower panel shows the corresponding peak CDF for an SSG84 kernel of 800′ FWHM. Similar results are obtained for the other component separation methods. Since the filter kernel normalization is free, and the prewhitened map to which the filter is applied is dimensionless, the plots are essentially in arbitrary units. 
The temperature peak distributions in Fig. 10 are consistent with Gaussian peak statistics, apart from a single anomalously cold peak on scales around 800′ FWHM. This corresponds to the previously reported Cold Spot. Although this exercise confirms that the Cold Spot is a rare cold feature, as already noted by Cruz et al. (2005) and confirmed in this paper, the most peculiar characteristic of the Cold Spot is not its coldness, but rather its size. A more detailed analysis of its nature is presented in Sect. 5.7.
The probability that the observed sky exceeds the value of the KS deviation for the adopted fiducial cosmology can be determined by counting the number of simulations with . The pvalues for the KS test comparing the CDF of the observed sky with the median peak CDF derived from simulations for several different kernels and representative scales are summarized in Table 10. The similarly derived pvalues for the total peak counts are summarized in Table 11. Most of the results indicate that the two distributions are highly consistent, with the exception of results for the SSG84 filter on scales of about 500′ FWHM. This deviation appears to be related to a hemispherical asymmetry in the peak CDFs, and will be discussed further in Sect. 5.6.
Modified upper tail probability (mUTP) for the KS test, comparing the data with the median peak CDF derived from simulations.
Modified upper tail probability (mUTP) for the total peak count, comparing the data with the peak count CDF derived from simulations.
One can also test whether the observed values of the parameters, σ and γ as defined in Eq. (27), are consistent with the simulation ensemble, under the assumption that the peak distributions in the Planck data are described by a Gaussian peak CDF. Figure 11 demonstrates the consistency of the bestfit values of these parameters, corresponding to the peak distributions in Fig. 10, with equivalent values derived from the simulations.
Inspired by the analysis of the WMAP firstyear data in Larson & Wandelt (2004) which found fewer extreme peaks than expected, we additionally evaluate whether the distributions of maxima and minima are separately consistent with simulations. The mean of all maxima, and the negative of the mean of all minima, are calculated for the filtered map, and the observed values are compared to the simulated distributions in Fig. 12. The observed minima/maxima means are found to be in good agreement with the fiducial values.
Fig. 13 Fraction of the Gaussian random field realizations in which the coldest peak is as cold as or colder than that observed, as a function of SMHW filter scale for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). 
The probability that the coldest peak seen on the sky is consistent with the adopted fiducial cosmology is evaluated as a function of both filter shape and size by counting the number of simulations with . The results obtained for the SMHW filter are summarized in Fig. 13. Consistent behaviour is seen when the GAUSS and SSG84 filters are applied. The error bars represent the sampling uncertainty due to the finite number of realizations, and are determined using a bootstrap method. As the filters overlap substantially, different points are highly correlated. The Planck CMB maps are consistent with the expectations of a statistically isotropic Gaussian model. The most significant deviation, found at an effective filter bandwidth given by ℓ = 20, is attributable to a single region on the sky – the Cold Spot.
4.5.4. Peak locations as a function of scale
Fig. 14 Peak positions and CDF rank summarized for all filtering scales. The three skyview panels in the top row show a Lambert projection of the north pole, the usual full sky Mollweide projection, and a Lambert projection of the south pole. The lower panel shows the peak heights (in percentile of the peak distribution on the horizontal axis) as a function of filter scale (on the vertical axis, in logarithmic scale), truncated to larger scales for clarity. Circles represent peaks (nodes of the graph) coloured according to their percentile level, and scaled according to kernel size. Black lines represent edges connecting peaks at different scales (according to a minimal distance measure). The components connected to the coldest and hottest peaks at any scale are highlighted by thicker edges, and are navy blue and dark red in colour. Note that there are thick lines that do not touch the 0 and 1 percentiles in the plot view. Those edges are connected to extreme percentile values, but at scales smaller than those shown in the plot. The Cold Spot is represented by the connected nodes that have the smallest percentiles except for the coarsest scale in the plot view. 
The application of a filter kernel of variable size to a map extends it into what can be considered a “multiscale space”, such that features on different scales are represented by a oneparameter family of smoothed maps. This technique is often used for feature detection and mathematical morphology analysis. Here, we introduce a morphological description of temperature maps based on the peak connectedness graph in multiscale space, and apply this technique to a statistical analysis of the Planck CMB data. Like most morphological analyses, it is equally applicable to intrinsically nonGaussian maps, but here we focus on the Gaussian random field statistics and attempt to understand what features of the CMB temperature map are responsible for the Cold Spot.
To construct a multiscale representation, we trace the location of the peaks in the smoothed, whitened CMB map as the smoothing scale is varied. As the smoothing scale increases, peaks merge and the total peak count decreases. Linking closest peak neighbours in positionscale space, from the finest to the coarsest resolution, produces an acyclic graph that encapsulates the peak “merger tree” history as the scale is varied. A summary of all the peak positions and CDF ranks for the SSG84 filter kernel on scales ranging from 120′ to 1200′ FWHM is shown in Fig. 14. The peaks are represented by discs of varying size (reflecting the filter scale) and colour (reflecting the peak temperature rank), with peaks at all scales projected onto a single map. The lower panel shows the peak linkage graph on the coarser scales; for the statistical analysis 81 filter scales are used, logspaced from 120′ to 1200′. Peaks of the same type (i.e., maxima to maxima and minima to minima) are linked to the closest peak on the coarser scale according to a distance measure, ds^{2} + df^{2}, where ds^{2} is the metric on the unit sphere, and df^{2} is the difference of peak temperature ranks (but only if that distance is within a predetermined fraction of the filter scale FWHM).
Fig. 15 Distribution of node degrees in the multiscale peak linkage graph determined for the SMICA map (cyan), compared with the median (red line), first to third quartile (blue box), and 95% (whiskers) derived from 1000 FFP8 simulations. 
The resulting peak linkage graph is then analysed for connectedness. The simplest quantifiable measure is the nodedegree distribution, shown in Fig. 15 for SMICA. The nodedegree distribution is highly peaked at 2; this population corresponds to a single peak being traced across multiple scales. Prewhitening effectively decorrelates the Gaussian map across different scales, so that the resulting node distribution shows a sizeable population of degree 0 and 1 nodes. When compared to the linkage graphs derived from the simulation ensemble, the nodedegree distribution of the peak linkage graph derived from Planck CMB data is consistent, with a slight excess in node counts of degrees 5 and 6.
5. Anomalies in the microwave sky
The previous section established the lack of evidence for significant nonGaussianity in the Planck data. Here we consider several important anomalies that were originally detected in the WMAP sky maps, and later confirmed in the analyses described in PCIS13. Many of these are connected to evidence for a violation of isotropy, or to a preferred direction, in the CMB. Tests that involve dipolar power asymmetry, either directly or via measures of directionality, are collected together in Sect. 6. In this section we consider only those anomalies not directly related to dipolar power asymmetry.
The microwave sky is intrinsically statistically anisotropic due to our motion with respect to the CMB rest frame. The resulting Doppler boosting effect, introduced in Sect. 1, was detected in the 2013 Planck data (Planck Collaboration XXVII 2014). For completeness, Appendix B repeats the analysis with the Planck full mission data set, though based only on the full velocity estimator (β), which is the sum of the modulation and the aberration contributions. However, since the effects of Doppler boosting are now included in the simulations used for that analysis, this constitutes a consistency check for this release. More importantly, since both the data and simulations now include the effect, it is not necessary to consider deboosted data in many of the studies reported here, unlike in PCIS13 (although one exception in Sect. 6.4 makes use of unboosted simulations to search for the frequencydependent signature of the effect in the SEVEM100, SEVEM143, and SEVEM217 sky maps). However, we note that some care must be taken due to the absence of the aberration contribution in the simulations. Indeed, this leads to the slightly, but not alarmingly, low PTE for β_{} in Appendix B. However, we not expect any impact on the results presented in this section.
Before presenting our results, we return to the issue of a posteriori correction, which particle physicists refer to as correcting for the “lookelsewhere effect” (LEE). Since there are many tests that can be performed on the data to look for a violation of statistical isotropy, we expect some to indicate detections at, for example, roughly 3σ levels, since even a statistically isotropic CMB sky is a realization of an underlying statistical process corresponding to many independent random variables. However, in the absence of an existing theoretical framework (i.e., a physical model) to predict such anomalies, it is difficult to interpret their significance. It is then necessary, and equally challenging, to address the question of how often such detections would be found for statistically isotropic Gaussian skies. Unfortunately, it is not always clear how to answer this question.
There will always be a degree of subjectivity when deciding exactly how to assess the significance of these types of features in the data. As an example, one could argue that the largescale dipole modulation signal we see is coming specifically from superHubble modes, in which case performing an LEE correction for dipole modulation that could have been seen on small scales (ℓ ≳ 100) would not make sense. Models for such a superHubble modulation exist and an example was examined in Planck Collaboration XX (2016), the conclusion being that the model could only explain part of the dipole modulation and that the allowed part was perfectly consistent with cosmic variance.
In this paper, we adopt a pragmatic approach. When there is a clear mechanism for doing so, we attempt to correct for the “multiplicity of tests”, or the possible ways in which an anomalous signal might have been detected but was not, as a consequence of any a posteriori (datadriven) choices made in searching for it. In such cases, a strong dependence of the significance on the correction would indicate that we should be cautious about the uncorrected result. When such an obvious correction is not possible, we clearly describe the methodology applied to the data and its limitations. With this approach, we also recognize that any statistical assessment is partially subjective, including those that purport to correct for the LEE.
Although many of the observed effects described in this and the next section may elude theoretical prediction today, we continue to highlight them since there is a real possibility that the significance of one or more might increase at a later date, perhaps when polarization data are included in the analysis, and lead to new insights into early Universe physics. Alternatively, such observations may directly motivate the construction of models that can make predictions for features that can be sought in new data sets. This is particularly the case for anomalies on the largest angular scales, which may have a specific connection to inflation.
5.1. Variance, skewness, kurtosis
Previous analyses of the WMAP data (Monteserín et al. 2008; Cruz et al. 2011; Gruppuso et al. 2013) have reported that the variance of the CMB sky is lower than that of simulations based on the ΛCDM model. PCIS13 confirmed this, and proposed a possible explanation of the apparent incompatibility of the observed variance with a fiducial cosmological model that has been determined from the same data set. Specifically, whilst the mapbased variance is dominated by contributions from large angular scales on the sky, the cosmological parameter fits are relatively insensitive to these loworder ℓmodes, and are instead largely dominated by scales corresponding to ℓ> 50. Therefore the variance of the map appears to be anomalous, since there is a dearth of largeangularscale power compared to the model.
In Sect. 4.1, we again confirmed the presence of low variance in the data. Here, we extend the analysis to investigate which angular scales are responsible for the low variance by applying the unit variance estimator to lower resolution componentseparated maps, specifically those from N_{side} = 1024 to N_{side} = 16, with the corresponding common masks, and then comparing the results with those determined from 1000 MC simulations. The results are shown in Fig. 16.
Fig. 16 Lower tail probability of the variance (top panel), skewness (centre panel), and kurtosis (bottom panel) obtained at different resolutions from the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) sky maps. 
All of the componentseparation methods that we consider yield very consistent results which indicate an increasingly anomalous low variance at lower resolutions, with the lowertail probability reaching a minimum value of 0.5% at N_{side} = 16. We then consider the impact of a possible lookelsewhere effect by evaluating the minimum lower tail probability of each simulation irrespective of the N_{side} resolution at which it occurs. By comparing the distribution of these values with that of the data, we infer that the probability is slightly weakened to a value of about 1%. These results are compatible with a lack of power on large angular scales. Since the variance estimator is heavily weighted towards lowℓ modes, this has an increasing impact on the observed variance when going from high to low resolution sky maps. Conversely, the skewness and kurtosis are consistent with the simulations, although there is some indication of a weak scaledependence (albeit at low significance).
We also investigate the stability of the results at N_{side} = 16 with respect to the possible presence of residual foregrounds by considering two additional masks obtained by extending the edge of the N_{side} = 16 common mask by 5° and 9°, reducing the usable sky fraction from 58% to 48% and 40%, respectively. We then reapply the unit variance estimator to the low resolution componentseparated maps with these masks and determine the variance, skewness, and kurtosis values (see Table 12).
Lowertail probability for the variance, skewness, and kurtosis of the low resolution N_{side} = 16 componentseparated maps obtained with the common mask and two extended versions thereof.
The results from 48% of the sky reveal that only 1 simulation in 1000 is found to be more anomalous (i.e., exhibit lower variance) than the observed map. In addition, both the skewness and kurtosis become more compatible with the ΛCDM model. With the more aggressive mask, the lowertail probability slightly increases again. However, given the limited number of pixels involved in the analysis, this shift may be related to the effects of sample variance.
Overall, our results may be explained by the presence of a lowvariance anomaly in the primordial CMB signal – the stability of the lowvariance significance argues against foreground contamination being responsible for the lack of observed power. This is reinforced by the decrease in variance when regions close to the common mask borders, where foreground residuals are most likely to be observed, are omitted from the analysis.
5.2. Npoint correlation function anomalies
5.2.1. Lack of largeangle correlations
We first reassess the lack of correlation seen in the 2point angular correlation function at large angular separations as reported in Sect. 4.3, and previously noted for both WMAP and the 2013 Planck temperature maps (Bennett et al. 2003; Copi et al. 2015). We attempt to quantify this lack of structure using the statistic proposed by Spergel et al. (2003): (30)where Ĉ_{2}(θ) is our estimate of the 2point correlation function. Generally, the upper limit on the integral has been taken to correspond to a separation angle of 60°, possibly (as noted by Copi et al. 2009) motivated by the COBEDMR 4year results (Hinshaw et al. 1996). Inspection of the top panel of Fig. 2 suggests that the Planck 2point function lies close to zero between 80° and 170°, but for consistency with previous work we compute the statistic S_{1 / 2}, for . The results are presented in Table 13. We find that the data indeed show a lack of correlations on large angular scales, with a significance consistent with that found by Copi et al. (2015) (although note that the sense of the pvalues differs between the papers).
Probabilities of obtaining values for the S_{1 / 2} and statistics for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2015 temperature CMB maps with resolution parameter N_{side} = 64, estimated using the Commander, NILC, SEVEM, and SMICA maps.
Possible criticisms of the S_{1 / 2} statistic include that it has been designed a posteriori to test for a lack of largeangle correlations, and that it does not account for the high degree of correlation between bins at different angular scales. We can address these concerns, at least in part, by considering a modified version of the commonly used and well understood χ^{2} statistic used in previous studies. In order to test the same hypothesis as the S_{1 / 2} statistic – that there are no correlations on scales larger than some angular cutoff – we do not subtract an averaged 2point function when computing the χ^{2}, i.e., we use a statistic defined as (31)where i_{min}, i_{max} denote the index of the bins corresponding to the minimum and maximum value of the separation angles θ_{min} and θ_{max}, respectively. In this analysis, we adopt θ_{min} = 60° and θ_{max} = 180°. M_{ij} is the covariance matrix given by Eq. (8), estimated using MC simulations corresponding to the fiducial ΛCDM model. The results are shown in Table 13. The significance level of the anomaly is slightly smaller for the statistic than that derived with S_{1 / 2}. We note that this statistic is closely related to the A(x) measure proposed by Hajian (2007).
A further potential criticism of the S_{1 / 2} statistic relates to the a posteriori choice of 60° for the separation angle that delineates the interesting region of behaviour of the correlation function. We therefore consider the generalized statistic S(x) and compute its value for all values of x, both for the data and for the simulations. Then, for each value of x, we determine the number of simulations with a higher value of S(x), and hence infer the most significant value of the statistic and the separation angle that it corresponds to. However, since such an analysis is sensitive to the LEE, we define a global statistic to evaluate the true significance of the result. Specifically, we repeat the procedure for each simulation, and search for the largest probability irrespective of the value of x at which it occurs. The fraction of these probabilities higher than the maximum probability found for the data defines a global pvalue. As seen in Table 13, this corresponds to values of order 98% for all of the CMB estimates.
The previous analyses essentially test how consistent the observed 2point correlation function data is with a lack of correlations on large angular scales, in particular for separation angles θ> 60°. A conventional χ^{2} statistic allows us to test the consistency of this quantity with the predictions of the ΛCDM model. In this case, the statistic is defined as in Eq. (7), except that we constrain the computations to those bins that correspond to the intervals defined by θ< 60° and θ> 60°. The results of these studies are shown in Table 14.
Probabilities of obtaining values for the χ^{2} statistic for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2015 temperature CMB maps with resolution parameter N_{side} = 64, estimated using the Commander, NILC, SEVEM, and SMICA maps.
The analysis for θ< 60° indicates that the observed 2point function is a good match to the mean 2point function predicted by the ΛCDM model. Moreover, for θ> 60° the results suggest that the problem is that the fit of the data to the model is too good, and this is even more pronounced for an analysis in the full separation angle range.
Overall, the tests indicate an unusually good fit of the observed 2point function both to zero and to the predictions of the ΛCDM model for angles above 60°. This problem may be related to the fact that the theoretical variance for the bestfit model is larger than the observed value at large scales, so that the simulations based on this model that have been used in all of the statistical tests may overestimate the variance of the 2point function.
5.2.2. Hemispherical asymmetry
We now turn to a reassessment of the asymmetry between the realspace Npoint correlation functions computed on hemispheres reported previously for the WMAP (Eriksen et al. 2005) and Planck 2013 temperature maps (PCIS13). We initially focus the analysis on the hemispheres determined in the ecliptic coordinate frame for which the largest asymmetry was observed. However, we also carry out the corresponding calculations in other relevant reference frames, such as those defined by the Doppler boost (DB, see Sect. 6.4, Appendix B, and Planck Collaboration XXVII 2014) and the dipole modulation (DM, see Sect. 6.2) directions. We use the same configurations of the Npoint functions as described in Sect. 4.3. However, here the functions are not averaged over the full sky and depend on a choice of specific direction, so they constitute tools for studying statistical isotropy rather than nonGaussianity (Ferreira & Magueijo 1997).
As in Sect. 4.3, we analyse the CMB estimates at a resolution of N_{side} = 64 and quantify their agreement with the fiducial cosmological model using a χ^{2} statistic. The results determined from the Planck 2015 temperature data for the ecliptic hemispheres are shown in Fig. 17. If we consider that the χ^{2} statistic itself can act as a measure of fluctuation level, then asymmetry between the two measured hemispheres can be quantified by the ratio of the corresponding χ^{2} values. The probabilities of obtaining values of the χ^{2} statistic or ratio for the Planck fiducial ΛCDM model at least as large as the observed values are given in Table 15. Since we do not have any predictions concerning the behaviour of a given hemisphere, in the case of the χ^{2} ratios we provide the complementary probabilities of the 2tailed statistic.
The significance levels of the 3 and 4point functions in the northern hemisphere are nominally very high, exceeding 99.9% for the pseudocollapsed 3point function. However, proper interpretation requires that one recognize that the analysis is affected by a posteriori choices for the smoothing scale and reference frame defining the hemispheres. This typically leads to an overestimation of the significance of the results. Accounting for such effects requires the repetition of the analysis for all possible reference directions and also for data at other resolutions. Unfortunately, because of computational limitations, such an extended analysis is not possible for these higherorder statistics. Nevertheless, the observed properties of the Planck data are consistent with a clear lack of fluctuations in a direction towards the north ecliptic pole. However, the χ^{2}ratio statistic indicates a slightly smaller significance level for the asymmetry, not exceeding 99% for any of the Npoint functions.
The results for the Npoint correlation functions determined in the DB and DM reference frames for the SMICA map are shown in Fig. 18 and the probabilities are presented in Table 16. Note that the positive hemisphere for the ecliptic reference frame corresponds to the southern hemisphere in the previous table. Whilst the largest asymmetry is seen in ecliptic coordinates, a substantial asymmetry is present also for the DM direction. This can be explained by the fact that the DM direction is more closely aligned with the south ecliptic pole (with a separation of around 47°) than the DB direction is. For the DB direction we do not find any significant asymmetry. The equivalent results for Commander, NILC, and SEVEM are consistent with those shown here.
Fig. 17 Difference of the Npoint correlation functions determined from the N_{side} = 64Planck CMB 2015 temperature estimates and the corresponding means estimated from 1000 simulations. Results are shown for the 2point, pseudocollapsed 3point (upper left and right panels, respectively), equilateral 3point, and connected rhombic 4point functions (lower left and right panels, respectively). Correlation functions are shown for the analysis performed on northern (blue) and southern (red) hemispheres determined in the ecliptic coordinate frame. The solid, dashed, dotdashed, and dotted lines correspond to the Commander, NILC, SEVEM, and SMICA maps, respectively. Note that the lines lie on top of each other. The shaded dark and light grey regions indicate, for reference, the 68% and 95% confidence regions, respectively, determined from the SMICA simulations. 
Probabilities of obtaining values for the χ^{2} statistic and ratio of χ^{2} of the Npoint functions shown in Fig. 17 for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2015 CMB maps estimated on northern and southern ecliptic hemispheres.
Fig. 18 Difference of the Npoint correlation functions determined from the N_{side} = 64PlanckSMICA CMB 2015 temperature estimates and the corresponding means estimated from 1000 simulations. Results are shown for the 2point, pseudocollapsed 3point (upper left and right panels, respectively), equilateral 3point, and connected rhombic 4point functions (lower left and right panels, respectively). Correlation functions are shown for the analysis performed on negative (blue) and positive (red) hemispheres determined in the ecliptic (solid lines), Doppler boost (DB, dashed lines), and dipole modulation (DM, dotdashed lines) coordinate frames. The shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively. 
Probabilities of obtaining values for the χ^{2} statistic and ratio of χ^{2} of the Npoint functions shown in Fig. 18 for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the SMICA map on hemispheres defined by the ecliptic (first column), Doppler boost (DB, second column), and dipolar modulation (DM, third column) reference frames.
In conclusion, the correlation functions for the Planck 2015 temperature data are consistent with the results presented in PCIS13. Specifically, we observe that 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 consistent with Gaussian simulations.
5.3. Constraints on quadrupolar modulation
The most natural extension of the class of statistically anisotropic models that we have considered previously involves the quadrupolar modulation of an initially statistically isotropic CMB sky map. No detection of a corresponding quadrupolar power asymmetry is currently claimed. An initial BipoSH analysis of the WMAP 7year data (Bennett et al. 2011) found evidence of corresponding nonzero spectra, and , in ecliptic coordinates. However, Hanson et al. (2010) demonstrated that the signal could be attributed to an incomplete treatment of beam asymmetries in the data, and this was subsequently confirmed in Bennett et al. (2013). The corresponding analysis of the Planck 2013 data indicated consistency with statistical isotropy (Planck Collaboration XXIII 2014).
Here, we proceed further and consider the quadrupolar modulation of the primordial power spectrum as suggested by Ackerman et al. (2007): (32)Given such a spectrum, the CMB temperature field is expected to exhibit a correlation between a_{ℓm} and with Δ = 0,2. Therefore, the BipoSH coefficients and are sensitive to g_{2M}. In the limit of weak anisotropy, Kim & Komatsu (2013) proposed an optimal estimator for g_{2M}: (33)where a is the CMB data vector in harmonic space and C is its covariance matrix, and (34)Here, ⟨ (C^{1}a^{∗})_{ℓm}(C^{1}a)_{ℓ′m′} ⟩ _{g2M = 0} is the mean field in the absence of the quadrupolar modulation. Observationspecific issues such as incomplete sky coverage, inhomogeneous noise, and asymmetric beams will result in a nonzero mean field, which can be estimated for the Planck data using simulations. Due to the otherwise prohibitive computational cost, we adopt a diagonal approximation for the inverse of the covariance matrix: (35)where C_{ℓ} and N_{ℓ} are the signal and noise power spectra respectively. Uncertainties are computed by applying the estimator to simulations.
Constraints on the quadrupolar modulation, determined from the Commander, NILC, SEVEM, and SMICA foregroundcleaned maps.
Table 17 presents results from an analysis of the Planck data using the extended common mask, UTA76, and limiting the range of multipoles to 2 ≤ ℓ ≤ 1200. When including data at higher ℓvalues, simulations show evidence for large statistical uncertainties in the recovered g_{2M} values that are a consequence of the many holes in the mask related to point sources. Therefore, imposing this limit ℓ ≤ 1200 does not significantly affect the constraining power of the analysis. We then estimate the amplitude of the quadrupolar modulation using the relation g_{2} = ^{(}1 / 5 ∑ _{M}  g_{2M}  ^{2}^{)}^{1 / 2}. Due to the nature of the estimator, which is necessarily positive, the estimation is biased. For an unbiased assessment, we estimate the mean and standard deviation of g_{2} from simulations. We find no evidence for quadrupolar modulation of the primordial power spectrum. However, the derived limits allow us to impose tight constraints on statistically anisotropic inflationary models, such as those including vector fields during inflation. A companion paper, Planck Collaboration XX (2016), contains a more complete discussion on the theoretical implications of this constraint.
5.4. Pointparity asymmetry
The CMB anisotropy field defined on the sky, , may be divided into symmetric, , and antisymmetric, , functions with respect to the centre of the sphere, as previously described in PCIS13. These functions have even and odd parity, and thus correspond to spherical harmonics with even and odd ℓmodes, respectively. On the very large scales corresponding to the SachsWolfe plateau of the temperature power spectrum (2 ≤ ℓ ≤ 30), the Universe should be parity neutral with no particular parity preference exhibited by the CMB fluctuations. However, an odd pointparity preference has previously been observed in the WMAP data releases (Land & Magueijo 2005a,b; Kim & Naselsky 2010a,b; Gruppuso et al. 2011) and the Planck 2013 results. Here, we investigate the parity asymmetry in the 2015 temperature maps at N_{side} = 32. We consider the following estimator: (36)where D_{+}(ℓ_{max}) and D_{−}(ℓ_{max}) are given by (37)is the total number of even (+) or odd (−) multipoles included in the sum up to ℓ_{max}, and is the temperature angular power spectrum computed using a quadratic maximum likelihood (QML) estimator (Gruppuso et al. 2011). The ℓ(ℓ + 1) / (2π) factor in Eq. (37) effectively flattens the spectrum across the ℓrange of the SachsWolfe plateau (up to ℓ = 50) in a ΛCDM model.
Figure 19 presents the ratio, R^{TT}(ℓ_{max}), for the 2015 componentseparated maps, together with the distribution determined from the SMICA MC simulations which serves as a reference for the expected behaviour of the statistic in a parityneutral Universe. The distributions for the other CMB maps are very similar. The four componentseparation products are in good agreement, indicating an oddparity preference at very large scales for the multipole range considered in this test.
Fig. 19 Ratio R^{TT}(ℓ_{max}) for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) determined at N_{side} = 32. The shaded grey regions indicate the distribution of the statistic derived from the SMICA MC simulations, with the dark, lighter, and light grey bands corresponding to the 1, 2, and 3σ confidence levels. 
Figure 20 shows the lowertail probability for the data as compared to simulations as a function of ℓ_{max}. The results are in good agreement with those in PCIS13. The cleaned CMB maps yield generally consistent profiles which signify an anomalous oddparity preference in the multipole region ℓ_{max} = 20–30. The minimum in the lowertail probability occurs at ℓ = 28 corresponding to a value of 0.2% for NILC, SEVEM, and SMICA, and 0.3% for Commander^{4}.
As a first attempt to quantify any a posteriori effects in the significance levels, we consider how many MC simulations appear in the lower tail of the MC distribution with a probability equal to, or lower than, 0.2%, for at least one ℓ_{max} value over a specific range. For ℓ_{max} in the range 3−50, the total number of simulated maps with this property is less than 20 over 1000 MC maps, implying that, even considering the LEE, an oddparity preference is observed with a lowertail probability of less than 2%.
Fig. 20 Lowertail probability of the pointparity estimator for Commander (red), NILC (orange), SEVEM, (green), and SMICA (blue). 
5.5. Mirrorparity asymmetry
For the Planck 2013 data release, we studied the properties of the temperature data at a resolution of N_{side} = 16 under reflection with respect to a plane to search for mirror symmetries. Such a symmetry might be connected to nontrivial topologies (Starobinsky 1993; Stevens et al. 1993; de OliveiraCosta et al. 1996). In Planck Collaboration XXIII (2014), we reported evidence for an antisymmetry plane, with a perpendicular direction given by (l,b) = (264°,−17°), However, the probability of the results was slightly dependent on the method of foreground cleaning, with a pvalue ranging from 0.5% for CommanderRuler to 8.9% for SMICA. The same direction was also found in the WMAP 7year data (Finelli et al. 2012), and is close to that determined for the dipole modulation in the Planck 2013 data release (PCIS13), suggesting possible connections between the two directional anomalies.
We now proceed to reanalyse the status of mirror symmetries using the Planck 2015 full mission temperature data at both N_{side} = 16 and N_{side} = 32. In order to avoid possible bias introduced by the use of the Galactic mask^{5} the results are derived from the fullsky Commander, NILC, and SMICA maps described in Sect. 2. For SEVEM, a customized map is first produced by inpainting about 3% of the map along the Galactic plane using a diffusive inpainting technique. This is then smoothed to the appropriate lower resolutions for further analysis. Following Finelli et al. (2012), we consider the estimators in the pixel domain given by: (38)where the sum is over all N_{pix}HEALPix pixels, is the CMB temperature anisotropy measured at the pixel defined by the unit vector , and is the opposite direction with respect to the plane defined by , i.e., (39)Note that we expect S^{+} to be small if the points on opposite sides of the mirror are negatives of each other, and S^{−} to be small when they are the same.
Lowertail probability for the S^{±} statistics of the componentseparated maps at N_{side} = 16 and N_{side} = 32.
Fig. 21 Histograms of the S^{+} (top panel) and S^{−} (bottom panel) statistic. The vertical lines show the minimum value for the estimator computed at N_{side} = 32 for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) maps. The grey area shows the same quantity computed from 1000 simulated SMICA maps. 
We compute these quantities for each of the 3072 (12288) directions defined at resolution N_{side} = 16 (32), and allow the j and k indices to run over all of the pixels of the lowresolution fullsky maps. We perform the same analysis on 1000 FFP8 simulations and store the minimum value of S^{±} for each of these to compute probabilities. The results are summarized in Table 18 and Fig. 21.
We confirm that the full mission Planck temperature data at N_{side} = 16 exhibits the most anomalous mirror antisymmetry in the direction (l,b) = (264°,−17°), consistent with the result from the 2013 nominal mission data, with a probability which ranges from 1.6% for SEVEM to 2.9% for Commander. This is within 40° of the preferred direction identified by the dipole modulation analysis in Sect. 6.2. The corresponding results at N_{side} = 32 yield approximately the same direction, (l,b) = (264°,−16°), with a slightly increased probability, ranging from 0.8% for SEVEM to 1.9% for Commander.
We also note that the CMB pattern exhibits a mirror symmetry in the direction (l,b) = (260°,48°), consistent with that found in the WMAP 7year data (Finelli et al. 2012), and close to that identified by the solar dipole (Planck Collaboration VIII 2016). However, the significance of the symmetry pattern is less than in the antisymmetric case.
This extension of the analysis to higher resolution than in our previous work shows that the antisymmetry property does not seem to be confined to the largest angular scales, although we have not attempted to correct for any a posteriori choices made in the analysis. The detailed connection of this antisymmetry property to the lowvariance and hemispherical asymmetry observed on these scales remains an open issue.
5.6. Local peak statistics
Fig. 22 KSdeviation of the peak distribution for 70° radius discs centred on the positive and negative asymmetry directions determined from the SMICA CMB temperature map in Sect. 6.2. From top to bottom, the plots correspond to maps filtered with a GAUSS kernel of 40′ FWHM, an SSG84 filter of 500′ FWHM, and an SSG84 filter of 800′ FWHM, respectively. 
Local extrema or peaks, as introduced in Sect. 4.5.3, can be employed to search for localized anomalies on the CMB sky by examining how their statistical properties vary in patches as a function of location.
Initially, we consider a further test for asymmetry by examining the differences in the peak distribution when divided according to orientation with respect to a previously specified asymmetry direction. In particular, we select the peaks both in a disc of radius 70° centred on (l,b) = (225°,−18°) (the positive direction of the dipole defined in Sect. 6.2 for SMICA) and in the corresponding antipodal disc, then construct the empirical peak height CDFs to be compared with the fullsky median FFP8 distribution, as shown in Fig. 22. For maps filtered with a 40′ FWHM GAUSS filter the distribution of the peaks for the positivedirection disc is in general agreement with the full sky result, while that for the negativedirection is marginally different. Moreover, this pattern of behaviour is seen over a number of filtering scales, both for the KS deviation from the median fullsky simulated CDFs, and the spread of extremal values when comparing positive and negative regions. We also find that the properties of the negative disc affect the pvalue results for a full sky KS test on data filtered with an SSG84 filter of 500′ FWHM, as seen in Sect. 4.5.3.
We can then extend the analysis for the 40′ GAUSSfiltered data by considering the variation in the peak statistical properties for a set of discs, each of which is centred on a pixel defined at N_{side} = 256. The simplest statistics to consider are the peak number counts. We therefore consider discs of 30° diameter and compute the peak counts for each disc. These are then compared to the corresponding peak count CDFs determined from simulations, and the upper and lowertail probabilities are assigned by counting the number of simulations above and below the observed counts at the same location. These quantities can then be visualized in the form of N_{side} = 256 sky maps. The derived −log _{10}(UTP) maps for each componentseparation method are shown in Fig. 23. While we find that the total counts of peaks for the sky coverage defined by the common mask is consistent with simulations, significant regional variation is seen. Indeed, the pvalue for certain disc locations drops to 0.1% (i.e., the sky counts exceed anything seen in simulations). However, one needs to account for the a posteriori selection of significant regions in the determination of the true significance. It should also be noted that regional variations of the UTP are seen at similar levels when inspecting the peakcount statistics maps derived for randomly selected realizations of the simulations. Moreover, the significance of such peakcounting anomalies is degraded with larger disc diameters, and becomes insignificant for counts on the full sky. Thus, no significant anomalies can be claimed for the peakcount statistics of the Planck data.
A powerful nonparametric test of statistical isotropy is provided by the twosample KSdeviation between the full sky empirical peak height CDF F_{n}(x) (see Eq. (28)) and an empirical peak height CDF F_{n′}(x) derived from a subsample of the distribution, again defined by the peaks within discs of 30° diameter as defined above. The twosample KSdeviation (40)for a partial sky region shares samples between the two CDFs, and can be calculated extremely efficiently using rank statistics according to (41)where r and r′ denote the ranks of a value with index i in the full set of n and restricted set of n′ samples, respectively. Maps of the upper tail probability are then determined by comparison with the equivalent quantities computed from simulations; −log _{10}(UTP) maps are shown in Fig. 24. The majority of the selected locations are consistent with the fullsky distribution, thus indicating the statistical isotropy of the Planck maps. The most prominent feature in each of the local KSdeviation maps appears south of the Galactic centre and may be associated with a cold region crossing the Galactic plane. However, as with the peak counts, it cannot be interpreted as statistically anomalous.
Fig. 23 Map of −log _{10}(UTP) for peak counts in the Planck 40′ GAUSSfiltered temperature data, where each pixel encodes the probability determined for a 30° diameter disc centred on it. 
Fig. 24 Map of −log _{10}(UTP) for the twosample KSdeviation where each pixel encodes the probability determined for a 30° diameter disc centred on it, as computed from the Planck 40′ GAUSSfiltered temperature data. 
5.7. The Cold Spot
Fig. 25 Top: temperature patch centred on the Cold Spot. Bottom: peak merger tree within the Cold Spot region. The figure shows a region centred on the Cold Spot location in gnomonic projection, with all the peaks in SSG84filtered maps with FWHM ranging from 80′ to 1200′ overlaid on the same plot. The size of the coloured circles is proportional to the filtering scale. The colour corresponds to the peak value, normalized in units of σ for a given filter scale. In both panels the data are from the SMICA CMB map at full resolution. 
Since its discovery in the WMAP firstyear data (Vielva et al. 2004), the Cold Spot, centred at Galactic coordinates (l,b) = (210°,−57°) has been one of the most extensively studied largescale CMB anomalies. In the 2013 release (Planck Collaboration XXIII 2014), Planck confirmed the apparently anomalous nature of this feature in temperature, in terms of the area of the SMHW coefficients on angular scales of ≈ 10° on the sky; the 2015 release has also confirmed this feature (see Sects. 4.5.2 and 4.5.3). The CMB temperature anisotropies around the Cold Spot as observed by Planck are shown in the top panel of Fig. 25. The peak merger tree within the Cold Spot region is presented in the lower panel of the figure and provides a multiscale view of its structure (see Sect. 4.5.4 for details).
The robustness of the detection of the anomalies discussed in this paper is a nontrivial issue. For the particular case of the Cold Spot, this has been reviewed by Vielva (2010), and addressed in detail by Cruz et al. (2006), paying specific attention to the impact of a posteriori choices. In particular, the latter study focused on the original test that indicated the presence of this feature on the sky, confirming a significance between 1% and 2%. An alternative analysis of the significance based on two statistical tests with different levels of conservativeness was made by McEwen et al. (2005), providing values of 0.1% and 4.7%, respectively. The statistical significance of the Cold Spot was questioned by Zhang & Huterer (2010) who found a low significance after performing a study based on different kernels. As discussed in more detail by Vielva (2010), this result can also be interpreted as evidence that not all kernels are necessarily suitable for the detection of arbitrary nonGaussian features.
The possibility that the Cold Spot arises from instrumental systematics (Vielva et al. 2004) or foreground residuals (Liu & Zhang 2005; Cruz et al. 2006) has been largely rejected. However, several nonstandard physical mechanisms have been proposed as possible explanations. These include the gravitational effect produced by a collapsing cosmic texture (Cruz et al. 2007), the linear and nonlinear ISW effect caused by a void in the largescale structure (e.g., Tomita 2005; Inoue & Silk 2006; Rudnick et al. 2007; Tomita & Inoue 2008; Finelli et al. 2016), a cosmic bubble collision within the eternal inflation framework (Czech et al. 2010; Feeney et al. 2011; McEwen et al. 2012), and a localized version of the inhomogeneous reheating scenario within the inflationary paradigm (Bueno Sanchez 2014).
Since the other scenarios lack additional evidence, the void hypothesis would seem to be the most plausible, depending on the sizes, density contrasts, and profiles assumed in the computations, some of which are not in agreement with either observation (Cruz et al. 2008) or current Nbody studies (Cai et al. 2010; Watson et al. 2014). However, Szapudi et al. (2015) have recently detected a large void in the WISE2MASS galaxy catalogue aligned with the Cold Spot, with an estimated radius of around 200 h^{1} Mpc, an averaged density contrast of , and centred on a redshift of z ≈ 0.15. Large voids with similar characteristics are not unusual in the standard ΛCDM model (Nadathur et al. 2014). In fact, Nbody simulations predict about 20 such voids in the local Universe (z< 0.5). However, Zibin (2014) and Nadathur et al. (2014) indicate that the expected signal due to the linear and nonlinear ISW effects caused by this structure is not large enough to explain the temperature decrement associated with the Cold Spot.
The new Planck data release allows us to further explore the statistical nature of the Cold Spot. Two previous studies (Zhao 2013; Gurzadyan et al. 2014) have claimed inconsistencies of the internal properties of the Cold Spot with the Gaussian hypothesis, which we readdress here. In particular, we consider the smallscale fluctuations within a disclike region of radius ≈ 25°.
Several statistical quantities are computed from the fullresolution temperature maps within the Cold Spot region. This is divided into a central disc of diameter 1° surrounded by a set of 13 concentric annuli with central radii spaced in steps of about 2°, thus allowing us to build angular profiles for the mean, variance, skewness, and kurtosis. These are then compared to specialized CMB realizations, generated as follows. A set of Gaussian CMB skies is simulated using the FFP8 reference spectrum, and convolved with a Gaussian beam of 5′ FWHM. As for the FFP8 simulations themselves, these maps are rescaled, as discussed previously. Only those that contain a spot as extreme as the Cold Spot at a scale R = 300′ in SMHW space are retained, and these are rotated such that each simulated cold spot is relocated to the actual position of the Cold Spot (this ensures that the noise properties are identical for both data and simulations). This selection criterion corresponds to the characteristic that originally indicated the presence of the Cold Spot in the observed sky. As a final step, for each remaining CMB simulation a noise realization is added, consistent with each componentseparation method.
The results are presented in Fig. 26. Focusing on the profile of the mean value, it is apparent that the largest deviations from the simulations appear on scales around 15°, which corresponds to a hot ring structure, as seen in Fig. 25 and previously discussed in Cayón et al. (2005) and Nadathur et al. (2014). Notice that on the smallest scales the mean profile is also somewhat deviant with respect to the simulations, but this may be connected to selection bias, since we are considering CMB simulations containing a spot that is at least as cold as the Cold Spot. However, if we consider the distribution of the profiles corresponding to the coldest spots instead of the spots as extreme as the Cold Spot (removing the bias at the smallest scales) then the results do not change substantially (see below).
In order to quantify possible deviations from Gaussianity, we determine the probability of finding a χ^{2} value larger than that of the data for each statistic, as summarized in Table 19. The χ^{2} value for the data is computed using an estimate of the covariance matrix between different radial scales determined from the Cold Spot simulations (1000 for each componentseparation method), and then compared to the theoretical χ^{2} distribution with 13 degrees of freedom. The results indicate that the angular profile for the mean is poorly described by the simulations, of which less than 1% are found to have a higher χ^{2} than the data (when considering the distribution corresponding to the coldest spot this probability becomes approximately 2%). We have checked that this deviation is not obviously associated with a particular subrange of angular scales, implying that the mean profile is anomalous over the full range considered. Conversely, the radial profiles of the higherorder moments are compatible with the Gaussian simulations. The latter results are then in contradiction with a similar analysis (using discs instead of rings) by Zhao (2013) for the WMAP 9year data. However, it appears that this may be related to the criteria applied for the selection of the Gaussian simulations used to define the null hypothesis. In particular, Zhao (2013) used the coldest pixel in real space as a means to identify those simulations that should be retained, as opposed to the existence of cold spots as extreme as the Cold Spot selected in the SMHW coefficient map at R = 300′. Since it is not implicit that such a temperature extremum is necessarily associated with an extended cold region, particularly one defined in wavelet space, the simulations used by Zhao (2013) did not contain features comparable to the nature of the Cold Spot. This explains why the Cold Spot seemed to be anomalous when looking at the smallscale fluctuations.
Fig. 26 From left to right: mean, variance, skewness, and kurtosis angular profiles computed for rings at radii θ centred on the Cold Spot position for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The expected value obtained from the simulations is denoted by the black dashed line and the dark and light grey regions represent the 1σ and 2σ intervals. 
In conclusion, it appears that only the mean temperature profile of the Cold Spot should be considered anomalous when compared to CMB cold spots that are as statistically extreme. All other measures of its internal structure are consistent with expectations.
As a final remark, we note that the highpass filtering currently applied to the Planck CMB polarization maps severely limits the possibility of conducting targeted analyses to discriminate between different possible origins of the Cold Spot. For example, no polarization signal would be expected in those models producing secondary anisotropies due to a gravitational effect, whereas a specific pattern might be expected in a bubble collision scenario (Czech et al. 2010). Appropriate tests will be pursued in future work, once the largescale polarization data are available.
6. Dipole modulation and directionality
In this section, we examine isotropy violation related to dipolar asymmetry, various forms of which have been noted since the early WMAP releases (Eriksen et al. 2004a). We perform a nonexhaustive series of tests in an attempt to narrow down the nature of the asymmetry (on the assumption that it is not simply a statistical fluke). First, we will briefly describe some similarities and differences between the tests that are important for making a proper comparison of the results.
All the tests in this section share in common the fitting of a dipole. This is done either by fitting for a dipole explicitly in a map of power on the sky (Sects. 6.1 and 6.5), by employing Bayesian techniques in pixel space for a specific model (Sect. 6.2), or by measuring the coupling of ℓ to ℓ ± 1 modes in the CMB covariance matrix (Sects. 6.3, 6.4, and 6.6). The differences arise from how the fitted dipoles are combined, which determines the specific form of asymmetry that the test is sensitive to.
The tests can be divided into two categories, amplitudebased and directionbased. Sections 6.1 to 6.4 are all sensitive to the amplitude of a dipole modulation. Specifically, Sect. 6.1 looks for dipole modulation in the pixeltopixel variance of the data, while Sects. 6.2−6.4 all search for dipole modulation of the angular power spectrum. The distinction between these two approaches is mainly one of ℓ weighting.
Sections 6.5 and 6.6 both examine aspects of directionality in the data, where the directions are extracted from dipole fits but combined in different ways. Section 6.5 fits for dipoles in band power (with similar results for variance) and only uses the direction information, while Sect. 6.6 weights each dipole equally across all scales and uses the amplitude information as well.
The differences between the approaches of these sections should be kept in mind when comparing their results. For example, although Sects. 6.5 and 6.6 both look for a directional signal in the data, they are optimized for different forms of deviations from statistical isotropy. It is therefore unsurprising that they arrive at different results. However, the signal found in Sect. 6.5, if not simply a statistical fluke, is constrained by the results of Sect. 6.6.
Regarding the impact on the dipolar modulation results of the lack of the aberration contribution to the simulations, we note the following. In general the analyses are either sensitive only to large angular scales, or only claim possible detections on such scales, where the effect of aberration will be negligible and hence the conclusions are unlikely to change. A possible exception is in relation to the results of Sect. 6.5, where claims are made about effects extending out to ℓ_{max} = 1500. It is plausible that the effects of aberration could start to become important on these scales.
6.1. Variance asymmetry
The study of power asymmetry via the local variance of the CMB fluctuations was first performed by Akrami et al. (2014) for the Planck 2013 and WMAP 9year temperature data. The approach was motivated by its conceptual and implementational simplicity, its directly intuitive interpretation, and by virtue of being defined in pixel space, a useful complementarity to other mostly harmonicbased methods. The statistic was computed over patches of different sizes and positions on the sky, and compared with the values obtained from statistically isotropic simulations. It was found that none of the 1000 available simulations had a larger variance asymmetry than that estimated from the data. This suggested the presence of asymmetry at a statistical significance of at least 3.3σ, with a preferred direction (l,b) ≈ (212°,−13°) in good agreement with other studies. In this section, we revisit the variance asymmetry and report the results of the analysis for the Planck 2015 temperature maps at full resolution, N_{side} = 2048.
The analysis proceeds as follows. We consider a set of discs of various sizes centred on the pixels of a HEALPix map defined by a specific N_{side} value. For each sky map, we first remove the monopole and dipole components from the masked sky and then compute the variance of the fluctuations on a given disc using only the unmasked pixels. This yields a localvariance map at the HEALPix resolution of interest. We also estimate the expected average and variance of the variance on each disc from the simulations and then subtract the resulting average variance map from both the observed and simulated localvariance maps. Finally, we define the amplitude and direction of the asymmetry by fitting a dipole to each of the localvariance maps, where each pixel is weighted by the inverse of the variance of the variances computed from the simulations at that pixel. At all stages, we use only the discs for which more than 10% of the area is unmasked, although our results are robust against the choice of this value. The computed localvariance amplitudes are then used to compare the data with statistically isotropic simulations. Note that we use only the dipole amplitudes of the localvariance maps to measure the significance of the asymmetry; the amplitudes of higher multipoles were shown by Akrami et al. (2014) to be consistent with statistically isotropic simulations and we therefore do not consider them in the present paper.
In Akrami et al. (2014), the sensitivity of the method to the disc size was assessed using both statistically isotropic and anisotropic simulations. The free parameters, i.e., the number and size of the discs, were then fixed by these simulations. It was found that for 3072 patches centred on the set of pixels defined at N_{side} = 16, the simulated asymmetry signals were not detected when either very small (r_{disc}< 4°) or very large (r_{disc}> 16°) discs were used.
Fig. 27 Upper panel: pvalues for variance asymmetry measured as the number of simulations with localvariance dipole amplitudes larger than those inferred from the data, as a function of disc radius for the four componentseparated maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), and for unfiltered and highpassfiltered cases. For the filtered case, the Commander curve is covered by the SMICA curve for small (r_{disk} ≤ 8) disks, and by the SEVEM curve for large disks (r_{disk}> 8). Lower panel: localvariance dipole directions for the SMICA map. The colours, as indicated by the colourbar, correspond to different values of the highpass filter central multipole ℓ_{0}. The size of a marker disc corresponds, from small to large, to the size of the disc used in the analysis, namely 4°, 12°, 20°, and 70°. The dipole directions from the Commander, NILC, and SEVEM componentseparation methods are consistent with the case shown here. The lowℓ and WMAP9 directions are identical to those in Fig. 35. 
The former effect is due to a combination of the low number of pixels per disc and an insufficient number of discs to cover the entire sky when N_{side} = 16 reference grids are used. However, it has recently been shown by Adhikari (2015) that using a larger number of small discs (by increasing N_{side} to 32, 64, 128, and 256, depending on the disc size) in order to cover the entire sky allows the localvariance method to detect the largescale anomalous asymmetry as well as the Doppler boost signal from the Planck 2013 data, at a significance of > 3.3σ. Fantaye (2014) has demonstrated that the Doppler boost signal can be detected at a similar level of significance using needlet bandpass filtering of the data, even with large discs, when simulations are deboosted. Here, in contrast to the 2013 analysis, we use maps which contain Doppler boosting, for both simulations and data, and therefore we do not detect any Doppler boost signal when using a large number of small discs.
The low observed significance levels when large discs are used is due to the cosmic variance associated with the largestscale modes. Motivated by the analysis of Fantaye (2014), and in order to address this issue, we also perform analyses using a Butterworth highpass filter, (42)centred at multipoles ℓ_{0} = 5, 10, 15, 20, and 30. In addition, the filtering of low multipoles allows us to establish the contribution of such modes to any detected asymmetry.
Here, based on the analysis of Akrami et al. (2014), we restrict our analysis to those disc sizes for which 3072 discs, corresponding to an N_{side} = 16 map, cover the entire sky, i.e., to the range 4°–90°. Consistent results can be obtained by choosing other values of N_{side} for a given disc size provided that the entire sky is covered by the discs. Here, for simplicity, we work with the same N_{side} (= 16) for all disc sizes.
Our results for the measured amplitude of the variance asymmetry, compared to the values from the simulations, as well as the corresponding dipole directions, are shown in Fig. 27. The pvalues are given for different disc sizes and in terms of the number of simulations with localvariance dipole amplitudes greater than the ones measured from the data. Note that since the discs with different sizes used in our analysis are correlated, the significance levels are also correlated. For this reason we choose to show the pvalues as a function of disc size instead of combining them into a single number. Moreover, it should be noted that the significance values we present here do not incorporate any corrections to account for the choice of parameters adopted during method calibration, specifically the dipole amplitudes and directions for the anisotropic simulations that were used to fix the range of disc sizes and number of patches.
It can be seen from the upper panel of Fig. 27 that for the unfiltered map the significance of the power asymmetry drops quickly when we increase the disc size to radii greater than 16°. This is no longer the case, however, when the lowest multipoles are filtered out. For example, when the filter scale is set to ℓ_{0} = 5, i.e., when the very low multipoles which are affected most by cosmic variance are suppressed, the variance asymmetry is detected at the 3σ level for all disc sizes, as shown in Fig. 27. Table 20 presents the pvalues of the variance asymmetry using 8° discs and for various values of ℓ_{0}. Our results show that variance asymmetry is detected with a remarkable significance for all disc sizes when very low multipoles are filtered out. In addition, the variance asymmetry amplitude slowly decreases with increasing ℓ_{0}, as seen in the upper panel of Fig. 28. For ℓ_{0} ≳ 20, the dipole amplitude becomes too small and we find no significant variance asymmetry. It is interesting to note, however, that the dipole directions found for large ℓ_{0} are closely aligned with those found for ℓ_{0}< 20.
pvalues for the variance asymmetry measured by 8° discs for the four componentseparated temperature maps and different highpass filter scales.
The lower panel of Fig. 27 shows the dipole directions we find using different disc sizes and different filter scales for SMICA. The dipole directions for the Commander, NILC, and SEVEM componentseparated maps are very similar to those shown. The asymmetry directions found here are consistent with those determined by other analyses in this paper.
In the upper panel of Fig. 28, we show the localvariance dipole amplitudes for the 8° discs as a function of the central multipole of the highpass filter, ℓ_{0}. In the lower panel of the same figure we show, as an example, the meansubtracted and inversevarianceweighted localvariance map using 8° discs for the Commander componentseparation method. The pixels of the map are given in terms of the lower and uppertail probabilities of the values from the data compared to the values from the simulations. The maps for NILC, SEVEM, and SMICA are very similar. The numerical values of the localvariance dipole amplitudes and directions for the Commander method are given in Table 21; the values for the NILC, SEVEM, and SMICA methods are similar.
Fig. 28 Upper panel: localvariance dipole amplitude for 8° discs as a function of the central multipole of the highpass filter, ℓ_{0}, for the four componentseparation methods, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The grey regions, from dark to light, correspond, respectively, to 1σ, 2σ, and 3σ percentiles from the 1000 FFP8 simulations processed by the Commander method. Lower panel: meansubtracted and inversevarianceweighted localvariance map for the 8° discs and for the Commander componentseparation method; each pixel is given in terms of the lower and uppertail probability of the measured value on that pixel compared to the values from the simulations. The pixels in grey correspond to the centres of the 8° discs on which the number of unmasked pixels in the full resolution map is lower than our threshold. The black curve superposed on the map indicates the boundary of the opposing hemispheres along the asymmetry axis. It is clear that the largest fraction of >95% outliers (red pixels) lie on the positive amplitude hemisphere of the local variance dipole, while the <5% outliers (blue pixels) are on the opposite hemisphere. The corresponding maps for NILC, SEVEM, and SMICA are very similar to the one shown here. 
6.2. Dipole modulation: pixelbased likelihood
Localvariance dipole amplitudes and directions.
In PCIS13 we presented an analysis of the apparent anisotropic distribution of largescale power in the Planck 2013 temperature data within the parametric framework defined by Gordon (2007) and Hoftuft et al. (2009), who introduced an explicit dipole modulation field to model potential hemispherical power asymmetry. The following is a direct update of that analysis using the Planck 2015 CMB data at N_{side} = 32, retaining the 2013 common mask to explicitly test for consistency with the earlier study. All results are found to be in excellent agreement. In the following, we therefore only consider a smoothing scale of 5° FWHM as a representative example. This is the highest angular resolution accessible for an N_{side} = 32 map.
Recall first the basic data model adopted in the dipole modulation approach: rather than assuming the CMB sky to be a statistically isotropic Gaussian field, we allow for an additional dipole modulation, resulting in a data model of the form d = BMs + n, where is an offset dipole field multiplying an intrinsically isotropic signal s with a dipole of amplitude α pointing towards some preferred direction . B denotes convolution with an instrumental beam, and n denotes instrumental noise. Additionally, we model the power spectrum of the underlying statistically isotropic field in terms of a twoparameter amplitude–tilt model of the form , where is the bestfit Planck 2015 ΛCDM spectrum (Planck Collaboration XI 2016). The two parameters q and n can accommodate a deficit in power at low ℓ as compared to the bestfit cosmology that would otherwise create a tension with the underlying statistically isotropic model and result in the analysis measuring a combination of both asymmetry and power mismatch.
Fig. 29 Top: marginal constraints on the dipole modulation amplitude, as derived from Planck 2015 temperature observations at a smoothing scale of 5° FWHM for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The plot corresponds directly to Fig. 32 of Planck Collaboration XXIII (2014). The Commander, SEVEM, and SMICA posteriors coincide almost perfectly both internally, and with the corresponding SMICA 2013 posterior, shown as a dashed black line. Bottom: corresponding marginal twodimensional constraints on the lowℓ power spectrum amplitude and tilt, (q,n), defined relative to the bestfit Planck 2015 ΛCDM model. 
In the absence of any dipole modulation, α = 0, the total data covariance matrix is given by C = BS_{iso}B^{T} + N, where S_{iso} is the standard statistically isotropic CMB covariance matrix given by the power spectrum, C_{ℓ}, N is the noise covariance matrix, and the corresponding likelihood is given by the usual expression for a multivariate Gaussian distribution. With dipole modulation, this generalizes straightforwardly to C = BMS_{iso}M^{T}B^{T} + N, with the likelihood given by (43)Figure 29 and Table 22 summarize this fivedimensional likelihood in terms of marginal parameters for each of the four Planck CMB maps, as evaluated over the common mask using the multidimensional gridbased Snake algorithm (Mikkelsen et al. 2013). All results correspond to a smoothing scale of 5° FWHM, the highest resolution supported by an N_{side} = 32HEALPix grid, but, as in 2013, we consider all smoothing scales between 5° and 10° FWHM, reaching similar conclusions in each case: the dipole modulation results derived from the Planck 2015 temperature maps are essentially identical to the 2013 results, with improved internal consistency between the four CMB maps due to better mitigation of systematic errors. The bestfit dipole modulation amplitude at 5° FWHM is 6–7% whilst the lowℓ power spectrum has an approximately 3–5% lower amplitude compared to the bestfit ΛCDM prediction. These results are fully consistent with expectations given that the Planck 2013 sky maps were already cosmicvariancelimited on these angular scales, and the 2015 maps differ from the 2013 maps at the level of only a few microkelvin (Planck Collaboration IX 2016).
6.3. Dipole modulation: QML analysis
In this section we use the QML estimator introduced in Moss et al. (2011) and described in Appendix C to assess the level of dipole modulation in our estimates of the CMB sky at N_{side} = 2048. The specific implementation is essentially identical to that used in Hanson & Lewis (2009), Planck Collaboration XVII (2014), and Planck Collaboration XXVII (2014), and exploits the fact that dipole modulation of any cosmological parameter is equivalent to coupling of ℓ to ℓ ± 1 modes in the CMB covariance matrix to leading order (see Appendix C). Planck Collaboration XX (2016) presents an alternate analysis for a specific isocurvature model.
Since we are interested in dipole modulation there are three independent estimators. For our particular approach, these are a realvalued m = 0 and a complexvalued m = 1 estimator, and take the form Here T_{ℓm} are Cinverse filtered data and . We adopt the inversevariance filter from Planck Collaboration XVII (2014), where the approximate filter functions are also specified. We define δC_{ℓℓ + 1} ≡ dC_{ℓ}/dX + dC_{ℓ + 1}/dX, where X is the parameter modulated, and A_{ℓm} and B_{ℓm} are numerical coefficients (details can be found in Appendix C). The factor f_{1m} corrects the normalization for errors introduced by masking: (46)where M(Ω) is the mask. Finally, we correct the direction for the effects of inhomogeneous noise which is not accounted for in the filtering process, by weighting the by the inverse of the variance derived from filtered and meanfield corrected simulations.
The physics is readily accessible in this estimator: the ℓdependence in modulation determined by the parameter X is expressed in the δC_{ℓℓ + 1} factor, and the relevant scales appear directly in the limits of the sum. We consider the estimator over the range ℓ_{min} = 2 ≤ ℓ ≤ ℓ_{max}. The modulation amplitude and direction are then given by It is worth reemphasizing that the quantities , , and are all dependent on the ℓ range considered.
Amplitude (A) and direction of the lowℓ dipole modulation signal determined from the QML analysis for the range ℓ ∈ [2,64].
As a consequence of the central limit theorem, for sufficiently large ℓ_{max} the s are Gaussiandistributed with mean zero, so that the amplitude parameter has a MaxwellBoltzmann distribution. We fit to this distribution for ℓ_{max} ≥ 10 when computing the pvalue, so as not to be influenced by Poisson noise in the tails of the empirical distribution (and we have determined that this is a good fit to the simulations by applying a KS test). For the case of scalar amplitude modulation (i.e., X = A_{s}), and ℓ_{min} = 2, the cosmicvariancelimited expectation for the modulation amplitude from statistically isotropic skies is (50)This is the cosmic variance for a scaleinvariant dipole modulation, and gives a more explicit expression than the scaling discussed in Hanson & Lewis (2009).
Fig. 30 Probability determined from the QML analysis for a Monte Carlo simulation to have a larger dipole modulation amplitude than the Commander (red), NILC (orange), SEVEM (green), or SMICA (blue) data sets, with (top panel) ℓ_{min} = 2 or (bottom panel) ℓ_{min} = 100. No significant modulation is found once the lowℓ signal is removed. We emphasize that the statistic here is cumulative and apparent trends in the curves can be misleading. 
The top panel of Fig. 30 presents results for the pvalue of the fitted modulation amplitude as a function of ℓ_{max}. Note that there are several peaks, at ℓ ≈ 40 and ℓ ≈ 67 (the focus of most attention in the literature), and ℓ ≈ 240. The latter peak, while not previously emphasized, is also present in the WMAP results (see Fig. 15 in Bennett et al. 2011). It is also interesting to note that a modulation amplitude is observed at ℓ_{max} ≈ 800 that is somewhat lower than what one would typically expect for a statistically isotropic sky. However, the significance is not at the level of the excess dipole modulation at low ℓ and will not be discussed further. The dip at ℓ_{max} ≈ 67, with a pvalue of 0.9–1.0%, corresponds to the wellknown lowℓ dipole modulation^{6}. Table 23 presents the corresponding dipole modulation parameters, which are seen to be consistent with previous studies. Note that the mean amplitude expected for a set of statistically isotropic simulations at this ℓ_{max} is 2.9% (in close agreement with the expected value due to cosmic variance, Eq. (50)).
We have therefore determined a phenomenological signature of modulation for ℓ = 2–67 with a pvalue of 0.9–1.0%. If such a signal had been predicted by a specific model, then we could claim a significance of about 3σ. However, in the absence of such an a priori model, we can assess how often we might find a 3σ effect by chance, given that it could have occurred over any ℓ range. Since we are looking for a largescale phenomenon, we assume that the analysis should include the corresponding lowℓ modes and start at ℓ = 2. In order to correct for a posteriori effects we then adopt the following scheme.

1.
We calculate the modulation of each simulation on the scales 2–ℓ, where ℓ ∈ [10,ℓ_{max}]. For each simulation we find the modulation that gives thesmallest probability, η (in the same way that was done for the data).

2.
With the distribution of ηs given by the simulations we then compare this to the data. That is, we calculate the probability that one would find oneself in a Hubble patch with a modulation amplitude up to ℓ ∈ [10,ℓ_{max}] that is as significant as (or more significant than) the modulation in the real data.
If ℓ_{max} = 132 (as chosen by Bennett et al. 2011), the probability of achieving a modulation as large as the Planck data in this range is higher than 10% (see Fig. 31). This is in agreement with the findings of the WMAP team (which found 10% and 13% in the same ℓrange, using two different masks). Here, we do not quote a specific PTE for the dipole modulation since it depends on the choice of both ℓ_{max} (albeit not so sensitively) and ℓ_{min} (which we have decided not to marginalize over). However, it appears to be the case that the dipole modulation that we observe is quite unremarkable. That is, Gaussian fluctuations in a statistically isotropic Universe will reasonably often result in a dipole modulation with a comparable level of significance to that presented here.
Fig. 31 Probability determined from the QML analysis for obtaining a dipole modulation amplitude at least as anomalous as the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) data sets, for the range ℓ ∈ [10,ℓ_{max}]. The vertical line corresponds to ℓ_{max} = 132 which was used as the search limit in Bennett et al. (2011). The probability grows approximately logarithmically with ℓ_{max}. This means that the adopted probability to exceed is fortunately not very sensitive to ℓ_{max}, and for any reasonable choice is above 10%. 
Beyond this, evidence for dipole modulation is found at ℓ ≈ 200–300, with a smaller dip at ℓ ≈ 500. Given that the dipole modulation estimator is a cumulative quantity, it is possible that these features are statistically enhanced by the usual lowℓ signal. To test this we analyse the dipole modulation as a function of ℓ_{max} again, with the restriction ℓ_{min} = 100 applied in order to completely remove any lowℓ influence. The outcome is presented in Fig. 30 (bottom). It is clear that even before introducing posterior corrections no significant modulation is found, indicating that the pvalues of the features at ℓ> 100 were indeed exaggerated by the lowℓ modulation.
6.4. Bipolar spherical harmonics
In the absence of the assumption of statistical isotropy, the CMB twopoint correlation function can be most generally expanded in the bipolar spherical harmonic (BipoSH) basis representation as follows: (51)The BipoSH basis functions, are tensor products of ordinary spherical harmonic functions, and the corresponding expansion coefficients are termed BipoSH coefficients (Hajian & Souradeep 2003; Hajian & Souradeep 2006). The BipoSH basis provides a complete representation of any form of statistical isotropy violation with the key advantage of separating the angular scaledependence of the signal in spherical harmonic multipoles, ℓ, from the nature of the violation indexed in the bipolar multipole space by L. Consequently, it is possible to simultaneously determine that such a signal is dipolar (L = 1), quadrupolar (L = 2), octopolar (L = 3), and so on, in nature and that the power is restricted to specific ranges of angular scales.
The estimation of BipoSH coefficients from CMB maps is a natural generalization of the more routinely undertaken estimation of the angular power spectrum C_{l}. To allow a direct connection to the angular power, we further introduce a set of BipoSH spectra at every bipolar harmonic moment, (L,M), labelled by a difference index d, defined as follows: (52)where are the ClebschGordon coefficients and for brevity the notation . BipoSH spectra, clearly, are then simply a generalized set of CMB angular power spectra, with the standard CMB angular power spectrum being one of them^{7}. While quantifies the properties of the statistically isotropic part of the CMB fluctuations, the additional BipoSH coefficients quantify the statistically anisotropic part of the CMB twopoint correlation function.
Amplitude (A) and direction of the dipole modulation in Galactic coordinates as estimated for the multipole range ℓ ∈ [2,64] using a BipoSH analysis.
Thus BipoSH provides a mathematically complete description of all possible violations of statistical isotropy in a Gaussian CMB sky map. It is then always possible to translate any specific model for such a signal into the language of BipoSH and provide a common approach for the multiple specialized tests that have been implemented previously in this paper and elsewhere. However, improving on the analysis of the 2013 Planck data, a new formalism is developed in order to reliably analyse a masked sky, as concisely described in Appendix D. Aluri et al. (2015) provides a more detailed description of the approach and includes an explicit demonstration of its validity using simulations.
Initially, we revisit the simple phenomenological model of dipole modulation of the CMB sky from Sect. 6.2, (53)where represents the modulated CMB sky, is the underlying (statistically isotropic) random CMB sky, and is a dipolar field. The BipoSH coefficients resulting from such a modulation are given by Here corresponds to the BipoSH coefficients of the unknown, but statistically isotropic, unmodulated CMB field, m_{1M} are the spherical harmonic coefficients of the modulation field, and C_{ℓ} is the bestfit CMB angular power spectrum.
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 ℓ in Eq. (55). This additional information is important for identifying the origin of the isotropybreaking signal, which could be either cosmological or due to systematic artefacts.
We perform the analysis for the N_{side} = 2048 component separated CMB maps with an apodized version of the common mask at that resolution and reconstruct the modulation signal in independent bins of width Δℓ = 64 up to ℓ_{max} = 512. The application of the common mask introduces a mean field bias in the BipoSH coefficients derived from the data. This bias is estimated from the FFP8 simulations and subtracted from the derived coefficients. The process of masking induces a coupling between the modulation field and the mask that results in a modification of the spectral shape of the modulation signal by the modified shape function (MSF; see Appendix D for details). Further, the covariance of the biassubtracted BipoSH coefficients is not easy to derive analytically in this case. To overcome this problem, we consider the diagonal approximation to the covariance matrix and estimate it from simulations.
The results presented in the top panel of Fig. 32 indicate that the dipole modulation signal is most significant in the lowest multipole window ℓ ∈ [2,64]. Note that the power in the dipole modulation field m_{1} = (  m_{11}  ^{2} +  m_{10}  ^{2} +  m_{1−1}  ^{2}) / 3 is related to the dipole amplitude by . The bestfit amplitude (A) and direction corresponding to the reconstructed dipole modulation field from this lowest multipole bin is quoted in Table 24 for each componentseparation method. Also shown are the corresponding results for the cleaned frequency maps SEVEM100, SEVEM143, and SEVEM217. As expected for signals with a cosmological origin, no evidence for frequency dependence is seen.
Fig. 32 Top: measured dipole modulation (L = 1) power in nonoverlapping CMB multipole bins for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) as determined from a BipoSH analysis of the data. The power in the dipole of the modulation field is a χ^{2}distributed variable with 3 degrees of freedom. The shaded regions in the plot depict, in darkgrey, grey, and lightgrey respectively, the 1, 2, and 3σ equivalent intervals of the distribution function derived from simulations, while the solid black line denotes its median. Significant power in the dipole modulation is seen to be limited to ℓ = 2–64 and does not extend to higher multipoles. Bottom: dipole modulation direction as determined from the SMICA map. The directions found from the other component separation maps are consistent with this analysis. The coloured circles denote the central value of the multipole bin used in the analysis, as specified in the colour bar. The lowℓ and WMAP9 directions are identical to those in Fig. 35. 
Since the amplitude of the dipole modulation field is consistent with zero within 2σ for all of the higher ℓbins considered, it is plausible that the simple modulation model in Eq. (53) is inadequate to describe the features seen in the BipoSH spectra 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. It is also intriguing to note that, although in most cases the amplitude of the modulation dipole is seen at low significance, the directions in the first four bins, ℓ_{32} ∈ [2,64], ℓ_{96} ∈ [65,128], ℓ_{160} ∈ [129,192], and ℓ_{224} ∈ [193,256], are seen to be clustered together, as shown in the bottom panel of Fig. 32. Note that the lower significance of the modulation for the multipole bins at ℓ> 64 results in larger errors for their respective directions than the value quoted for the ℓ ∈ [2,64] bin recorded in Table 24.
We extend our analysis to carry out the dipole modulation reconstruction in cumulative bins up to ℓ_{max} = 512, making cumulative increments in the multipole in steps of Δℓ = 64. The results of this analysis are summarized in Fig. 33.
Fig. 33 Top: measured dipole modulation power in cumulative CMB multipole bins for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) as determined from a BipoSH analysis of the data.. Colour coding as in Fig. 32. Note that the measurements in cumulative bins indicate a power in excess of 2σ up to multipole ℓ_{max} ~ 320. The value on the horizontal axis denotes the maximum multipole used in the analysis, with ℓ_{min} = 2. Bottom: modulation dipole direction as recovered from the SMICA map. The directions found from the other componentseparation maps are consistent with these directions. The colourcoded points represent the directions recovered for the specific ℓ_{max} used in the analysis, with ℓ_{min} = 2. The lowℓ and WMAP9 directions are identical to those in Fig. 35. 
As noted previously, as a consequence of our motion with respect to the CMB rest frame, the observed CMB map is expected to be statistically anisotropic, as has been demonstrated in Planck Collaboration XXVII (2014) and Appendix B. Reassuringly, in PCIS13 it was established that such a signal would not contaminate a dipole modulation signal up to ℓ_{max} ≈ 700. We now confirm the Doppler boost signal using the BipoSH methodology.
An equivalent description of the Doppler boost in terms of BipoSH coefficients is given by where , β = v/c denotes the peculiar velocity of our local rest frame with respect to the CMB, and b_{ν} is the frequencydependent boost factor, as discussed in more detail in Planck Collaboration XXVII (2014).
Since the Doppler boost signal has a frequency dependence, we perform our analysis on the SEVEM100, SEVEM143, and SEVEM217 maps at N_{side} = 2048, and adopt values of b_{ν} = 1.51,1.96, and 3.07, respectively. A minimum variance estimator for β_{1M}, as discussed in Appendix D, is adopted with the shape function replaced by the corresponding Doppler boost term given in Eq. (56). Corresponding unboosted CMB simulations were also used, in particular to correct for the mean field bias. However, we use a set of Dopplerboosted simulations in order to estimate the error on the reconstructed Doppler boost vector.
Since it is expected that the low multipole modes of the spectrum are contaminated by the dipolar signal reported previously, in order to monitor the impact of this anomalous signal on the Doppler reconstruction we implement a cumulative analysis using multipoles with a varying ℓ_{min} from 2 to 640 in increments of Δℓ_{min} = 128 and a fixed ℓ_{max} = 1024^{8}. The recovered Doppler amplitudes from the three SEVEM frequency cleaned maps as a function of ℓ_{min} are shown in the top panel of Fig. 34, while the lower panel indicates the corresponding direction in Galactic coordinates determined from the SEVEM217 data. Table 25 records the bestfit amplitudes and directions for ℓ ∈ [640,1024].
Fig. 34 Top: amplitude  β  of the Doppler boost from the SEVEM100, SEVEM143, and SEVEM217 maps for different multipole bins determined using a BipoSH analysis. The maximum multipole of each bin is fixed at ℓ_{max} = 1024, while ℓ_{min} is incremented from ℓ = 2 to ℓ = 640 in steps of Δℓ = 128. The dashed line corresponds to the actual dipole boost amplitude,  β  = 1.23 × 10^{3}. Bottom: Doppler boost direction measured in Galactic coordinates from SEVEM217. The coloured circles denote ℓ_{min} used in the analysis, while ℓ_{max} = 1024 is held fixed. The lowℓ and WMAP9 directions are identical to those in Fig. 35. 
Doppler boost amplitude ( β ) and direction in Galactic coordinates derived over the multipole range ℓ ∈ [640,1024] as evaluated from a BipoSH analysis.
6.5. Angular clustering of the power distribution
In the Planck 2013 data release we reported a possible deviation from statistical isotropy in the multipole range ℓ = 2–600, thus confirming earlier findings based on the WMAP data (Hansen et al. 2009; Axelsson et al. 2013). This claim of asymmetry extending to higher multipoles was made only on the basis of the alignment of preferred directions as determined from maps of the power distribution on the sky for specific multipole ranges. In particular, it was found that the directions of the dipoles fitted to such maps in the multipole range ℓ = 2–600 were significantly more aligned than in simulations. In addition, we showed that the ratio of the power spectra in the two opposite hemispheres defined by the asymmetry axis for ℓ = 2−600 was not statistically anomalous (as later confirmed over the extended multipole range ℓ = 2−2000 by Quartin & Notari 2015).
Fig. 35 Dipole directions for independent 100multipole bins of the local power spectrum distribution from ℓ = 2 to 1500 in the SMICA map with the common mask applied. We also show the preferred dipolar modulation axis (labelled as “lowℓ”) derived in Sect. 6.2, as well as the total direction for ℓ_{max} = 600 determined from WMAP9 (Axelsson et al. 2013). The average directions determined from the two multipole ranges ℓ ∈ [2,300] and ℓ ∈ [750,1500] are shown as blue and red rings, respectively. The error on the derived direction that results from masking the data is about 60°, with only small variations related to bin size. 
Here, we test for the alignment in the Planck 2015 data set. We adopt the approach for the estimation of the dipole alignment that was described in detail in PCIS13, a brief summary of which follows.

1.
Local power spectra are estimated from the data at N_{side} = 2048 for 12 patchesof the sky corresponding to the N_{side} = 1HEALPix base pixels. Only thosehighresolution pixels surviving the application of the commonmask are included in the analysis^{9}. As a consequenceof this masking, when patches based on HEALPix pixels with N_{side}> 1are used, the available sky fraction for those patches close to theGalactic plane is too small for powerspectrum estimation. Formost of the analysis, we use the crossspectra determined fromhalfmission data sets^{10}. Due to a mismatch betweenthe noise level in the data and the simulated maps, the results basedon autospectra are less reliable and also more prone to othersystematic effects than the crossspectra. We therefore do notconsider such results here. The spectra are binned over various binsizes between Δℓ = 8 and Δℓ = 32.

2.
For each power spectrum multipole bin, an N_{side} = 1HEALPix map with the local power distribution is constructed.

3.
The bestfit dipole amplitude and direction are estimated from this map using inversevariance weighting, where the variance is determined from the local spectra computed from the simulations. We do not compute error bars for the direction, but expect this to be accounted for in part by the use of equivalently treated simulations in the clustering analysis.

4.
A measure of the alignment of the different multipole blocks is then constructed. In PCIS13, we considered the mean angle between all possible pairs of dipole directions up to a given ℓ_{max}. Here, for greater consistency with Sect 6.6, we use the mean of the cosine of the angles, rather than of the angles themselves, between all pairs of dipoles. This effectively corresponds to the Rayleigh statistic (RS) introduced formally in Sect. 6.6, and we will refer to it as such, although it differs by ignoring all amplitude information. Clearly, smaller values of the RS correspond to less clustering.

5.
The clustering as a function of ℓ_{max} is then assessed using pvalues determined as follows. We first construct the RS using all multipoles up to ℓ_{max}. The pvalue is then given by the fraction of simulations with a higher RS than for the data for this ℓ_{max}. A small pvalue therefore means that there are few simulations that exhibit as strong clustering as the data. Note that the pvalues are highly correlated as the RS is a cumulative function of ℓ_{max}.

6.
We then define two measures of significance. To achieve this, it is necessary to reduce the 1499 different pvalues determined for ℓ_{max} ∈ [2,1500] to a single measure of clustering. We do this in two different ways, using the mean of these pvalues, and by finding the minimum of the pvalues, for both the data and for each available simulation. We then determine the percentage of simulations with (i) a lower mean pvalue and (ii) a lower minimum pvalue than the data. Note that these two measures of significance take into account different aspects of the data. Note further that since the RS is cumulative and the pvalues therefore correlated, different scales are weighted unequally and a detection in the mean and/or minimum pvalue may be difficult to interpret and to correct for the multiplicity of tests effect (LEE).
Note that the statistics defined in step 6 above correspond to two choices of what were referred to as “global statistics” in PCIS13 in order to assess the degree to which the significance of the results depends on a specific choice for ℓ_{max}. The mean pvalue over all available ℓ_{max} measures the degree to which clustering is present over large multipole ranges independently of whether the clustering is strongly focused in one given direction. Clearly the pvalues for different ℓ_{max} are strongly correlated, but if the clustering is present only over a small multipole range, the RS will drop and the corresponding pvalues will eventually rise. By comparing this value to simulations, we test not only whether the dipole alignment in the data is stronger than in statistically isotropic random simulations, but also whether it is present over larger ranges of multipoles than expected. The minimum pvalue will give strong detections if there is a strong asymmetry over a limited multipole range or weaker clustering over larger multipole ranges when the clustering is strongly focused in a given direction.
For Commander, NILC, and SEVEM, only 500 simulations are available. However, 5000 simulations are available for SMICA, which allows a better estimate of significance to be determined when the probabilities obtained are very low. In this case, we use half of the 5000 simulations to calibrate the statistic (obtain pvalues following step 5 above) and the remaining half to determine significance levels (compute the mean and minimum over these pvalues as a function of ℓ_{max} following step 6). When using 500 simulations, it is necessary to use the same set of simulations to calibrate as well as to obtain probabilities. A related issue with these results is that this set of simulations (corresponding to the first 500 out of the 5000 available for SMICA) are observed to yield higher pvalues for the clustering angle due to a statistical fluctuation. Another 9 sets of 500 simulations that can be obtained from partitioning the 5000 available SMICA simulations all result in lower pvalues. As a consequence, we observe that results based on the larger number of simulations often give lower pvalues than when only 500 simulations are used.
In Fig. 35 we show the dipole directions of the 15 lowest 100multipole bins for the SMICA map. Here, the binning has been chosen for visualization purposes; in further analysis of the Planck data we use finer ℓintervals. The preferred lowℓ modulation direction determined in Sect. 6.2 is also indicated, along with the WMAP9 result determined over the range ℓ = 2 to 600 (Axelsson et al. 2013). The observed clustering of the dipole directions is similar to that shown in figure 27 of PCIS13. Note that differences in masking, foreground subtraction, and residual systematic effects will displace the direction of a given dipole with respect to the previous analysis. Similar behaviour is seen for all of the Planck componentseparated maps.
Fig. 36 Derived pvalues for the angular clustering of the power distribution as a function of ℓ_{max}, determined for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), based on 500 simulations. For SMICA, the pvalues based on 2500 simulations are also shown (black). The pvalues are based on the fraction of simulations with a higher RS, determined over the ℓrange up to the given ℓ_{max}, compared to the data. The results shown here have been marginalized over bin sizes in the range Δℓ = 8 to Δℓ = 32. 
In PCIS13, we calculated the mean angle between all possible pairs of dipole directions determined from maps of the local power in multipole bins of size Δℓ = 16. Here we test the possible bias arising from such a choice by considering bin sizes between Δℓ = 8 and Δℓ = 32 in steps of 2. The lower limit avoids significant bintobin coupling in the power spectra for smaller binnings, whilst the upper limit excludes cases where there are an insufficient number of derived dipoles from which the mean angle can be calculated, this leading to poor statistics. In addition to showing results for each bin size, we also calculate the varianceweighted mean of the power spectra over all bin sizes (the C_{ℓ} for a given bin size is weighted by where N_{b} is the bin size). In this way, we marginalize over bin sizes to obtain local power spectra and thereby the RS for each single multipole.
Figure 36 shows the pvalues for the different componentseparated maps, derived as described in step 5 above. We see that the results based on 500 simulations for NILC, SEVEM, and SMICA are in good agreement. The Commander results are less consistent, but this may be related to the fact that component separation was performed independently for the halfmission solutions, in contrast to the other methods, where componentseparation solutions were obtained from the full mission data only. For SMICA, we also show pvalues based on 2500 simulations. These more accurate results show lower pvalues, and may indicate that those determined from only 500 simulations are not sufficiently stable. Note also that for ℓ< 100 the pvalues are not consistent with the detection of a lowℓ asymmetry/modulation, as seen by other methods in this paper. However, for ℓ< 100, there are very few bins and the variance of the RS might therefore be too high for this effect to be visible.
In agreement with the conclusions in PCIS13, a large degree of alignment is seen at least to ℓ_{max} ≈ 600. However, in contrast to the earlier results where the pvalues started increasing systematically for ℓ_{max}> 1000, the current pvalues remain low for ℓ_{max}> 750. The full componentseparated maps which have higher resolution and sensitivity are used for the current analysis, instead of the singlefrequency foregroundcleaned map (SEVEM143) used in PCIS13. We note that the results for the updated SEVEM143 map are consistent with the earlier analysis, both with and without correction for the Doppler modulation. Note also that the SMICA results with improved statistics (based on 2500 simulations) generally show lower pvalues than the corresponding results based on 500 simulations.
Significance of the angular clustering of the power distribution.
Table 26 presents the fraction of simulations with a lower mean/minimum pvalue than in the data for a number of different cases. The table shows probabilities for SMICA with different bin sizes (showing only every second bin size since these are correlated), as well as for the results marginalized over bin sizes. We also show results for the different componentseparated maps, results based on halfring crossspectra instead of halfmission crossspectra, and results using a different ℓweighting scheme, specifically (2ℓ + 1)C_{ℓ} instead of ℓ(ℓ + 1)C_{ℓ}, the former being a measure of the variance of the temperature fluctuations. The table indicates probabilities of approximately 0–2% for most of these cases, although results for the smallest bin size show much less significant results. This could be due to the strong anticorrelations between adjacent bins found for this bin size in those Galactic N_{side} = 1 patches with very small available sky fraction. For the other bin sizes, these correlations are much weaker. Note that many of the significances based on minimum pvalue are only upper limits. This is due to the fact that the limited number of simulations in some cases results in the lowest minimum pvalue being zero. When the minimum pvalue in the data is zero, we show the percentage of simulations which also have zero as the minimum pvalue. Clearly this fraction is only an upper limit on the real significance.
Fig. 37 Derived pvalues for the angular clustering analysis as a function of ℓ_{max}, determined from SMICA based on 2500 simulations. The pvalues are based on the fraction of simulations with a higher Rayleigh statistic up to the given ℓ_{max} than in the data. The RS here is calculated over all pairs of dipole directions where one dipole in each pair is computed in the range [ℓ_{lim},ℓ_{max}], and the other is determined in the range [2,ℓ_{lim}]. The plot shows pvalues for ℓ_{lim} = 300 (purple), ℓ_{lim} = 400 (yellow), ℓ_{lim} = 500 (pink), and ℓ_{lim} = 700 (cyan). The results have been marginalized over bin sizes in the range Δℓ = 8 to Δℓ = 32. 
In order to further investigate the ℓdependence of the asymmetry, we follow two approaches from PCIS13. Firstly, we restrict the analysis to multipoles above a minimum multipole ℓ_{min}. Table 26 indicates that clustering at the < 1% significance level is still found when considering only those multipoles with ℓ_{min} greater than 100. However, when this limit is increased to 200, no significant clustering is found. We then calculate the RS between pairs of dipoles where one dipole is determined from an ℓrange above a certain limiting multipole ℓ_{lim}, and the other dipole below this limit. Figure 37 shows the RS as a function of ℓ_{max} for some selected values of ℓ_{lim}. The ℓ_{lim} = 300 curve (purple) indicates that dipole directions for ℓ> 1000 are significantly aligned with dipoles for ℓ< 300. Similarly, the ℓ_{lim} = 700 curve (cyan) indicates that the dipole directions for ℓ = 700–1000 are strongly correlated with the dipole directions for ℓ< 700.
Combining these results, we note that when using only multipoles with (i) ℓ> 200; or (ii) ℓ< 200, no significant clustering is found. The strong clustering significance shown to persist to high multipoles in Fig. 36 must therefore be the result of clustering of the dipole directions between low and high multipoles as supported by Fig. 37. The low pvalues can be explained by the alignment of dipole directions for multipoles extending all the way to ℓ = 1500 correlated with directions for ℓ< 200. The observed asymmetry is therefore not consistent with a model based on dipole modulation or power asymmetry located in one specific multipole range or for one given direction, but rather as a correlation of the dipole directions between ℓ< 200 and ℓ> 200. This correlation with lower multipoles is found to persist all the way to ℓ_{max} = 1500.
An advantage of the directional analysis performed here is that it focuses on a central issue for tests of deviation from isotropy – whether there is a preferred direction. Indeed, Bunn & Scott (2000) noted that the CMB may exhibit a pattern that cannot be identified from the power spectrum, but which would indicate some nontrivial largescale structure. Evidence for the close correlation and alignment of directions on different angular scales may present a signature of broken statistical isotropy, since in the standard model, these directions should all be independent random variables. In this context, we do not quote a specific direction for such asymmetry here since our results indicate a clustering of angles between different multipoles, but not necessarily that all multipoles are clustered about one specific direction. However, crucially we have shown that the measured clustering is driven by the correlations of directions between higher and lower multipoles.
Some of the analyses in other sections of the paper focus on dipolar modulation, a specific model for a dipolar power enhancement of the statistically isotropic CMB field towards a preferred direction of the sky, and use methods optimized for the detection of such a signal. While the results of Sect. 6.6 show no detection of the clustering of directions, there is no clear contradiction with the results presented here, since they are based on tests for a_{ℓm} correlations between different multipoles as expected in the dipolar modulation model. The clustering analysis presented here is a modelindependent test for deviations from statistical isotropy which could induce very different correlation structure. It is therefore sensitive to other forms of asymmetry, such as the addition of power in one part of the sky or more general phase correlations.
6.6. Rayleigh statistic: QML analysis
Results from Sect. 6.5 and in PCIS13 suggest that, beyond a dipole modulation of power on large angular scales, some form of directional asymmetry continues to small scales. There are also indications from Sect. 6.5 that the directions of dipolar asymmetry are correlated between large and small angular scales. Since the nature of the asymmetry is unknown we use the RS, a generic test for directionality that makes minimal assumptions about the nature of the asymmetry. This statistic has been used both in previous CMB studies (Stannard & Coles 2005) and other areas of cosmology (Scott 1991). In our context, for a statistically isotropic sky this statistic is identical to a threedimensional random walk. The implementation here incorporates all information pertaining to modulation, not just the direction. The approach in this section differs from that of Sect. 6.5 in the method of reconstructing power, the choice of binning, and the choice of how to weight directions in each bin. Another important difference is that Sect. 6.5 only considers the direction of dipolar asymmetry and does not take into account its amplitude.
The statistic is cumulative and thus narrowing down the specific scales from which a signal may be originating is a nontrivial task. However, it is the case that all statistics that measure this form of asymmetry (dipole modulation or largescale clustering of power) are in some way cumulative and so we will not worry about this issue any further. Another disadvantage of this approach is that it will generally be less powerful than a test that uses a specific model for the directionality. Again, this is a distinction shared when one compares any nonparametric versus parametric statistic.
The construction of the statistic is as follows.

1.
Beginning with the estimator from Eqs. (44)and (45), compute the followingbinned quantities for the data and simulation:For each ℓ this computes the coupling of ℓ to ℓ + 1. We emphasize that this is a very natural choice of binning the estimator, since any parameter that is dipole modulated will lead to coupling of ℓ to ℓ ± 1 modes, albeit with different ℓweightings (below we describe why this is not an important issue).

2.
Construct a threedimensional vector out of the three estimators for both the data and the simulations^{11}, as defined by Eqs. (47)−(49).

3.
Compute the mean amplitude from simulations and divide all vectors (data and simulations) by this amplitude. This choice ensures that each vector is treated equally, since we have no a priori reason to weight some scales more than others.

4.
Add this new vector to the previous vector. If this is the first time going through this process the previous vector is the zero vector.

5.
Repeat with ℓ → ℓ + 1. Note that the statistics of this process are identical to a three dimensional random walk.
Given that a dipole modulation amplitude of roughly 3σ significance is known to exist at low ℓ (before a posteriori correction), one would expect a similar level of detection of asymmetry to be determined by the RS. Indeed, we find that asymmetry is present out to ℓ ≈ 240. Figure 38 (top) presents the pvalues derived when the RS is computed as a function of ℓ_{max} from ℓ = 2. The minimum pvalue obtained by the data is 0.1–0.2%, to be compared to the value of 0.9–1.0% obtained for the dipole modulation amplitude at ℓ_{max} = 67. The direction preferred by the data for ℓ_{max} ≈ 240 is (l,b) = (208°,−29°), which is approximately 20° away from the dipole modulation direction determined to ℓ ≈ 64.
Fig. 38 Rayleigh statistic pvalues determined from the QML analysis as a function of ℓ_{max} for the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) data sets, with (top panel) ℓ_{min} = 2 and (bottom panel) ℓ_{min} = 100. The general pattern of peaks is very similar to that in Fig. 30. We emphasize that the statistic here is cumulative and as such trends in the curves can be misleading. 
Fig. 39 Probability to exceed (PTE) the pvalue of the signal from the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) data at ℓ = 230–240 (which is the multipole range with the most significant deviation) when searching over a range of multipoles up to ℓ_{max}, for the RS determined from the QML analysis. Much like the equivalent curve for dipole modulation, the PTE appears to grow approximately logarithmically with ℓ_{max}. 
We correct for a posteriori statistics using the same procedure as in Sect. 6.3. Specifically, we count how often simulations find asymmetry in the range 10 ≤ ℓ ≤ ℓ_{max} that is more significant than that found for the data. From Fig. 39 it is clear that generic asymmetry at the significance level found in our CMB sky occurs about 6% or 8% of the time (depending on the range of ℓ one decides to search over).
While the PTE here is not very low, it is nevertheless somewhat lower than for the usual dipole modulation test. Hence, it seems worth exploring whether any of this signal comes from higher multipoles. Therefore we compute the RS starting at ℓ_{min} = 100, to avoid the influence of asymmetry at lower ℓ. The lower panel of Fig. 38 presents the corresponding pvalues as a function of ℓ_{max}. There is a striking similarity with the lower panel of Fig. 30. It is clear that, even in the absence of a posteriori correction, we find no significant asymmetry at larger ℓ. Hence most of the signal we are seeing in Fig. 38 (top) is due to the usual lowℓ asymmetry.
We would like to stress that the results here are very similar to the results of the previous section. For each of the statistics used we are simply asking whether there is significant coupling of ℓ with ℓ ± 1 modes. The details of how to optimally combine these couplings for a given ℓ range depends on whether we are talking about dipole modulation or directionality (or some other related test, e.g., variance asymmetry). These details will change the range of scales over which the strongest signal in the data is found.
7. Sensitivity of anomalies to enhanced sky coverage
One of the critical aspects in searching for anomalous features in sky maps is to ensure that the region being investigated constitutes a fair and unbiased sample. Since many of the claimed anomalies are on large angular scales, this implies that minimal masking should be applied to the data. However, residual foregrounds then become a significant consideration. The masks applied to the four componentseparated maps studied in the bulk of this paper have been defined at high resolution, and then conservatively degraded for lower resolution studies. Such a procedure inevitably reduces the sky coverage available for analysis, and can be particularly problematic if significant structures are aligned by chance with the masked regions. Indeed, the WMAP team (Bennett et al. 2011) have drawn attention to several such features in their ILC reconstruction of the CMB sky, and these are clearly also present in the PlanckCommander, NILC, SEVEM, and SMICA sky maps. A large cold spot is seen near to the Galactic centre, a significant fraction of which lies within the common mask at any resolution. However, despite its location and visual impression, the feature is neither likely to be attributable to residual foreground emission, nor is it inconsistent with the ΛCDM model (Gott et al. 2007). In addition, four elongated cold fingers stretching from near the Galactic equator to the south Galactic pole are seen, although no equivalent features are evident in the northern sky. Bennett et al. (2011) have noted that the alignment of the ℓ = 2 and ℓ = 3 multipoles (Tegmark et al. 2003) seems to be intimately connected with these largescale cool fingers and the intervening warm regions. One of the latter also corresponds to the wellknown “Bianchi VII_{h}” main lobe originally found in Jaffe et al. (2005).
Although we would ideally pursue full sky analyses, we prefer to remain mindful of the influence of residual foregrounds, but still seek to minimize the extent of any mask applied for analysis. In this context, and specifically for largeangularscale studies, we consider the properties of an additional estimate of the CMB sky, also generated using the Commander component separation methodology. In particular, we note that the Planck lowℓ likelihood analysis (Planck Collaboration XI 2016) uses the temperature solution from this study, degraded to a resolution of N_{side} = 16. The LklCommander map, as we now refer to it, is initially derived from input data sets (32 bands) at 1° FWHM resolution and N_{side} = 256. This includes Planck individual detector and detector set maps from 30–857 GHz, the 9year WMAP observations between 23 and 94 GHz, and the 408 MHz sky survey (Haslam et al. 1982), whereas the Commander map described in Planck Collaboration IX (2016) includes Planck data alone. It is believed that the 32band solution is better (on large angular scales) than the Planckonly map, because the larger number of input frequencies allows more detailed foreground modelling, and in particular the separation of the lowfrequency foregrounds into synchrotron, freefree, and spinning dust components. An associated confidence mask (hereafter LklT_{256}93) is then defined based on a goodnessoffit measure per pixel, corresponding to a rejection of 7.3% of the pixels on the sky. A detailed discussion of these results can be found in Planck Collaboration X (2016).
Lowertail probability for the variance, skewness, and kurtosis of the LklCommander map.
We now consider the implications of using the LklCommander map for studies of several largeangularscale anomalies observed in previous sections, in particular since the larger sky coverage permitted by this data set should constitute a better sample of the Universe. Note that, at the resolutions of interest for the following analyses, the noise level is negligible (even accounting for the WMAP contribution) and should not have significant impact on the results. The exact details of the noise contribution to simulations is therefore unimportant.
7.1. Variance, skewness, and kurtosis
We begin by estimating the variance, skewness, and kurtosis of the CMB. We apply the unit variance estimator to the LklCommander map, and specifically to the version used in the lowℓ likelihood analysis, which is smoothed to 440′ FWHM at a resolution of N_{side} = 16. A corresponding lowℓ mask is generated by a simple degrading of the mask at N_{side} = 256, then setting those N_{side} = 16 pixels with a value less than 0.5 to zero and all others to unity. The resultant lowℓ likelihood mask rejects only 6.4% of the sky. We compare the results for both this mask (also to be referred to as LklT_{16}94), and the standard common mask at this resolution (UT_{16}58). The results are summarized in Table 27 and show that, when using the lowℓ likelihood mask, the lower tail probability for the variance is 7.0%. This value is higher than the corresponding values for the component separated maps as shown in Table 12. In addition the skewness and kurtosis are less consistent with Gaussianity than the component separated maps. However, when using the standard common mask at N_{side} = 16, the lower tail probability of the variance, skewness, and kurtosis become more compatible with those derived earlier.
There are two possible explanations for this behaviour. Either the variance of the CMB in the region close to the Galactic plane is intrinsically high, perhaps due to the presence of the various features noted above, or the presence of residual foregrounds increases the variance of the map. In order to attempt to distinguish between these options, we again apply the unit variance estimator to the standard componentseparated maps^{12}, but this time utilising the lowℓ mask. Although the componentseparated maps are likely to contain some foreground contamination in the regions omitted by application of the UT_{16}58 mask, it is appropriate to recall that this was constructed in a conservative way, and may also mask parts of the sky where the level of residual foregrounds can be considered negligible. In addition, we investigate the cleaned frequency maps produced by the SEVEM algorithm in order to test for the presence of frequencydependent residual foregrounds. The results of the unit variance estimator analysis are summarized in Table 28.
Lowertail probability for the variance, skewness, and kurtosis of the LklCommander map compared to the component separated maps, obtained using the lowℓ likelihood mask LklT_{16}94.
All of the component separated maps show an increase in the lower tail probability from about 0.5% when the UT_{16}58 mask is applied to roughly 7% for the LklT_{16}94 mask. The small variations in results for the different maps may be attributable to the presence of residual foregrounds close to the Galactic plane. However, the increased probabilities can also be explained by the presence of CMB structures with higher variance within that region which is not rejected by the less conservative mask. Indeed, since the componentseparated maps are affected by different residual foregrounds, if the source of the changes in probabilities is due only to the residual foregrounds, then we would expect a larger dispersion than what is observed. We also note that when we apply the lowℓ likelihood mask the skewness and kurtosis values are shifted towards more extreme values. This implies that the sky signal is less Gaussian for the larger sky fraction, despite the results remaining compatible with the ΛCDM model assumed for the null tests. Both Commander maps are noteworthy in this regard.
An important issue is whether the changes in the statistics can simply be attributed to differences in the masks. We determine how many simulations show an increase in variance at least as large as that seen for the LklCommander map when comparing the values derived for the UT_{16}58 and lowℓ likelihood masks. Similarly, we determine how many simulations have increased skewness or kurtosis values with shifts at least as large as observed. When the three statistics are considered separately, the fraction of simulations that indicate such changes are 7.6%, 4.3%, and 13.9% for the variance, skewness, and kurtosis, respectively. Of course, such subsets of the simulations also include cases where a large shift in the statistic is observed, but the statistic would not be considered anomalous for either mask. If we also impose the requirement that the simulations have these shifts for all three quantities simultaneously, then only 2 maps from 1000 are found. Of course, such a requirement is rather strong, and at this stage we are likely to be approaching the limits of what can be said based on modelindependent null tests. Indeed, in order to assess whether these results are sensitive to a posteriori choices, we repeat the analysis but successively take each simulation as the reference. Thus, for each simulation the shift in the variance, skewness, and kurtosis is computed and then we determine how many times we find a case in which two or less of the remaining simulations simultaneously show larger shifts for the three moments. We find that 48 maps from 1000 satisfy these conditions. Given this, it is difficult to draw strong conclusions about the significance or otherwise of the maskrelated changes in variance.
7.2. Npoint correlation functions
The connection between sky coverage and the observed structure of the 2point correlation function for large angular separations has previously been discussed in the literature, in particular in connection with the S_{1 / 2} statistic discussed in Sect. 5.2. Bennett et al. (2011) consider that the use of a Galactic mask when computing these quantities is suboptimal, and note that a fullsky computation of the 2point correlation function from the 7year WMAP ILC map lies within the 95% confidence region determined by simulations of their bestfit ΛCDM model over all angular separations. However, Copi et al. (2009) suggest that the origin of the inconsistencies between the fullsky and cutsky largescale angular correlations remains unknown, and that the observed discrepancies may indicate that the Universe is not statistically isotropic on these scales. We therefore consider the Npoint correlation functions, and related statistics, of the LklCommander map to contribute to this debate.
Fig. 40 Npoint correlation functions determined from the N_{side} = 64Planck CMB 2015 temperature estimates. Results are shown for the 2point, pseudocollapsed 3point (upper left and right panels, respectively), equilateral 3point, and connected rhombic 4point functions (lower left and right panels, respectively). The brown three dotdashed, purple dashed, and red dotdashed lines correspond to the LklCommander map analysed using the lowℓ and common masks and the Commander map analysed using the common mask, respectively. Note that the dashed and dotdashed lines lie on top of each other. The black solid line indicates the mean for 1000 MC simulations. The shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively, estimated using 1000 Commander simulations. See Sect. 4.3 for the definition of the separation angle θ. 
Probabilities to exceed the observed values of the χ^{2} statistics for the LklCommander and Commander maps at N_{side} = 64.
We compare results computed for both the LklCommander and Commander maps at N_{side} = 64 after smoothing to a FWHM of 160′. A mask is constructed for the LklCommander map by degrading the LklT_{256}93 mask to N_{side} = 64 and setting all resulting pixels with a value less than 0.5 to zero, with the remainder set to unity. The LklT_{64}92 mask retains 92% of the sky, to be compared to the 67% usable sky coverage allowed by the UT_{64}67 common mask at this resolution.
The results are presented in Fig. 40 where we compare the Npoint functions for the data and the mean values estimated from 1000 Commander simulations. The probabilities for obtaining values of the χ^{2} statistic for the Planck fiducial ΛCDM model at least as large as the observed values are provided in Table 29. For the estimation of the probabilities, we use the same set of 1000 Commander simulations for both versions of the Commander data. As noted previously, the details of the simulations for such highly smoothed data is essentially unimportant. We also provide an analysis of the LklCommander map using the common mask to enable a direct comparison with the analysis of the Commander map. In this latter case, the results for both maps are in excellent agreement. However, the LklCommander map is more consistent with simulations when the LklT_{64}92 mask is adopted for the 2point and pseudocollapsed 3point functions, but less consistent for the equilateral 3point and rhombic 4point function results. Nevertheless, the results are generally in agreement with expectations for a Gaussian, statistically isotropic model of the CMB fluctuations.
Probabilities for obtaining values of the S_{1 / 2} and statistics for the simulations at least as large as the observed values of the statistic estimated from the LklCommander and Commander maps using the LklT_{64}92 and UT_{64}67 masks, respectively. We also show the corresponding estimation of the global pvalue for the S(x) statistic.
The increased consistency of the 2point function with simulations when analysing a larger sky fraction is consistent with the observations in Copi et al. (2009). We therefore quantify this further by determining the statistical quantities introduced in Sect. 5.2 for the LklCommander map. In particular, we reassess the lack of correlation determined previously for large angular scales. It is evident from Table 30 that the results for the S_{1 / 2} and statistics are less anomalous when the lowℓ mask is applied. Moreover, the global pvalue for the S(x) statistic is substantially smaller.
We also repeat the conventional χ^{2} analysis but constraining the computations to the two separate ranges defined by θ< 60° and θ> 60°. The results of these studies are shown in Table 31. The analysis for seperation angles θ> 60° indicates that the unusually good fit of the observed 2point function to the mean 2point function determined for the ΛCDM model is independent of the mask used in the analysis. Conversely, the results for the angles θ< 60° indicate a strong dependence on the mask. It appears that the decreased significance of the χ^{2} statistic for the 2point function of the LklCommander map reported in Table 29 is related mainly to correlations in the data for separation angles smaller than 60°.
Probabilities for obtaining values of the χ^{2} statistic for the simulations at least as large as the observed values of the statistic estimated from the LklCommander and Commander maps using the LklT_{64}92 and UT_{64}67 masks, respectively.
Our results do appear to indicate that computations made on larger sky fractions increase the consistency of the 2point function with simulations. We therefore also test how the hemispherical asymmetry observed previously is affected. The results for the ecliptic frame are presented in Table 32. We find that the asymmetry is larger for the LklCommander map than for the Commander map in the case of the 2point function, but does not change substantially for the 3point and 4point functions.
Probabilities for obtaining values of the χ^{2} statistic and ratio of χ^{2} of the Npoint functions for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic on the northern and southern ecliptic hemispheres estimated from the LklCommander and Commander maps using the LklT_{64}92 and UT_{64}67 masks, respectively.
7.3. Dipole modulation and directionality
7.3.1. Variance asymmetry
pvalues for the variance asymmetry measured by different discs from the Planck 2015 LklCommander and Commander temperature solutions using the LklT_{256}93 and UT_{256}73 masks, respectively.
Localvariance dipole directions measured by 8° discs for the Planck 2015 LklCommander and Commander temperature solutions.
Here we apply the localvariance analysis of Sect. 6.1 to the LklCommander map and compare the results with those of the Commander map. Contrary to the analysis of Sect. 6.1, where fullresolution (N_{side} = 2048) maps were used, here the Commander map is downgraded to N_{side} = 256 in order to consistently compare the results for both maps. The simulations used for estimating the significance levels are also downgraded to the same resolution, and convolved with the corresponding beam function. Otherwise, the procedure is identical to the one described in Sect. 6.1, e.g., the same number of discs has been used to construct the localvariance maps. Here we only present the results when no highpass filtering has been applied to the maps; this is to avoid confusion as our objective in this section is only to compare the general properties of the LklCommander map to those of the standard componentseparated maps.
Table 33 summarizes the significance levels measured by our variance asymmetry analysis using discs of different radii, for the Planck 2015 Commander and LklCommander temperature maps. The pvalues represent the fraction of simulations with localvariance dipole amplitudes larger than those inferred from the data. We in addition present in Table 34 the preferred variance asymmetry directions for both maps using 8° discs.
Our results show consistency between the two maps. The small change in the preferred direction is expected from the change in the mask, and agrees specifically with the directions found by the analysis of the QML dipole modulation analysis in Sect. 7.3.3. One interesting observation is that the large variance asymmetry significance is now extended to cases where larger discs are used. Note that no highpass filtering has been applied in the present analysis, and therefore pvalues inferred from the Commander map increase with the disc size. As explained in Sect. 6.1, the low observed significance levels for larger discs is due to the cosmic variance associated with the largestscale modes. The observed increase in the significance levels for the LklCommander map is therefore interestingly consistent with this picture; the mask in this case is smaller and therefore a larger fraction of the sky is available. This in turn provides more data on the largest scales, and therefore lowers the impact of the cosmic variance.
7.3.2. Dipole modulation: pixelbased likelihood
Table 35 presents constraints on the dipole modulation model as derived from the LklCommander map and the LklT_{32}93 mask that includes 93% of the sky, updating the results from Sect. 6.2 for the Commander map. We find that all previously reported results are robust with respect to data selection and sky coverage. In particular, the bestfit dipole modulation amplitude at 5° FWHM is 5.9% in the LklCommander map, and is thus stable to within about 0.3σ when increasing the sky fraction from 78% to 93%. Likewise, the marginal lowℓ power spectrum amplitude, q, shifts upward by 0.4σ, and the power spectrum tilt, n, downward by 0.3σ, for the same sky fraction increases.
To assess the statistical significance of these shifts, we compare with Gaussian statistics, creating two Gaussian random vectors with 78 and 93 elements, respectively, where the first 78 elements of the latter vector are identical to the first vector. From these, we compute the difference between the two means, after normalizing each so that their individual errors in the mean are unity. Repeating this simple calculation 10^{5} times, we find that 48% of all Gaussian realizations observe shifts larger than 0.3σ, and 34% observe shifts larger than 0.4σ. Thus, the parameter differences due to the different data selection and sky fractions reported above are consistent with expectations from random Gaussian statistics.
7.3.3. Dipole modulation: QML analysis
We also repeat the QML dipole modulation analysis of Sect. 6.3 for the LklCommander map and corresponding mask. Table 36 summarizes the results of the lowℓ dipole modulation for the LklCommander temperature solution, compared with the Commander map.
The bestfit modulation amplitude for LklCommander is 5.8% and the small 0.5% shift from the Commander bestfit amplitude corresponds to a decrease of approximately 0.4σ. These results mirror very closely the results found above for the pixelbased likelihood approach to dipole modulation, as expected, and the observed shifts are perfectly consistent with those expected from the change in the mask.
7.3.4. Bipolar spherical harmonics
We next perform a dipole modulation analysis on the LklCommander temperature map using the BipoSH formalism from Sect. 6.4. The dipole modulation amplitude inferred from the analysis is smaller that that deduced from analysing the Commander map as seen in Table 37. However, it should be noted that the probability for simulations to yield a dipole modulation amplitude equal to or greater than the amplitude inferred from data is 0.4%, which is smaller by a factor of approximately 2.4 as compared to the pvalue inferred from analysis on Commander. The reduction in the dipole amplitude and the enhanced significance can both be attributed to the reduced power bias which is a result of the increased sky coverage.
Amplitude (A) and direction of the dipole modulation in Galactic coordinates as estimated for the multipole range ℓ ∈ [2,64] using the BipoSH analysis on LklCommander and Commander maps. The former results were derived using the LklT_{256}93 mask; the latter are those determined previously in Sect. 6.4.
7.4. Summary
Using a larger sky fraction in our analyses leads to small changes in the results related to largeangularscale anomalies, but these are essentially consistent with expectations from random Gaussian statistics. In particular, the asymmetry in power on the sky, as parameterized by a dipole modulation model, is robust to mask changes.
8. Polarization analysis
As previously discussed in Sect. 2, large angularscale CMB fluctuations in the Planck polarization data have been suppressed by a postprocessing highpass filter to minimize the impact of systematic artefacts. Therefore, no polarization results concerning CMB statistical anomalies on such scales are presented in this paper. In addition, a noise mismatch between simulations and data also limits our ability to study polarization more generally. Nevertheless, a local analysis of the polarization data for stacked patches of the sky can still be performed, in order to test the statistical properties of the CMB anisotropies. In this case, the stacking procedure mitigates the impact of the smallscale noise and potential systematic effects.
Traditionally, the Stokes parameters Q and U are used to describe the CMB polarization anisotropies (e.g., Zaldarriaga & Seljak 1997). Such quantities are not rotationally invariant, thus for the stacking analysis it is convenient to consider a local rotation of the Stokes parameters, resulting in quantities denoted by Q_{r} and U_{r}, as described in Sect. 8.1. Additionally, several other related quantities can be defined.
The polarization amplitude and polarization angle , are commonly used quantities in, for example, Galactic astrophysics. However, unbiased estimators of these quantities in the presence of anisotropic and/or correlated noise are hard to define (Plaszczynski et al. 2014). Of course, a direct comparison of the observed (noisebiased) quantity to simulations analysed in the same manner is possible, but we elect here to defer the study of this representation of the polarization signal, using maps of the polarization amplitude only to define peaks around which stacking can be applied.
The rotationally invariant quantities referred to as E and B modes are commonly used for the global analysis of CMB data. Although the Emode maps are not analysed in detail here, they are considered qualitatively, so that it is appropriate to recall their construction. Since the quantities Q ± iU, defined relative to the direction vectors , transform as spin2 variables under rotations around the axis, they can be expanded as (62)where are the spinweighted spherical harmonics and are the corresponding harmonic coefficients. If we define then the invariant quantities are given by
8.1. Stacking around temperature hot and cold spots
The stacking of CMB anisotropies around peaks (hot and cold spots) on the sky yields characteristic temperature and polarization patterns that contain valuable information about the physics of recombination (Komatsu et al. 2011). Statistical analysis of stacked images differs from the other tests in this paper in several respects. First, peakrelated new physics may be revealed that is difficult to find in a global analysis, for example, the nonGaussian CMB cold spots predicted by a modulated preheating model (Bond et al. 2009). Secondly, stacking is a local operation, which naturally avoids maskinduced complications. Thus stacking can be used as a transparent and intuitive method to test the robustness of anomalies found with other methods. Alternatively, it can be applied as a quality indicator of the data at the map level.
Our stacking procedure is as follows. Hot (or cold) peaks are selected in the temperature map as local extrema with negative (or positive) second derivatives, and classified relative to a given threshold ν (in rms units of the temperature map). Since the spinorial components Q and U are expressed in a local coordinate system, we employ a configuration in which the Stokes parameters around a peak at the direction can be superposed (Kamionkowski et al. 1997). In particular, we use a locally defined rotation of the Stokes parameters that is written as: where φ is the angle between the axis aligned along a meridian (pointing to the south by convention) in the local coordinate system centred on a peak at and the great circle connecting this peak to a position . This definition decomposes the linear polarization into radial (Q_{r}> 0) and tangential (Q_{r}< 0) contributions around the peaks. This definition of Q_{r} is equivalent to the “tangential shear” used in weak lensing studies.
For visualization purposes, a flat patch around each peak is then extracted, and the average stacked image computed from the subset. A position on the sky at an angular distance θ from the central peak is labelled with the flatsky coordinates (69)Here ϖ = 2sin(θ/ 2) ≈ θ is the effective flatsky radius. For the angular scales of a few degrees considered in the stacking analyses the difference between ϖ and θ is negligible. We use ϖ for analyses in the flatsky approximation, and θ for analyses directly on the sphere.
The stacking process tends to provide an image with azimuthal symmetry about its centre, due to the almost uncorrelated orientations of the temperature peaks. The stacked images of temperature patches around hot spots selected above the null threshold for both the Commander data and a corresponding simulation are shown in the top row of Fig. 41. The observed patterns are in excellent agreement. Stacking around cold spots yields similar patterns but with flipped sign. Given the symmetry, it is often useful to consider the radial profile obtained by averaging the stacked image over the azimuthal angle φ. Figure 42 shows such a profile determined from the stacked temperature image.
Fig. 41 From top to bottom, T, Q, U, Q_{r}, and U_{r} stacked images (in μK units) extracted around temperature hot spots selected above the null threshold (ν = 0) in the Commander sky map for data (left column) and an equivalent simulation (right column). The horizontal and vertical axes of the flatsky projection are labelled in degrees. 
At this point, it is useful to consider the underlying physics represented by the various patterns in the stacked images. During recombination, the sound horizon extends an angle θ_{s} = r_{s}/D_{A} ≈ 0.011 (0.61°), where r_{s} ≈ 0.15 Gpc is the size of the sound horizon at recombination and D_{A} ≈ 14 Gpc is the angulardiameter distance to the last scattering surface. To understand the ring patterns in the stacking image, projection effects must be taken into account. Firstly, all 3D modes with wavenumber k ≥ ℓ/D_{A} contribute to a 2D ℓmode. More modes contribute to, and therefore enhance, the power at lower ℓ. For the first acoustic peak, the net effect is a π/ 4 phase shift towards lower ℓ, such that ℓ_{s} ≈ (π−π/ 4) /θ_{s} ≈ 220. The projected acoustic scale on the temperature map is of order (0.81°). Secondly, the stacked 2D modes around peaks interfere with each other. The first dark ring appears at (1.0°). The factor 1.22 is the ratio of the first minimum of the projection kernel, the Bessel function J_{0}, to the first minimum of the unprojected cosine wave.
Fig. 42 Radial profile μ_{T}(ϖ) derived from the stacked temperature image (see Fig. 41 or 45). The denominators σ_{0} and σ_{2} are the theoretical rms values of CMB T and ∇^{2}T, respectively. The theoretical ⟨ μ_{T}(ϖ) ⟩ is a linear combination of ⟨ T(ϖ)(T(0) /σ_{0}) ⟩ (green dashdotted line) and ⟨ T(ϖ)(−∇^{2}T(0)) /σ_{2}) ⟩ (blue dotted line). For all four componentseparated maps, the deviation of μ_{T} from the ensemble mean ⟨ μ_{T} ⟩ of the fiducial model (here the Planck 2015 ΛCDM best fit) is consistent with cosmic variance, and can be related to the lowℓ power deficit. The example powerdeficit ⟨ μ_{T} ⟩ (purple dashed line) is the theoretical prediction of ⟨ μ_{T} ⟩ if the fiducial model C_{ℓ}s are reduced by 10% in the range 2 ≤ ℓ ≤ 50. 
The dark ring can also be regarded as a consequence of the correlation between T and −∇^{2}T. At the temperature maxima −∇^{2}T is positive, with an amplitude of order . Thus, the quadratic terms in the local expansion of the temperature field have a negative contribution that grows as . At the quadratic terms dominate and the T(− ∇^{2}T) correlation becomes negative. Meanwhile, the T(− ∇^{2}T) correlation tends to zero on the scale , where the temperature autocorrelation becomes weak and the local quadratic expansion starts to fail. As shown in Fig. 42, the dark ring appears at the critical point where the T(− ∇^{2}T) correlation reaches its minimum and turns back toward zero.
We have discussed the projection effects that make the projected radial acoustic scale on a stacked T image larger than θ_{s}. For Q_{r}, the most striking patterns in the image have more intuitive simple explanations, since the stacking is essentially the realspace equivalent of the temperature polarization correlation. The projection function contains an extra ℓ^{2} factor, which enhances the highℓ power and reduces the projected radial acoustic scale, coincidentally, back to ≈ θ_{s}. The quadrupole responsible for the polarization around peaks is induced by gravity on angular scales larger than twice the size of the horizon at decoupling. In the case of an overdensity, this causes a flow of photons towards the gravitational well on these scales, inducing a quadrupolar pattern (see, e.g., Coulson et al. 1994). The spherical symmetry of the gravitational interaction causes a rotation of the quadrupole in the vicinity of the well, resulting in a radial configuration in polarization. This radial polarization pattern implies Q_{r}> 0 and an overdensity implies T< 0 by the SachsWolfe formulae, which leads to anticorrelation on these scales. Similarily, an underdensity leads to an outward flow and induces a tangential polarization pattern, once again leading to anticorrelation on these scales. At smaller scales, the polarized contribution is dominated by the dynamics of the photon fluid. The acoustic oscillations modulate the polarization pattern, leading to the different rings in the stacked images. The most noticeable rings in the stacked Q_{r} image are approximately at θ_{s} and 2θ_{s}. Thanks to the ℓ^{2} enhancement, multiple acoustic peaks in the TE power spectrum may be captured and projected into ring patterns in the stacked polarization images. As photons flow towards the overdensity, they are compressed and the temperature increases, slowing the fluid descent into the well. Eventually, the radiation pressure becomes large enough to reverse the photon flow. This expansion cools the photons until they fall back towards the well. Note that the resulting inner ring was not observed in the WMAP analysis (Komatsu et al. 2011), since the resolution was too low.
Fig. 43 Mean radial profiles of T, Q_{r}, and U_{r} in μK obtained for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). Each individual panel contains (top) the mean radial profiles and (bottom) the differences (denoted “Diff”) between the mean profiles of the data and those computed from the ensemble mean of the simulations. Results based on stacks around temperature hot and cold spots are shown in the left and right columns, respectively. Upper plots present results for peaks selected above the null threshold, while lower plots show the equivalent results for peak amplitudes above (hot spots) or below (cold spots) 3 times the dispersion of the temperature map. The black dots (connected by dashed lines) depict the mean profiles and the shaded regions correspond to the 1σ (68%) and 2σ (95%) error bars. The mean profiles and error bars are determined from SEVEM simulations. Note that the Diff curves for each componentseparation method are computed using the corresponding ensemble average, although only the ensemble average from SEVEM is shown here. 
Figure 41 clearly reveals all of the features described above. The two bright rings at θ_{s} ≈ 0.011 (0.6°) and 2θ_{s} ≈ 0.021 (1.2°) are the predicted patterns associated with the first acoustic peak at ℓ ≈ 310, while the two faint rings are a striking illustration of the detection of multiple acoustic peaks in the TE power spectrum. The largescale anticorrelation is suppressed due to the scaledependent bias which results from the fact that peaks are defined by the second derivatives of the temperature field (e.g., Desjacques 2008).
We are now in a position to discuss the consistency of the Planck results with the predictions of a ΛCDM cosmology. For simplicity, further analysis is focused on the angular profiles, and specifically the mean, μ(θ), estimated as the average of the angular profiles around all hot (cold) peaks above (below) a certain threshold ν. This analysis is performed directly on the sphere to avoid any repixelization error. Note that the expected value of the mean temperature angular profile is proportional to , whilst the expected values of the Q_{r} and U_{r} mean angular profiles are approximately proportional to and , respectively. Since T has even parity and B has odd parity, the expectation value for is zero, and the U_{r} mean angular profile is therefore expected to vanish.
A χ^{2} estimator is used to quantify the differences between the profiles obtained from the data and the expected values estimated with simulations: (70)with the covariance matrix defined as (71)where the sum is over the N simulations used to estimate this matrix and is the ensemble average. Note that although the profiles in Fig. 41 are derived from data at a resolution N_{side} = 1024, faster convergence of the χ^{2} statistic is achieved using maps at a lower resolution. We have verified that the results remain unchanged when adopting data with N_{side} = 512.
Figure 43 presents a comparison between the profiles obtained from the componentseparated data and the mean value estimated from simulations processed through the SEVEM pipeline. Note that the error bars for the temperature profiles are asymmetric due to a bias in the selection of the peaks above a given threshold. Results for hot and cold spots are shown for two different thresholds, ν = 0 and ν = 3. There is generally excellent agreement between the different componentseparation methods. A systematic deviation between the data and the simulations for the hot peaks in temperature (ν = 0) is seen at a level greater than 1σ. This discrepancy increases at higher thresholds, reaching values of about 2σ for the ν = 3 case. Similar behaviour is seen for the cold spots. For the Q_{r} angular profiles, the most striking differences appear around θ = 2° in the ν = 3 case for hot peaks, and around for the cold peaks. For the U_{r} angular profiles, where a null signal is expected (i.e., only noise is expected to be present), deviations at similar levels are seen.
Fig. 44 χ^{2} distributions obtained from the T (left column), Q_{r} (middle column), and U_{r} (right column) mean radial profiles centred on temperature hot spots selected above the null threshold (upper row) and three times the dispersion of the map (bottom row). The black lines correspond to the theoretical χ^{2} distribution with 12 degrees of freedom, whilst the histograms show the distributions determined from 100 simulations computed through the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) pipelines. The vertical lines represent the χ^{2} values obtained from the data. 
pvalues of the T, Q_{r}, and U_{r} angular profiles computed from the stacking of hot and cold spots selected above the ν = 0 and ν = 3 thresholds.
Table 38 presents the corresponding pvalues for this comparison. A theoretical χ^{2} distribution is used to determine the probability that a sky drawn from the ΛCDM cosmology has a value larger than that derived from the data. We have verified this approach by comparing the empirical χ^{2} distribution estimated from 100 simulations (in which the mean value and the covariance matrix are computed from a further 900 simulations) with the theoretical distribution with the corresponding degrees of freedom (see Fig. 44). The χ^{2} value of the data is then estimated using the mean value and the covariance matrix determined from simulations. Although some differences are found among the componentseparation methods, a general consistency between model and data is found.
Although the χ^{2} test has the advantage of being sensitive to different types of deviations between model and data, does not assume prior knowledge about possible departures from the model, and can account for correlations between the various tests from which it is constructed, it can nevertheless be suboptimal under certain conditions. This appears to be the case when considering the systematic shift between data and simulations seen in the temperature profiles μ_{T} – the χ^{2} statistic is not particularly sensitive to systematic deviations of constant sign. We therefore consider an alternative quantity, the integrated profile deviation, defined as (72)where R, the size of stacking patches, is taken to be 3.̊5 in this case. The weighting function is chosen to be proportional to the expected profile, but the results are robust for other choices, e.g., W = 1. The pvalues obtained in this case are given in Table 39. These are consistent with what might be expected from visual inspection of the plots, i.e., the deviations are typically close to 2σ. These deviations are likely to be connected to the deficit in the observed power spectrum at low multipoles, as may be seen in Fig. 42. Here, the purple dashed line indicates the reduction in if the theoretical C_{ℓ} values are reduced by 10% over the range 2 ≤ ℓ ≤ 50.
8.2. Generalized stacking
In this section, a much wider class of stacking methods is introduced, with particular emphasis on oriented stacking, a novel approach that has not previoulsy been explored in the literature. We regard the stacking as oriented if the orientation of the local coordinate frame, and in particular the φ = 0 axis, is correlated with the map that is being stacked. Thus, the stacking methodology in Sect. 8.1 is considered unoriented, because the orientation is defined relative to the local meridian pointing towards the Galactic south, rather than any property of the data themselves. Alternative approaches to unoriented stacking can also be considered. In this subsection, the orientation of each patch is chosen randomly from a uniform distribution in [0,2π). The unoriented T and Q_{r} images can then be directly compared with previous sections.
For unoriented stacking, the ensemble average of stacked fields cannot result in any intrinsic φdependence, as this would be averaged out by the uncorrelated orientation choices. The φdependence due to a specific choice of representation can always be removed via a local rotation. For example, the ensemble averages of Q + iU around unoriented temperature peaks are proportional to e^{2iφ}. A local rotation (Q,U) → (Q_{r},U_{r}) (Kamionkowski et al. 1997) removes the e^{2iφ} factor and compresses the information into a single real map Q_{r}. For oriented stacking, the φdependence can be a mixture of a few Fourier modes (e^{imφ}, for integer m). Each m mode corresponds to a radial (ϖdependent) function.
pvalues of Δμ_{T} computed from the stacking of hot and cold spots selected above the ν = 0 and ν = 3 thresholds from the Commander, NILC, SEVEM, and SMICA maps.
In what follows, we use the N_{side} = 1024 componentseparated maps at a resolution of 10′ FWHM. The use of this higher resolution as compared to the N_{side} = 512 data used in Sect. 8.1 is motivated by the smallerscale features that are expected to result from the oriented stacking.
We also introduce the concept of the noisefree ensemble average (NFEA), which is defined as the ensemble average of stacked CMBonly maps for a fiducial cosmology. Recall that the fiducial model for the simulated sky maps, the Planck 2013 bestfit ΛCDM model (Planck Collaboration XVI 2014), differs from the updated Planck 2015 bestfit ΛCDM model (Planck Collaboration XIII 2016). In previous sections, this mismatch was partially accommodated by rescaling the CMB signal by a fixed scale factor. Here, we instead specifically adopt the 2015 best fit as a fiducial model for the data. When comparing the data to the simulations, we subtract the corresponding NFEA to minimize any bias resulting from cosmology dependence.
In the context of random Gaussian fields, the NFEA can be computed straightforwardly following Bardeen et al. (1986): (73)where M is the map (around the central peak) to be stacked, and w is the collection of Gaussian variables (on the central peak) that are related to peak selection and orientation determination. Equation (73) is only valid for Gaussian random variables. If the patch is rotated before stacking, the field value evaluated at a dynamic coordinate is, in general, not a random Gaussian variable. However, statistical isotropy guarantees that the rotation of patches is equivalent to an orientation constraint on the nonzerospin field. For example, orienting each patch in the direction where U = 0 and Q> 0 is equivalent to the unoriented stacking case where only peaks satisfying the additional constraint −ϵ/ 2 < arg(Q + iU) <ϵ/ 2 (ϵ → 0^{+}) are selected.
A further source of statistical bias can arise from noise mismatch between the simulations and the data. Since the effect of noise is to introduce random shifts in the peaks and hence suppress patterns in the stacked images, any noise mismatch can lead to pattern mismatch between the data and simulations. For the temperature data, the contribution due to noise mismatch is estimated to be at the subpercent level, lower than the cosmic variance. For stacking on polarization peaks, the impact of the noise mismatch cannot be safely ignored. Thus, for quantitative comparisons in this paper, we only consider stacking on temperature peaks.
8.2.1. Oriented temperature stacking
The most straightforward way to orient a patch centred on a temperature peak is to align the horizontal axis with the major axis defined by a local quadratic expansion of the temperature field around the peak. The disadvantage of doing so is that the orientation is dominated by smallscale fluctuations that are noisesensitive. A better choice is to use the major axis of the inverse Laplacian ∇^{2}T that filters out the smallscale power. The inverse Laplacian is defined as: (74)where are the harmonic coefficients of the masked temperature map. Spin2 maps Q_{T}, U_{T} are then defined by: (75)In the flatsky limit, and U_{T} ≈ −2∂_{x}∂_{y}∇^{2}T. To align the ∇^{2}T axes of the patches, we rotate each patch so that U_{T} vanishes and Q_{T} ≥ 0 for the central peak.
Figure 45 presents the stacked images of SMICA temperature patches centred on temperature hot spots selected above the threshold ν = 0, in both unoriented and oriented forms. These are seen to be in excellent agreement with their accompanying NFEAs, and, in the case of the unoriented stacks, with the results shown in Fig. 41, despite the different stacking methodologies adopted (and component separation method selected for visualization purposes).
Fig. 45 Comparison between unoriented stacking (upper panels) and oriented stacking (lower panels) of temperature patches around temperature hot spots selected above the null threshold (ν = 0). The left panels are the stacked SMICA maps, and the right panels their corresponding NFEAs. The image units are μK. 
The oriented T image is notably different from the unoriented one. The alignment between the major axis (of ∇^{2}T) and the horizontal axis results in an ellipse elongated along the horizontal axis, rather than a central disc. Moreover, the quadraticterm contribution is suppressed along the horizontal axis where the temperature profile is smoother, and enhanced along the vertical axis where the temperature profile is sharper. As a consequence, the dark ring visible in the upper panel at 1° splits into two cold blobs along the vertical axis.
To proceed with a quantitative analysis, we extract Fourier modes T_{m}(ϖ) from the stacked map T_{stack}(ϖ,φ) as follows: (76)where δ_{m0} is the Kronecker delta function. For odd m, the NFEA ⟨T_{m}⟩ vanishes due to statistical isotropy. For even m, a straightforward calculation shows that only T_{0}(ϖ), which is equivalent to μ_{T}(ϖ), and T_{2}(ϖ) have nonzero NFEAs.
As discussed previously in Sect. 8.1, there are some shortcomings of the standard χ^{2} procedure that is generally used to assess the consistency of the data with simulations. The problem is simplified by studying the statistics of an integrated profile deviation: (77)where R, the size of the stacking patches, is taken to be 2° in our examples. The purpose of removing the NFEA, ⟨T_{m}(ϖ)⟩, which differs for the data and the simulations, is to minimize the impact of the cosmology dependence. A natural choice for the filter is ⟨T_{m}(ϖ)⟩ itself with a proper normalization: (78)For the filter given by Eq. (78), the integrated profile deviation describes the relative deviation from the NFEA. If ΛCDM is the correct model, the deviation is due to cosmic variance and noise. The distribution of is obtained from simulations.
Table 40 presents a comparison of the values derived from the Planck data and the FFP8 simulations. No inconsistencies in excess of the 3σ level have been found, although tensions around 2σ are seen.
The m = 0 projection kernel J_{0} [(ℓ + 1 / 2)ϖ] peaks at low ℓ. Thus is cosmicvariance sensitive and the apparent discrepancy in it could be related to a lowℓ power deficit. An example is shown in Fig. 42 for illustration. To test the robustness of this result, we have tried three additional filters: a tophat filter W = 1, a linear filter W = ϖ, and a Gaussian filter with σ_{g} = 1°. In all cases, the power deficit remains at about the 2σ level.
Although the deficit is not significant enough to falsify the ΛCDM model, further investigation of its properties may still be interesting and help us understand the other anomalies discussed in this paper. We consider two possibilities. Firstly the amplitude of the deficit is of order 5–10%, which coincides with the level of hemispherical power asymmetry discussed in Sect. 6.1. To test whether the deficit is localized on one hemisphere, we define the “north” direction to be aligned with the power asymmetry direction at (l,b) = (212°,−13°) (Akrami et al. 2014) and compute on the northern and southern hemispheres separately. The results are presented in Table 41. Although the deficit is more significant for the southern hemisphere, it remains consistent with the ΛCDM prediction. Secondly, it is of interest to determine whether the deficit is related to the Cold Spot discussed in Sect. 5.7. We therefore mask out the Cold Spot using a disc of radius 6° and repeat the calculation. The impact of this region on the deficit is insignificant.
Tensions at the 2σ level are also seen for . However, due to the additional ℓ^{2} factor in the projection kernel, the oriented (m = 2) component is more sensitive to highℓ power where the cosmic variance is small, and an understanding of the noise properties of the data is more important. The former implies that the related uncertainty in is, in general, smaller than that in . However, a mismatched cosmology, perhaps arising from a different primordial power amplitude A_{s}, can then lead to significant tension between the data and the simulations. Indeed, we find that without application of our cosmology calibration (i.e., the subtraction of the NFEA in Eq. (77)) the tension between the data and simulations increases by about 0.5σ, whereas the variation of the tension is ≲ 0.2σ. The highℓ sensitivity of also requires the use of an accurate noise model, and it is possible that the 1–2σ tension in may be alleviated once improved noise simulations are available.
8.2.2. Oriented polarization stacking
The stacked Q and U images can be decomposed into Fourier modes, . For unoriented Q + iU stacking on temperature peaks, only P_{2}(ϖ) has a nonzero NFEA, and it can be linked to the conventional Q_{r} stacking via P_{2} = −Q_{r}. Figure 46 shows that the stacked Q_{r} image is in excellent agreement with its NFEA and the corresponding stacked image (fourth panel) in Fig. 41, despite the different stacking methodologies adopted (and componentseparation method selected for visualization purposes). The length and orientation of the headless vectors represent the polarization amplitude, , and direction.
Fig. 46 Stacked Q_{r} image around temperature hot spots selected above the null threshold (ν = 0) in the SMICA sky map. The left panel corresponds to the observed data and the right panel shows the NFEA. The image units are μK. The headless vectors (black solid lines) are the polarization directions for stacked Q_{stack}, U_{stack}. The lengths of the headless vectors are proportional to the polarization amplitude P_{stack}. 
We next consider oriented stacking of the polarization maps, again using Q_{T}, U_{T} to define the orientation of the patches. The stacked polarization images around temperature peaks have m = 0,2,4 Fourier components. We can also choose to stack the polarization maps on P_{T} peaks, where . This picks up m = 0,4 Fourier modes with no circularly symmetric (Q_{r}, m = 2) mode. In Fig. 47 we compare the (Q,U) images stacked centred either on T peaks (top panel) or on P_{T} peaks (bottom panel) with their corresponding NFEAs, and find excellent agreement.
Fig. 47 Oriented stacking of polarization fields (Q, U) on temperature maxima (upper panels) and P_{T} maxima (lower panels). In both cases the threshold ν = 0 is used and the orientation is chosen such that U_{T} = 0 and Q_{T} ≥ 0 on the central peak. The image units are μK. The left panels are the stacked SMICA maps, and the right panels their NFEAs. See Fig. 46 for the meaning of the headless vectors (black dashed lines). 
For a quantitative comparison, we only consider stacking on temperature peaks and define the polarization integrated profile deviation (79)where by default the filter is (80)The comparison of (m = 0,2,4) between the data and the simulations is shown in Table 42, where the results are seen to be in excellent agreement.
Finally, we note that the peak selection does not have to be made from the temperature map. In Fig. 48 we show a few examples of stacking on polarization peaks using the N_{side} = 512 maps. The higherresolution polarization data are too noisy for peak selection. In the upper panels, we compare stacked images of the Emode map centred around Emode peaks with the corresponding NFEA. We find that the noise impact is relatively minor for 20′ FWHM maps and the plots are in qualitatively good agreement. Another possibility, shown in the lower panels, is to stack polarization maps centred on peaks determined from the corresponding polarization amplitude map. In this case the peaks are strongly biased by the quadratic noise contribution and quite visible deviation from the NFEA is observed in the stacked image.
Fig. 48 Top: Emode maps stacked on the unoriented Emode maxima computed above the null threshold ν = 0. Bottom: Q stacked around oriented polarization amplitude (P) maxima. In this case, no threshold is used and the orientation is chosen such that U = 0 and Q ≥ 0 on the central peak. The left panels are the stacked SMICA maps, and the right panels their corresponding NFEAs. See Fig. 46 for the meaning of the headless vectors (black dashed lines). The image units are μK. 
9. Conclusions
In this paper, we have presented a study of the statistical isotropy and Gaussianity of the CMB using the Planck 2015 data, including the full mission for temperature. We do not claim that our results support or refute any particular physical model. Rather, we focus on nullhypothesis testing: a number of tests are performed, then pvalues are calculated and reported. It is in the very nature of such a modelindependent approach to leave the detailed interpretation to the reader. However, we do address the important subject of a posteriori correction where possible.
The statistical tests are performed on maps of the CMB anisotropy that result from the application of the four componentseparation methods described in Planck Collaboration IX (2016). All of the results presented here are robust with respect to the choice of componentseparated CMB map. This is important since it demonstrates the high quality and equivalence of the Planck componentseparated data products rendered by different methodologies under varying assumptions.
We find that the CMB is largely consistent with statistical isotropy, although there are a few indications of anomalies with respect to the expectations of ΛCDM. Some of the tests we have performed are the same as those in PCIS13, in which case the results are consistent. Since many of these anomalies were also observed in the WMAP temperature data, we reemphasize explicitly the statement we made in 2013 – that the agreement between the two independent experiments effectively rules out the possibility that the origin of these features can be found in residual systematic artefacts present in either data set. We have also performed a number of new tests, in order to try to narrow down the nature of the apparent violations of statistical isotropy. In addition, although the componentseparated polarization maps contained in the Planck 2015 release are highpass filtered, we have performed a stacking analysis that tests some aspects of the polarized sky while mitigating the impacts of noise and systematic effects.
In Sect. 4, we examined aspects of the Gaussianity of the CMB fluctuations. Tests of skewness, kurtosis, multinormality, Npoint functions, and Minkowski functionals yielded no indications of significant departures from Gaussianity, while the variance of the CMB map was found to be low, in agreement with previous studies (PCIS13). Firstorder moments of filtered maps also exhibit the lowvariance anomaly, as well as a kurtosis excess on certain scales associated with the Cold Spot. A new study of peak statistics finds results consistent with the expectations for a Gaussian random field, although the Cold Spot is again detected.
Section 5 provides an updated study of several previously known peculiarities. We study in detail the low variance anomaly, which appears to be associated with the known lowℓ deficit in the angular power spectrum. We confirm the lack of largescale angular correlations, relatively featureless northern ecliptic hemisphere 3 and 4point functions, and indications of violations of point and mirrorparity symmetry, although we make little or no attempt to correct these for a posteriori effects. We place tight constraints on a quadrupolar power modulation. The Cold Spot is examined further, and, while we find variance, skewness, and kurtosis angular profiles consistent with the expectations of statistically isotropic simulations, the mean temperature profile is anomalous at roughly the 1% level, apparently due to the surrounding hot ring – the feature that deviates most from the Gaussian model.
In Sect. 6 we perform a series of tests probing the wellknown largescale dipolar power asymmetry. We detect the asymmetry via pixeltopixel variance, as well as by measuring power explicitly or indirectly via ℓ to ℓ ± 1 mode coupling. The latter approach lends itself to a posteriori correction, which reduces the significance of the asymmetry substantially when no model for the anomaly is assumed. In addition, we perform two independent but related tests of directionality. One finds suggestions of anomalous clustering of directions out to relatively small scales while the other does not, evidently due to being optimized for slightly different forms of directionality.
Section 7 demonstrates that the significances of several largeangularscale anomalies are robust to the use of larger sky coverage, with the observed small changes being consistent with expectations from random Gaussian statistics.
Finally, Sect. 8 presents the results of the stacking of temperature and polarization peaks. We find results that are largely consistent with statistically isotropic simulations, both for oriented and unoriented stacking. The exception is a low unoriented temperature profile, which seems to be yet another reflection of the largescale power deficit.
With the Planck 2015 release, we are probably near the limit of our ability to probe the CMB anomalies with temperature fluctuations alone. The use of largeangularscale polarization, expected for the final Planck release, should enable independent tests of these peculiar features. Importantly, this will reduce or eliminate the subjectivity and ambiguity in interpreting their statistical significance. It is a tantalizing possibility that some of the anomalies described in this paper will take us beyond the standard model of cosmology.
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 and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).
The digits 8 and 4 denote the order of the spherical SavitzkyGolay kernel and the smoothing weight, described in Appendix A.
The Galactic mask induces a preferred direction in the analysis of the MC simulation ensemble, which affects the significance of the results determined from the data. See BenDavid & Kovetz (2014) for a discussion.
Actually only SEVEM and SMICA achieve their minimum at ℓ_{max} = 67, whereas NILC and Commander achieve theirs at ℓ_{max} = 14 and 240, respectively. Such scatter is expected when searching over a large number of possible ℓ ranges. The reconstructed amplitudes for each componentseparation method are well within the error budgets of the estimator.
The BipoSH spectra, as defined in Eq. (52), restrict us to working with only evenparity BipoSH coefficients (L + d is even) due to the vanishing of otherwise. While most known isotropyviolating phenomena like weak lensing, Doppler boost, noncircular beams, etc., can only produce evenparity BipoSH spectra, measurement of oddparity BipoSH spectra can be used to test for systematic effects, or to search for the signatures of exotic effects such as the lensing of CMB photons by tensor metric perturbations.
Departing from the analysis in PCIS13, we do not use an apodized version of the common mask. Simulations indicate that the error on the power spectrum for those multipoles in the range 300 to 500 where the significance is highest is up to 20% larger in this case, with the corresponding error on preferred direction being typically 8% larger.
Note that here we have not specified what δC_{ℓℓ + 1} is (it is fully specified by choosing a parameter X to modulate). This is because we have decided to weight each ℓ equally and thus any strictly positive choice for δC_{ℓℓ + 1} will be equivalent, since in step 3 we force the mean length of the vectors at each ℓ to be equal.
Note that the SEVEM maps used in this section have been inpainted within 3% of the sky towards the Galactic centre using a simple diffusive inpainting algorithm. This prevents residual foreground contamination from propagating to neighbouring regions when downgrading the map. The other componentseparated maps are not preprocessed in this way since some form of inpainting of the most contaminated regions was already implemented as part of the component separation algorithms.
Acknowledgments
The Planck Collaboration acknowledges the support of: ESA; CNES and CNRS/INSUIN2P3INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, 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); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planckcollaboration. Some of the results in this paper have been derived using the HEALPix package.
References
 Ackerman, L., Carroll, S. M., & Wise, M. B. 2007, Phys. Rev. D, 75, 083502 [NASA ADS] [CrossRef] [Google Scholar]
 Adhikari, S. 2015, MNRAS, 446, 4232 [NASA ADS] [CrossRef] [Google Scholar]
 Akrami, Y., Fantaye, Y., Shafieloo, A., et al. 2014, ApJ, 784, L42 [NASA ADS] [CrossRef] [Google Scholar]
 Aluri, P. K., Pant, N., Rotti, A., & Souradeep, T. 2015, Phys. Rev. D, 92, 083015 [NASA ADS] [CrossRef] [Google Scholar]
 Axelsson, M., Fantaye, Y., Hansen, F. K., et al. 2013, ApJ, 773, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Baldi, P., Kerkyacharian, G., Marinucci, D., & Picard, D. 2009, Annals of Statistics, 37, 1150 [Google Scholar]
 Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15 [NASA ADS] [CrossRef] [Google Scholar]
 BenDavid, A., & Kovetz, E. D. 2014, MNRAS, 445, 2116 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, C. L., Halpern, M., Hinshaw, G., et al. 2003, ApJS, 148, 1 [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]
 Bond, J. R., & Efstathiou, G. 1987, MNRAS, 226, 655 [NASA ADS] [Google Scholar]
 Bond, J. R., Frolov, A. V., Huang, Z., & Kofman, L. 2009, Phys. Rev. Lett., 103, 071301 [NASA ADS] [CrossRef] [Google Scholar]
 Bueno Sanchez, J. 2014, Phys. Lett. B, 739, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Bunn, E. F., & Scott, D. 2000, MNRAS, 313, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Cai, Y.C., Cole, S., Jenkins, A., & Frenk, C. S. 2010, MNRAS, 407, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Cayón, L., Jin, J., & Treaster, A. 2005, MNRAS, 362, 826 [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., Schwarz, D. J., & Starkman, G. D. 2007, Phys. Rev. D, 75, 023507 [Google Scholar]
 Copi, C. J., Huterer, D., Schwarz, D. J., & Starkman, G. D. 2009, MNRAS, 399, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Copi, C. J., Huterer, D., Schwarz, D. J., & Starkman, G. D. 2015, MNRAS, 451, 2978 [NASA ADS] [CrossRef] [Google Scholar]
 Coulson, D., Crittenden, R. G., & Turok, N. G. 1994, Phys. Rev. Lett., 73, 2390 [NASA ADS] [CrossRef] [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., Turok, N., Vielva, P., MartínezGonzález, E., & Hobson, M. 2007, 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. 2011, 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. 2011, MNRAS, 412, 1038 [NASA ADS] [Google Scholar]
 Czech, B., Kleban, M., Larjo, K., Levi, T. S., & Sigurdson, K. 2010, J. Cosmol. Astropart. Phys., 12, 23 [NASA ADS] [CrossRef] [Google Scholar]
 de OliveiraCosta, A., Smoot, G. F., & Starobinsky, A. A. 1996, ApJ, 468, 457 [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]
 Desjacques, V. 2008, Phys. Rev. D, 78, 103503 [NASA ADS] [CrossRef] [Google Scholar]
 Efstathiou, G. 2004, MNRAS, 348, 885 [NASA ADS] [CrossRef] [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., Novikov, D. I., Lilje, P. B., Banday, A. J., & Górski, K. M. 2004b, 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]
 Fantaye, Y. 2014, ArXiv eprints [arXiv:1409.1114] [Google Scholar]
 Fantaye, Y., Marinucci, D., Hansen, F., & Maino, D. 2015, Phys. Rev. D, 91, 063501 [NASA ADS] [CrossRef] [Google Scholar]
 Feeney, S. M., Johnson, M. C., Mortlock, D. J., & Peiris, H. V. 2011, Phys. Rev. Lett., 107, 071301 [NASA ADS] [CrossRef] [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, J. Cosmol. Astropart. Phys., 1207, 049 [Google Scholar]
 Finelli, F., GarciaBellido, J., Kovacs, A., Paci, F., & Szapudi, I. 2016, MNRAS, 455, 1246 [NASA ADS] [CrossRef] [Google Scholar]
 Fixsen, D. J., Cheng, E. S., Gales, J. M., et al. 1996, ApJ, 473, 576 [NASA ADS] [CrossRef] [Google Scholar]
 Gay, C., Pichon, C., & Pogosyan, D. 2012, Phys. Rev. D, 85, 023011 [NASA ADS] [CrossRef] [Google Scholar]
 Gordon, C. 2007, ApJ, 656, 636 [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]
 Gott, J. R., Colley, W. N., Park, C.G., Park, C., & Mugnolo, C. 2007, MNRAS, 377, 1668 [NASA ADS] [CrossRef] [Google Scholar]
 Gruppuso, A., Finelli, F., Natoli, P., et al. 2011, MNRAS, 411, 1445 [NASA ADS] [CrossRef] [Google Scholar]
 Gruppuso, A., Natoli, P., Paci, F., et al. 2013,J. Cosmol. Astropart. Phys., 7, 47 [Google Scholar]
 Gurzadyan, V. G., Kashin, A. L., Khachatryan, H., et al. 2014, A&A, 566, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hajian, A. 2007, ArXiv eprints [arXiv:astroph/0702723] [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., Eriksen, H. K., & Lilje, P. B. 2009, ApJ, 704, 1448 [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]
 Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, A&AS, 47, 1 [NASA ADS] [Google Scholar]
 Hikage, C., Matsubara, T., Coles, P., et al. 2008, MNRAS, 389, 1439 [NASA ADS] [CrossRef] [Google Scholar]
 Hinshaw, G., Banday, A. J., Bennett, C. L., et al. 1996, ApJ, 464, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Hinshaw, G., Weiland, J. L., Hill, R. S., et al. 2009, ApJS, 180, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Hoftuft, J., Eriksen, H. K., Banday, A. J., et al. 2009, ApJ, 699, 985 [NASA ADS] [CrossRef] [Google Scholar]
 Hou, Z., Banday, A. J., & Górski, K. M. 2009, MNRAS, 396, 1273 [NASA ADS] [CrossRef] [Google Scholar]
 Inoue, K. T., & Silk, J. 2006, ApJ, 648, 23 [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]
 Kamionkowski, M., & Knox, L. 2003, Phys. Rev. D, 67, 063001 [NASA ADS] [CrossRef] [Google Scholar]
 Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. D, 55, 7368 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J., & Komatsu, E. 2013, Phys. Rev. D, 88, 101301 [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]
 Knutsson, H., & Westin, C.F. 1993, Conf. IEEE Comput. Soc., 515 [Google Scholar]
 Kogut, A., Spergel, D. N., Barnes, C., et al. 2003, ApJS, 148, 161 [NASA ADS] [CrossRef] [MathSciNet] [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]
 Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 [NASA ADS] [CrossRef] [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]
 Larson, D. L., & Wandelt, B. D. 2004, ApJ, 613, L85 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, D. L., & Wandelt, B. D. 2005, ArXiv eprints [arXiv:astroph/0505046] [Google Scholar]
 Liu, X., & Zhang, S. N. 2005, ApJ, 633, 542 [NASA ADS] [CrossRef] [Google Scholar]
 Marinucci, D., Pietrobon, D., Balbi, A., et al. 2008, MNRAS, 383, 539 [NASA ADS] [CrossRef] [Google Scholar]
 MartinezGonzá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., Vielva, P., Wiaux, Y., et al. 2007, J. Fourier Analysis and Applications, 13, 495 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, J. D., Feeney, S. M., Johnson, M. C., & Peiris, H. V. 2012, Phys. Rev. D, 85, 103502 [NASA ADS] [CrossRef] [Google Scholar]
 Mecke, K. R., Buchert, T., & Wagner, H. 1994, A&A, 288, 697 [NASA ADS] [Google Scholar]
 Mikkelsen, K., Næss, S. K., & Eriksen, H. K. 2013, ApJ, 777, 172 [NASA ADS] [CrossRef] [Google Scholar]
 Mitra, S., Rocha, G., Górski, K. M., et al. 2011, ApJS, 193, 5 [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]
 Nadathur, S., Lavinto, M., Hotchkiss, S., & Räsänen, S. 2014, Phys. Rev. D., 90, 103510 [NASA ADS] [CrossRef] [Google Scholar]
 Notari, A., & Quartin, M. 2015, J. Cosmol. Astropart. Phys., 6, 47 [Google Scholar]
 Planck Collaboration I. 2014, A&A, 571, A1 [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 XIX. 2014, A&A, 571, A19 [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 XXVII. 2014, A&A, 571, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2016, A&A, 594, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2016, A&A, 594, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2016, A&A, 594, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2016, A&A, 594, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2016, A&A, 594, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2016, A&A, 594, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2016, A&A, 594, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2016, A&A, 594, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2016, A&A, 594, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2016, A&A, 594, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2016, A&A, 594, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2016, A&A, 594, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIV. 2016, A&A, 594, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XV. 2016, A&A, 594, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVI. 2016, A&A, 594, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVII. 2016, A&A, 594, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2016, A&A, 594, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2016, A&A, 594, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XX. 2016, A&A, 594, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXI. 2016, A&A, 594, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2016, A&A, 594, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIII. 2016, A&A, 594, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIV. 2016, A&A, 594, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXV. 2016, A&A, 594, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2016, A&A, 594, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2016, A&A, 594, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVIII. 2016, A&A, 594, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Plaszczynski, S., Montier, L., Levrier, F., & Tristram, M. 2014, MNRAS, 439, 4048 [NASA ADS] [CrossRef] [Google Scholar]
 Pogosyan, D., Gay, C., & Pichon, C. 2009, Phys. Rev. D, 80, 081301 [NASA ADS] [CrossRef] [Google Scholar]
 Pontzen, A., & Peiris, H. V. 2010, Phys. Rev. D, 81, 103008 [NASA ADS] [CrossRef] [Google Scholar]
 Quartin, M., & Notari, A. 2015, J. Cosmol. Astropart. Phys., 1, 008 [Google Scholar]
 Rudnick, L., Brown, S., & Williams, L. R. 2007, ApJ, 671, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Savitzky, A., & Golay, M. J. E. 1964, J. Anal. Chem., 36, 1627 [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]
 Scott, D. 1991, A&A, 242, 1 [NASA ADS] [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]
 Stannard, A., & Coles, P. 2005, MNRAS, 364, 929 [NASA ADS] [CrossRef] [Google Scholar]
 Starobinsky, A. A. 1993, JETP Lett., 57, 622 [NASA ADS] [Google Scholar]
 Stevens, D., Scott, D., & Silk, J. 1993, Phys. Rev. Lett., 71, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Szapudi, I., Kovács, A., Granett, B. R., et al. 2015, MNRAS, 450, 288 [NASA ADS] [CrossRef] [Google Scholar]
 Tegmark, M., de OliveiraCosta, A., & Hamilton, A. 2003, Phys. Rev. D, 68, 123523 [NASA ADS] [CrossRef] [Google Scholar]
 Tomita, K. 2005, Phys. Rev. D, 72, 103506 [NASA ADS] [CrossRef] [Google Scholar]
 Tomita, K., & Inoue, K. T. 2008, Phys. Rev. D, 77, 103522 [NASA ADS] [CrossRef] [Google Scholar]
 Vanmarcke, E. 1983, in Random Fields (Cambridge, Massachusetts, USA: The MIT Press), March (Paper), ed. E. Vanmarcke, 372 [Google Scholar]
 Vielva, P. 2007, in SPIE Conf. Ser., 6701 [Google Scholar]
 Vielva, P. 2010, Adv. Astron., 2010, 592094 [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]
 Watson, W. A., Diego, J. M., Gottlöber, S., et al. 2014, MNRAS, 438, 412 [NASA ADS] [CrossRef] [Google Scholar]
 Wick, G. C. 1950, Phys. Rev., 80, 268 [NASA ADS] [CrossRef] [Google Scholar]
 Zaldarriaga, M., & Seljak, U. 1997, Phys. Rev. D, 55, 1830 [Google Scholar]
 Zhang, R., & Huterer, D. 2010, Astropart. Phys., 33, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, W. 2013, MNRAS, 433, 3498 [NASA ADS] [CrossRef] [Google Scholar]
 Zibin, J. P. 2014, ArXiv eprints [arXiv:1408.4442] [Google Scholar]
Appendix A: Generalized SavitzkyGolay polynomials
In the construction of optimal linear filters, one needs to combine information about the (statistically isotropic) CMB signal, anisotropic instrumental noise, masking to be applied for the elimination of foreground contributions, and a model for any nonGaussian signal for matched filtering. These can be combined in a general framework of normalized convolutions (Knutsson & Westin 1993), where the filtered field is defined as (A.1)where B is the (multiscale) filtering beam function, T is the temperature, a and w their respective weights, and ⋆ denotes the usual convolution operation (A.2)In the absence of a specific model for the nonGaussian signal, the beam functions can be taken to be orthogonal polynomials on a disc, weighted by some smoothing function, while the weights applied to the temperature maps are determined by the CMB and noise covariance.
In a simple approach, the information about the CMB signal can be utilized by prewhitening the map by convolving it with an isotropic beam function derived from the isotropic bestfit CMB power spectrum combined with a diagonal approximation to the instrumental noise covariance. After the componentseparated CMB maps are prewhitened, and the corresponding mask is applied to the resulting map, the multiscale filtering kernel b_{ℓ} is applied at various scales.
Fig. A.1 Generalized SavitzkyGolay polynomials are orthogonal to polynomials up to degree n on a disc, with smoothing weight applied. Upper left panel shows a few representative polynomial kernels (SSG21 in red, SSG42 in dark green, SSG84 in blue) and Gaussian (in black) as a function of radius (scaled to the same FWHM of 800′), lower left shows their harmonic space representation. Right column shows the prewhitening kernel for the SMICA temperature map on the top (in light blue), and the corresponding composite kernels (WHITE*SSG21, etc.) on the bottom (in the same colours). 
In this paper, the maps are prewhitened with the 2013 bestfit cosmological parameter CMB spectrum (Planck Collaboration XV 2014), coadded to an isotropic noise power spectrum derived from the halfmission, halfdifference noise maps appropriate for each componentseparation method. No adjustment is made either for the recalibration of the 2015 data relative to the nominal results that the cosmological spectrum is derived from, or for the mismatch in noise level between the halfmission, halfdifference and fullmission maps. This implies that the filtering is suboptimal, but the data and simulations are treated consistently so there should be no significant impact on the results. The resulting prewhitening beam function w_{ℓ} for the SMICA temperature map is shown in Fig. A.1.
The peak detector wavelets are taken to be SavitzkyGolay polynomials (Savitzky & Golay 1964), generalized to be defined on a disc with a polynomial smoothing weight function applied, as shown in Fig. A.1. A generalized spherical SavitzkyGolay kernel of order n and smoothing weight k (referred to as SSGnk in the text) is defined by a polynomial function of a radial variable x = sin(θ/ 2) / sin(θ_{max}/ 2), (A.3)which is normalized to have unit mean on a disc and is orthogonal to all nonconstant polynomials up to order n, (A.4)These are essentially highorder lowpass filters in harmonic space, but have compact support on the sphere. A few representative SavitzkyGolay polynomials are compared to a Gaussian kernel in Fig. A.1. Combined with prewhitening, the total effect of the filters applied is described by the composite beam functions shown in Fig. A.1.
One should note a slight ℓspace bandwidth mismatch between differently shaped kernels with the same FWHM value in real space, which is clear from the lower left panel of Fig. A.1. While not a problem in general, some care should be exercised when directly comparing results for different shape kernels. In particular, the ℓ value at which the filter kernel coefficient reaches b_{ℓ} = b_{max}/ 2 differs by a factor of 1.58 between the GAUSS and SSG84 kernels of the same FWHM.
Appendix B: Doppler boosting
The main effect of our relative motion with respect to the CMB rest frame is a dominant contribution to the CMB dipole (C_{1}); this is boosting of the monopole and has been detected previously (Kogut et al. 2003; Fixsen et al. 1996; Hinshaw et al. 2009). A subtler consequence of our motion is the boosting of all other multipoles. In fact, there are really two effects at work. The first is a modulation effect which increases power by approximately 0.25% in the direction of our motion and decreases it by the same amount in the opposite direction. This can equivalently be thought of as coupling between the multipoles ℓ and ℓ ± 1. The second is an aberration effect which shifts the apparent direction in which CMB photons arrive at our detectors toward the velocity direction.
Planck Collaboration XXVII (2014) reported a detection of this Doppler boosting, and an associated measurement of its velocity signature of 384 ± 78 (statistical) ± 115 (systematic) km s^{1} in the known dipole direction, (l,b) = (264°,48°). Here, we demonstrate that the Planck 2015 data release remains in agreement with this result, by considering the angular scales 500 ≤ ℓ ≤ 2000. However, since the simulations employed in the analysis partially contain the effects of Doppler boosting (as noted in Sect. 3 the aberration contribution was erroneously omitted), we report a consistency check rather than a detection.
It is useful to perform a harmonic transform on the peculiar velocity vector, (B.1)where only the L = 1 modes are nonzero. Following the convention in Planck Collaboration XXVII (2014), we rotate to an orthonormal basis, labelled β_{} (along the expected velocity direction), β_{×} (parallel to the Galactic plane), and β_{⊥} (the remaining vector).
Significance measures for the β estimates for the 143 × 217 data set.
The peculiar velocity is detected using estimators that pick out the offdiagonal components of the CMB covariance matrix (B.2)The weight function W^{βv} is a sum of the modulation (b_{v}W^{τ}) and aberration (W^{φ}) effects. We quote results based on orthogonalized weight matrices, Due to the clear connection between the velocity estimators and those used for the lensing analysis, we adopt the same data (143 GHz and 217 GHz sky maps, with dust foregrounds removed using the 857 GHz data as a template) and mask as used in Planck Collaboration XV (2016). The results summarized in Table B.1 show a slight excess of signal in the dipole direction of the data compared to simulations. This is due to the simulations used containing the modulation, but not aberration, part of the Doppler boost signal.
Appendix C: Generalized modulation estimator
Consider a parameter X that the (primary) CMB power spectrum is dependent on. Let X have a dipolar dependence of the form (this could correspond to a gradient in X across our observable volume), where X_{0} is the average value, is the direction to the last scattering surface, and is the gradient direction. To linear order in ΔX/X, the measured spherical harmonics coefficients are given by (C.1)where the are the unmodulated statistically isotropic modes. The are coupling coefficients given by where From Eq. (C.1) we can find the covariance matrix to first order in the components ΔX_{M}: (C.6)where δC_{ℓℓ + 1} = dC_{ℓ}/ dX + dC_{ℓ + 1}/ dX. To determine the bestfit parameters, we proceed by maximizing the CMB likelihood function (C.7)where d is the CMB temperature data. Equation (C.7) is maximized for the ΔX_{M} that satisfy (C.8)From Eq. (C.6) it is clear that the CMB covariance can be decomposed into an isotropic part (C_{ℓ}) and a small anisotropic part proportional to ΔX_{M}. By inverting Eq. (C.6) and using the orthogonality of the , we can determine the bestfit parameters and , to first order in the anisotropy. These estimators are the fullsky, nonoise versions of Eqs. (44) and (45).
Errors can easily be found by expanding the loglikelihood about the bestfit parameters. The Fisher matrix is defined as (C.11)Upon switching bases, we find We can then assign the standard errors, .
Appendix D: Weightedvariance modified shape function estimator
The BipoSH representation characterizes the offdiagonal elements in the covariance matrix and is a generalization of the angular power spectrum, C_{ℓ}, (D.1)In general, it is not possible to analyse the full sky even for componentseparated maps, due to the presence of residual contributions from diffuse Galactic emission and point sources. However, the application of a mask leads to coupling between the spherical harmonic modes. Hence, the correlation function is no longer described only by C(θ) or the power spectrum C_{ℓ}, and other quantities are required to completely quantify the statistical field.
We obtain an analytic expression for the observed BipoSH coefficients after the application of a mask in terms of the corresponding coefficients of the unmasked sky, and those of the mask itself, (D.2)where , are the BipoSH coefficients of the masked sky map, correspond to the BipoSH coefficients of the unmasked sky, are the BipoSH coefficient of the mask itself, are the ClebschGordon coefficients, and the term {} in Eq. (D.2) is the 9j −symbol. This quantifies the coupling between the BipoSH coefficients of the CMB sky map and those of the mask itself.
The underlying CMB sky may have deviations from statistical isotropy, as discussed in Sect. 6.4, due either to a dipole modulation (L = 1) of unknown origin, or to Doppler boosting (L = 1) of the temperature field. The BipoSH coefficients of such statistical isotropyviolating fields can be given by (D.3)Here corresponds to the BipoSH coefficients of the unknown but statistically isotropic CMB field. This couples with BipoSH coefficients of the mask to introduce a mean field linear bias , which is estimated from simulations and subtracted from the BipoSH coefficients obtained from the masked sky. The φ_{LM} are the spherical harmonic coefficients of the field that breaks statistical isotropy, and is the shape function. Shape functions for dipole modulation and Doppler boosting are given in Eqs. (54) and (56), respectively.
Due to symmetries of the mask, which is largely defined by foreground residuals towards the Galactic plane, the dominant BipoSH modes of the mask correspond to J = { 0,2 } ,K = 0. Hence, for all practical purposes, signal is retained in the L = 1 mode itself, although masking modifies the shape function, now defined as the modified shape funtion in the rest of the text. A weighted variance modified shape function is defined as (D.4)where and the weights are chosen such that .
Here is the MSF, which can be evaluated as (D.5)The weights are then given by (D.6)
All Tables
Lowertail probabilities for the variance, skewness, and kurtosis of the componentseparated maps.
Lowertail probabilities for the Npdf χ^{2} statistics derived from the Planck 2015 componentseparated maps at N_{side} = 16 and 32.
Probabilities of obtaining values for the χ^{2} statistic of the Npoint functions for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2015 temperature CMB maps with resolution parameter N_{side} = 64, estimated using the Commander, NILC, SEVEM, and SMICA methods.
Probability as a function of resolution for the unnormalized, classical Minkowski functionals.
Modified upper tail probability (mUTP ) for the cold (top) and hot (bottom) areas.
Modified upper tail probability (mUTP) for the KS test, comparing the data with the median peak CDF derived from simulations.
Modified upper tail probability (mUTP) for the total peak count, comparing the data with the peak count CDF derived from simulations.
Lowertail probability for the variance, skewness, and kurtosis of the low resolution N_{side} = 16 componentseparated maps obtained with the common mask and two extended versions thereof.
Probabilities of obtaining values for the S_{1 / 2} and statistics for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2015 temperature CMB maps with resolution parameter N_{side} = 64, estimated using the Commander, NILC, SEVEM, and SMICA maps.
Probabilities of obtaining values for the χ^{2} statistic for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2015 temperature CMB maps with resolution parameter N_{side} = 64, estimated using the Commander, NILC, SEVEM, and SMICA maps.
Probabilities of obtaining values for the χ^{2} statistic and ratio of χ^{2} of the Npoint functions shown in Fig. 17 for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2015 CMB maps estimated on northern and southern ecliptic hemispheres.
Probabilities of obtaining values for the χ^{2} statistic and ratio of χ^{2} of the Npoint functions shown in Fig. 18 for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the SMICA map on hemispheres defined by the ecliptic (first column), Doppler boost (DB, second column), and dipolar modulation (DM, third column) reference frames.
Constraints on the quadrupolar modulation, determined from the Commander, NILC, SEVEM, and SMICA foregroundcleaned maps.
Lowertail probability for the S^{±} statistics of the componentseparated maps at N_{side} = 16 and N_{side} = 32.
pvalues for the variance asymmetry measured by 8° discs for the four componentseparated temperature maps and different highpass filter scales.
Amplitude (A) and direction of the lowℓ dipole modulation signal determined from the QML analysis for the range ℓ ∈ [2,64].
Amplitude (A) and direction of the dipole modulation in Galactic coordinates as estimated for the multipole range ℓ ∈ [2,64] using a BipoSH analysis.
Doppler boost amplitude ( β ) and direction in Galactic coordinates derived over the multipole range ℓ ∈ [640,1024] as evaluated from a BipoSH analysis.
Lowertail probability for the variance, skewness, and kurtosis of the LklCommander map.
Lowertail probability for the variance, skewness, and kurtosis of the LklCommander map compared to the component separated maps, obtained using the lowℓ likelihood mask LklT_{16}94.
Probabilities to exceed the observed values of the χ^{2} statistics for the LklCommander and Commander maps at N_{side} = 64.
Probabilities for obtaining values of the S_{1 / 2} and statistics for the simulations at least as large as the observed values of the statistic estimated from the LklCommander and Commander maps using the LklT_{64}92 and UT_{64}67 masks, respectively. We also show the corresponding estimation of the global pvalue for the S(x) statistic.
Probabilities for obtaining values of the χ^{2} statistic for the simulations at least as large as the observed values of the statistic estimated from the LklCommander and Commander maps using the LklT_{64}92 and UT_{64}67 masks, respectively.
Probabilities for obtaining values of the χ^{2} statistic and ratio of χ^{2} of the Npoint functions for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic on the northern and southern ecliptic hemispheres estimated from the LklCommander and Commander maps using the LklT_{64}92 and UT_{64}67 masks, respectively.
pvalues for the variance asymmetry measured by different discs from the Planck 2015 LklCommander and Commander temperature solutions using the LklT_{256}93 and UT_{256}73 masks, respectively.
Localvariance dipole directions measured by 8° discs for the Planck 2015 LklCommander and Commander temperature solutions.
Amplitude (A) and direction of the dipole modulation in Galactic coordinates as estimated for the multipole range ℓ ∈ [2,64] using the BipoSH analysis on LklCommander and Commander maps. The former results were derived using the LklT_{256}93 mask; the latter are those determined previously in Sect. 6.4.
pvalues of the T, Q_{r}, and U_{r} angular profiles computed from the stacking of hot and cold spots selected above the ν = 0 and ν = 3 thresholds.
pvalues of Δμ_{T} computed from the stacking of hot and cold spots selected above the ν = 0 and ν = 3 thresholds from the Commander, NILC, SEVEM, and SMICA maps.
All Figures
Fig. 1 Variance, skewness, and kurtosis for the four different componentseparation methods – Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) – compared to the distributions derived from 1000 Monte Carlo simulations. 

In the text 
Fig. 2 Npoint correlation functions determined from the N_{side} = 64Planck CMB 2015 temperature maps. Results are shown for the 2point, pseudocollapsed 3point (upper left and right panels, respectively), equilateral 3point, and connected rhombic 4point functions (lower left and right panels, respectively). The red dotdotdotdashed, orange dashed, green dotdashed, and blue long dashed lines correspond to the Commander, NILC, SEVEM, and SMICA maps, respectively. Note that the lines lie on top of each other. The black solid line indicates the mean determined from 1000 SMICA simulations. The shaded dark and light grey regions indicate the corresponding 68% and 95% confidence regions, respectively. See Sect. 4.3 for the definition of the separation angle θ. 

In the text 
Fig. 3 Needlet space MFs for Planck 2015 data using the four componentseparated maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue); the grey regions, from dark to light, correspond, respectively, to 1, 2, and 3σ confidence regions estimated from the 1000 FFP8 simulations processed by the Commander method. The columns from left to right correspond to the needlet parameters j = 4,6, and 8, respectively; the jth needlet parameter has compact support over multipole ranges [2^{j−1},2^{j + 1}]. The ℓ_{c} = 2^{j} value indicates the central multipole of the corresponding needlet map. Note that to have the same range at all the needlet scales, the vertical axis has been multiplied by a factor that takes into account the steady decrease of the variance of the MFs as a function of scale. 

In the text 
Fig. 4 Histograms of χ^{2} for the Planck 2015 Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) foregroundcleaned maps analysed with the common mask. The χ^{2} is obtained by combining the three MFs in needlet space with an appropiate covariance matrix. The histograms are for the FFP8 simulations, while the vertical lines are for the data. The figures from left to right are for the needlet scales j = 4,6, and 8, with the central multipoles ℓ_{c} = 2^{j} shown in each panel. 

In the text 
Fig. 5 Comparison of the window functions (normalized to have equal area) for the SMHW (blue), GAUSS (yellow), and SSG84 (magenta) filters. The scales shown are 25′ (top) and 250′ (bottom). 

In the text 
Fig. 6 Modified upper tail probabilities (mUTP) obtained from the analyses of the filter coefficients as a function of the filter scale R for the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) sky maps. From left to right, the panels correspond to standard deviation, skewness, and kurtosis results, when determined using the SMHW (top), GAUSS (middle), and SSG84 (bottom) filters. The squares represent UTP values above 0.5, whereas circles represent UTP values below 0.5. 

In the text 
Fig. 7 Modified upper tail probabilities (mUTP) obtained from the analyses of the SMHW coefficients as a function of the wavelet scale R for the SEVEM100 (blue), SEVEM143 (yellow), SEVEM217 (magenta), and SEVEM (green) maps. From left to right, the panels correspond to the standard deviation, skewness, and kurtosis. 

In the text 
Fig. 8 Modified upper tail probabilities (mUTP) obtained from the analyses of the SMHW coefficients as a function of the wavelet scale R for the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) halfring halfdifference noise estimates. From left to right, the panels correspond to the standard deviation, skewness, and kurtosis. 

In the text 
Fig. 9 Cold and hot areas for thresholds ν> 3.0 as determined from the SEVEM temperature map. From top to bottom, the maps are for SMHW scales of R = 200′, R = 250′, R = 300′, and R = 400′. 

In the text 
Fig. 10 Cumulative density function of the peak distribution for the SMICA CMB temperature map. The top row shows the peak CDF for maps filtered with a GAUSS kernel of 40′ FWHM. The bottom row shows the corresponding peak CDF for an SSG84 kernel of 800′ FWHM. The spectral shape parameter γ (see Eq. (27)) is the bestfit value for the simulated ensemble, as indicated by the cyan circle in Fig. 11. Similar results are obtained for the other componentseparation methods. 

In the text 
Fig. 11 Distribution of bestfit Gaussian peak CDF spectral shape parameters, σ and γ (as defined in Eq. (27)), recovered from 1000 simulations, as indicated by the black dots and the smoothed density map and compared to those derived for the observed sky (shown by the red star). The blue contours enclose 68% and 95% of the parameter distribution, and the cyan circle represents the bestfit parameters for the median peak CDF determined from simulations. The upper panel shows the peak CDF parameters for the SMICA map filtered with a GAUSS kernel of 40′ FWHM. The lower panel shows the corresponding peak CDF for an SSG84 kernel of 800′ FWHM. Similar results are obtained for the other componentseparation methods. 

In the text 
Fig. 12 Cumulative density function of the mean amplitude of all extrema, maxima (red) and minima (blue), derived from simulations, compared to the equivalent values observed for the SMICA CMB temperature map. The upper panel shows the peak mean amplitudes for maps filtered with a GAUSS kernel of 40′ FWHM. The lower panel shows the corresponding peak CDF for an SSG84 kernel of 800′ FWHM. Similar results are obtained for the other component separation methods. Since the filter kernel normalization is free, and the prewhitened map to which the filter is applied is dimensionless, the plots are essentially in arbitrary units. 

In the text 
Fig. 13 Fraction of the Gaussian random field realizations in which the coldest peak is as cold as or colder than that observed, as a function of SMHW filter scale for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). 

In the text 
Fig. 14 Peak positions and CDF rank summarized for all filtering scales. The three skyview panels in the top row show a Lambert projection of the north pole, the usual full sky Mollweide projection, and a Lambert projection of the south pole. The lower panel shows the peak heights (in percentile of the peak distribution on the horizontal axis) as a function of filter scale (on the vertical axis, in logarithmic scale), truncated to larger scales for clarity. Circles represent peaks (nodes of the graph) coloured according to their percentile level, and scaled according to kernel size. Black lines represent edges connecting peaks at different scales (according to a minimal distance measure). The components connected to the coldest and hottest peaks at any scale are highlighted by thicker edges, and are navy blue and dark red in colour. Note that there are thick lines that do not touch the 0 and 1 percentiles in the plot view. Those edges are connected to extreme percentile values, but at scales smaller than those shown in the plot. The Cold Spot is represented by the connected nodes that have the smallest percentiles except for the coarsest scale in the plot view. 

In the text 
Fig. 15 Distribution of node degrees in the multiscale peak linkage graph determined for the SMICA map (cyan), compared with the median (red line), first to third quartile (blue box), and 95% (whiskers) derived from 1000 FFP8 simulations. 

In the text 
Fig. 16 Lower tail probability of the variance (top panel), skewness (centre panel), and kurtosis (bottom panel) obtained at different resolutions from the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) sky maps. 

In the text 
Fig. 17 Difference of the Npoint correlation functions determined from the N_{side} = 64Planck CMB 2015 temperature estimates and the corresponding means estimated from 1000 simulations. Results are shown for the 2point, pseudocollapsed 3point (upper left and right panels, respectively), equilateral 3point, and connected rhombic 4point functions (lower left and right panels, respectively). Correlation functions are shown for the analysis performed on northern (blue) and southern (red) hemispheres determined in the ecliptic coordinate frame. The solid, dashed, dotdashed, and dotted lines correspond to the Commander, NILC, SEVEM, and SMICA maps, respectively. Note that the lines lie on top of each other. The shaded dark and light grey regions indicate, for reference, the 68% and 95% confidence regions, respectively, determined from the SMICA simulations. 

In the text 
Fig. 18 Difference of the Npoint correlation functions determined from the N_{side} = 64PlanckSMICA CMB 2015 temperature estimates and the corresponding means estimated from 1000 simulations. Results are shown for the 2point, pseudocollapsed 3point (upper left and right panels, respectively), equilateral 3point, and connected rhombic 4point functions (lower left and right panels, respectively). Correlation functions are shown for the analysis performed on negative (blue) and positive (red) hemispheres determined in the ecliptic (solid lines), Doppler boost (DB, dashed lines), and dipole modulation (DM, dotdashed lines) coordinate frames. The shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively. 

In the text 
Fig. 19 Ratio R^{TT}(ℓ_{max}) for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) determined at N_{side} = 32. The shaded grey regions indicate the distribution of the statistic derived from the SMICA MC simulations, with the dark, lighter, and light grey bands corresponding to the 1, 2, and 3σ confidence levels. 

In the text 
Fig. 20 Lowertail probability of the pointparity estimator for Commander (red), NILC (orange), SEVEM, (green), and SMICA (blue). 

In the text 
Fig. 21 Histograms of the S^{+} (top panel) and S^{−} (bottom panel) statistic. The vertical lines show the minimum value for the estimator computed at N_{side} = 32 for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) maps. The grey area shows the same quantity computed from 1000 simulated SMICA maps. 

In the text 
Fig. 22 KSdeviation of the peak distribution for 70° radius discs centred on the positive and negative asymmetry directions determined from the SMICA CMB temperature map in Sect. 6.2. From top to bottom, the plots correspond to maps filtered with a GAUSS kernel of 40′ FWHM, an SSG84 filter of 500′ FWHM, and an SSG84 filter of 800′ FWHM, respectively. 

In the text 
Fig. 23 Map of −log _{10}(UTP) for peak counts in the Planck 40′ GAUSSfiltered temperature data, where each pixel encodes the probability determined for a 30° diameter disc centred on it. 

In the text 
Fig. 24 Map of −log _{10}(UTP) for the twosample KSdeviation where each pixel encodes the probability determined for a 30° diameter disc centred on it, as computed from the Planck 40′ GAUSSfiltered temperature data. 

In the text 
Fig. 25 Top: temperature patch centred on the Cold Spot. Bottom: peak merger tree within the Cold Spot region. The figure shows a region centred on the Cold Spot location in gnomonic projection, with all the peaks in SSG84filtered maps with FWHM ranging from 80′ to 1200′ overlaid on the same plot. The size of the coloured circles is proportional to the filtering scale. The colour corresponds to the peak value, normalized in units of σ for a given filter scale. In both panels the data are from the SMICA CMB map at full resolution. 

In the text 
Fig. 26 From left to right: mean, variance, skewness, and kurtosis angular profiles computed for rings at radii θ centred on the Cold Spot position for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The expected value obtained from the simulations is denoted by the black dashed line and the dark and light grey regions represent the 1σ and 2σ intervals. 

In the text 
Fig. 27 Upper panel: pvalues for variance asymmetry measured as the number of simulations with localvariance dipole amplitudes larger than those inferred from the data, as a function of disc radius for the four componentseparated maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), and for unfiltered and highpassfiltered cases. For the filtered case, the Commander curve is covered by the SMICA curve for small (r_{disk} ≤ 8) disks, and by the SEVEM curve for large disks (r_{disk}> 8). Lower panel: localvariance dipole directions for the SMICA map. The colours, as indicated by the colourbar, correspond to different values of the highpass filter central multipole ℓ_{0}. The size of a marker disc corresponds, from small to large, to the size of the disc used in the analysis, namely 4°, 12°, 20°, and 70°. The dipole directions from the Commander, NILC, and SEVEM componentseparation methods are consistent with the case shown here. The lowℓ and WMAP9 directions are identical to those in Fig. 35. 

In the text 
Fig. 28 Upper panel: localvariance dipole amplitude for 8° discs as a function of the central multipole of the highpass filter, ℓ_{0}, for the four componentseparation methods, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The grey regions, from dark to light, correspond, respectively, to 1σ, 2σ, and 3σ percentiles from the 1000 FFP8 simulations processed by the Commander method. Lower panel: meansubtracted and inversevarianceweighted localvariance map for the 8° discs and for the Commander componentseparation method; each pixel is given in terms of the lower and uppertail probability of the measured value on that pixel compared to the values from the simulations. The pixels in grey correspond to the centres of the 8° discs on which the number of unmasked pixels in the full resolution map is lower than our threshold. The black curve superposed on the map indicates the boundary of the opposing hemispheres along the asymmetry axis. It is clear that the largest fraction of >95% outliers (red pixels) lie on the positive amplitude hemisphere of the local variance dipole, while the <5% outliers (blue pixels) are on the opposite hemisphere. The corresponding maps for NILC, SEVEM, and SMICA are very similar to the one shown here. 

In the text 
Fig. 29 Top: marginal constraints on the dipole modulation amplitude, as derived from Planck 2015 temperature observations at a smoothing scale of 5° FWHM for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The plot corresponds directly to Fig. 32 of Planck Collaboration XXIII (2014). The Commander, SEVEM, and SMICA posteriors coincide almost perfectly both internally, and with the corresponding SMICA 2013 posterior, shown as a dashed black line. Bottom: corresponding marginal twodimensional constraints on the lowℓ power spectrum amplitude and tilt, (q,n), defined relative to the bestfit Planck 2015 ΛCDM model. 

In the text 
Fig. 30 Probability determined from the QML analysis for a Monte Carlo simulation to have a larger dipole modulation amplitude than the Commander (red), NILC (orange), SEVEM (green), or SMICA (blue) data sets, with (top panel) ℓ_{min} = 2 or (bottom panel) ℓ_{min} = 100. No significant modulation is found once the lowℓ signal is removed. We emphasize that the statistic here is cumulative and apparent trends in the curves can be misleading. 

In the text 
Fig. 31 Probability determined from the QML analysis for obtaining a dipole modulation amplitude at least as anomalous as the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) data sets, for the range ℓ ∈ [10,ℓ_{max}]. The vertical line corresponds to ℓ_{max} = 132 which was used as the search limit in Bennett et al. (2011). The probability grows approximately logarithmically with ℓ_{max}. This means that the adopted probability to exceed is fortunately not very sensitive to ℓ_{max}, and for any reasonable choice is above 10%. 

In the text 
Fig. 32 Top: measured dipole modulation (L = 1) power in nonoverlapping CMB multipole bins for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) as determined from a BipoSH analysis of the data. The power in the dipole of the modulation field is a χ^{2}distributed variable with 3 degrees of freedom. The shaded regions in the plot depict, in darkgrey, grey, and lightgrey respectively, the 1, 2, and 3σ equivalent intervals of the distribution function derived from simulations, while the solid black line denotes its median. Significant power in the dipole modulation is seen to be limited to ℓ = 2–64 and does not extend to higher multipoles. Bottom: dipole modulation direction as determined from the SMICA map. The directions found from the other component separation maps are consistent with this analysis. The coloured circles denote the central value of the multipole bin used in the analysis, as specified in the colour bar. The lowℓ and WMAP9 directions are identical to those in Fig. 35. 

In the text 
Fig. 33 Top: measured dipole modulation power in cumulative CMB multipole bins for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) as determined from a BipoSH analysis of the data.. Colour coding as in Fig. 32. Note that the measurements in cumulative bins indicate a power in excess of 2σ up to multipole ℓ_{max} ~ 320. The value on the horizontal axis denotes the maximum multipole used in the analysis, with ℓ_{min} = 2. Bottom: modulation dipole direction as recovered from the SMICA map. The directions found from the other componentseparation maps are consistent with these directions. The colourcoded points represent the directions recovered for the specific ℓ_{max} used in the analysis, with ℓ_{min} = 2. The lowℓ and WMAP9 directions are identical to those in Fig. 35. 

In the text 
Fig. 34 Top: amplitude  β  of the Doppler boost from the SEVEM100, SEVEM143, and SEVEM217 maps for different multipole bins determined using a BipoSH analysis. The maximum multipole of each bin is fixed at ℓ_{max} = 1024, while ℓ_{min} is incremented from ℓ = 2 to ℓ = 640 in steps of Δℓ = 128. The dashed line corresponds to the actual dipole boost amplitude,  β  = 1.23 × 10^{3}. Bottom: Doppler boost direction measured in Galactic coordinates from SEVEM217. The coloured circles denote ℓ_{min} used in the analysis, while ℓ_{max} = 1024 is held fixed. The lowℓ and WMAP9 directions are identical to those in Fig. 35. 

In the text 
Fig. 35 Dipole directions for independent 100multipole bins of the local power spectrum distribution from ℓ = 2 to 1500 in the SMICA map with the common mask applied. We also show the preferred dipolar modulation axis (labelled as “lowℓ”) derived in Sect. 6.2, as well as the total direction for ℓ_{max} = 600 determined from WMAP9 (Axelsson et al. 2013). The average directions determined from the two multipole ranges ℓ ∈ [2,300] and ℓ ∈ [750,1500] are shown as blue and red rings, respectively. The error on the derived direction that results from masking the data is about 60°, with only small variations related to bin size. 

In the text 
Fig. 36 Derived pvalues for the angular clustering of the power distribution as a function of ℓ_{max}, determined for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), based on 500 simulations. For SMICA, the pvalues based on 2500 simulations are also shown (black). The pvalues are based on the fraction of simulations with a higher RS, determined over the ℓrange up to the given ℓ_{max}, compared to the data. The results shown here have been marginalized over bin sizes in the range Δℓ = 8 to Δℓ = 32. 

In the text 
Fig. 37 Derived pvalues for the angular clustering analysis as a function of ℓ_{max}, determined from SMICA based on 2500 simulations. The pvalues are based on the fraction of simulations with a higher Rayleigh statistic up to the given ℓ_{max} than in the data. The RS here is calculated over all pairs of dipole directions where one dipole in each pair is computed in the range [ℓ_{lim},ℓ_{max}], and the other is determined in the range [2,ℓ_{lim}]. The plot shows pvalues for ℓ_{lim} = 300 (purple), ℓ_{lim} = 400 (yellow), ℓ_{lim} = 500 (pink), and ℓ_{lim} = 700 (cyan). The results have been marginalized over bin sizes in the range Δℓ = 8 to Δℓ = 32. 

In the text 
Fig. 38 Rayleigh statistic pvalues determined from the QML analysis as a function of ℓ_{max} for the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) data sets, with (top panel) ℓ_{min} = 2 and (bottom panel) ℓ_{min} = 100. The general pattern of peaks is very similar to that in Fig. 30. We emphasize that the statistic here is cumulative and as such trends in the curves can be misleading. 

In the text 
Fig. 39 Probability to exceed (PTE) the pvalue of the signal from the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) data at ℓ = 230–240 (which is the multipole range with the most significant deviation) when searching over a range of multipoles up to ℓ_{max}, for the RS determined from the QML analysis. Much like the equivalent curve for dipole modulation, the PTE appears to grow approximately logarithmically with ℓ_{max}. 

In the text 
Fig. 40 Npoint correlation functions determined from the N_{side} = 64Planck CMB 2015 temperature estimates. Results are shown for the 2point, pseudocollapsed 3point (upper left and right panels, respectively), equilateral 3point, and connected rhombic 4point functions (lower left and right panels, respectively). The brown three dotdashed, purple dashed, and red dotdashed lines correspond to the LklCommander map analysed using the lowℓ and common masks and the Commander map analysed using the common mask, respectively. Note that the dashed and dotdashed lines lie on top of each other. The black solid line indicates the mean for 1000 MC simulations. The shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively, estimated using 1000 Commander simulations. See Sect. 4.3 for the definition of the separation angle θ. 

In the text 
Fig. 41 From top to bottom, T, Q, U, Q_{r}, and U_{r} stacked images (in μK units) extracted around temperature hot spots selected above the null threshold (ν = 0) in the Commander sky map for data (left column) and an equivalent simulation (right column). The horizontal and vertical axes of the flatsky projection are labelled in degrees. 

In the text 
Fig. 42 Radial profile μ_{T}(ϖ) derived from the stacked temperature image (see Fig. 41 or 45). The denominators σ_{0} and σ_{2} are the theoretical rms values of CMB T and ∇^{2}T, respectively. The theoretical ⟨ μ_{T}(ϖ) ⟩ is a linear combination of ⟨ T(ϖ)(T(0) /σ_{0}) ⟩ (green dashdotted line) and ⟨ T(ϖ)(−∇^{2}T(0)) /σ_{2}) ⟩ (blue dotted line). For all four componentseparated maps, the deviation of μ_{T} from the ensemble mean ⟨ μ_{T} ⟩ of the fiducial model (here the Planck 2015 ΛCDM best fit) is consistent with cosmic variance, and can be related to the lowℓ power deficit. The example powerdeficit ⟨ μ_{T} ⟩ (purple dashed line) is the theoretical prediction of ⟨ μ_{T} ⟩ if the fiducial model C_{ℓ}s are reduced by 10% in the range 2 ≤ ℓ ≤ 50. 

In the text 
Fig. 43 Mean radial profiles of T, Q_{r}, and U_{r} in μK obtained for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). Each individual panel contains (top) the mean radial profiles and (bottom) the differences (denoted “Diff”) between the mean profiles of the data and those computed from the ensemble mean of the simulations. Results based on stacks around temperature hot and cold spots are shown in the left and right columns, respectively. Upper plots present results for peaks selected above the null threshold, while lower plots show the equivalent results for peak amplitudes above (hot spots) or below (cold spots) 3 times the dispersion of the temperature map. The black dots (connected by dashed lines) depict the mean profiles and the shaded regions correspond to the 1σ (68%) and 2σ (95%) error bars. The mean profiles and error bars are determined from SEVEM simulations. Note that the Diff curves for each componentseparation method are computed using the corresponding ensemble average, although only the ensemble average from SEVEM is shown here. 

In the text 
Fig. 44 χ^{2} distributions obtained from the T (left column), Q_{r} (middle column), and U_{r} (right column) mean radial profiles centred on temperature hot spots selected above the null threshold (upper row) and three times the dispersion of the map (bottom row). The black lines correspond to the theoretical χ^{2} distribution with 12 degrees of freedom, whilst the histograms show the distributions determined from 100 simulations computed through the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) pipelines. The vertical lines represent the χ^{2} values obtained from the data. 

In the text 
Fig. 45 Comparison between unoriented stacking (upper panels) and oriented stacking (lower panels) of temperature patches around temperature hot spots selected above the null threshold (ν = 0). The left panels are the stacked SMICA maps, and the right panels their corresponding NFEAs. The image units are μK. 

In the text 
Fig. 46 Stacked Q_{r} image around temperature hot spots selected above the null threshold (ν = 0) in the SMICA sky map. The left panel corresponds to the observed data and the right panel shows the NFEA. The image units are μK. The headless vectors (black solid lines) are the polarization directions for stacked Q_{stack}, U_{stack}. The lengths of the headless vectors are proportional to the polarization amplitude P_{stack}. 

In the text 
Fig. 47 Oriented stacking of polarization fields (Q, U) on temperature maxima (upper panels) and P_{T} maxima (lower panels). In both cases the threshold ν = 0 is used and the orientation is chosen such that U_{T} = 0 and Q_{T} ≥ 0 on the central peak. The image units are μK. The left panels are the stacked SMICA maps, and the right panels their NFEAs. See Fig. 46 for the meaning of the headless vectors (black dashed lines). 

In the text 
Fig. 48 Top: Emode maps stacked on the unoriented Emode maxima computed above the null threshold ν = 0. Bottom: Q stacked around oriented polarization amplitude (P) maxima. In this case, no threshold is used and the orientation is chosen such that U = 0 and Q ≥ 0 on the central peak. The left panels are the stacked SMICA maps, and the right panels their corresponding NFEAs. See Fig. 46 for the meaning of the headless vectors (black dashed lines). The image units are μK. 

In the text 
Fig. A.1 Generalized SavitzkyGolay polynomials are orthogonal to polynomials up to degree n on a disc, with smoothing weight applied. Upper left panel shows a few representative polynomial kernels (SSG21 in red, SSG42 in dark green, SSG84 in blue) and Gaussian (in black) as a function of radius (scaled to the same FWHM of 800′), lower left shows their harmonic space representation. Right column shows the prewhitening kernel for the SMICA temperature map on the top (in light blue), and the corresponding composite kernels (WHITE*SSG21, etc.) on the bottom (in the same colours). 

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.