Unexpected Topology of the Temperature Fluctuations in the Cosmic Microwave Background

TL;DR

This study investigates the topology of temperature fluctuations in the CMB using Betti numbers, revealing deviations from Gaussian simulations and suggesting potential implications for primordial non-Gaussianity or non-trivial topology models.

Contribution

It introduces the use of Betti numbers and relative homology to analyze CMB topology, highlighting significant deviations from Gaussian LCDM simulations at specific scales.

Findings

Significant deviations between observed and simulated Betti numbers at 2-7 degrees.

Detection of mildly unusual Euler characteristic behavior at 3.66 degrees.

Strong differences in non-parametric tests across most scales.

Abstract

We study the topology generated by the temperature fluctuations of the Cosmic Microwave Background (CMB) radiation, as quantified by the number of components and holes, formally given by the Betti numbers, in the growing excursion sets. We compare CMB maps observed by the Planck satellite with a thousand simulated maps generated according to the LCDM paradigm with Gaussian distributed fluctuations. The survey of the CMB over $\mathbb{S}^2$ is incomplete due to obfuscation effects by bright point sources and other extended foreground objects like our own galaxy. To deal with such situations, where analysis in the presence of "masks" is of importance, we introduce the concept of relative homology. The parametric $\chi^2$-test shows differences between observations and simulations, yielding $p$-values at per-cent to less than per-mil levels roughly between 2 to 7 degrees. The highest observed deviation for $b_0$ and $b_1$ is approximately between $3\sigma$-4$\sigma$ at scales of 3 to 7 degrees. There are reports of mildly unusual behaviour of the Euler characteristic at 3.66 degrees in the literature, computed from independent measurements of the CMB temperature fluctuations by Planck's predecessor WMAP satellite. The mildly anomalous behaviour of Euler characteristic is related to the strongly anomalous behaviour of components and holes. These are also the scales at which the observed maps exhibit low variance compared to the simulations. Non-parametric tests show even stronger differences at almost all scales. Regardless, beyond the trivial possibility that this may still be a manifestation of an extreme Gaussian case, these observations, along with the super-horizon scales involved, may motivate to look at primordial non-Gaussianity. Alternative scenarios worth exploring may be models with non-trivial topology.

Figures

Tables

Table 1. Table 1: Percentage of sky area covered by the unmasked regions for the various degraded resolutions between N = 1024 𝑁 1024 N=1024 and 16 16 16 , for mask binarization threshold 0.9 0.9 0.9 .

Resolution	% unmasked
1024	77.19
512	76.52
256	75.50
128	73.37
64	72.41
32	69.39
16	66.24

Table 2. Table 2: Table displaying the two-tailed p 𝑝 p -values for relative homology obtained from parametric (Mahalanobis distance) and non-parametric (Tukey depth) tests, for four mask binarization thresholds. The last entry for each threshold is the summary statistic computed across all resolutions. Marked in boldface are p 𝑝 p -values 0.05 0.05 0.05 or smaller.

Relative homology
	Mahalanobis			Tukey Depth
resolution	$b_0$	$b_1$	${\mathrm{𝖤𝖢}}_{\mathrm{rel}}$	$b_0$	$b_1$	${\mathrm{𝖤𝖢}}_{\mathrm{rel}}$
threshold = 0.70
1024	0.236	0.244	0.472			0.302
512	0.492	0.325	0.666	0.134		0.685
256	0.602	0.513	0.760	0.201	0.211	0.579
128	0.260	0.441	0.627		0.171	0.327
64	0.335	0.278	0.528	0.171
32	0.082	0.302	0.442
16	0.018	0.030	0.120
8	0.373		0.012	0.430
summary	0.002	0.001	0.002
threshold = 0.80
1024	0.225	0.278	0.472			0.410
512	0.499	0.348	0.686	0.130		0.649
256	0.559	0.481	0.679	0.139	0.218	0.538
128	0.295	0.484	0.633			0.331
64	0.250	0.269	0.382
32	0.082	0.406	0.538
16	0.024	0.043	0.082
8	0.202		0.013	0.142		0.220
summary	0.001	0.001	0.002		0.032
threshold = 0.90
1024	0.225	0.278	0.472			0.410
512	0.526	0.340	0.661	0.264		0.641
256	0.601	0.584	0.738	0.150	0.306	0.589
128	0.259	0.525	0.486		0.160
64	0.611	0.248	0.444	0.547		0.428
32	0.053	0.305	0.300		0.243	0.340
16	0.009	0.054	0.081
8	0.408		0.007	0.401
summary	0.010		0.001		0.002
threshold = 0.95
1024	0.225	0.278	0.472			0.410
512	0.531	0.340	0.654	0.207		0.445
256	0.602	0.550	0.702	0.139	0.345	0.434
128	0.309	0.524	0.554		0.165	0.446
64	0.420	0.157	0.430	0.276
32	0.076	0.246	0.174
16	0.436	0.016	0.027	0.383
8	0.871	0.047	0.188	0.846		0.517
summary	0.213	0.003	0.010		0.009	0.001

Table 3. Table 3: Table displaying the two-tailed p 𝑝 p -values for absolute homology obtained from parametric (Mahalanobis distance) and non-parametric (Tukey depth) tests, for four mask binarization thresholds. The last entry for each threshold is the summary statistic computed across all resolutions. Marked in boldface are p 𝑝 p -values 0.05 0.05 0.05 or smaller.

Absolute homology
	$Mahalanobis$			Tukey Depth
resolution	${\beta }_0$	${\beta }_1$	$\mathrm{𝖤𝖢}$	${\beta }_0$	${\beta }_1$	$\mathrm{𝖤𝖢}$
threshold = 0.70
1024	0.416	0.318	0.629			0.501
512	0.492	0.371	0.573	0.241	0.133	0.437
256	0.535	0.308	0.488	0.131
128	0.490	0.316	0.542	0.163		0.445
64	0.413	0.184	0.396	0.180
32	0.012	0.024	0.035
16	0.220	0.017	0.135	0.164
8		0.584	0.004		0.647
summary	0.064	0.276	0.040	0.152	0.502
threshold = 0.80
1024	0.416	0.318	0.629			0.501
512	0.492	0.371	0.573	0.241	0.133	0.437
256	0.535	0.308	0.488	0.131
128	0.490	0.316	0.542	0.163		0.445
64	0.413	0.184	0.396	0.180
32	0.012	0.024	0.035
16	0.220	0.017	0.135	0.164
8	0.001	0.584	0.004		0.647
summary		0.204	0.001		0.010	0.001
threshold = 0.90
1024	0.281	0.339	0.642			0.301
512	0.523	0.389	0.697	0.135		0.520
256	0.402	0.329	0.422			0.307
128	0.542	0.216	0.342	0.153		0.298
64	0.276	0.228	0.421			0.448
32	0.045	0.172	0.231		0.175
16	0.058	0.019	0.037
8	0.223	0.475	0.071	0.199	0.510
summary	0.009	0.386	0.024	0.371	0.631	0.206
threshold = 0.95
1024	0.281	0.339	0.642			0.301
512	0.523	0.407	0.698	0.250		0.512
256	0.380	0.339	0.302	0.132
128	0.456	0.226	0.248
64	0.491	0.449	0.794	0.235	0.346	0.431
32	0.261	0.262	0.553		0.308
16	0.137	0.066	0.050	0.155
8	0.155	0.024	0.008
summary	0.049	0.033	0.036

Equations

$$0$$

$$\to {\mathsf{H}}_1(E,M)\to {\mathsf{H}}_0(M)\to {\mathsf{H}}_0(E)\to {\mathsf{H}}_0(E,M)\to 0;$$

$$f_{final}$$

$$a_{lm}^{out}=\frac{b_l^{out}{p}_l^{out}}{b_l^{in}{p}_l^{in}}{a}_{lm}^{in},$$

$$b_0$$

$$b_1$$

$$+\#\{[b,d)\in Dg{m}_1(\mathbb{E}\cup \mathbb{M})\mid +\infty >b\ge \nu >d\};$$

$$b_2$$

$$+\#\{[b,d)\in Dg{m}_2(\mathbb{E}\cup \mathbb{M})\mid +\infty >b\ge \nu >d\}.$$

$$d_{\mathrm{Mahal}}^2(\mathbf{y})=(\mathbf{y}-\bar{\mathbf{x}}{)}^{\mathtt{T}}{\mathbf{S}}^{-1}(\mathbf{y}-\bar{\mathbf{x}}).$$

$$p_{\mathrm{Mahal}}(\mathbf{y})=P[{\chi }_m^2>d_{\mathrm{Mahal}}^2(\mathbf{y})].$$

$$p_{\mathrm{dep}}(\mathbf{y})=\#\{j\mid d_j>d_{dep}(y)\}/n$$

Full text

1 Introduction

Cosmological background. The Lambda Cold Dark Matter (or LCDM) standard paradigm of cosmology postulates that the Universe consists primarily of cold non-relativistic dark matter, which reveals its presence only through gravitational interactions, and the Universe is currently driven by dark energy, causing accelerated volume expansion in this model. The Cosmic Microwave Background (or CMB) radiation, which originates at the epoch of recombination, is the most important observational probe into the validity of the standard paradigm today (Jones, 2017). It is the earliest visible light and offers a glimpse into the processes during the nascent stage of the Universe. Fluctuations about the mean in the temperature field of the CMB correspond to the fluctuations in the distribution of matter in the early Universe, thus understanding the CMB is crucial to understanding the primordial Universe.

The LCDM paradigm together with the inflationary theories in their simplest forms, predict the primordial perturbations to be realizations of a homogeneous and isotropic Gaussian random field (Guth & Pi, 1982). This hypothesis is supported experimentally by CMB observations (Smoot et al., 1992; Bennett et al., 2003; Spergel et al., 2007; Komatsu et al., 2011; Planck Collaboration et al., 2015) and theoretically by the Central Limit Theorem. While it has largely been agreed upon that the CMB exhibits characteristics of a homogeneous and isotropic Gaussian field, there are lingering doubts. The pioneering works of Eriksen et al. (2004a) and Park (2004) challenge the assumption of homogeneity, and the alignment of low multipoles (Copi et al., 2015) challenges the assumption of isotropy. Other noted anomalies include the vanishing correlation function at large scales, and the unusually low variance at approximately $3$ degrees; see Schwarz et al. (2016) for a review and possible interpretations. Planck Collaboration et al. (2014b) independently confirms these anomalies.

The primordial non-Gaussianity remains a topic of ongoing debate. Deviations from Gaussianity, if found, will point to new physics driving the Universe in its nascent stages. The consensus is biased towards its absence (Komatsu et al., 2011; Planck Collaboration et al., 2014b; Matsubara, 2010; Bartolo et al., 2010; Planck Collaboration et al., 2014c); see also (Buchert et al., 2017) for a review and a model-independent route of analysis. Despite mildly unusual behaviour of the Euler characteristic, pointed out in Eriksen et al. (2004b) and Park (2004), the methods employed until today have not provided compelling evidence of non-Gaussianity in the CMB. In contrast, the topological methods of this paper find the observed CMB maps (Planck Collaboration et al., 2016a) to be significantly different from the Full Focal Plane 8 (FFP8) simulations (Planck Collaboration et al., 2014b, 2016b) that assume the initial perturbations to be Gaussian.

Topological methods. Topology is the branch of mathematics concerned with properties of shapes and spaces preserved under continuous deformations, such as stretching and bending, but not tearing and gluing. It is related to, but different from geometry, which measures size and shape. Both geometry and topology have been used in the past to study the structure of the CMB radiation and other cosmic fields. Historically, the predominant tools in this endeavour were the Minkowski functionals, which for a $2$-manifold embedded in the $3$-dimensional space are related to the enclosed volume, the area, the total mean curvature, and the total Gaussian curvature. By the Gauss–Bonnet Theorem, for $2$-manifolds, the latter is $2\pi$ times the Euler characteristic (Euler, 1758), thus providing a bridge between geometry and topology. Early topological studies of the cosmic mass distribution were based on the Euler characteristic of the iso-density surfaces, which generically are $2$-manifolds (Doroshkevich, 1970; Bardeen et al., 1986; Gott et al., 1986; Park et al., 2013). The full set of Minkowski functionals was later introduced to cosmology in (Mecke et al., 1994; Schmalzing & Buchert, 1997; Schmalzing & Gorski, 1998). For Gaussian, and Gaussian-related random fields, the expected values of the Minkowski functionals of excursion sets have known analytic expressions (Adler, 1981; Adler & Taylor, 2010), which is one of the main reasons they have played a key role in the study of real valued fields arising in cosmology and other disciplines.

While the Minkowski functionals have been instructive, the topological information contained in them is limited and convolved with geometric information. Moreover, they are not equipped to address the hierarchical aspects of the matter distribution directly, although partial Minkowski functionals (Schmalzing et al., 1999) may be useful in certain settings. We, therefore, analyse CMB fluctuations in terms of the purely topological concepts of homology (Munkres, 1984), as quantified by Betti numbers (Betti, 1871) and persistence (Edelsbrunner et al., 2002; Edelsbrunner & Harer, 2010). By the Euler-Poincaré formula, the Euler characteristic is the alternating sum of the Betti numbers, implying that the latter provide a finer description of topology (Munkres, 1984). A broad exposition of these concepts, in a cosmological setting, is given in (Pranav et al., 2017; Pranav, 2015; van de Weygaert et al., 2011); also see (Park et al., 2013; Sousbie, 2011; Shivashankar et al., 2016; Adler et al., 2017; Makarenko et al., 2018; Cole & Shiu, 2018) for some applications. Related but slightly different methodologies used for the analysis of cosmological datasets, emanating from concepts in Morse theory maybe found in (Colombi et al., 2000; Novikov et al., 2006; Sousbie et al., 2008).

Results. Our main result is an anomaly of the observed CMB radiation when compared with simulations based on Gaussian prescriptions. The $\chi^{2}$-test yields a significant difference between the number of components and holes in the observed sky compared to the simulations, with $p$-values at per-cent to less than per-mil levels at scales of roughly 2 to 7 degrees. The differences peak sporadically at more than $3\sigma$ at these scales. Non-parametric tests see an even more significant difference between the observation and the simulations at almost all scales. The $\chi^{2}$-test shows the anomaly at roughly the same scales at which the power spectrum exhibits a dip. Eriksen et al. (2004b) reports a mildly unusual Euler characteristic at approximately 3 degrees in the earlier measurements of the CMB radiation by the Wilkinson Microwave Anisotropy Probe (WMAP) satellite, which is related to the anomalous behaviour of components and holes. The noted anomaly motivates a closer look at the standard paradigm. Possible scenarios include, but are not limited to, primordial non-Gaussianity, as well as models with non-trivial topology (Aurich & Steiner, 2001; Aurich et al., 2007; Bernui et al., 2018).

Overview. The workflow in this paper is straightforward: topological descriptors are computed from cosmology data, and statistical tests based on these descriptors are used to compare the observations with simulations. Seection 2 gives a summary of the topological concepts. Section 3 describes the data, the computational pipeline, and a brief account of the statistical tests employed. Section 4 presents the main results of the paper, followed by a summary and conclusions in Section 5.

2 Topological Background

Since the CMB radiation is observed as a scalar field on the $2$-dimensional sphere, the topological concepts needed in this paper are elementary, namely the components and the holes of subsets of this sphere. To count them in the presence of regions with unreliable data, we compute the ranks of the homology groups relative to the mask that covers these regions.

2.1 Excursion sets and absolute homology

Writing ${{\mathbb{S}}}^{2}$ for the $2$-dimensional sphere and $f\colon{{\mathbb{S}}}^{2}\to{{\mathbb{R}}}$ for the temperature field of the CMB, the excursion set at a temperature $\nu$ is the subset of the sphere in which the temperature is $\nu$ or larger: ${{\mathbb{E}}}(\nu)=\{x\in{{\mathbb{S}}}^{2}\mid f(x)\geq\nu\}$. It is a closed set, and we write ${{\beta}_{0}}(\nu)$ for its number of components. A hole is a component of the complement, ${{\mathbb{S}}}^{2}\setminus{{\mathbb{E}}}(\nu)$. Assuming there is at least one hole, we write ${{\beta}_{1}}(\nu)+1$ for the number of holes, and we set ${{\beta}_{2}}(\nu)=0$, because ${{\mathbb{E}}}(\nu)$ does not cover the entire sphere. On the other hand, if there is no hole, we set ${{\beta}_{1}}(\nu)=0$ and ${{\beta}_{2}}(\nu)=1$; see the left panel of Figure 1 for an illustration. These definitions are motivated by the more general theory (Munkres, 1984) in which the $p$-th Betti number is the rank of the $p$-th homology group: ${{\beta}_{p}}={{\rm rank\,}{{{\sf H}_{p}}}}$ for $p=0,1,2$. These are the basic objects of homology. The Euler characteristic of the excursion set is the alternating sum of the Betti numbers: ${\sf EC}(\nu)={{\beta}_{0}}(\nu)-{{\beta}_{1}}(\nu)+{{\beta}_{2}}(\nu)$.

The Euler characteristic has a long history in the CMB literature, largely due to the fact that a simple analytic formula for its expected value is known when the CMB is modeled as a Gaussian random field (Adler, 1981; Adler & Taylor, 2010). While such formulas are not known for the Betti numbers, the information they carry is richer. Pranav et al. (2018) present a numerical study of the Betti numbers of Gaussian random fields, and compare them to the Euler characteristic and Minkowski functionals, and find that Betti numbers present a more detailed account of the topological properties of the field compared to the Euler characteristic (also see Park et al. (2013)). In general, near the mean level of $\nu$, one expects the components and holes of the excursion set to be of similar size and number. Accordingly, one expects ${{\beta}_{0}}(\nu)$ and ${{\beta}_{1}}(\nu)$ to be of similar magnitude, combining to give an Euler characteristic close to zero. Such an Euler characteristic tells us nothing about the individual Betti numbers, beyond the fact that they are similar.

2.2 Masks and relative homology

We define the mask to be the region in which the data is not reliable, and denote it by ${{\mathbb{M}}}\subseteq{{\mathbb{S}}}^{2}$. In our application, it includes a belt around the equator corresponding to the thickened disk of the Milky Way, along with other galactic and extra-galactic bright foreground objects that interfere with the observation of the CMB radiation. In an effort to exclude the mask from our computations, we consider the reduced excursion set: ${{\mathbb{E}}}(\nu)\setminus{{\mathbb{M}}}$. Treating ${{\mathbb{M}}}$ as a closed set, this difference is not necessarily closed. An appropriate topological measure is the relative homology of a pair of closed spaces, $(E,M)$, with the second being contained in the first. In our setting, the pair is $E={{\mathbb{E}}}(\nu)$ and $M={{\mathbb{M}}}\cap{{\mathbb{E}}}(\nu)$. Just as in the absolute case, we get relative homology groups in dimensions [math], $1$, and $2$, and we use their ranks for quantification. It is tempting to refer to these ranks as relative Betti numbers, but this is not the traditional terminology, and we just write ${{b}_{p}}={{\rm rank\,}{{{\sf H}_{p}}(E,M)}}$ for $p=0,1,2$. If $M=\emptyset$, then ${{b}_{p}}={{\beta}_{p}}$, for all three choices of $p$, but if the mask overlaps with the excursion set, then there are differences. We explain some of these differences with reference to Figure 1: If ${{\mathbb{M}}}$ overlaps a component of $E$, this component is no longer counted because every vertex in it bounds a path connecting it to the mask. If ${{\mathbb{M}}}$ overlaps the component in two disconnected pieces, we count a new loop, namely the path connecting these two pieces. If ${{\mathbb{M}}}$ covers part of a hole, this hole is still counted because the part of its boundary curve outside the mask is open, with endpoints in ${{\mathbb{M}}}$. If a hole of ${{\mathbb{M}}}$ is contained in the excursion set, we get a surface without boundary. The relation between absolute and relative homology is compactly expressed by the exact sequence of the pair $M\subseteq E$ (Munkres, 1984):

[TABLE]

Without going into details, this means that we can assign non-negative integers to the arrows so that the rank of each group is the sum of integers assigned to its incoming and outgoing arrows. For example in Figure 1, we have $0\to 0\to 0\to 1\to 1\to 1\to 2\to 4\to 3\to 0\to 0$, and it is easy to find the assignment of integers that satisfies the stated property.

2.3 Variationally maximal loops

When we count ${{\beta}_{1}}+1$ holes in absolute homology, we really count ${{\beta}_{1}}$ loops needed to separate them. In relative homology, the connection is not as intuitive because we also have open loops, whose endpoints lie in the mask; see Figure 2. Generally, there are uncountably many ways to draw a loop, and in homology they are all considered equivalent. The set of equivalent loops is called a homology class, and any one of the loops in the class is a representative. These classes are the elements of the $1$-dimensional homology group, which is a vector space. The rank of this group counts the classes that are needed to span the vector space.

For visualization, it is desirable to have a unique representative for each class. Similar to the intuitive notion of the rim of a crater, we choose this representative as high as possible, alternating between peaks and saddles of $f$ which it connects via ridges within the reduced excursion set. We refer to this loop as the variationally maximal representative of its class; see Figure 2 for an example. While constructing variational maxima for smooth scalar fields may be problematic, the persistence algorithm applied to a piecewise linear scalar field produces them as a side effect of reducing the boundary matrix; see 3.2.

3 Data and methods

3.1 Data

Decades after the accidental discovery of the CMB, its first space based observational probe was carried out by the Cosmic Microwave Background Explorer (COBE) satellite (Smoot et al., 1992), establishing that the CMB is perfect black-body radiation. Later, the Wilkinson Microwave Anisotropy Probe was launched to study the temperature fluctuations in greater detail (Spergel et al., 2007). Most recently, the high precision Planck mission was launched, measuring temperature fluctuations to an accuracy of $10^{-5}$ degrees (Planck Collaboration et al., 2014a), and at a resolution of five arc-minutes, giving the most detailed and precise measurement of CMB temperature fluctuations currently available. We use the Planck maps for our analyses (Planck Collaboration et al., 2016a, b).

Format. The CMB sky maps are presented in the HealPix format (Górski et al., 2005), which is an equal-area pixelisation of the sphere, which we denote as ${{\mathbb{S}}}^{2}$; see Figure  3. Using the faces of the rhombic dodecahedron, we start by decomposing the sphere into twelve patches. Fine resolution is achieved by dividing these patches into $N^{2}$ equal area pixels each, so that the total number of pixels at this resolution is $12\times N^{2}$. At maximum resolution $N=2048$, the maps have about 50 million pixels.

Observed sky. The Planck satellite observes the sky at seven different frequency bands, leading to component-separated maps using four different techniques: Commander-Ruler (C-R), NILC, SEVEM, and SMICA; cf.  Planck Collaboration et al. (2016a). We use these component-separated CMB maps throughout. These maps are contaminated by noise from various sources, including inherent detector noise, and efforts by the Planck team to denoise the data have not been completely successful. Consequently, our analysis is performed on the maps produced by combining the CMB and noise maps for each realization of the simulation. This is a fairly simple task, given that the map making exercise is linear in nature:

[TABLE]

Simulations. In addition to the observed data, the Planck team released a set of Full Focal Plane 8 (or FFP8) simulations (Planck Collaboration et al., 2016b) of both the CMB and noise. We use 1000 NILC simulations for our computational experiments. These simulations assume that the CMB is a Gaussian random field, consistent with the null-hypothesis of Gaussianity for the CMB, which is what we wish to check, and we use them to estimate the error-bars for testing the significance of differences between observed and simulated maps.

Degradation. In order to perform a scale-dependent analysis of the CMB maps, we degrade them to resolutions between $N=1024$ and $8$, dividing $N$ by $2$ from one level to the next. The process of degradation amounts to decomposing them into spherical harmonics on the full sky at the input resolution. The spherical harmonics coefficients $a_{lm}$ are then convolved to the new resolution using the formula (Planck Collaboration et al., 2014b)

[TABLE]

where $b_{l}$ is the beam transfer function, $p_{l}$ is the pixel window function, and in and out denote the input and the output functions at the different resolutions respectively. They are then synthesized into maps at the output resolution directly.

Masks. The observation of the CMB by the Planck satellite is incomplete in some regions of the sky, typically as a result of interference from bright foreground objects, such as our own galactic disk and bright point sources. In these regions, the CMB sky map is reconstructed as a constrained Gaussian field. In order to avoid these areas in the analysis, we use the most conservative UT78 mask released by the Planck team, see Figure 4. It is a combination of all foreground objects with least sky coverage and so provides for a conservative analysis. The mask is a binary map, where reliable pixels of the CMB map are marked by the value 1, and the unreliable parts by 0.

For the scale dependent analysis, we also degrade the masks, so that the map and the mask have the same resolution. Degrading the original binary UT78 mask converts it into a non-binary mask in a thickened zone at the boundary separating the reliable part of the mask from the non reliable part. Figure 5 presents these yet to be binarized masks. To re-convert it to a binary mask, we set a range of binarization thresholds for our experiments: 0.7, 0.8, 0.9, 0.95. Pixels with values above or equal to the binarization threshold are marked as 1, and the rest as 0. Figure 6 presents the binarized maps at various degraded resolution, for binarization threshold 0.9. Table 1 presents the percentage of sky that is useable for analysis after masking at various resolutions, for this threshold. The percentage of usable area drops with decreasing resolution, with only $66\%$ for $N=16$. Figure 7 presents a visualization of the degraded and masked maps for all the resolutions analysed in this paper in the Mollweide projection view.

3.2 Computational Pipeline

The computational pipeline is tailored specifically to the Planck data. The preprocessing step involves converting the CMB maps given in absolute units to a dimensionless unit, corrected for mean and scaled by the standard deviation (computed using the non-masked pixels only). We use the HealPix package for the preprocessing step. The output of this operation is the normalized temperature values on $12{N}^{2}$ pixels, along with their coordinates on the sphere, which is the input to subsequent steps, which we discuss in five subsections: (i) triangulating the surface of the sphere with the pixel centers as vertices, (ii) sorting the vertices, edges, triangles to form an upper-star filteration, (iii) computing the persistence in terms of a reduced boundary matrix, (iv) computing the ranks of the relative homology groups ${{b}_{p}},p=0,1,2$ and (v) the variationally optimal loops from the reduced matrix. The software is written in C++ and designed to handle very large data sets.

Triangulation. The HealPix format stores the data in twelve square arrays of ${N}^{2}$ pixels each, with ${N}=2048$ at the finest resolution. The centers of these pixels are points on the faces of a rhombic dodecahedron. With central projection, these points are mapped to the $2$-sphere and distorted to achieve an approximately equal-area decomposition. Taking the convex hull of these points in ${{\mathbb{R}}}^{3}$, we get a convex polytope whose boundary is homeomorphic to the sphere. Most of the faces will be triangles, and the occasional faces with $k\geq 4$ sides can be subdivided into $k-2$ triangles to get a triangulation of the sphere. This triangulation is the input to all the down-stream computations. Consisting of $V=12{N}^{2}$ vertices, $3V-6$ edges, and $2V-4$ triangles, it represents the temperature field, $f\colon{{\mathbb{S}}}^{2}\to{{\mathbb{R}}}$, by storing the temperature value at every vertex. We implicitly assume a piecewise linear interpolation along the edges and the triangles. Figure 8 illustrates such a triangulation, using colours to visualize the temperature field. We use the CGAL   library (The CGAL Project, 2018) to implement the triangulation.

Upper-star filtration. Given a triangulation $K$ of ${{\mathbb{S}}}^{2}$, let $K(\nu)\subseteq K$ contain all simplices (vertices, edges, triangles) whose temperature values are $\nu$ or larger. We use $K(\nu)$ as a proxy for ${{\mathbb{E}}}(\nu)$, the corresponding excursion set. Indeed, because of the linear interpolation, there is a deformation retraction from ${{\mathbb{E}}}(\nu)$ to $K(\nu)$ (Edelsbrunner & Harer, 2010), which implies that the two have corresponding components and holes. To process the sequence of excursion sets, it thus makes sense to sort the vertices of $K$ in the order of decreasing temperature value. More precisely, we order the simplices of $K$ such that ${\sigma}$ precedes ${\tau}$ if (i) $f({\sigma})>f({\tau})$ or (ii) $f({\sigma})=f({\tau})$ and ${\rm dim\,}{{\sigma}}<{\rm dim\,}{{\tau}}$, in which $f({\sigma})$ is the minimum temperature value of the one, two, or three vertices of ${\sigma}$. The remaining ties are broken arbitrarily. Assuming any two vertices have different temperature values, then the edges and triangles that immediately follow a vertex are exactly the ones in the upper star of that vertex. We therefore call any ordering that satisfies (i) and (ii) an upper-star filter of $K$ and $f$. The corresponding upper-star filtration consists of all prefixes of the filter, each representing an excursion set. This filtration is instrumental in computing the persistence of components and holes.

Computing persistence. Given an upper-star filter of the piecewise linear temperature field, there is optimized software available to compute its persistence (Bauer et al., 2014). We base our persistence computation on an adaptation of the software. This software is a sophisticated implementation of the basic algorithm, which we now describe and modify to get the variationally optimal loops. We write ${\sigma}_{1},{\sigma}_{2},\ldots,{\sigma}_{n}$ for the simplices in the triangulation of the sphere, sorted into an upper-star filter. Let $\partial[1..n,1..n]$ be the corresponding ordered boundary matrix, with $\partial[i,j]=1$, if ${\sigma}_{i}$ is a face of ${\sigma}_{j}$ and ${\rm dim\,}{{\sigma}_{i}}={\rm dim\,}{{\sigma}_{j}}-1$, and $\partial[i,j]=0$, otherwise. This matrix is sparse and stored as such. The standard persistence algorithm reduces the matrix from left to right. To reduce column $j$, we subtract columns to the left of $j$ with the goal to move the lowest $1$ in column $j$ higher or eliminate it altogether. We use modulo $2$ arithmetic, so subtracting is the same as adding: $1-1=1+1=0$. We call column $j$ reduced if it is zero or its lowest $1$ has only [math]s in the same row to its left. We modify the standard algorithm by continuing the reduction even if the lowest $1$ cannot be changed any more, calling the final result fully reduced. To be unambiguous, we explain this algorithm in pseudo-code, where we write ${\rm pivot}{({j})}$ for the row index of the lowest $1$ in column $j$.

[TABLE]

The running time of this algorithm is cubic in the number of simplices in the worst case, but the available optimized software is typically much faster.

Ranks of relative homology groups. For computing the ranks of the homology groups relative to the mask, we set the vertices belonging to the mask at $+\infty$, and consider the complex ${{\mathbb{M}}}$, induced by the union of these vertices. This mask is closed by definition. We then compute the filtration and persistence diagram corresponding to absolute homology of ${{\mathbb{E}}}\cup{{\mathbb{M}}}$. Writing $Dgm_{p}({{\mathbb{E}}}\cup{{\mathbb{M}}})$ for the $p$-dimensional persistence diagram, and recalling that each diagram consists of intervals with real birth- and death-value, $b>d$, we get the ranks of homology groups relative to the mask:

[TABLE]

For computing absolute homology, we set the mask pixels at $-\infty$ and consider the union of such vertices as the mask, which is open by definition.

3.3 Statistical Tests

The data consists of topological summaries (${{b}_{0}},{{b}_{1}},{\sf EC_{\rm rel}}$) obtained from $1000$ simulations, as well as of the observed CMB field, processed according to NILC scheme. The goal is to estimate the probability that the physical model that produced the simulations would produce quantities consistent with those from the observed CMB field. Let $\mathbf{x}_{i}\in\mathbb{R}^{m}$, $i=1,\ldots,n$, be a sample of i.i.d. $m$-dimensional vectors, drawn from a distribution $F$. Let $\mathbf{y}\in\mathbb{R}^{m}$ be another sample point, assumed to be drawn from a distribution $G$. We wish to test the (null) hypothesis that $F=G$, and shall give the test results in terms of $p$-values, that compute the probability that $\mathbf{y}$ is ‘consistent’ with this hypothesis. We consider two methods of testing for statistical consistency. The first is a parametric test based on the Mahalanobis distance (Mahalanobis, 1936), also known as the $\chi^{2}$-test. The second is a non-parametric test based on the Tukey depth (Tukey, 1975). The $\chi^{2}$-test is more standard but has the disadvantage of assuming that the compared quantities follow a Gaussian distribution, while the Tukey depth works without any assumption on the distribution.

Mahalanobis distance or $\chi^{2}$-test. Let $\bar{\mathbf{x}}=\sum_{i=1}^{n}\mathbf{x}_{i}/n$ and $\mathbf{S}=\sum_{i=1}^{n}(\mathbf{x}_{i}-\bar{\mathbf{x}})(\mathbf{x}_{i}-\bar{\mathbf{x}})^{\tt T}/(n-1)$ the sample mean and covariance matrix of the sample $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$. The squared Mahalanobis distance of $\mathbf{y}$ to $\bar{\mathbf{x}}$ is then

[TABLE]

If $F$ is assumed to be Gaussian and $n$ is large, then under the hypothesis that $G=F$ the squared Mahalanobis distance (5) is approximately distributed as a $\chi^{2}$ distribution with $m$ degrees of freedom. Thus the corresponding $p$-value is

[TABLE]

Tukey depth. As will be shown in the data analysis, the distribution $F$ does not always conform to elliptical contours and therefore may not be Gaussian. In such a setting, $p$-values computed using the Mahalanobis distance may not be reliable.

The Tukey half-space depth provides a general metric for identifying outliers in a flexible manner and in a non-parametric setting. Take $\mathbf{x}_{i}$, $i=1,\dots,n$ and $\mathbf{y}$ as above, making no assumptions on the structure of $F$ and $G$, and let $\mathbf{z}$ be any point in $\mathbb{R}^{m}$. Then the half-space depth $d_{\rm dep}(\mathbf{z};\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$ of $\mathbf{z}$ within the sample of the $\mathbf{x}_{i}$ is the smallest fraction of the $n$ points $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ to either side of any hyperplane passing through $\mathbf{z}$. By definition, the half-space depth is a number between 0 and 0.5. Points that have the same depth constitute a non-parametric estimate of the isolevel contour of the distribution $F$.

To evaluate a $p$-value for $\mathbf{y}$, we first compute $d_{j}=d_{\rm dep}(\mathbf{x}_{j};\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$ for every point $\mathbf{x}_{j}$, $j=1,\ldots,n$, yielding an empirical distribution of depth. Then the $p$-value is computed as the proportion of points whose depth is lower than that of $\mathbf{y}$:

[TABLE]

Note that, by construction, the depth $p$-value increases in units of $1/n$. For computing half-space depths below, we use the R package depth.

4 Results

We use the Planck maps for our analyses, which measure fluctuations about the mean in the CMB temperature to an accuracy of $10^{-5}$K (Planck Collaboration et al., 2014a). We compare observed maps with $1000$ Full Focal Plane 8 (FFP8) simulations cleaned using the NILC technique (Planck Collaboration et al., 2016b). The simulations are based on the LCDM paradigm and assume that the temperature fluctuations have a Gaussian distribution. In addition, we also compare the observed maps cleaned using the C-R, SEVEM and SMICA techniques with the NILC simulations, and note that they exhibit similar characteristics. However, since these maps arise from different processing techniques, the results, while broadly reproducing and corroborating those based on the NILC technique, are subject to various interpretations, and hence we abstain from discussing them. We perform our analyses for a range of scales between $0.05$ and $7.33$ degrees, which correspond to resolutions between ${N}=1024$ and ${N}=8$ in the HealPix format (Górski et al., 2005). Further degradation of the maps destabilizes the statistics due to the low number of data points in these cases. We do this for a range of mask binarization thresholds: 0.7, 0.8, 0.9 and 0.95. See Section 3 for the details of degradation and masking.

We present our analyses in terms of the ranks of relative homology groups, ${{b}_{p}}$ for $0\leq p\leq 1$. The relative components and loops are quantified by the relative component function, ${{b}_{0}}:{{\mathbb{R}}}\to{{\mathbb{R}}}$, and the relative loop function, ${{b}_{1}}:{{\mathbb{R}}}\to{{\mathbb{R}}}$. We present the graphs of ${{b}_{0}}$, ${{b}_{1}}$, and of the (relative) Euler characteristic, ${\sf EC_{\rm rel}}$, followed by statistical tests that estimate the significance of results. If $f(x):{{\mathbb{S}}}^{2}\to{{\mathbb{R}}}$ is the absolute temperature at a location $x$, and $f_{0}$ the mean temperature of the distribution, the dimensionless temperature is given by: $\nu(x)=(f(x)-f_{0})/\sigma(f)$, where $\sigma(f)$ is the standard deviation computed from the non-masked pixels. We then obtain the ranks of relative homology groups as functions of the normalized temperature.

4.1 Ranks of relative homology groups

To carry out omnibus tests, we choose 13 a-priori levels, $\ell_{-6},\dots,\ell_{6}$, where $\ell_{k}=k/2$, so that the normalized temperature thresholds run from $-3$ to $+3$ in steps of $0.5$, and consider collections of random variables ${{b}_{0}}(\ell_{k})$, ${{b}_{1}}(\ell_{k})$, and ${\sf EC_{\rm rel}}(\ell_{k})$, for $-6\leq k\leq 6$.

Figures 9–11 show the curves of ${{b}_{0}}$, ${{b}_{1}}$, and ${\sf EC_{\rm rel}}$ for resolutions between ${N}=1024$ and ${N}=8$, for mask threshold 0.9. The graphs present the average curve (black) from $1000$ NILC simulations, with error-bands drawn up to $3\sigma$. The individual curves from simulations are drawn in dotted black, a few of which escape the $3\sigma$ band. Also plotted are curves from NILC observed maps in yellow. ${{b}_{0}}(\nu)$, ${{b}_{1}}(\nu)$ and ${\sf EC_{\rm rel}}(\nu)$ show a difference from simulations peaking near $2\sigma$ for some temperature levels for all resolutions. Additionally, ${{b}_{0}}$ and ${{b}_{1}}$ show differences peaking between $3-4\sigma$ sporadically between $N=32$ and $N=8$.

4.2 Experimental evidence of Euler characteristic suppression

As noted before, the Euler-Poincaré formula states that the Euler characteristic is the alternating sum of the Betti numbers. As a consequence, the signals in Euler characteristic are suppressed by design, due to the cancellation of the constituent Betti numbers. Our experiments provide evidence for such suppressions of the topological signals emanating from the Euler characteristic. As an example, consider the quantities at the degraded resolution $N=16$, and temperature threshold value $\nu=0.5$. ${b}_{0}$ at this resolution and threshold has a significant difference between observations and simulations at $3.7\sigma$ (Figure 9), but the corresponding value for the Euler characteristic is $2.4\sigma$ (Figure 11). This is because of the cancellation effects between ${{b}_{0}}$ and ${{b}_{1}}$ in determining the Euler characteristic. One may find more instances of such cancellation effects in the graphs.

4.3 Statistical significance of the results

We consider the two methods detailed in Section 3.3, and present $p$-values of the observed maps for both. We consider the variables ${{b}_{0}}(\ell_{k=0,\ldots,6})$, ${{b}_{1}}(\ell_{k=-6,\ldots,0})$, and ${\sf EC_{\rm rel}}(\ell_{k=-6,\ldots,6})$ for estimating the statistical significance of the results. The choice of regions is determined by the fact that ${{b}_{0}}(\nu)$ tends to be small, and carries little information for $\nu<0$, ${{b}_{1}}(\nu)$ tends to be small for $\nu>0$, while the Euler characteristic is informative over the full range of levels. We perform summary and specific tests for mask binarization threshold values corresponding to 0.7, 0.8, 0.9 and 0.95. For the summary tests, we take the full set of quantities for all the resolutions between $N=1024$ and $N=8$, and the relevant thresholds, as vectors separately for each of the three topological quantities. The last entry for each mask threshold in Table 2 presents the summary $\chi^{2}$ and depth $p$-values. Overall, there is ample indication that the observations differ from the simulations. This is followed by specific tests for each resolution. The rest of the entries in Table 2 present the Mahalanobis and the depth $p$-values for each resolution, for different mask thresholds. The Mahalanobis distances are particularly small for ${{b}_{1}}$ at $N=16$, and very significant at $N=8$, across all binarization thresholds. ${{b}_{0}}$ and ${\sf EC_{\rm rel}}$ also show significance, but are not stable across binarization thresholds. The depth $p$-values are very significant for ${{b}_{1}}$ for $N=16$ and $N=8$, while ${{b}_{0}}$ shows high significance at $N=32$ consistently across the range of binarization thresholds. The depth $p$-values also show high significance at higher resolutions, more often for ${{b}_{1}}$ than for ${b}_{0}$, but not at all for ${\sf EC_{\rm rel}}$, presumably because of cancellation effects. ${{b}_{1}}$ shows significance more often than ${{b}_{0}}$ and ${\sf EC_{\rm rel}}$ when considering Tukey depth, and is an order of magnitude more significant compared to ${{b}_{0}}$ and ${\sf EC_{\rm rel}}$, when considering the Mahalanobis distance.

Regardless of the choice of test, it is evident that the model and observation disagree significantly at least in the number of loops on a range of scales approximately between 1 to 7 degrees. For the Mahalanobis values, the general trend of significance increases up to N = 8, providing additional evidence the deviations are not purely due to chance. For both tests, the $p$-values for the summary tests tend to be more significant than the for individual resolutions. Another general trend is that the non-parametric test shows the difference between the observed and the simulated maps to be starker than the parametric test. To help interpret this difference, Figures 12 – 15 shows plots that visualize to what extent the assumption of a Gaussian distribution for the compared quantities is justified. See also Table 3 for a comparison with $p$-values for the absolute homology. The results indicate similar trends as in the relative homology case.

4.4 Principal component graphs

Figure 12 presents the projection onto the first two principal components for the summary tests, which include results from all resolutions. Mahalanobis and depth contours corresponding to $p$-values of 0.1, 0.01 and 0.001 are shown in blue (top) and purple (bottom). Observed CMB points are in red. Examining the diagrams corresponding to the Mahalanobis distance, the hypothesis that the distribution conforms to elliptical contours is questionable.

Figures 13 – 15 present the projection onto the first two principal components respectively for ${{b}_{0}}$, ${{b}_{1}}$, and ${\sf EC_{\rm rel}}$ for specific resolutions. Also drawn are the Mahalanobis (top two rows) and Tukey depth (bottom two rows) contours. In general, it is the case that the symmetric Mahalanobis contours do not always fit the data. However, as the resolution decreases, the Mahalanobis contours, which are Gaussian in nature, seem to fit the data well, and may be a reasonable approximation after all.

5 Summary and Conclusions

We provide evidence for the deviation of the observed Planck CMB maps from the Gaussian predictions of the standard LCDM model. Specifically, we find an over-abundance of loops in the observed maps, deviating from the simulations at per cent to less than per mil levels. This is in terms of $p$-values computed using $\chi^{2}$ statistics, between the resolutions $N=32$ and $N=8$. The difference in the number of components and loops peaks sporadically at more than $3\sigma$ from the predictions between $N=32$ and $N=8$. Results based on smoothed maps corroborate with those based on degraded maps in terms of approximate scales at which the anomaly is observed. We also compute the absolute homology for the dataset, and confirm that the results are consistent with those from relative homology. External evidence that these deviations are not a result of overanalysing the data comes from the fact that the variance of the observed CMB is anomalous with respect to the standard model at $N=16$ (Planck Collaboration et al., 2014b), and the computed power spectrum exhibits a dip roughly at this range of scales. In addition, there are reports of a mildly significant Euler characteristic at $3.66$ degrees ($N=16$) (Eriksen et al., 2004b), computed from independent measurements of the CMB by Planck’s predecessor – Wilkinson Microwave Anisotropy Probe (WMAP) satellite. This can be explained by the significantly high number of loops and components, together with cancellation effects that the Euler characteristic suffers from, by definition. Similar observations by independent satellites makes it unlikely that the source of the anomaly has its origin in instrumental noise or systematic effects. Moreover the medium super-horizon scales at we observe it, could possibly point to a cosmological origin. The non-parametric Tukey depth test shows the observations to be different from the simulations at almost all resolutions. Regardless of the preferred test, evidently the topological structure of the CMB is deviant from the simulations, at least on some scale.

We can rule out this anomaly being the effect of the cold spot in the CMB sky, or any previously detected directional anomalies. Our statistics are based on a large numbers of loops surrounding the low density regions, to which the loop generated by the cold spot may contribute at most only a few, and often only one. Moreover, to support this claim, we visually confirm that these loops are scattered all over the sky (see Figure 16). We also test and confirm that simulations that are based on Gaussian prescriptions and match the characteristics of the observed “dipped” power spectrum cannot resolve this anomaly. Additionally, we present topological methods that are suitable in the presence of obfuscating masks. As such, the results presented in this paper are robust despite lacking full sky coverage, as well as being model-independent.

In conclusion, we reiterate that we present clear evidence of departure of the observed CMB maps with respect to the simulations based on the LCDM paradigm, but make no attempt to address the issue of the physical mechanism behind this phenomenon; a question we leave to the wider cosmological community. Nevertheless, our analysis demonstrates the existence of unexpected topology in the CMB. Possible, but non-exhaustive scenarios worth exploring may be primordial non-Gaussianity, as well as models with non-trivial topology.

Acknowledgements

PP and RA acknowledge the support of ERC advanced grant Understanding Random Systems through Algebraic Topology (URSAT) (no: 320422, PI: RA). This work is also part of a project that has received funding for PP and TB from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement ERC advanced grant 740021– Advances in Research on THeories of the dark UniverSe (ARTHUS), PI: TB). HE and HW acknowledge the support by the Office of Naval Research, through grant N62909-18-1-2038, and by the DFG Collaborative Research Center TRR 109, ‘Discretization in Geometry and Dynamics’, through grant I02979-N35 of the Austrian Science Fund (FWF). PP is grateful to Julian Borill from the Planck consortium for providing the data, and for the illuminating discussions and inputs. PP also thanks Hans Kristen Eriksen and Anne Ducout for numerous discussions in the early stages.