Planck 2018 results. VII. Isotropy and Statistics of the CMB

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 $\Lambda$CDM cosmological model, yet also confirm the presence of several so-called"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 $E$-mode 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 non-Gaussianity 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., $\ell \lesssim 400$). In this regime, no unambiguous detections of cosmological non-Gaussianity, 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 $\Lambda$CDM cosmological model, and also gives a clear indication of how Planck provides state-of-the-art measurements of CMB temperature and polarization on degree scales.


Introduction
This paper, one of a set associated with the 2018 release of data from the Planck 1 mission (Planck Collaboration I 2018), 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 2018), 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 full-sky measurements of the CMB temperature fluctuations of the presence of such "anomalies," and indicating that a modest degree of deviaflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA). 1 arXiv:1906.02552v1 [astro-ph.CO] 6 Jun 2019 tion 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 "look-elsewhere 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 cosmic-variance 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 full-sky CMB polarization signal have the potential to provide an improvement in significance of the detection of large-scale 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 Sachs-Wolfe (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 full-mission Planck full-sky 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 high-pass filtering of the component-separated 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 2018), means that such a procedure is no longer necessary. Nevertheless, our studies remain limited both by the relatively low signal-to-noise 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 non-Gaussianity, as might arise from non-standard inflationary models, are provided in a companion paper (Planck Collaboration IX 2018). 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 full-mission 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 non-Gaussianity. 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 B-mode maps from the Stokes Q and U data.

Polarization analysis preamble
Traditionally, the Stokes parameters Q and U are used to describe CMB polarization anisotropies (e.g., Zaldarriaga & Seljak 1997a). However, unlike intensity, Q and U are not scalar quantities, but rather components of the rank-2 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 (noise-biased) 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 cross-estimator 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 pointn ref so that Q r (n;n ref ) = −Q (n) cos (2φ) − U (n) sin (2φ), U r (n;n ref ) = Q (n) sin (2φ) − U (n) cos (2φ), 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 positionn. 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 vectorsn, transform as spin-2 variables under rotations around then axis, they can be expanded as where ±2 Y m (n) denotes the spin-weighted spherical harmonics and a (±2) m 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 component-separation approaches. Nevertheless, the presence of residual foregrounds mandates the use of a mask, the application of which during the generation of E-and B-mode maps results in E/B mixing (Lewis et al. 2002;. In Appendix A, we describe the method adopted in this paper to reduce such mixing.

Data description
In this paper, we use data from the Planck 2018 full-mission 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, thermally-unstable phase of the mission. The release includes sky maps at nine frequencies in temperature, and seven in polarization, provided in HEALPix format , 3 with a pixel size defined by the N side parameter. 4 For polarization studies, the 353-GHz maps are based on polarization-sensitive bolometer (PSB) observations only (see Planck Collaboration III 2018, for details).
Estimates of the instrumental noise contribution and limits on time-varying systematic artefacts can be inferred from maps that are generated by splitting the full-mission data sets in various ways. For LFI, half-ring 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, odd-ring (O) and even-ring (E) maps are constructed using alternate pointing periods, i.e., either odd or even numbered rings, to avoid the correlations observed previously in the half-ring data sets. However, for convenience and consistency, we will refer to both of these ring-based splits as "odd-even" (OE), in part as recognition of the signal-to-noise ratios of the LFI and HFI maps and their relative contributions to the component-separated maps described below. Half-mission (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 full-mission data set in the case of HFI. Note that important information on the level of noise and systematic-effect residuals can be inferred from maps constructed from half-differences of the half-mission (HMHD) and odd-even (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 2018;Planck Collaboration III 2018).
As in previous studies, we base our main results on estimates of the CMB from four component-separation algorithms, -Commander, NILC, SEVEM, and SMICA -as described in Planck Collaboration IV (2018). 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. Foreground-cleaned versions of the 70-, 100-, 143-, and 217-GHz sky maps generated by the SEVEM algorithm, hereafter referred to as SEVEM-070, SEVEM-100, SEVEM-143, and SEVEM-217, 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 sky-coverage masks.
The CMB temperature maps are derived using all channels, from 30 to 857 GHz, and provided at a common angular resolution of 5 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 "odd-even" data split that combines the LFI half-ring 1 with the HFI odd-ring data, and the LFI half-ring 2 with the HFI even-ring data.
Lower-resolution versions of these data sets are also used in the analyses presented in this paper. The downgrading procedure is as follows. The full-sky 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,  Fig. 1. Component-separated CMB maps at 80 resolution. Columns show temperature T , and E-and B-mode maps, respectively, while rows show results derived with different component-separation methods. The temperature maps are inpainted within the common mask, but are otherwise identical to those described in Planck Collaboration IV (2018). The E-and B-mode 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.
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 component-separated Q and U data. Figure 1 presents the E-and B-mode maps for the four component-separation methods at a resolution of N side = 128, with the corresponding common masks overplotted. Planck Collaboration IV (2018) notes that some broad large-scale features aligned with the Planck scanning strategy are observed in the Q and U data. The detailed impact on E-and B-mode 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 large-scale residuals in the pairwise differences of the component-separated maps. Finally, we note that the B-mode polarization is strongly noise dominated on all scales, therefore, although shown here for complete-ness, we do not present a comprehensive statistical analysis of these maps.
In general, we make use of standardized masks made available for temperature and polarization analysis, as described in detail in Planck Collaboration IV (2018). These masks are then downgraded for lower-resolution 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 5 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

Fig. 2.
Pairwise differences between maps from the four CMB component-separation pipelines, smoothed to 80 resolution. Columns show temperature, T , and E-and B-mode 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.
pre-processing is applied to them before the application of the component-separation algorithms. Specifically, the value of any unobserved pixel is replaced by the value and 2018 data sets, due to a change in the condition number threshold at the map-making stage.
of the corresponding N side = 64 parent pixel. Analysis of the component-separated 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 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 E-mode polarization data, at a resolution N side = 128. From left to right, full-mission, HM, and OE masks are shown. Note that the masks for E-and B-mode analysis are extended relative to those derived for Q and U studies, in order to reduce the reconstruction residuals.
the component-separation 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 full-mission standardized masks.
The masks for E-and B-mode 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 CMB-plus-noise realizations, as propagated through all four Planck component-separation pipelines, then evaluating the residuals with respect to the input fullsky maps. The final mask is specified by requiring that the maximal rms level of the residuals observed in the simulations is less than 0.5 µK, significantly below the cosmological E-mode signal.
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 full-mission, 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.

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 data-reduction pipeline over the full-mission period. These are extensions of the "full focal-plane" 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 (2018) 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 frequency-dependent. The signal realizations in- Table 1. Standardized data sets used in this paper. The resolutions of the sky maps used are defined in terms of the N side parameter and corresponding FWHM of the Gaussian beam with which they are convolved. The fraction of unmasked pixels in the corresponding common masks for the full-mission (Full), as well as the HM and OE data splits, are also specified. clude the frequency-specific beam properties of the LFI and HFI data sets implemented by the FEBeCoP (Mitra et al. 2011) beam-convolution approach. Given that the instrumental noise properties of the Planck data are complex, we make use of a set of socalled "end-to-end" 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 map-making 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 (2018). A similar approach is followed by LFI to generate noise MCs that capture important characteristics of the scanning strategy, detector response, and data-reduction pipeline over the fullmission period (Planck Collaboration II 2018). A total of 300 realizations are generated at each Planck frequency, for the full-mission, 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 full-focal plane simulations.
Finally, the CMB signal and noise simulations are propagated through the various component-separation pipelines using the same weights as derived from the Planck fullmission data analysis (Planck Collaboration IV 2018). The signal and noise realizations are then permuted to generate 999 simulations 7 for each component-separation 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 analysis. 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 2018) 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 signal-tonoise 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 full-mission data set. Note that detailed comparisons have been undertaken using the power spectra of the individual frequency maps. Figure 18 of Planck Collaboration II (2018) compares half-ring halfdifference (HRHD) spectra for the LFI 30-, 44-and 70-GHz data with simulations, finding good agreement over most angular scales. Figure 17 of Planck Collaboration III (2018) makes a similar comparison for the HFI 100-, 143-, 217-and 353-GHz half-mission HMHD and OEHD data.
Of more importance to this paper, however, is the consistency of the data and simulations after various component-separation methods have been applied. As established in Planck Collaboration IV (2018), the corresponding end-to-end 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 2018). 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 point-source residual contribution not included in FFP10. It should, therefore, be apparent that systematic shifts over some ranges of angular scale could contribute to p-value 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 component-separation 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.

Tests of non-Gaussianity
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 2018). Searching for departures from this scenario is crucial for its validation, yet there is no unique signature of non-Gaussianity. 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 non-exhaustive 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 signal-to-noise remains relatively low. The analyses are therefore applied to all four component-separation products (Commander, NILC, SEVEM, and SMICA) at a given resolution with the accompanying common mask, and significance levels are determined by comparison with the corresponding results derived from the FFP10 simulations. The consistency of the results derived from the various component-separation 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 non-Gaussianity depending on the statistical test applied, given the significant level of the systematic effects modelled therein.

One-dimensional 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 component-separated maps using the unit-variance estimator . This method requires a normalized variance sky map, u X defined as: where X i is the observed temperature at pixel i, σ 2 X,0 is the variance of the CMB signal, and σ 2 i,N is the variance of the noise for that pixel, estimated using the FFP10 MC simulations. The CMB variance is then determined by finding thê σ 2 X,0 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 lower-tail probability of the variance, skewness, and kurtosis determined at different resolutions from the four component-separated maps (left columns) and from the SEVEM frequency-cleaned maps (right columns), after applying the appropriate common mask. There is good agreement between the maps, although the NILC results indicate a slightly lower p-value 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 figure 15 of Planck Collaboration IV (2018). We note that Planck showing a decreasing lower-tail probability with decreasing resolution. This lower-tail 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 (Planck Collaboration XVI 2016;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. In polarization we follow a different approach, as a consequence of the lower signal-to-noise ratio. Specifically, we subtract the noise contribution to the total variance of the polarization maps and define the estimator Signal-to-noise 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 (right-hand term of Eq. (7)). Note that the same colour scheme for distinguishing the four component-separation maps is used throughout this paper.σ where Q and U are the Stokes parameters of the observed polarization maps, and Q 2 N + U 2 N MC are noise estimates determined from MC simulations. Planck Collaboration IV (2018) 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 signal-dominated polarization results, we also define a cross-variance 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 OE-cleaned data split, and Q N 1 Q N 2 + U N 1 U N 2 MC is the corresponding noise contribution to the total variance in polarization estimated from the corresponding simulations. Note that a cross-estimator 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 signal-to-noise ratio of the polarization variance for the component-separated 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, Q 2 N + U 2 N MC . All of the methods show similar behaviour, with a maximum signal-to-noise 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 signal-to-noise ratio increases again, but with an amplitude that depends noticeably on the component-separation method considered . At very high resolutions the signal-to-noise ratio drops, as expected.
In Fig. 6 we show the lower-tail probabilities of the variance determined from the full-mission and the HM and OE data splits using the appropriate common mask, compared to the corresponding results from MC simulations. At high resolutions, the lower-tail probability determined from the variance of the full-mission 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 lower-tail probability becomes more compatible with the MC simulations when we consider the cross-variance 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 tem-  perature data. This behaviour is common to all of the component-separated methods and also to the SEVEM-143 frequency-cleaned data, although with different probabilities at a given resolution. The SEVEM-070 frequency-cleaned 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 component-separation methods and resolutions, even when a cross-estimator is considered, suggests that correlated noise and residual systematics in the data may not be sufficiently well described by the current set of simulations.
In an attempt to minimize the impact of correlated noise, we consider the cross-variance estimated between pairs of frequencies from the SEVEM frequency-cleaned maps, using full-mission, HM and OE data sets. The results are shown in Fig. 7. Note that the combination of the SEVEM-070 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 70-GHz cleaned data. In addition, the cross-variance determined between the SEVEM-070 and SEVEM-217 maps yields a particular low probability. Since the 217-GHz data are used to clean the 70-GHz 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 cross-variance data.
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 p-values found, depending on the specific combinations of data considered. This indicates that the lower signal-to-noise 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.

N-point correlation functions
In this section, we present tests of the non-Gaussianity of the Planck 2018 temperature and polarization CMB data using real-space N -point correlation functions.
An N -point correlation function is defined as the average product of N observables, measured in a fixed relative orientation on the sky, where the unit vectorsn 1 , . . . ,n N span an N -point polygon. If statistical isotropy is assumed, these functions do not depend on the specific position or orientation of the N -point 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 directionn. 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 coordinatesystem-independent N -point 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, n ref , is specified by the centre of mass of the polygon (Gjerløw et al. 2010). In the case of the 2-point 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 w i 1 , · · · , w i N 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 2-point function, as well as the pseudo-collapsed and equilateral configurations for the 3-point function. Following Eriksen et al. (2005), the pseudo-collapsed 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 is the difference between the observed,Ĉ N (θ i ), and the corresponding average from the MC simulation ensemble, C N (θ i ) , of the N -point 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 ∆ (k) N (θ i ) is the kth simulated N -point 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 N sim = N sim −1. However, due to degeneracies in the covariance matrix resulting from an overdetermined system and a precision in estimation of the matrix elements of order ∆M ij ∼ 2/N sim , the inversion of the matrix is unstable.
To avoid this, a singular-value decomposition (SVD) of the matrix is performed, and only those modes that have singular values larger than 2/N sim 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 (N sim − N bin − 1)/N sim (Hartlap et al. 2007).
We analyse the CMB estimates at a resolution of N side = 64 due to computational limitations. The results for the 2-point correlation functions of the CMB maps are presented in Fig. 8, while in Fig. 9 the 3-point 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 3-point functions derived from the simulations indicates the presence of small non-Gaussian contributions, presumably associated with modelled systematic effects that are included in the simulations. Furthermore, both the mean and associated confidence regions vary between component-separation 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 (2018). 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.
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 N -point 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 goodness-of-fit between the data and the simulations than a simple inspection of the figures can reveal.
The N -point function results show excellent consistency between the CMB temperature maps estimated using the different component-separation 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 pseudo-inverse 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 signal-to-noise ratio of the polarized data on large angular scales. Interestingly, a tendency towards very high probability values is observed for the pseudo-collapsed T T Q r 3-point functions for all methods, and for the equilateral T Q r Q r functions in the case of Commander and SEVEM.
As an alternative to the Stokes parameters, we also consider N -point functions computed from the temperature and E-mode 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 T T E 3-point functions for all component-separation methods.
Nevertheless, we conclude that no strong evidence is found for statistically significant deviations from Gaussianity of the CMB temperature and polarization maps using N -point correlation functions.
Finally, we note that the results for the T T correlation function confirm the lack of structure at large separation angles, noted in the WMAP first-year 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 T Q r correlation function.

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 non-Gaussianity 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 E-mode representation, since the MFs computed from the spin-2 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 MF-based 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 1983;Pogosyan et al. 2009;Gay et al. 2012;Matsubara 2010;Fantaye et al. 2015). This factorization is valid in the weakly non-Gaussian 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 configu-  ration. 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 scale-dependent 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 non-Gaussian features detected with the MFs.
First, we undertake a real-space analysis by computing the three normalized functionals described above at different resolutions and smoothing scales for each of the four component-separation 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 simulations: whereȳ sim ≡ y sim 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 E-mode 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 variance-weighted difference with respect to the mean of the simulations for the three MFs. The right-most 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 p-values determined for the different component-separation 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 E-mode.
For the temperature results, the χ 2 values computed for the different component-separation methods are more consistent than was the case for the Planck 2015 analysis, for all scales. In the case of the E-mode results, we find no significant discrepancy between the Planck data and the FFP10 simulations. The striking variation in the p-values for the four component-separation methods, is also observed when considering individual realizations in the set of simulations.
As a complement to the pixel-based 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 pixel-space 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 high-frequency; and the doublelocalization properties of needlets (in real and harmonic space) allow a much more precise, scale-by-scale, interpretation of any possible anomalies. While the behaviour of standard all-scale (pixel-based) 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 all-scale MFs. This is because if one changes the shape of the power spectrum while still keeping [2 + 1]C constant, the pixel-space 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 non-Gaussianity detected by the MF analysis.
The needlet components of a scalar CMB field are defined by Marinucci et al. (2008) and Baldi et al. (2009) are given by where j on the left-hand side is the needlet index and j on the right-hand side is a power. Here, X (n) denotes the component at multipole of the CMB map X(n) (corresponding to temperature or the polarization E or B modes), i.e., wheren ∈ S 2 denotes the pointing direction, B 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 ( /B 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 Euler-Poincaré characteristic, boundary length, and area of the full sphere, respectively. The quantities v k are the normalized MFs given in Eq. (17), while the needlet-scale amplitudes A j k 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 B = 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 needlet-scale 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 pixel-space analysis without alteration.
Once the needlet maps are obtained, we follow the identical procedure as in the pixel-space case to compute the three MFs. The results derived from the Commander, NILC, SEVEM, and SMICA component-separated temperature maps for needlet scales B = 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 pixel-space MFs case, the results indicate a high degree of consistency among the four component-separation methods.
The results for the E-mode polarization data are shown in Fig. 13. Compared to the temperature results, the different component-separation 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.
In summary, both the pixel-and needlet-space MF analyses show that the 2018 Planck temperature and polarization maps are consistent with the Gaussian simulations over angular scales corresponding to max = 2048.

Peak statistics
In this section, we present non-Gaussianity 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 component-separated 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 sphericalharmonics-based methods. Peak statistics are particularly sensitive to non-Gaussian 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 and Hou et al. (2009), and Planck data in Planck Collaboration XVI (2016).
We consider peak statistics from the Planck componentseparated temperature and polarization maps at N side = 1024, with E-mode maps reconstructed by the purified inpainting method described in Appendix A. The maps are pre-whitened by convolving them with an isotropic function derived from the isotropic best-fit CMB power spectrum, combined with a diagonal approximation to the instrumen- tal noise covariance. Then, a confidence mask is applied, and weighted convolution is performed with a 2D-Gaussian smoothing kernel (that we label as "GAUSS"), as described in Appendix A of Planck Collaboration XVI (2016). 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 density-function (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 CDFF (x) 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 p-values we derive its CDF directly from the simulations. The peak distributions for T and E-mode 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 CDFF (x) 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.
A further statistical test of isotropic Gaussian random field expectations is to check if the best-fit parameters, σ and γ, to the observed empirical peak CDF, agree with those derived from individual simulated realizations. The distribution of best-fit 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.
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 To summarize, the temperature and E-mode peak statistics determined from the Planck component-separated maps show no significant anomalies, except perhaps for the previously known Cold Spot discussed in Section 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.

Non-oriented 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 Marcos-Caballero 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,n ref , is specified by the centre of each extremumn 0 . 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 foreground-cleaned map at a HEALPix pixel resolution of N side = 1024 convolved with a Gaussian beam of 10 , and a noise-weighted 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 signal-to-noise 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 Cumulative density-function of the peak distributions for the SMICA temperature T (left) and reconstructed E-mode polarization (right) maps. The top row shows the peak CDF filtered with a GAUSS kernel of 120 FWHM, the bottom row shows the peak CDF filtered with the same kernel of 600 FWHM. The spectral shape parameter γ (see Eq. 29) is the best-fit 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 component-separation methods.
Planck data allows additional inner rings to be observed compared to WMAP.
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 component-separation pipelines, compared to the simulations. The corresponding p-values for the comparison are presented in Table 7. We see that all component-separation methods yield consistent results. No significant differences for the intensity profiles µ T (θ) are observed with respect to the results found in the Planck-2015 analysis. As in the latter case, a systematic deviation between the data and the mean value of simulations is present. This was previously interpreted (PCIS15) as an effect connected with the deficit in the observed power spectrum at low multipoles.
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 p-values. We therefore consider Distribution of best-fit 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 E-mode polarization (right) maps. The blue contours enclose 68 % and 95 % of the parameter distribution, and the cyan circle represents the best-fit parameters for the median peak CDF determined from simulations. The upper panel shows the peak CDF parameters for the SMICA map filtered with a GAUSS kernel of 120 FWHM, while the lower panel shows the corresponding peak CDF using the same kernel with 600 FWHM. Similar results are obtained for the other component-separation methods.
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 p-values obtained in this case are given in Table 8, and are consistent with the deviations shown in Fig. 19. Turning to polarization, while Table 7 generally shows consistent results between the data and simulations, somewhat low p-values 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 p-values 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 component-separation methods. In addition, it was shown in Planck Collaboration Int. XLIX (2016) Table 7. 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. In this table, small p-values would be considered anomalous. 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 non-oriented stacking are consistent with the predictions of the ΛCDM model.

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 Planck Collaboration XVI (2016), 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 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 cross-correlation power spectra; these can be used to characterize non-Gaussian 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 signal-to-noise 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 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 spin-2 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 flat-sky approximation, the temperature-derived 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 equal-area azimuthal Lambert projection to represent the stacks, which introduces the radial variable which is almost identical to the usual angular variable θ in the flat-sky 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 non-trivial angular modes:  Fig. 19. Mean radial profiles of T , Q r , and U r in micro-kelvin obtained for Commander (red), NILC (orange), SEVEM (green), and SMICA (blue) at N side = 1024. Each individual panel contains the mean radial profiles (top part) and the differences between the mean profiles of the data and those computed from the ensemble mean of the simulations (denoted "Diff," lower part). 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 component-separation method are computed using the corresponding set of ensemble averages, although only the ensemble average from SEVEM is shown here.
where δ m0 is the Kronecker delta and X stands for T or E. Statistical isotropy of the stacked field X would imply that ensemble-averaged radial functions for odd ms vanish. Similar to non-oriented 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 profileX m ( ).
To assess foreground contamination in componentseparated CMB temperature and polarization data, we present unoriented and oriented stacks of temperature T and full-sky reconstructions of E and B modes for Planck CMB maps compared to low-and high-frequency 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 component-separated CMB maps at N side = 1024, with a 10 FWHM beam. The middle rows show stacks of the LFI 30-GHz 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 beam-convolved with the LFI 30-effective beam, thus representing the total lowfrequency foreground emission. The bottom rows show the HFI 353-GHz polarization-sensitive bolometer maps convolved with a 10 FWHM Gaussian beam and degraded to N side = 1024, then CMB-corrected using the SMICA map further beam-convolved with the HFI effective beam, thereby representing the total high-frequency emission. The polarization maps are transformed to E and B modes via the full-sky spherical harmonic transform, while the Q T and U T signals used to define patch orientations are derived from full-sky 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 • .
The CMB data show a characteristic ringing pattern in the temperature stacks, associated with the first acoustic peak, a high signal-to-noise pattern in the E-mode stack (as expected from the CMB T E correlation), and no evidence of detectable patterns in the B-mode stack (reaching , reconstructed E-mode (centre) and B-mode (right) polarization stacked on intensity maxima for the CMB (top row, SMICA map at N side =1024 and 10 FWHM beam), low-frequency foregrounds (middle row, LFI 30-GHz map smoothed using a 10 FWHM beam, with subtraction of the SMICA CMB map at 10 FWHM resolution smoothed by the 30-GHz LFI beam), and high-frequency foregrounds (bottom row, 353-GHz HFI polarization-sensitive bolometer map smoothed using a 10 FWHM beam, with subtraction of the SMICA CMB map at 10 FWHM resolution smoothed by the 353-GHz 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 B-mode stacks.
a noise floor of about 0.05 µK), which serves as a null test for foreground contamination (given the fact that both lowand high-frequency foregrounds display T B correlations, as described below To investigate the quantitative agreement of the Planck component-separated maps with fiducial cosmological and noise models, we compare the radial profiles (Eq. 39) of oriented temperature and polarization stacks for all component-separation 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 , reconstructed E-mode (centre) and B-mode (right) polarization stacked on intensity maxima for the CMB (top row, SMICA map at N side =1024 and 10 FWHM beam), low-frequency foregrounds (middle row, LFI 30-GHz map smoothed using a 10 FWHM beam, with subtraction of the SMICA CMB map at 10 FWHM resolution smoothed by the 30-GHz LFI beam), and high-frequency foregrounds (bottom row, 353-GHz HFI polarization-sensitive bolometer map smoothed using a 10 FWHM beam, with subtraction of the SMICA CMB map at 10 FWHM resolution smoothed by the 353-GHz HFI beam). The orientation of the stacked patch is chosen so that 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 B-mode stacks.   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 B-mode 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 spin-2 orientation operator.
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 minimum-to-maximum 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 upper-tailed p-values of the χ 2 deviation (Eq. 32) of oriented radial profiles, while Table 10 summarizes the two-tailed p-values of the integrated integrated profile deviations (Eq. 40). None of the p-values appear anomalous.
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.

28
Planck Collaboration: Isotropy and statistics of the CMB

Anomalies in the microwave sky
The previous section established the lack of evidence for significant non-Gaussianity 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 large-angle 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 odd-parity 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 E-mode 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.

Lack of large-angle correlations
We assess the lack of correlation in the 2-point 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. 2013a;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. (2013b) 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Ĉ XY 2 (θ) is an estimate of the corresponding 2-point correlation function. This is a generalization of the S 1/2 statistic (Spergel et al. 2003) computed from the temperature auto-correlation 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. (2013a) (although note that the sense of the p-values differs between the papers). The p-values 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. Copi et al. (2013b) demonstrated that the signal-tonoise 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 σ S 1/2 < ∼ 5 × 10 −5 µK 4 for the S T T 1/2 statistic, σ S 1/2 < ∼ 7 × 10 −6 µK 4 for the remaining statistics with estimated values at least of order 10 −5 µK 4 , and σ S 1/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 best-fit 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 S XY 1/2 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 component-separation methods. No equivalent excursions are seen for any statistics involving polarization data, 10 We fit the quadrupole to the masked maps using a modified version of the HEALPix routine remove_dipole.  which also show notable variations between componentseparation methods. We provide a more quantitative assessement 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 look-elsewhere effect, we define a global statistic to evaluate the true signifi- cance of the result. Specifically, we repeat the procedure for each simulation, and search for the largest probability irre- spective 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 p-value. As seen in Table 13, this corresponds to values of order 99 % for all of the component-separated temperature maps, and much smaller global p-values 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 (PCIS15) (i.e., 99 % compared to 98 %). More importantly, the p-values for the temperature maps after removing the best-fit quadrupole drop to around 86 %, indicating again that the anomalous value of the S T T statistic is connected to the observed quadrupole. As pointed out by Copi et al. (2013b), the crosscorrelation between temperature and polarization can be used to determine whether the S T T value is low due to a statistical fluke in the context of a ΛCDM cosmology. Specifically, if this hypothesis were true, the S T Q statistic would also likely be small on large angular scales, given that the polarization signal is partially correlated with temperature. Copi et al. (2013b) have determined that the S T Q 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 seven-year 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 T Q (48 • , 120 • ) statistic measured with the Planck data are significantly smaller, so that we cannot rule out the statistical fluke hypothesis.  Muir et al. (2018) have noted that the low value of the large-angle correlation statistic S T T 1/2 can be connected to the temperature 2-point correlation function for antipodean points, C T T 2 (180 • ). This is motivated by fact that, as seen in Fig. 8, the otherwise largely flat C T T 2 (θ) drops to negative values close to θ = 180 • . Here, we extend the analysis of the temperature correlation function to the functions including polarization information, C XY 2 (180 • ). Table 15 presents the probabilities of finding larger values of the 2-point 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 p-values differs between the papers). When the best-fit 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.

Hemispherical asymmetry
In this section, we reassess the asymmetry between the real-space N -point 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 2-point and 3-point 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 N -point 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 non-Gaussianity (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 dipole-modulation (DM) analysis in Sect. 7.2. In the latter case, the positive and negative hemispheres are defined as those for which the 31 Planck Collaboration: Isotropy and statistics of the CMB 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 component-separated maps in Table 24, for the data combination T T, T E, EE. We do not include an additional analysis of the hemispheric split defined by the Doppler-boost 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 N -point functions for the data and the mean values estimated from the FFP10 MC simulations.
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 asymmetry between the two measured hemispheres can be quantified by the ratio of the corresponding χ 2 values. The probabilities of obtaining values of the χ 2 statistic 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 2-tailed statistic.
Looking at the results in Tables 16-19, we no longer observe the high significance level for the pseudo-collapsed 3-point function of the temperature map in the northern ecliptic hemisphere, as reported for 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 T T Q r equilateral 3-point 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 p-values for the SMICA T U r 2-point function and the NILC Q r U r U r pseudo-collapsed 3-point 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 correlation-function configurations considered will increase the likelihood of finding some apparently anomalous behaviour, even for statistically isotropic data.
In the DM reference frame, the 3-point 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 component-separated maps for the Q r Q r Q r , T T U r and T Q r Q r pseudo-collapsed 3-point functions, and for a few configurations of the equilateral 3-point func-     tion. Of course, the inconsistency in p-values for different component-separated maps argues against results of cosmological significance, but rather indicates the presence of different levels of foreground or systematic-effect residuals, depending on the component-separation method.

Point-parity 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 parity-symmetric T + (n) and parity-antisymmetric T − (n) functions defined as where for a given directionn, the antipodal point is −n.
Under the parity transformation T (P ) (n) = P [T (n)] = T (−n), it can be shown that the corresponding a m coefficients transform such that a (P ) m = (−1) a m . As a conse- quence, T + (n) and T − (n) 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 B-mode angular power spectra with only even or odd -modes, since under the parity transformation a B, m = (−1) +1 a B, m . On the largest angular scales, 2 < < 30, corresponding to the Sachs-Wolfe plateau of the temperature power spectrum, we expect the angular power spectra of T + (n) and T − (n) to be of comparable amplitude. However, an odd point-parity preference has been observed both in various WMAP data releases (Land & Magueijo 2005b;Kim & Naselsky 2010;Gruppuso et al. 2011), and confirmed in the Planck 2013 and 2015 studies. Furthermore, Land & Magueijo (2005a) have shown that this cannot be easily ascribed to the presence of residual Galactic foreground emission because of its even point-parity 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 C TT + ( max ) and C TT − ( max ) are given by with +,− tot being the total number of even (+) or odd (−) multipoles included in the sum up to max , and C TT , is the temperature angular power spectrum computed using a quadratic maximum-likelihood (QML) estimator (Gruppuso et al. 2011;Molinari et al. 2014).
The top-left 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 large-scale parity. We additionally consider an additional estimate of the CMB sky determined using the Commander component-separation methodology that has been optimized for large angular scale analyses. This map is also used in the Planck lowlikelihood analysis (Planck Collaboration V 2018), thus we refer to it as Lkl-Commander. This corresponding result is shown in the top right panel of the figure. The results from the component-separation products are in good agreement, and indicate an odd-parity preference for the multiple range considered in this test. The lower-left panel of Fig. 25 shows the lower-tail 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 lower-tail 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 lower-tail 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.
The Lkl-Commander map shows good agreement with the standard component-separated 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 lower-tail 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 Lkl-Commander analysis, we count the number of MC simulations with a lower-tail 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 look-elsewhere effect, an odd-parity preference is observed with a lower-tail probability of about 1.6 %.
In the lower-right panel of Fig. 25, the lower-tail probability as a function of min is presented when max is constrained to be 24, corresponding to the lowest probability found for the Lkl-Commander 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 point-parity analysis to include polarization data, specifically considering the T E 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 signal-tonoise 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 T E or EE spectra, and C X +,− ( max ) 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 2018), we determine the power spectrum using the cross-quasi-QML estimator that forms the basis of the HFI low-likelihood in polarization (Planck Collaboration V 2018). Here, the signal covariance matrix is constructed according to the FFP10 fiducial model, while the noise is specified using the FFP8 11 143-GHz noise-covariance matrix. This is clearly a suboptimal description of the noise present in the component-separated 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 noise-covariance matrix. To verify this, we also analyse the SEVEM 143-GHz 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 cross-estimator to both the HM and OE data sets.
Examination of the EE and T E power spectra determined from the component-separated 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 2018) 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 point-parity analysis results in probabilities that reach values of 98.5 % for certain component-separated 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 lower-tail probability as a function of max for the HM and OE data. Results from the component-separated 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. 25. Lower panels: The lower-tail 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.
maps. As in the T E analysis, the data and MC simulations are consistent, although the p-values decrease over the max range of 20 to 30, reaching minima of about 6 % at max = 27 for the HM data; however, the OE results do not indicate such a trend.
It is interesting to note that the C XY 2 (180 • ) 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., C T Qr 2 (180 • ) = C T Ur 2 (180 • ) = 0. The factor (−1) splits the sum into contributions from even (parity-symmetric) and odd (parity-antisymmetric) multipoles with opposite signs. The functions are therefore differences between parity-even and parity-odd multipole band powers, 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 lower-tail probability is found. However, the low signal-to-noise 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.

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 SMICA T T T EEE QML estimator (see Table 24). Separate peak CDFs are then constructed for these directions, then compared to the full-sky 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 E-mode 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 The Cold Spot was first detected in the WMAP temperature data as an anomalously cold region in terms of the spherical Mexican-hat 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;Planck Collaboration XVI 2016). Recently, Marcos-Caballero et al. (2017a) analysed the local properties of the Cold Spot (and other large-scale 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) andFernández-Cobos 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 large-scale 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 Marcos-Caballero et al. 2017b,a, 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 to 4 in Fig. 29, all located in the southern Galactic hemisphere. Given their locations, they are expected to make a significant contribution to the hemispheric-power asymmetry. In addition, the peaks should induce a particular pattern in the polarization field, due to the correlation between the temperature and E-mode polarization.
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 T B 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 Marcos-Caballero et al. 2016, for technical details).
In order to analyse the large-scale 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 Marcos-Caballero et al. (2016) and generate 2000 CMB simulations that are constrained to have a peak located at the same position as every large-scale 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 large-scale 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 Marcos-Caballero 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 (Marcos-Caballero 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 (Marcos-Caballero 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: HereX(θ i ) 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 large-scale 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 Fernández-Cobos 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 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 Marcos-Caballero 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 signal-to-noise 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.

Dipole modulation and directionality
In this section, we examine dipolar violations of isotropy. In Planck Collaboration XVI (2016), 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 dipolar-modulation models to the data. Aluri & Shafieloo (2017) applied a local-variance estimator to low-resolution E-mode 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 pixel-based method to maps of the squared polarization amplitude, but found no evidence for the presence of a dipolar-modulation 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 pixel-to-pixel 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 band-power 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.

Variance asymmetry
We first study the Planck 2018 temperature and polarization maps using the local-variance method introduced by Akrami et al. (2014). Despite its relative simplicity, the local-variance estimator serves as a powerful method for detecting violations of statistical isotropy if they are of a type corresponding to a large-scale 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 E-mode 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 large-scale power asymmetry in CMB temperature data, and allow us to focus on the question of whether a corresponding large-scale 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 local-variance 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 local-variance maps. Finally, we fit a dipole to each of the local-variance 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 local-variance 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 local-variance maps here. However, although we do not consider higher-multipole asymmetries, Table 20. Amplitudes A estimated from Eq. (51) of the T and Q r profiles for the large-scale peaks considered here. The Cold Spot corresponds to Peak 5.

Amplitudes Peak
Comm such features may exist in local-variance maps estimated from the polarization data. Table 21 presents the results for a local-variance analysis of the full-resolution, N side = 2048, Planck 2018 temperature data over a range of disc radii, 4 • ≤ r disc ≤ 20 • . The p-values are measured as the fractional number of simulations with local-variance dipole amplitudes larger than those inferred from the data. The results are consistent for all the component-separated 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 local-variance 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 p-values as the disc radius is increased to values larger than 8 • .
Since our focus is on large angular scales, we repeat the local-variance analysis for low-resolution, N side = 64, temperature maps. Figure 31 (upper panel) shows the pvalues for the component-separated 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 component-separation methods. The overall level of significance is lower compared to the full-resolution results, but high significance (low p-value) results are found for small disc radii. The lowest p-values correspond to the smallest disc radius, i.e., 4 • discs, and the p-values 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 local-variance dipole directions for the four component-separation methods are provided in Table 22 for both the high-resolution and lowresolution cases.
The analysis of the polarization data is significantly more subtle than for temperature because of the inherently low signal-to-noise. We first validate the technique by applying it to polarization simulations. The analysis shows that the direct application of the method to isotropic simulations of E-mode polarization signal returns local-variance Table 21. p-values for variance asymmetry measured using different discs for the full-resolution (N side = 2048) Planck 2018 Commander, NILC, SEVEM, and SMICA temperature solutions. The values represent the fraction of simulations with local-variance dipole amplitudes larger than those inferred from the data, and hence small p-values correspond to anomalously large variance asymmetry.  Table 22. Local-variance dipole directions for the variance asymmetry of the four component-separated temperature maps. All directions quoted here are for 4 • discs, and for full-resolution (N side = 2048), as well as downgraded, lowresolution (N side = 64) maps. dipole directions that are not uniformly distributed. This arises from the strongly anisotropic, correlated, and non-Gaussian 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 signal-to-noise ratio is large.
The following procedure is adopted to minimize the effects of the anisotropic noise. We first preprocess the E- Fig. 31. Upper panel: p-values for variance asymmetry measured as the fraction of simulations with local-variance dipole amplitudes larger than those inferred from the data, for low-resolution (N side = 64) temperature maps. The pvalues are given as a function of disc radius and for the four component-separated temperature maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue). Lower panel: Corresponding local-variance dipole directions for the four component-separation temperature maps, and for 4 • discs with the lowest p-values. 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. mode 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 local-variance dipole directions for an analysis of the HMHD E-mode simulation maps when unmodified, and after transformation according to Eq. (53) with p = 1 or 2. Since we are interested in large-scale 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 non-uniformity 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. Figure 33 (upper panel) presents the p-values obtained by applying the local-variance estimator to the four component-separated E-mode polarization maps. Even though we typically see an increase in p-values 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 component-separated 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.
Nevertheless, the preferred directions we have found for the E-mode 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 E-mode and temperature-variance dipole directions for the data as well as for all the simulations, and computed p-values 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 local-variance dipole directions and corresponding p-values for the four component-separation methods, together with the separation angles between temperature and E-mode dipole directions and their associated p-values.
Despite the bias correction and weighting procedure employed for reducing the non-uniformity 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.

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 dipole-modulated skies and skies consistent with ΛCDM (although this statement is somewhat model dependent). Furthermore, the Fig. 32. Impact of bias correction and weighting on the uniformity of dipole directions for E-mode polarization HMHD simulation maps. Upper panel: local-variance dipole directions for the Commander simulations of E-mode polarization maps, obtained for 4 • discs, when no bias correction and weighting have been applied to the maps. Middle panel: The 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 local-variance estimator. Lower panel: The same as in the upper and middle panels, but when the transformation X i → (X i − M i )/σ 2 i has been applied to the maps (i.e., with a square in the denominator). The results for the other component-separation 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.
analysis we carry out here is a purely phenomenological one performed in multipole space, with no attempt to connect to any real-space modulation; however, several physical kspace models are considered in a companion paper Planck Collaboration X (2018).

Fig. 33.
Upper panel: p-values for variance asymmetry measured as the number of simulations with local-variance dipole amplitudes larger than those inferred from the data, as a function of disc radius for the four componentseparated E-mode polarization maps, Commander (red), NILC (orange), SEVEM (green), and SMICA (blue), at the HEALPix resolution of N side = 64. Lower panel: corresponding local-variance dipole directions for the four componentseparation E-mode 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.
We employ a version of the estimator proposed in Planck Collaboration XVI (2016), which was further developed in  and Contreras et al. (2017). In particular, we usẽ with Here the δC +1 are determined by the model of modulation; for this case we have assumed scale-invariant power modulation and thus δC +1 = 2A s . The remaining terms are:X M , the spherical harmonic transform of the modulation amplitude and direction; W Z = T T , T E, or EE; W m and Z m , which are inverse-covariance filtered data; Table 23. Local-variance dipole directions and p-values for the variance asymmetry of the four component-separated E-mode polarization maps at the HEALPix resolution of N side = 64. The p-values represent the fraction of simulations with local-variance dipole amplitudes larger than those inferred from the data. All directions quoted here are for 4 • discs. We also show the separation angles and degrees of alignment between the preferred directions inferred from temperature, T , and E-mode polarization data, with the p-values measured as the fraction of simulations with separation angles smaller than those inferred from the data.
Commander . . . F W ≈ W m W * m ; and the last term on the right-handside of Eq. (56) denotes the mean-field correction (details, including the precise form of F W 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 B m = ( + m + 1)( + m + 2) 2(2 + 1)(2 + 3) .
TheX M 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 E-mode polarization, these coupling coefficients neglect corrections on the very largest scales ( 10) due to the non-local definition of E-modes. This gives rise to a slightly different coupling induced in E as compared to T by a dipole modulation of the primordial fluctuations ). 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 inverse-variance weighting over all data combinations (T T , T E, EE) to obtain a combined minimum-variance 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 scale-invariant 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 . For example, modulated fluctuations that exhibit a 7 % dipolar asymmetry in T to ≈ 65 are not expected to produce an E-mode asymmetry of the same amplitude and over the same scales. Therefore, strictly speaking, the "model" being tested when we combine T T , T E, and EE has no physical basis; a more physical, k-space approach is considered in the companion paper Planck Collaboration X (2018).
In Fig. 34 we show the p-values of the full-mission data compared to the FFP10 simulations, considering T T (top left), EE (top right), or T E (bottom left) data independently. While the T T results are largely consistent with our previous analysis (differences are within the expected scatter given the different analysis masks, PCIS15), our EE and T E results are new. The polarization-only 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 cross-correlation results are rather featureless in the full range of max considered. In the max ≈ 250 region NILC shows systematically lower p-values compared to the other methods, although the modulation amplitudes are still statistically consistent. This difference could be attributable to a percent level p-value 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.
In the bottom right panel of Fig. 34 we combine temperature and polarization (including T E) data and show the p-values of the data compared to simulations. Only the p-value 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 k-space, as previously mentioned (also see , and references therein). Table 24. Amplitude and direction of the low-dipoleasymmetry signal determined from the QML analysis for the range = 2-64. The errors are calculated from the FFP10 simulations. Recall that the expectation for noisefree Gaussian skies (the "cosmic variance of dipole modulation," see equation 50 in PCIS15) corresponds to an amplitude of approximately 3 %. The polarization data on their own show no evidence of modulation, and the addition of polarization has very little effect on the temperature asymmetry signal. In Table 24 we show the measured amplitude and direction of dipolar modulation in the oft-quoted = 2-64 range, with and without polarization data. Polarization clearly has very little effect on the modulation parameters in this region, and T E in particular has a nearly negligible effect at this scale. It is worth recalling that the expectation for purely Gaussian skies is a dipole-modulation amplitude of approximately 3 % (see equation 50 in PCIS15). In Table 25   Table 25. Amplitude and direction of the low-dipole modulation signal determined from the QML analysis for the range = 2-220. The errors are calculated from the FFP10 simulations. This range corresponds to the lowest p-value for the T T , T E, EE data. we show the amplitude and direction of dipolar modulation for the combined T T , T E, and EE data in the = 2-220 range, which is the range with the lowest p-value. 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 lowanomaly on the larger scales where the dipole-modulation amplitude is larger. Removing angular scales < 100 eliminates the significance of the modulation entirely.
In PCIS15 we applied a look-elsewhere 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.

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 sub-sections, in PCIS15, we confirmed the apparent presence of a deviation from statistical isotropy in the Planck data using an angular-clustering 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 inde- p-values here correspond to anomalously large dipole-modulation amplitudes. We emphasize that the statistic here is cumulative and apparent trends in the curves can be misleading. pendent 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. 2018). Nevertheless, regarding this as a purely empirical question, it is important to repeat the directional-clustering analysis and broaden its scope to include polarization. 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.
modes due to incomplete sky coverage is removed in the power-spectrum estimation during the inversion of the full T EB 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. 16 2. For each power-spectrum multipole bin, an N side = 1 HEALPix map with the local power distribution is constructed. 3. The best-fit dipole amplitude and direction is estimated from this map using inverse-variance weighting, where the variance is determined from the local spectra computed from the simulations. As in previous papers (PCIS15 and PCIS13), 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 p-values determined as follows. We first construct the RS using all multipoles up to max . The p-value is then given by the fraction of simulations with a higher RS than for the data for this max . A small p-value therefore means that there are few simulations that exhibit 16 See footnote 15.
as strong clustering as the data. Note that the p-values are highly correlated because the RS as a function of max is cumulative. 6. After calculating results for each bin size, we calculate the variance-weighted mean of the power spectra over all bin sizes (the C for a given bin size is weighted by 1/ √ N b 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 100-multipole 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 inverse-variance weighting in the latter case. Deviation from a uniform distribution is seen for the polarization directions, most notably for the small-scale 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.
The three panels on the left of Fig. 36 show the T T , T E, and EE dipole directions determined for fifteen successive 100-multipole 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 north-south 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 p-values 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 to 1000, the p-values are below 1 % for all values of the maximum multipole for all methods. However, from max ≈ 1000, the p-values increase rapidly. Figure 37 presents the equivalent results for the HM data set.
The application of the OE-or HM-unobserved 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 per-pixel inverse-noise 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 p-values 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 p-values for max > 1000 and check if a similar trend can be seen. As a trade-off between ensuring a low p-value and having a sufficient number of simulations to obtain reliable statistics, we select those simulations with a mean p-value p < 10 % for > 1000. In Fig. 38 we show the 68 % spread of these p-values obtained from OE Commander simulations with the full-mission 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 p-values increase for smaller scales as a result of the higher uncertainty on direction caused by the unobserved pixel masks.
In Fig. 39 we show the results for the Commander OE analysis when the full-mission 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.
We note that for max < 100 the temperature p-values 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 large-angular-scale features. Note also that the dip in p-values for the EE polarization at ≈ 200 is still seen even with the restriction min > 100.
It is also apparent from Fig. 39 that some p-values are close to 100 %, in particular for T E. A high p-value 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 p-value of the RS statistic for the T E 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 T T 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 T T . In this way, the statistic measures the cross-correlation between T T and EE directions. In Fig. 40 we show the p-values for the cross-correlation statistic as a function of multipole. Note that a strong correlation between the position of the T T and EE dipoles are detected for the same multipole range where the EE polarization dipoles (and T T dipoles) are clustered.
We have seen that the EE directions appear clustered for max ≈ 200. This is also the case for T T , 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 T T and EE dipoles clustered among themselves, but that they also appear to be clustered towards the same direction. If T T and EE were independent, this would be highly unexpected, even if T T and EE were separately clustered. However, we know that the T E spectrum is non-zero, giving rise to a correlation between T and E. In order to test whether such a T T -EE directional correlation is expected in the case when both T T and EE are clustered individually, we perform the following test. We examine all simulations having a minimum p-value of less than 1 % for both T T and EE in overlapping multipole ranges (similar to what we observed in the data). Only two simulations have overlapping T T and EE multipole ranges using this criterion, but neither of them has a significant correlation between T T 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 T T (and the T E 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 T T , T E, and EE. There is some hint here of the correlation between T T 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 T E directions. Clearly the T E directions are more scattered than T T and EE, as expected from Fig. 39. The directions for several multipole blocks in the range < 100 coincide quite well with the T T and EE directions, but most of the 10-multipole blocks point in different directions. To further compare the angular clustering of T E dipoles with T T and EE dipoles, we construct the statistic measuring the correlation between T E and T T directions, as well as between T E and EE directions, in the same way as described above for the correlation between T T and EE. We find no significant correlations between the T T and T E directions, but the correlation between T E and EE again shows some sign of similar multipole ranges where EE is aligned and the T T and EE directions correlate. This is shown as a magenta line in Fig. 40. The correlation between T E and EE directions at ≈ 200 may arise as a result of the fact that T T and EE are aligned, as well as the fact that T T and EE directions are also correlated at exactly these multipoles. However, since none of the simulations show a high Note that the maps have been rotated about the Galactic north-south axis, so that Galactic longitude l = 180 • is in the centre of the map. The top panel shows the directions for the T T power spectrum, the middle panel for the T E 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 p-values 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 p-values 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 p-values would correspond to anomalously aligned dipole directions. The results shown here have been marginalized over bin sizes in the range ∆ = 8 to ∆ = 32. 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 T T and T E if this was really the case.
In PCIS15, we made a more detailed study of the total significance of the alignment based on the p-values in figures like Fig. 39. Specifically, we compared the mean and minimum p-values (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 sim-ulations have lower mean p-values than the data. This is valid for most masks or foreground-subtraction methods. Because the data have a minimum p-value of 0 for all masks and methods, only an upper bound on significance can be set. Considering only min > 100, the p-values are still low.
Turning now to polarization, for T E, the numbers are completely consistent with simulations, while for EE, the results are unstable to the choice of masks and componentseparation methods. The mean p-value is always consistent with simulations, which is not unexpected given the short multipole range where the p-values are low (as discussed above). However, for the minimum p-value, 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 signifi- Commander simulations with a mean p-value less than 10 % for > 1000. The dark green band (dark green boundaries) represents results computed from the OE data set using the full-mission 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, full-mission common mask; middle line, OE common mask; and upper line, HM common mask. cance can be set. A similar conclusion can be made for the T T /EE angular correlation, but the results are in this case more stable with choice of mask and component-separation method.

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 large-angular-scale 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 model-independent approach and focus on null-hypothesis testing by calculating and reporting p-values for a number of statistical tests. These tests are performed on maps of the CMB anisotropy that originate from the four component-separation methods, Commander, NILC, SEVEM, and SMICA, described in Planck Collaboration IV (2018). For polarization studies, we consider both maps of the Stokes parameters, Q and U , and of the E-mode signal generated using a novel method described in Appendix A. The consistency of the results determined from the component-separated maps is an important indicator of the cosmological origin of any significant p-values, 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 sig- natures 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, N -point correlation functions, Minkowski functionals, peak statis-  Fig. 36, but for the Commander OE maps with the full-mission 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 p-values 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 full-mission 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. tics, 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 large-scale angular correlations, a hemispherical asymmetry in the behaviour of N -point functions or peak distributions, a violation of point-parity 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 well-known large-scale dipolar power asymmetry. Neither investigations using a variance estimator nor via to ± 1 mode coupling find strong evidence of this asymmetry. However, an interesting align- Fig. 40. Derived p-values for the angular correlation of T T 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 full-mission 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 T E and EE directions for the Commander HM split with the common mask. The p-values 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 p-values would correspond to anomalously aligned dipole directions. The results shown here have been marginalized over bin sizes in the range ∆ = 8 to ∆ = 32. ment 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 component-separated map in question. The mode-coupling estimator indicates that the polarization-only 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 p-value dips would increase in significance when even statistically isotropic polarization data were added. Secondly, all of these dipolemodulation and hemispheric-asymmetry 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-T T dipole-modulation 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 large-scale and lownoise polarization experiments.

Fig. 41.
Dipole directions for independent 10-multipole bins of the local power-spectrum distribution from = 2 to 200 in the Commander HM maps with the full-mission 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 north-south 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 T T , T E, 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.
A notable feature of all of the polarization analyses is the variation in p-values for a given test between the four component-separated maps. This is a consequence of the fact that different component-separation 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 2018), thus allowing more robust measurement of the optical-depth-to-reionization τ , residual systematics, and our ability to simulate them, can limit the kind of statistical tests of non-Gaussianity 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 E-modes at the cosmic-variance 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. and therefore not rotationally invariant. Neverthless, scalar E and pseudo-scalar B maps can be determined from the measured quantities, as described in Sect. 2. Such maps offer definite advantages for any map-based analyses, such as those presented in the main body of this paper. Since they are generated via the application of non-local spherical harmonic transformations, the E and B estimators are non-trivial 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 1997b;, defined by They are related to the standard E and B estimators in harmonic space by a factor of ( + 2)( + 1) ( − 1) = ( + 1) ( + 1) − 2 , (A.2) which corresponds to an application of the bi-Laplacian operator such that E, B = −∇ 2 (∇ 2 + 2) ψ E,B in real space. Unlike E and B, E and B maps can be derived from Q and U by a (local) second-derivative 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 spin-0 nor spin-2 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 B-mode components, although numerically this would be implemented using spherical harmonic transforms rather than matrix multiplications.
Mode mixing is a well-known 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 full-sky quantities. For an isotropic Gaussian random field with an ensemble-average spectrum C , the masked sky mode averages C are given by where the coupling matrices K 1 2 depend on the method used, as illustrated in Fig. A.1. For reconstruction of low-multipoles, maximumlikelihood estimators are widely used (Efstathiou 2004;Bielewicz et al. 2004;de Oliveira-Costa & 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 2011;Bunn & Wandelt 2017). Some of the approaches discussed in the literatire require solving large linear algebra problems (typically via iterative solvers), and could be expensive on high-resolution 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 B-mode maps in the case of incomplete sky coverage in Sects. A.2, A.3, and A.4. The suitability of these approaches for our purposes depends upon the uniformity of the reconstruction, since the method-specific 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 B-modes from the Q and U data after the application of a mask that zeros the problematic pixels. As a consequence, E-and B-mode mixing does result, with mode coupling matrices expressible analytically in terms of Wigner-3j symbols and the power spectrum of the mask W = m |W m | 2 /(2 + 1) defined by a window function Mode-coupling matrices up to max = 512 for the common polarization mask with point sources at N side = 1024 for masking (top), inpainting (middle), and purified-inpainting (bottom) methods. The coloured shading represents normalized matrix elements K 1 2 / K 1 1 K 2 2 on a logarithmic scale, with values spanning from 10 −8 (deep blue) to 1 (dark brown).

3
(2 3 + 1) The resultant mode-coupling 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 E and B maps instead, there is no mode mixing between E and B. Unfortunately, second-order derivative operators enhance the noise power, which results in large artefacts in the reconstruction due to high-to-low-mode-coupling in the masking operator, unless extremely strong apodization is applied to the mask (as described in Smith 2006).  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 localx 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 double-arrow 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, second-order finite-difference schemes need more than just nearest neighbours to work.

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 non-Gaussianity 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 Sect. A.3.2, and we further develop the details of the finite-differencing miltigrid approach in Sect. 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 multi-grid methods of Brandt (1977) adapted to the spherical HEALPix pixelization, as described in Sect. A.3.3.
Mode-coupling 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 mode-coupling matrices are not symmetric and result in excellent suppression of the low-to-high-mode mixing, at the expense of increased high-to-low-mode mixing. Similarly to the case of pure E-and B-mode masking, high-to-low-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. Finite-difference stencils
Finite-difference approximations of differential operators are non-trivial to evaluate on the HEALPix grid, especially for polarization. In this section we discuss first-and secondorder-accurate 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 complexc i for tensor values, and h 2 = π/3N 2 side being the average area of a pixel. We will derive the stencil weights in a flat-sky 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, n ≡ sin θ cos φ, sin θ sin φ, cos θ , θ ≡ ∂n ∂θ = cos θ cos φ, cos θ sin φ, − sin θ , φ ≡ 1 sin θ ∂n ∂φ = − sin φ, cos φ, 0 , (A.7) associated with it, which form a basis on a tangent plane (x,ŷ) ≡ (θ,φ) and a normal directionẑ ≡n. Nearby points can be projected onto a selected tangent plane via gnomonic projection, which maps the radius vector r into ρ = r r ·n −n, ρ ·n ≡ 0, (A.8) and thus introduces local coordinates on a tangent plane, (x, y) = (x · ρ,ŷ · ρ). (A.9) Second-order-accurate 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 . Under rotation of the basis (ê 1 ,ê 2 ) used to define the Stokes parameters Q and U , i.e., 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 (x,ŷ) basis when projecting onto a tangent plane. The appropriate angle of rotation can be computed as tan ψ =x ·φ −ŷ ·θ x ·θ +ŷ ·φ , (A.12) with transformation of Stokes parameters most conveniently implemented as a complex phase rotation Q + iU = e −2iψ (Q + iU ). (A.13) A linear-order shift of the average pixel position in the HEALPix grid breaks the symmetry of the local Taylor expansion, and there is a unique second-order-accurate nearest-neighbour discretization for the Laplace operator, as opposed to a one-parameter family on the rectangular grid. Unfortunately, it turns out to be unconditionally unstable for diffusion-type problems, so we will not discuss it here. Instead, we will use an approximate first-order stencil based on the isotropic weighting c i ≡ c(ρ i /h), with coefficients normalized by L[const] = 0 and L[ρ 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 byc i = e −2iψi c i , 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 first-order forward in time discretization of the flow equation φ(t + δt) − φ(t) = (δt/h 2 )L[φ] updates a pixel according to using the weighted sum of itself and its neighbours, where α = δt/h 2 , and the largest timestep δt that can be taken is determined by Courant-Lewy 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: as detailed in Brandt (1977). Recursive application of coarse-grid 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. 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 "w-stroke," 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.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,p, 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 full-sky map, or "ambiguous" attribution where some mode mixing is allowed: p ≡ a + e + b; M ·p = M · p. (A.22) Note that the polarization map components e and b are not pure in the sense of being orthogonal to the entire Eand B-mode linear spaces on the masked sky, as in ; 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: a = p; e = 0; b = 0.
(A.23) Then the method proceeds by peeling off non-ambiguous E and B modes one by one through the explicit construction of a Krylov subspace (Krylov 1931) generated by the inverse bi-Laplacian, K, acting on the masked maps. The obvious starting point that contains most of the CMB power is the E-mode projection of the masked polarization map: w = EM · a. (A.24) To prevent the constructed Krylov subspace basis from becoming degenerate, we use Lanczos bi-orthogonalization (Lanczos 1950). To do this we keep two recent normalized basis vectors, u and v, initially set to β = w · w 1 2 , u = 0, v = w β , (A.25) and project them out (from the next E-mode generated using the inverse bi-Laplacian operator, K) to obtain a new mode from the previous ones via w = EKM·v, w → w−β u, w → w− w·v v. (A.26) Once the new mode is constructed, it is normalized and the most recent basis set is updated: Lanczos bi-orthogonalization 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 E-mode map: α = a · v ; a → a − α v; e → e + α v. (A.28) 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 B-mode 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 ambiguous-mode map.
Inpainting is a linear operation, and can be represented by a matrix operator F. A suitable method for inpainting high-resolution maps is the multi-grid approach, as discussed above, which has computational costs of O(n log n). The final reconstructed polarization map is generated by inpainting the remaining ambiguous modes as a → FM · a. (A.29) 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 multi-grid approach. The low-to-high-mode-mixing is much improved over the equivalent case with masking alone, while the increase in the high-to-low-mode-mixing is substantially less than if a simple inpainting approach were applied. Purified inpainting therefore represents a good compromise for the E-and B-mode reconstruction of CMB polarization maps.
The computational costs of the purified-inpainting method are dominated by the two round-trip full-resolution spherical transforms required for the initial projection (Eq. A.24) of the E and B modes, and the multi-grid inpainting of the ambiguous modes. The Krylov subspace construction does eventually produce a complete basis for the E-and B-mode subspaces, but is too expensive to run to completion for high-resolution 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 bi-Laplacian 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 B-mode maps for a given method can be evaluated from realistic MC simulations of the CMB signal plus noise, where the true full-sky E * and B * maps are known, and can be directly compared to the reconstructed onesẼ andB. Various metrics of performance can be assigned to the masked residual maps δE ≡ M · (Ẽ − E * ) and δB ≡ M · (B − B * ). 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 purified-inpainting methods. The latter seems to perform on par or better than other direct reconstruction methods published so far, and is competitive with the maximum-likelihood estimators, at substantially lower computational cost. Apodization and inpainting methods do not perform as well because of the E-and B-mode mixing arising from the mask, while the pure projection method of Smith (2006) actually does much worse as a consequence of the high--to-low-aliasing of the noise power. The purified-inpainting method, presented here, offers a balance between these two considerations, and yields an order of magnitude lower residuals for N side = 1024 Planck maps.
Extending the mask beyond that used in the recontruction 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 CMB-plus-noise 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 component-separation pipelines are shown in Fig. A  along with the confidence mask that admits 65 % of the sky for further analysis, as derived from the iso-levels of the total rms reconstruction residual δE 2 + δB 2 smoothed by an 80 FWHM Gaussian beam. Figure A.6 shows the pseudo-C anisotropy power spectra of masked residuals for reconstructed E and B modes for all four component-separation pipelines.
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 E-mode signal, thus allowing sensitive tests  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 component-separation methods (i.e., the average of the square of each of the four residuals). The grey area represents the common polarization mask, while the semi-transparent 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. Pseudo-C 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 component-separation 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. Semi-transparent areas filled with grey show 68 %, 95 %, and minimum-tomaximum bounds for individual realizations.
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 purified-inpainting 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).