© Planck Collaboration 2020

1. Introduction

This paper, one of a set associated with the 2018 release of data from the Planck¹ mission (Planck Collaboration I 2020), describes a compendium of studies undertaken to determine the statistical properties of both the temperature and polarization anisotropies of the cosmic microwave background (CMB).

The ΛCDM model explains the structure of the CMB in detail (Planck Collaboration VI 2020), yet it remains entirely appropriate to look for hints of departures from, or tensions with, the standard cosmological model, by examining the statistical properties of the observed radiation. Indeed, in recent years, tantalizing evidence has emerged from the WMAP and Planck full-sky measurements of the CMB temperature fluctuations of the presence of such “anomalies”, and indicating that a modest degree of deviation from global isotropy exists. Such features appear to exert a statistically mild tension against the mainstream cosmological models that themselves invoke the fundamental assumptions of global statistical isotropy and Gaussianity.

A conservative explanation for the temperature anomalies is that they are simply statistical flukes. This is particularly appealing given the generally modest level of significance claimed, and the role of a posteriori choices (also referred to as the “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 2020), 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 2020). The current paper is organized as follows. Section 2 provides a brief introduction to the study of polarized CMB data. Section 3 summarizes the Planck 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.

2. Polarization analysis preamble

Traditionally, the Stokes parameters Q and U are used to describe CMB polarization anisotropies (e.g., Zaldarriaga & Seljak 1997). However, unlike intensity, Q and U are not scalar quantities, but rather components of the 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,

(1)

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² and applies a correction for noise bias determined from simulations. A cross-estimator based on polarization observations from two maps, P² = Q₁Q₂ + U₁U₂ is also considered.

In addition, a local rotation of the Stokes parameters, resulting in quantities denoted by Q_r and U_r, is employed in Sects. 5.2 and 5.5. In this case, a local frame is defined with respect to a reference point so that

(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 position .

Finally, the rotationally invariant quantities referred to as E and B modes are commonly used for the global analysis of CMB data. Since the quantities Q ± iU, defined relative to the direction vectors , transform as spin-2 variables under rotations around the axis, they can be expanded as

(3)

where denotes the spin-weighted spherical harmonics and are the corresponding harmonic coefficients. If we define

(4)

then the invariant quantities are given by

(5)

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; Bunn et al. 2003). In Appendix A, we describe the method adopted in this paper to reduce such mixing.

3. Data description

In this paper, we use data from the Planck 2018 full-mission data release (“PR3”) that are made available on the Planck Legacy Archive (PLA²). 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 (Górski et al. 2005)³, with a pixel size defined by the N_side parameter⁴. For polarization studies, the 353-GHz maps are based on polarization-sensitive bolometer (PSB) observations only (see Planck Collaboration III 2020, 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 2020; Planck Collaboration III 2020).

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 (2020). These provide data sets determined from combinations of the Planck raw frequency maps with minimal Galactic foreground residuals, although some contributions from unresolved extragalactic sources are present in the temperature solutions. 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 frequency-dependent 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′ full width at half maximum (FWHM) and N_side = 2048. In contrast to the 2013 and 2015 releases, these do not contain a contribution from the second order temperature quadrupole (Kamionkowski & Knox 2003). An additional window function, applied in the harmonic domain, smoothly truncates power in the maps over the range ℓ_min  ≤ ℓ  ≤ ℓ_max, such that the window function is unity at ℓ_max = 3400 and zero at ℓ_max = 4000. The polarization solutions include information from all channels sensitive to polarization, from 30 to 353 GHz, at the same resolution as the temperature results, but only including contributions from harmonic scales up to ℓ_max = 3000. In the context of these CMB maps, we refer to an “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, 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 (2020) 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 completeness, 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 (2020). 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⁵ at a given frequency, some 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 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 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 full-sky maps. The final mask is specified by requiring that the maximal root mean square (rms) level of the residuals observed in the simulations is less than 0.5 μK, significantly below the cosmological E-mode signal. Examples of common masks are shown in Fig. 3.

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.

4. Simulations

The results presented in this paper are derived using Monte Carlo (MC) simulations. These provide both the reference set of sky maps used for the null tests employed here, and form the basis of any debiasing in the analysis of the real data, as required by certain statistical methods. The simulations include Gaussian CMB signals and instrumental noise realizations that capture important characteristics of the Planck scanning strategy, telescope, detector responses, and 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 (2020) is slightly lower, at τ = 0.054 ± 0.007. 1000 realizations of the CMB sky are generated including lensing, Rayleigh scattering, and Doppler boosting⁶ effects, the latter two of which are frequency-dependent. The signal realizations include 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 so-called “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 time-ordered 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 (2020). 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 full-mission period (Planck Collaboration II 2020). 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 full-mission data analysis (Planck Collaboration IV 2020). The signal and noise realizations are then permuted to generate 999 simulations⁷ 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 p-value. This is the probability of obtaining a test statistic at least as extreme as the observed one, under the assumption that the null hypothesis (i.e., primordial Gaussianity and isotropy of the CMB) as represented by the simulations is true. However, this also requires that the simulated reference data set adopts a cosmological model that is sufficiently consistent with that preferred by the data. We have noted above that the τ-value used in the FFP10 simulations is high relative to that preferred by the latest cosmological. As a preliminary assessment, we have considered the predicted variance of the CMB signal for two values of τ, specifically 0.060 (as adopted by the FFP10 simulations), and 0.052 (which is representative of the value determined by an analysis of the HFI data in Planck Collaboration VI 2020) with A_s set to an appropriate value. We find that the latter reduces the polarization variance by approximately 20% at N_side = 16 and 32. It may be necessary to take this effect into account when interpreting the polarization results in what follows.

Similar considerations apply to the simulated noise and residual systematic effects, particularly given the signal-to-noise 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 (2020) compares half-ring half-difference (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 (2020) 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 (2020), 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 2020). Alternatively, biases can arise due to the lack of foreground residuals in the simulations. On small angular scales, the power observed in the temperature data exceeds that of the simulations due to a 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.

5. 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 2020). 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⁸ 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.

5.1. 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 (Cruz et al. 2011). This method requires a normalized variance sky map, u^X defined as:

(6)

where X_i is the observed temperature at pixel i, is the variance of the CMB signal, and is the variance of the noise for that pixel, estimated using the FFP10 MC simulations. The CMB variance is then determined by finding the value for which the variance of the normalized map u^X is unity. The skewness and kurtosis are then subsequently computed from the appropriately normalized map.

In Fig. 4 we show the 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 Fig. 15 of Planck Collaboration IV (2020). We note that Planck Collaboration IV (2020) has demonstrated the presence of a noise mismatch between the observed data and simulations, as traced by the HMHD and OEHD maps. However, this is not relevant for analysis of the temperature data, given its very high signal-to-noise ratio. The results for 1D moments presented here are in very good agreement with the Planck 2015 analysis (PCIS15), 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 (PCIS15; Gruppuso et al. 2013) has shown that the low variance anomaly becomes less significant with increasing sky coverage. Indeed, when we apply the 2016 common mask to the current data set, the probability decreases to 0.7–0.8%, in better agreement with the previous results. The skewness and kurtosis results do not show any anomalous behaviour, in agreement with earlier analyses.

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

(7)

where Q and U are the Stokes parameters of the observed polarization maps, and are noise estimates determined from MC simulations. Planck Collaboration IV (2020) indicates a mismatch between the noise in the data and that in simulations for map resolutions above N_side = 256. This corresponds to a few percent of the theoretical CMB variance up to N_side = 1024, while it is much larger at the highest resolution. Since the noise mismatch is likely to affect the less 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:

(8)

where Q₁, Q₂, U₁, and U₂ are the Stokes parameters of the two maps from either the HM- or OE-cleaned data split, and 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, . 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 temperature 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.

5.2. 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,

(9)

where the unit vectors 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 direction . Following the standard CMB convention, Q and U are defined with respect to the local meridian of the spherical coordinate system of choice. To obtain coordinate-system-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, , 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 θ₁, …, θ_2N − 3 characterizing the shape and size of the polygon,

(10)

Here, pixel weights represent masking and are set to 1 or 0 for included or excluded pixels, respectively.

The shapes of the polygons selected for the analysis are not more optimal for testing Gaussianity than other configurations, but are chosen because of ease of implementation and for comparison of the results with those for the 2013 and 2015 Planck data sets. In particular, we consider the 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 χ² statistic to quantify the agreement between the observed data and simulations. This is defined by

(11)

Here, is the difference between the observed, , 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 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

(12)

where . However, due to degeneracies in the covariance matrix resulting from an overdetermined system and a precision in estimation of the matrix elements of order , the inversion of the matrix is unstable. To avoid this, a singular-value decomposition (SVD) of the matrix is performed, and only those modes that have singular values larger than are used in the computation of the χ² statistic (Gaztañaga & Scoccimarro 2005). We note that this is a modification of the procedure used in previous Planck analyses (PCIS13; PCIS15). Finally, we also correct for bias in the inverse covariance matrix by multiplying it by a factor (Hartlap et al. 2007).

We analyse the CMB estimates at a resolution of N_side = 64 due to computational limitations. The results for the 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 N-point 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 (2020). To avoid biases, it is essential to compare the statistical properties of a given map with the associated simulations. Comparing with simulations without systematic effects could lead to incorrect conclusions.

The probabilities of obtaining values of the χ² 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 χ² 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 TTQ_r 3-point functions for all methods, and for the equilateral TQ_rQ_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 χ² statistic for the Planck fiducial ΛCDM model at least as large as the observed temperature and polarization values are provided in Table 4. Here, we see that the most significant deviations between the data and the simulations occur for the TTE 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 TT 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 TQ_r correlation function.

5.3. Minkowski functionals

In this section, we present a morphological analysis of the Planck 2018 temperature and polarization CMB maps using Minkowski functionals. The Minkowski functionals (hereafter MFs) describe the morphology of fields in any dimension and have long been used to investigate 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, σ₀. The three MFs, namely the area V₀(ν) = A(ν), the perimeter V₁(ν) = C(ν), and the genus V₂(ν) = G(ν), are defined respectively as

(13)

(14)

(15)

where N_ν is the number of pixels with |ΔT|/σ₀ >  |ν|, 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 σ₀. 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 σ₀ in our analysis.

Furthermore, the MFs can be written as a product of a function A_k (k = 0, 1, 2), which depends only on the Gaussian power spectrum, and v_k, which is a function only of the threshold ν (see e.g., Vanmarcke et al. 1983; Pogosyan et al. 2009; Gay et al. 2012; Matsubara 2010; Fantaye et al. 2015). This factorization is valid in the weakly 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 configuration⁹. The analytical expressions are

(16)

(17)

with H_n, the Hermite function,

(18)

The amplitude A_k depends only on the shape of the power spectrum C_ℓ through the parameters σ₀ and σ₁, the rms of the field and its first derivative, respectively:

(19)

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 σ₀ units, providing a total of 36 different statistics y = {v₀, v₁, v₂}. A χ² 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 :

(20)

(21)

where is the mean of the statistics y computed on the simulations, i, j are the threshold indices from the combined MFs, and

(22)

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 χ² 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 component-separation methods. The grey shaded regions in the MFs plot and the histogram in the χ² 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 χ² 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 double-localization 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), and are given by

(23)

(24)

where j on the left-hand side is the needlet index and j on the right-hand side is a power. Here, denotes the component at multipole ℓ of the CMB map (corresponding to temperature or the polarization E or B modes), i.e.,

(25)

where denotes the pointing direction, 𝔹 is a fixed parameter that controls the needlet’s band width (usually taken to be between 1 and 3), and b(.) is a smooth function such that ∑_jb²(ℓ/𝔹^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

(26)

where t₀ = 2, t₁ = 0, and t₂ = 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 have a similar form to A_k, but with the variances of the map and its first derivative given by

(27)

(28)

We adopted the needlet parameters 𝔹 = 2, j = 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 for this analysis. Note that the jth needlet scale has compact support over the multipole range [2^j − 1, 2^j + 1]. For clarity in all the figures, we refer to the different needlet scales by their central multipole ℓ_c = 2^j.

To obtain the needlet maps at different scales, we initially decompose into spherical harmonics the temperature and polarization maps, inpainted using diffusive and purified inpainting (see Appendix A), respectively, at N_side = 1024. In all cases, we set the maximum multipole to ℓ_max = 2048, which is twice the maximum resolution considered for the needlet MF analysis. We then obtain the jth 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 𝔹 = 2, j = 2, 6, 9 are shown in Fig. 12. As can be seen from both the MF and χ² 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 component-separated 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.

5.4. 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 spherical-harmonics-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

(29)

where σ is the rms random field amplitude, and γ is the shape parameter, dependent on the spectrum of the Gaussian random field. Peak locations and amplitudes, and various derived quantities, such as their correlation functions, have previously been used to characterize the WMAP maps in Larson & Wandelt (2004, 2005) and Hou et al. (2009), and Planck data in PCIS15.

We consider peak statistics from the Planck component-separated 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 instrumental 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 PCIS15. The mask is further extended by rejecting pixels with an effective convolution weight that differs from unity by more than 12%, and peaks within it are extracted and analysed. The empirical cumulative density-function (CDF) of peak values x, defined for a set of n peaks at values {X_i},

(30)

is generated by sorting the peak values {X_i} extracted from the map in ascending order, and comparison to the median CDF derived from an identical analysis of the simulations. A statistical measure of the difference of the two distributions is provided by the Kolmogorov–Smirnov deviation

(31)

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 CDF from simulations. The upper panels show the difference between the observed and median simulated CDF values, , with the grey bands representing the 68.3%, 95.4%, and 99.7% regions of the simulated CDF distributions. Other component separation methods produce similar results, and no significant deviations from Gaussian expectations are observed in the polarization peak statistics.

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 the negative of the mean of all minima, are calculated for the filtered map, and the observed values are compared to the simulated distributions in Fig. 17. The observed minima/maxima counts and means are not significantly different from the fiducial model.

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 Sect. 6.5. They are consistent with the predictions of the ΛCDM model, and our understanding of the instrument and noise properties of the 2018 Planck data.

5.5. Stacking of CMB peaks

5.5.1. 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 small-scale noise and residual systematic effects, thus minimizing the impact of inconsistencies in these properties between the data and simulations.

Hot (or cold) peaks are selected in the CMB intensity map as local maxima (or minima) by comparison with their nearest neighbour pixels, and grouped into different ranges above (below) a given threshold ν (in rms units of the intensity map). In order to facilitate the comparison of the polarization signal between different peaks, we use transformed Stokes parameters, Q_r and U_r, defined by Eq. (2), where the reference point, , is specified by the centre of each extremum . The Q_r component then traces the linear polarization in terms of radial (Q_r >  0) and tangential (Q_r <  0) contributions with respect to the centre of the peak. This stacking method is referred to as “non-oriented”, 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 Planck data allows additional inner rings to be observed compared to WMAP.

We now consider the consistency of the Planck non-oriented 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 χ² estimator is used to quantify the differences between the μ(θ) profiles obtained from the data and the expected values estimated with simulations:

(32)

with the covariance matrix defined as

(33)

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 in PCIS15 as an effect connected with the deficit in the observed power spectrum at low multipoles.

PCIS15 also demonstrated that the χ² 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 an alternative quantity, the integrated profile deviation Δμ_T(W), to evaluate the consistency between the data and the model. This is defined as

(34)

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) that there is a systematic effect in the U_r component due to uncertainty in the orientation of the polarization sensitive detectors.

In summary, given our understanding of the Planck data and simulations, the results from non-oriented stacking are consistent with the predictions of the ΛCDM model.

5.5.2. Oriented stacking

The stacking method of the previous subsection can be generalized by orienting the local coordinate frame of the patch to be stacked in a way that is correlated with the map being stacked. This approach, which we call “oriented stacking”, first introduced in PCIS15, allows extra information to be extracted from the stacked data. For unoriented stacking, the ensemble average cannot result in any intrinsic angular dependence, since it would be averaged by the uncorrelated orientation choices. For example, the ensemble average of the combination of polarization Stokes parameters Q + iU around unoriented temperature peaks has overall angular dependence e^2iϕ, which can be removed by a local rotation (Eq. (2)), as was carried out in the previous section. However, for oriented stacking, the angular dependence is a linear combination of a few Fourier modes, e^imϕ, with the exact mode content determined by the spin of the field being stacked and the spin of the orientation operator. For scalar fields stacked using orientations determined by the spin-2 operator, only m = 0 and m = 2 modes are present.

Choices of what to stack, where to centre the patch, and how to orient it, provide a multitude of statistical tests that are complementary to the auto- and 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 ∇⁻²T, namely the tensor ∇_i∇_j∇⁻²T.

The inverse Laplacian of the temperature ∇⁻²T is readily computed using the spherical harmonic transform

(35)

and its covariant derivatives can be computed using standard HEALPix routines. Perhaps a simpler interpretation of the tensor ∇_i∇_j∇⁻²T is suggested by its projection onto the Stokes parameters using spin-2 spherical harmonics:

(36)

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

(37)

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

(38)

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:

(39)

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,

(40)

can be used to evaluate the consistency between the data and the model. Once again, the weighting function W is chosen to be proportional to the expected profile .

To assess foreground contamination in component-separated 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 low-frequency 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 TE correlation), and no evidence of detectable patterns in the B-mode stack (reaching a noise floor of about 0.05 μK), which serves as a null test for foreground contamination (given the fact that both low- and high-frequency foregrounds display TB correlations, as described below). Polarized dust foregrounds in particular were previously investigated in Planck Collaboration XI (2020). Both low- and high-frequency foreground stacks show no acoustic peak ringing, with oriented temperature stacks being notably different from the CMB. The high-frequency stacks show strong TE and pronounced TB correlations, likely driven by the properties of the polarized dust emission that is dominant at 353 GHz. Low-frequency foreground stacks show weaker, but still noticeable, TE and TB correlations, with no prominent pattern in the E mode, as is the case for high-frequency foregrounds. All foreground stacks are quite distinct from CMB ones, and foreground leakage should result in extraneous correlations in CMB component-separated stacks, which are not observed. Our results indicate that the residual foreground contamination in Planck component-separated maps is below levels that are measurable by this method.

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 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 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 χ² deviation (Eq. (32)) of oriented radial profiles, while Table 10 summarizes the two-tailed p-values of the 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.

6. 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 component-separation effects that does not exist for the component-separated data that we analyse here. Therefore our initial approach is simply to test for evidence of anomalies in the polarization data, in addition to verifying once again those seen previously in the temperature maps.

6.1. Lack of 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. 2015; PCIS15). In particular, we extend the analysis to polarization data, which were previously too noisy and/or contaminated by residual systematic artefacts on such scales to use them for verification of the temperature anomaly. We consider the statistic proposed by Copi et al. (2013) in an analysis of the WMAP data:

(41)

where X, Y can denote the temperature anisotropy T or the two Stokes parameters Q_r and U_r (as defined in Sect. 5.2), and is an estimate of the corresponding 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, θ₁, and upper, θ₂, limits of the separation angle range considered are 60° and 180°, respectively.

To check the consistency of the current data set with previous analyses, we again adopt this range of separation angles and present the results in Tables 11 and 12. We find that the temperature data show a lack of correlation on large angular scales, with a significance consistent with that found by Copi et al. (2015) (although note that the sense of the 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. (2013) demonstrated that the signal-to-noise ratio in the WMAP polarization maps was insufficient to allow meaningful estimates of S^XY to be made. For the Planck polarization data, we have estimated errors on the statistics from the standard deviation of these values computed from the corresponding FFP10 HMHD maps. We find that the errors are σ_S1/2 ≲ 5 × 10⁻⁵ μK⁴ for the statistic, σ_S1/2 ≲ 7 × 10⁻⁶ μK⁴ for the remaining statistics with estimated values at least of order 10⁻⁵ μK⁴, and σ_S1/2 ≲ 5 × 10⁻⁷ μK⁴ for the statistics with estimated values of order 10⁻⁶ μK⁴. 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¹⁰ from the temperature maps. The corresponding probabilities, estimated from similarly-corrected 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 large-angle correlations. Copi et al. (2009) argue that all modes below ℓ ≤ 5 contribute to this, by cancellation with each other and with higher order modes.

A potential criticism of the statistic relates to the a posteriori choice of the range of separation angles to delineate the interesting region of behaviour of the correlation function. As in PCIS15, we consider the generalized statistic S^XY(θ, 180° ) and compute it for all values of θ both for the data and simulations. Figure 23 indicates that the only excursion outside of the 95% confidence regions determined from simulations is observed for the temperature data when the lower bound of the integral is θ ≈ 60°, for all component-separation methods. No equivalent excursions are seen for any statistics involving polarization data, which also show notable variations between component-separation methods. We provide a more quantitative assessment as follows. For each value of θ, we determine the number of simulations with a higher value of S^XY(θ, 180° ), and hence infer the most significant value of the statistic and the separation angle that it corresponds to. However, since such an analysis is sensitive to the look-elsewhere effect, we define a global statistic to evaluate the true significance of the result. Specifically, we repeat the procedure for each simulation, and search for the largest probability irrespective of the value of θ at which it occurs. The fraction of these probabilities lower than the maximum probability found for the data defines a global 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 (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^TT statistic is connected to the observed quadrupole.

As pointed out by Copi et al. (2013), the cross-correlation between temperature and polarization can be used to determine whether the S^TT value is low due to a statistical fluke in the context of a ΛCDM cosmology. Specifically, if this hypothesis were true, the S^TQ statistic would also likely be small on large angular scales, given that the polarization signal is partially correlated with temperature. Copi et al. (2013) have determined that the S^TQ statistic would be smaller than 1.403 μK⁴ 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^TQ(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 can be connected to the temperature 2-point correlation function for antipodean points, . This is motivated by fact that, as seen in Fig. 8, the otherwise largely flat drops to negative values close to θ = 180°. Here, we extend the analysis of the temperature correlation function to the functions including polarization information, .

Table 15 presents the probabilities of finding larger values of the 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.

6.2. 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 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 TT,  TE,  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 χ² statistic defined by Eq. (11), where the correlation function, mean, and covariance matrix are computed for a corresponding hemisphere. If we consider that the χ² statistic itself can act as a measure of fluctuation level, then between the two measured hemispheres can be quantified by the ratio of the corresponding χ² values. The probabilities of obtaining values of the χ² 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 χ² 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 χ² 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 the 2013 and 2015 Planck data sets (PCIS13; PCIS15). This may be a consequence of the use of different masks in the various analyses, or of the improved treatment of poorly determined modes in the estimated correlation matrix used for the computation of the χ² statistic, as described in Sect. 5.2. The largest asymmetry, at a significance of around 97–99%, is observed for the TTQ_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 SMICATU_r 2-point function and the NILCQ_rU_rU_r pseudo-collapsed 3-point function are very similar for the two hemispheres; as a consequence, the probabilities for the corresponding χ² 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 pseudo-collapsed case. Similar behaviour is also observed for some of the component-separated maps for the Q_rQ_rQ_r, TTU_r and TQ_rQ_r pseudo-collapsed 3-point functions, and for a few configurations of the equilateral 3-point function. 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.

6.3. 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 and parity-antisymmetric functions defined as

(42)

where for a given direction , the antipodal point is . Under the parity transformation , it can be shown that the corresponding a_ℓm coefficients transform such that . As a consequence, and are comprised of spherical harmonics with only even or odd ℓ-modes, respectively.

For the polarization field, a similar approach can be followed by expanding the Q and U stokes parameter maps as

(43)

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 and .

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 and to be of comparable amplitude. However, an odd point-parity preference has been observed both in various WMAP data releases (Land & Magueijo 2005a; Kim & Naselsky 2010; Gruppuso et al. 2011), and confirmed in the Planck 2013 and 2015 studies. Furthermore, Land & Magueijo (2005b) have shown that this cannot be easily ascribed to the presence of residual Galactic foreground emission because of its even 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:

(44)

where and are given by

(45)

with being the total number of even (+) or odd (−) multipoles included in the sum up to ℓ_max, and , is the temperature angular power spectrum computed using a quadratic 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 component-separated 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 scale analyses. This map is also used in the Planck low-ℓ likelihood analysis (Planck Collaboration V 2020), 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 TE and EE power spectra. As in the temperature analysis, the data are compared with MC simulations that describe the expected parity preference for the ΛCDM model. However, the ratio estimator cannot be used in this case, since the signal-to-noise 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:

(46)

where X corresponds to either the TE or EE spectra, and is defined analogously to the temperature case in Eq. (45).

In order to minimize the effect of incomplete knowledge of the instrumental noise, and given that some noise mismatch has been observed between the data and simulations (Planck Collaboration IV 2020), we determine the power spectrum using the cross-quasi-QML estimator that forms the basis of the HFI low-ℓ likelihood in polarization (Planck Collaboration V 2020). Here, the signal covariance matrix is constructed according to the FFP10 fiducial model, while the noise is specified using the FFP8¹¹ 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 component-separated 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 TE 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 2020) and their correlation with foregrounds. These effects are not fully captured by the current set of MC simulations. As a consequence, the inclusion of the ℓ = 2 multipole in the 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.

Figure 27 shows equivalent results for the D^EE(ℓ_max) estimator, though with Commander taken as a representative example of the behaviour of the component-separated maps. As in the TE analysis, the data and MC simulations are consistent, although the p-values decrease over the ℓ_max range of 20–30, reaching minima of about 6% at ℓ_max = 27 for the HM data; however, the OE results do not indicate such a trend.

It is interesting to note that the statistic introduced in Sect. 6.1 is related to the measurement of parity asymmetry (see also Muir et al. 2018). In particular, when the functions are expressed in terms of the angular power spectra, it can be seen that the following relationships hold:

(47)

(48)

(49)

(50)

The expected values for the remaining correlation functions are zero, i.e., . 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.

6.4. Peak distribution asymmetry

Localized anomalies on the CMB sky can be searched for by testing how the statistical properties of local extrema (or peaks) vary in patches as a function of location. Since we expect the asymmetry measured by the QML estimator in Sect. 7.2 to be mirrored by a slight difference in the temperature and polarization peak distributions in the corresponding positive and negative hemispheres, we present and examine these quantities here. Note that the test is not powerful enough in itself to determine a modulation direction, but relies on the dipole orientation results presented in Table 24.

We extract peaks from the filtered maps discussed in Sect. 5.4 within the 70° radius discs centred on the positive and negative asymmetry directions determined by the SMICATT,  TE,  EE QML estimator (see Table 24). Separate peak CDFs are then constructed for these directions, then compared to the 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 observed in the negative direction, as previously observed for the temperature results in PCIS15.

6.5. The Cold Spot and other large-scale peaks in temperature and polarization

The Cold Spot was first detected in the WMAP temperature data as an anomalously cold region in terms of the spherical 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; PCIS15). 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¹². Its internal structure was then subsequently studied by Chiang (2018) in the context of an analysis of the variation of the acoustic peak positions in different parts of the sky. It was found that the Cold Spot region shows a large synchronous shift of the peaks towards smaller multipoles.

In this section, we analyse the polarization pattern of the Cold Spot, considered as a minimum in the temperature field when smoothed with a Gaussian of standard deviation R = 5°, to search for evidence as to whether its origin is primordial or otherwise. Vielva et al. (2011) and Ferná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. 2017a,b, for more details). Since these peaks dominate the temperature field on large angular scales, the particular value of the smoothing scale used in their identification is not critical for locating the peaks. The scale R = 10° is chosen in order to highlight these structures while preserving their geometrical properties, such as curvature and eccentricity. These features correspond to two maxima and two minima, shown as Peaks 1–4 in Fig. 29, all located in the southern Galactic hemisphere. Given their locations, they are expected to make a significant contribution to the 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 TB correlation and neglecting the eccentricity of the peaks, contain all of the temperature and polarization information that we are interested in.

In order to have a complete characterization of the peaks, we have to calculate derivatives of the fields up to second order. Since we are interested in the azimuthally symmetric signal, only the peak height and the Laplacian are relevant to the Q_r profiles. Once these two quantities are calculated from the temperature smoothed at a given scale R, their values are used to construct the expected theoretical profile in the polarization field (see 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 cosmic-variance 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:

(51)

(52)

Here and X(θ_i) are the measured profile and the corresponding theoretical expectation for the angular bin centred at θ_i, respectively. In these equations, X represents the T or the Q_r profile, depending on the case we are analysing. Since we are interested in 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 (Vielva et al. 2011; 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 model in Fernández-Cobos et al. (2013) was calculated using only the Laplacian (i.e., SMHW coefficient). Moreover, the Cold Spot amplitude is characterized with respect to the observed standard deviation, not using the theoretical expectation, as in our. This difference in the theoretical models results in different significance levels for the forecasts. However, the normalization of the theoretical profiles is not relevant when the significance of the hypothesis (A/σ_A) is calculated from the observed amplitude A.

The T and Q_r profiles determined from the four component separation methods are shown in Fig. 30. The corresponding profile amplitudes are given in Table 20. The temperature results are compatible with the analysis of the same peaks performed in 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.

7. Dipole modulation and directionality

In this section, we examine dipolar violations of isotropy. In PCIS15, we performed a non-exhaustive 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¹³. 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.

7.1. 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, 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 p-values 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 full-resolution 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 low-resolution 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 simulations. The analysis shows that the direct application of the method to isotropic simulations of E-mode polarization signal returns local-variance 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-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:

(53)

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 component-separation 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 component-separation 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 component-separation 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.

7.2. Dipole modulation: QML analysis

In this section, we use the QML estimator originally introduced in Moss et al. (2011) to assess the level of dipole asymmetry in the CMB sky, and further extend the analysis to polarization data. We note, however, that Contreras et al. (2017) found that Planck polarization data are unlikely to be able to distinguish between dipole-modulated skies and skies consistent with ΛCDM (although this statement is somewhat model dependent). Furthermore, the analysis we carry out here is a purely phenomenological one performed in multipole space, with no attempt to connect to any real-space modulation; however, several physical k-space models are considered in a companion paper (Planck Collaboration X 2020).

We employ a version of the estimator proposed in PCIS15, which was further developed in Zibin & Contreras (2017) and Contreras et al. (2017). In particular, we use

(54)

(55)

with

(56)

Here the δC_ℓℓ′ are determined by the model of modulation; for this case we have assumed scale-invariant power modulation out to some ℓ_max and thus δC_ℓℓ′ = 2(C_ℓ + C_ℓ′). The remaining terms are: , the spherical harmonic transform of the modulation amplitude and direction; WZ= TT, TE, or EE; W_ℓm and Z_ℓm, which are inverse-covariance filtered data; ; and the last term on the right-hand-side of Eq. (56) denotes the mean-field correction (details, including the precise form of can be found in Appendix A.1 of Planck Collaboration XV 2016 and Appendix B of Contreras et al. 2017). The parentheses in the superscripts indicate symmetrization over the enclosed variables. The coupling coefficients are given by

(57)

(58)

The can be transformed into amplitudes of modulation along each Cartesian axis, since they are simply the spherical harmonic decomposition of a dipole. Likewise we can write the modulation in terms of its spherical coordinates (amplitude and direction) as

(59)

(60)

(61)

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 (Contreras et al. 2017). This has little effect on the comparison of the data with simulations, however, since the simulations are treated in the same way. A potentially more significant effect is a mismatch in the optical depth between data and simulations; below we estimate this to have essentially negligible influence.

The estimators of Eqs. (54) to (55) can be combined with inverse-variance weighting over all data combinations (TT, TE, EE) to obtain a combined minimum-variance estimator, given by

(62)

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 (Contreras et al. 2017). 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 TT, TE, and EE has no physical basis; a more physical, k-space approach is considered in the companion paper Planck Collaboration X (2020).

In Fig. 34 we show the p-values of the full-mission data compared to the FFP10 simulations, considering TT (top left), EE (top right), or TE (bottom left) data independently. While the TT results are largely consistent with our previous analysis (differences are within the expected scatter given the different analysis masks, PCIS15), our EE and TE results are new. The 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 temperature-polarization 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 TE) 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 Zibin & Contreras 2017, and references therein).

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 TE 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 Eq. (50) in PCIS15). In Table 25 we show the amplitude and direction of dipolar modulation for the combined TT, TE, and EE data in the ℓ = 2 − 220 range, which is the range with the lowest 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 low-ℓ anomaly 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 p-values (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 temperature-modulation signature.

7.3. Angular clustering of the power distribution

Despite the lack of evidence for any strong anomaly in the amplitude of dipole modulation discussed in the previous 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 band-power estimates, it was found that the directions were aligned at the 2–3σ level up to ℓ  =  1500 when compared to simulations. In the standard cosmological model, although such maps are expected to exhibit dipolar distributions of power due to Gaussian random fluctuations, the associated directions should be independent random variables. Evidence for the close correlation and alignment of directions on different angular scales then appears to be a signature of broken statistical isotropy. Since we do not observe a significant amplitude of dipole modulation over similar angular scales, then the result of finding clustering of directions seems mysterious – it is hard to imagine a concrete (e.g., inflationary) model that would cluster the directions without affecting the amplitude of the modulation (but see Hansen et al. 2019). Nevertheless, regarding this as a purely empirical question, it is important to repeat the directional-clustering analysis and broaden its scope to include polarization.

Here, we repeat the analysis using the Planck 2018 data and extend it to include polarization measurements. The local power spectra, TT, TE, and EE, are estimated directly from maps of the temperature and Stokes Q and U parameters. Since the analysis is sensitive to differences between the noise properties of the data and simulations, we consider only cross-spectra computed between the two HM or OE maps for each component-separation method.

We adopt the same approach for the estimation of the dipole alignment as described in detail in PCIS15, a brief summary of which follows.

1. Local TT, TE and EE power spectra are estimated (using the MASTER approach, Hivon et al. 2002) directly from the T, Q, and U maps at N_side = 2048 for the 12 patches of the sky corresponding to the N_side = 1 HEALPix base pixels. Leakage between the E and B modes due to incomplete sky coverage is removed in the power-spectrum estimation during the inversion of the full TEB kernel. The same results would be expected using the E maps, but the extended masks required for these maps would increase the uncertainty of the results. The spectra are binned over various (even) bin sizes between Δℓ = 8 and Δℓ = 32¹⁴.

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 (PCIS13; PCIS15), the fitted dipole amplitudes are found to be consistent with those determined from simulations.

4. A measure of the alignment of the different multipole blocks is then constructed. We use the mean of the cosine of the angles between all pairs of dipoles. This is essentially the Rayleigh statistic (RS) and we will refer to it as such, although it differs by not including any amplitude information. Smaller values of the RS correspond to less clustering.

5. The clustering as a function of ℓ_max is then assessed using 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 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 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 TT, TE, 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 − 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 component-separated 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 p-values 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 TE. 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 TE data is exceeded in 20–40% of the simulations (for the example of Commander, for the various masks and data splits), and hence does not appear to be anomalous.

In order to further study the alignment in EE polarization for ℓ <  300, we tested whether the directions of the EE dipoles in this range are correlated with the directions for the TT dipoles. Here we made a small change in the statistic as already described: in point 4 in the above description of the method, we instead use the mean of the cosine of the angles between all pairs of dipoles, where one dipole is always taken from EE and the other always from TT. In this way, the statistic measures the cross-correlation between TT 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 TT and EE dipoles are detected for the same multipole range where the EE polarization dipoles (and TT dipoles) are clustered.

We have seen that the EE directions appear clustered for ℓ_max ≈ 200. This is also the case for TT, and so we would like to test whether the two preferred directions are also related. What Fig. 40 shows is that not only are the TT and EE dipoles clustered among themselves, but that they also appear to be clustered towards the same direction. If TT and EE were independent, this would be highly unexpected, even if TT and EE were separately clustered. However, we know that the TE spectrum is non-zero, giving rise to a correlation between T and E. In order to test whether such a TT–EE directional correlation is expected in the case when both TT and EE are clustered individually, we perform the following test. We examine all simulations having a minimum p-value of less than 1% for both TT and EE in overlapping multipole ranges (similar to what we observed in the data). Only two simulations have overlapping TT and EE multipole ranges using this criterion, but neither of them has a significant correlation between TT and EE directions. While this may seem to suggest that the effect in EE is not just due to the already known directional clustering in TT (and the TE correlation), it is hard to draw firm conclusions without many more simulations.

We investigate this further in Fig. 41 where we show the dipole directions for blocks of 10 multipoles for ℓ <  200 for TT, TE, and EE. There is some hint here of the correlation between TT and EE directions seen in Fig. 40. Note that in Fig. 40 this angular correlation appears stronger for the HM split than for the OE split; given the large errors on direction for polarization, differences in the noise properties (including systematic effects) for the HM and OE split may be responsible. The middle panel in Fig. 41 shows the TE directions. Clearly the TE directions are more scattered than TT and EE, as expected from Fig. 39. The directions for several multipole blocks in the range ℓ <  100 coincide quite well with the TT and EE directions, but most of the 10-multipole blocks point in different directions. To further compare the angular clustering of TE dipoles with TT and EE dipoles, we construct the statistic measuring the correlation between TE and TT directions, as well as between TE and EE directions, in the same way as described above for the correlation between TT and EE. We find no significant correlations between the TT and TE directions, but the correlation between TE and EE again shows some sign of similar multipole ranges where EE is aligned and the TT and EE directions correlate. This is shown as a magenta line in Fig. 40. The correlation between TE and EE directions at ℓ ≈ 200 may arise as a result of the fact that TT and EE are aligned, as well as the fact that TT and EE directions are also correlated at exactly these multipoles. However, since none of the simulations show a high significance for all three of these statistics in one common multipole range for one given simulation, we are unable to test this with simulations. Nevertheless, we might have expected a similar correlation between TT and TE if this was really the case.

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 simulations 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 TE, the numbers are completely consistent with simulations, while for EE, the results are unstable to the choice of masks and component-separation 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 significance can be set. A similar conclusion can be made for the TT/EE angular correlation, but the results are in this case more stable with choice of mask and component-separation method.

8. Conclusions

In this paper, we present a study of the statistical isotropy and Gaussianity of the CMB using the Planck 2018 temperature and polarization data. The Planck 2015 release essentially corresponded to the limit of our ability to probe CMB anomalies with temperature fluctuations alone. The use of 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 (2020). 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 signatures are well known; thus, in this summary, we focus on the properties of the polarization data.

In Sect. 5, we examine aspects of the Gaussianity of the polarized CMB fluctuations. Tests of 1D moments, N-point correlation functions, Minkowski functionals, peak statistics, and the oriented and unoriented stacking of peaks yield no indications of significant departures from Gaussianity. In addition, no evidence is found for a low variance of the polarized sky signal.

Section 6 provides an updated study of several previously known peculiarities. We find no evidence in the polarization data of a lack of 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 alignment of the preferred directions of the temperature and E-mode 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 dipole-modulation 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–TT 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 low-noise polarization experiments.

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 2020), 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.

---

Acknowledgments

The Planck Collaboration acknowledges the support of: ESA; CNES, and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, JA, and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planck-collaboration. This work has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement numbers 687312, 776282, and 772253.

References

Appendix A: E- and B-mode map reconstruction

A.1. Introduction and notation

Polarization of the CMB is usually measured on the sky in terms of the Stokes parameters Q and U (see e.g., Hu & White 1997). However, these are not scalar quantities and therefore not rotationally invariant. Nevertheless, scalar E and 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 1997; Bunn et al. 2003), defined by

(A.1)

They are related to the standard E and B estimators in harmonic space by a factor of

(A.2)

which corresponds to an application of the bi-Laplacian operator such that ℰ, ℬ = −∇²(∇² + 2) ψ_E, B in real space. Unlike E and B, ℰ and ℬ 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 B-mode maps from partial sky coverage is that neither spin-0 nor spin-2 spherical harmonics are orthogonal under the masked inner product,

(A.3)

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 are given by

(A.4)

where the coupling matrices K_ℓ1ℓ2 depend on the method used, as illustrated in Fig. A.1.

For reconstruction of low-ℓ multipoles, maximum-likelihood 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 et al. 2003; Bunn 2011; Bunn & Wandelt 2017). Some of the approaches discussed in the literature require solving large linear algebra problems (typically via iterative solvers), and could be expensive on 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 Appendices A.2, A.3, and A.4. The suitability of these approaches for our purposes depends upon the uniformity of the reconstruction, since the 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 defined by a window function :

(A.5)

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 ℰ and ℬ 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).

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 nearest-neighbour 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 finite-difference “Laplacian stencils”, as detailed in Appendix A.3.2, and we further develop the details of the finite-differencing multigrid approach in Appendix A.3.3.

One particular aspect of our improved method is that when computing the average over the nearest neighbours of a given masked pixel, their contributions are Gaussian-weighted 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 ∇²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 Appendix 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 ℰ- and ℬ- 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 second-order-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:

(A.6)

with real c_i for scalar and complex for tensor values, and being the average area of a pixel. We will derive the stencil weights in a 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,

(A.7)

associated with it, which form a basis on a tangent plane and a normal direction . Nearby points can be projected onto a selected tangent plane via gnomonic projection, which maps the radius vector r into

(A.8)

and thus introduces local coordinates on a tangent plane,

(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 used to define the Stokes parameters Q and U, i.e.,

(A.10)

their values transform as

(A.11)

In the HEALPix polarization convention, the Stokes parameters Q and U of a pixel are always defined with respect to the basis of the pixel, and must be rotated to the basis when projecting onto a tangent plane. The appropriate angle of rotation can be computed as

(A.12)

with transformation of Stokes parameters most conveniently implemented as a complex phase rotation

(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 ℒ[const.] = 0 and ℒ[ρ²] = 4, which leads to

(A.14)

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.,

(A.15)

where the width σ can be tuned. We chose it to be σ = 1.61, which gives near perfect residual cancelation at the poles. This is the scalar Laplacian stencil we will adopt in what follows, while the complex Laplacian stencil for polarization is defined by , as explained above.

A.3.3. Multigrid methods

Inpainting a map ϕ can be viewed as a diffusive flow ∂_tϕ = ∇²ϕ 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 updates a pixel according to

(A.16)

using the weighted sum of itself and its neighbours, where α = δt/h², and the largest time step δ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₀) becomes negative, so in practice the fastest diffusion is often achieved by replacing a pixel by a weighted sum of its neighbours:

(A.17)

N successive iterations will diffuse the solution across roughly N^1/2 pixels, so running diffusion flow like this for full convergence to the static solution ∇²ϕ = 0 is very expensive on high-resolution maps. In previous studies (Planck Collaboration XXIV 2014; Planck Collaboration XVII 2016), a finite number of steps like Eq. (A.17) was applied to inpaint the masked areas.

Convergence of diffusion flow can be accelerated if one notes that a coarser spatial grid allows bigger time steps, and thus the solution of a linear elliptic boundary value problem on a grid with spacing h,

(A.18)

can be found iteratively from the approximate solution by downgrading the residual

(A.19)

to a coarser grid with spacing 2h, thus obtaining the coarse-grid correction

(A.20)

which can be upgraded back to a fine grid with spacing h and used to improve the solution, i.e.,

(A.21)

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 w-stroke 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 w-stroke are enough to reach the static solution to double precision. This requires only 288 diffusion steps at full resolution, and is substantially faster than the 2000 iterations at full resolution that were used in Planck Collaboration XXIV (2014). Improved performance of the inpainting routine allows for further applications and testing to be carried out with the same computing resources.

A.4. Purified inpainting

A.4.1. Method description

Here, we demonstrate that the advantages of masking and inpainting can be combined into a single method.

The motivation for the approach is to construct a full-sky polarization map, , from a masked one, M ⋅ p, such that it agrees identically with p on the unmasked portion of the sky, and assigns modes either through strict attribution to the E and B subspaces of the full-sky map, or “ambiguous” attribution where some mode mixing is allowed:

(A.22)

Note that the polarization map components e and b are not pure in the sense of being orthogonal to the entire E- and B-mode linear spaces on the masked sky, as in Bunn et al. (2003); that requirement results in a large linear algebra problem, which is expensive to solve, although efficient methods for that have been developed recently (Bunn & Wandelt 2017). Instead, their purpose is to “purify” the ambiguous mode map a, ensuring that most of the pure modes are projected out, thus minimizing the power leakage between E and B in the ambiguous mode.

Initially, the entire polarization map p is assigned to an ambiguous mode:

(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:

(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

(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

(A.26)

Once the new mode is constructed, it is normalized and the most recent basis set is updated:

(A.27)

Lanczos bi-orthogonalization guarantees that the dot product of the constructed basis vectors ⟨v_i, v_j⟩ forms a tri-diagonal 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.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, ∇²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 𝒪(n log n). The final reconstructed polarization map is generated by inpainting the remaining ambiguous modes as

(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 and . Various metrics of performance can be assigned to the masked residual maps and . Figure A.4 compares the residual maps for one SMICA realization from the FFP10 simulations, as reconstructed via each of the masking, simple inpainting, and 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 reconstruction process itself can help to reduce the residuals levels yet further, in particular given that these residuals tend to be localized near the mask boundary. To systematically define the optimal confidence mask, we performed purified-inpainting 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.5, 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² + δB² 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 to be undertaken concerning the isotropy and statistical properties of the Planck polarization data.

All Tables

All Figures