Issue 
A&A
Volume 641, September 2020
Planck 2018 results



Article Number  A7  
Number of page(s)  61  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201935201  
Published online  11 September 2020 
Planck 2018 results
VII. Isotropy and statistics of the CMB
^{1}
AIM, CEA, CNRS, Université ParisSaclay, Université ParisDiderot, Sorbonne Paris Cité, 91191 GifsurYvette, France
^{2}
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
^{3}
African Institute for Mathematical Sciences, 68 Melrose Road, Muizenberg, Cape Town, South Africa
^{4}
Aix Marseille Univ., CNRS, CNES, LAM, Marseille, France
^{5}
Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{6}
Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZuluNatal, Westville Campus, Private Bag X54001, Durban 4000, South Africa
^{7}
CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada
^{8}
CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
^{9}
Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA
^{10}
California Institute of Technology, Pasadena, CA, USA
^{11}
Centre for Theoretical Cosmology, DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
^{12}
Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, CA, USA
^{13}
Département de Physique Théorique, Université de Genève, 24 quai E. Ansermet, 1211 Genève 4, Switzerland
^{14}
Département de Physique, École normale supérieure, PSL Research University, CNRS, 24 rue Lhomond, 75005 Paris, France
^{15}
Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{16}
Departamento de Física, Universidad de Oviedo, C/ Federico García Lorca, 18, Oviedo, Spain
^{17}
Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{18}
Department of Mathematics, University of Stellenbosch, Stellenbosch 7602, South Africa
^{19}
Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC, Canada
^{20}
Department of Physics & Astronomy, University of the Western Cape, Cape Town 7535, South Africa
^{21}
Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, Helsinki, Finland
^{22}
Department of Physics, Princeton University, Princeton, NJ, USA
^{23}
Department of Physics, University of California, Santa Barbara, CA, USA
^{24}
Department of Physics, University of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, IL, USA
^{25}
Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, Via Marzolo 8, 35131 Padova, Italy
^{26}
Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
^{27}
Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 2, Roma, Italy
^{28}
Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria, 16, Milano, Italy
^{29}
Dipartimento di Fisica, Università degli Studi di Trieste, Via A. Valerio 2, Trieste, Italy
^{30}
Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy
^{31}
European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, Madrid, Spain
^{32}
European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{33}
Gran Sasso Science Institute, INFN, Viale F. Crispi 7, 67100 L’Aquila, Italy
^{34}
Haverford College Astronomy Department, 370 Lancaster Avenue, Haverford, PA, USA
^{35}
Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland
^{36}
INAF – OAS Bologna, Istituto Nazionale di Astrofisica – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Area della Ricerca del CNR, Via Gobetti 101, 40129 Bologna, Italy
^{37}
INAF – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, Padova, Italy
^{38}
INAF – Osservatorio Astronomico di Trieste, Via G.B. Tiepolo 11, Trieste, Italy
^{39}
INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy
^{40}
INAF/IASF Milano, Via E. Bassini 15, Milano, Italy
^{41}
INFN – CNAF, Viale Berti Pichat 6/2, 40127 Bologna, Italy
^{42}
INFN, Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
^{43}
INFN, Sezione di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
^{44}
INFN, Sezione di Milano, Via Celoria 16, Milano, Italy
^{45}
INFN, Sezione di Roma 2, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy
^{46}
Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK
^{47}
Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay, Bât. 121, 91405 Orsay Cedex, France
^{48}
Institut d’Astrophysique de Paris, CNRS (UMR7095), 98bis boulevard Arago, 75014 Paris, France
^{49}
Institute Lorentz, Leiden University, PO Box 9506, 2300 RA Leiden, The Netherlands
^{50}
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{51}
Institute of Theoretical Astrophysics, University of Oslo, Blindern, Oslo, Norway
^{52}
Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, La Laguna, Tenerife, Spain
^{53}
Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, Santander, Spain
^{54}
Istituto Nazionale di Fisica Nucleare, Sezione di Padova, Via Marzolo 8, 35131 Padova, Italy
^{55}
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, USA
^{56}
Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester M13 9PL, UK
^{57}
Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA
^{58}
Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{59}
Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), UTIAS, The University of Tokyo, Chiba 277 8583, Japan
^{60}
Laboratoire de Physique Subatomique et Cosmologie, Université GrenobleAlpes, CNRS/IN2P3, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{61}
Laboratoire de Physique Théorique, Université ParisSud 11 & CNRS, Bâtiment 210, 91405 Orsay, France
^{62}
Low Temperature Laboratory, Department of Applied Physics, Aalto University, Espoo 00076 Aalto, Finland
^{63}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{64}
Mullard Space Science Laboratory, University College London, Surrey RH5 6NT, UK
^{65}
NAOCUKZN Computational Astrophysics Centre (NUCAC), University of KwaZuluNatal, Durban 4000, South Africa
^{66}
National Centre for Nuclear Research, ul. A. Soltana 7, 05400 Otwock, Poland
^{67}
Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, 00716 Warsaw, Poland
^{68}
SISSA, Astrophysics Sector, Via Bonomea 265, 34136 Trieste, Italy
^{69}
San Diego Supercomputer Center, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA
^{70}
School of Chemistry and Physics, University of KwaZuluNatal, Westville Campus, Private Bag X54001, Durban 4000, South Africa
^{71}
School of Physical Sciences, National Institute of Science Education and Research, HBNI, Jatni 752050, Odissa, India
^{72}
School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff CF24 3AA, UK
^{73}
School of Physics and Astronomy, Sun Yatsen University, 2 Daxue Rd, Tangjia, Zhuhai, PR China
^{74}
School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{75}
School of Physics, Indian Institute of Science Education and Research Thiruvananthapuram, Maruthamala PO, Vithura, Thiruvananthapuram 695551, Kerala, India
^{76}
School of Physics, The University of New South Wales, Sydney, NSW 2052, Australia
^{77}
Simon Fraser University, Department of Physics, 8888 University Drive, Burnaby, BC, Canada
^{78}
Sorbonne Université, Observatoire de Paris, Université PSL, École normale supérieure, CNRS, LERMA, 75005 Paris, France
^{79}
Sorbonne UniversitéUPMC, UMR7095, Institut d’Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France
^{80}
Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, Moscow 117997, Russia
^{81}
Space Science Data Center – Agenzia Spaziale Italiana, Via del Politecnico snc, 00133 Roma, Italy
^{82}
Space Sciences Laboratory, University of California, Berkeley, CA, USA
^{83}
The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, 106 91 Stockholm, Sweden
^{84}
UPMC Univ Paris 06, UMR7095, 98bis boulevard Arago, 75014 Paris, France
^{85}
Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex 4, France
^{86}
Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
Received:
13
November
2018
Accepted:
24
January
2019
Analysis of the Planck 2018 data set indicates that the statistical properties of the cosmic microwave background (CMB) temperature anisotropies are in excellent agreement with previous studies using the 2013 and 2015 data releases. In particular, they are consistent with the Gaussian predictions of the ΛCDM cosmological model, yet also confirm the presence of several socalled “anomalies” on large angular scales. The novelty of the current study, however, lies in being a first attempt at a comprehensive analysis of the statistics of the polarization signal over all angular scales, using either maps of the Stokes parameters, Q and U, or the Emode signal derived from these using a new methodology (which we describe in an appendix). Although remarkable progress has been made in reducing the systematic effects that contaminated the 2015 polarization maps on large angular scales, it is still the case that residual systematics (and our ability to simulate them) can limit some tests of nonGaussianity and isotropy. However, a detailed set of null tests applied to the maps indicates that these issues do not dominate the analysis on intermediate and large angular scales (i.e., ℓ ≲ 400). In this regime, no unambiguous detections of cosmological nonGaussianity, or of anomalies corresponding to those seen in temperature, are claimed. Notably, the stacking of CMB polarization signals centred on the positions of temperature hot and cold spots exhibits excellent agreement with the ΛCDM cosmological model, and also gives a clear indication of how Planck provides stateoftheart measurements of CMB temperature and polarization on degree scales.
Key words: cosmology: observations / cosmic background radiation / polarization / methods: data analysis / methods: statistical
© Planck Collaboration 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
This paper, one of a set associated with the 2018 release of data from the Planck^{1} mission (Planck Collaboration I 2020), describes a compendium of studies undertaken to determine the statistical properties of both the temperature and polarization anisotropies of the cosmic microwave background (CMB).
The ΛCDM model explains the structure of the CMB in detail (Planck Collaboration VI 2020), yet it remains entirely appropriate to look for hints of departures from, or tensions with, the standard cosmological model, by examining the statistical properties of the observed radiation. Indeed, in recent years, tantalizing evidence has emerged from the WMAP and Planck fullsky measurements of the CMB temperature fluctuations of the presence of such “anomalies”, and indicating that a modest degree of deviation from global isotropy exists. Such features appear to exert a statistically mild tension against the mainstream cosmological models that themselves invoke the fundamental assumptions of global statistical isotropy and Gaussianity.
A conservative explanation for the temperature anomalies is that they are simply statistical flukes. This is particularly appealing given the generally modest level of significance claimed, and the role of a posteriori choices (also referred to as the “lookelsewhere effect”), i.e., 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. However, determining whether this is the case, or alternatively whether the anomalies are due to real physical features of the cosmological model, cannot be determined by further investigation of the temperature fluctuations on the angular scales of interest, since those data are already cosmicvariance limited.
Polarization fluctuations also have their origin in the primordial gravitational potential, and have long been recognized as providing the possibility to independently study the anomalies found in the temperature data, given that they are largely sourced by different modes. The expectation, then, is that measurements of the fullsky CMB polarization signal have the potential to provide an improvement in significance of the detection of largescale anomalies. Specifically, it is important to determine in more detail whether any anomalies are observed in the CMB polarization maps, and if so, whether they are related to existing features in the CMB temperature field. Conversely, the absence of corresponding features in polarization might imply that the temperature anomalies (if they are not simply statistical excursions) could be due to a secondary effect such as the integrated SachsWolfe (ISW) effect (Planck Collaboration XIX 2014; Planck Collaboration XXI 2016), or alternative scenarios in which the anomalies arise from physical processes that do not correlate with the temperature, e.g., texture or defect models. 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.
In this paper, we present a first comprehensive attempt at assessing the isotropy of the Universe via an analysis of the fullmission Planck fullsky polarization data. Analysis of the 2015 data set in polarization (Planck Collaboration XVI 2016, hereafter PCIS15) was limited on large angular scales by the presence of significant residual systematic artefacts in the High Frequency Instrument (HFI) data (Planck Collaboration VII 2016; Planck Collaboration VIII 2016) that necessitated the highpass filtering of the componentseparated maps. This resulted in the suppression of structure on angular scales larger than approximately 5°. However, the identification, modelling, and removal of previously unexplained systematic effects in the polarization data, in combination with new mapmaking and calibration procedures (Planck Collaboration Int. XLVI 2016; Planck Collaboration III 2020), means that such a procedure is no longer necessary. Nevertheless, our studies remain limited both by the relatively low signaltonoise ratio of the polarization data, and the presence of residual systematic artefacts that can be significant with respect to detector sensitivity and comparable to the cosmological signal. A detailed understanding of the latter, in particular, have a significant impact on our ability to produce simulations that are needed to allow a meaningful assessment of the data. These issues will be subsequently quantified and the impact on results discussed.
The current work covers all relevant aspects related to the phenomenological study of the statistical isotropy and Gaussian nature of the CMB measured by the Planck satellite. Constraints on isotropy or nonGaussianity, as might arise from nonstandard inflationary models, are provided in a companion paper (Planck Collaboration IX 2020). The current paper is organized as follows. Section 2 provides a brief introduction to the study of polarized CMB data. Section 3 summarizes the Planck fullmission data used for the analyses, and important limitations of the polarization maps that are studied. Section 4 describes the characteristics of the simulations that constitute our reference set of Gaussian sky maps representative of the null hypothesis. In Sect. 5 the null hypothesis is tested with a number of standard tests that probe different aspects of nonGaussianity. This includes tests of the statistical nature of the polarization signal observed by Planck using a local analysis of stacked patches of the sky. Several important anomalous features of the CMB sky are studied in Sect. 6, using both temperature and polarization data. Aspects of the CMB fluctuations specifically related to dipolar asymmetry are examined in Sect. 7. Section 8 provides the main conclusions of the paper. Finally, in Appendix A a detailed description is provided of the novel method, called “purified inpainting”, used to generate E and Bmode maps from the Stokes Q and U data.
2. Polarization analysis preamble
Traditionally, the Stokes parameters Q and U are used to describe CMB polarization anisotropies (e.g., Zaldarriaga & Seljak 1997). However, unlike intensity, Q and U are not scalar quantities, but rather components of the rank2 polarization tensor in a specific coordinate basis associated with the map. Such quantities are not rotationally invariant, thus in many analyses it is convenient to consider alternate, but related, polarization quantities.
The polarization amplitude P and polarization angle Ψ, defined as follows,
are commonly used quantities in, for example, Galactic astrophysics. However, completely unbiased estimators of these quantities in the presence of anisotropic and/or correlated noise are difficult to determine (Plaszczynski et al. 2014). Of course, it is still possible to take the observed (noisebiased) quantity and directly compare it to simulations analysed in the same manner. As an alternative, Sect. 5.1 works with the quantity P^{2} and applies a correction for noise bias determined from simulations. A crossestimator based on polarization observations from two maps, P^{2} = Q_{1}Q_{2} + U_{1}U_{2} is also considered.
In addition, a local rotation of the Stokes parameters, resulting in quantities denoted by Q_{r} and U_{r}, is employed in Sects. 5.2 and 5.5. In this case, a local frame is defined with respect to a reference point so that
where ϕ denotes the angle between the axis aligned along a meridian in the local coordinate system centred on the reference point and the great circle connecting this point to a position .
Finally, the rotationally invariant quantities referred to as E and B modes are commonly used for the global analysis of CMB data. 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
where denotes the spinweighted spherical harmonics and are the corresponding harmonic coefficients. If we define
then the invariant quantities are given by
In practice, the Q and U data sets that are analysed are the end products of sophisticated componentseparation approaches. Nevertheless, the presence of residual foregrounds mandates the use of a mask, the application of which during the generation of E and Bmode maps results in E/B mixing (Lewis et al. 2002; Bunn et al. 2003). In Appendix A, we describe the method adopted in this paper to reduce such mixing.
3. Data description
In this paper, we use data from the Planck 2018 fullmission data release (“PR3”) that are made available on the Planck Legacy Archive (PLA^{2}). The raw data are identical to those used in 2015, except that the HFI omits 22 days of observations from the final, thermallyunstable phase of the mission. The release includes sky maps at nine frequencies in temperature, and seven in polarization, provided in HEALPix format (Górski et al. 2005)^{3}, with a pixel size defined by the N_{side} parameter^{4}. For polarization studies, the 353GHz maps are based on polarizationsensitive bolometer (PSB) observations only (see Planck Collaboration III 2020, for details).
Estimates of the instrumental noise contribution and limits on timevarying systematic artefacts can be inferred from maps that are generated by splitting the fullmission data sets in various ways. For LFI, halfring maps are generated from the first and second half of each stable pointing period, consistent with the approach in the 2013 and 2015 Planck papers. For HFI, oddring (O) and evenring (E) maps are constructed using alternate pointing periods, i.e., either odd or even numbered rings, to avoid the correlations observed previously in the halfring data sets. However, for convenience and consistency, we will refer to both of these ringbased splits as “oddeven” (OE), in part as recognition of the signaltonoise ratios of the LFI and HFI maps and their relative contributions to the componentseparated maps described below. Halfmission (HM) maps are generated from a combination of Years 1 and 3, and Years 2 and 4 for LFI, or the first and second half of the fullmission data set in the case of HFI. Note that important information on the level of noise and systematiceffect residuals can be inferred from maps constructed from halfdifferences of the halfmission (HMHD) and oddeven (OEHD) combinations. In particular, the OE differences trace the instrumental noise, but filter away any component fluctuating on timescales longer than the pointing period, whereas the HM differences are sensitive to the time evolution of instrumental effects. A significant number of consistency checks are applied to this set of maps. Full details are provided in two companion papers (Planck Collaboration II 2020; Planck Collaboration III 2020).
As in previous studies, we base our main results on estimates of the CMB from four componentseparation algorithms, – Commander, NILC, SEVEM, and SMICA – as described in Planck Collaboration IV (2020). These provide data sets determined from combinations of the Planck raw frequency maps with minimal Galactic foreground residuals, although some contributions from unresolved extragalactic sources are present in the temperature solutions. Foregroundcleaned versions of the 70, 100, 143, and 217GHz sky maps generated by the SEVEM algorithm, hereafter referred to as SEVEM070, SEVEM100, SEVEM143, and SEVEM217, respectively, allow us to test the frequency dependence of the cosmological signal, either to verify its cosmological origin, or to search for specific frequencydependent effects. In all cases, possible residual emission is then mitigated in the analyses by the use of skycoverage masks.
The CMB temperature maps are derived using all channels, from 30 to 857 GHz, and provided at a common angular resolution of 5′ full width at half maximum (FWHM) and N_{side} = 2048. In contrast to the 2013 and 2015 releases, these do not contain a contribution from the second order temperature quadrupole (Kamionkowski & Knox 2003). An additional window function, applied in the harmonic domain, smoothly truncates power in the maps over the range ℓ_{min} ≤ ℓ ≤ ℓ_{max}, such that the window function is unity at ℓ_{max} = 3400 and zero at ℓ_{max} = 4000. The polarization solutions include information from all channels sensitive to polarization, from 30 to 353 GHz, at the same resolution as the temperature results, but only including contributions from harmonic scales up to ℓ_{max} = 3000. In the context of these CMB maps, we refer to an “oddeven” data split that combines the LFI halfring 1 with the HFI oddring data, and the LFI halfring 2 with the HFI evenring data.
Lowerresolution versions of these data sets are also used in the analyses presented in this paper. The downgrading procedure is as follows. The fullsky maps are decomposed into spherical harmonics at the input HEALPix resolution, these coefficients are then convolved to the new resolution using the appropriate beam and pixel window functions, then the modified coefficients are used to synthesize a map directly at the output HEALPix resolution.
Specific to this paper, we consider polarization maps determined via the method of “purified inpainting” (described in Appendix A) from the componentseparated Q and U data. Figure 1 presents the E and Bmode maps for the four componentseparation methods at a resolution of N_{side} = 128, with the corresponding common masks overplotted. Planck Collaboration IV (2020) notes that some broad largescale features aligned with the Planck scanning strategy are observed in the Q and U data. The detailed impact on E and Bmode map generation is unclear, thus some caution should be exercised in the interpretation of the largest angular scales in the data. Figure 2 does indicate the presence of largescale residuals in the pairwise differences of the componentseparated maps. Finally, we note that the Bmode polarization is strongly noise dominated on all scales, therefore, although shown here for completeness, we do not present a comprehensive statistical analysis of these maps.
Fig. 1. Componentseparated CMB maps at 80′ resolution. Columns show temperature T, and E and Bmode maps, respectively, while rows show results derived with different componentseparation methods. The temperature maps are inpainted within the common mask, but are otherwise identical to those described in Planck Collaboration IV (2020). The E and Bmode maps are derived from the Stokes Q and U maps following the method described in Appendix A. The dark lines indicate the corresponding common masks used for analysis of the maps at this resolution. Monopoles and dipoles have been subtracted from the temperature maps, with parameters fitted to the data after applying the common mask. 
Fig. 2. Pairwise differences between maps from the four CMB componentseparation pipelines, smoothed to 80′ resolution. Columns show temperature, T, and E and Bmode maps, respectively, while rows show results for different pipeline combinations. The grey regions correspond to the appropriate common masks. Monopoles and dipoles have been subtracted from the temperature difference maps, with parameters fitted to the data after applying the common mask. 
In general, we make use of standardized masks made available for temperature and polarization analysis, as described in detail in Planck Collaboration IV (2020). These masks are then downgraded for lowerresolution studies as follows. The binary mask at the starting resolution is first downgraded in the same manner as a temperature map. The resulting 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 again generate a binary mask.
In the case of the data cuts, some additional care must be taken with masking. Since the HFI HM and OE maps contain many unobserved pixels^{5} at a given frequency, some preprocessing is applied to them before the application of the componentseparation algorithms. Specifically, the value of any unobserved pixel is replaced by the value of the corresponding N_{side} = 64 parent pixel. Analysis of the componentseparated maps derived from the data cuts then requires masking of these pixels. However, a simple merge of the unobserved pixel masks at each frequency for a given data cut is likely to be insufficient, since the various convolution and deconvolution processes applied by the componentseparation algorithms will cause leakage of the inpainted values into neighbouring pixels. The masks are therefore extended as follows. Starting with the initial merge of the unobserved pixels over all frequencies, the unobserved pixels are selected and their neighbouring pixels are also masked. This is repeated three times. Lower resolution versions are generated by degrading the binary mask to the target resolution, then setting all pixels with values less than a threshold of 0.95 to zero, while all other pixels have their values set to unity. Masks appropriate for the analysis of the HM and OE maps are generated by combining the unobserved pixel masks with the fullmission standardized masks.
The masks for E and Bmode analysis are extensions of the those applied to the Q and U maps before executing the purified inpainting technique. Specifically, an optimal confidence mask is defined by performing reconstructions on simulated CMBplusnoise realizations, as propagated through all four Planck componentseparation pipelines, then evaluating the residuals with respect to the input fullsky maps. The final mask is specified by requiring that the maximal root mean square (rms) level of the residuals observed in the simulations is less than 0.5 μK, significantly below the cosmological Emode signal. Examples of common masks are shown in Fig. 3.
Fig. 3. Examples of common masks. From top to bottom, the masks correspond to those used for analysing temperature maps, polarization represented by the Stokes Q and U parameters, and Emode polarization data, at a resolution N_{side} = 128. From left to right, fullmission, HM, and OE masks are shown. Note that the masks for E and Bmode analysis are extended relative to those derived for Q and U studies, in order to reduce the reconstruction residuals. 
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. Often, we will simply refer to such a mask as the “common mask”, irrespective of the resolution or data split in question. However, in the latter case, we will refer to fullmission, HM or OE common masks, where appropriate, to avoid confusion. Table 1 lists the N_{side} and FWHM values defining the resolution of these maps, together with the different masks and their sky coverage fractions that accompany the signal maps.
Standardized data sets used in this paper.
4. Simulations
The results presented in this paper are derived using Monte Carlo (MC) simulations. These provide both the reference set of sky maps used for the null tests employed here, and form the basis of any debiasing in the analysis of the real data, as required by certain statistical methods. The simulations include Gaussian CMB signals and instrumental noise realizations that capture important characteristics of the Planck scanning strategy, telescope, detector responses, and datareduction pipeline over the fullmission period. These are extensions of the “full focalplane” simulations described in Planck Collaboration XII (2016), with the latest set being known as “FFP10”.
The fiducial CMB power spectrum corresponds the cosmology described by the parameters in Table 2. Note that the preferred value of τ in Planck Collaboration VI (2020) is slightly lower, at τ = 0.054 ± 0.007. 1000 realizations of the CMB sky are generated including lensing, Rayleigh scattering, and Doppler boosting^{6} effects, the latter two of which are frequencydependent. The signal realizations include the frequencyspecific beam properties of the LFI and HFI data sets implemented by the FEBeCoP (Mitra et al. 2011) beamconvolution approach.
Cosmological parameters for the FFP10 simulations, used to make the simulated maps in this paper, and throughout the Planck 2018 papers.
Given that the instrumental noise properties of the Planck data are complex, we make use of a set of socalled “endtoend” simulations. For HFI, residual systematics must be accounted for in the scientific analysis of the polarized sky signal, thus the simulations include models of all systematic effects, together with noise and sky signal (a fixed CMB plus foregrounds fiducial sky). Realistic timeordered information for all HFI frequencies are then generated and subsequently propagated through the mapmaking algorithm to produce frequency maps. Finally, the sky signal is removed and the resulting maps of noise and residual systematics can be added to the set of CMB realizations. More details can be found in Planck Collaboration Int. XLVI (2016) and Planck Collaboration III (2020). A similar approach is followed by LFI to generate noise MCs that capture important characteristics of the scanning strategy, detector response, and datareduction pipeline over the fullmission period (Planck Collaboration II 2020). A total of 300 realizations are generated at each Planck frequency, for the fullmission, HM and OE data splits. In what follows, we will often refer to simulations of the noise plus systematic effects simply as “noise simulations”. The noise and CMB realizations are then considered to form the FFP10 fullfocal plane simulations.
Finally, the CMB signal and noise simulations are propagated through the various componentseparation pipelines using the same weights as derived from the Planck fullmission data analysis (Planck Collaboration IV 2020). The signal and noise realizations are then permuted to generate 999 simulations^{7} for each componentseparation method to be compared to the data.
In the analyses presented in this paper, we often 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) as represented by the simulations is true. However, this also requires that the simulated reference data set adopts a cosmological model that is sufficiently consistent with that preferred by the data. We have noted above that the τvalue used in the FFP10 simulations is high relative to that preferred by the latest cosmological. As a preliminary assessment, we have considered the predicted variance of the CMB signal for two values of τ, specifically 0.060 (as adopted by the FFP10 simulations), and 0.052 (which is representative of the value determined by an analysis of the HFI data in Planck Collaboration VI 2020) with A_{s} set to an appropriate value. We find that the latter reduces the polarization variance by approximately 20% at N_{side} = 16 and 32. It may be necessary to take this effect into account when interpreting the polarization results in what follows.
Similar considerations apply to the simulated noise and residual systematic effects, particularly given the signaltonoise regime of the polarized data. In order to quantify the agreement of the noise properties and systematic effects in the data and simulations, we use differences computed from various subsets of the fullmission data set. Note that detailed comparisons have been undertaken using the power spectra of the individual frequency maps. Figure 18 of Planck Collaboration II (2020) compares halfring halfdifference (HRHD) spectra for the LFI 30, 44 and 70GHz data with simulations, finding good agreement over most angular scales. Figure 17 of Planck Collaboration III (2020) makes a similar comparison for the HFI 100, 143, 217 and 353GHz halfmission HMHD and OEHD data.
Of more importance to this paper, however, is the consistency of the data and simulations after various componentseparation methods have been applied. As established in Planck Collaboration IV (2020), the corresponding endtoend simulations exhibit biases at the level of several percent with respect to the observations on intermediate and small scales, with reasonable agreement on larger scales. These discrepancies in part originate from the individual frequency bands. For example, the power in the 100–217 GHz HFI simulations underestimates the noise in the data (Planck Collaboration III 2020). Alternatively, biases can arise due to the lack of foreground residuals in the simulations. On small angular scales, the power observed in the temperature data exceeds that of the simulations due to a pointsource residual contribution not included in FFP10. It should, therefore, be apparent that systematic shifts over some ranges of angular scale could contribute to pvalue uncertainties in subsequent studies.
We attempt to verify that the analyses presented in this paper are not sensitive to the differences between the simulations and data. In particular, the comparison of the HMHD and OEHD maps for each componentseparation method with those computed from the ensemble of FFP10 simulations allows us to define the angular scales over which the various statistical tests applied to the data can be considered reliable. These may vary depending on the analysis being undertaken.
5. Tests of nonGaussianity
A key prediction of the standard cosmological model is that an early phase of accelerated expansion, or inflation, gave rise to fluctuations that correspond to a homogeneous and isotropic Gaussian field, and that the corresponding statistical properties were imprinted directly on the primordial CMB (Planck Collaboration XXII 2014; Planck Collaboration XX 2016; Planck Collaboration X 2020). Searching for departures from this scenario is crucial for its validation, yet there is no unique signature of nonGaussianity. Nevertheless, the application of a variety of tests^{8} over a range of angular scales allows us to probe the data for inconsistencies with the theoretically motivated Gaussian statistics.
In previous work (PCIS13; PCIS15), we demonstrated that the Planck temperature anisotropies are indeed consistent with Gaussianity, except for a few apparent anomalies discussed further in the following section. Here, we again apply a nonexhaustive set of tests to the temperature fluctuations in order to confirm previous results, then extend the studies to the polarization data. Of course, significant evidence of deviation from Gaussianity in the statistics of the measured CMB anisotropies is usually considered to be an indicator of the presence of residual foregrounds or systematic artefacts in the data. It is important to be able to mitigate against such possibilities, particularly in the case of polarization anisotropies, where the signaltonoise remains relatively low. The analyses are therefore applied to all four componentseparation 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 FFP10 simulations. The consistency of the results derived from the various componentseparation techniques then provides a strong argument against significant contamination of the data. However, the fidelity of the simulations is limited by the accuracy with which the systematic effects can be modelled; therefore we use HMHD and OEHD null tests to evaluate the agreement of the data and simulations over the scales of interest. It is plausible that the simulations of the polarized signal show evidence of a small level of nonGaussianity depending on the statistical test applied, given the significant level of the systematic effects modelled therein.
5.1. Onedimensional moments
In this section we consider simple tests of Gaussianity based on moments of the CMB temperature and polarization anisotropy maps.
For the temperature analysis, we repeat the study performed in PCIS15 and measure the variance, skewness, and kurtosis of the Planck 2018 componentseparated maps using the unitvariance estimator (Cruz et al. 2011). This method requires a normalized variance sky map, u^{X} defined as:
where X_{i} is the observed temperature at pixel i, is the variance of the CMB signal, and is the variance of the noise for that pixel, estimated using the FFP10 MC simulations. The CMB variance is then determined by finding the value for which the variance of the normalized map u^{X} is unity. The skewness and kurtosis are then subsequently computed from the appropriately normalized map.
In Fig. 4 we show the lowertail probability of the variance, skewness, and kurtosis determined at different resolutions from the four componentseparated maps (left columns) and from the SEVEM frequencycleaned maps (right columns), after applying the appropriate common mask. There is good agreement between the maps, although the NILC results indicate a slightly lower pvalue for the variance at intermediate and high resolutions. This may be related to the small relative power deficit observed between NILC and the other component separation methods over the multipole range ℓ = 100 − 300, as shown in Fig. 15 of Planck Collaboration IV (2020). We note that Planck Collaboration IV (2020) has demonstrated the presence of a noise mismatch between the observed data and simulations, as traced by the HMHD and OEHD maps. However, this is not relevant for analysis of the temperature data, given its very high signaltonoise ratio. The results for 1D moments presented here are in very good agreement with the Planck 2015 analysis (PCIS15), showing a decreasing lowertail probability with decreasing resolution. This lowertail probability is related to the presence of the well known lack of power on large angular scales. However, in the previous analysis we found a minimum value for the probability of 0.5% at N_{side} = 16 for all the maps considered, compared to a probability of roughly 1% here. The difference can be explained by the fact that the 2018 common mask rejects less of the sky than the 2016 common mask, and previous work (PCIS15; Gruppuso et al. 2013) has shown that the low variance anomaly becomes less significant with increasing sky coverage. Indeed, when we apply the 2016 common mask to the current data set, the probability decreases to 0.7–0.8%, in better agreement with the previous results. The skewness and kurtosis results do not show any anomalous behaviour, in agreement with earlier analyses.
Fig. 4. Lowertail probabilities of the variance (top), skewness (centre), and kurtosis (bottom), determined from the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) componentseparated temperature maps (left) and the SEVEM070 (light green), SEVEM100 (dark blue), SEVEM143 (yellow), and SEVEM217 (magenta) frequencycleaned maps (right) at different resolutions. 
In polarization we follow a different approach, as a consequence of the lower signaltonoise ratio. Specifically, we subtract the noise contribution to the total variance of the polarization maps and define the estimator
where Q and U are the Stokes parameters of the observed polarization maps, and are noise estimates determined from MC simulations. Planck Collaboration IV (2020) indicates a mismatch between the noise in the data and that in simulations for map resolutions above N_{side} = 256. This corresponds to a few percent of the theoretical CMB variance up to N_{side} = 1024, while it is much larger at the highest resolution. Since the noise mismatch is likely to affect the less signaldominated polarization results, we also define a crossvariance estimator that determines the variance from the two maps available for each data split, HM or OE, respectively:
where Q_{1}, Q_{2}, U_{1}, and U_{2} are the Stokes parameters of the two maps from either the HM or OEcleaned data split, and is the corresponding noise contribution to the total variance in polarization estimated from the corresponding simulations. Note that a crossestimator should be less affected by noise mismatch, although correlated noise remains an issue. However, it is impossible to assess if the latter is well described by the simulations.
In Fig. 5 we show the expected signaltonoise ratio of the polarization variance for the componentseparated maps, determined by comparing the theoretical variance of the signal at different resolutions (as evaluated from the Planck FFP10 fiducial model, including beam and pixel window function effects) to the corresponding MC estimate of the noise, . All of the methods show similar behaviour, with a maximum signaltonoise ratio of about 0.8 on intermediate scales, N_{side} = 512. The minimum ratio is observed at N_{side} = 64, as explained by the fact that the EE angular power spectrum exhibits a low amplitude over the multipole range ℓ = 10 − 100. At very large scales, N_{side} = 16, the signaltonoise ratio increases again, but with an amplitude that depends noticeably on the componentseparation method considered. At very high resolutions the signaltonoise ratio drops, as expected.
Fig. 5. Signaltonoise ratio for the variance estimator in polarization for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), obtained by comparing the theoretical variance from the Planck FFP10 fiducial model with an MC noise estimate (righthand term of Eq. (7)). Note that the same colour scheme for distinguishing the four componentseparation maps is used throughout this paper. 
In Fig. 6 we show the lowertail probabilities of the variance determined from the fullmission and the HM and OE data splits using the appropriate common mask, compared to the corresponding results from MC simulations. At high resolutions, the lowertail probability determined from the variance of the fullmission data approaches zero. As previously noted, this is due to the poor agreement between the noise properties of the data and the MC simulations, in particular at high resolution. This explanation is further supported by the fact that the lowertail probability becomes more compatible with the MC simulations when we consider the crossvariance analyses. However, given the uncertainties in the properties of the correlated noise in the simulations, we prefer to focus on the intermediate and large angular scales, N_{side} ≤ 256. We note that there is a trend towards lower probabilities as the resolution decreases from N_{side} = 256, similar to what is observed with the temperature data. This behaviour is common to all of the componentseparated methods and also to the SEVEM143 frequencycleaned data, although with different probabilities at a given resolution. The SEVEM070 frequencycleaned map is not compatible with the MC simulations for resolutions lower than N_{side} = 128. This may be due to the presence of either residual foregrounds in the data or systematic effects that are not sufficiently well represented by the MC simulations. Although the compatibility of the data with the MC simulations is generally adequate, the large variation of probabilities seen for different componentseparation methods and resolutions, even when a crossestimator is considered, suggests that correlated noise and residual systematics in the data may not be sufficiently well described by the current set of simulations.
Fig. 6. Lowertail probabilities of the variance determined from the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) componentseparated polarization maps (left) and the SEVEM070 (light green), SEVEM100 (dark blue), SEVEM143 (yellow), and SEVEM217 (magenta) frequencycleaned maps (right) at different resolutions. The top, middle and bottom rows correspond to results evaluated with the fullmission, HM and OEcrossvariance estimates, respectively. In this figure, small pvalues would correspond to anomalously low variance. 
In an attempt to minimize the impact of correlated noise, we consider the crossvariance estimated between pairs of frequencies from the SEVEM frequencycleaned maps, using fullmission, HM and OE data sets. The results are shown in Fig. 7. Note that the combination of the SEVEM070 data with higher frequency maps is consistent with the MC simulations, supporting the idea that residual systematic effects in the former, which are not well described by the corresponding simulations, bias results computed only with the 70GHz cleaned data. In addition, the crossvariance determined between the SEVEM070 and SEVEM217 maps yields a particular low probability. Since the 217GHz data are used to clean the 70GHz map, it is probable that this particular combination is more affected by correlated residuals than elsewhere. However, the effect disappears when considering the HM or OE crossvariance data.
Fig. 7. Lowertail probabilities of the crossvariance determined between the SEVEM frequencycleaned polarization maps (top) or from the HM (centre) or OE (bottom) SEVEM frequencycleaned maps at different resolutions. In this figure, small pvalues would correspond to anomalously low variance. 
In summary, we confirm previous results based on the analysis of the temperature anisotropy (PCIS13; PCIS15), indicating that the data are consistent with Gaussianity, although exhibiting low variance on large angular scales, with a probability of about 1% as compared to our fiducial cosmological model. In polarization we find reasonable consistency with MC simulations on intermediate and large angular scales, but there is a considerable range of pvalues found, depending on the specific combinations of data considered. This indicates that the lower signaltonoise ratio of the Planck data in polarization, and, more specifically, the uncertainties in our detailed understanding of the noise characterization (both in terms of amplitude and correlations between angular scales) limits our ability to pursue further investigate the possible presence of anomalies in the 1D moments.
5.2. Npoint correlation functions
In this section, we present tests of the nonGaussianity of the Planck 2018 temperature and polarization CMB data using realspace Npoint correlation functions.
An Npoint correlation function is defined as the average product of N observables, measured in a fixed relative orientation on the sky,
where the unit vectors span an Npoint polygon. If statistical isotropy is assumed, these functions do not depend on the specific position or orientation of the Npoint polygon on the sky, but only on its shape and size. In the case of the CMB, the fields, X, correspond to the temperature, T, and the two Stokes parameters, Q and U, which describe the linearly polarized radiation in direction . Following the standard CMB convention, Q and U are defined with respect to the local meridian of the spherical coordinate system of choice. To obtain coordinatesystemindependent Npoint correlation functions, we define Stokes parameters in a radial system, denoted by Q_{r} and U_{r}, according to Eq. (2), where the reference point, , is specified by the centre of mass of the polygon (Gjerløw et al. 2010). In the case of the 2point function, this corresponds to defining a local coordinate system in which the local meridian passes through the two points of interest (see Kamionkowski et al. 1997).
The correlation functions are estimated by simple product averages over all sets of N pixels fulfilling the geometric requirements set by the 2N − 3 parameters θ_{1}, …, θ_{2N − 3} characterizing the shape and size of the polygon,
Here, pixel weights represent masking and are set to 1 or 0 for included or excluded pixels, respectively.
The shapes of the polygons selected for the analysis are not more optimal for testing Gaussianity than other configurations, but are chosen because of ease of implementation and for comparison of the results with those for the 2013 and 2015 Planck data sets. In particular, we consider the 2point function, as well as the pseudocollapsed and equilateral configurations for the 3point function. Following Eriksen et al. (2005), the pseudocollapsed configuration corresponds to an (approximately) isosceles triangle, where the length of the baseline falls within the second bin of the separation angles and the length of the longer edge of the triangle, θ, parametrizes its size. Analogously, in the case of the equilateral triangle, the size of the polygon is parametrized by the length of the edge, θ.
We use a simple χ^{2} statistic to quantify the agreement between the observed data and simulations. This is defined by
Here, is the difference between the observed, , and the corresponding average from the MC simulation ensemble, ⟨C_{N}(θ_{i})⟩, of the Npoint correlation function for the bin with separation angle θ_{i}, normalized by the standard deviation of the difference, σ_{N}(θ_{i}), and N_{bin} is the number of bins used for the analysis. If is the kth simulated Npoint correlation function difference and N_{sim} is the number of simulations, then the covariance matrix (normalized to unit variance) M_{ij} is estimated by
where . However, due to degeneracies in the covariance matrix resulting from an overdetermined system and a precision in estimation of the matrix elements of order , the inversion of the matrix is unstable. To avoid this, a singularvalue decomposition (SVD) of the matrix is performed, and only those modes that have singular values larger than are used in the computation of the χ^{2} statistic (Gaztañaga & Scoccimarro 2005). We note that this is a modification of the procedure used in previous Planck analyses (PCIS13; PCIS15). Finally, we also correct for bias in the inverse covariance matrix by multiplying it by a factor (Hartlap et al. 2007).
We analyse the CMB estimates at a resolution of N_{side} = 64 due to computational limitations. The results for the 2point correlation functions of the CMB maps are presented in Fig. 8, while in Fig. 9 the 3point functions for the Commander maps are shown. In the figures, the Npoint functions for the data are compared with the mean values estimated from the FFP10 MC simulations. Note that the mean behaviour of the 3point functions derived from the simulations indicates the presence of small nonGaussian contributions, presumably associated with modelled systematic effects that are included in the simulations. Furthermore, both the mean and associated confidence regions vary between componentseparation methods, which reflects the different weightings given to the individual frequency maps that contribute to the CMB estimates, and the systematic residuals contained therein. Some evidence for this behaviour can also be found in the analysis of HMHD and OEHD maps in the companion paper Planck Collaboration IV (2020). To avoid biases, it is essential to compare the statistical properties of a given map with the associated simulations. Comparing with simulations without systematic effects could lead to incorrect conclusions.
Fig. 8. 2point correlation functions determined from the N_{side} = 64 Planck CMB 2018 temperature and polarization maps. Results are shown for the Commander, NILC, SEVEM, and SMICA maps (first, second, third, and fourth rows, respectively). The solid lines correspond to the data, while the black three dotsdashed lines indicate the mean determined from the corresponding FFP10 simulations, and the shaded dark and light grey areas indicate the corresponding 68% and 95% confidence regions, respectively. 
Fig. 9. 3point correlation functions determined from the N_{side} = 64 Planck CMB Commander 2018 temperature and polarization maps. Results are shown for the pseudocollapsed 3point (upper panel) and equilateral 3point (lower panel) functions. The red solid line corresponds to the data, while the black three dotsdashed line indicates the mean determined from the FFP10 Commander simulations, and the shaded dark and light grey regions indicate the corresponding 68% and 95% confidence areas, respectively. See Sect. 5.2 for the definition of the separation angle θ. 
The probabilities of obtaining values of the χ^{2} statistic for the Planck fiducial ΛCDM model at least as large as the observed temperature and polarization values are provided in Table 3. It is worth noting that the values of the Npoint functions and their associated errors are strongly correlated between different angular separations. The estimated probabilities, which take into account such correlations, therefore provide more reliable information on the goodnessoffit between the data and the simulations than a simple inspection of the figures can reveal.
Probabilities of obtaining values for the χ^{2} statistic of the Npoint functions determined from the Planck fiducial ΛCDM model at least as large as those obtained from the Commander, NILC, SEVEM, and SMICA temperature and polarization (Q and U) maps at N_{side} = 64 resolution.
The Npoint function results show excellent consistency between the CMB temperature maps estimated using the different componentseparation methods. Some differences between results for the 2015 and 2018 temperature data sets are caused by the use of different masks in the analysis, and the adoption of the pseudoinverse matrix in the computation of the χ^{2} statistic, as described by Eq. (11). In the case of polarization, some scatter is observed between the functions computed for different methods, which is a consequence of the relatively low signaltonoise ratio of the polarized data on large angular scales. Interestingly, a tendency towards very high probability values is observed for the pseudocollapsed TTQ_{r} 3point functions for all methods, and for the equilateral TQ_{r}Q_{r} functions in the case of Commander and SEVEM.
As an alternative to the Stokes parameters, we also consider Npoint functions computed from the temperature and Emode polarizations maps. The probabilities of obtaining values of the χ^{2} statistic for the Planck fiducial ΛCDM model at least as large as the observed temperature and polarization values are provided in Table 4. Here, we see that the most significant deviations between the data and the simulations occur for the TTE 3point functions for all componentseparation methods.
Probabilities of obtaining values for the χ^{2} statistic of the Npoint functions determined from the Planck fiducial ΛCDM model at least as large as the those obtained from the Commander, NILC, SEVEM, and SMICA temperature and polarization (Emode) maps at N_{side} = 64 resolution.
Nevertheless, we conclude that no strong evidence is found for statistically significant deviations from Gaussianity of the CMB temperature and polarization maps using Npoint correlation functions.
Finally, we note that the results for the TT correlation function confirm the lack of structure at large separation angles, noted in the WMAP firstyear data by Bennett et al. (2003) and in previous Planck analyses (PCIS13; PCIS15). We will discuss this issue further in Sect. 6.1, where we also consider the behaviour of the TQ_{r} correlation function.
5.3. Minkowski functionals
In this section, we present a morphological analysis of the Planck 2018 temperature and polarization CMB maps using Minkowski functionals. The Minkowski functionals (hereafter MFs) describe the morphology of fields in any dimension and have long been used to investigate nonGaussianity and anisotropy in the CMB (see Planck Collaboration XXIII 2014, and references therein). They are additive for disjoint regions of the sky and invariant under rotations and translations. For the polarization data, we analyse the scalar Emode representation, since the MFs computed from the spin2 Q and U Stokes parameters are no longer invariant under rotation after the application of a mask (Chingangbam et al. 2017).
We compute MFs for the regions colder and hotter than a given threshold ν, usually defined in units of the sky rms amplitude, σ_{0}. The three MFs, namely the area V_{0}(ν) = A(ν), the perimeter V_{1}(ν) = C(ν), and the genus V_{2}(ν) = G(ν), are defined respectively as
where N_{ν} is the number of pixels with Δ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 or cold spot. There are two approaches to the calculation of σ_{0}. The first possibility is to use a population rms, which can be inferred from the average variance of the simulations. Using this estimator provides robust results for low resolutions. An alternative is to use the sample rms, estimated directly from the map in question. Cammarota & Marinucci (2016) have shown that this approach increases the sensitivity of MFbased tests, and thus we adopt this definition of σ_{0} in our analysis.
Furthermore, the MFs can be written as a product of a function A_{k} (k = 0, 1, 2), which depends only on the Gaussian power spectrum, and v_{k}, which is a function only of the threshold ν (see e.g., Vanmarcke et al. 1983; Pogosyan et al. 2009; Gay et al. 2012; Matsubara 2010; Fantaye et al. 2015). This factorization is valid in the weakly nonGaussian case. In this paper, we use the normalized MFs, v_{k}, to focus on deviations from Gaussianity, with reduced sensitivity to the cosmic variance of the Gaussian power spectrum. However, we have verified that the results derived using both normalized or unnormalized MFs are consistent in every configuration^{9}. The analytical expressions are
with H_{n}, the Hermite function,
The amplitude A_{k} depends only on the shape of the power spectrum C_{ℓ} through the parameters σ_{0} and σ_{1}, the rms of the field and its first derivative, respectively:
with ω_{k} ≡ π^{k/2}/Γ(k/2 + 1).
In order to characterize the MFs, we consider two approaches for the scaledependent analysis of the temperature and polarization sky maps: in real space via a standard Gaussian smoothing and degradation of the maps; and in harmonic space by using needlets. Such a complete investigation should provide insight regarding the harmonic and spatial nature of possible nonGaussian features detected with the MFs.
First, we undertake a realspace analysis by computing the three normalized functionals described above at different resolutions and smoothing scales for each of the four componentseparation methods. The appropriate common mask is applied for a given scale. The MFs are evaluated for 12 thresholds ranging between −3 and 3 in σ_{0} units, providing a total of 36 different statistics y = {v_{0}, v_{1}, v_{2}}. A χ^{2} value is then computed by combining these, assuming a Gaussian likelihood for the MFs at every threshold, taking into account their correlations (Ducout et al. 2012) using a covariance matrix computed from the FFP10 :
where is the mean of the statistics y computed on the simulations, i, j are the threshold indices from the combined MFs, and
is the covariance matrix estimated from the FFP10 simulations. The covariance matrix is well converged for this low number of statistics (i.e., 36).
The results for temperature and Emode polarization data are shown in Figs. 10 and 11, respectively. The first three columns of panels in these figures show the normalized MFs together with their varianceweighted difference with respect to the mean of the simulations for the three MFs. The rightmost column of panels in Figs. 10 and 11 presents the χ^{2} obtained when the three MFs are combined with an appropriate covariance matrix derived using the FFP10 simulations. The vertical lines in these figures represent the data, with different colours for the different componentseparation methods. The grey shaded regions in the MFs plot and the histogram in the χ^{2} plot are determined from the FFP10 simulations. Table 5 presents the corresponding pvalues determined for the different componentseparation techniques and map resolutions, between N_{side} = 16 and N_{side} = 2048 for temperature, and between N_{side} = 16 and N_{side} = 1024 for the polarization Emode.
Fig. 10. Realspace normalized MFs determined from the Planck 2018 temperature data using the four componentseparated maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The grey region corresponds to the 99th percentile area, estimated from the FFP10 simulations processed by the SMICA method, while the dashed curves with matching colours outline the same interval for the other componentseparation methods. Results are shown for analyses at N_{side} = 32, 256, and 1024. The rightmost column shows the χ^{2} obtained by combining the three MFs in real space with an appropriate covariance matrix derived from FFP10 simulations. The vertical lines correspond to values from the Planck data. 
Fig. 11. As in Fig. 10, but for the Planck 2018 Emode polarization data. The Planck data are consistent with the Gaussian FFP10 simulations, but variations between the different componentseparation methods are evident. 
Probability as a function of resolution determined using normalized realspace MFs for the temperature and polarization Emode data.
For the temperature results, the χ^{2} values computed for the different componentseparation methods are more consistent than was the case for the Planck 2015 analysis, for all scales. In the case of the Emode results, we find no significant discrepancy between the Planck data and the FFP10 simulations. The striking variation in the pvalues for the four componentseparation methods, is also observed when considering individual realizations in the set of simulations.
As a complement to the pixelbased analysis, we also determine the MFs of needlet coefficient maps on various scales (see Table 6). Measuring the MFs in needlet space, as compared to the usual pixelspace case, has two clear advantages: the needlet maps are minimally affected by masked regions due to the localization of the needlet filter in pixel space, especially at highfrequency; and the doublelocalization properties of needlets (in real and harmonic space) allow a much more precise, scalebyscale, interpretation of any possible anomalies. While the behaviour of standard allscale (pixelbased) 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. This is because if one changes the shape of the power spectrum while still keeping ∑_{ℓ}[2ℓ+1]C_{ℓ} constant, the pixelspace MFs will not change but the needlet MFs are affected. This sensitivity to the shape of the power spectrum can be used to understand in more detail the nature of any possible nonGaussianity detected by the MF analysis.
Probability as a function of needlet scale.
The needlet components of a scalar CMB field are defined by Marinucci et al. (2008) and Baldi et al. (2009), and are given by
where j on the lefthand side is the needlet index and j on the righthand side is a power. Here, denotes the component at multipole ℓ of the CMB map (corresponding to temperature or the polarization E or B modes), i.e.,
where denotes the pointing direction, 𝔹 is a fixed parameter that controls the needlet’s band width (usually taken to be between 1 and 3), and b(.) is a smooth function such that ∑_{j}b^{2}(ℓ/𝔹^{j}) = 1 for all ℓ. In Fantaye et al. (2015), it is shown that a general analytical expression for MFs at a given needlet scale j can be written as
where t_{0} = 2, t_{1} = 0, and t_{2} = 4π, are the EulerPoincaré characteristic, boundary length, and area of the full sphere, respectively. The quantities v_{k} are the normalized MFs given in Eq. (17), while the needletscale amplitudes have a similar form to A_{k}, but with the variances of the map and its first derivative given by
We adopted the needlet parameters 𝔹 = 2, j = 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 for this analysis. Note that the jth needlet scale has compact support over the multipole range [2^{j − 1}, 2^{j + 1}]. For clarity in all the figures, we refer to the different needlet scales by their central multipole ℓ_{c} = 2^{j}.
To obtain the needlet maps at different scales, we initially decompose into spherical harmonics the temperature and polarization maps, inpainted using diffusive and purified inpainting (see Appendix A), respectively, at N_{side} = 1024. In all cases, we set the maximum multipole to ℓ_{max} = 2048, which is twice the maximum resolution considered for the needlet MF analysis. We then obtain the jth needletscale map, by computing Eq. (24) using the HEALPix map2alm routine at the appropriate N_{side}. Specifically, we use N_{side} = 16 for j = 1, 2, 3, 4, and N_{side} = 2^{j} for the remaining needlet scales. These choices allow us to adopt the same masks used for the pixelspace analysis without alteration.
Once the needlet maps are obtained, we follow the identical procedure as in the pixelspace case to compute the three MFs. The results derived from the Commander, NILC, SEVEM, and SMICA componentseparated temperature maps for needlet scales 𝔹 = 2, j = 2, 6, 9 are shown in Fig. 12. As can be seen from both the MF and χ^{2} plots, the Planck 2018 temperature data are consistent with the Gaussian FFP10 simulations. This is true for all the needlet scales considered. Similar to the pixelspace MFs case, the results indicate a high degree of consistency among the four componentseparation methods.
Fig. 12. Needletspace normalized MFs of the Planck 2018 temperature data using the four componentseparated maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The grey region corresponds to the 99th percentile area, estimated from the FFP10 simulations processed by the SMICA method, while the dashed curves with matching colours outline the same interval for the other componentseparation methods. The needlet MFs are denoted by the central multipole of the needlet filter ℓ_{c} = 2^{j}; the jth needlet parameter has compact support over the multipole range [2^{j − 1}, 2^{j + 1}]. The rightmost column shows the χ^{2} obtained by combining the three MFs in needlet space with an appropriate covariance matrix derived from FFP10 simulations. The vertical lines correspond to values from the Planck data. 
The results for the Emode polarization data are shown in Fig. 13. Compared to the temperature results, the different componentseparation methods show greater variation, although no significant deviations between data and simulations are observed. A similar degree of scatter in the results determined for the OEHD and HMHD componentseparated maps suggests that noise can play a significant role in explaining this variation.
Fig. 13. As in Fig. 12, but for the Planck 2018 Emode polarization data. The Planck data are consistent with the Gaussian FFP10 simulations, but some variations between the different componentseparation methods are evident. 
In summary, both the pixel and needletspace MF analyses show that the 2018 Planck temperature and polarization maps are consistent with the Gaussian simulations over angular scales corresponding to ℓ_{max} = 2048.
5.4. Peak statistics
In this section, we present nonGaussianity tests of the Planck 2018 temperature and polarization data using the statistical properties of local extrema (both minima and maxima, to be referred to collectively as “peaks”) determined from the componentseparated CMB maps. The peaks, defined as pixels whose amplitudes are either higher or lower than the corresponding values for all of their nearest neighbours, compress the information contained in the map and provide tests complementary to sphericalharmonicsbased methods. Peak statistics are particularly sensitive to nonGaussian features localized in real space.
The statistical properties of peaks for an isotropic Gaussian random field were derived in Bond & Efstathiou (1987). In particular, the fraction of peaks with amplitudes x above a certain threshold x/σ > ν is given by
where σ is the rms random field amplitude, and γ is the shape parameter, dependent on the spectrum of the Gaussian random field. Peak locations and amplitudes, and various derived quantities, such as their correlation functions, have previously been used to characterize the WMAP maps in Larson & Wandelt (2004, 2005) and Hou et al. (2009), and Planck data in PCIS15.
We consider peak statistics from the Planck componentseparated temperature and polarization maps at N_{side} = 1024, with Emode maps reconstructed by the purified inpainting method described in Appendix A. The maps are prewhitened by convolving them with an isotropic function derived from the isotropic bestfit CMB power spectrum, combined with a diagonal approximation to the instrumental noise covariance. Then, a confidence mask is applied, and weighted convolution is performed with a 2DGaussian smoothing kernel (that we label as “GAUSS”), as described in Appendix A of PCIS15. The mask is further extended by rejecting pixels with an effective convolution weight that differs from unity by more than 12%, and peaks within it are extracted and analysed. The empirical cumulative densityfunction (CDF) of peak values x, defined for a set of n peaks at values {X_{i}},
is generated by sorting the peak values {X_{i}} extracted from the map in ascending order, and comparison to the median CDF derived from an identical analysis of the simulations. A statistical measure of the difference of the two distributions is provided by the Kolmogorov–Smirnov deviation
Although the Kolmogorov–Smirnov deviation has a known limiting distribution, to evaluate pvalues we derive its CDF directly from the simulations.
The peak distributions for T and Emode peaks of the SMICA CMB map filtered at two different scales (120′ and 600′ FWHM) are shown in Fig. 14. The lower panels show empirical peak CDFs and total peak counts, compared to the Gaussian random field peak CDFs derived from Eq. (29) by fitting parameters σ and γ to the median CDF from 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. Other component separation methods produce similar results, and no significant deviations from Gaussian expectations are observed in the polarization peak statistics.
Fig. 14. Cumulative densityfunction of the peak distributions for the SMICA temperature T (left) and reconstructed Emode polarization (right) maps. Top row: peak CDF filtered with a GAUSS kernel of 120′ FWHM, bottom row: peak CDF filtered with the same kernel of 600′ FWHM. The spectral shape parameter γ (see Eq. (29)) is the bestfit value for the simulated ensemble, as indicated by the cyan circle in Fig. 15. No significant deviations from Gaussian expectations are observed. Similar results are obtained for other componentseparation methods. 
A further statistical test of isotropic Gaussian random field expectations is to check if the bestfit parameters, σ and γ, to the observed empirical peak CDF, agree with those derived from individual simulated realizations. The distribution of bestfit values of σ and γ from simulations is compared to the observed value in Fig. 15 for the same data as presented in Fig. 14. Once again, the polarization results are consistent with Gaussian expectations.
Fig. 15. Distribution of bestfit Gaussian peak CDF spectral shape parameters, σ and γ (as defined in Eq. (29)), recovered from FFP10 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) for the SMICA temperature T (left) and the reconstructed Emode polarization (right) maps. 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. Upper panel: peak CDF parameters for the SMICA map filtered with a GAUSS kernel of 120′ FWHM, lower panel: corresponding peak CDF using the same kernel with 600′ FWHM. Similar results are obtained for the other componentseparation methods. 
Extending the analysis of Larson & Wandelt (2004) to polarization data, we also evaluate whether the distributions of maxima and minima are separately consistent with simulations. Counts of maxima and minima in the filtered maps are compared to the distributions determined from simulations in Fig. 16. 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. 17. The observed minima/maxima counts and means are not significantly different from the fiducial model.
Fig. 16. Cumulative densityfunction of the number of all extrema, maxima (red) and minima (blue), derived from simulations, compared to the equivalent values observed for the SMICA temperature T (left) and the reconstructed Emode polarization (right). Upper panel: peak counts for maps filtered with a GAUSS kernel of 120′ FWHM. Lower panel: corresponding peak count CDF for the same kernel of 600′ FWHM. Similar results are obtained for the other componentseparation methods. 
Fig. 17. Cumulative densityfunction of the mean amplitude of all extrema, maxima (red) and minima (blue), derived from simulations, compared to the equivalent values observed for the SMICA temperature T (left) and the reconstructed Emode polarization (right). Upper panel: peak mean amplitudes for maps filtered with a GAUSS kernel of 120′ FWHM. Lower panel: corresponding peak CDF for the same kernel of 600′ FWHM. Similar results are obtained for the other component separation methods. The amplitude values are shown in arbitrary (dimensionless) units determined by map prewhitening. 
To summarize, the temperature and Emode peak statistics determined from the Planck componentseparated maps show no significant anomalies, except perhaps for the previously known Cold Spot discussed in Sect. 6.5. They are consistent with the predictions of the ΛCDM model, and our understanding of the instrument and noise properties of the 2018 Planck data.
5.5. Stacking of CMB peaks
5.5.1. Nonoriented stacking
The stacking of CMB anisotropies in both temperature and polarization around the locations of extrema generates characteristic patterns that connect to the physics of recombination and anisotropy power spectra, as discussed in detail in PCIS15 and MarcosCaballero et al. (2016). Comparison of the results from the Planck CMB maps with the predictions of the Planck fiducial ΛCDM model acts as both a test of their consistency, and an assessment of the quality of the data at the map level. Furthermore, the stacking procedure is expected to mitigate the impact of smallscale noise and residual systematic effects, thus minimizing the impact of inconsistencies in these properties between the data and simulations.
Hot (or cold) peaks are selected in the CMB intensity map as local maxima (or minima) by comparison with their nearest neighbour pixels, and grouped into different ranges above (below) a given threshold ν (in rms units of the intensity map). In order to facilitate the comparison of the polarization signal between different peaks, we use transformed Stokes parameters, Q_{r} and U_{r}, defined by Eq. (2), where the reference point, , is specified by the centre of each extremum . The Q_{r} component then traces the linear polarization in terms of radial (Q_{r} > 0) and tangential (Q_{r} < 0) contributions with respect to the centre of the peak. This stacking method is referred to as “nonoriented”, because the orientation is defined relative to the local meridian rather than any property of the data themselves.
Figure 18 compares the patterns seen in the Planck and WMAP data when averaging over patches centred on CMB intensity maxima above ν = 0 and 3. To enhance the visualization, a random rotation of the patch is performed around each maximum before stacking (in particular, this allows a residual pattern in U_{r} due to pixelization effects to be removed). Specifically, we consider the SEVEM foregroundcleaned map at a HEALPix pixel resolution of N_{side} = 1024 convolved with a Gaussian beam of 10′, and a noiseweighted combination of the WMAP V and W bands at a HEALPix pixel resolution of N_{side} = 512 convolved with a Gaussian beam of 30′. Note that when the signaltonoise is sufficient, the stacking procedure tends to provide an image with azimuthal symmetry about its centre, due to the almost uncorrelated orientations of the temperature peaks. Indeed, the Planck analysis for maxima above ν = 3 clearly reveals two rings, while the WMAP data are noise dominated at the same resolution. Furthermore, when maxima are selected above ν = 0, the enhanced resolution of the Planck data allows additional inner rings to be observed compared to WMAP.
Fig. 18. Upper panels: patches of T, Q_{r}, and U_{r} in microKelvin for maxima above ν = 0 from the WMAP V+W data at a HEALPix resolution of N_{side} = 512 convolved with a Gaussian beam of 30′ (first row), and from the SEVEM data at a HEALPix resolution of N_{side} = 1024 convolved with a Gaussian beam of 10′ (second row). Lower panels: the same thing, but only for maxima above ν = 3. The polarization direction is represented as a headless vector. Note that, in the bottom panel, the length of the headless vectors is divided by 3 with respect to the scale used in the upper panel. 
We now consider the consistency of the Planck nonoriented results with the predictions of ΛCDM by focussing on the mean value of the angular profiles μ(θ) estimated as the average of the profiles around all hot (cold) peaks above (below) a certain threshold ν. Although the analysis is performed on data at a resolution N_{side} = 1024, the profiles are only sampled with 16 bins to ensure that the covariance matrix can be well estimated from the simulations.
A χ^{2} estimator is used to quantify the differences between the μ(θ) profiles obtained from the data and the expected values estimated with simulations:
with the covariance matrix defined as
where the index k denotes a simulation, N is the total number of simulations used to estimate this matrix, and is the ensemble average.
Figure 19 presents the results for maxima and minima selected at thresholds of ν = 0 and 3 for the CMB maps provided by the four componentseparation pipelines, compared to the simulations. The corresponding pvalues for the comparison are presented in Table 7. We see that all componentseparation methods yield consistent results. No significant differences for the intensity profiles μ_{T}(θ) are observed with respect to the results found in the Planck2015 analysis. As in the latter case, a systematic deviation between the data and the mean value of simulations is present. This was previously interpreted in PCIS15 as an effect connected with the deficit in the observed power spectrum at low multipoles.
Fig. 19. Mean radial profiles of T, Q_{r}, and U_{r} in microkelvin obtained for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) at N_{side} = 1024. The top panel in each plot shows the radial profile compared to the simulation average (grey line), while the lower panel shows the difference of the radial profile and simulation average on smaller scales. Results based on stacks around temperature maxima and minima are shown in the upper and lower rows, respectively. The left three columns present results for peaks selected above the null threshold, while the right three columns show the equivalent results for peak amplitudes above (maxima) or below (minima) 3 times the dispersion of the temperature map. The black dots (connected by dashed lines) show the mean value from simulations and the shaded regions correspond to the ±1σ (68%) and ±2σ (95%) error bars estimated from SEVEM simulations. Note that the “Diff” curves for each componentseparation method are computed using the corresponding set of ensemble averages, although only the ensemble average from SEVEM is shown here. 
Fraction of simulations with higher values of χ^{2} than the observed ones for the T, Q_{r}, and U_{r} angular profiles, computed from the stacking of hot and cold extrema selected above or below the ν = 0 and ν = 3 thresholds.
PCIS15 also demonstrated that the χ^{2} statistic is suboptimal when considering the systematic shift between data and simulations, as seen in the intensity profiles μ_{T}, and may lead to misleading pvalues. We therefore consider an alternative quantity, the integrated profile deviation Δμ_{T}(W), to evaluate the consistency between the data and the model. This is defined as
where R represents the size of stacking patches (3° in this case), and the weighting function W is chosen to be proportional to the expected profile. The pvalues obtained in this case are given in Table 8, and are consistent with the deviations shown in Fig. 19.
Fraction of simulations with higher values of χ^{2} than the observed ones for Δμ_{T}, computed from the stacking of hot and cold spots selected above the ν = 0 and ν = 3 thresholds.
Turning to polarization, while Table 7 generally shows consistent results between the data and simulations, somewhat low pvalues for the minima below ν = 0 are observed. However, as is apparent from Fig. 19, the deviation between the data and the mean value of simulations at this threshold is less significant for the maxima. We note that no evidence of asymmetry is found between the number of maxima and minima in the data when compared to simulations. In addition, the sum of the profiles from maxima and minima is consistent with zero for the data with respect to the simulations. Since similar pvalues are found when comparing the mean angular profiles from the HMHD and OEHD data splits to corresponding simulations, we consider that it is unlikely that this discrepancy for the minima is cosmological in origin.
We find that the noise level traced by the U_{r} profile seems to be slightly higher than was seen in PCIS15. In addition, for ν = 0 and low values of θ(< 0.° 5), the mean determined from the simulations does not tend to a null value, especially for the NILC and SMICA maps (although this is not observed in the HMHD and OEHD analyses). Both effects may be related to noise correlations introduced by changes in the raw data processing, and variations in the weights ascribed to the frequency maps by the different componentseparation methods. In addition, it was shown in Planck Collaboration Int. XLIX (2016) that there is a systematic effect in the U_{r} component due to uncertainty in the orientation of the polarization sensitive detectors.
In summary, given our understanding of the Planck data and simulations, the results from nonoriented stacking are consistent with the predictions of the ΛCDM model.
5.5.2. Oriented stacking
The stacking method of the previous subsection can be generalized by orienting the local coordinate frame of the patch to be stacked in a way that is correlated with the map being stacked. This approach, which we call “oriented stacking”, first introduced in PCIS15, allows extra information to be extracted from the stacked data. For unoriented stacking, the ensemble average cannot result in any intrinsic angular dependence, since it would be averaged by the uncorrelated orientation choices. For example, the ensemble average of the combination of polarization Stokes parameters Q + iU around unoriented temperature peaks has overall angular dependence e^{2iϕ}, which can be removed by a local rotation (Eq. (2)), as was carried out in the previous section. However, for oriented stacking, the angular dependence is a linear combination of a few Fourier modes, e^{imϕ}, with the exact mode content determined by the spin of the field being stacked and the spin of the orientation operator. For scalar fields stacked using orientations determined by the spin2 operator, only m = 0 and m = 2 modes are present.
Choices of what to stack, where to centre the patch, and how to orient it, provide a multitude of statistical tests that are complementary to the auto and crosscorrelation power spectra; these can be used to characterize nonGaussian data (such as polarized foreground emission) and have the advantage of being easy to visualize. Here we focus on patch positions and orientations determined by the highest signaltonoise channel, namely temperature T, and present stacks of temperature and polarization data on temperature peaks (either all peaks, or just those selected above a certain threshold ν). The orientation of the patch can be random, or determined by the second derivative of the temperature (since gradients vanish at the peak). While the most straightforward way to orient a patch centred on a temperature peak is to align the axes with principal directions defined by a local quadratic expansion of the temperature field around the peak, this would be susceptible to noise, and a better choice is to use the principal directions of a local quadratic expansion of the inverse Laplacian of the temperature ∇^{−2}T, namely the tensor ∇_{i}∇_{j}∇^{−2}T.
The inverse Laplacian of the temperature ∇^{−2}T is readily computed using the spherical harmonic transform
and its covariant derivatives can be computed using standard HEALPix routines. Perhaps a simpler interpretation of the tensor ∇_{i}∇_{j}∇^{−2}T is suggested by its projection onto the Stokes parameters using spin2 spherical harmonics:
This highlights the fact that the direction defined by the tensor relates to the temperature T in the same way as the direction of polarization relates to the E mode. In the flatsky approximation, the temperaturederived Stokes parameters Q_{T} and U_{T} are
We orient the patch so that U_{T} vanishes and Q_{T} is positive for the central peak. We use an equalarea azimuthal Lambert projection to represent the stacks, which introduces the radial variable
which is almost identical to the usual angular variable θ in the flatsky limit, but allows for less deformation if large angular sizes are considered. Rectangular coordinates on the patch are defined as usual by x = ϖ cos ϕ and y = ϖ sin ϕ. Given the stacked field X_{stack}(ϖ, ϕ), information can be further compressed by extracting radial profiles of the nontrivial angular modes:
where δ_{m0} is the Kronecker delta and X stands for T or E. Statistical isotropy of the stacked field X would imply that ensembleaveraged radial functions for odd ms vanish. Similar to nonoriented stacking (Eq. (34)), the integrated profile deviation,
can be used to evaluate the consistency between the data and the model. Once again, the weighting function W is chosen to be proportional to the expected profile .
To assess foreground contamination in componentseparated CMB temperature and polarization data, we present unoriented and oriented stacks of temperature T and fullsky reconstructions of E and B modes for Planck CMB maps compared to low and highfrequency foregrounds in Figs. 20 and 21, respectively. The stacking procedure is identical, except for the data sets used. The top rows of these figures show stacks of SMICA componentseparated CMB maps at N_{side} = 1024, with a 10′ FWHM beam. The middle rows show stacks of the LFI 30GHz maps at N_{side} = 1024 (with leakage correction applied), further convolved with a 10′ FWHM Gaussian beam, and corrected for the CMB contribution by subtracting the SMICA CMB map (at 10′ FWHM resolution) itself beamconvolved with the LFI 30effective beam, thus representing the total lowfrequency foreground emission. The bottom rows show the HFI 353GHz polarizationsensitive bolometer maps convolved with a 10′ FWHM Gaussian beam and degraded to N_{side} = 1024, then CMBcorrected using the SMICA map further beamconvolved with the HFI effective beam, thereby representing the total highfrequency emission. The polarization maps are transformed to E and B modes via the fullsky spherical harmonic transform, while the Q_{T} and U_{T} signals used to define patch orientations are derived from fullsky temperature maps. Temperature peaks are selected within the common intensity mask above the ν = 0 threshold, and patches either rotated randomly, as shown in Fig. 20, or oriented so that U_{T} = 0 is at the peak, as in Fig. 21. The maps to be stacked are masked with the common intensity mask for the temperature channel, and with the confidence mask of Appendix A for the polarization channels; monopole and dipole patterns in unmasked regions are fitted and subtracted out from all maps (a procedure that is more applicable to CMB maps, but is also applied to foregrounds as well for consistency), and the result stacked on temperature peaks selected as described above. All stacks are presented in units of μK_{CMB}, and cover square patches of size 4° by 4°.
Fig. 20. Nonoriented stacks of intensity T (left), reconstructed Emode (centre) and Bmode (right) polarization stacked on intensity maxima for the CMB (top row, SMICA map at N_{side} = 1024 and 10′ FWHM beam), lowfrequency foregrounds (middle row, LFI 30GHz map smoothed using a 10′ FWHM beam, with subtraction of the SMICA CMB map at 10′ FWHM resolution smoothed by the 30GHz LFI beam), and highfrequency foregrounds (bottom row, 353GHz HFI polarizationsensitive bolometer map smoothed using a 10′ FWHM beam, with subtraction of the SMICA CMB map at 10′ FWHM resolution smoothed by the 353GHz HFI beam). The common temperature mask was used for temperature peak selection, and the polarization confidence mask of Appendix A was used for E and Bmode stacks. 
Fig. 21. Oriented stacks of intensity T (left), reconstructed Emode (centre) and Bmode (right) polarization stacked on intensity maxima for the CMB (top row, SMICA map at N_{side} = 1024 and 10′ FWHM beam), lowfrequency foregrounds (middle row, LFI 30GHz map smoothed using a 10′ FWHM beam, with subtraction of the SMICA CMB map at 10′ FWHM resolution smoothed by the 30GHz LFI beam), and highfrequency foregrounds (bottom row, 353GHz HFI polarizationsensitive bolometer map smoothed using a 10′ FWHM beam, with subtraction of the SMICA CMB map at 10′ FWHM resolution smoothed by the 353GHz HFI beam). The orientation of the stacked patch is chosen so that the principal directions of the tensor ∇_{i}∇_{j}∇^{−2}T are aligned with the horizontal and vertical axes of the stack. The common temperature mask was used for temperature peak selection, and the polarization confidence mask of Appendix A was used for the E and Bmode stacks. 
The CMB data show a characteristic ringing pattern in the temperature stacks, associated with the first acoustic peak, a high signaltonoise pattern in the Emode stack (as expected from the CMB TE correlation), and no evidence of detectable patterns in the Bmode stack (reaching a noise floor of about 0.05 μK), which serves as a null test for foreground contamination (given the fact that both low and highfrequency foregrounds display TB correlations, as described below). Polarized dust foregrounds in particular were previously investigated in Planck Collaboration XI (2020). Both low and highfrequency foreground stacks show no acoustic peak ringing, with oriented temperature stacks being notably different from the CMB. The highfrequency stacks show strong TE and pronounced TB correlations, likely driven by the properties of the polarized dust emission that is dominant at 353 GHz. Lowfrequency foreground stacks show weaker, but still noticeable, TE and TB correlations, with no prominent pattern in the E mode, as is the case for highfrequency foregrounds. All foreground stacks are quite distinct from CMB ones, and foreground leakage should result in extraneous correlations in CMB componentseparated stacks, which are not observed. Our results indicate that the residual foreground contamination in Planck componentseparated maps is below levels that are measurable by this method.
To investigate the quantitative agreement of the Planck componentseparated maps with fiducial cosmological and noise models, we compare the radial profiles (Eq. (39)) of temperature and polarization stacks for all componentseparation methods with identically processed FFP10 simulations. We inpaint temperature T within the common intensity mask, and reconstruct E and B modes using purified inpainting within the common polarization mask, as described in Appendix A. Monopoles and dipoles are fitted and subtracted from the resulting maps, and once again we stack on temperature peaks above a threshold ν = 0 within the common intensity mask oriented so that U_{T} = 0 at the peak. The confidence mask of Appendix A is used for the E and Bmode stacks, and the common intensity mask for T stacks. Radial profiles (Eq. (39)) for m = 0, 2, and 4 are extracted from the stacked data, and are compared to simulations. The results for m = 0 and 2 are shown in Fig. 22; the results for m = 4 are consistent with zero signal, as should be the case for scalar maps stacked with the spin2 orientation operator.
Fig. 22. Radial profiles of oriented stacks of intensity T (left two columns) and reconstructed E mode (right two columns) for Commander, NILC, SEVEM, and SMICA pipelines. The top panel in each plot shows the radial profile compared to the simulation average (grey line), while the lower panel shows the difference of the radial profile and simulation average on smaller scales. The shaded grey regions represent 68%, 95%, and minimumtomaximum bounds for individual realizations. 
The top panels of Fig. 22 show the observed radial profiles of oriented stacks compared to the ensemble average of the FFP10 simulations, while the lower panels show the differences on a magnified scale. The shaded grey regions represent 68% and 95% confidence, and minimumtomaximum bounds using profiles obtained from individual realizations. All observed radial profiles in oriented stacking agree with simulations within the variance bounds. Table 9 summarizes the uppertailed pvalues of the χ^{2} deviation (Eq. (32)) of oriented radial profiles, while Table 10 summarizes the twotailed pvalues of the integrated profile deviations (Eq. (40)). None of the pvalues appear anomalous.
pvalues of χ^{2} computed from the oriented stacking of hot and cold spots selected above the ν = 0 threshold.
pvalues of ΔX_{m} computed from the oriented stacking of hot and cold spots selected above the ν = 0 threshold.
To summarize this section, oriented stacking results reinforce the conclusion of the previous section that the stacked Planck CMB data are consistent with the predictions of the standard ΛCDM model, given our understanding of the instrument and simulations.
6. Anomalies in the microwave sky
The previous section established the lack of evidence for significant nonGaussianity in the Planck temperature and polarization data. Here we reconsider several noteworthy features detected both in the WMAP temperature sky maps, and later confirmed in the Planck analyses described in PCIS13 and PCIS15. These include a lack of largeangle correlations, a hemispherical power asymmetry (either a simple excess of power in one hemisphere or a continuous dipolar modulation of the CMB anistropy over the sky), a preference for oddparity modes in the angular power spectrum, and an unexpectedly large temperature decrement in the southern hemisphere. Tests that involve dipolar power asymmetry, either directly or via measures of directionality, are collected together in Sect. 7, but in this section we consider tests of other kinds of anomaly.
The existence of these features is uncontested, but, given the modest significances at which they deviate from the standard ΛCDM cosmological model, and the a posteriori nature of their detection, the extent to which they provide evidence for a violation of isotropy in the CMB remains unclear. It is plausible that they are indeed simply statistical fluctuations. Nevertheless, if any one of them has a physical origin, it would be extremely important, and hence further investigation is certainly worthwhile. However, given that the Planck temperature data are cosmic variance limited on both large and intermediate angular scales, new information is required to determine if a real physical effect on the primordial fluctuations is indicated, or otherwise. This can be achieved by the analysis of the Planck polarization data.
Of course, the Emode polarization is partially correlated with the temperature anisotropy, so that it is not a fully statistically independent probe of the anomalies. Indeed, this correlation could result in a polarized feature due to the presence of a chance fluctuation in the temperature map, or could modify any intrinsic polarization anomaly. In principle, it is possible to split the polarization (temperature) signal into two parts, one that is correlated with the temperature (polarization), and one that is uncorrelated. Such an approach has already been applied to the WMAP data by Frommert & Enßlin (2010). However, the methodology requires a mathematical description of the noise plus residuals of the systematic and componentseparation effects that does not exist for the componentseparated data that we analyse here. Therefore our initial approach is simply to test for evidence of anomalies in the polarization data, in addition to verifying once again those seen previously in the temperature maps.
6.1. Lack of largeangle correlations
We assess the lack of correlation in the 2point angular correlation function at large angular separations, as previously noted for both the WMAP and Planck temperature maps (Bennett et al. 2003; Copi et al. 2015; PCIS15). In particular, we extend the analysis to polarization data, which were previously too noisy and/or contaminated by residual systematic artefacts on such scales to use them for verification of the temperature anomaly. We consider the statistic proposed by Copi et al. (2013) in an analysis of the WMAP data:
where X, Y can denote the temperature anisotropy T or the two Stokes parameters Q_{r} and U_{r} (as defined in Sect. 5.2), and is an estimate of the corresponding 2point correlation function. This is a generalization of the S_{1/2} statistic (Spergel et al. 2003) computed from the temperature autocorrelation function where the lower, θ_{1}, and upper, θ_{2}, limits of the separation angle range considered are 60° and 180°, respectively.
To check the consistency of the current data set with previous analyses, we again adopt this range of separation angles and present the results in Tables 11 and 12. We find that the temperature data show a lack of correlation 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). The pvalues for the temperature maps are slightly larger than those determined from the 2015 data set (PCIS15). This could be caused either by changes in the mask, or the inclusion of systematic effects in the FFP10 simulations. However, in temperature, the latter are relatively small compared to the cosmological signal on large angular scales.
Values for the statistic (in units of μK^{4}) for the Planck 2018 data with resolution parameter N_{side} = 64.
Probabilities of obtaining values for the statistic for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2018 CMB maps with resolution parameter N_{side} = 64, estimated using the Commander, NILC, SEVEM, and SMICA maps.
Copi et al. (2013) demonstrated that the signaltonoise ratio in the WMAP polarization maps was insufficient to allow meaningful estimates of S^{XY} to be made. For the Planck polarization data, we have estimated errors on the statistics from the standard deviation of these values computed from the corresponding FFP10 HMHD maps. We find that the errors are σ_{S1/2} ≲ 5 × 10^{−5} μK^{4} for the statistic, σ_{S1/2} ≲ 7 × 10^{−6} μK^{4} for the remaining statistics with estimated values at least of order 10^{−5} μK^{4}, and σ_{S1/2} ≲ 5 × 10^{−7} μK^{4} for the statistics with estimated values of order 10^{−6} μK^{4}. Similar errors are obtained from the OEHD maps.
A possible explanation for the lack of correlation in the temperature maps is due to the low observed value of the quadrupole. We therefore repeat the analysis after removing the bestfit quadrupole^{10} from the temperature maps. The corresponding probabilities, estimated from similarlycorrected FFP10 simulations, are recorded in the second row of Table 12 and indicate that the low power in the quadrupole alone does contribute to the absence of largeangle correlations. Copi et al. (2009) argue that all modes below ℓ ≤ 5 contribute to this, by cancellation with each other and with higher order modes.
A potential criticism of the statistic relates to the a posteriori choice of the range of separation angles to delineate the interesting region of behaviour of the correlation function. As in PCIS15, we consider the generalized statistic S^{XY}(θ, 180° ) and compute it for all values of θ both for the data and simulations. Figure 23 indicates that the only excursion outside of the 95% confidence regions determined from simulations is observed for the temperature data when the lower bound of the integral is θ ≈ 60°, for all componentseparation methods. No equivalent excursions are seen for any statistics involving polarization data, which also show notable variations between componentseparation methods. We provide a more quantitative assessment as follows. For each value of θ, we determine the number of simulations with a higher value of S^{XY}(θ, 180° ), 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 lookelsewhere effect, 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 θ at which it occurs. The fraction of these probabilities lower than the maximum probability found for the data defines a global pvalue. As seen in Table 13, this corresponds to values of order 99% for all of the componentseparated temperature maps, and much smaller global pvalues for quantities involving the polarization data. Note that the results for the temperature data are slightly more significant than those for the 2015 data set (i.e., 99% compared to 98%). More importantly, the pvalues for the temperature maps after removing the bestfit quadrupole drop to around 86%, indicating again that the anomalous value of the S^{TT} statistic is connected to the observed quadrupole.
Fig. 23. S^{XY}(θ, 180° ) statistic for the N_{side} = 64 Planck CMB 2018 temperature and polarization maps. Results are shown for the Commander, NILC, SEVEM, and SMICA maps (first, second, third, and fourth rows, respectively). The solid lines correspond to the data, while the black three dotsdashed lines indicate the median determined from the corresponding FFP10 simulations, and the shaded dark and light grey areas indicate the corresponding 68% and 95% confidence areas, respectively. 
Global pvalue for the S^{XY}(θ, 180° ) statistic for the Planck fiducial ΛCDM model at most as large as the observed values of the statistic for the Planck 2018 CMB maps with resolution parameter N_{side} = 64.
As pointed out by Copi et al. (2013), the crosscorrelation between temperature and polarization can be used to determine whether the S^{TT} value is low due to a statistical fluke in the context of a ΛCDM cosmology. Specifically, if this hypothesis were true, the S^{TQ} statistic would also likely be small on large angular scales, given that the polarization signal is partially correlated with temperature. Copi et al. (2013) have determined that the S^{TQ} statistic would be smaller than 1.403 μK^{4} over the angular separation range [48° ,120° ] in 99% of their constrained realizations based on the properties of the WMAP sevenyear data and assuming the statistical fluke hypothesis to be true. A value of the measured statistic exceeding this limit would then allow the hypothesis to be ruled out. Table 14 indicates that the values of the S^{TQ}(48° ,120° ) statistic measured with the Planck data are significantly smaller, so that we cannot rule out the statistical fluke hypothesis.
Values for the S^{TQ}(48° ,120° ) statistic (in μK^{4}) and the probability of obtaining values for the Planck fiducial ΛCDM model at least as large as the observed values for the Planck 2018 CMB maps with resolution parameter N_{side} = 64.
Muir et al. (2018) have noted that the low value of the largeangle correlation statistic can be connected to the temperature 2point correlation function for antipodean points, . This is motivated by fact that, as seen in Fig. 8, the otherwise largely flat drops to negative values close to θ = 180°. Here, we extend the analysis of the temperature correlation function to the functions including polarization information, .
Table 15 presents the probabilities of finding larger values of the 2point correlation functions at 180° from the data than those estimated using the FFP10 simulations. The results determined from the temperature maps are consistent with those reported by Muir et al. (2018), i.e., around 89% (although note that the sense of the pvalues differs between the papers). When the bestfit quadrupole is removed, the significance falls to around 60%, indicating that the low value of the quadrupole can impact the statistic. We do not observe any anomalous behaviour for the functions evaluated at this angular separation that include polarization fields.
Probabilities of getting the statistic for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2018 CMB maps with resolution parameter N_{side} = 64.
6.2. Hemispherical asymmetry
In this section, we reassess the asymmetry between the realspace Npoint correlation functions computed on hemispheres, reported previously for the WMAP (Eriksen et al. 2005) and Planck temperature maps (PCIS13; PCIS15). We pay special attention to testing this asymmetry using the Planck polarization data. Several different 2point and 3point functions drawn from permutations of the T, Q, and U maps are tested.
As in Sect. 5.2, we analyse the CMB estimates at a resolution of N_{side} = 64 and use the same configurations of the Npoint functions. 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).
We test the asymmetry for two separate cases: (i) the hemispheres determined in the ecliptic coordinate frame (i.e., those hemispheres separated by the ecliptic equator); and (ii) those determined by the dipolemodulation (DM) analysis in Sect. 7.2. In the latter case, the positive and negative hemispheres are defined as those for which the DM amplitude is positive or negative, respectively. The DM direction adopted here, (l, b) = (221° , − 20° ), then corresponds to the pole of the positive hemisphere, taken to be the average over the results determined from the four componentseparated maps in Table 24, for the data combination TT, TE, EE. We do not include an additional analysis of the hemispheric split defined by the Dopplerboost direction as presented in PCIS15, since the observed asymmetry described there was less significant than observed for the ecliptic and DM frames. The results for cases (i) and (ii) are presented in Fig. 24, where, as in PCIS15, we compute differences between the Npoint functions for the data and the mean values estimated from the FFP10 MC simulations.
Fig. 24. Difference of the Npoint correlation functions determined from the N_{side} = 64 Planck CMB Commander 2018 temperature and polarization maps and the corresponding means estimated from FFP10 MC simulations. Results are shown for the 2point (upper panels), pseudocollapsed 3point (middle panels), and equilateral 3point (lower panels) functions. The blue and red dashed lines correspond to the functions computed on the negative and positive hemispheres, respectively, determined in the dipolemodulation coordinate frame. The blue and red solid lines correspond to the functions computed on the northern and southern ecliptic hemispheres, respectively. The shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively, determined for the negative DM hemisphere. See Sect. 5.2 for the definition of the separation angle θ. 
As in Sect. 5.2, we quantify the agreement of the correlation functions of the CMB estimates with the fiducial cosmological model using a χ^{2} statistic defined by Eq. (11), where the correlation function, mean, and covariance matrix are computed for a corresponding hemisphere. If we consider that the χ^{2} statistic itself can act as a measure of fluctuation level, then 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 for the Planck fiducial ΛCDM model at least as large as the observed values are given in Tables 16 and 17 for the ecliptic and DM reference frames, respectively. The corresponding probabilities for the ratio of the χ^{2} values are given in Tables 18 and 19. 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.
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 for the Commander, NILC, SEVEM, and SMICA temperature and polarization Q and U maps estimated on northern and southern ecliptic hemispheres.
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 statistics for the Commander, NILC, SEVEM, and SMICA temperature and polarization Q and U maps estimated on negative and positive hemispheres defined by the DM reference frame.
Looking at the results in Tables 16–19, we no longer observe the high significance level for the pseudocollapsed 3point function of the temperature map in the northern ecliptic hemisphere, as reported for the 2013 and 2015 Planck data sets (PCIS13; PCIS15). This may be a consequence of the use of different masks in the various analyses, or of the improved treatment of poorly determined modes in the estimated correlation matrix used for the computation of the χ^{2} statistic, as described in Sect. 5.2. The largest asymmetry, at a significance of around 97–99%, is observed for the TTQ_{r} equilateral 3point function (although not for the SEVEM map). This comes from the very high probabilities found for the northern ecliptic hemisphere. It is also worth noting that the pvalues for the SMICATU_{r} 2point function and the NILCQ_{r}U_{r}U_{r} pseudocollapsed 3point function are very similar for the two hemispheres; as a consequence, the probabilities for the corresponding χ^{2} ratios are very small, typically of order 1%. We cannot offer any explanation as to the similarity of these correlation functions in the two hemispheres. However, one should recognize that a posteriori choices for the smoothing scale and reference frames defining the hemispheres will lead to an overestimation of the significance of the results. In addition, the large number of correlationfunction configurations considered will increase the likelihood of finding some apparently anomalous behaviour, even for statistically isotropic data.
Probabilities of obtaining values for the ratio of χ^{2} of the Npoint functions for the Planck fiducial ΛCDM model at least as large as the observed values of the statistics for the Commander, NILC, SEVEM, and SMICA temperature and polarization Q and U maps estimated on northern and southern ecliptic hemispheres.
Probabilities of obtaining values for the ratio of χ^{2} of the Npoint functions for the Planck fiducial ΛCDM model at least as large as the observed values of the statistics for the Commander, NILC, SEVEM, and SMICA temperature and polarization Q and U maps estimated on negative and positive hemispheres defined by the DM reference frame.
In the DM reference frame, the 3point temperature correlation functions in the negative hemisphere are somewhat significant, reaching a level of 98–99% for the pseudocollapsed case. Similar behaviour is also observed for some of the componentseparated maps for the Q_{r}Q_{r}Q_{r}, TTU_{r} and TQ_{r}Q_{r} pseudocollapsed 3point functions, and for a few configurations of the equilateral 3point function. Of course, the inconsistency in pvalues for different componentseparated maps argues against results of cosmological significance, but rather indicates the presence of different levels of foreground or systematiceffect residuals, depending on the componentseparation method.
6.3. Pointparity asymmetry
Now we turn to another way in which statistical isotropy can be broken at large scales. The CMB temperature anisotropy field on the sky can be written as the sum of paritysymmetric and parityantisymmetric functions defined as
where for a given direction , the antipodal point is . Under the parity transformation , it can be shown that the corresponding a_{ℓm} coefficients transform such that . As a consequence, and are comprised of spherical harmonics with only even or odd ℓmodes, respectively.
For the polarization field, a similar approach can be followed by expanding the Q and U stokes parameter maps as
It can then be shown that these functions are described by E and Bmode angular power spectra with only even or odd ℓmodes, since under the parity transformation and .
On the largest angular scales, 2 < ℓ < 30, corresponding to the SachsWolfe plateau of the temperature power spectrum, we expect the angular power spectra of and to be of comparable amplitude. However, an odd pointparity preference has been observed both in various WMAP data releases (Land & Magueijo 2005a; Kim & Naselsky 2010; Gruppuso et al. 2011), and confirmed in the Planck 2013 and 2015 studies. Furthermore, Land & Magueijo (2005b) have shown that this cannot be easily ascribed to the presence of residual Galactic foreground emission because of its even pointparity nature. Here we investigate the parity asymmetry of the Planck 2018 temperature data at N_{side} = 32 using the same estimator adopted in PCIS15, which is defined as follows:
where and are given by
with being 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 maximumlikelihood (QML) estimator (Gruppuso et al. 2011; Molinari et al. 2014).
The topleft panel of Fig. 25 presents the ratio, R^{TT}(ℓ_{max}), determined from the 2018 componentseparated maps after applying the common mask, together with the distribution of the SMICA MC simulations that are representative of the expected behaviour of the statistic in a Universe with no preferred largescale parity. We additionally consider an additional estimate of the CMB sky determined using the Commander componentseparation methodology that has been optimized for large scale analyses. This map is also used in the Planck lowℓ likelihood analysis (Planck Collaboration V 2020), thus we refer to it as LklCommander. This corresponding result is shown in the top right panel of the figure. The results from the componentseparation products are in good agreement, and indicate an oddparity preference for the multiple range considered in this test. The lowerleft panel of Fig. 25 shows the lowertail probability as a function of ℓ_{max} for the R^{TT}(ℓ_{max}) ratios shown in the upper panels. The profiles of the cleaned CMB maps are very consistent, and indicate a lowertail probability of about 1% over the range of multipoles ℓ_{max} = 20 − 30. This is higher than the typical value found in PCIS15, a difference that we ascribe to changes in the common mask. Indeed, if the 2015 maps are analysed after applying the 2018 common mask, then a lowertail probability of about 1% is also found. Similarly, if the 2015 common mask is applied to the 2018 data, then a probability of 0.3% is found at ℓ_{max} = 27, in good agreement with the earlier results.
Fig. 25. Upper panels: ratio R^{TT}(ℓ_{max}) of the Planck 2018 data 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, respectively. Left panel: ratio computed from the componentseparated maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (cyan), after application of the common mask. Note that these results substantially overlap each other. Right panel: ratio computed when the LklCommander map (magenta) is analysed with the corresponding mask. Lower panels: lowertail probability estimator as a function of ℓ_{max} (left), and ℓ_{min} with a fixed ℓ_{max} = 24 (right). For these lower panels small probabilities would correspond to anomalously odd parity. 
The LklCommander map shows good agreement with the standard componentseparated maps when the 2018 common mask is applied. However, if the confidence mask specifically associated with the lowℓ likelihood is used, then, as seen in the top right and lower panels of Fig. 25, lower probabilities are seen over the ℓ_{max} range of 20–30, with a minimum lowertail probability of 0.2% determined for ℓ_{max} = 24. This behaviour is in agreement with previous studies (see, for example, Gruppuso et al. 2018), which indicate that considering a reduced sky coverage, where the analysed regions are at higher Galactic latitudes, diminishes the odd point parity preference. In addition, we also observe probabilities of about 0.9% in the multipole region 7–9. At the spectral level, the different parity preference with respect to the 2018 common mask and the lowℓ mask seems to be associated with shifts in the amplitudes of multipoles at ℓ = 3 and 7.
In order to quantify the impact of any a posteriori effects on the significance levels of the LklCommander analysis, we count the number of MC simulations with a lowertail probability equal to, or lower than, 0.2%, for at least one ℓ_{max} value over a specific multipole range. For ℓ_{max} in the range 3–64, 16 out of the 999 MC maps have this property, implying that, even considering the lookelsewhere effect, an oddparity preference is observed with a lowertail probability of about 1.6%.
In the lowerright panel of Fig. 25, the lowertail probability as a function of ℓ_{min} is presented when ℓ_{max} is constrained to be 24, corresponding to the lowest probability found for the LklCommander analysis. It is apparent that the probability increases as the lowest multipoles are omitted, demonstrating that the anomaly is mostly driven by the largest scales.
We then extend the pointparity analysis to include polarization data, specifically considering the TE and EE power spectra. As in the temperature analysis, the data are compared with MC simulations that describe the expected parity preference for the ΛCDM model. However, the ratio estimator cannot be used in this case, since the signaltonoise is lower than for the temperature maps; this means that estimates of the power spectra may become negative, so that numerical problems can arise if the denominator of Eq. (44) approaches zero when summing values up to ℓ_{max}. To avoid this, a less sensitive but more robust estimator is considered:
where X corresponds to either the TE or EE spectra, and is defined analogously to the temperature case in Eq. (45).
In order to minimize the effect of incomplete knowledge of the instrumental noise, and given that some noise mismatch has been observed between the data and simulations (Planck Collaboration IV 2020), we determine the power spectrum using the crossquasiQML estimator that forms the basis of the HFI lowℓ likelihood in polarization (Planck Collaboration V 2020). Here, the signal covariance matrix is constructed according to the FFP10 fiducial model, while the noise is specified using the FFP8^{11} 143GHz noisecovariance matrix. This is clearly a suboptimal description of the noise present in the componentseparated maps, but only affects the variance of the estimated power spectra. However, since our statistical tests do not require this quantity, but are based on the dispersion of the MC simulations, the results are unaffected by the choice of the noisecovariance matrix. To verify this, we also analyse the SEVEM 143GHz cleaned map, the noise properties of which should be more closely matched by this covariance matrix. The consistency of results between the various componentseparated maps indicate that no bias arises from this choice. Finally, to test consistency, we apply the crossestimator to both the HM and OE data sets.
Examination of the EE and TE power spectra determined from the componentseparated maps indicates that the largest angular scales, in particular ℓ = 2, are likely to be affected by the presence of residual systematics, probably arising mainly from the ADC nonlinearity (Planck Collaboration III 2020) and their correlation with foregrounds. These effects are not fully captured by the current set of MC simulations. As a consequence, the inclusion of the ℓ = 2 multipole in the pointparity analysis results in probabilities that reach values of 98.5% for certain componentseparated products. In what follows, we limit the multipole range of interest to ℓ = 3 − 30.
Figure 26 presents the results for the D^{TE}(ℓ_{max}) estimator. The upper panels show, as examples, the distributions determined from the SMICA data and corresponding MC simulations. Results from the other component separated maps behave similarly. In the lower panels, we show the lowertail probability as a function of ℓ_{max} for the HM and OE data. Results from the componentseparated maps show similar trends, but with a range of probabilities. Nevertheless, they are compatible with the MC simulations and show a mild tendency towards decreasing probability with higher ℓ_{max}, although the lowest value observed is of order 20%. Results from the HM and OE data splits are in good agreement.
Fig. 26. Upper panels: D^{TE}(ℓ_{max}) for the SMICA HM (left panel) and OE (right panel) data. The shaded grey regions indicate the distribution of the statistic derived from the corresponding MC simulations as in Fig. 25. Lower panels: lowertail probability of the polarization estimators as a function of ℓ_{max} for the HM data (left panel) and OE data (right panel) for all componentseparated methods. Results for the SEVEM 143GHz cleaned maps are also shown (magenta line). For these lower panels small probabilities would correspond to anomalously odd parity. 
Figure 27 shows equivalent results for the D^{EE}(ℓ_{max}) estimator, though with Commander taken as a representative example of the behaviour of the componentseparated maps. As in the TE analysis, the data and MC simulations are consistent, although the pvalues decrease over the ℓ_{max} range of 20–30, reaching minima of about 6% at ℓ_{max} = 27 for the HM data; however, the OE results do not indicate such a trend.
Fig. 27. Upper panels: D^{EE}(ℓ_{max}) for the Commander HM (left panel) and OE (right panel) data. The shaded grey regions indicate the distribution of the statistic derived from the corresponding MC simulations as in Fig. 25. Lower panels: lowertail probability of the polarization estimators as a function of ℓ_{max}, as described in Fig. 26. For these lower panels small probabilities would correspond to anomalously odd parity. 
It is interesting to note that the statistic introduced in Sect. 6.1 is related to the measurement of parity asymmetry (see also Muir et al. 2018). In particular, when the functions are expressed in terms of the angular power spectra, it can be seen that the following relationships hold:
The expected values for the remaining correlation functions are zero, i.e., . The factor ( − 1)^{ℓ} splits the sum into contributions from even (paritysymmetric) and odd (parityantisymmetric) multipoles with opposite signs. The functions are therefore differences between parityeven and parityodd multipole bandpowers, and can be compared to Eq. (46). However, the expressions are not identical, since, in the latter case, the power spectrum is weighted by the factor ℓ(ℓ + 1), as compared to the (2ℓ+1) weighting used in the former case. Nevertheless, the results presented in Table 15 seem to be consistent with those in this section.
In summary, no anomalous lowertail probability is found. However, the low signaltonoise of the Planck polarization data over this range of scales is a limiting factor for this analysis. Future higher sensitivity observations are certainly required in order to investigate the point parity asymmetry in polarization.
6.4. Peak distribution asymmetry
Localized anomalies on the CMB sky can be searched for by testing how the statistical properties of local extrema (or peaks) vary in patches as a function of location. Since we expect the asymmetry measured by the QML estimator in Sect. 7.2 to be mirrored by a slight difference in the temperature and polarization peak distributions in the corresponding positive and negative hemispheres, we present and examine these quantities here. Note that the test is not powerful enough in itself to determine a modulation direction, but relies on the dipole orientation results presented in Table 24.
We extract peaks from the filtered maps discussed in Sect. 5.4 within the 70° radius discs centred on the positive and negative asymmetry directions determined by the SMICATT, TE, EE QML estimator (see Table 24). Separate peak CDFs are then constructed for these directions, then compared to the fullsky peak distribution, and the median CDF of the simulations. The results are presented in Fig. 28, which indicates the Kolmogorov–Smirnov deviation from the median simulation CDF of the full sky. No significant difference between the Emode peak distributions for the positive and negative directions is observed when the data are filtered on a 120′ scale. For a 600′ filtering scale, a marginal difference between the data and simulations is observed in the negative direction, as previously observed for the temperature results in PCIS15.
Fig. 28. Kolmogorov–Smirnov deviation of the peak distribution for the SMICA temperature T (left) and the reconstructed Emode polarization (right) maps in 70° radius discs centred on the positive and negative asymmetry directions as determined by the SMICATT, TE, EE QML estimator in Sect. 7.2. Top panels: maps filtered with GAUSS kernels of 120′ and bottom panels: filtering with 600′ FWHM. 
6.5. The Cold Spot and other largescale peaks in temperature and polarization
The Cold Spot was first detected in the WMAP temperature data as an anomalously cold region in terms of the spherical Mexicanhat wavelet (SMHW) coefficients (Vielva et al. 2004; Cruz et al. 2005) and later confirmed with Planck data (PCIS13; PCIS15, and references therein). The shape and the local properties of the Cold Spot have been studied using a variety of statistical approaches (e.g., Cayón et al. 2005; PCIS15). Recently, MarcosCaballero et al. (2017a) analysed the local properties of the Cold Spot (and other largescale peaks) via multipolar profiles that characterize the local shape in terms of the discrete Fourier transform of the azimuthal angle. In that paper, the anomalous nature of the Cold Spot is identified by being one of the most prominent peaks in curvature, as quantified by comparison to the predictions of ΛCDM^{12}. Its internal structure was then subsequently studied by Chiang (2018) in the context of an analysis of the variation of the acoustic peak positions in different parts of the sky. It was found that the Cold Spot region shows a large synchronous shift of the peaks towards smaller multipoles.
In this section, we analyse the polarization pattern of the Cold Spot, considered as a minimum in the temperature field when smoothed with a Gaussian of standard deviation R = 5°, to search for evidence as to whether its origin is primordial or otherwise. Vielva et al. (2011) and FernándezCobos et al. (2013) provide forecasts on whether given experimental configurations can distinguish between these hypotheses, and at what level of significance. In addition to the Cold Spot, four additional largescale peaks are considered, selected as the most anomalous structures on large angular scales, here taken to be extrema with a filter of R = 10° (see MarcosCaballero et al. 2017a,b, for more details). Since these peaks dominate the temperature field on large angular scales, the particular value of the smoothing scale used in their identification is not critical for locating the peaks. The scale R = 10° is chosen in order to highlight these structures while preserving their geometrical properties, such as curvature and eccentricity. These features correspond to two maxima and two minima, shown as Peaks 1–4 in Fig. 29, all located in the southern Galactic hemisphere. Given their locations, they are expected to make a significant contribution to the hemisphericpower asymmetry. In addition, the peaks should induce a particular pattern in the polarization field, due to the correlation between the temperature and Emode polarization.
Fig. 29. Locations of the largescale maxima or minima, numbered 1 to 5, considered in the analysis and indicated on an inpainted version of the Commander map at N_{side} = 2048. In particular, Peak 5 corresponds to the Cold Spot. 
Our analysis of the peaks is based on the calculation of the T and Q_{r} angular profiles. Since these are computed by averaging the signal over azimuthal angle, the analysis only characterizes the circularly symmetric part of each peak. However, the T and Q_{r} profiles, in the absence of any TB correlation and neglecting the eccentricity of the peaks, contain all of the temperature and polarization information that we are interested in.
In order to have a complete characterization of the peaks, we have to calculate derivatives of the fields up to second order. Since we are interested in the azimuthally symmetric signal, only the peak height and the Laplacian are relevant to the Q_{r} profiles. Once these two quantities are calculated from the temperature smoothed at a given scale R, their values are used to construct the expected theoretical profile in the polarization field (see MarcosCaballero et al. 2016, for technical details).
In order to analyse the largescale peaks, the profiles are calculated by averaging the temperature and polarization fields over azimuthal angle in different bins of the radial distance θ. We use intervals in θ from 0 to 30°, with a width of 1°. In the case of polarization, the Q and U components are rotated to obtain the Q_{r} and U_{r} Stokes parameters. The mean profile and the covariance matrix are discretized in the same way as the observed peak profiles. Assuming that the CMB field is Gaussian, the statistical distribution of the profiles, which is obtained by conditioning the values of the peak height and the Laplacian, is also Gaussian. In the absence of noise, the covariance matrix corresponds to the cosmic variance, modified to account for the fact that the derivatives at the centre of the peak are fixed to the measured values. In the case of the temperature profiles, this contribution introduces large correlations between different angular bins, which leads to a problem when a mask is applied to the data. Incomplete sky coverage modifies the correlation between the bins of the data, producing a disagreement between the data and the theoretical predictions. In order to properly characterize the masked data, we adopt the methodology developed in MarcosCaballero et al. (2016) and generate 2000 CMB simulations that are constrained to have a peak located at the same position as every largescale peak under consideration. Since the covariance matrix does not depend on the value of the peak height and curvature, we use the same simulations for each of the largescale peaks with R = 10°. However, the covariance matrix depends on the scale of the peak, and thus separate simulations are employed for the Cold Spot. Since the mean profiles are not affected by the mask, we use the analytical expressions in MarcosCaballero et al. (2016) to calculate the theoretical models instead of the average value of the simulations. This procedure allows us to reduce the sample variance due to the finite set of realizations.
In polarization, the contribution of the peak to the covariance is small compared with the cosmic variance (MarcosCaballero et al. 2016) and is ignored in subsequent calculations. The polarization analysis is then carried out at N_{side} = 512 with the common mask applied. The cosmicvariance part of the covariance matrix is calculated theoretically from the angular power spectrum (MarcosCaballero et al. 2016, 2017a), whereas the noise part is estimated from the FFP10 noise simulations. In order to take into account the possible contribution of the anisotropic noise, the covariances are calculated at the specific locations of each peak. The consistency of the noise between data and simulations is verified by computing the profiles in the OEHD and HMHD maps. The noise values in the profiles at each peak location are compatible with those obtained from the simulations. Since we are analysing only large scales, the noise can be ignored in the case of the temperature profiles.
The temperature and polarization signals induced by the peaks can be estimated by rescaling the theoretical profiles by a factor A. The value of A is obtained using the maximum likelihood estimator assuming that the profiles are Gaussian. In this case, it is possible to show that the distribution of A is also Gaussian with the following mean and standard deviation:
Here and X(θ_{i}) are the measured profile and the corresponding theoretical expectation for the angular bin centred at θ_{i}, respectively. In these equations, X represents the T or the Q_{r} profile, depending on the case we are analysing. Since we are interested in largescale peaks, angles up to 60° are considered in the analysis, with a bin width of 1°. The matrix C_{ij} is the covariance of the angular bins, which includes both cosmic variance and noise.
Forecasts regarding the ability of polarization measurements to distinguish between the standard model and other alternatives make use of the Fisher discriminant (Vielva et al. 2011; FernándezCobos et al. 2013). This method is based on the overlapping of the probability distributions of the amplitude A for different hypotheses. In the case of the Cold Spot, the standard model and alternatives are compared, assuming zero correlation between temperature and polarization. The predicted significance levels for detection of the polarization signal of the Cold Spot in these papers is higher than presented here, since a different characterization of the theoretical profile is used; while the model we assume is based on the temperature and the Laplacian observed in the data, the mean value of the model in FernándezCobos et al. (2013) was calculated using only the Laplacian (i.e., SMHW coefficient). Moreover, the Cold Spot amplitude is characterized with respect to the observed standard deviation, not using the theoretical expectation, as in our. This difference in the theoretical models results in different significance levels for the forecasts. However, the normalization of the theoretical profiles is not relevant when the significance of the hypothesis (A/σ_{A}) is calculated from the observed amplitude A.
The T and Q_{r} profiles determined from the four component separation methods are shown in Fig. 30. The corresponding profile amplitudes are given in Table 20. The temperature results are compatible with the analysis of the same peaks performed in MarcosCaballero et al. (2017a). Regarding the polarization analysis, the covariance of Q_{r} is not affected by the presence of a peak in the field, which results in larger uncertainties for polarization than for temperature. Indeed, the signaltonoise ratio for the Q_{r} amplitude factors are only about 1, and the polarization signal of the peaks cannot be detected; more sensitive polarization data would be required to further test the model. The values of the amplitude factors measured from the peak profiles are consistent with the standard ΛCDM model.
Fig. 30. Temperature (left column) and Q_{r} (right column) profiles of the largescale Peaks 1 (top panel) through 5 (bottom panel), as shown in Fig. 29. The black solid lines represent the profiles obtained from the SEVEM CMB map. The red curves correspond to the expected theoretical profiles, while the blue lines represent the theoretical profiles rescaled by the estimated amplitude A for each profile. The shaded regions represent the ±1, ±2 and ±3σ confidence levels. 
7. Dipole modulation and directionality
In this section, we examine dipolar violations of isotropy. In PCIS15, we performed a nonexhaustive series of tests to try to elucidate the nature of the observed asymmetry in temperature, and we follow this approach again here, both to reconfirm the previous temperature results, and to search for evidence of equivalent signatures in the polarization data.
Despite warnings about the use of the 2015 Planck polarization maps for the study of isotropy and statistics, several papers have attempted to fit dipolarmodulation models to the data. Aluri & Shafieloo (2017) applied a localvariance estimator to lowresolution Emode maps derived from the 2015 Commander solution, and found a power asymmetry at the level of around 3% over the range ℓ = 20 − 240, with a preferred direction broadly aligned with the CMB dipole^{13}. The ℓ range was selected on the basis that an apparent noise mismatch between the data and the FFP8 simulations could be minimized by the application of a simple scaling factor. However, we note that this mismatch is almost certainly due to the absence of modelled systematic effects in the FFP8 simulations, rather than an actual miscalibration of the instrumental noise. Conversely, Ghosh & Jain (2018) applied a pixelbased method to maps of the squared polarization amplitude, but found no evidence for the presence of a dipolarmodulation signal. They suggested that this may be due to residual systematics masking a real effect, since strong clustering of the dipole directions inferred from the FFP8 simulations was also observed.
These results indicate the necessity to use improved simulations that characterize the data more completely, and to adapt the estimators to eliminate bias resulting from the anisotropic noise and systematic effect residuals present in both data and simulations. As usual, we test for inconsistencies between the data and simulations using null tests based on the HMHD and OEHD data splits.
All the tests in this section involve the fitting of a dipole. In Sects. 7.1 and 7.3, we fit a dipole explicitly to a map of power on the sky. On the other hand, in Sect. 7.2 we measure the coupling of ℓ to ℓ ± 1 modes in the CMB multipole covariance matrix. There are differences in how we combine the fitted dipoles, which determine the particular form of dipolar asymmetry that the test will be sensitive to.
The tests can also be distinguished by whether they are based on amplitude or direction. The approaches of Sects. 7.1 and 7.2 are both sensitive to the dipole modulation amplitude. In particular, in Sect. 7.1 we examine the data for dipolar modulation of the pixeltopixel variance, while in Sect. 7.2 we search for modulation of the angular power spectra. These approaches differ mainly in terms of their ℓ weighting. Directionality in the data is the subject of Sect. 7.3, with the directions determined by bandpower dipole fits.
Given these differences in the approaches in this section, it is important to keep in mind that the results cannot usually be directly compared, even though all probe some aspect of dipolar asymmetry.
7.1. Variance asymmetry
We first study the Planck 2018 temperature and polarization maps using the localvariance method introduced by Akrami et al. (2014). Despite its relative simplicity, the localvariance estimator serves as a powerful method for detecting violations of statistical isotropy if they are of a type corresponding to a largescale power asymmetry, as observed in the Planck 2015 temperature data (PCIS15).
Here, we closely follow the previous temperature analysis, but extend its application to the scalar Emode polarization data. The method can briefly be described as follows. We define a set of discs of various sizes uniformly distributed on the sky. The centres of the discs are defined to be the pixel centroids of a HEALPix map at some specific low resolution, here taken to be N_{side} = 16 for both temperature and polarization analyses. Since it is important to cover the entire sky, we do not work with discs of radii smaller than 4° (given our choice of N_{side} for the centroids). These combinations of N_{side} and disc radii have been shown to be adequate for detecting largescale power asymmetry in CMB temperature data, and allow us to focus on the question of whether a corresponding largescale anomalous power asymmetry exists in the polarization data. For each sky map (data and simulations), we first remove the monopole and dipole components from the masked map and then compute the variance of the fluctuations for each disc of a given size using only the unmasked pixels. This results in localvariance maps at the HEALPix resolution of N_{side} = 16. We also estimate the expected average and variance of the variances on each disc from the simulations, and then subtract the average variance map from both the observed and simulated localvariance maps. Finally, we fit a dipole to each of the localvariance maps using the HEALPix remove_dipole routine, where each pixel is weighted by the inverse of the variance of the variances that we have computed from the simulations at that pixel. Note that, at all stages of the analysis, we work only with those discs for which more than 10% of the area is unmasked. In this way, we define an amplitude and direction for the variance asymmetry of each map. We then compare the localvariance amplitudes for the observed data to those of the simulations, containing statistically isotropic CMB realizations, and assess the level of anisotropy in the data. Since Akrami et al. (2014) have shown that the amplitudes of higher multipoles in such fits to temperature data are consistent with statistically isotropic simulations, we work only with the dipole amplitudes of the localvariance maps here. However, although we do not consider highermultipole asymmetries, such features may exist in localvariance maps estimated from the polarization data.
Table 21 presents the results for a localvariance analysis of the fullresolution, N_{side} = 2048, Planck 2018 temperature data over a range of disc radii, 4° ≤r_{disc} ≤ 20°. The pvalues are measured as the fractional number of simulations with localvariance dipole amplitudes larger than those inferred from the data. The results are consistent for all the componentseparated maps, and also with those of the Planck 2015 analysis, although a higher level of significance is seen here. In particular, none of the simulations yield localvariance dipole amplitudes larger than those determined from the data for the 4°, 6°, and 8° discs. This illustrates that the Planck temperature sky is asymmetric if one chooses to focus on variance over this range of angular scale. We also see the expected increase in pvalues as the disc radius is increased to values larger than 8°.
pvalues for variance asymmetry measured using different discs for the fullresolution (N_{side} = 2048) Planck 2018 Commander, NILC, SEVEM, and SMICA temperature solutions.
Since our focus is on large angular scales, we repeat the localvariance analysis for lowresolution, N_{side} = 64, temperature maps. Figure 31 (upper panel) shows the pvalues for the componentseparated maps as a function of disc radius, over the range 4° ≤r_{disc} ≤ 90°. The preferred lowℓ modulation direction determined from the temperature data in Sect. 7.2 is also indicated. Excellent agreement is found between the componentseparation methods. The overall level of significance is lower compared to the fullresolution results, but high significance (low pvalue) results are found for small disc radii. The lowest pvalues correspond to the smallest disc radius, i.e., 4° discs, and the pvalues increase with disc size. Figure 31 (lower panel) also shows the preferred directions for the four maps when 4° discs are used. The observed directions are in excellent agreement with each other, with the results of the fullresolution analysis, and with the findings of the Planck 2015 analysis. The values of the 4°disc localvariance dipole directions for the four componentseparation methods are provided in Table 22 for both the highresolution and lowresolution cases.
Fig. 31. Upper panel: pvalues for variance asymmetry measured as the fraction of simulations with localvariance dipole amplitudes larger than those inferred from the data, for lowresolution (N_{side} = 64) temperature maps. The pvalues are given as a function of disc radius and for the four componentseparated temperature maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). Lower panel: corresponding localvariance dipole directions for the four componentseparation temperature maps, and for 4° discs with the lowest pvalues. Note that the four directions match almost perfectly, so that the symbols essentially overlap. For reference, we also show the CMB dipole direction, the north ecliptic pole (NEP), the south ecliptic pole (SEP), and the preferred dipolar modulation axis (labelled as “lowℓ”) derived from the temperature data in Sect. 7.2. 
Localvariance dipole directions for the variance asymmetry of the four componentseparated temperature maps.
The analysis of the polarization data is significantly more subtle than for temperature because of the inherently low signaltonoise. We first validate the technique by applying it to simulations. The analysis shows that the direct application of the method to isotropic simulations of Emode polarization signal returns localvariance dipole directions that are not uniformly distributed. This arises from the strongly anisotropic, correlated, and nonGaussian structure of the Planck polarization noise, which needs to be corrected before any statistical method of anisotropy detection is applied to the data. This is not an issue for the temperature maps, since the signaltonoise ratio is large.
The following procedure is adopted to minimize the effects of the anisotropic noise. We first preprocess the Emode polarization maps (both the data and simulations) pixel by pixel and apply a bias correction and weighting. We consider the following types of transformation to the maps:
Here, X_{i} is the value of the map at pixel i, M_{i} is the mean at pixel i computed from the simulations, σ_{i} is the standard deviation at pixel i (again computed from the simulations), and p ≥ 1 is an integer.
Figure 32 shows the distribution of the localvariance dipole directions for an analysis of the HMHD Emode simulation maps when unmodified, and after transformation according to Eq. (53) with p = 1 or 2. Since we are interested in largescale anomalies, we restrict the study to a resolution of N_{side} = 64, and provide results for the Commander data and 4° discs only; the results for the other componentseparation methods are very similar. The untransformed case clearly shows that the noise does not yield uniformly distributed dipole directions, while by applying the transformation given by Eq. (53) a significant improvement in the uniformity is observed. Since the nonuniformity for simulated maps including signal arises from the noise, this transformation also improves the uniformity of the recovered dipole directions. Additionally, we see slightly more improvement for p = 2, and adopt this weighting for the analysis of the polarization data.
Fig. 32. Impact of bias correction and weighting on the uniformity of dipole directions for Emode polarization HMHD simulation maps. Upper panel: localvariance dipole directions for the Commander simulations of Emode polarization maps, obtained for 4° discs, when no bias correction and weighting have been applied to the maps. Middle panel: same as in the upper panel, but when the transformation X_{i} → (X_{i} − M_{i})/σ_{i} has been applied to the polarization maps before applying the localvariance estimator. Lower panel: same as in the upper and middle panels, but when the transformation has been applied to the maps (i.e., with a square in the denominator). The results for the other componentseparation methods are very similar to the ones presented here. For reference, we also show the CMB dipole direction, the north ecliptic pole (NEP), the south ecliptic pole (SEP), and the preferred dipolar modulation axis (labelled as “lowℓ”) derived from the temperature data in Sect. 7.2. 
Figure 33 (upper panel) presents the pvalues obtained by applying the localvariance estimator to the four componentseparated Emode polarization maps. Even though we typically see an increase in pvalues between smaller and larger disc radii, as observed for the temperature data, the detailed behaviour differs. Moreover, the curves seem to be divided into two groups, one including Commander and SEVEM, and the other consisting of SMICA and NILC. This might reflect differences in the componentseparation approaches, particularly given that the former methods operate in the pixel domain, and the latter in the harmonic domain. There is also a likely connection with a variation in residual systematic effects present in the componentseparated maps, depending on the componentseparation methodology applied. Furthermore, these differences, and particularly the relatively high values for the SMICA and NILC maps, argue against a detection of cosmological power asymmetry in the polarization data.
Fig. 33. 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 Emode polarization maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), at the HEALPix resolution of N_{side} = 64. Lower panel: corresponding localvariance dipole directions for the four componentseparation Emode polarization maps, and for 4° discs. For reference, we also show the CMB dipole direction, the north ecliptic pole (NEP), the south ecliptic pole (SEP), and the preferred dipolar modulation axis (labelled as “lowℓ”) derived from the temperature data in Sect. 7.2. 
Nevertheless, the preferred directions we have found for the Emode polarization data, as shown in the lower panel of Fig. 33, are intriguingly close to those determined for the temperature data. In order to quantify the significance of this alignment, we have computed the separation angles between the Emode and temperaturevariance dipole directions for the data as well as for all the simulations, and computed pvalues defined as the fraction of simulations with T − E separation angles smaller than those inferred from the data. Again, the results cannot be interpreted as evidence of power asymmetry in polarization, and the SEVEM result is a notable outlier. Table 23 summarizes these results, providing the coordinates of the preferred localvariance dipole directions and corresponding pvalues for the four componentseparation methods, together with the separation angles between temperature and Emode dipole directions and their associated pvalues.
Localvariance dipole directions and pvalues for the variance asymmetry of the four componentseparated Emode polarization maps at the HEALPix resolution of N_{side} = 64.
Despite the bias correction and weighting procedure employed for reducing the nonuniformity of the dipole directions in simulations, it is clear that the treatment of anisotropic noise plays an important role in our analysis of the polarization data, and its impact on the variance asymmetry results may need further consideration. While the close alignment of the temperature and polarization directions could simply be a coincidence, future data sets may offer additional insight.
7.2. Dipole modulation: QML analysis
In this section, we use the QML estimator originally introduced in Moss et al. (2011) to assess the level of dipole asymmetry in the CMB sky, and further extend the analysis to polarization data. We note, however, that Contreras et al. (2017) found that Planck polarization data are unlikely to be able to distinguish between dipolemodulated skies and skies consistent with ΛCDM (although this statement is somewhat model dependent). Furthermore, the analysis we carry out here is a purely phenomenological one performed in multipole space, with no attempt to connect to any realspace modulation; however, several physical kspace models are considered in a companion paper (Planck Collaboration X 2020).
We employ a version of the estimator proposed in PCIS15, which was further developed in Zibin & Contreras (2017) and Contreras et al. (2017). In particular, we use
with
Here the δC_{ℓℓ′} are determined by the model of modulation; for this case we have assumed scaleinvariant power modulation out to some ℓ_{max} and thus δC_{ℓℓ′} = 2(C_{ℓ} + C_{ℓ′}). The remaining terms are: , the spherical harmonic transform of the modulation amplitude and direction; WZ= TT, TE, or EE; W_{ℓm} and Z_{ℓm}, which are inversecovariance filtered data; ; and the last term on the righthandside of Eq. (56) denotes the meanfield correction (details, including the precise form of can be found in Appendix A.1 of Planck Collaboration XV 2016 and Appendix B of Contreras et al. 2017). The parentheses in the superscripts indicate symmetrization over the enclosed variables. The coupling coefficients are given by
The can be transformed into amplitudes of modulation along each Cartesian axis, since they are simply the spherical harmonic decomposition of a dipole. Likewise we can write the modulation in terms of its spherical coordinates (amplitude and direction) as
Note that for Emode polarization, these coupling coefficients neglect corrections on the very largest scales (ℓ ≲ 10) due to the nonlocal definition of Emodes. This gives rise to a slightly different coupling induced in E as compared to T by a dipole modulation of the primordial fluctuations (Contreras et al. 2017). This has little effect on the comparison of the data with simulations, however, since the simulations are treated in the same way. A potentially more significant effect is a mismatch in the optical depth between data and simulations; below we estimate this to have essentially negligible influence.
The estimators of Eqs. (54) to (55) can be combined with inversevariance weighting over all data combinations (TT, TE, EE) to obtain a combined minimumvariance estimator, given by
We calculate the variance from the scatter of simulations, although they agree closely with the Fisher errors given in PCIS15.
Here we test for an ℓspace asymmetry in temperature and polarization. The model considered is a scaleinvariant modulation from ℓ_{min} = 2 out to a variable maximum multipole, ℓ_{max}. It is important to stress that, when we combine temperature and polarization, the phenomenological, ℓspace approach we adopt here is very simplistic and does not capture the behaviour of the modulated fluctuations arising from the different k–ℓ kernels for T and E modes (Contreras et al. 2017). For example, modulated fluctuations that exhibit a 7% dipolar asymmetry in T to ℓ ≈ 65 are not expected to produce an Emode asymmetry of the same amplitude and over the same scales. Therefore, strictly speaking, the “model” being tested when we combine TT, TE, and EE has no physical basis; a more physical, kspace approach is considered in the companion paper Planck Collaboration X (2020).
In Fig. 34 we show the pvalues of the fullmission data compared to the FFP10 simulations, considering TT (top left), EE (top right), or TE (bottom left) data independently. While the TT results are largely consistent with our previous analysis (differences are within the expected scatter given the different analysis masks, PCIS15), our EE and TE results are new. The polarizationonly results show mildly significant asymmetry to ℓ_{max} ≈ 250, and are featureless elsewhere. This is minimally affected by the mismatch in τ between the data and simulations. The effect of excluding the ℓ ≲ 10 data (where τ is most relevant) from the analysis modifies the amplitude estimation by of order 10% for scales ℓ ≳ 50. Thus a 10% mismatch in τ would (at most) correspond to a 1% error in the amplitude, which decreases at smaller scales. This error is negligible compared to the noise contribution on these scales. The temperaturepolarization crosscorrelation results are rather featureless in the full range of ℓ_{max} considered. In the ℓ_{max} ≈ 250 region NILC shows systematically lower pvalues compared to the other methods, although the modulation amplitudes are still statistically consistent. This difference could be attributable to a percent level pvalue excess of asymmetry in the OEHD data observed in the same region; however, such an excess does not appear in the HMHD data. No other combinations of OEHD and HMHD data show any significant excess of asymmetry.
Fig. 34. Probability determined from the QML analysis for a Monte Carlo simulation to have a larger dipolemodulation amplitude than the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) data sets, with ℓ_{min} = 2. The topleft, topright, and bottomleft panels use TT, EE, and TE, while the final panel is for all data combined. Small pvalues here correspond to anomalously large dipolemodulation amplitudes. We emphasize that the statistic here is cumulative and apparent trends in the curves can be misleading. 
In the bottom right panel of Fig. 34 we combine temperature and polarization (including TE) data and show the pvalues of the data compared to simulations. Only the pvalue dip at ℓ_{max} ≈ 250 falls lower than the corresponding dip in temperature alone. Moreover, it is important to stress that the combined significance increasing at that scale is not in itself evidence that the asymmetry has a genuine, physical origin. As pointed out in Contreras et al. (2017), the combination of a statistically isotropic polarization signal with temperature data will, simply due to Gaussian statistics, spuriously increase the significance of a ≈3σ temperature signal with 30% probability for Planck. Furthermore, our phenomenological model does not properly combine temperature and polarization in 3D kspace, as previously mentioned (also see Zibin & Contreras 2017, and references therein).
In Table 24 we show the measured amplitude and direction of dipolar modulation in the oftquoted ℓ = 2 − 64 range, with and without polarization data. Polarization clearly has very little effect on the modulation parameters in this region, and TE in particular has a nearly negligible effect at this scale. It is worth recalling that the expectation for purely Gaussian skies is a dipolemodulation amplitude of approximately 3% (see Eq. (50) in PCIS15). In Table 25 we show the amplitude and direction of dipolar modulation for the combined TT, TE, and EE data in the ℓ = 2 − 220 range, which is the range with the lowest pvalue. The fact that the amplitude for this ℓrange is smaller than that for ℓ = 2 − 64 is consistent with the prediction for statistically isotropic skies (as noted in PCIS15, the amplitude typical of modulation should decrease as 1/ℓ_{max}). The proximity of the directions observed for the two scales, ℓ ≤ 64 and ℓ ≤ 220, was also noted in PCIS15, where tests of directionality were performed. There we showed that the two directions are correlated at a slightly higher level than seen in simulations, but that this can be traced to the lowℓ anomaly on the larger scales where the dipolemodulation amplitude is larger. Removing angular scales ℓ < 100 eliminates the significance of the modulation entirely.
Amplitude and direction of the lowℓ dipoleasymmetry signal determined from the QML analysis for the range ℓ = 2 − 64.
Amplitude and direction of the lowℓ dipole modulation signal determined from the QML analysis for the range ℓ = 2–220.
In PCIS15 we applied a lookelsewhere correction to the corresponding results and demonstrated that if we allow the range of ℓ_{max} to vary (as we have done here) then small pvalues (of order 0.1–1%) occur of order 10% of the time. Repeating that analysis here would yield a similar result. We have found that, as expected, the Planck polarization data have not been able to refute or confirm the original signal found in temperature. The polarization data alone also appear to be consistent with statistical isotropy. Better polarization data are required to further test whether there might be a physical origin for the original temperaturemodulation signature.
7.3. Angular clustering of the power distribution
Despite the lack of evidence for any strong anomaly in the amplitude of dipole modulation discussed in the previous subsections, in PCIS15, we confirmed the apparent presence of a deviation from statistical isotropy in the Planck data using an angularclustering analysis, as previously seen in PCIS13. In particular, some alignment of preferred directions determined from maps of the temperature power distribution on the sky was observed over a wide range of angular scales.
Specifically, by calculating the power spectrum locally in patches, for various multipole ranges, and then fitting dipoles to maps of these bandpower estimates, it was found that the directions were aligned at the 2–3σ level up to ℓ = 1500 when compared to simulations. In the standard cosmological model, although such maps are expected to exhibit dipolar distributions of power due to Gaussian random fluctuations, the associated directions should be independent random variables. Evidence for the close correlation and alignment of directions on different angular scales then appears to be a signature of broken statistical isotropy. Since we do not observe a significant amplitude of dipole modulation over similar angular scales, then the result of finding clustering of directions seems mysterious – it is hard to imagine a concrete (e.g., inflationary) model that would cluster the directions without affecting the amplitude of the modulation (but see Hansen et al. 2019). Nevertheless, regarding this as a purely empirical question, it is important to repeat the directionalclustering analysis and broaden its scope to include polarization.
Here, we repeat the analysis using the Planck 2018 data and extend it to include polarization measurements. The local power spectra, TT, TE, and EE, are estimated directly from maps of the temperature and Stokes Q and U parameters. Since the analysis is sensitive to differences between the noise properties of the data and simulations, we consider only crossspectra computed between the two HM or OE maps for each componentseparation method.
We adopt the same approach for the estimation of the dipole alignment as described in detail in PCIS15, a brief summary of which follows.
1. Local TT, TE and EE power spectra are estimated (using the MASTER approach, Hivon et al. 2002) directly from the T, Q, and U maps at N_{side} = 2048 for the 12 patches of the sky corresponding to the N_{side} = 1 HEALPix base pixels. Leakage between the E and B modes due to incomplete sky coverage is removed in the powerspectrum estimation during the inversion of the full TEB kernel. The same results would be expected using the E maps, but the extended masks required for these maps would increase the uncertainty of the results. The spectra are binned over various (even) bin sizes between Δℓ = 8 and Δℓ = 32^{14}.
2. For each powerspectrum multipole bin, an N_{side} = 1 HEALPix map with the local power distribution is constructed.
3. The bestfit dipole amplitude and direction is estimated from this map using inversevariance weighting, where the variance is determined from the local spectra computed from the simulations. As in previous papers (PCIS13; PCIS15), the fitted dipole amplitudes are found to be consistent with those determined from simulations.
4. A measure of the alignment of the different multipole blocks is then constructed. We use the mean of the cosine of the angles between all pairs of dipoles. This is essentially the Rayleigh statistic (RS) and we will refer to it as such, although it differs by not including any amplitude information. 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 because the RS as a function of ℓ_{max} is cumulative.
6. After calculating results for each bin size, we 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 RS values for each single multipole.
Since the test is based on the angular correlation between power dipoles (dipoles of the angular distribution of power as described above) and not on the absolute direction, the results should be unaffected by any preferred directions in the data caused by its noise properties, as long as these are matched by the noise properties of the simulations. Nevertheless, we have tested the uniformity of the estimated dipole directions in the simulated data. We find that the EE polarization directions are strongly influenced by the high noise level, with a strong bias towards the ecliptic plane. In order to reduce the noise bias on direction, we weight the Q and U maps by their inverse noise variance before estimating the local spectra. Figure 35 shows the histogram of the ecliptic latitudes of the power dipoles determined from simulated maps constructed in 100multipole bins. The horizontal black lines indicate the 95% confidence interval obtained from uniformly distributed directions. The coloured lines show the distribution of directions from temperature and EE polarization dipoles estimated on small and large scales, using inversevariance weighting in the latter case. Deviation from a uniform distribution is seen for the polarization directions, most notably for the smallscale EE results (magenta line). The ecliptic longitude values fall consistently within the 95% confidence interval determined from uniformly distributed directions for both temperature and polarization, although some modulation is seen.
Fig. 35. Distribution of directions for power dipoles estimated from independent 100multipole bins for simulated HM Commander maps with the common mask applied. The fraction of dipoles pointing towards a given ecliptic latitude direction is shown, where the data are binned by the cosine of the direction, since it is this quantity that is uniformly distributed on the sky for a Gaussian random field. The horizontal black lines show the expected 95% deviation taken from a uniform distribution. The coloured lines show the distributions for temperature at ℓ < 800 (red line), temperature at ℓ > 800 (green line), EE polarization at ℓ < 800 (blue line), and EE polarization at ℓ > 800 (magenta line). The temperature results are determined using uniform weighting, while the EE results have inversevariance weighting applied. 
The three panels on the left of Fig. 36 show the TT, TE, and EE dipole directions determined for fifteen successive 100multipole bins from the OE split of the Commander data with the OE common mask applied. This particular binning has been chosen for visualization purposes. The temperature results (top panel) are consistent with those of PCIS15 (although note that, in order to highlight the clustering, the plots in the current paper are rotated 180° about the Galactic northsouth axis, so that the Galactic longitude of l = 180° is at the centre of the image). The preferred lowℓ modulation direction determined from the temperature data in Sect. 7.2 is also indicated. The right panels of Fig. 36 presents the corresponding pvalues for the power asymmetry as a function of ℓ_{max}. The significance of the temperature alignment as a function of ℓ_{max} is consistent with earlier results up to ℓ_{max} ≈ 1000. We see that from ℓ_{max} ≈ 150 − 1000, the pvalues are below 1% for all values of the maximum multipole for all methods. However, from ℓ_{max} ≈ 1000, the pvalues increase rapidly. Figure 37 presents the equivalent results for the HM data set.
Fig. 36. Left panels: dipole directions for independent 100multipole bins of the local powerspectrum distribution from ℓ = 2 to 1500 in the Commander OE maps, with the OE common mask applied. The preferred directions for maps of specific multipole bins are colour coded according to the central value of the bin, ℓ_{central}, as shown in the colour bar. Note that the maps have been rotated about the Galactic northsouth axis, so that Galactic longitude l = 180° is in the centre of the map. The top panel shows the directions for the TT power spectrum, the middle panel for the TE power spectrum and the bottom panel for the EE power spectrum. In all panels, we also show the CMB dipole direction, the north ecliptic pole (NEP), the south ecliptic pole (SEP), and the preferred dipolar modulation axis (labelled as “lowℓ”) derived from the temperature data in Sect. 7.2. The average directions determined from the two multipole ranges ℓ = 2 − 300 and ℓ = 750 − 1500 are shown as blue and dark red (open) rings, respectively. The mean error on the derived directions that results from masking the data is 44° in the range ℓ < 800 and 62° in the range ℓ > 800 for temperature, but 78° in the range ℓ < 800 and 91° in the range ℓ > 800 for EE polarization. Right panels: derived pvalues for the angular clustering of the power distribution in OE maps as a function of ℓ_{max}, determined for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), based on FFP10 simulations using the OE common mask. These are the same colours used throughout the paper (e.g., see Fig. 5), while the grey line shows the Commander results when excluding the first 100 multipoles in the analysis. These 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, hence small pvalues would correspond to anomalously aligned dipole directions. The results shown here have been marginalized over bin sizes in the range Δℓ = 8 to Δℓ = 32. 
The application of the OE or HMunobserved pixel masks is necessary for the analysis of the componentseparated maps in order to avoid complications related to inpainting in these pixels. However, the Commander data set is the exception here, since it applies perpixel inversenoise weighting per frequency channel, so that unobserved pixels in a given channel are simply given zero weight in the parametric fits. Thus a valid CMB estimate is provided for such pixels, at the expense of higher noise. This is particularly relevant when comparing results to the 2015 analysis, since the number of unobserved pixels has increased significantly between the 2015 and 2018 data sets. Indeed, it seems apparent that the application of the unobserved pixel masks has significantly increased the error on direction for ℓ_{max} > 1000, resulting in a corresponding change in pvalues for these multipoles.
To test the change in significance for ℓ_{max} > 1000 due to the inclusion of the unobserved pixel masks, we investigate all simulations with low pvalues for ℓ_{max} > 1000 and check if a similar trend can be seen. As a tradeoff between ensuring a low pvalue and having a sufficient number of simulations to obtain reliable statistics, we select those simulations with a mean pvalue ⟨p⟩< 10% for ℓ > 1000. In Fig. 38 we show the 68% spread of these pvalues obtained from OE Commander simulations with the fullmission and OE common masks applied, and the HM Commander simulations with the HM common mask applied. The red solid lines show the corresponding results for the data. We see a similar trend in the simulations as observed in the data: the pvalues increase for smaller scales as a result of the higher uncertainty on direction caused by the unobserved pixel masks.
Fig. 38. Modulation results at high ℓ. The shaded bands show the 68% spread of pvalues for asymmetry taken from Commander simulations with a mean pvalue less than 10% for ℓ > 1000. The dark green band (dark green boundaries) represents results computed from the OE data set using the fullmission common mask, the grey band (with grey boundaries) to the same data set analysed with the OE common mask applied, and the light green band (with light green boundaries) corresponds to the HM data set analysed using the HM common mask. The red solid lines show the same three cases for the data: lower line, fullmission common mask; middle line, OE common mask; and upper line, HM common mask. 
In Fig. 39 we show the results for the Commander OE analysis when the fullmission common mask has been applied. We see that, for ℓ < 1000, the results are largely consistent with the previous analysis when the unobserved pixel mask was also applied. From ℓ ≈ 150 to 1150, the pvalues are consistently below 1% for all multipoles. This is in good agreement with the Commander results in PCIS15, although we note that the latter results were computed using the 2015 common mask that did not include unobserved pixels. For the EE polarization signal, some alignment seems to be indicated, reaching p < 1% at ℓ ≈ 150, which corresponds to the multipole range where the alignment in temperature also starts to be seen. The EE alignment will be discussed in more detail below.
Fig. 39. Left panels: as for Fig. 36, but for the Commander OE maps with the fullmission common mask applied. The mean error on the derived direction that results from masking the data is 39° in the range ℓ < 800 and 50° in the range ℓ > 800 for temperature, but 78° in the range ℓ < 800 and 91° in the range ℓ > 800 for EE polarization. Right panels: derived pvalues for the angular clustering of the power distribution in OE maps as a function of ℓ_{max}, determined for Commander (red line) based on simulations with the fullmission common mask applied. The grey line shows the Commander results when excluding the first 100 multipoles in the analysis, the black solid line shows results for the Commander HM split, and the black dashed line corresponds to the grey line for the Commander HM split, in all cases with the common mask applied. 
We note that for ℓ_{max} < 100 the temperature pvalues are not consistent with the detection of a lowℓ asymmetry/modulation, as seen in Sect. 7.1. However, over this ℓ range, there are very few bins and the variance of the RS might therefore be too high for the effect to be visible. We further test the multipole dependence of the alignment by restricting the analysis to multipoles above ℓ_{min} = 100. The grey lines in Fig. 39 show the Commander results in this case, and indicate that clustering at the p < 1% level is still found for temperature. The clustering significance can therefore not be solely attributed to largeangularscale features. Note also that the dip in pvalues for the EE polarization at ℓ ≈ 200 is still seen even with the restriction ℓ_{min} > 100.
It is also apparent from Fig. 39 that some pvalues are close to 100%, in particular for TE. A high pvalue means a low value for the RS statistic, and hence it seems worthwhile asking if such low values are unusual. In fact we find that scanning over the range of ℓ_{max}, the maximum pvalue of the RS statistic for the TE data is exceeded in 20–40% of the simulations (for the example of Commander, for the various masks and data splits), and hence does not appear to be anomalous.
In order to further study the alignment in EE polarization for ℓ < 300, we tested whether the directions of the EE dipoles in this range are correlated with the directions for the TT dipoles. Here we made a small change in the statistic as already described: in point 4 in the above description of the method, we instead use the mean of the cosine of the angles between all pairs of dipoles, where one dipole is always taken from EE and the other always from TT. In this way, the statistic measures the crosscorrelation between TT and EE directions. In Fig. 40 we show the pvalues for the crosscorrelation statistic as a function of multipole. Note that a strong correlation between the position of the TT and EE dipoles are detected for the same multipole range where the EE polarization dipoles (and TT dipoles) are clustered.
Fig. 40. Derived pvalues for the angular correlation of TT and EE dipole directions in Commander maps as a function of ℓ_{max}, determined using the HM split with the fullmission common mask (black), the HM split with the HM common mask (green), the OE map with the fullmission common mask (red), and the OE map with the OE common mask (blue), based on FFP10 simulations. The magenta line shows the corresponding correlation between TE and EE directions for the Commander HM split with the common mask. 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, hence small pvalues would correspond to anomalously aligned dipole directions. The results shown here have been marginalized over bin sizes in the range Δℓ = 8 to Δℓ = 32. 
We have seen that the EE directions appear clustered for ℓ_{max} ≈ 200. This is also the case for TT, and so we would like to test whether the two preferred directions are also related. What Fig. 40 shows is that not only are the TT and EE dipoles clustered among themselves, but that they also appear to be clustered towards the same direction. If TT and EE were independent, this would be highly unexpected, even if TT and EE were separately clustered. However, we know that the TE spectrum is nonzero, giving rise to a correlation between T and E. In order to test whether such a TT–EE directional correlation is expected in the case when both TT and EE are clustered individually, we perform the following test. We examine all simulations having a minimum pvalue of less than 1% for both TT and EE in overlapping multipole ranges (similar to what we observed in the data). Only two simulations have overlapping TT and EE multipole ranges using this criterion, but neither of them has a significant correlation between TT and EE directions. While this may seem to suggest that the effect in EE is not just due to the already known directional clustering in TT (and the TE correlation), it is hard to draw firm conclusions without many more simulations.
We investigate this further in Fig. 41 where we show the dipole directions for blocks of 10 multipoles for ℓ < 200 for TT, TE, and EE. There is some hint here of the correlation between TT and EE directions seen in Fig. 40. Note that in Fig. 40 this angular correlation appears stronger for the HM split than for the OE split; given the large errors on direction for polarization, differences in the noise properties (including systematic effects) for the HM and OE split may be responsible. The middle panel in Fig. 41 shows the TE directions. Clearly the TE directions are more scattered than TT and EE, as expected from Fig. 39. The directions for several multipole blocks in the range ℓ < 100 coincide quite well with the TT and EE directions, but most of the 10multipole blocks point in different directions. To further compare the angular clustering of TE dipoles with TT and EE dipoles, we construct the statistic measuring the correlation between TE and TT directions, as well as between TE and EE directions, in the same way as described above for the correlation between TT and EE. We find no significant correlations between the TT and TE directions, but the correlation between TE and EE again shows some sign of similar multipole ranges where EE is aligned and the TT and EE directions correlate. This is shown as a magenta line in Fig. 40. The correlation between TE and EE directions at ℓ ≈ 200 may arise as a result of the fact that TT and EE are aligned, as well as the fact that TT and EE directions are also correlated at exactly these multipoles. However, since none of the simulations show a high significance for all three of these statistics in one common multipole range for one given simulation, we are unable to test this with simulations. Nevertheless, we might have expected a similar correlation between TT and TE if this was really the case.
Fig. 41. Dipole directions for independent 10multipole bins of the local powerspectrum distribution from ℓ = 2 to 200 in the Commander HM maps with the fullmission common mask applied. This is a finer binning in order to investigate the directional clustering of different multipole ranges. Note that the maps have been rotated about the Galactic northsouth axis, such that Galactic longitude l = 180° is in the centre of the map. The top, middle, and bottom plots show maps based on the TT, TE, and EE spectra, respectively. In all panels, we also show the CMB dipole direction, the north ecliptic pole (NEP), the south ecliptic pole (SEP), and the preferred dipolar modulation axis (labelled as “lowℓ”) derived from the temperature data in Sect. 7.2. 
In PCIS15, we made a more detailed study of the total significance of the alignment based on the pvalues in figures like Fig. 39. Specifically, we compared the mean and minimum pvalues (over a range of ℓ_{max}) computed from the data with the value from simulations. Here we repeat this study of “global statistics” for selected cases of particular interest (see Bennett et al. 2011 and PCIS15 for discussion of how this is an attempt to assess the overall significance without making choices a posteriori). The results are shown in Table 26. For temperature, we see that none of the simulations have lower mean pvalues than the data. This is valid for most masks or foregroundsubtraction methods. Because the data have a minimum pvalue of 0 for all masks and methods, only an upper bound on significance can be set. Considering only ℓ_{min} > 100, the pvalues are still low.
Significance of the angular clustering of the power distribution.
Turning now to polarization, for TE, the numbers are completely consistent with simulations, while for EE, the results are unstable to the choice of masks and componentseparation methods. The mean pvalue is always consistent with simulations, which is not unexpected given the short multipole range where the pvalues are low (as discussed above). However, for the minimum pvalue, the lowest values seen in the multipole range ℓ = 100–300 are zero for some masks and some methods. This also occurs in some simulations and hence only an upper bound on the significance can be set. A similar conclusion can be made for the TT/EE angular correlation, but the results are in this case more stable with choice of mask and componentseparation method.
8. Conclusions
In this paper, we present a study of the statistical isotropy and Gaussianity of the CMB using the Planck 2018 temperature and polarization data. The Planck 2015 release essentially corresponded to the limit of our ability to probe CMB anomalies with temperature fluctuations alone. The use of largeangularscale polarization measurements enables largely independent tests of these peculiar features. In principle, this can reduce or eliminate the subjectivity and ambiguity in interpreting their statistical significance.
As in previous work, we follow a modelindependent approach and focus on nullhypothesis testing by calculating and reporting pvalues for a number of statistical tests. These tests are performed on maps of the CMB anisotropy that originate from the four componentseparation methods, Commander, NILC, SEVEM, and SMICA, described in Planck Collaboration IV (2020). For polarization studies, we consider both maps of the Stokes parameters, Q and U, and of the Emode signal generated using a novel method described in Appendix A. The consistency of the results determined from the componentseparated maps is an important indicator of the cosmological origin of any significant pvalues, or otherwise.
The temperature results are consistent with previous findings in PCIS13 and PCIS15. Specifically, the observed fluctuations are largely compatible with Gaussian statistics and statistical isotropy, with some indications of departures from the expectations of ΛCDM in a few cases. Such signatures are well known; thus, in this summary, we focus on the properties of the polarization data.
In Sect. 5, we examine aspects of the Gaussianity of the polarized CMB fluctuations. Tests of 1D moments, Npoint correlation functions, Minkowski functionals, peak statistics, and the oriented and unoriented stacking of peaks yield no indications of significant departures from Gaussianity. In addition, no evidence is found for a low variance of the polarized sky signal.
Section 6 provides an updated study of several previously known peculiarities. We find no evidence in the polarization data of a lack of largescale angular correlations, a hemispherical asymmetry in the behaviour of Npoint functions or peak distributions, a violation of pointparity symmetry, or a polarization signature associated with the Cold Spot.
In Sect. 7 we perform a series of tests searching for the signature in polarization of the wellknown largescale dipolar power asymmetry. Neither investigations using a variance estimator nor via ℓ to ℓ ± 1 mode coupling find strong evidence of this asymmetry. However, an interesting alignment of the preferred directions of the temperature and Emode dipolar modulation is found using the variance asymmetry estimator at a modest significance, depending on the componentseparated map in question. The modecoupling estimator indicates that the polarizationonly results show some apparent asymmetry over scales up to ℓ_{max} ≈ 250, a range that overlaps the scales of interest for the variance asymmetry. Similarly, an independent, but related, test of directionality finds suggestions of some alignment of directions in the EE polarization signal beginning at ℓ_{max} ≈ 150 and extending to ℓ_{max} ≈ 250.
There are some caveats worth pointing out here. Firstly, as described in Contreras et al. (2017), one could predict of order 30% chance that any of the pvalue dips would increase in significance when even statistically isotropic polarization data were added. Secondly, all of these dipolemodulation and hemisphericasymmetry statistics are just measuring slightly different weightings of the ℓ to ℓ ± 1 couplings on the same sky, and hence they cannot be considered to be independent of each other. Lastly, we are only testing phenomenological models here, rather than physical modulation models where there is a prediction for how scales in temperature and in polarization might be separately modulated. Hence it is unclear if the hints of EE–TT dipolemodulation alignment are what we would expect if there was a physical mechanism responsible for some modulation. Whether the hint of alignment between the temperature and the polarization dipolar power asymmetry is more than a coincidence can only be addressed once new data are available from forthcoming largescale and lownoise polarization experiments.
A notable feature of all of the polarization analyses is the variation in pvalues for a given test between the four componentseparated maps. This is a consequence of the fact that different componentseparation methods respond to noise and residual systematic effects in different ways. However, it may also indicate an incomplete understanding of the noise properties of the data, both in terms of amplitude and correlations between angular scales. This should not be considered surprising, given that Planck was not optimized for polariation measurements. Although remarkable progress has been made in reducing the systematic effects that contaminated the 2015 polarization maps on large angular scales, particularly for the HFI instrument (Planck Collaboration Int. XLVI 2016; Planck Collaboration III 2020), thus allowing more robust measurement of the opticaldepthtoreionization τ, residual systematics, and our ability to simulate them, can limit the kind of statistical tests of nonGaussianity and isotropy that can be applied to the data. Nevertheless, a detailed set of null tests applied to the maps indicates that these issues do not dominate the analysis on intermediate and large angular scales, particularly for the statistical tests presented in this paper.
Future experiments that can measure the cosmological Emodes at the cosmicvariance limit are required in order to unambiguously test for the presence of anomalies in the polarized sky. However, given the amplitude of the effects seen in the Planck temperature data, it may still remain difficult to claim high significance (> 3σ) detections in polarization (although detailed forecasts related to this are highly model dependent). Nevertheless, this should not prevent us from undertaking such searches, since any detection of anomalies in the polarized sky signal will inevitably 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).
These are pixels that were either never seen by any of the bolometers present at a given frequency, or for which the polarization angle coverage is too poor to support a reliable decomposition into the three Stokes parameters. Note that the number of unobserved pixels has increased significantly between the 2015 and 2018 data sets, due to a change in the condition number threshold at the mapmaking stage.
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, where both effects are aligned with the CMB dipole (Challinor & van Leeuwen 2002; Planck Collaboration XXVII 2014). Both contributions are present in the FFP10 simulations.
One of the more important tests in the context of inflationary cosmology is related to the analysis of the bispectrum. This is explored thoroughly in Planck Collaboration IX (2020), and is therefore not discussed further in this paper.
FFP8 is the name given to the previous version of full focalplane simulations utilized for the analysis of the 2015 Planck maps (Planck Collaboration XII 2016).
Referred to as the “Solar dipole” in accompanying Planck papers (see also Sect. 2.1 of Planck Collaboration I 2020).
In order to maintain consistency with previous Planck Collaboration papers, some of the figures will show results using 100multipole bins. These 100ℓ bins were constructed from the 16ℓ bins using the prescription described in PCIS13 and Hansen et al. (2009). Although this means that there are not exactly 100 multipoles in each bin, they will nevertheless be referred to as “100ℓ bins”. Note that these bins are used for illustrative purposes in the figures only, while the analysis always uses smaller bins.
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. This work has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement numbers 687312, 776282, and 772253.
References
 Abrial, P., Moudden, Y., Starck, J. L., et al. 2008, Stat. Meth., 5, 289 [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., & Shafieloo, A. 2017, ArXiv eprints [arXiv:1710.00580] [Google Scholar]
 Baldi, P., Kerkyacharian, G., Marinucci, D., & Picard, D. 2009, Ann. Stat., 37, 1150 [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]
 Bielewicz, P., Gorski, K. M., & Banday, A. J. 2004, MNRAS, 355, 1283 [NASA ADS] [CrossRef] [Google Scholar]
 Bond, J. R., & Efstathiou, G. 1987, MNRAS, 226, 655 [NASA ADS] [Google Scholar]
 Bowyer, J., Jaffe, A. H., & Novikov, D. I. 2011, ArXiv eprints [arXiv:1101.0520] [Google Scholar]
 Brandt, A. 1977, Math. Comput., 31, 333 [CrossRef] [MathSciNet] [Google Scholar]
 Bucher, M., & Louis, T. 2012, MNRAS, 424, 1694 [NASA ADS] [CrossRef] [Google Scholar]
 Bunn, E. F. 2011, Phys. Rev. D, 83, 083003 [NASA ADS] [CrossRef] [Google Scholar]
 Bunn, E. F., & Wandelt, B. 2017, Phys. Rev. D, 96, 043523 [NASA ADS] [CrossRef] [Google Scholar]
 Bunn, E. F., Zaldarriaga, M., Tegmark, M., & de OliveiraCosta, A. 2003, Phys. Rev. D, 67, 023501 [NASA ADS] [CrossRef] [Google Scholar]
 Cammarota, V., & Marinucci, D. 2016, ArXiv eprints [arXiv:1603.09588] [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]
 Chiang, L.Y. 2018, ApJ, 861, 8 [CrossRef] [Google Scholar]
 Chingangbam, P., Ganesan, V., Yogendran, K. P., & Park, C. 2017, Phys. Lett. B, 771, 67 [CrossRef] [Google Scholar]
 Contreras, D., Zibin, J. P., Scott, D., Banday, A. J., & Górski, K. M. 2017, Phys. Rev. D, 96, 123522 [NASA ADS] [CrossRef] [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. 2013, MNRAS, 434, 3590 [CrossRef] [Google Scholar]
 Copi, C. J., Huterer, D., Schwarz, D. J., & Starkman, G. D. 2015, MNRAS, 451, 2978 [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., Vielva, P., MartínezGonzález, E., & Barreiro, R. B. 2011, MNRAS, 412, 2383 [NASA ADS] [CrossRef] [Google Scholar]
 de OliveiraCosta, A., & Tegmark, M. 2006, Phys. Rev. D, 74, 023005 [NASA ADS] [CrossRef] [Google Scholar]
 Ducout, A., Bouchet, F. R., Colombi, S., Pogosyan, D., & Prunet, S. 2012, MNRAS, 423 [Google Scholar]
 Efstathiou, G. 2004, MNRAS, 348, 885 [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., Marinucci, D., Hansen, F., & Maino, D. 2015, Phys. Rev. D, 91, 063501 [NASA ADS] [CrossRef] [Google Scholar]
 Feeney, S. M., Peiris, H. V., & Pontzen, A. 2011, Phys. Rev. D, 84, 103002 [NASA ADS] [CrossRef] [Google Scholar]
 FernándezCobos, R., Vielva, P., MartínezGonzález, E., Tucci, M., & Cruz, M. 2013, MNRAS, 435, 3096 [NASA ADS] [CrossRef] [Google Scholar]
 Ferreira, P. G., & Magueijo, J. 1997, Phys. Rev. D, 56, 4578 [NASA ADS] [CrossRef] [Google Scholar]
 Frommert, M., & Enßlin, T. A. 2010, MNRAS, 403, 1739 [NASA ADS] [CrossRef] [Google Scholar]
 Gay, C., Pichon, C., & Pogosyan, D. 2012, Phys. Rev. D, 85, 023011 [NASA ADS] [CrossRef] [Google Scholar]
 Gaztañaga, E., & Scoccimarro, R. 2005, MNRAS, 361, 824 [NASA ADS] [CrossRef] [Google Scholar]
 Ghosh, S., & Jain, P. 2018, ArXiv eprints [arXiv:1807.02359v1] [Google Scholar]
 Gjerløw, E., Eriksen, H. K., Banday, A. J., Górski, K. M., & Lilje, P. B. 2010, ApJ, 710, 689 [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]
 Gruetjen, H. F., Fergusson, J. R., Liguori, M., & Shellard, E. P. S. 2017, Phys. Rev. D, 95, 043532 [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, JCAP, 7, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Gruppuso, A., Kitazawa, N., Lattanzi, M., et al. 2018, Phys. Dark Universe, 20, 49 [CrossRef] [Google Scholar]
 Hansen, F. K., Banday, A. J., Górski, K. M., Eriksen, H. K., & Lilje, P. B. 2009, ApJ, 704, 1448 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, F. K., Trombetti, T., Bartolo, N., et al. 2019, A&A, 626, A13 [CrossRef] [EDP Sciences] [Google Scholar]
 Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Hou, Z., Banday, A. J., & Górski, K. M. 2009, MNRAS, 396, 1273 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., & White, M. 1997, New Astron., 2, 323 [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., & Naselsky, P. 2010, ApJ, 714, L265 [NASA ADS] [CrossRef] [Google Scholar]
 Krylov, A. N. 1931, Otdel. mat. i estest. nauk, VII, 491 [Google Scholar]
 Lanczos, C. 1950, J. Res. Nat. Bur. Std., 45, 255 [CrossRef] [Google Scholar]
 Land, K., & Magueijo, J. 2005a, Phys. Rev. D, 72, 101302 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Land, K., & Magueijo, J. 2005b, Phys. Rev. Lett., 95, 071301 [NASA ADS] [CrossRef] [PubMed] [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]
 Lewis, A., Challinor, A., & Turok, N. 2002, Phys. Rev. D, 65, 023505 [NASA ADS] [CrossRef] [Google Scholar]
 MarcosCaballero, A., FernándezCobos, R., MartínezGonzález, E., & Vielva, P. 2016, JCAP, 4, 058 [CrossRef] [Google Scholar]
 MarcosCaballero, A., MartínezGonzález, E., & Vielva, P. 2017a, JCAP, 5, 023 [CrossRef] [Google Scholar]
 MarcosCaballero, A., MartínezGonzález, E., & Vielva, P. 2017b, JCAP, 2, 026 [CrossRef] [Google Scholar]
 Marinucci, D., Pietrobon, D., Balbi, A., et al. 2008, MNRAS, 383, 539 [NASA ADS] [CrossRef] [Google Scholar]
 Matsubara, T. 2010, Phys. Rev. D, 81, 083505 [NASA ADS] [CrossRef] [Google Scholar]
 Mitra, S., Rocha, G., Górski, K. M., et al. 2011, ApJS, 193, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Molinari, D., Gruppuso, A., Polenta, G., et al. 2014, MNRAS, 440, 957 [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]
 Muir, J., Adhikari, S., & Huterer, D. 2018, Phys. Rev. D, 98, 023521 [NASA ADS] [CrossRef] [Google Scholar]
 Nishizawa, A. J., & Inoue, K. T. 2016, MNRAS, 462, 588 [CrossRef] [Google Scholar]
 Patra, M., & Karttunen, M. 2006, Differ. Equ., 22, 936 [Google Scholar]
 Planck Collaboration XIX. 2014, A&A, 571, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2014, A&A, 571, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIII. 2014, A&A, 571, A23 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Planck Collaboration XXIV. 2014, A&A, 571, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2014, A&A, 571, A27 [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 XII. 2016, A&A, 594, A12 [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 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 I. 2020, A&A, 641, A1 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2020, A&A, 641, A2 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2020, A&A, 641, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2020, A&A, 641, A4 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2020, A&A, 641, A5 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2020, A&A, 641, A7 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2020, A&A, 641, A8 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2020, A&A, 641, A9 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2020, A&A, 641, A10 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2020, A&A, 641, A11 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2020, A&A, 641, A12 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLVI. 2016, A&A, 596, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLIX. 2016, A&A, 596, A110 [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]
 Rocha, G., Contaldi, C. R., Bond, J. R., & Gorski, K. M. 2011, MNRAS, 414, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, K. M. 2006, Phys. Rev. D, 74, 083002 [NASA ADS] [CrossRef] [Google Scholar]
 Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Starck, J. L., Fadili, M. J., & Rassat, A. 2013, A&A, 550, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vanmarcke, E. 1983, in Random Fields, ed. E. Vanmarcke (Cambridge, MA, USA: The MIT Press), 372 [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., MartínezGonzález, E., Cruz, M., Barreiro, R. B., & Tucci, M. 2011, MNRAS, 410, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Zaldarriaga, M., & Seljak, U. 1997, Phys. Rev. D, 55, 1830 [Google Scholar]
 Zibin, J. P., & Contreras, D. 2017, Phys. Rev. D, 95, 063011 [CrossRef] [Google Scholar]
Appendix A: E and Bmode map reconstruction
A.1. Introduction and notation
Polarization of the CMB is usually measured on the sky in terms of the Stokes parameters Q and U (see e.g., Hu & White 1997). However, these are not scalar quantities and therefore not rotationally invariant. Nevertheless, scalar E and pseudoscalar B maps can be determined from the measured quantities, as described in Sect. 2. Such maps offer definite advantages for any mapbased analyses, such as those presented in the main body of this paper. Since they are generated via the application of nonlocal spherical harmonic transformations, the E and B estimators are nontrivial to construct in the typical case of partial sky coverage, resulting from the need to mask out strong foreground contributions in the data.
Alternative sets of related scalars are occasionally used in the literature (Zaldarriaga & Seljak 1997; Bunn et al. 2003), defined by
They are related to the standard E and B estimators in harmonic space by a factor of
which corresponds to an application of the biLaplacian operator such that ℰ, ℬ = −∇^{2}(∇^{2} + 2) ψ_{E, B} in real space. Unlike E and B, ℰ and ℬ maps can be derived from Q and U by a (local) secondderivative operator, although the noise power at high ℓ is then significantly enhanced.
The central problem for the reconstruction of E and Bmode maps from partial sky coverage is that neither spin0 nor spin2 spherical harmonics are orthogonal under the masked inner product,
and thus mode mixing generally occurs. Here, we introduce a vector notation for the polarization map, p ≡ (Q_{i}, U_{i}), and the masking operator, M ≡ diag(M_{i}), which multiplies an individual map pixel i by a mask value M_{i}. To keep the notation concise, we will denote all other linear operators similarly, e.g., E and B correspond to the projection of the polarization map onto its E and Bmode components, although numerically this would be implemented using spherical harmonic transforms rather than matrix multiplications.
Mode mixing is a wellknown problem for the estimation of power spectra on a masked sky (see, for example, Rocha et al. 2011). In this case, the effect of masking results in an estimated power spectrum that is a linear combination of the fullsky quantities. For an isotropic Gaussian random field with an ensemble average spectrum ⟨C_{ℓ}⟩, the masked sky mode averages are given by
where the coupling matrices K_{ℓ1ℓ2} depend on the method used, as illustrated in Fig. A.1.
Fig. A.1. Modecoupling matrices up to ℓ_{max} = 512 for the common polarization mask with point sources at N_{side} = 1024 for masking (top), inpainting (middle), and purifiedinpainting (bottom) methods. The coloured shading represents normalized matrix elements on a logarithmic scale, with values spanning from 10^{−8} (deep blue) to 1 (dark brown). 
For reconstruction of lowℓ multipoles, maximumlikelihood estimators are widely used (Efstathiou 2004; Bielewicz et al. 2004; de OliveiraCosta & Tegmark 2006; Feeney et al. 2011). Filling in the missing data in the CMB maps can be done using various statistical priors (Abrial et al. 2008; Bucher & Louis 2012; Starck et al. 2013; Nishizawa & Inoue 2016), in particularly using constrained Gaussian realizations. Methods targeting decomposition of polarization into pure E and B modes plus ambiguous components have already been developed (Bunn et al. 2003; Bunn 2011; Bunn & Wandelt 2017). Some of the approaches discussed in the literature require solving large linear algebra problems (typically via iterative solvers), and could be expensive on highresolution maps. With the large number of simulations that need to be processed for Planck data analysis, numerical performance becomes a very important issue.
We consider three direct approaches to the computation of E and Bmode maps in the case of incomplete sky coverage in Appendices A.2, A.3, and A.4. The suitability of these approaches for our purposes depends upon the uniformity of the reconstruction, since the methodspecific residuals are generally quite inhomogeneous and dependent on the mask.
A.2. Masking
The simplest method to implement is the direct computation of E and Bmodes from the Q and U data after the application of a mask that zeros the problematic pixels. As a consequence, E and Bmode mixing does result, with mode coupling matrices expressible analytically in terms of Wigner3j symbols and the power spectrum of the mask defined by a window function :
The resultant modecoupling matrices are symmetric and are shown in Fig. A.1 for the common polarization mask at N_{side} = 1024 resolution. In practice, it is often faster to evaluate these matrices using Monte Carlo methods rather than explicit summation involving Wigner symbols.
As an alternative, if one masks ℰ and ℬ maps instead, there is no mode mixing between E and B. Unfortunately, secondorder derivative operators enhance the noise power, which results in large artefacts in the reconstruction due to hightolowℓ modecoupling in the masking operator, unless extremely strong apodization is applied to the mask (as described in Smith 2006).
A.3. Simple inpainting
A.3.1. Overview
The application of a diffusive inpainting procedure to the masked pixels of input sky maps has proven to be a satisfactory approach to handle incomplete sky coverage when searching for evidence of primordial nonGaussianity in the Planck temperature and polarization data (Planck Collaboration XXIV 2014; Planck Collaboration XVII 2016). To be more explicit, the procedure works by replacing each masked pixel by the average of all nearestneighbour pixels, then the process is repeated over a large number of iterations. Such an approach is straightforward to implement, but the convergence rate of the inpainted solution is slow for the largest scales.
To address this, we adopt slightly improved finitedifference “Laplacian stencils”, as detailed in Appendix A.3.2, and we further develop the details of the finitedifferencing multigrid approach in Appendix A.3.3.
One particular aspect of our improved method is that when computing the average over the nearest neighbours of a given masked pixel, their contributions are Gaussianweighted by the distance to their positions in the tangent plane via gnomonic projection. In addition, basis orientation differences between the Q and U polarization components are properly projected (via parallel transport in the tangent plane) to form a tensor Laplacian stencil. We improve on the speed of the inpainting algorithm by noting that an infinite number of relaxation iterations converges to the solution of an elliptical Laplace equation ∇^{2}T = 0 with the Dirichlet boundary conditions given by the unmasked pixels; this can be solved in O(n log n) operations using the standard multigrid methods of Brandt (1977) adapted to the spherical HEALPix pixelization, as described in Appendix A.3.3.
Modecoupling matrices for inpainted temperature maps have been presented in Gruetjen et al. (2017). Here, we extend the results to polarization, as shown in Fig. A.1, where the matrices have been computed for the common polarization mask at N_{side} = 1024 resolution. Unlike for the case of masking, the inpainting modecoupling matrices are not symmetric and result in excellent suppression of the lowtohighℓ mode mixing, at the expense of increased hightolowℓ mode mixing. Similarly to the case of pure ℰ and ℬ mode masking, hightolowℓ mixing renders the inpainting of Q and U maps susceptible to the transfer of noise artefacts into lowℓ patterns, albeit to a lesser extent.
A.3.2. Finitedifference stencils
Finitedifference approximations of differential operators are nontrivial to evaluate on the HEALPix grid, especially for polarization. In this section we discuss first and secondorderaccurate stencils for the Laplace operator using finite differences. In general, they can be represented as a weighted sum of the intensity and polarization for the pixel and its nearest neighbours:
with real c_{i} for scalar and complex for tensor values, and being the average area of a pixel. We will derive the stencil weights in a flatsky approximation by projecting onto a tangent plane through the pixel where the derivative operator is being evaluated.
A point on the unit sphere with spherical coordinates θ and ϕ has a local set of three orthonormal vectors,
associated with it, which form a basis on a tangent plane and a normal direction . Nearby points can be projected onto a selected tangent plane via gnomonic projection, which maps the radius vector r into
and thus introduces local coordinates on a tangent plane,
Secondorderaccurate discretizations of the Laplacian operator on rectangular grids are well known and easy to derive (Patra & Karttunen 2006); however, HEALPix pixels are placed differently, as is illustrated in Fig. A.2. Deformation of the grid is never small and does not scale down with increasing N_{side}. In addition, differences in orientation of the local bases are never small around the poles (as is obvious from the left column of Fig. A.2), and care must be taken in the parallel transport of the polarization tensor represented by the Stokes parameters Q and U.
Fig. A.2. Examples of the finitedifference stencils for HEALPix pixels with eight (left and centre) and seven (right) neighbours for N_{side} = 8 (top) and N_{side} = 1024 (bottom) in ring ordering (see also Bowyer et al. 2011). Red circles represent positions of HEALPix pixel centres in a gnomonic projection onto a plane tangent to the central pixel (i.e., looking straight down at the tangent plane), with the dotted grid aligned with local and directions, illustrating the average pixel pitch h = (π/3)^{1/2}/N_{side}. The white bars represent the directions of polarization, specifically the direction of the polarization ellipse for the +Q polarization mode. Single and doublearrow vectors show projections of and directions, respectively, for neighbouring pixels onto the tangent plane. The black cross corresponds to the average pixel position, while the blue dotted ellipse represents the pixel position covariance. Although for some positions on the sky, the polarization directions are aligned, this is not at all true near the poles (pixel 0); hence just adding Q and U does not make sense. Additionally we can see that because the grid is distorted, secondorder finitedifference schemes need more than just nearest neighbours to work. 
Under rotation of the basis used to define the Stokes parameters Q and U, i.e.,
their values transform as
In the HEALPix polarization convention, the Stokes parameters Q and U of a pixel are always defined with respect to the basis of the pixel, and must be rotated to the basis when projecting onto a tangent plane. The appropriate angle of rotation can be computed as
with transformation of Stokes parameters most conveniently implemented as a complex phase rotation
A linearorder shift of the average pixel position in the HEALPix grid breaks the symmetry of the local Taylor expansion, and there is a unique secondorderaccurate nearestneighbour discretization for the Laplace operator, as opposed to a oneparameter family on the rectangular grid. Unfortunately, it turns out to be unconditionally unstable for diffusiontype problems, so we will not discuss it here. Instead, we will use an approximate firstorder stencil based on the isotropic weighting c_{i} ≡ c(ρ_{i}/h), with coefficients normalized by ℒ[const.] = 0 and ℒ[ρ^{2}] = 4, which leads to
In previous studies (Planck Collaboration XXIV 2014; Planck Collaboration XVII 2016) equal weighting was widely used, but discretization residuals can be significantly improved at little expense by using Gaussian weighting, i.e.,
where the width σ can be tuned. We chose it to be σ = 1.61, which gives near perfect residual cancelation at the poles. This is the scalar Laplacian stencil we will adopt in what follows, while the complex Laplacian stencil for polarization is defined by , as explained above.
A.3.3. Multigrid methods
Inpainting a map ϕ can be viewed as a diffusive flow ∂_{t}ϕ = ∇^{2}ϕ applied to the masked areas, subject to the Dirichlet boundary conditions provided by the unmasked data. A firstorder forward in time discretization of the flow equation updates a pixel according to
using the weighted sum of itself and its neighbours, where α = δt/h^{2}, and the largest time step δt that can be taken is determined by CourantLewy stability analysis. The scheme is almost guaranteed to be unstable if the coefficient (1 + αc_{0}) becomes negative, so in practice the fastest diffusion is often achieved by replacing a pixel by a weighted sum of its neighbours:
N successive iterations will diffuse the solution across roughly N^{1/2} pixels, so running diffusion flow like this for full convergence to the static solution ∇^{2}ϕ = 0 is very expensive on highresolution maps. In previous studies (Planck Collaboration XXIV 2014; Planck Collaboration XVII 2016), a finite number of steps like Eq. (A.17) was applied to inpaint the masked areas.
Convergence of diffusion flow can be accelerated if one notes that a coarser spatial grid allows bigger time steps, and thus the solution of a linear elliptic boundary value problem on a grid with spacing h,
can be found iteratively from the approximate solution by downgrading the residual
to a coarser grid with spacing 2h, thus obtaining the coarsegrid correction
which can be upgraded back to a fine grid with spacing h and used to improve the solution, i.e.,
as detailed in Brandt (1977). Recursive application of coarsegrid correction (Eq. (A.20)) interlaced with diffusion steps (Eq. (A.17)), as illustrated in Fig. A.3, achieves convergence to the static solution in O(log N_{side}) iterations, most of which are on coarser grids, and thus very fast. This is the method we use to diffusively inpaint intensity and polarization maps in this paper.
Fig. A.3. Multigrid inpainting schedule for an N_{side} = 128 map. Filled nodes represent successive iterations of diffusion steps (Eq. (A.17)), downward strokes represent calculation of the residual (Eq. (A.19)) and its downgrade to the coarser grid, while upward strokes represent the upgrade of the correction (Eq. (A.20)) to the finer grid, and solid arrows represent the merge of the correction into the solution (Eq. (A.21)). Eight diffusion steps are used at all grid levels except the lowest, where 64 steps are used, which is enough to find the static solution there. 
The algorithm proceeds by first constructing the multigrid structure, which will contain temporary maps and residuals to be corrected. The inpainting mask is recursively degraded, and Laplacian stencils are precomputed for pixels to be inpainted at all grid levels. Note that only the strict interior of the masked region should be inpainted at coarse levels to avoid boundary effects that degrade the convergence rate. The finest grid is initialized with the map to be inpainted, while the coarse grids will contain corrections and are initialized to zero. The inpainting is carried out by repeated application of the recursive “wstroke”, as illustrated in Fig. A.3 for an N_{side} = 128 example. The wstroke takes a small number of diffusion steps (Eq. (A.17)) on the region to be inpainted, computes and downgrades the residual (Eq. (A.19)), recursively calls itself on a coarse grid twice to obtain the correction (Eq. (A.20)), upgrades the correction and merges it with the solution (Eq. (A.21)), then takes a small number of diffusion steps (Eq. (A.17)) again. The entire pattern is repeated until the desired convergence accuracy is reached. For the masks used in temperature and polarization analysis in this paper, 18 iterations of the wstroke are enough to reach the static solution to double precision. This requires only 288 diffusion steps at full resolution, and is substantially faster than the 2000 iterations at full resolution that were used in Planck Collaboration XXIV (2014). Improved performance of the inpainting routine allows for further applications and testing to be carried out with the same computing resources.
A.4. Purified inpainting
A.4.1. Method description
Here, we demonstrate that the advantages of masking and inpainting can be combined into a single method.
The motivation for the approach is to construct a fullsky polarization map, , from a masked one, M ⋅ p, such that it agrees identically with p on the unmasked portion of the sky, and assigns modes either through strict attribution to the E and B subspaces of the fullsky map, or “ambiguous” attribution where some mode mixing is allowed:
Note that the polarization map components e and b are not pure in the sense of being orthogonal to the entire E and Bmode linear spaces on the masked sky, as in Bunn et al. (2003); that requirement results in a large linear algebra problem, which is expensive to solve, although efficient methods for that have been developed recently (Bunn & Wandelt 2017). Instead, their purpose is to “purify” the ambiguous mode map a, ensuring that most of the pure modes are projected out, thus minimizing the power leakage between E and B in the ambiguous mode.
Initially, the entire polarization map p is assigned to an ambiguous mode:
Then the method proceeds by peeling off nonambiguous E and B modes one by one through the explicit construction of a Krylov subspace (Krylov 1931) generated by the inverse biLaplacian, K, acting on the masked maps. The obvious starting point that contains most of the CMB power is the Emode projection of the masked polarization map:
To prevent the constructed Krylov subspace basis from becoming degenerate, we use Lanczos biorthogonalization (Lanczos 1950). To do this we keep two recent normalized basis vectors, u and v, initially set to
and project them out (from the next Emode generated using the inverse biLaplacian operator, K) to obtain a new mode from the previous ones via
Once the new mode is constructed, it is normalized and the most recent basis set is updated:
Lanczos biorthogonalization guarantees that the dot product of the constructed basis vectors ⟨v_{i}, v_{j}⟩ forms a tridiagonal matrix, and thus avoids any stability issues associated with the basis vectors becoming nearly linearly dependent. Whenever a new basis vector, v, is constructed following Eqs. (A.25) or (A.27), it is projected out from the ambiguous mode map and assigned to the Emode map:
Once a sufficient number of E modes are extracted from the map, B modes are peeled off next in exactly the same way, via the Bmode projection operator, B, instead of E. The Krylov subspace construction (Eqs. (A.24) through (A.27)) can be restarted several times on the remaining ambiguous mode map, but it reaches a point of diminishing returns rather quickly. The exact number of Krylov modes to peel depends on a compromise between performance and quality of reconstruction, and is discussed in the next section.
After the ambiguous modes are purified, they still need to be reconstructed on the full sky. The correct approach is to inpaint the masked region by solving the elliptical partial differential equation, ∇^{2}a = 0, in the mask interior, subject to fixed values of a at the mask boundary, which minimizes the total extra amount of power introduced to the reconstructed ambiguousmode map.
Inpainting is a linear operation, and can be represented by a matrix operator F. A suitable method for inpainting highresolution maps is the multigrid approach, as discussed above, which has computational costs of 𝒪(n log n). The final reconstructed polarization map is generated by inpainting the remaining ambiguous modes as
A.4.2. Performance considerations
Figure A.1 presents the coupling matrices corresponding to the purified inpainting of the polarization data after application of the common polarization mask at N_{side} = 1024. The temperature maps are inpainted following the usual multigrid approach. The lowtohighℓ modemixing is much improved over the equivalent case with masking alone, while the increase in the hightolowℓ modemixing is substantially less than if a simple inpainting approach were applied. Purified inpainting therefore represents a good compromise for the E and Bmode reconstruction of CMB polarization maps.
The computational costs of the purifiedinpainting method are dominated by the two roundtrip fullresolution spherical transforms required for the initial projection (Eq. (A.24)) of the E and B modes, and the multigrid inpainting of the ambiguous modes. The Krylov subspace construction does eventually produce a complete basis for the E and Bmode subspaces, but is too expensive to run to completion for highresolution maps. An investigation of the computational cost versus quality of the reconstruction suggests that a drastic truncation works well. In fact, we extract only 32 modes in the code used here. The Krylov basis maps generated by an inverse biLaplacian operator generally have very red spectra, due to the inverse (ℓ−1)ℓ(ℓ + 1)(ℓ+2) factor in Eq. (A.2), hence do not require spherical transforms at full ℓ_{max}.
A.5. Reconstruction residuals and confidence mask
The accuracy of reconstruction of E and Bmode maps for a given method can be evaluated from realistic MC simulations of the CMB signal plus noise, where the true fullsky E^{*} and B^{*} maps are known, and can be directly compared to the reconstructed ones and . Various metrics of performance can be assigned to the masked residual maps and . Figure A.4 compares the residual maps for one SMICA realization from the FFP10 simulations, as reconstructed via each of the masking, simple inpainting, and purifiedinpainting methods. The latter seems to perform on par or better than other direct reconstruction methods published so far, and is competitive with the maximumlikelihood estimators, at substantially lower computational cost. Apodization and inpainting methods do not perform as well because of the E and Bmode mixing arising from the mask, while the pure projection method of Smith (2006) actually does much worse as a consequence of the highℓtolowℓ aliasing of the noise power. The purifiedinpainting method, presented here, offers a balance between these two considerations, and yields an order of magnitude lower residuals for N_{side} = 1024 Planck maps.
Fig. A.4. E and Bmode reconstruction residuals, shown in the left and right columns respectively, for the methods of mask multiplication (top), simple inpainting (middle), and purified inpainting (bottom). Residuals are shown for a single realization of a SMICA componentseparated simulation containing both CMB signal and noise at a resolution N_{side} = 1024. The grey area represents the common polarization mask applied during the reconstruction. 
Extending the mask beyond that used in the reconstruction process itself can help to reduce the residuals levels yet further, in particular given that these residuals tend to be localized near the mask boundary. To systematically define the optimal confidence mask, we performed purifiedinpainting reconstructions for the FFP10 CMBplusnoise realizations, as propagated through all four Planck component separation pipelines and subsequently provided at a number of resolutions, then evaluated the rms signal in the residual maps. Combined reconstruction residuals for all four componentseparation pipelines are shown in Fig. A.5, along with the confidence mask that admits 65% of the sky for further analysis, as derived from the isolevels of the total rms reconstruction residual δE^{2} + δB^{2} smoothed by an 80′ FWHM Gaussian beam. Figure A.6 shows the pseudoC_{ℓ} anisotropy power spectra of masked residuals for reconstructed E and B modes for all four componentseparation pipelines.
Fig. A.5. Rms of the combined reconstruction residuals for E and B modes (left and right columns, respectively) determined from the purified inpainting of the FFP10 simulations for all four componentseparation methods (i.e., the average of the square of each of the four residuals). The grey area represents the common polarization mask, while the semitransparent grey area indicates an expansion of the confidence mask, which together admit 65% of the sky for further analysis. The mask increment is determined from thresholding the total residual δE^{2} + δB^{2}, smoothed by an 80′ FWHM Gaussian. The rms reconstruction accuracy is better than 0.5 μK, with the largest deviations mainly localized near the mask boundary. 
Fig. A.6. PseudoC_{ℓ} spectra of residuals for E and B modes (left and right columns, respectively) determined from the purified inpainting of the FFP10 simulations for all four componentseparation methods. Solid curves show averaged spectra, with red corresponding to the residual masked with the reconstruction mask and blue corresponding to the residual masked with the confidence mask. Semitransparent areas filled with grey show 68%, 95%, and minimumtomaximum bounds for individual realizations. 
The choice of the confidence mask threshold is motivated by Fig. A.7, which shows the maximal rms reconstruction residuals versus the sky fraction. Note that a slight reduction in the sky fraction admitted by the confidence mask results in more than an order of magnitude decrease in the residuals. Smoothing the residuals in order to simplify the detailed geometry of the confidence mask results in a slight increase of the maximal residual, but is nevertheless beneficial for many analyses, given the dependence of the mode coupling on this structure. The final confidence mask is selected, such that the rms reconstruction residuals are less than 0.5 μK, significantly below the cosmological Emode signal, thus allowing sensitive tests to be undertaken concerning the isotropy and statistical properties of the Planck polarization data.
Fig. A.7. Maximal rms of reconstruction residuals for E modes (green), B modes (orange), and combined E and B modes (blue) determined from purifiedinpainting simulations as a function of sky fraction. The red circle and the thin dashed lines show the reconstruction quality within the common polarization mask used in reconstruction (77.7% of the sky admitted, maximum rms combined residual of 9.0 μK). The blue star and the thin solid lines show the reconstruction quality within the confidence mask (64.8% of the sky admitted, with a maximum rms combined residual of 0.5 μK). 
All Tables
Cosmological parameters for the FFP10 simulations, used to make the simulated maps in this paper, and throughout the Planck 2018 papers.
Probabilities of obtaining values for the χ^{2} statistic of the Npoint functions determined from the Planck fiducial ΛCDM model at least as large as those obtained from the Commander, NILC, SEVEM, and SMICA temperature and polarization (Q and U) maps at N_{side} = 64 resolution.
Probabilities of obtaining values for the χ^{2} statistic of the Npoint functions determined from the Planck fiducial ΛCDM model at least as large as the those obtained from the Commander, NILC, SEVEM, and SMICA temperature and polarization (Emode) maps at N_{side} = 64 resolution.
Probability as a function of resolution determined using normalized realspace MFs for the temperature and polarization Emode data.
Fraction of simulations with higher values of χ^{2} than the observed ones for the T, Q_{r}, and U_{r} angular profiles, computed from the stacking of hot and cold extrema selected above or below the ν = 0 and ν = 3 thresholds.
Fraction of simulations with higher values of χ^{2} than the observed ones for Δμ_{T}, computed from the stacking of hot and cold spots selected above the ν = 0 and ν = 3 thresholds.
pvalues of χ^{2} computed from the oriented stacking of hot and cold spots selected above the ν = 0 threshold.
pvalues of ΔX_{m} computed from the oriented stacking of hot and cold spots selected above the ν = 0 threshold.
Values for the statistic (in units of μK^{4}) for the Planck 2018 data with resolution parameter N_{side} = 64.
Probabilities of obtaining values for the statistic for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2018 CMB maps with resolution parameter N_{side} = 64, estimated using the Commander, NILC, SEVEM, and SMICA maps.
Global pvalue for the S^{XY}(θ, 180° ) statistic for the Planck fiducial ΛCDM model at most as large as the observed values of the statistic for the Planck 2018 CMB maps with resolution parameter N_{side} = 64.
Values for the S^{TQ}(48° ,120° ) statistic (in μK^{4}) and the probability of obtaining values for the Planck fiducial ΛCDM model at least as large as the observed values for the Planck 2018 CMB maps with resolution parameter N_{side} = 64.
Probabilities of getting the statistic for the Planck fiducial ΛCDM model at least as large as the observed values of the statistic for the Planck 2018 CMB maps with resolution parameter N_{side} = 64.
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 for the Commander, NILC, SEVEM, and SMICA temperature and polarization Q and U maps estimated on northern and southern ecliptic hemispheres.
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 statistics for the Commander, NILC, SEVEM, and SMICA temperature and polarization Q and U maps estimated on negative and positive hemispheres defined by the DM reference frame.
Probabilities of obtaining values for the ratio of χ^{2} of the Npoint functions for the Planck fiducial ΛCDM model at least as large as the observed values of the statistics for the Commander, NILC, SEVEM, and SMICA temperature and polarization Q and U maps estimated on northern and southern ecliptic hemispheres.
Probabilities of obtaining values for the ratio of χ^{2} of the Npoint functions for the Planck fiducial ΛCDM model at least as large as the observed values of the statistics for the Commander, NILC, SEVEM, and SMICA temperature and polarization Q and U maps estimated on negative and positive hemispheres defined by the DM reference frame.
pvalues for variance asymmetry measured using different discs for the fullresolution (N_{side} = 2048) Planck 2018 Commander, NILC, SEVEM, and SMICA temperature solutions.
Localvariance dipole directions for the variance asymmetry of the four componentseparated temperature maps.
Localvariance dipole directions and pvalues for the variance asymmetry of the four componentseparated Emode polarization maps at the HEALPix resolution of N_{side} = 64.
Amplitude and direction of the lowℓ dipoleasymmetry signal determined from the QML analysis for the range ℓ = 2 − 64.
Amplitude and direction of the lowℓ dipole modulation signal determined from the QML analysis for the range ℓ = 2–220.
All Figures
Fig. 1. Componentseparated CMB maps at 80′ resolution. Columns show temperature T, and E and Bmode maps, respectively, while rows show results derived with different componentseparation methods. The temperature maps are inpainted within the common mask, but are otherwise identical to those described in Planck Collaboration IV (2020). The E and Bmode maps are derived from the Stokes Q and U maps following the method described in Appendix A. The dark lines indicate the corresponding common masks used for analysis of the maps at this resolution. Monopoles and dipoles have been subtracted from the temperature maps, with parameters fitted to the data after applying the common mask. 

In the text 
Fig. 2. Pairwise differences between maps from the four CMB componentseparation pipelines, smoothed to 80′ resolution. Columns show temperature, T, and E and Bmode maps, respectively, while rows show results for different pipeline combinations. The grey regions correspond to the appropriate common masks. Monopoles and dipoles have been subtracted from the temperature difference maps, with parameters fitted to the data after applying the common mask. 

In the text 
Fig. 3. Examples of common masks. From top to bottom, the masks correspond to those used for analysing temperature maps, polarization represented by the Stokes Q and U parameters, and Emode polarization data, at a resolution N_{side} = 128. From left to right, fullmission, HM, and OE masks are shown. Note that the masks for E and Bmode analysis are extended relative to those derived for Q and U studies, in order to reduce the reconstruction residuals. 

In the text 
Fig. 4. Lowertail probabilities of the variance (top), skewness (centre), and kurtosis (bottom), determined from the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) componentseparated temperature maps (left) and the SEVEM070 (light green), SEVEM100 (dark blue), SEVEM143 (yellow), and SEVEM217 (magenta) frequencycleaned maps (right) at different resolutions. 

In the text 
Fig. 5. Signaltonoise ratio for the variance estimator in polarization for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), obtained by comparing the theoretical variance from the Planck FFP10 fiducial model with an MC noise estimate (righthand term of Eq. (7)). Note that the same colour scheme for distinguishing the four componentseparation maps is used throughout this paper. 

In the text 
Fig. 6. Lowertail probabilities of the variance determined from the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) componentseparated polarization maps (left) and the SEVEM070 (light green), SEVEM100 (dark blue), SEVEM143 (yellow), and SEVEM217 (magenta) frequencycleaned maps (right) at different resolutions. The top, middle and bottom rows correspond to results evaluated with the fullmission, HM and OEcrossvariance estimates, respectively. In this figure, small pvalues would correspond to anomalously low variance. 

In the text 
Fig. 7. Lowertail probabilities of the crossvariance determined between the SEVEM frequencycleaned polarization maps (top) or from the HM (centre) or OE (bottom) SEVEM frequencycleaned maps at different resolutions. In this figure, small pvalues would correspond to anomalously low variance. 

In the text 
Fig. 8. 2point correlation functions determined from the N_{side} = 64 Planck CMB 2018 temperature and polarization maps. Results are shown for the Commander, NILC, SEVEM, and SMICA maps (first, second, third, and fourth rows, respectively). The solid lines correspond to the data, while the black three dotsdashed lines indicate the mean determined from the corresponding FFP10 simulations, and the shaded dark and light grey areas indicate the corresponding 68% and 95% confidence regions, respectively. 

In the text 
Fig. 9. 3point correlation functions determined from the N_{side} = 64 Planck CMB Commander 2018 temperature and polarization maps. Results are shown for the pseudocollapsed 3point (upper panel) and equilateral 3point (lower panel) functions. The red solid line corresponds to the data, while the black three dotsdashed line indicates the mean determined from the FFP10 Commander simulations, and the shaded dark and light grey regions indicate the corresponding 68% and 95% confidence areas, respectively. See Sect. 5.2 for the definition of the separation angle θ. 

In the text 
Fig. 10. Realspace normalized MFs determined from the Planck 2018 temperature data using the four componentseparated maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The grey region corresponds to the 99th percentile area, estimated from the FFP10 simulations processed by the SMICA method, while the dashed curves with matching colours outline the same interval for the other componentseparation methods. Results are shown for analyses at N_{side} = 32, 256, and 1024. The rightmost column shows the χ^{2} obtained by combining the three MFs in real space with an appropriate covariance matrix derived from FFP10 simulations. The vertical lines correspond to values from the Planck data. 

In the text 
Fig. 11. As in Fig. 10, but for the Planck 2018 Emode polarization data. The Planck data are consistent with the Gaussian FFP10 simulations, but variations between the different componentseparation methods are evident. 

In the text 
Fig. 12. Needletspace normalized MFs of the Planck 2018 temperature data using the four componentseparated maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). The grey region corresponds to the 99th percentile area, estimated from the FFP10 simulations processed by the SMICA method, while the dashed curves with matching colours outline the same interval for the other componentseparation methods. The needlet MFs are denoted by the central multipole of the needlet filter ℓ_{c} = 2^{j}; the jth needlet parameter has compact support over the multipole range [2^{j − 1}, 2^{j + 1}]. The rightmost column shows the χ^{2} obtained by combining the three MFs in needlet space with an appropriate covariance matrix derived from FFP10 simulations. The vertical lines correspond to values from the Planck data. 

In the text 
Fig. 13. As in Fig. 12, but for the Planck 2018 Emode polarization data. The Planck data are consistent with the Gaussian FFP10 simulations, but some variations between the different componentseparation methods are evident. 

In the text 
Fig. 14. Cumulative densityfunction of the peak distributions for the SMICA temperature T (left) and reconstructed Emode polarization (right) maps. Top row: peak CDF filtered with a GAUSS kernel of 120′ FWHM, bottom row: peak CDF filtered with the same kernel of 600′ FWHM. The spectral shape parameter γ (see Eq. (29)) is the bestfit value for the simulated ensemble, as indicated by the cyan circle in Fig. 15. No significant deviations from Gaussian expectations are observed. Similar results are obtained for other componentseparation methods. 

In the text 
Fig. 15. Distribution of bestfit Gaussian peak CDF spectral shape parameters, σ and γ (as defined in Eq. (29)), recovered from FFP10 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) for the SMICA temperature T (left) and the reconstructed Emode polarization (right) maps. 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. Upper panel: peak CDF parameters for the SMICA map filtered with a GAUSS kernel of 120′ FWHM, lower panel: corresponding peak CDF using the same kernel with 600′ FWHM. Similar results are obtained for the other componentseparation methods. 

In the text 
Fig. 16. Cumulative densityfunction of the number of all extrema, maxima (red) and minima (blue), derived from simulations, compared to the equivalent values observed for the SMICA temperature T (left) and the reconstructed Emode polarization (right). Upper panel: peak counts for maps filtered with a GAUSS kernel of 120′ FWHM. Lower panel: corresponding peak count CDF for the same kernel of 600′ FWHM. Similar results are obtained for the other componentseparation methods. 

In the text 
Fig. 17. Cumulative densityfunction of the mean amplitude of all extrema, maxima (red) and minima (blue), derived from simulations, compared to the equivalent values observed for the SMICA temperature T (left) and the reconstructed Emode polarization (right). Upper panel: peak mean amplitudes for maps filtered with a GAUSS kernel of 120′ FWHM. Lower panel: corresponding peak CDF for the same kernel of 600′ FWHM. Similar results are obtained for the other component separation methods. The amplitude values are shown in arbitrary (dimensionless) units determined by map prewhitening. 

In the text 
Fig. 18. Upper panels: patches of T, Q_{r}, and U_{r} in microKelvin for maxima above ν = 0 from the WMAP V+W data at a HEALPix resolution of N_{side} = 512 convolved with a Gaussian beam of 30′ (first row), and from the SEVEM data at a HEALPix resolution of N_{side} = 1024 convolved with a Gaussian beam of 10′ (second row). Lower panels: the same thing, but only for maxima above ν = 3. The polarization direction is represented as a headless vector. Note that, in the bottom panel, the length of the headless vectors is divided by 3 with respect to the scale used in the upper panel. 

In the text 
Fig. 19. Mean radial profiles of T, Q_{r}, and U_{r} in microkelvin obtained for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) at N_{side} = 1024. The top panel in each plot shows the radial profile compared to the simulation average (grey line), while the lower panel shows the difference of the radial profile and simulation average on smaller scales. Results based on stacks around temperature maxima and minima are shown in the upper and lower rows, respectively. The left three columns present results for peaks selected above the null threshold, while the right three columns show the equivalent results for peak amplitudes above (maxima) or below (minima) 3 times the dispersion of the temperature map. The black dots (connected by dashed lines) show the mean value from simulations and the shaded regions correspond to the ±1σ (68%) and ±2σ (95%) error bars estimated from SEVEM simulations. Note that the “Diff” curves for each componentseparation method are computed using the corresponding set of ensemble averages, although only the ensemble average from SEVEM is shown here. 

In the text 
Fig. 20. Nonoriented stacks of intensity T (left), reconstructed Emode (centre) and Bmode (right) polarization stacked on intensity maxima for the CMB (top row, SMICA map at N_{side} = 1024 and 10′ FWHM beam), lowfrequency foregrounds (middle row, LFI 30GHz map smoothed using a 10′ FWHM beam, with subtraction of the SMICA CMB map at 10′ FWHM resolution smoothed by the 30GHz LFI beam), and highfrequency foregrounds (bottom row, 353GHz HFI polarizationsensitive bolometer map smoothed using a 10′ FWHM beam, with subtraction of the SMICA CMB map at 10′ FWHM resolution smoothed by the 353GHz HFI beam). The common temperature mask was used for temperature peak selection, and the polarization confidence mask of Appendix A was used for E and Bmode stacks. 

In the text 
Fig. 21. Oriented stacks of intensity T (left), reconstructed Emode (centre) and Bmode (right) polarization stacked on intensity maxima for the CMB (top row, SMICA map at N_{side} = 1024 and 10′ FWHM beam), lowfrequency foregrounds (middle row, LFI 30GHz map smoothed using a 10′ FWHM beam, with subtraction of the SMICA CMB map at 10′ FWHM resolution smoothed by the 30GHz LFI beam), and highfrequency foregrounds (bottom row, 353GHz HFI polarizationsensitive bolometer map smoothed using a 10′ FWHM beam, with subtraction of the SMICA CMB map at 10′ FWHM resolution smoothed by the 353GHz HFI beam). The orientation of the stacked patch is chosen so that the principal directions of the tensor ∇_{i}∇_{j}∇^{−2}T are aligned with the horizontal and vertical axes of the stack. The common temperature mask was used for temperature peak selection, and the polarization confidence mask of Appendix A was used for the E and Bmode stacks. 

In the text 
Fig. 22. Radial profiles of oriented stacks of intensity T (left two columns) and reconstructed E mode (right two columns) for Commander, NILC, SEVEM, and SMICA pipelines. The top panel in each plot shows the radial profile compared to the simulation average (grey line), while the lower panel shows the difference of the radial profile and simulation average on smaller scales. The shaded grey regions represent 68%, 95%, and minimumtomaximum bounds for individual realizations. 

In the text 
Fig. 23. S^{XY}(θ, 180° ) statistic for the N_{side} = 64 Planck CMB 2018 temperature and polarization maps. Results are shown for the Commander, NILC, SEVEM, and SMICA maps (first, second, third, and fourth rows, respectively). The solid lines correspond to the data, while the black three dotsdashed lines indicate the median determined from the corresponding FFP10 simulations, and the shaded dark and light grey areas indicate the corresponding 68% and 95% confidence areas, respectively. 

In the text 
Fig. 24. Difference of the Npoint correlation functions determined from the N_{side} = 64 Planck CMB Commander 2018 temperature and polarization maps and the corresponding means estimated from FFP10 MC simulations. Results are shown for the 2point (upper panels), pseudocollapsed 3point (middle panels), and equilateral 3point (lower panels) functions. The blue and red dashed lines correspond to the functions computed on the negative and positive hemispheres, respectively, determined in the dipolemodulation coordinate frame. The blue and red solid lines correspond to the functions computed on the northern and southern ecliptic hemispheres, respectively. The shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively, determined for the negative DM hemisphere. See Sect. 5.2 for the definition of the separation angle θ. 

In the text 
Fig. 25. Upper panels: ratio R^{TT}(ℓ_{max}) of the Planck 2018 data 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, respectively. Left panel: ratio computed from the componentseparated maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (cyan), after application of the common mask. Note that these results substantially overlap each other. Right panel: ratio computed when the LklCommander map (magenta) is analysed with the corresponding mask. Lower panels: lowertail probability estimator as a function of ℓ_{max} (left), and ℓ_{min} with a fixed ℓ_{max} = 24 (right). For these lower panels small probabilities would correspond to anomalously odd parity. 

In the text 
Fig. 26. Upper panels: D^{TE}(ℓ_{max}) for the SMICA HM (left panel) and OE (right panel) data. The shaded grey regions indicate the distribution of the statistic derived from the corresponding MC simulations as in Fig. 25. Lower panels: lowertail probability of the polarization estimators as a function of ℓ_{max} for the HM data (left panel) and OE data (right panel) for all componentseparated methods. Results for the SEVEM 143GHz cleaned maps are also shown (magenta line). For these lower panels small probabilities would correspond to anomalously odd parity. 

In the text 
Fig. 27. Upper panels: D^{EE}(ℓ_{max}) for the Commander HM (left panel) and OE (right panel) data. The shaded grey regions indicate the distribution of the statistic derived from the corresponding MC simulations as in Fig. 25. Lower panels: lowertail probability of the polarization estimators as a function of ℓ_{max}, as described in Fig. 26. For these lower panels small probabilities would correspond to anomalously odd parity. 

In the text 
Fig. 28. Kolmogorov–Smirnov deviation of the peak distribution for the SMICA temperature T (left) and the reconstructed Emode polarization (right) maps in 70° radius discs centred on the positive and negative asymmetry directions as determined by the SMICATT, TE, EE QML estimator in Sect. 7.2. Top panels: maps filtered with GAUSS kernels of 120′ and bottom panels: filtering with 600′ FWHM. 

In the text 
Fig. 29. Locations of the largescale maxima or minima, numbered 1 to 5, considered in the analysis and indicated on an inpainted version of the Commander map at N_{side} = 2048. In particular, Peak 5 corresponds to the Cold Spot. 

In the text 
Fig. 30. Temperature (left column) and Q_{r} (right column) profiles of the largescale Peaks 1 (top panel) through 5 (bottom panel), as shown in Fig. 29. The black solid lines represent the profiles obtained from the SEVEM CMB map. The red curves correspond to the expected theoretical profiles, while the blue lines represent the theoretical profiles rescaled by the estimated amplitude A for each profile. The shaded regions represent the ±1, ±2 and ±3σ confidence levels. 

In the text 
Fig. 31. Upper panel: pvalues for variance asymmetry measured as the fraction of simulations with localvariance dipole amplitudes larger than those inferred from the data, for lowresolution (N_{side} = 64) temperature maps. The pvalues are given as a function of disc radius and for the four componentseparated temperature maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). Lower panel: corresponding localvariance dipole directions for the four componentseparation temperature maps, and for 4° discs with the lowest pvalues. Note that the four directions match almost perfectly, so that the symbols essentially overlap. For reference, we also show the CMB dipole direction, the north ecliptic pole (NEP), the south ecliptic pole (SEP), and the preferred dipolar modulation axis (labelled as “lowℓ”) derived from the temperature data in Sect. 7.2. 

In the text 
Fig. 32. Impact of bias correction and weighting on the uniformity of dipole directions for Emode polarization HMHD simulation maps. Upper panel: localvariance dipole directions for the Commander simulations of Emode polarization maps, obtained for 4° discs, when no bias correction and weighting have been applied to the maps. Middle panel: same as in the upper panel, but when the transformation X_{i} → (X_{i} − M_{i})/σ_{i} has been applied to the polarization maps before applying the localvariance estimator. Lower panel: same as in the upper and middle panels, but when the transformation has been applied to the maps (i.e., with a square in the denominator). The results for the other componentseparation methods are very similar to the ones presented here. For reference, we also show the CMB dipole direction, the north ecliptic pole (NEP), the south ecliptic pole (SEP), and the preferred dipolar modulation axis (labelled as “lowℓ”) derived from the temperature data in Sect. 7.2. 

In the text 
Fig. 33. 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 Emode polarization maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), at the HEALPix resolution of N_{side} = 64. Lower panel: corresponding localvariance dipole directions for the four componentseparation Emode polarization maps, and for 4° discs. For reference, we also show the CMB dipole direction, the north ecliptic pole (NEP), the south ecliptic pole (SEP), and the preferred dipolar modulation axis (labelled as “lowℓ”) derived from the temperature data in Sect. 7.2. 

In the text 
Fig. 34. Probability determined from the QML analysis for a Monte Carlo simulation to have a larger dipolemodulation amplitude than the Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) data sets, with ℓ_{min} = 2. The topleft, topright, and bottomleft panels use TT, EE, and TE, while the final panel is for all data combined. Small pvalues here correspond to anomalously large dipolemodulation amplitudes. We emphasize that the statistic here is cumulative and apparent trends in the curves can be misleading. 

In the text 
Fig. 35. Distribution of directions for power dipoles estimated from independent 100multipole bins for simulated HM Commander maps with the common mask applied. The fraction of dipoles pointing towards a given ecliptic latitude direction is shown, where the data are binned by the cosine of the direction, since it is this quantity that is uniformly distributed on the sky for a Gaussian random field. The horizontal black lines show the expected 95% deviation taken from a uniform distribution. The coloured lines show the distributions for temperature at ℓ < 800 (red line), temperature at ℓ > 800 (green line), EE polarization at ℓ < 800 (blue line), and EE polarization at ℓ > 800 (magenta line). The temperature results are determined using uniform weighting, while the EE results have inversevariance weighting applied. 

In the text 
Fig. 36. Left panels: dipole directions for independent 100multipole bins of the local powerspectrum distribution from ℓ = 2 to 1500 in the Commander OE maps, with the OE common mask applied. The preferred directions for maps of specific multipole bins are colour coded according to the central value of the bin, ℓ_{central}, as shown in the colour bar. Note that the maps have been rotated about the Galactic northsouth axis, so that Galactic longitude l = 180° is in the centre of the map. The top panel shows the directions for the TT power spectrum, the middle panel for the TE power spectrum and the bottom panel for the EE power spectrum. In all panels, we also show the CMB dipole direction, the north ecliptic pole (NEP), the south ecliptic pole (SEP), and the preferred dipolar modulation axis (labelled as “lowℓ”) derived from the temperature data in Sect. 7.2. The average directions determined from the two multipole ranges ℓ = 2 − 300 and ℓ = 750 − 1500 are shown as blue and dark red (open) rings, respectively. The mean error on the derived directions that results from masking the data is 44° in the range ℓ < 800 and 62° in the range ℓ > 800 for temperature, but 78° in the range ℓ < 800 and 91° in the range ℓ > 800 for EE polarization. Right panels: derived pvalues for the angular clustering of the power distribution in OE maps as a function of ℓ_{max}, determined for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), based on FFP10 simulations using the OE common mask. These are the same colours used throughout the paper (e.g., see Fig. 5), while the grey line shows the Commander results when excluding the first 100 multipoles in the analysis. These 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, hence small pvalues would correspond to anomalously aligned dipole directions. The results shown here have been marginalized over bin sizes in the range Δℓ = 8 to Δℓ = 32. 

In the text 
Fig. 37. As for Fig. 36, but for the HM maps, with the HM common mask applied. 

In the text 
Fig. 38. Modulation results at high ℓ. The shaded bands show the 68% spread of pvalues for asymmetry taken from Commander simulations with a mean pvalue less than 10% for ℓ > 1000. The dark green band (dark green boundaries) represents results computed from the OE data set using the fullmission common mask, the grey band (with grey boundaries) to the same data set analysed with the OE common mask applied, and the light green band (with light green boundaries) corresponds to the HM data set analysed using the HM common mask. The red solid lines show the same three cases for the data: lower line, fullmission common mask; middle line, OE common mask; and upper line, HM common mask. 

In the text 
Fig. 39. Left panels: as for Fig. 36, but for the Commander OE maps with the fullmission common mask applied. The mean error on the derived direction that results from masking the data is 39° in the range ℓ < 800 and 50° in the range ℓ > 800 for temperature, but 78° in the range ℓ < 800 and 91° in the range ℓ > 800 for EE polarization. Right panels: derived pvalues for the angular clustering of the power distribution in OE maps as a function of ℓ_{max}, determined for Commander (red line) based on simulations with the fullmission common mask applied. The grey line shows the Commander results when excluding the first 100 multipoles in the analysis, the black solid line shows results for the Commander HM split, and the black dashed line corresponds to the grey line for the Commander HM split, in all cases with the common mask applied. 

In the text 
Fig. 40. Derived pvalues for the angular correlation of TT and EE dipole directions in Commander maps as a function of ℓ_{max}, determined using the HM split with the fullmission common mask (black), the HM split with the HM common mask (green), the OE map with the fullmission common mask (red), and the OE map with the OE common mask (blue), based on FFP10 simulations. The magenta line shows the corresponding correlation between TE and EE directions for the Commander HM split with the common mask. 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, hence small pvalues would correspond to anomalously aligned dipole directions. The results shown here have been marginalized over bin sizes in the range Δℓ = 8 to Δℓ = 32. 

In the text 
Fig. 41. Dipole directions for independent 10multipole bins of the local powerspectrum distribution from ℓ = 2 to 200 in the Commander HM maps with the fullmission common mask applied. This is a finer binning in order to investigate the directional clustering of different multipole ranges. Note that the maps have been rotated about the Galactic northsouth axis, such that Galactic longitude l = 180° is in the centre of the map. The top, middle, and bottom plots show maps based on the TT, TE, and EE spectra, respectively. In all panels, we also show the CMB dipole direction, the north ecliptic pole (NEP), the south ecliptic pole (SEP), and the preferred dipolar modulation axis (labelled as “lowℓ”) derived from the temperature data in Sect. 7.2. 

In the text 
Fig. A.1. Modecoupling matrices up to ℓ_{max} = 512 for the common polarization mask with point sources at N_{side} = 1024 for masking (top), inpainting (middle), and purifiedinpainting (bottom) methods. The coloured shading represents normalized matrix elements on a logarithmic scale, with values spanning from 10^{−8} (deep blue) to 1 (dark brown). 

In the text 
Fig. A.2. Examples of the finitedifference stencils for HEALPix pixels with eight (left and centre) and seven (right) neighbours for N_{side} = 8 (top) and N_{side} = 1024 (bottom) in ring ordering (see also Bowyer et al. 2011). Red circles represent positions of HEALPix pixel centres in a gnomonic projection onto a plane tangent to the central pixel (i.e., looking straight down at the tangent plane), with the dotted grid aligned with local and directions, illustrating the average pixel pitch h = (π/3)^{1/2}/N_{side}. The white bars represent the directions of polarization, specifically the direction of the polarization ellipse for the +Q polarization mode. Single and doublearrow vectors show projections of and directions, respectively, for neighbouring pixels onto the tangent plane. The black cross corresponds to the average pixel position, while the blue dotted ellipse represents the pixel position covariance. Although for some positions on the sky, the polarization directions are aligned, this is not at all true near the poles (pixel 0); hence just adding Q and U does not make sense. Additionally we can see that because the grid is distorted, secondorder finitedifference schemes need more than just nearest neighbours to work. 

In the text 
Fig. A.3. Multigrid inpainting schedule for an N_{side} = 128 map. Filled nodes represent successive iterations of diffusion steps (Eq. (A.17)), downward strokes represent calculation of the residual (Eq. (A.19)) and its downgrade to the coarser grid, while upward strokes represent the upgrade of the correction (Eq. (A.20)) to the finer grid, and solid arrows represent the merge of the correction into the solution (Eq. (A.21)). Eight diffusion steps are used at all grid levels except the lowest, where 64 steps are used, which is enough to find the static solution there. 

In the text 
Fig. A.4. E and Bmode reconstruction residuals, shown in the left and right columns respectively, for the methods of mask multiplication (top), simple inpainting (middle), and purified inpainting (bottom). Residuals are shown for a single realization of a SMICA componentseparated simulation containing both CMB signal and noise at a resolution N_{side} = 1024. The grey area represents the common polarization mask applied during the reconstruction. 

In the text 
Fig. A.5. Rms of the combined reconstruction residuals for E and B modes (left and right columns, respectively) determined from the purified inpainting of the FFP10 simulations for all four componentseparation methods (i.e., the average of the square of each of the four residuals). The grey area represents the common polarization mask, while the semitransparent grey area indicates an expansion of the confidence mask, which together admit 65% of the sky for further analysis. The mask increment is determined from thresholding the total residual δE^{2} + δB^{2}, smoothed by an 80′ FWHM Gaussian. The rms reconstruction accuracy is better than 0.5 μK, with the largest deviations mainly localized near the mask boundary. 

In the text 
Fig. A.6. PseudoC_{ℓ} spectra of residuals for E and B modes (left and right columns, respectively) determined from the purified inpainting of the FFP10 simulations for all four componentseparation methods. Solid curves show averaged spectra, with red corresponding to the residual masked with the reconstruction mask and blue corresponding to the residual masked with the confidence mask. Semitransparent areas filled with grey show 68%, 95%, and minimumtomaximum bounds for individual realizations. 

In the text 
Fig. A.7. Maximal rms of reconstruction residuals for E modes (green), B modes (orange), and combined E and B modes (blue) determined from purifiedinpainting simulations as a function of sky fraction. The red circle and the thin dashed lines show the reconstruction quality within the common polarization mask used in reconstruction (77.7% of the sky admitted, maximum rms combined residual of 9.0 μK). The blue star and the thin solid lines show the reconstruction quality within the confidence mask (64.8% of the sky admitted, with a maximum rms combined residual of 0.5 μK). 

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.