# Loops abound in the cosmic microwave background: A $4\sigma$ anomaly on
super-horizon scales

Pratyush Pranav  Univ Lyon, ENS de Lyon, Univ Lyon1, CNRS, Centre de Recherche
Astrophysique de Lyon UMR5574, F–69007, Lyon, France

###### Abstract

We present a topological analysis of the temperature fluctuation maps from the
_Planck 2020_ Data release 4 (DR4) based on the NPIPE data processing
pipeline. For comparison, we also present the topological characteristics of
the maps from _Planck 2018_ Data release 3 (DR3). We perform our analysis in
terms of the homology characteristics of the maps, invoking relative homology
to account for analysis in the presence of masks. We perform our analysis for
a range of smoothing scales spanning sub- and super-horizon scales
corresponding to
$FWHM=5^{\prime},10^{\prime},20^{\prime},40^{\prime},80^{\prime},160^{\prime},320^{\prime},640^{\prime}$.
Our main result indicates a significantly anomalous behavior of the loops in
the observed maps compared to simulations that are modeled as isotopic and
homogeneous Gaussian random fields. Specifically, we observe a $4\sigma$
deviation between the observation and simulations in the number of loops at
$FWHM=320^{\prime}$ and $FWHM=640^{\prime}$, corresponding to super-horizon
scales of $5$ degrees and larger. In addition, we also notice a mildly
significant deviation at $2\sigma$ for all the topological descriptors for
almost all the scales analyzed. Our results show a consistency across
different data releases, and therefore, the anomalous behavior deserves a
careful consideration regarding its origin and ramifications. Disregarding the
unlikely source of the anomaly being instrumental systematics, the origin of
the anomaly may be genuinely astrophysical – perhaps due to a yet unresolved
foreground, or truly primordial in nature. Given the nature of the topological
descriptors, that potentially encodes information of all orders, non-
Gaussianities, of either primordial or late-type nature, may be potential
candidates. Alternate possibilities include the Universe admitting a non-
trivial global topology, including effects induced by large-scale topological
defects.

##  1 Introduction

At the epoch of recombination, matter and radiation separate, allowing
radiation to stream freely in the Universe. This free-streaming radiation
permeating the Universe, that we observe as the _Cosmic Microwave Background_
(CMB) radiation, encodes a treasure trove of information about the initial
conditions in the Universe (Ryden, 2003; Jones, 2017). Despite having a
remarkably consistent average temperature, the CMB still exhibits tiny
deviations of the order of $10^{-5}$ from the background average. The
temperature fluctuations in the CMB trace the fluctuations in the underlying
mater distribution in the infant Universe, that are linked to the spontaneous
quantum fluctuations generated in an otherwise homogeneous medium (Harrison,
1970; Peebles & Yu, 1970). Thus studying the properties of the temperature
fluctuations in the CMB is essential towards understanding the properties of
the primordial matter field.

The Lambda Cold Dark Matter (LCDM) paradigm is the standard paradigm of
cosmology. Together with the inflationary models in their simplest form
(Starobinsky, 1982; Guth & Pi, 1982), the standard model of cosmology predicts
the nature of the primordial stochastic matter distribution field to be that
of an isotropic and homogeneous Gaussian random field (Harrison, 1970; Guth,
1981). This prediction finds allies theoretically in the Central Limit
Theorem, and observationally in the various measurements of the CMB
temperature anisotropy field via ground- and space-based probes such as the
BOOMERanG balloon-based experiment Masi (2002) (Jaffe et al., 2001), and the
Wilkinson Microwave Anisotropy Probe (WMAP) satellite (Bennett et al., 2013).
The latest endeavor of measuring the CMB temperature anisotropies has
materialized through the launch of the _Planck_ satellite, which boasts of the
highest resolution in measurements till date, at scales of a few arc-minutes
(Planck Collaboration et al., 2020a). Despite the general consensus that the
CMB exhibits the characteristics of an isotropic and homogeneous Gaussian
random field, there is a growing body of evidence regarding anomalies in the
observed CMB field with respect to the base model. Notable among them are the
observed hemispherical asymmetry in the CMB power spectrum (Eriksen et al.,
2004a), as well as the alignment of low multipoles (Copi et al., 2015); see
Schwarz et al. (2016) for a review. These observed anomalies raise doubts
about the assumption of statistical isotropy and homogeneity respectively.

Testing the assumption of Gaussianity requires tools which encode information
about higher orders. Traditional endeavor in this direction has focused on
higher order correlation functions, which are generally extremely resource
intensive computationally (Planck Collaboration et al., 2016b). Recently,
attention has turned towards developing alternative tools beyond the
correlation functions and multi-spectra, which may potentially encode
information of all orders. Principal such tools have arisen from integral
geometry, and involve computing the _Minkowski functionals_ or the _Lifshitz-
Killing curvatures_ (Adler, 1981; Mecke et al., 1994; Schmalzing & Gorski,
1998; Sahni et al., 1998; Codis et al., 2013; Ducout et al., 2013; Matsubara,
2010; Chingangbam et al., 2017; Pranav et al., 2019b; Telschow et al., 2019).
The $j$-th Minkowski functional and $(D-j)$-th Lifshitz-Killing curvature of a
$D$-dimensional manifold ${{\mathbb{M}}}$ are related by
$Q_{j}({{\mathbb{M}}})\ =\ j!\omega_{j}{\cal L}_{D-j}({{\mathbb{M}}})$, where
$j=0,\dots,D,$ and $\omega_{j}$ is the volume of the $j$-dimensional unit
ball. There are $d$ such quantifiers for a $D$-dimensional set, where
$d=0,\ldots,D$. All but one are purely geometrical quantities, related to the
$d$-dimensional volume of the manifold. The exception is the $0$-th Lifshitz-
Killing curvature, or equivalently the $D$-th Minkowski functional, which is
related to a purely topological quantity, the _Euler characteristic_ (Euler,
1758; Adler, 1981), via Gauss’s _Theorema Egrerium_ (Gauss, 1900; Adler, 1981;
Pranav et al., 2019b). The Minkowski functional computations of the CMB have
consistently shown the observations to be congruent with the standard model
(Planck Collaboration et al., 2016a).

More recently, developments in computational topology have paved way for
extracting topological information from datasets, at the level of _homology_
(Munkres, 1984; Edelsbrunner & Harer, 2010; van de Weygaert et al., 2011;
Pranav, 2015; Pranav et al., 2017), and its hierarchical extension _persistent
homology_ (Edelsbrunner & Harer, 2010; Pranav et al., 2017). _Topological data
analysis_ (TDA) involving homology and persistent homology has started finding
application in astrophysical disciplines recently, for example in the context
of structure identification (Shivashankar et al., 2016; Xu et al., 2019) and
quantification of large-scale structures (van de Weygaert et al., 2011; Kono
et al., 2020; Wilding et al., 2020), including detection and quantification of
non-Gaussianities (Feldbrugge et al., 2019; Biagetti et al., 2020). Homology
describes the topology of a space by identifying the holes and the topological
cycles that bound them. A $d$-dimensional space may contain topological cycles
of $0$ up to $d$-dimensions. The cycles and holes are associated with the
_homology groups_ of the space. The $p$-th _Betti number_ , ${{\beta}_{p}}$,
is the rank of the $p$-th homology group, $\mathbb{H}_{p,p=0\ldots d}$. While
itself a purely topological quantity, the Euler characteristic is also the
alternating sum of the Betti numbers of all ambient dimensions of a manifold,
as denoted by the Euler-Poincaré formula (Adler, 1981; Pranav et al., 2017,
2019b). The Euler characteristic has a long history in the analysis of
cosmological fields (Gott et al., 1986; Pogosyan et al., 2009; Park et al.,
2013; Appleby et al., 2020).

Building on and refining existing tools from computational topology, in the
context of analyzing the CMB field, this paper presents the homology
characteristics of the temperature fluctuation maps of the cosmic microwave
background obtained by the _Planck_ satellite (Planck Collaboration et al.,
2020a). We perform our experiments on the fourth and the final data release
_Planck 2020_ Data release 4 (DR4), which is based on the NPIPE data
processing pipeline (Planck Collaboration et al., 2020b). The NPIPE dataset
represents a natural evolution of the Planck data processing pipeline,
integrating the best practices from the LFI and HFI pipelines separately. The
result is an overall amplification of signal and reduction in the associated
systematic, noise and residuals at almost all angular scales (Planck
Collaboration et al., 2020b). For comparison, we also present results for the
_Planck 2018_ Data release 3 (DR3) (Planck Collaboration et al., 2020a), which
is based on the _Full Focal Plane_ (FFP) data processing pipeline (Planck
Collaboration et al., 2016c) resulting in the FFP10 simulations (Planck
Collaboration et al., 2020d). The paper follows the spirit of Pranav et al.
(2019a) in methods and analysis, where we present results for the intermediate
_Planck 2015_ Data release 2 (DR2) (Planck Collaboration et al., 2016a). The
novel aspect of the methodologies presented here and in Pranav et al. (2019a)
is an analysis pipeline that takes into account regions with unreliable data
on ${{\mathbb{S}}}^{2}$. In the case of CMB, this is reflected in the
obfuscation effects that the measurements suffer from due to the foreground
objects such as our own galaxy, as well as other extra- and intra-galactic
foreground sources. We mask such regions and compute the homology of the
excursion sets relative to the mask.

Our main results indicate an anomalous behavior of loops when comparing the
observational maps to the simulations. Specifically, we detect a $4\sigma$
deviation in the number of loops between the observational maps and
simulations at scales of approximately $5$ degrees and larger. This is on top
of a generally deviant behavior of both components and loops at around
$2\sigma$ for almost all scales analyzed, when comparing observations to
simulations. The Euler characteristic, being influenced by both the components
and the loops, shows commensurate deviations. Even though differing in
details, the results show a generally consistent trend across datasets. The
results merit a serious consideration in view of the twin facts that the the
data processing pipeline employed by the Planck team has evolved consistently,
resulting in increasingly more accurate temperature maps, and that the
temperature maps are in excellent agreement across data releases, not least
due to a high signal-to-noise ratio a-priori.

We present a brief description of the topological background in Section 2,
followed by the results in Section 3, which are based on the NPIPE dataset. We
discuss the ramifications of the results and conclude the main body of the
paper in Section 4. The appendices present a brief description of the datasets
and the computational pipeline, as well as additional results based on the
FFP10 dataset.

(a)  Figure 1: A visualization of the temperature fluctuations in the CMB sky.
The survey surface ${{\mathbb{S}}}^{2}$ is distorted at each point in the
direction of the surface normal. The distortion is proportional to the
fluctuation in direction and magnitude. The visualization is based on the
observed CMB sky cleaned by the NPIPE pipeline and smoothed at $5$ degrees.

(a)

(b)

Figure 2: The CMB sky thresholded at moderately positive (left) and negative
(right) levels. For high thresholds, the excursion set is dominated by
isolated components, while at low thresholds it gives the appearance of a
single connected surface indented by numerous holes. For sufficiently low
thresholds, the holes fill up, and the excursion set covers the entire sphere,
composed of a connected surface without boundary, enclosing a single void (cf.
Figure 1). (a)  Figure 3: A visualization of the temperature fluctuations in
the CMB sky in the presence of masks (painted in grey). It consists of an
equatorial belt and numerous patches on the northern and the southern cap,
corresponding to our galaxy and other bright foreground objects. The
visualization is based on the observed CMB sky cleaned by the NPIPE pipeline
and smoothed at $5$ degrees.

##  2 Topological background

Commensurate with our intention of analyzing the topology of the CMB
temperature fluctuations, we restrict ourselves to topological definitions on
${{\mathbb{S}}}^{2}$ (Pranav et al., 2019a), and invoke _relative homology_ to
account for analysis in the presence of masked regions. Standard reference for
this section is Edelsbrunner & Harer (2010); also see Pranav et al. (2019a)
for discussion in the context of CMB.

###  2.1 Homology characteristics of excursion sets of ${{\mathbb{S}}}^{2}$

Denoting the CMB temperature fluctuations on ${{\mathbb{S}}}^{2}$ as
$f\colon{{\mathbb{S}}}^{2}\to{{\mathbb{R}}}$, we define the excursion set at a
temperature $\nu$ as the subset of ${{\mathbb{S}}}^{2}$ where the temperature
is larger than or equal to $\nu$:

${{\mathbb{E}}}(\nu)=\\{x\in{{\mathbb{S}}}^{2}\mid f(x)\geq\nu\\}.$ (1)

If ${{\mathbb{E}}}(\nu)$ does not cover entire ${{\mathbb{S}}}^{2}$, it may be
composed of isolated components and holes. Figure 2 presents excursion sets
corresponding to two different thresholds. For high thresholds, presented in
the left panel, the excursion set is dominated by components, while for low
thresholds, presented in the right panel, the excursion set is dominated by a
few large connected objects indented with holes, that are bounded by loops.
The _Betti numbers_ ${{\beta}_{0}}$ and ${{\beta}_{1}}$ count the number of
independent components and loops of the excursion set respectively. In
general, for a $d$-dimensional topological space, ${{\beta}_{p}}$ is the rank
of the $p$-th _homology group_ , ${{\sf H}_{p}};p=0,\ldots,d$, and counts the
number of independent $p$-dimensional cycles (Munkres, 1984; Edelsbrunner &
Harer, 2010; Pranav et al., 2017). If ${{\mathbb{E}}}(\nu)$ does not cover
entire ${{\mathbb{S}}}^{2}$, the number of independent loops is one less than
the total number of loops. If ${{\mathbb{E}}}(\nu)$ covers entire
${{\mathbb{S}}}^{2}$, there are no loops, and ${{\beta}_{2}}=1$, because of
the void enclosed by the boundary-less surface of the sphere. A related
quantity that has a long history of usage in cosmological analyses is the
_Euler characteristic_ , or alternatively the _genus_ (Gott et al., 1986; Park
et al., 2013), which is the alternating sum of the Betti numbers of the
excursion set:

${\sf EC}(\nu)={{\beta}_{0}}(\nu)-{{\beta}_{1}}(\nu)+{{\beta}_{2}}(\nu).$ (2)

The Euler characteristic also has a geometric interpretation as one of the
_Lifshitz-Killing_ curvatures of the manifold (Adler, 1981; Pranav et al.,
2019a).

###  2.2 Masks and relative homology

The measurement of CMB signal is unreliable in certain parts of the sky due to
interference from bright foreground objects. These include extended objects
such as our galaxy, as well as bright point sources. We mask such regions, and
compute the homology characteristics of the excursion set relative to the
mask. Figure 3 presents a visualization of the masked CMB sky. Letting
${{\mathsf{M}}}\subseteq{{\mathbb{S}}}^{2}$ be the mask, and
${{\mathbb{E}}}(\nu)$ the excursion set, we consider the _relative homology_
of the pair of closed spaces, $(E,M)$, where $E={{\mathbb{E}}}(\nu)$ and
$M={{\mathsf{M}}}\cap{{\mathbb{E}}}(\nu)$. Note that $M$ is contained in $E$.
We denote the rank of the _relative homology groups_ of the pair $(E,M)$ by
${{b}_{p}}={{\rm rank\,}{{{\sf H}_{p}}(E,M)}};p=0,1,2$. The Betti numbers
computed considering the pair $(E,M)$ are different from the Betti numbers of
the excursion set without a mask. For a more detailed discussion about
relative homology in the context of masked CMB sky see Pranav et al. (2019a).
The _relative Euler characteristic_ , as in the case of absolute homology, is
the alternating sum of the rank of relative homology groups:

${\sf EC_{\rm rel}}(\nu)={{b}_{0}}(\nu)-{{b}_{1}}(\nu)+{{b}_{2}}(\nu).$ (3)

##  3 Results

We present our results in terms of the ranks of relative homology groups,
${{b}_{p}}$ for $0\leq p\leq 1$. 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. The
main paper presents results from the DR4 NPIPE dataset, based on $600$
simulations, obtained using the SEVEM component separation pipeline. Similar
results for the the DR3 FFP10 dataset, based on $300$ simulations, obtained
using the SMICA component separation pipeline, are presented in the appendix
for comparison.

(a)

(b)

(c)

Figure 4: Significance of difference for ${{b}_{0}}$, ${{b}_{1}}$ and ${\sf
EC_{\rm rel}}$. The observational curve is presented in red, and the curves
corresponding to simulations are presented in gray. The plotted quantity is
the difference between the observation and the mean, scaled by the standard
deviation, the latter two computed from the simulations. We notice a $4\sigma$
deviation between observation and simulations in the number of loops at
${N}=32$, and ${N}=16$, corresponding to $FWHM=320^{\prime}$ and
$FWHM=640^{\prime}$. This is on top of a $2\sigma$ deviation that we notice
generally for both components and loops across almost all resolutions. The
Euler characteristic reflects the deviations in components and loops, albeit
at smaller significance, due to cancellation effects (Pranav et al., 2019a).
Relative homology
---
|  | Mahalanobis | Tukey Depth
Res | FWHM | ${{b}_{0}}$ | ${{b}_{1}}$ | ${\sf EC_{\rm rel}}$ | ${{b}_{0}}$ | ${{b}_{1}}$ | ${\sf EC_{\rm rel}}$
threshold = 0.90
2048 | 5 | 0.768 | 0.570 | 0.858 | 0.630 | 0.318 | 0.912
1024 | 10 | 0.780 | 0.599 | 0.383 | 0.672 | 0.352 | 0.000
512 | 20 | 0.280 | 0.549 | 0.203 | 0.332 | 0.308 | 0.000
256 | 40 | 0.585 | 0.746 | 0.772 | 0.362 | 0.560 | 0.642
128 | 80 | 0.120 | 0.353 | 0.563 | 0.000 | 0.312 | 0.730
64 | 160 | 0.401 | 0.671 | 0.768 | 0.465 | 0.718 | 0.755
32 | 320 | 0.040 | 0.004 | 0.001 | 0.000 | 0.000 | 0.000
16 | 640 | 0.987 | 0.005 | 0.203 | 0.983 | 0.000 | 0.440
summary | NA | 0.478 | 0.023 | 0.045 | 0.803 | 0.000 | 0.000
(a)  Table 1: Table displaying the two-tailed $p$-values for relative homology
obtained from parametric (Mahalanobis distance) and non-parametric (Tukey
depth) tests, for different resolutions and smoothing scales for the NPIPE
dataset. The last entry is the $p$-value for the summary statistic computed
across all resolutions. Marked in boldface are $p$-values $0.05$ or smaller.
(a)
Figure 5: A visualization of the loops surrounding the low density regions at
a moderately negative threshold. The gaps in the manifold correspond to the
mask that has been removed. We clearly notice the equatorial belt
corresponding to our galaxy, and more patches in the visible cap. Some loops
live fully in the excursion region, not influenced by the mask. We color them
in red. We also depict some representative relative loops, that have their end
points in mask. We draw closed loops for this category as well, coloring the
portions in white where they overlap with the masked regions. The
visualization is based on the observed CMB map from the NPIPE data set
smoothed at $5$ degrees – approximately the scale at which the loops start
showing deviant behavior.

###  3.1 Graphs of ${{b}_{0}}$, ${{b}_{1}}$ and ${\sf EC_{\rm rel}}$

We choose a number of a-priori levels, $\ell_{k}$, where $\ell_{k}=k/2$,
setting $k_{{{b}_{0}}}=[0:6]$, $k_{{{b}_{1}}}=[-6:0]$, and $k_{{\sf EC_{\rm
rel}}}=[-6:6]$. We do so to restrict ourselves to analyzing ${{b}_{0}}$ for
positive thresholds, ${{b}_{1}}$ for negative thresholds, and ${\sf EC_{\rm
rel}}$ across the full threshold range. 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
consider collections of random variables ${{b}_{0}}(\ell_{k})$,
${{b}_{1}}(\ell_{k})$, and ${\sf EC_{\rm rel}}(\ell_{k})$. For each level, we
compute the mean $\mu_{sim}$ and the standard deviation $\sigma_{sim}$ from
the simulations, for each of the topological descriptor. If the observed value
is $f_{obs}$, the normalized difference, $\Delta$ is given by:

$\Delta=\frac{f_{obs}-\mu_{sim}}{\sigma_{sim}}.$ (4)

Figure 4 presents the curves for the normalized differences for ${{b}_{0}}$
(top two rows), ${{b}_{1}}$ (middle two rows), and ${\sf EC_{\rm rel}}$
(bottom two rows). For each quantity, from top-left to bottom-right, we
present plots corresponding to the $8$ different resolutions and smoothing
scales enumerated in Table 2. The normalized difference curves for the
observational map is presented in red, overlapped with similar curves for the
simulations, presented in gray. The first thing we notice is that the
observational curves for both ${{b}_{0}}$ and ${{b}_{1}}$ show mildly
significant deviation at approximately $2\sigma$ across resolutions, for at
least some thresholds. A more interesting observation is the highly
significant deviation of the observational curves corresponding to relative
loops at $4\sigma$ for ${{b}_{1}}$, at ${N}=32$ and ${N}=16$, corresponding to
$FWHM=320^{\prime}$ and $FWHM=640^{\prime}$. ${{b}_{1}}$ quantifies the number
of independent loops in the manifold. Figure 5 presents a visualization of
some such loops at moderately negative thresholds, for observational maps
smoothed at $5$ degrees. The smoothing scale is approximately the scale at
which the loops start exhibiting a deviant behavior from the base model. The
relative Euler characteristic also shows deviations for these resolutions,
albeit at a slightly smaller significance, owing to cancellation effects
(Pranav et al., 2019a). However, the Euler characteristic is not strictly an
independent quantity, due to it being an alternating sum of the ranks of the
relative homology groups. As such, it merely reflects the deviations in the
contributing Betti numbers. For comparison, we also show graphs of the
topological quantities computed from the DR3 FFP10 dataset in Figure 8. We
notice similar trends in the behavior of the three topological descriptors as
the case presented in the main body.

###  3.2 Statistical significance of the results

We consider the two methods detailed in Section 2, and present $p$-values of
the observed maps for both. Considering all the three topological quantities
${{b}_{0}}$, ${{b}_{1}}$ and ${\sf EC_{\rm rel}}$ we first compute the
statistics combining all resolutions, in order to probe for non-random
deviations. Thereafter, we compute the statistics for all resolutions
separately, to ascribe a scale-dependence to the signals.

Summary tests. To test for the evidence of a non-random discrepancy, we take
the full set of measurements for the relevant thresholds, combined across
resolutions, for the three variables. This results in a data vector of $m=56$
entries for ${{b}_{0}}$ and ${{b}_{1}}$, and $m=104$ for ${\sf EC_{\rm rel}}$.
The length of the vector is commensurate with the fact that we analyze
${{b}_{0}}$ and ${{b}_{1}}$ for positive and negative thresholds respectively,
while we analyze ${\sf EC_{\rm rel}}$ for the full range of levels. The last
entry in the Table 1 presents the $p$-values for the variables for both the
parametric and the non-parametric tests. Overall, there is an indication that
the observations are deviant from the simulations, at least in the the number
of loops, which is an independent quantity. By association, Euler
characteristic also shows deviations. We also note that the non-parametric
test shows higher significance compared to the parametric test. Table 3 in the
appendix presents the $p$-values for the FFP10 dataset. The trends are
similar, albeit with higher significance. The Mahalanobis test indicates
stronger deviation for the FFP10 dataset compared to the NPIPE dataset.
Similarly, the Tukey depth test exhibits more instances of highly significant
deviations for all the three variables, where the number of simulations are
not enough to resolve the $p$-values.

Tests at specific resolutions. All entries before the last in Table 1 present
$p$-values for the variables, computed from maps degraded at specific
resolutions detailed in Table 2. Considering the Mahalanobis distance,
${{b}_{1}}$ shows highly significant difference between the observational maps
and the simulations for ${N}=32$ and ${N}=16$. Additionally, ${{b}_{0}}$ also
shows mildly significant difference at ${N}=32$, but the differences are an
order smaller than what is exhibited by ${{b}_{1}}$. The Euler characteristic
shows significant difference at ${N}=32$. This is due to the significant
difference shown by the contributing ${{b}_{0}}$ and ${{b}_{1}}$ at this
resolution. In general, due to the Euler characteristic being an alternating
sum of the Betti numbers, its behavior is influenced by both the Betti
numbers. As an example, for the next lower resolution, ${N}=16$, the Euler
characteristic shows no significant difference, even though ${{b}_{1}}$
exhibits significant difference. This is because of the highly non-significant
behavior of ${{b}_{0}}$, whose contribution ends up canceling the contribution
of ${{b}_{1}}$ towards the Euler characteristic. Results for the FFP10 dataset
are presented in Table 3 in the appendix. As in the summary test case, the
Mahalanobis test indicates stronger deviation for the FFP10 dataset compared
to the NPIPE dataset in general. The Tukey depth test exhibits more instances
of highly significant deviations across different resolutions for all the
three variables, where the number of simulations are not enough to resolve the
$p$-values.

##  4 Discussions and conclusion

We present a topological analysis of the temperature fluctuations in the CMB
in terms of homology. To account for regions with unreliable data, we compute
the homology of the excursion sets relative to the mask covering these areas.
We perform our analysis on the fourth and final NPIPE data release from the
Planck team. The pipeline represents a natural evolution of the data
processing pipeline, commensurate with better understanding of systematics,
residuals and noise over a period of time across successive data releases. It
incorporates the best strategies from both the LFI and HFI processing
pipeline, so that there is an overall reduced level of noise and residuals
across all scales (Planck Collaboration et al., 2020b). We also investigate
the _Planck 2018_ Data release 3 (DR3), accompanied by the FFP10 simulations
(Planck Collaboration et al., 2020d), for comparison and completeness. The
present paper is a successor to Pranav et al. (2019a), where we investigate
the topology of the temperature fluctuation maps based on the intermediate
_Planck 2015_ Data release 2 (DR2), accompanied by FFP8 simulations. Between
the two papers, we have investigated topological characteristics of the CMB
temperature fluctuation maps for the latest three data releases by the Planck
team, and the overall trends in the results across the datasets inform on the
consistency of the data processing pipeline, as well as our own methodologies.

Our main result indicates an anomaly in the behavior of the loops in the
observed sky compared to the simulations based on the standard LCDM model that
predicts the nature of the primordial perturbations to be that of an isotropic
and homogeneous Gaussian random field. The number of loops in the
observational maps shows significantly anomalous deviations from the
simulations at $4\sigma$, at scales of $5$ degrees and larger. We also notice
a general mild flaring of all the topological quantities computed from the
observational maps at approximately $2\sigma$ for almost all resolutions. In
order to test for non-random discrepancies, we compute the $p$-values using
the parametric $\chi^{2}$ test, as well as the non-parametric Tukey depth
test. For the parametric test, the number of loops shows per-mil deviations
between observations and simulation, which is an order of magnitude larger
than the deviation shown by the number of components. Non-parametric tests
exhibit trends that are starker, more often for the loops compared to the
other quantities, in the sense that the number of simulations turn out to be
inadequate to resolve the $p$-values reliably. The trend in our results show a
consistency across the latest data releases NPIPE and PR3 analyzed in this
paper, as well as the intermediate data release PR2 analyzed in (Pranav et
al., 2019a). The trend is commensurate with the general observation by the
Planck team that the quality of the temperature maps has been consistent since
the second data release PR2 (Planck Collaboration et al., 2016a). The
consistency of our results across the three datasets independently confirms
this.

The origin and repercussions of the observed anomaly deserve careful
consideration. Concerning the origin of the anomaly, the least interesting,
and a highly unlikely possibility is instrumental systematics as source. This
is considering the twin facts that the results are consistent across data
releases and different processing pipelines, and that the Euler characteristic
computed from the CMB maps obtained by _Planck_ ’s predecessor Wilkinson
Microwave Anisotropy Probe (WMAP) satellite also exhibits mildly significant
deviations between observations and simulations (Eriksen et al., 2004b).
Disregarding systematics, a more interesting possibility is the source of the
anomaly being a genuine astrophysical signal, perhaps truly primordial in
nature, or due to a yet unresolved foreground effect. If the signals are
primordial, it opens up the possibility of the Universe admitting a non-
trivial global topology (Aurich & Steiner, 2001; Roukema et al., 2004; Di
Valentino et al., 2019), possibly induced by large-scale topological defects
(Durrer, 1999; Vilenkin & Shellard, 2000; Bouchet et al., 2001; Gangui, 2002),
as well as that of primordial non-Gaussianity (Verde et al., 2001; Planck
Collaboration et al., 2016b). The latter scenario has justification based on
the nature of the method used to detect anomalies. Methods based on comparing
the power spectrum of different patches of the sky, as example hemispheres,
inform about homogeneity properties (Eriksen et al., 2004a; Adler et al.,
2017), where as investigating the alignment of multipoles informs about
isotropy properties (Copi et al., 2015). Results based on these methods do not
encode information about higher orders, and hence cannot shed light on non-
Gaussianities. Computing higher order correlations is expensive, however $3$
-point correlation functions show the observations and simulations to be
consistent (c.f. Planck Collaboration et al. (2016a, 2020c). Alternative
methods, potentially encoding information of all orders emerge from geometry
and topology. Principal such tools from integral geometry are the Minkowski
functionals (Adler, 1981; Mecke et al., 1994; Matsubara, 2010; Ducout et al.,
2013), that also show consistent behavior between simulations and observations
of CMB (Planck Collaboration et al., 2016a, 2020c). In this scenario, the
purely topological tools employed in this paper present, for the first time,
an anomaly that may have contributing influence from higher orders of
correlation, thereby potentially pointing to a non-Gaussian signal. Future
research will involve a coordinated effort in all the aforementioned
directions to establish a deeper understanding of the origin and nature of the
anomaly.

## Acknowledgements

I am greatly indebted to Robert Adler, Thomas Buchert, Herbert Edelsbrunner,
Bernard Jones, Armin Schwarzman, Gert Vegter, and Rien van de Weygaert for
encouraging this solo venture, and for acting as mentors over years. My
gratitude also to Julian Borrill and Reijo Keskitalo for their patience in
clarifying doubts, and their constructive comments on the draft. I also thank
Tal Eliezri for insightful comments on the artwork. This work is supported by
the ERC advanced grant ARThUs (grant no: 740021; PI: TB), with contributing
influence from ERC advanced grant URSAT (grant no: 320422; PI: RA). I
gratefully acknowledge the support of PSMN (Pôle Scientifique de Modélisation
Numérique) of the ENS de Lyon, and the Department of Energy’s National Energy
Research Scientific Computing Center (NERSC) at Lawrence Berkeley National
Laboratory, operated under Contract No. DE-AC02-05CH11231, for the use of
computing resources.

## References

  * Adler (1981) Adler, R. 1981, The Geometry of Random Fields, Classics in applied mathematics (Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104))
  * Adler et al. (2017) Adler, R. J., Agami, S., & Pranav, P. 2017, Proceedings of the National Academy of Sciences, 114, 11878
  * Appleby et al. (2020) Appleby, S., Park, C., Hong, S. E., Hwang, H. S., & Kim, J. 2020, Astrophysical Journal, 896, 145
  * Aurich & Steiner (2001) Aurich, R. & Steiner, F. 2001, Mon. Not. Royal Astro. Soc., 323, 1016
  * Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, Astrophys. J. Suppl., 208, 20
  * Biagetti et al. (2020) Biagetti, M., Cole, A., & Shiu, G. 2020, arXiv e-prints, arXiv:2009.04819
  * Bouchet et al. (2001) Bouchet, F. R., Peter, P., Riazuelo, A., & Sakellariadou, M. 2001, Phys. Rev. D, 65, 021301
  * Chingangbam et al. (2017) Chingangbam, P., Yogendran, K. P., Joby, P. K., et al. 2017, Jour. Cos. and Part. Phys., 12, 023
  * Codis et al. (2013) Codis, S., Pichon, C., Pogosyan, D., Bernardeau, F., & Matsubara, T. 2013, Mon. Not. Royal Astro. Soc., 435, 531
  * Copi et al. (2015) Copi, C., Huterer, D., Schwarz, D., & Starkman, G. 2015, Monthly Notices of the Royal Astronomical Society, 449, 3458
  * Di Valentino et al. (2019) Di Valentino, E., Melchiorri, A., & Silk, J. 2019, Nature Astronomy
  * Ducout et al. (2013) Ducout, A., Bouchet, F. R., Colombi, S., Pogosyan, D., & Prunet, S. 2013, MNRAS, 429, 2104
  * Durrer (1999) Durrer, R. 1999, New. Astro. Review, 43, 111
  * Edelsbrunner & Harer (2010) Edelsbrunner, H. & Harer, J. 2010, Computational Topology - an Introduction (American Mathematical Society), I–XII, 1–241
  * Eriksen et al. (2004a) Eriksen, H. K., Hansen, F. K., Banday, A. J., Górski, K. M., & Lilje, P. B. 2004a, Astrophysical Journal, 605, 14
  * Eriksen et al. (2004b) Eriksen, H. K., Novikov, D. I., Lilje, P. B., Banday, A. J., & Górski, K. M. 2004b, Astrophysical Journal, 612, 64
  * Euler (1758) Euler, L. 1758, Novi Commentarii academiae scientiarum Petropolitanae, 4, 140
  * Feldbrugge et al. (2019) Feldbrugge, J., van Engelen, M., van de Weygaert, R., Pranav, P., & Vegter, G. 2019, Journal of Cosmology and Astroparticle Physics, 2019, 052–052
  * Gangui (2002) Gangui, A. 2002, International Journal of Modern Physics A, 17, 4273
  * Gauss (1900) Gauss, C. F. 1900, K. Gesellschaft Wissenschaft, 8
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, Astrophysical Journal, 622, 759
  * Gott et al. (1986) Gott, III, J. R., Dickinson, M., & Melott, A. L. 1986, Astrophysical Journal, 306, 341
  * Guth (1981) Guth, A. H. 1981, Physical Review D, 23, 347
  * Guth & Pi (1982) Guth, A. H. & Pi, S.-Y. 1982, Physical Review Letters, 49, 1110
  * Harrison (1970) Harrison, E. R. 1970, Physical Review D, 1, 2726
  * Jaffe et al. (2001) Jaffe, A. H., Ade, P. A. R., Balbi, A., et al. 2001, Phys. Rev. Lett., 86, 3475
  * Jones (2017) Jones, B. J. T. 2017, Precision Cosmology: The First Half Million Years (Cambridge University Press)
  * Kono et al. (2020) Kono, K. T., Takeuchi, T. T., Cooray, S., Nishizawa, A. J., & Murakami, K. 2020, arXiv e-prints, arXiv:2006.02905
  * Mahalanobis (1936) Mahalanobis, P. C. 1936, in Proceedings National Institute of Science, India, Vol. 2, 49–55
  * Masi (2002) Masi, S. 2002, Progress in Particle and Nuclear Physics, 48, 243
  * Matsubara (2010) Matsubara, T. 2010, Physical Review D, 81, 083505
  * Mecke et al. (1994) Mecke, K. R., Buchert, T., & Wagner, H. 1994, Astronomy & Astrophysics, 288, 697
  * Munkres (1984) Munkres, J. 1984, Elements of Algebraic Topology, Advanced book classics (Perseus Books)
  * Park et al. (2013) Park, C., Pranav, P., Chingangbam, P., et al. 2013, Journal of Korean Astronomical Society, 46, 125
  * Peebles & Yu (1970) Peebles, P. J. E. & Yu, J. T. 1970, Astrophysical Journal, 162, 815
  * Planck Collaboration et al. (2016a) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016a, Astron. & Astrophysics, 594, A16
  * Planck Collaboration et al. (2016b) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016b, Astron. & Astrophysics, 594, A17
  * Planck Collaboration et al. (2016c) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016c, Astron. & Astrophysics, 594, A12
  * Planck Collaboration et al. (2020a) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020a, Astron. & Astrophysics, 641, A1
  * Planck Collaboration et al. (2020b) Planck Collaboration, Akrami, Y., Andersen, K. J., et al. 2020b, Astron. & Astrophysics, 643, A42
  * Planck Collaboration et al. (2020c) Planck Collaboration, Akrami, Y., Ashdown, M., et al. 2020c, Astron. & Astrophysics, 641, A7
  * Planck Collaboration et al. (2020d) Planck Collaboration, Akrami, Y., Ashdown, M., et al. 2020d, Astron. & Astrophysics, 641, A4
  * Pogosyan et al. (2009) Pogosyan, D., Gay, C., & Pichon, C. 2009, Physical Review D, 80, 081301
  * Pranav (2015) Pranav, P. 2015, Persistent Holes in the Universe: A Hierarchical Topology of the Cosmic Mass Distribution (University of Groningen)
  * Pranav et al. (2019a) Pranav, P., Adler, R. J., Buchert, T., et al. 2019a, Astron. & Astrophysics, 627, A163
  * Pranav et al. (2017) Pranav, P., Edelsbrunner, H., van de Weygaert, R., et al. 2017, Mon. Not. Royal Astro. Soc., 465, 4281
  * Pranav et al. (2019b) Pranav, P., van de Weygaert, R., Vegter, G., et al. 2019b, Mon. Not. Royal Astro. Soc., 485, 4167
  * Roukema et al. (2004) Roukema, B. F., Lew, B., Cechowska, M., Marecki, A., & Bajtlik, S. 2004, Astron. & Astrophysics, 423, 821
  * Ryden (2003) Ryden, B. 2003, Introduction to Cosmology (Addison-Wesley)
  * Sahni et al. (1998) Sahni, V., Sathyprakash, B., & Shandarin, S. 1998, Astrophysical Journal, 507, L109
  * Schmalzing & Gorski (1998) Schmalzing, J. & Gorski, K. M. 1998, Mon. Not. Royal Astro. Soc., 297, 355
  * Schwarz et al. (2016) Schwarz, D. J., Copi, C. J., Huterer, D., & Starkman, G. D. 2016, Classical and Quantum Gravity, 33, 184001
  * Shivashankar et al. (2016) Shivashankar, N., Pranav, P., Natarajan, V., et al. 2016, IEEE Trans. Vis. Comput. Graph., 22, 1745
  * Starobinsky (1982) Starobinsky, A. A. 1982, Physics Letters B, 117, 175
  * Telschow et al. (2019) Telschow, F., Schwartzman, A., Cheng, D., & Pranav, P. 2019, arXiv e-prints, arXiv:1908.02493
  * Tukey (1975) Tukey, J. W. 1975, in Proceedings of the 1974 international congress of mathematicians, Vol. 2, 523–531
  * van de Weygaert et al. (2011) van de Weygaert, R., Pranav, P., Jones, B. J. T., et al. 2011, ArXiv e-prints, arXiv:1110.5528
  * van de Weygaert et al. (2011) van de Weygaert, R., Vegter, G., Edelsbrunner, H., et al. 2011, Transactions on Computational Science, 14, 60
  * Verde et al. (2001) Verde, L., Jimenez, R., Kamionkowski, M., & Matarrese, S. 2001, Mon. Not. Royal Astro. Soc., 325, 412
  * Vilenkin & Shellard (2000) Vilenkin, A. & Shellard, E. P. S. 2000, Cosmic Strings and Other Topological Defects (Cambridge University Press)
  * Wilding et al. (2020) Wilding, G., Nevenzeel, K., van de Weygaert, R., et al. 2020, arXiv e-prints, arXiv:2011.12851
  * Xu et al. (2019) Xu, X., Cisewski-Kehe, J., Green, S., & Nagai, D. 2019, Astronomy and Computing, 27, 34–52

##  Appendix A Data and methods

In this Section, we briefly describe the dataset used in the experiments. We
also present a short account of the computational pipeline, referring the
reader to Pranav et al. (2019a) for details.

###  A.1 Datasets

Over a period of time, the Planck team has invested significant effort in
understanding and calibrating the source of noise as well as systematics that
affect observations. This has led to the release of a series of datasets over
a period of time, with successive data releases achieving better calibration
and accuracy. In Pranav et al. (2019a), we performed our experiments on the
second data release (DR2 hereafter). which is one of the intermediate data
releases. In this paper, we concentrate on two data releases simultaneously,
which we briefly describe below.

_Planck 2020_ Data Release 4 NPIPE dataset. The primary results are based on
the observational maps and simulations based on the NPIPE (Planck
Collaboration et al. 2020b) analysis pipeline (hereafter just NPIPE), which
employs the SEVEM component separation technique. The NPIPE dataset is the
final data release from the Planck team, and the pipeline is the most evolved
and sophisticated of all data generation pipelines. It combines some of the
most powerful features of the separate LFI and HFI analysis piplelines,
resulting in frequency and component maps that have lower levels of noise and
systematics at essentially all angular scales (c.f. Planck Collaboration et
al. (2020b)). The observational maps are accompanied by a suite of $600$
simulations modeled as isotropic, homogeneous Gaussian random fields.

_Planck 2018_ Data Release 3 FFP10 dataset. For comparison, we also analyze
the temperature maps from _Planck 2018_ full-mission data release (DR3
hereafter), which is the third data release by the Planck team, and employs
the SMICA component separation technique. The noise and systematics in this
release are better constrained than in the previous data releases. The
component separated observational maps are accompanied by a suite of $300$
simulations generated using the _Full Focal Plane Plane_ pipeline (hereafter
referred to as the FFP10 simulations) (Planck Collaboration et al. 2020d),
also based on the standard model.

###  A.2 Computational pipeline

(a)
Figure 6: A visualization of the mask employed for the analyses in this paper.
It is the common mask released in Data release 3 (2018), and masks our galaxy
as well as other bright foreground sources, including point sources.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 7: A visualization of the degraded, smoothed and masked observational maps employed in this paper. The maps are degraded at ${N}=2048,1024,512,256,128,64,32,16$, and smoothed at $FWHM=5^{\prime},10^{\prime},20^{\prime},40^{\prime},80^{\prime},160^{\prime},320^{\prime},640^{\prime}$ respectively. The mask is degraded and smoothed at the same scales as the map, and subsequently thresholded at $0.9$ for re-binarization. Resolution | FWHM(arcmin) | % unmasked
---|---|---
2048 | 5 | 77.9
1024 | 10 | 76.9
512 | 20 | 75.6
256 | 40 | 74.7
128 | 80 | 73.6
64 | 160 | 71.3
32 | 320 | 68.8
16 | 640 | 64.5
Table 2: Percentage of sky area covered by the unmasked regions for the
various degraded resolutions and smoothing scales. The re-binarization
threshold is for mask is set at $0.9$.

The computational pipeline is composed of the steps described below in brief;
see Pranav et al. (2019a) for a substantially detailed description.

Pre-processing. We perform a range of pre-processing steps to the maps using
the HealPix package (Górski et al. 2005), which is also the format of the
initial data. First, the CMB and the noise maps are added pixel-wise for each
realization of the simulations. Subsequently, the observational and noise-
added simulation maps are degraded and smoothed at a range of scales for a
scale-dependent analysis. The mask, presented in Figure 6, which is a binary
map, is also degraded and smoothed at the same resolution as the CMB maps, and
subsequently thresholded at a value of $0.9$ to make it binary again. Table 2
presents the percentage of sky covered by the unmasked region for the degraded
resolution and associated smoothing scales analyzed in this paper. Figure 7
presents a visualization of the degraded and smoothed maps, after applying the
mask. Next, we compute the mean and standard deviation from the unmaksed
pixels, for each realization and resolution, and rescale the maps, pixel-wise,
by the standard deviation after subtracting the mean. As a final step, we
assign the value infinity to the masked pixels (numerically simulated by a
very large number).

Triangulation. As a first step, we project pixels of the maps projected onto
${{\mathbb{S}}}^{2}$, and we triangulate this set of points in
${{\mathbb{R}}}^{3}$. Taking the convex hull of this triangulation produces a
triangulation of the point-set on ${{\mathbb{S}}}^{2}$. This triangulation
consists of $V=12{N}^{2}$ vertices, $3V-6$ edges, and $2V-4$ triangles, where
${N}$ is the resolution parameter in HealPix format. It is the input to all
downstream computations, and represents the temperature field,
$f\colon{{\mathbb{S}}}^{2}\to{{\mathbb{R}}}$, by storing the temperature value
at each vertex (Pranav et al. 2019a). We assume a piece-wise linear
interpolation along higher dimensional simplices.

Upper-star filtration. Given the triangulation $K$ constructed in the previous
step, we order its simplices 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 vertices of ${\sigma}$. Any ordering that satisfies
(i) and (ii) is called 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 of $f$.

Persistence computation. We construct a boundary matrix from the upper-star
filtration of $K$. Writing ${\sigma}_{1},{\sigma}_{2},\ldots,{\sigma}_{n}$ for
the sorted simplices of the upper-star filtration, the boundary matrix
$\partial[1..n,1..n]$ is defined by $\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. Computing the
persistence birth death pair from this ordered boundary matrix involves
reduction of the columns to the _lowest-j_ form. We resort to reduction from
left to right, and a column of the matrix is reduced if if it is zero or its
lowest $1$ has only $0$s in the same row to its left. Each column with a
unique _lowest-j_ contributes to the birth death pair of the persistence
diagram, where the index of the birth and death simplices are precisely the
row and column indices of the _lowest-j_. Ranks of homology groups relative to
the mask are inferred from the persistence diagrams by setting the vertices
belonging to the mask at $+\infty$ (Pranav et al. 2019a):

${b}_{0}$ $\displaystyle=\\#\\{[b,d)\in
Dgm_{0}({{\mathbb{E}}}\cup{{\mathsf{M}}})\mid+\infty>b\geq\nu>d\\};$ (5)
$\displaystyle{{b}_{1}}$ $\displaystyle=\\#\\{[b,d)\in
Dgm_{0}({{\mathbb{E}}}\cup{{\mathsf{M}}})\mid+\infty=b>d\geq\nu\\}$
$\displaystyle+\\#\\{[b,d)\in
Dgm_{1}({{\mathbb{E}}}\cup{{\mathsf{M}}})\mid+\infty>b\geq\nu>d\\};$
$\displaystyle{{b}_{2}}$ $\displaystyle=\\#\\{[b,d)\in
Dgm_{1}({{\mathbb{E}}}\cup{{\mathsf{M}}})\mid+\infty=b>d\geq\nu\\}$
$\displaystyle+\\#\\{[b,d)\in
Dgm_{2}({{\mathbb{E}}}\cup{{\mathsf{M}}})\mid+\infty>b\geq\nu>d\\}.$

Tests for determining statistical significance. Our main aim is comparing the
observational CMB maps with the null-hypothesis simulation maps, which assume
isotropy, homogeneity and Gaussianity. 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 $G$. Let $\mathbf{y}\in\mathbb{R}^{m}$ be another sample point,
assumed to be drawn from a distribution $F$. 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 employ two different statistical tests, which we briefly
describe below.

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

$p_{\rm Mahal}(\mathbf{y})=P[\chi^{2}_{m}>d^{2}_{\rm Mahal}(\mathbf{y})].$ (6)

_Tukey depth._ The second method is a non-parametric test based on the _Tukey
depth_ (Tukey 1975). It a general metric for identifying outliers in a
flexible manner and in a non-parametric setting, making no assumptions on the
structure of $F$ and $G$. 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}$. Points that have the same depth constitute a
non-parametric estimate of the isolevel contour of the distribution $F$. 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.
The $p$-value of $\mathbf{y}$ is computed as the proportion of points whose
depth is lower than that of $\mathbf{y}$:

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

##  Appendix B FFP10 results

In this section, we present the graphs and the table of $p$-values
corresponding to the _Planck 2018_ Data release 3 (DR3) dataset. We base the
significance of our results on the accompanying $300$ FFP10 simulations.

(a)

(b)

(c)

Figure 8: Significance of difference for ${{b}_{0}}$, ${{b}_{1}}$ and ${\sf
EC_{\rm rel}}$ for the FFP10 dataset. Relative homology
---
|  | Mahalanobis | Tukey Depth
Resolution | FWHM | ${{b}_{0}}$ | ${{b}_{1}}$ | ${\sf EC_{\rm rel}}$ | ${{b}_{0}}$ | ${{b}_{1}}$ | ${\sf EC_{\rm rel}}$
threshold = 0.90
2048 | 5 | 0.454 | 0.573 | 0.887 | 0.000 | 0.367 | 0.617
1024 | 10 | 0.728 | 0.532 | 0.690 | 0.580 | 0.000 | 0.000
512 | 20 | 0.180 | 0.641 | 0.416 | 0.000 | 0.407 | 0.000
256 | 40 | 0.398 | 0.420 | 0.516 | 0.417 | 0.000 | 0.000
128 | 80 | 0.075 | 0.295 | 0.112 | 0.000 | 0.000 | 0.000
64 | 160 | 0.442 | 0.200 | 0.811 | 0.483 | 0.000 | 0.720
32 | 320 | 0.032 | 0.192 | 0.140 | 0.000 | 0.353 | 0.000
16 | 640 | 0.853 | 0.001 | 0.000 | 0.870 | 0.000 | 0.000
summary | NA | 0.001 | 0.002 | 0.000 | 0.000 | 0.000 | 0.000
(a)
Table 3: Table displaying the two-tailed $p$-values for relative homology
obtained from parametric (Mahalanobis distance) and non-parametric (Tukey
depth) tests, for different resolutions and smoothing scales for the FFP10
dataset. The last entry is the $p$-value for the summary statistic computed
across all resolutions. Marked in boldface are $p$-values $0.05$ or smaller.