INTRODUCTION

The cosmic microwave background (CMB) anisotropy on large angular scales has been measured with increasing precision over the past few decades by 3 fullsky experiments, namely, the Cosmic Background Explorer (COBE), the Wilkinson Microwave Anisotropy Probe (WMAP), and the Planck satellite. The CMB data from these experiments largely agree with the six-parameter standard ΛCDM model. However, several notable deviations at large angular scales in the temperature data have been identified. These include the lack of large-scale correlations (Hinshaw et al. 1996;Bennett et al. 2003;Copi et al. 2009), alignment of low multipoles (Tegmark et al. 2003;de Oliveira-Costa et al. 2004;Copi et al. 2004), point-parity asymmetry (Land & Magueijo 2005;Kim & Naselsky 2010), hemispherical power asymmetry (Eriksen et al. 2004;Monteserín et al. 2008;Akrami et al. 2014), and cold spots (Vielva et al. 2004;Cruz et al. 2005). While the lack of large scale correlations with the associated low quadrupole were ob-served by COBE (Hinshaw et al. 1996), the majority of these temperature deviations were first found in WMAP data (Bennett et al. 2011), and then were confirmed by Planck (Planck Collaboration VII 2020). Although the significance of these large-scale deviations is challenging to pin down (e.g, Bennett et al. 2011), throughout this paper we will refer to the features by their common name: "CMB anomalies". For a review of CMB temperature anomalies, see Schwarz et al. (2016).

For each anomaly, various statistical measures have been formulated to quantify its significance in a standard ΛCDM universe. Most of the statistical tests show mild (2-3σ) deviations compared to the ΛCDM model, and the a-posteriori nature of the observations requires that care be taken to avoid overstating the significance level. Consequently, there are three most commonly accepted postulations as to why these anomalies are present in the data. First, they could have cosmological origins. For example, there have been several attempts to use "justenough" slow-roll inflation to explain the lack of correlations on large scales (Ramirez & Schwarz 2012;Cicoli et al. 2014), and introducing modulation fields to the primordial curvature fluctuation spectrum P (k) can explain the hemispherical asymmetry and the quadrupoleoctupole alignment (Gordon et al. 2005;Dvorkin et al. 2008). The cosmological origin hypothesis is the most exciting explanation because it implies new physics beyond the standard ΛCDM model. However, these mechanisms have their corresponding weaknesses: it might not be worth introducing more parameters to explain a single or a few 2-3 sigma level anomalies; the modulation fields are also hard to physically motivate (Schwarz et al. 2016). A second postulate is that these anomalies could result from foreground effects due to the solar system, the Milky Way, and the local supercluster; or systematic effects in the data (see e.g. Hansen et al. 2004;Abramo et al. 2006;Copi et al. 2007;Groeneboom & Eriksen 2009;Hansen et al. 2012;Frejsel 2015). The primary challenge to these hypotheses is that the same anomalies are observed in both Planck and WMAP, which have different sky-scanning strategies, and rely primarily on different frequency bands to estimate the CMB. The third prevailing view is that the anomalies are merely unlikely fluctuations in our universe (Bennett et al. 2011). A better understanding of the anomalies will either support this "fluke hypothesis" or deepen our understanding of the universe. But since CMB temperature measurements have reached the cosmic-variancelimit at large angular scales, information from additional observations is needed.

To extend this investigation, studies have explored a range of cosmological probes, including 21-cm obser-vations (Shiraishi et al. 2016;Li et al. 2019), gravitational lensing (Yoho et al. 2014;Mukherjee & Souradeep 2016;Zibin & Contreras 2017), the integrated Sachs-Wolfe effect (Muir & Huterer 2016;Copi et al. 2016), and the Sunyaev-Zel'dovich effect (Cayuso & Johnson 2020). In this paper, however, we focus on CMB polarization measurements. In particular, the dominant Emode polarization component promises to provide the next most significant signal after the CMB temperature anisotropies. The Planck collaboration recently searched for the existence of such potential anomalies using their polarization data, but due to the low signalto-noise at large angular scales, different componentseparated maps resulted in different significance levels (Planck Collaboration VII 2020). Progress has been made by using single-frequency Planck polarization maps (Chiocchetta et al. 2021). Higher sensitivity observations can provide more details about these phenomena.

The situation may improve in the near future with new data from CMB polarization experiments that targeting the large angular scales over greater than 30% of the sky (e.g., Essinger-Hileman et al. 2014;Lee et al. 2020;Pérez-de-Taoro et al. 2016;Addamo et al. 2021;Gandilo et al. 2016;Hazumi et al. 2020;Kogut et al. 2011;Hanany et al. 2019). Recent studies forecast the constraining power of future polarization experiments on individual anomalies (e.g., Billi et al. 2019;O'Dwyer et al. 2020;Chiocchetta et al. 2021). It is also possible to separate polarization data into parts that are fully correlated or uncorrelated with temperature, as explored in Frommert & Enßlin (2010). Exploration of multiple anomalies and their correlations has been done in temperature (Muir et al. 2018), but not in polarization.

In this paper, we evaluate the significance of anomalies in temperature and polarization taking into consideration the correlation among anomalies in temperature and polarization. We first analyze Planck data to establish baseline constraints and validate our pipeline using existing results. Next we forecast constraints from the ground-based Cosmology Large Angular Scale Surveyor (CLASS) (Harrington et al. 2016;Dahal et al. 2022) and from LiteBIRD (Hazumi et al. 2020). Both CLASS and LiteBIRD anticipate sample-variance-limited E-mode measurements. Therefore, the primary difference between CLASS and LiteBIRD in this paper is the partial sky coverage of CLASS (∼ 70% of the sky fraction of LiteBIRD). The four CMB large-scale anomalies we considered in our work are: the lack of large-angular correlations, the alignment of the quadrupole and octupole, the point-parity asymmetry, and the hemispherical power asymmetry.

The structure of this paper is as follows: in Section 2, we describe the maps, masks, and simulations we used. Section 3 introduces the different large-scale anomalies and their statistics. We present our results in Section 4 and conclude in Section 5.

METHODS

Most analyses in this work are based on the temperature map T ( n) and maps of the linear polarization Stokes parameters Q( n) and U ( n). They can be expanded with spherical harmonics Y ℓm ( n) and spin-2 spherical harmonics ±2 Y ℓm ( n) as

where a T ℓm 's and a

ℓm 's are the expansion coefficients (Zaldarriaga & Seljak 1997).

For the polarization analysis of the hemispherical power asymmetry, we constructed the rotationally invariant E-mode map, which can be expanded with spherical harmonics as

The expansion coefficients a E ℓm 's are related to a

ℓm 's by

Denoting temperature and E-mode maps as X( n) where X ∈ {T, E}, the variance and covariance of their expansion coefficients give the power spectra C XY ℓ . In the ΛCDM model, this can be written as

where • refers to taking the ensemble average and XY ∈ {T T, T E, EE}.

The temperature two-point correlation function is defined as C T T (θ) ≡ T ( n1 )T ( n2 ) with n1 • n2 = cos θ, and is related to the power spectra (using equations defined above) as

where P ℓ 's are the Legendre Polynomials. For the Emode, as pointed out in Yoho et al. (2015), the physical interpretation of the two-point correlation function defined based on Equation 2 is ambiguous because an integral over the full sky is required to extract the a E ℓm information. We followed the path in Yoho et al. (2015) and defined the correlation function with local E-modes ( E( n)) that can be calculated from the Stokes Q and U using local spin raising and lowering operators (Zaldarriaga & Seljak 1997;Kamionkowski et al. 1997):

The corresponding E-mode correlation function

The estimation of power spectra will be explained in Section 2.3.

Data and Masks

We used several data products from the Planck 2018 data release.foot_0 We used the Planck SMICA component separation maps, both full mission and half mission (HM) mapsfoot_1 , in our analyses. For analyses based on other component separation methods (Commander, NILC, and SEVEM), we refer the reader to Planck Collaboration VII (2020). In this paper, we used Planck SMICA. Using different component separated maps does not affect our main conclusions. The Planck maps we used in the paper were originally created at HEALPix (Górski et al. 2005) resolution NSIDE = 2048, at an approximate angular resolution of 5 full width half maximum (FWHM). Since we are interested in large angular-scale structures, the original versions were downgraded to lower HEALPix resolutions at NSIDE = 64 and 16, with pixel resolutions being ∼ 55 and ∼ 220 , respectively. The NSIDE = 64 maps were used for map-based analyses of the quadrupole-octupole alignment and the hemispherical power asymmetry. For computational efficiency, NSIDE = 16 maps were used for computing power spectra at low multipoles. Power spectra were used in evaluating statistics for anomalies associated with the lack of large scale correlations and the point-parity asymmetry.

The original maps were downweighted in harmonic space according to (Muir et al. 2018;Planck Collaboration VII 2020) where b ℓ 's are Gaussian smoothing functions and p ℓ 's are HEALPix window functions. Superscripts (in) and (out) stand for input and output maps, respectively. For the Planck data, the function b

was the Gaussian approximation to the Planck beam transform. The Gaussian transform b (out) ℓ set the output smoothing scale and prevented aliasing of small-scale (high-ℓ) noise to larger scales when downgrading the resolution of the map. To downgrade maps to NSIDE = 64 and NSIDE = 16, the output Gaussian smoothing scale was chosen to be FWHM = 160 and FWHM = 640 , respectively. For power spectrum analyses of Planck polarization data, instead of using Equation 8, we smoothed maps with the filter from the low-ℓ Planck polarization analysis Planck Collaboration Int. XLVI (2016); Planck Collaboration V (2020):

Finally the downgraded resolution (NSIDE) map is constructed from the set of smoothed a

(out)

ℓm . The E-mode maps are created from NSIDE=64 maps directly, following Equation 2.

We used different masks in our work for different purposes. For the full-sky analysis, we used the 2018 version of the Planck common masks 3 for temperature and Q/U polarization. These masks were originally generated at NSIDE=2048. To obtain binary masks at lower resolutions, we first smoothed and downgraded them according to Equation 8, and then set pixel values which were smaller than 0.9 to be 0 (masked), otherwise 1 (unmasked). For the partial-sky coverage case, we adopted the survey limits for CLASS. Specifically, we used the same common masks and survey declination limits of -70 • < δ < 30 • . When dealing with Planck HM data and noise simulations, we additionally combined the Planck HM missing pixel masksfoot_3 with the Planck common masks. We first combined the HM missing pixel masks with common masks at NSIDE=2048 then processed the combinations in the same way as for the common masks. To reduce impacts from E/B mixing induced by constructing E-mode maps from a set of a ℓm 's computed on the masked sky, we created E-mode masks for the hemispherical power asymmetry analyses by setting thresholds on the residual E/B signal. Details can be found in Appendix A.

The temperature and Q/U polarization maps together with the masks we used in this work are shown in Fig- ure 1. E-mode masks are given in Appendix A. The sky fractions of all the masks used in this work are listed in Table 1. b For Q-O alignment, we did not mask out the galaxy for the full and HM case, and only adopted the declination limit for partial sky case. The corresponding sky fraction is 70. More details can be found in Section 3.2. c We combined the corresponding Planck HM missing pixel masks with the common masks when processing the masks.

Simulations

We used 4 different levels of noise simulations to estimate the 95% confidence intervals for the statistical measures of anomalies. To indicate the constraining power of Planck data, we used 300 realizations of Planck's end-to-end noise simulationsfoot_4 corresponding to the SMICA foreground separation method. These noise simulations were downgraded to NSIDE=64 (16) using the same method as applied to the data. To estimate the constraining ability of CLASS and LiteBIRD, we used white noise at a level of 15 µK • arcmin (Essinger-Hileman et al. 2014)foot_5 , and 2 µK • arcmin (Hazumi et al. 2020), respectively, and generated 300 simulations for both. Despite its design to mitigate spurious polarization signals from the atmosphere and the instrument (Harrington et al. 2021), the sensitivity of the CLASS maps at large scales will deviate from the 15 µK • arcmin white noise assumption, the extent of which is still being determined. Moreover, we present CLASS forecasts for temperature with the same noise prescription despite inaccessibility of the largest scales in temperature from the ground due to the correlated brightness fluctuations from the atmosphere. The CLASS temperature results should be taken as a toy case to indicate the impacts of partial sky coverage. Finally, to compare with the cosmic-variance-limited situation, we added 2 µK • arcmin white noise to temperature maps and 0.02 µK • arcmin white noise to polarization maps. Noise is needed for taking the xQML power spectra (Section 2.3), and the impact on anomaly estimators is negligible at these noise levels. For these cosmic-variance-limited simulations we do not apply any masks, and this specific case is referred to as the "Ideal" case. All white noise simulations were generated at NSIDE=64 (16) and smoothed with a 160 (640 ) FWHM Gaussian smoothing function to reproduce the same filtering as in Equation 8.

We simulated 10 4 Gaussian random CMB realizations. T QU map sets were produced at NSIDE=64 and 16 using the theoretical power spectra based on Planck's bestfit ΛCDM model 7 . The maps were smoothed according to Equation 8 or 9 so that the properties of the simulated CMB maps are consistent with that of the real data products and/or the noise simulations. Finally, we randomly combined the 300 noise simulations for Planck, CLASS, and LiteBIRD with the 10 4 CMB realizations, forming 3 sets of 10 4 signal+noise realizations for our analyses. The combinations are not completely inde-pendent from one another as the number of the noise simulations is much smaller than that of the CMB realizations. We denote the analyses in which the simulated CMB temperature and polarization signals are allowed to fluctuate freely according to the Planck's bestfit model as the 'unconstrained universe' study (Section 4.1). For studies based on constrained simulations, see e.g., Copi et al. (2015a); Chiocchetta et al. (2021).

In addition to the unconstrained universe, we considered two more simulation sets. The first is the special CMB realization (Special CMB). The Special CMB was selected from 10 5 unconstrained CMB simulations such that each E-mode anomaly estimator's probability-toexceed (PTE) with respect to the Ideal case simulations is either above ∼ 95% or below ∼ 5%. Tests based on the Special CMB distinguish the contribution of survey noise from that of cosmic variance in the posterior distributions of polarization anomaly estimators. More details about the Special CMB can be found in Section 4.2 and Appendix B. The second alternative simulation considered is the constrained universe, in which we fix the temperature simulations to be what has been measured from Planck SMICA data and only sample E and B components as Gaussian random fields according to Planck's best-fit ΛCDM model. When generating the constrained simulations, the temperature a ℓm 's were derived from the full sky data. We checked the impact from foreground residuals by comparing the anomaly estimator distributions (Figure C.1) constrained by either the full-sky SMICA temperature (our baseline analysis) or SMICA temperature with the Galactic regions inpainted as in Planck Collaboration IV (2020). The differences are subtle, which means the existing residual foreground won't significantly impact the conclusions drawn with constrained simulations. Looking into distributions from the constrained simulations is useful in making forecasts given the existing temperature measurements. We found that the predictions from the constrained universe (Appendix C) are largely consistent with those from the unconstrained universe.

xQML Power Spectra

We used xQML to estimate the power spectra. xQML is a quadratic, maximum-likelihood-based power spectrum estimator based on cross-correlation between maps (Vanneste et al. 2018). A primary advantage of crossspectra between maps with uncorrelated noise is that such spectra contain no noise bias. Compared to pseudo-C ℓ methods (e.g., Wandelt et al. 2001;Chon et al. 2004), quadratic maximum-likelihood estimators have the advantage of being optimal at large angular scales (lowℓ) (Tegmark 1997;Bond et al. 1998 Oliveira-Costa 2001). One disadvantage of xQML spectra is their computational inefficiency limits the ℓ-range.

All cross spectra in our work were computed based on NSIDE=16 maps, and the maximum ℓ used in our analyses is 30. Quadratic maximum-likelihood estimators require an estimate of the noise covariance matrix (NCVM). For Planck, the NCVM for temperature was computed based on 300 SMICA HM noise simulations smoothed and downgraded to NSIDE=16 in the same way as for other Planck data. The NCVM for estimating EE power spectra was computed based on the 300 SMICA HM noise simulations. We found constructing the estimator with only the diagonal component of the NCVMs produced more self-consistent and numerically stable results across the thousands of simulations performed. This estimator was then used to compute the cross spectra between Planck HM1&2 maps for both the data and the simulations. Although the diagonal-only NCVMs are not optimal, using them does not bias the resulting spectra, which are still more precise than their pseudo-C ℓ counterparts. Similarly, we used only the diagonal part of the covariance of the 300 white noise simulations for CLASS and LiteBIRD. We also adopted Planck HM masks to eliminate the impact of missing pixels in Planck High Frequency Instrument (HFI). As a check that the simulations were representative of the data, we computed PTE values, defined as the percentage of simulations with power greater than that from the data at each multipole used in this work. Figure 2 shows that the resulting PTE distribution is largely uniform, meaning that the xQML power spectra estimation of the simulations are consistent with that from the data. We did not include the quadrupole for EE analyses because Planck MC simulations are not able to characterize the residual systematics at ℓ = 2 (Planck Collaboration VII 2020).

Throughout the paper we used the constructed estimator to estimate the cross spectra between two maps with the same signal but different noise components.

ANOMALIES STUDIED

In this section, we give a brief introduction to the four anomalies considered in this work. Tables 2, 3, and 4 summarize masks and noise simulations we used for each anomaly, and are provided in the corresponding subsections.

Lack of Large-Angular Correlations

We begin with the lack of temperature correlation on large angular scales: the observed temperature twopoint correlation at large angular scales is much closer to zero than predicted from the ΛCDM model (Figure 3). This anomaly was first measured by COBE (Hinshaw et al. 1996), then confirmed by WMAP (Bennett et al. 2003) and Planck (Planck Collaboration VII 2020). Similar tests for, e.g., T Q, QQ, U U components can be found in, e.g., Copi et al. (2013); Yoho et al. (2015); Chiocchetta et al. (2021).

To quantify this anomaly, we chose the statistic in Spergel et al. (2003); Planck Collaboration VII (2020):

where C XX (θ) is the two-point correlation function, and X ∈ {T, E}. The integration lower and upper bounds were chosen to be -1 (θ max = 180 • ) and 1/2 (θ min = 60 • ), respectively, because this combination approximately reflects the greatest possible deviation for the temperature data (Planck Collaboration VII 2020).

The same integration bounds were adopted for EE for consistency. Instead of doing numerical integration on the two-point correlation function directly, we construct S XX 1/2 based on the power spectra by using Equations 5 and 7, which can be expressed as (Yoho et al. 2015;Chiocchetta et al. 2021):

[deg] where I ℓℓ is defined as

The exclusion of ℓ = 2 does not affect the statistics of S EE 1/2 significantly because the correlation functions do not strongly rely on C EE 2 (Chiocchetta et al. 2021). We used different values of ℓ max for temperature and polarization when computing the estimator for the lack of white noise 2 µK a We used the same configuration for these two anomalies.

b For both intensity and polarization if not mentioned explicitly. c HM refers to Planck half mission missing pixel masks; common refers to Planck common masks; and CLS refer to declination limit for CLASS. See Figure 1.

correlation. For temperature we used ℓ max = 30: this provides a good estimation because C T T ∼ ℓC T T ℓ ∼ 1/ℓ. Therefore, the contribution from ℓ > 30 is negligible for our purpose. For the E-mode correlation, we used ℓ max = 10. Due to the presence of (ℓ + 2)!/(ℓ -2)! in Equation 7, the contribution from higher multipoles dominates over that from lower multipoles, and the sum has been found to not converge through ℓ max = 1500 (Yoho et al. 2015). We chose ℓ max = 10 based on the cumulative signal-to-noise ratio S/N) as was done in (Chiocchetta et al. 2021): we found a similar plateau for the cumulative S/N at ℓ > 10 of Planck SMICA data, which indicates noise starts to dominate. Therefore, we chose ℓ max = 10 for the E-mode estimator. LiteBIRD has higher S/N and can go to higher ℓ, but we chose the same ℓ max to compare with Planck result. Similar tests based on two-point correlation can be found in, e.g., Copi et al. (2007); Planck Collaboration VII (2020), and these studies provide corresponding checks in the angular space. It is also possible to understand the lack of correlation via the amplitude of low multipoles as the relationship between them was noticed in Copi et al. (2009); Gruppuso ( 2014 We hereafter define the PTE as the percentage of simulations with statistic values -S XX 1/2 in the present case -greater than the value found in the Planck SMICA data (Section 3 and 4.1) or the Special CMB (Section 4.2) to reflect statistical significance for the anomalies. With this prescription, the PTE's for S XX 1/2 in the Planck case (i.e., using Planck simulations) are 97.0% for temperature and 24.0% for E-modes. A high PTE value here means the observed power on large angular scales is smaller than expected from the Planck simulations. The PTE value for the temperature estimator is consis-

T ( = 2) T ( = 3) T ( = 2 and = 3) E ( = 2) E ( = 3) E ( = 2 and = 3) tent with the p-value obtained with Planck public QML spectra in Muir et al. (2018) (p-value ∼ 6%, which corresponds to a PTE value ∼ 94%), but lower than reported in Planck Collaboration VII (2020) (> 99.9%). The inconsistency could be due to the fact that we used a different method to compute the S 1/2 estimator. The PTE for the E-mode is different from that in Chiocchetta et al. (2021) mainly because we used a different dataset. See Section 4 for more comments on this. 3.2. Quadrupole-Octupole Alignment Tegmark et al. (2003) first pointed out that the CMB temperature's quadrupole and octupole are planar and align well with each other using WMAP data. This was later verified in Frommert & Enßlin (2010); Copi et al. (2015a); Muir et al. (2018) and others. We show the quadrupole, octupole, and their superposition of SMICA temperature (left) and E-mode (right) maps in Figure 4. It can be seen clearly that the temperature quadrupole and octupole are both planar and aligned.

Probably the most mathematically complete tools to describe directions of multipoles are the Maxwell multipole vectors (MMVs), and there are several studies that apply MMVs to the CMB quadrupole-octupole (Q-O) alignment (Copi et al. 2004(Copi et al. , 2015a;;Schwarz et al. 2016;Muir et al. 2018). But here we chose a physically more intuitive method to show the alignment of the quadrupole and octupole following de Oliveira-Costa et al. ( 2004) and Dvorkin et al. (2008). The relation of this method to the MMVs is discussed in Copi et al. (2006).

Treating the CMB fluctuations as a wave function (using temperature as an example),

the direction for any multipole is defined to be the direction n of the axis that maximizes the angular momentum dispersion:

where a ℓm ( n) are the spherical harmonic coefficients of the CMB map in a rotated coordinate system with its z axis along the n direction (de Oliveira-Costa et al. 2004). By scaling the angular momentum dispersion as

with X ∈ {T, E}, the value of L 2 ℓ can be used to indicate how 'planar' a multipole is, or mathematically, how large m = ±ℓ are compared with the |m| < ℓ components for multipole ℓ. 8 Finally, the direction n * that maximizes

is treated as the 'joint' direction of quadrupole and octupole. The closer L 2 23 (n * ) is to 1, the more planar and aligned the quadrupole and octupole are, and n * is normal to that plane.

In practice, we estimated a ℓm 's without applying the common masks because otherwise the a ℓm 's can be biased (Copi et al. 2006). To understand the impacts from foreground residuals in temperature data, we computed the L

2(T T ) 23

estimator for all-sky temperature maps from multiple component separation methods: Planck PR3

8 The m = ±ℓ components have only azimuthal variation around n whereas m < ℓ components have polar variation as well.

Commander, NILC, SEVEM and SMICA. We also computed L

2(T T ) 23

for versions of these maps with regions of the strongest Galactic emission inpainted as in Planck Collaboration IV (2020). The directions and PTEs (with respect to synfast + Planck SMICA noise simulations) of the L

2(T T ) 23

for 8 different maps were found to be largely consistent with each other, which suggests that the temperature data are not strongly affected by foreground residuals, unless the foreground is impacting the 8 different maps in a similar way. We still adopted the declination limitation for the CLASS case, and the impact of this limitation can be found in Section 4. Once the alignment direction for a temperature map is settled, we use that direction for calculating L

2(EE) 23

directly to test whether the same axis and planarity is preferred by the E-modes. The alignment estimator relies on the phase information, hence would be more sensitive to issues related to residual foregrounds and systematics. However, given no existing reliable way to debias the Emode quadrupole phase, we chose to use the ℓ = 2 phase from Planck data directly. The E-mode estimator value could be biased and is not well-constrained in presence of Planck noise, which is reflected in the second panel in Figure 8.

The PTE values for these statistics comparing to Planck simulations are 0.3% for the temperature and 84.1% for the E-mode. A low PTE value here means the observed quadrupole and octupole are more planar and aligned in the Planck data than in most of the Planck simulations. The PTE value for the temperature estimator is consistent with the p-value (∼ 0.4%) reported in Muir et al. (2018).

Point-parity asymmetry

One of the earliest tests of whether temperature anisotropies have a preference for antisymmetric parity on large angular scales was carried out by Land & Magueijo (2005), and they did not find a strong preference based on the estimator and data they used. Equipped with a different estimator, Kim & Naselsky (2010) found the significance level was higher when they applied it to WMAP 7-yr temperature anisotropy data. More work has been done to better understand this anomaly since then (e.g., Gruppuso et al. 2011;Kim et al. 2012;Cheng et al. 2016;Aluri et al. 2017;Panda et al. 2021). The corresponding result in Planck Collaboration VII (2020) based on the Planck 2018 data (significance above 3σ for temperature) also indicates that this anomaly deserves investigation.

We quantified this point-parity asymmetry and the corresponding E-mode one using the same estimators as in Planck Collaboration VII (2020)

where C XX ± (ℓ max ) are sums over even (+) or odd (-) multipoles between ℓ min and ℓ max ,

and ℓ ± tot is the total number of even (+) or odd (-) multipoles included in the sum. C XX,± ℓ is the T-or E-mode angular power spectrum at even (+) or odd (-) multipoles. It follows that R T T ℓmax and D EE ℓmax stand respectively for the ratio-of and difference-between the average of band-powers for odd and even multipoles. The less sensitive D EE ℓmax was used for E-mode data because the low signal-to-noise ratio in Planck SMICA data resulted in negative estimation of C EE ℓ , which causes numerical problems as the denominator of the ratio approaches zero. Both estimators depend on the cut-off multipoles ℓ min and ℓ max . The ℓ min was chosen to be 2 for TT and 3 for EE, and the ℓ max was chosen to be 24 for both TT and EE because in Planck Collaboration VII (2020) the minimum lower-tail PTE for TT was found at ℓ = 24. This ℓ-range is motivated by the temperature anomaly and chosen for consistency with past studies. However, we note that, unlike with temperature, the E-mode estimator D EE 24 will be dominated by the reionization bump at the lowest multipoles. Future studies may well consider modifying the estimator, e.g., starting at higher ℓ min or using weights, to better capture even-odd variations up to higher E-mode multipoles. We show power spectra for the first 30 multipoles in Figure 5, and see Table 2 for a summary of masks and noise simulations used for this anomaly.

The PTE values for R T T 24 and D EE 24 , comparing to Planck simulations, are 97.7% for temperature and 83.2% for E-mode. A high PTE value here means the observed even multipoles are smaller than the odd multipoles in Planck data when compared to the Planck simulations. Although cut off at a different multipole, the PTE value for the temperature is consistent with the results for R T T 27 computed with Planck public QML spectra reported in Muir et al. (2018) (p-value ∼ 3%, which corresponds to a PTE value ∼ 97%, with respect to synfast simulations) and slightly less prominent than reported in Planck Collaboration VII (2020) (PTE value ∼ 99%). Additionally, in Muir et al. (2018) R T T 27 was found to be correlated with S T T 1/2 due to their dependence on the temperature quadrupole, and our results show similar 2-dimensional marginalized distributions with the R T T 24 estimator (see discussion in Section 4.1.2). The PTE for the E-mode is less anomalous than reported in Planck Collaboration VII (2020) (minimum p-value ∼ 6% at ℓ max = 27 for HM data), but is sta- tistically consistent given the large uncertainty on D EE 24 introduced by Planck noise (Figure 8).

Hemispherical Power Asymmetry

Using either COBE or WMAP data, Eriksen et al. (2004) and Hansen et al. (2004) first pointed out that the CMB temperature power spectrum is significantly stronger in the southern ecliptic hemisphere than in the northern hemisphere. Instead of using the power spectrum, one can also quantify the asymmetry in terms of variance in the map (e.g., Monteserín et al. 2008;Akrami et al. 2014). We adopted the latter approach and refer to this anomaly as "the hemispherical power asymmetry". While the asymmetry is most visually apparent on the largest angular scales, it even persists at scales below 10 • .

To quantify this anomaly, we estimated the dipole amplitude and direction in local-variance maps (LVMs) (Akrami et al. 2014;Planck Collaboration VII 2020). An LVM is a low-resolution map with pixel values capturing the localized variance information of neighbouring pixels in a higher resolution map. Therefore, the dipole amplitude of an LVM indicates the preference for the localized variances to be anomalously high in one hemisphere. We followed the procedure in Section 7.1 of Planck Collaboration VII (2020) for making LVMs with NSIDE=16 (θ pix = 3.66 • ) resolution from temperature and E-mode maps with NSIDE=64 (θ pix = 0.92 • ) resolution. Each pixel in the LVM map was assigned the variance of temperature or E-mode data within a disk of θ = 4 • radius, centered on the LVM pixel. (Therefore adjacent LVM pixels, being less than 4 • in width, are correlated.) Masks with resolution NSIDE=64, identified in Table 4, were used to exclude temperature or E-mode pixels from the variance calculations. Correspondingly, LVM pixels with more than 90% of the temperature or E-mode data masked in their associated 4 • -radius disks were excluded from further analysis.

ΛCDM and noise simulations, listed in Table 4, were used to compute the average variance and variance-on-the-variance for each pixel in the LVM. We emphasize that these quantities depend both on the ΛCDM model and the noise properties of each experiment. We then fit the monopole and the dipole for each LVM by

where v p is the LVM pixel value at pixel p; vp and σ 2 vp are the variance and variance-on-the-variance from each LVM pixel as determined by the simulations; d 0 captures the monopole component of the LVM; d ≡ (d x , d y , d z ) is the dipole component; and rp is the unit direction vector pointing to the center of pixel p. The sum is over all LVM pixels that were not excluded due to masking as described in the previous paragraph. Subtracting vp from the LVM removes the variance bias arising from ΛCDM fluctuations and the experimental noise. As with a standard Gaussian-likelihood analysis, σ 2 vp down-weights the noisier LVM data. This method was applied to both the temperature and E-mode LVMs. We tested with our simulations that best fit values of d x , d y and d z using this likelihood are unbiased (Section 4.3).

Figure 6 shows the LVMs from Planck SMICA data and corresponding dipole information.

The fitted dipole orientation (Galactic coordinates) and amplitude (d 4 , with 4 representing the radius, 4 • , of the disk on which we computed local variances) from Planck SMICA are (l, b) = (209

4

= 0.008 ± 0.004 µK 2 for E-mode. The errors are the 1σ ranges of the dipole orientations obtained from the LVM simulations with dipole estimated from Planck SMICA inserted. The dipole directions we obtained are consistent with the results in Planck Collaboration VII (2020) for the SMICA component-separated maps, where (209 • , -15 • ) for temperature and (219 • , -16 • ) for E-mode was reported. The PTE values for the dipole amplitudes are 0.1% for temperature and 23.7% for E-mode. A low PTE value here means the variance dipole from Planck is stronger than the expectation from Planck simulations. The temperature PTE we obtained is consistent with the finding in Planck Collaboration VII (2020) (none of the 1000 simulations had a dipole amplitude larger than obtained from data), but the E-mode one is higher than the reported 5.5%. We discuss more about our dipole estimator in Section 4.3.

FORECAST RESULTS

Unconstrained universe

We first checked the distributions of estimators including the effects of the cosmic variance, i.e., the dis- tributions from combinations of 10 4 CMB simulations and 300 noise simulations for each experiment. We also studied the correlations between estimators caused by the combination of cosmic variance and noise.

Constraints on the anomaly estimators

Table 5 shows the 95% confidence intervals derived from simulations for the 8 anomaly statistics log 10 S XX 1/2 , L

2(XX) 23

, R T T 24 , D EE 24 , d XX

4

(XX ∈ {T T, EE}) defined in the previous section. The second column in Table 5 shows the estimators with corresponding units, and the third column gives the values of estimators from Planck 2018 SMICA data, with PTEs computed based on simulations with Planck noise displayed in the brackets. In the last column we show the medians and 95% confidence intervals derived from the Ideal case (cosmicvariance-limited noise level with no mask applied). The three double columns in between show the unitless bias and error of the distributions with Planck, CLASS and LiteBIRD noise relative to those shown in the last column. The error is the 95% confidence interval width of the Planck/CLASS/LiteBIRD simulation normalized by the same width of the Ideal simulation. The bias is normalized in the same way and is computed as the difference between medians of the Planck/CLASS/LiteBIRD and the Ideal case.

For temperature, all the biases are below 0.15, meaning that the results from the three experiments are consistent with each other and the Ideal case. The errors are increased for the CLASS case for the odd-parity and hemispherical power asymmetry anomaly due to the limited sky coverage. Otherwise, no significant change from the current Planck constraints was found for any temperature anomaly estimator mainly because, on these scales, the uncertainty caused by the instrumental noise is smaller than that due to the cosmic variance.

For E-mode polarization, the errors on the estimator for lack of correlation (log S EE 1/2 ) are similar between the Ideal, CLASS, and LiteBIRD experiments. Furthermore, the bias is modest for CLASS and negligible for LiteBIRD. The significant bias for Planck (1.27) implies that the current Planck value is more than 10 times greater than would be expected from the Ideal E-mode measurement. This is because the S EE 1/2 estimated from Planck data is dominated by instrument noise, at least for the SMICA reduction. Instrument noise also explains the higher S EE 1/2 uncertainty from Planck. Similar analysis has been done in Chiocchetta et al. (2021) with the Planck HFI 100 × 143 power spectra, which has less noise contribution than the SMICA-based spectra used here. Table 5 shows no improvement for CLASS and LiteBIRD compared to Planck for the Q-O alignment estimator (L

2(EE) 23

) due to the contribution from cosmic variance, which (being isotropic) trivially saturates the error budget. Given a fixed Q-O alignment in the simulations, LiteBIRD can determine the corresponding Q-O alignment estimator to much higher precision than Planck. This will be discussed more in Section 4.2. The distribution of D EE 24 from Planck SMICA simulations is not significantly biased (0.27) from the Ideal one, but the error will be improved from 2.8 to 1.1 and 1.0 for CLASS and LiteBIRD, respectively. For the hemispherical power asymmetry, both the bias and error shrink significantly from the Planck case to CLASS and LiteBIRD. In particular, the high bias (1.37) and large error (3.6) of Planck SMICA indicate that the current SMICA-based constraints on the E-mode hemispherical power asymmetry are measurement-noise dominated.

Based on the observations listed above, we conclude that the current SMICA-based constraints are noise dominated for the commonly used E-mode anomaly estimators of lack of correlation, odd-parity asymmetry and hemispherical power asymmetry. The situation will improve with results from LiteBIRD-like experiments.

Correlation

Figure 7 shows the marginalized distributions of the 8-dimensional space spanned by the estimators. We use blue to denote contours and histograms from Planck, pink from CLASS, green from LiteBIRD, and black from Ideal simulations. The diagonal panels display the posterior distribution for each estimator from the four types of simulations. In these panels, the blue vertical lines mark the estimator values from Planck SMICA data, the black dashed vertical lines show the E-mode estimators for the Special CMB (Section 4.2), and the numbers on the top-right are the widths of the 95% intervals for the corresponding estimators. The off-diagonal panels show the 2-dimensional marginalized distributions, in which we used black unfilled contours to show distributions from the Ideal case and filled contours for those from the three experiments. The blue crosses mark measurements from SMICA data and black crosses mark the E-mode measurements from the Special CMB. The numbers on the top-right are the Pearson correlation coefficients.

Consistent with the conclusions in Section 4.1.1, the contours and distributions of Planck simulations are significantly biased relative to the Ideal case for log 10 S EE 1/2 and d EE 4 , and the contours for LiteBIRD are similar to the Ideal contours. The most significant correlation (∼ 0.3) is between log 10 S T T 1/2 and R T T 24 . (The CLASS simulations show a smaller correlation as the uncertainty in power spectra estimation was larger due to the limited sky coverage.) This correlation is consistent with the result in Muir et al. (2018) and is mostly due to the correlation between S T T 1/2 and the quadrupole C T T 2 , as a larger S T T 1/2 implies a stronger C T T 2 hence a larger R T T 24 . A similar correlation is not found between log 10 S EE 1/2 and D EE 24 because D EE 24 is largely determined by the octupole

1/2 is not. The second most significant correlation (∼ 0.17) is between R T T 24 and D EE 24 (except for the Planck simulations in which it is diminished by noise). This is mostly due to the correlation between C T T ℓ and C EE ℓ at low multipoles. The correlation coefficients in almost all the other off-diagonal panels are at or below the ∼ 0.1 level for all four cases including the Ideal one 9 , which means the ΛCDM model does not predict strong correlations for those pairs of estimators. Therefore, if an E-mode estimator reinforces the anomaly identified by its temperature counterpart, then this result would not be due to correlation between estimators but would present an independent challenge to the statistical-fluke hypothesis.

Tests based on the Special CMB

In Section 4.1, we computed the distribution of the anomaly estimators resulting from ΛCDM cosmic vari-9 The 0.2 correlation between D EE 24 and log 10 S EE 1/2 for Planck simulations can be blamed to the low signal-to-noise ratio.

1.864 1.809 1.867 1.868 1 0 1 2 log 10 S EE 1/2 +0.015 +0.003 +0.001 +0.001 1.412 1.053 0.989 0.983 0 . 6 0 0 . 7 5 0 . 9 0 L 2(TT) 23 +0.005 +0.059 +0.007 +0.007 +0.010 -0.012 -0.008 -0.008 0.321 0.322 0.321 0.321 0 . 2 5 0 . 5 0 0 . 7 5 L 2(EE) 23 +0.067 +0.088 +0.110 +0.110 -0.008 +0.009 -0.000 -0.003 +0.071 +0.073 +0.098 +0.098 0.564 0.587 0.576 0.575 0 . 8 1 . 2 1 . 6 R TT 24 +0.282 +0.235 +0.296 +0.312 +0.010 +0.015 +0.006 +0.012 +0.002 +0.008 +0.006 +0.005 +0.004 +0.016 +0.016 +0.018 0.632 0.767 0.594 0.564 0 . 0 0 5 0 . 0 0 0 0 . 0 0 5 0 . 0 1 0 D EE 24 +0.013 -0.033 -0.038 -0.037 +0.237 +0.032 +0.029 +0.026 +0.023 +0.008 +0.013 +0.013 -0.006 -0.041 -0.043 -0.044 +0.016 +0.166 +0.168 +0.172 0.032 0.011 0.011 0.011 4 0 8 0 1 2 0 1 6 0 d TT 4 -0.013 -0.018 -0.012 -0.008 -0.004 -0.003 -0.001 +0.007 -0.008 -0.032 -0.008 -0.008 +0.000 -0.007 +0.001 -0.004 -0.005 -0.025 -0.001 +0.007 +0.010 +0.000 -0.010 -0.017 92.19 117.34 92.30 86.78 3 4 5 6 log 10 S TT 1/2 0 . 0 0 4 0 . 0 0 8 0 . 0 1 2 d EE 4 +0.004 +0.011 +0.005 +0.004 1 0 1 2 log 10 S EE 1/2 -0.014 -0.001 +0.011 +0.008 0 . 6 0 0 . 7 5 0 . 9 0 L 2(TT) 23 -0.005 +0.010 +0.009 -0.004 0 . 2 5 0 . 5 0 0 . 7 5 L 2(EE) 23 +0.019 +0.014 +0.001 +0.007 0 . 8 1 . 2 1 . 6 R TT 24 +0.010 +0.002 +0.011 +0.003 0 . 0 0 5 0 . 0 0 0 0 . 0 0 5 0 . 0 1 0 D EE 24 +0.033 +0.008 -0.000 -0.002 4 0 8 0 1 2 0 1 6 0 d TT 4 +0.002 +0.023 +0.018 +0.011 0 . 0 0 4 0 . 0 0 8 0 . 0 1 2 d EE 4 0.011 0.006 0.004 0.003 Planck-like noise level CLASS-like noise level LiteBIRD-like noise level Cosmic-variance-limited ance and noise from Planck, CLASS and LiteBIRD. As reflected in Table 5 and Figure 7, if moving from Planck to LiteBIRD, the constraining power would increase for three of the E-mode anomaly estimators but no obvious change was found for the Q-O alignment estimator distribution. This is because the fluctuations caused by cosmic variance are maximal without additional noise.

In this section, we fix the CMB realization and focus on the bias and spread of the E-mode anomaly estimator distributions caused by noise only. Therefore, we generated the Special CMB such that the E-mode signal had all four anomalies, and we tested how different experiments will constrain the anomaly estimators for this specific realization.

The way we produced the Special CMB was by first generating 10 5 Gaussian CMB simulations based on Planck 2018 best-fit model and computing their anomaly estimator values in the cosmic-variance-limited case. Then the Special CMB was selected by setting thresholds on the PTEs of the E-mode anomaly estimators. The second column in Table 6 lists the PTEs of the CMB realization we picked. This realization has a more than 95% PTE (or less than a 5% PTE) for the lack of correlation and Q-O alignment, and an almost 5% PTE for the hemispherical power asymmetry. The odd-parity estimator is the least aberrant because we are seeking one realization that has all of these E-mode anomalies, which means that we are using the conditional probability instead of the marginal one. The black triangle-shape contours in the D EE 24 -log 10 S EE 1/2 panel of Figure 7 implies that a smaller log 10 S EE 1/2 corresponds to a smaller spread of D EE 24 . Consequently, when a realization lacks correlation, the probability for its D EE 24 to reach the 2σ level of the marginal distribution is extremely small. For the realization that we selected, the PTE of D EE 24 among simulations that have PTEs for log 10 S EE 1/2 above 95% is 98.5%, which is anomalous enough for our purpose. The visualization of the 4 E-mode anomalies for the Special CMB in analogy to Figures 3-6 can be found in Appendix B.

Our main analysis results based on the Special CMB are shown in Figure 8 and Table 6. In Figure 8, filled histograms are the posterior distributions from the Ideal simulations from Section 4.1 (i.e., its spread reflects cosmic variance). We use blue, pink and green unfilled histograms to denote distributions for the Special CMB with Planck, CLASS and LiteBIRD noise simulations, respectively. The black dashed vertical lines mark the estimator value of the Special CMB we picked and blue vertical lines are those from Planck SMICA, staying consistent with Figure 7. Unsurprisingly, we found that the spread of the unfilled histograms reduces as one goes from using Planck noise to LiteBIRD for all Emode anomaly estimators. The biases are different in different cases: in the first panel for log 10 (S EE ), we noticed that all three unfilled histograms are positively biased, which makes sense because noise contributes in S T T,EE 1/2 positively by definition as displayed in Equation 10, and the bias becomes smaller as the noise level becomes smaller. In the second panel, corresponding to the Q-O alignment, the bias for the Planck histogram is mainly due to the noise while, for CLASS, bias results from the incomplete sky coverage. This emphasizes the importance of having both low noise level and full sky coverage for constraining L 2(EE) 23

. Nearly no bias was found for all three cases for the odd-parity estimator in the third panel. In the last panel for the variance dipole amplitude, the behavior of the Planck and CLASS histograms are as expected: the positive bias goes down as noise level decreases. The bias of the LiteBIRD histogram is mainly related to the shape of the mask: we found that using LiteBIRD noise simulations but with full sky coverage only resulted in an -0.005 bias on the dipole amplitude. The unitless bias and error listed in Table 6 are consistent with our observations of Figure 8. Bias and error are still normalized by the error on the Ideal simulations as in Section 4.1.1, but the bias is now computed relative to the Special CMB estimator values. We conclude that, compared to the bias and error of Planck noise, LiteBIRD is promising in providing much stronger (more than an order of magnitude) constraints on the 4 E-mode anomaly estimators, while observation from CLASS is promising in constraining D EE 24 . While the sky coverage of CLASS limits its constraining power relative to the all-sky LiteBIRD measurement, the largescale E-mode data from CLASS will be valuable as a high-S/N check of LiteBIRD.

Constraints on the hemispherical power asymmetry orientation

We did two extra analyses to understand the recovery of the LVM dipole orientation from our method under the impacts of noise and the E-mode mask. To investigate the impacts of the noise, we added the dipole from the LVM of the Special CMB to the LVM of Gaussian CMB+Noise simulations for all three experiments. We then fitted for the dipole of these new LVMs using Equation 20, with vp and σ 2 vp computed for each experiment from the simulations without the Special CMB dipole. The left and middle panels in Figure 9 show 1, 2σ confidence regions of dipole orientations on NSIDE=16 maps for each experiment. In the left panels we used the Galactic coordinate system. To better show the shape of the contours, in the middle panels we ro- ), point-parity asymme-(D EE 24 ) and hemispherical power asymmetry (d EE 4 ). The gray filled histograms represent the spread caused by cosmic variance; other unfilled histograms represent the spread caused by Planck (blue), CLASS (pink) and LiteBIRD (green) noise, respectively. The black dashed vertical lines mark the estimator values of the Special CMB, and blue solid are those from the Planck SMICA measurements. tated the coordinate system so that the inserted dipole is located at the origin. To obtain the 1, 2σ contours, we first binned the recovered dipole orientations from the 10 4 simulations into the pixels of an NSIDE=16 map.

In the resulting spherical histogram of dipole directions, the 1, 2σ regions correspond to the highest-valued pixels whose sum is equal to (or slightly less than) 68%, 95% of the total number of simulations (10 4 ) contributing to the histogram.

The recovered dipole components d x , d y , and d z were unbiased. However, the projection of the elliptical distribution of the recovered dipoles onto the celestial sphere can skew the angular distribution, especially in the case of high noise. Nevertheless, for the LiteBIRD experiment, we found the distribution of dipole orientations from these simulations centers well on the input dipole direction (cyan stars). For CLASS and Planck, the apparent bias of the 2σ region is caused by radially projecting the 3-dimensional ellipsoidal normal distribution corresponding to symmetric distributions for d x , d y , and d z . The solid angle subtended by the 1σ (2σ) region in Planck SMICA is 12180 deg 2 (28415 deg 2 ), whereas it is 6191 deg 2 (18720 deg 2 ) for CLASS (a ∼ 2× reduction) and 3773 deg 2 (11870 deg 2 ) for LiteBIRD (a ∼ 3× reduction).

We then investigated how the dipole fit depends on the shape of the mask by inserting dipoles aligned with the centers of all the NSIDE=4 map pixels. We used 10 4 simulations in each direction, and the amplitude was fixed to be that from the Special CMB. The maps in the right panels of Figure 9 show the value of Ω(1σ), the fraction of the whole sky occupied by the 1σ orientation confidence region, for each direction. Ω(1σ) varies across the map mainly due to the presence of the mask, and we noticed that for LiteBIRD experiment, the pattern aligns well with the shape of the mask: pixels within the mask seem to have smaller Ω(1σ) values. This implies the fitting is more effective if the gradient of the dipole is unmasked. The alignment is less obvious for the CLASS experiment but follows the same general trend with smaller confidence regions within the Galactic mask or outside of the survey footprint. The Planck experiment result has a quite different pattern due to Center panels: Maps contain the same results as in the left column but are rotated such that the inserted dipole is located at the center. Right panels: Maps show the solid angles, Ω(1σ), of the 1σ dipole-orientation confidence region for additional simulations where the dipole from the Special CMB was recentered on each NSIDE=4 pixel. This result suggests that orientation localization is better when the maximum gradient of the dipole (versus its pole) lies outside the Galactic mask. We used different colormaps to highlight the different ranges of the Ω(1σ).

the lower signal-to-noise ratio. Instead of lying in the mask, the best constraints align with the high-noise regions of the Planck data. We found the mean recovered dipole amplitudes (d x , d y , d z ) matched the simulation input to better than 0.0003 times the standard error of the simulation ensemble for all the 192 input directions (NSIDE=4 map pixel centers). This means our estimation for d x , d y and d z are unbiased. The medians of Ω(1σ) are 0.285, 0.151, 0.096 for Planck, CLASS and LiteBIRD respectively, also suggesting a ∼ 2× (∼ 3×) reduction for CLASS (LiteBIRD) comparing to Planck SMICA. Finally, we emphasize that the uncertainty on the orientation depends on the amplitude of the input dipole, and the one we used from the special realization may not reflect the true case. Therefore, these tests are merely demonstrations based on a dipole with amplitude larger than that found in ∼ 95% of ideal ΛCDM realizations.

CONCLUSIONS

In this paper, we explored how CMB polarization data will improve our understanding of CMB anomalies. The four anomalies we studied are: the lack of correlation on large angular scales, the alignment of quadrupole and octupole, the point-parity asymmetry, and the hemispherical power asymmetry. The definitions for estimators of each of the anomalies are in Section 3. We forecast constraints from future experiments on the temperature and polarization estimators and explored the correlation between the estimators of different anomalies. Our main analyses were based on the combination of 10 4 Gaussian CMB simulations and 3 different types of 300 noise simulations based on Planck SMICA, CLASS and LiteBIRD. We also include tests based on a special CMB realization that was selected from Gaussian CMB simulations, of which all 4 E-mode anomaly estimators are at ∼ 2σ significance.

We found that Planck SMICA does not significantly constrain the four anomalies, but future E-mode measurements look promising. Our unitless bias and error in Table 5 show the constraining power on the lack of correlation S EE 1/2 (log 10 S EE 1/2 ) will be improved by factors of ∼ 20 (∼ 1.3) and ∼ 26 (∼ 1.4) moving from Planck to CLASS and LiteBIRD, respectively. The improvements on S EE 1/2 are high compared to the result in Chiocchetta et al. (2021) as they used a reduction of Planck data that has lower noise. Future E-mode studies with Planck data could yield similar gains for the other anomalies. While no analogous improvement in constraining power was found for the Q-O alignment based on the main simulation set, the constrained CMB simulation test (with quadrupole and octupole moments fixed to those of the Special CMB) showed significant improvement due to reduced instrumental noise (∼ 100 for LiteBIRD) as reflected in Table 6. The same test shows that the limited sky fraction makes it difficult for CLASS to constrain the Q-O alignment estimator. A factor-of-3 improvement was found for the point-parity asymmetry for both CLASS and LiteBIRD, while a factor of 2 (3) was found for the amplitude of the hemispherical power asymmetry for CLASS (LiteBIRD). The localization of the variance dipole extracted from the Special CMB improves from ∼ 30% of the sky (1σ) for Planck SMICA to ∼ 15% for CLASS and ∼ 9% for LiteBIRD. The localization depends on both the dipole amplitude and orientation with respect to masks and high-noise regions.

The correlation between different anomalies is negligible for most of the estimators except for two pairs: one is R T T 24 and log 10 S T T 1/2 , which has a Pearson correlation coefficient r ∼ 0.3 in the cosmic-variance-limited result; the other is R T T 24 and D EE 24 , with r ∼ 0.17. We note that the low correlations between this particular set of estimators does not imply that there are no correlations between potential underlying physical models as each anomaly could have several different estimators. The general lack of correlation found between temperature and polarization estimators implies that if estimators from future polarization experiments reproduce anomalies in temperature, then the statistical-fluke hypothesis will be challenged.

Software: HEALPix (Górski et al. 2005); matplotlib (Hunter 2007; Caswell et al. 2019); numpy (van der Walt et al. 2011); APPENDIX A. E-MODE MASK CONSTRUCTION

We constructed the E-mode mask following the method in Appendix A.5 of Planck Collaboration VII (2020).

The reconstruction of E and B-mode maps using masks described in section 2.2 introduces E/B-mixing, so we need to mask out extra regions of the E-mode map to reduce the impacts of compromised modes. To do this, we generated the root mean square (RMS) residual maps from CMB + noise simulations for the three different experiments. The RMS residual map we used was defined to be δ ) where STD means taking the standard deviation per pixel, Ẽ and B are the E-and B-mode maps reconstructed beginning with the masked sky, and E * and B * are those based on the full sky. We constructed residual maps δ E and δ B from all the simulations we generated (10 4 for each experiment). The E-mode masks were created by first requiring the pixel values of the averaged RMS residual map to not exceed 0.2 µK. The thresholds were chosen based on the upper panels of Figure A.1, which shows the maximal RMS residuals for E, B, and the combination of E and B modes. The masks were then smoothed using a Gaussian smoothing function with FWHM=160 arcmin. We set pixels in the smoothed maps with values smaller than 0.9 to be 0 (masked), otherwise 1 (unmasked). After smoothing and setting thresholds, we removed a few islands of isolated non-zero pixels (with radius < 4 • ) that remained in our masks. Due to the mask apodization and removal of small isolated unmasked regions, the triangles which represent the maximum δ EB and sky fraction of the Emode masks are not exactly on the δ EB curves. E ( = 2) E ( = 3) E ( = 2 and = 3) Figure B.2. The quadrupole (left), octupole (middle), and their superposition (right) of the Special CMB E-mode map. The PTE of L 2(EE) 23 estimator is 1.4%. C. CONSTRAINED UNIVERSE We followed the approach in Copi et al. (2013) for generating the constrained CMB realizations and hereby walk through some of the details. According to the ΛCDM theory, the spherical harmonic coefficients a T ℓm , a E ℓm , and a B ℓm of the CMB maps follow a complex normal distribution with zero mean and covariance matrix (Σ ℓ ) Σ ℓ =     C T T ℓ C T E ℓ 0 C T E ℓ C EE ℓ 0 0 0 C BB ℓ     , (C2) where C XY ℓ (XY in {T T, T E, EE, BB}) are from Planck's best-fit ΛCDM model. Its Cholesky decomposition Σ ℓ ≡ L ℓ • L T ℓ gives L ℓ =      C T T ℓ 0 0 C T E ℓ √ C T T ℓ C EE ℓ -(C T E ℓ ) 2 C T T ℓ 0 0 0 C BB ℓ      . (C3) To generate a unconstrained realization of âℓm = (â T ℓm , âE ℓm , âB ℓm ) T that have the desired covariance matrix, one needs to multiply a complex standard normal random vector ζ = (ζ 1 , ζ 2 , ζ 3 ) T with the L ℓ matrix. The constrained CMB realizations may be generated by setting âT ℓm to be what has been measured in the Planck SMICA map (denoted with a T data ℓm ). The constrained CMB âE ℓm = C T E ℓ C T T ℓ a T data ℓm + C EE ℓ -(C T E ℓ ) 2 C T T ℓ ζ 2 , âB ℓm = C BB ℓ ζ 3 .

(C4)

We followed the convention in HEALPix (Górski et al. 2005) by requiring that a X ℓm = (-1) m a X ℓ-m * where m ≥ 0 and X ∈ {E, B}. We generated 10 5 constrained CMB realizations and looked into the distributions of E-mode anomaly estimators after combining these CMB realizations with different noise simulations as has been done in the main text. The results are displayed in Table 7 and Figure C.1. We found that these new results are largely consistent with those for the unconstrained universe (Table 5 and Fig Confidence-curve matrix of polarization anomaly estimators based on the constrained simulation study, with contours and histograms from Planck (blue), CLASS (pink), LiteBIRD (green), and the cosmic-variancelimited Ideal case (black); see Table 5 for a summary of anomaly names and estimators. The numbers on the topright of off-diagonal panels are the Pearson correlation coefficients, and contours show 1σ and 2σ significance levels. The numbers on the top-right of diagonal panels are the 95% interval widths. The blue crosses and vertical lines represent Planck SMICA measurements, and black crosses and dashed vertical lines represent E-mode measurements from the Special CMB (see Section 4.2 for more about the Special CMB). Copi, C. J., Huterer, D., Schwarz, D. J., & Starkman, G. D. 2006, MNRAS, 367, 79, doi: 10.1111/j.1365-2966.2005.09980.x -. 2007, PhRvD, 75, 023507, doi: 10.1103/PhysRevD.75.023507 -. 2009, MNRAS, 399, 295, doi: 10.1111/j.1365-2966.2009.15270.x -. 2013, MNRAS, 434, 3590, doi: 10.1093/mnras/stt1287 -. 2015a, MNRAS, 449, 3458, doi: 10.1093/mnras/stv501 -. 2015b, MNRAS, 451, 2978, doi: 10.1093/mnras/stv1143 Copi, C. J., Huterer, D., & Starkman, G. D. 2004, PhRvD,