# Improved Galactic Foreground Removal for B-Modes Detection with Clustering
Methods

Giuseppe Puglisi1,2,3, Gueorgui Mihaylov4, Georgia V. Panopoulou5, Davide
Poletti6,7, Josquin Errard8, Paola A. Puglisi10, Giacomo Vianello11
1Computational Cosmology Center, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA
2Space Sciences Laboratory at University of California, 7 Gauss Way, Berkeley,
CA 94720
3Department of Physics, University of California, Berkeley, CA, USA 94720
4GSK and King’s College London, Department of Mathematics , Strand, London UK
5California Institute of Technology, MC350-17, 1200 East California Boulevard,
Pasadena, CA 91125, USA
6Universitá degli Studi di Milano-Bicocca, Piazza della Scienza 3, 20126
Milano, Italy
7 INFN - Sezione di Milano Bicocca, Piazza della Scienza 3, 20126 Milano,
Italy
8 AstroParticule et Cosmologie, Univ. Paris Diderot, CNRS/IN2P3, CEA/Irfu, Obs
de Paris, Sorbonne Paris Cité , France
9Royal Mail, 185 Farringdon Rd, London EC1A 1AA, UK
10Cape Analytics, 100 W Evelyn Ave UNIT 220, Mountain View, CA 94041, USA
E-mail<EMAIL_ADDRESS>Fellow

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Characterizing the sub-mm Galactic emission has become increasingly critical
especially in identifying and removing its polarized contribution from the one
emitted by the Cosmic Microwave Background (CMB). In this work, we present a
parametric foreground removal performed onto sub-patches identified in the
celestial sphere by means of spectral clustering. Our approach takes into
account efficiently both the geometrical affinity and the similarity induced
by the measurements and the accompanying errors. The optimal partition is then
used to parametrically separate the Galactic emission encoding thermal dust
and synchrotron from the CMB one applied on two nominal observations of
forthcoming experiments from the ground and from the space. Moreover, the
clustering is performed on tracers that are different from the data used for
component separation, e.g. the spectral index maps of dust and synchrotron.
Performing the parametric fit singularly on each of the clustering derived
regions results in an overall improvement: both controlling the bias and the
uncertainties in the CMB $B-$mode recovered maps. We finally apply this
technique using the map of the number of clouds along the line of sight,
$\mathcal{N}_{c}$, as estimated from HI emission data and perform parametric
fitting onto patches derived by clustering on this map. We show that adopting
the $\mathcal{N}_{c}$ map as a tracer for the patches related to the thermal
dust emission, results in reducing the $B-$mode residuals post-component
separation. The code is made publicly available .

###### keywords:

keyword1 – keyword2 – keyword3

††pubyear: 2021††pagerange: Improved Galactic Foreground Removal for B-Modes
Detection with Clustering Methods–D

##  1 Introduction

In the past few decades, the anisotropies of the Cosmic Microwave Background
(CMB) have been measured to constrain the parameters of the standard
cosmological model, the $\Lambda$-Cold-Dark-Matter ($\Lambda$CDM), with
unprecedented precision (Planck Collaboration et al., 2016d, 2020b). In recent
years, the polarization of the CMB has received significant attention. The
linear polarization of the CMB arises primarily from Thompson scattering of
photons with free electrons in the photon-baryon plasma of the epoch of
recombination. The linearly polarized emission anisotropies can be decomposed
into a curl-less field, commonly referred to as the $E-$modes, and a
divergence-less one, known as $B-$modes (Seljak & Zaldarriaga, 1997; Hu &
White, 1997). The $E-$modes are linked to the primordial scalar perturbations,
whereas $B-$modes at the degree scale can be produced only by tensor
perturbations of the space-time metric emitted at the _inflationary_ era
(Guth, 1981; Starobinsky, 1982). $B-$modes are expected to be observed at
degree angular scales, making them an interesting scientific target to
validate several theoretical models describing the early universe as their
amplitude, commonly parametrized by the tensor-to-scalar ratio $r$, is
proportional to the energy scale when inflation occurred.

At smaller angular scales ($\sim$arcmin), $B-$modes can be sourced by the
direct distortion of $E-$modes due to gravitational lensing of the intervening
large scale structure. Lensing $B-$modes have been detected in the past years
by many ground-based experiments, (e.g. The Polarbear Collaboration et al.,
2017, 2019; Choi et al., 2020; Bianchini et al., 2020). On the other hand,
primordial B-modes have not been detected yet and the best constraints have
recently been set to $r<0.044,\,95\%$ confidence level, by Tristram et al.
(2021) combining data from the BICEP2/Keck Array and Planck experiments.

Together with the instrumental sensitivity, the major challenge to detect
primordial $B-$modes is the polarized foreground radiation emitted at the same
CMB frequencies ($\nu\sim 100$ GHz) from our own Galaxy. Electrons accelerated
by spiralling along the Galactic magnetic field lines emit synchrotron
emission mostly dominating the foreground polarized emission at low
frequencies, $\nu<70$ GHz. On the other hand, thermal dust emission arising
from grains aligned with the Galactic magnetic field dominates the high
frequency end of the sub-mm foreground radiation, $\nu>100$ GHz. Polarization
from other processes, (e.g. from molecular lines Puglisi et al., 2017), is
expected to be subdominant with respect to that of Galactic synchrotron and
dust emission (Planck Collaboration et al., 2016b, 2020a).

Since the Galactic components have different frequency dependences, they can
be disentangled by means of component separation or foreground cleaning
techniques (Planck Collaboration et al., 2014, 2016a, 2016c). On one hand,
_blind_ algorithms, (e.g. Internal Linear Combination) are usually aimed at
removing emission that is different than the CMB, exploiting either the
assumption of statistical independence for most of the sky components (see
Remazeilles et al., 2018; Delabrouille et al., 2003; Maino et al., 2007) or
the maximum entropy principle (Stolyarov et al., 2005). On the other hand,
_non-blind_ approaches e.g. internal template subtraction or parametric
fitting (Dunkley et al., 2009; Eriksen et al., 2008; Leach et al., 2008;
Stompor et al., 2008; Bobin et al., 2008; Hansen et al., 2006; Bennett et al.,
1992) are employed in order to distinguish several Galactic components. These
come at the cost of making assumptions to model the frequency scaling of all
the emission involved, commonly accounting for a very large number of free
parameters. Moreover, what accentuates the difficulties of foreground removal
is the fact that the properties of the interstellar medium (ISM) change not
only _spatially_ , i.e. in different locations in the sky (Planck
Collaboration et al., 2016b, 2020a) but also along the same line-of-sight
(Tassis & Pavlidou, 2015; Clark, 2018; Panopoulou & Lenz, 2020), eventually
leading to frequency dependencies that are more complex than typically assumed
(Chluba et al., 2017; Pelgrims et al., 2021; Mangilli et al., 2021;
Azzoni_2021). As a result, algorithms relying on parametric component
separation would tend to reconstruct CMB maps with large residual foreground
bias as the spatial variability of foregrounds is mismodelled by fitting a
single parameter across the whole observed sky.

A partition of the sky, where the parameters are fit independently on multiple
regions, results in a mitigation of the bias but might increase the
statistical uncertainties due to instrumental noise. This is mainly due to the
fact that the fit is performed on a smaller number of pixels, encoding a lower
signal-to-noise ratio (SNR) than the case with all the observed pixels.
Moreover, performing a pixel-by-pixel foreground removal is clearly unfeasible
for experiments involving large sky footprints and/or high resolution because
fitting for spectral indices in each sky pixel is a very costly process.

Recently, several techniques to encompass the spatial variability of the
foregrounds have been proposed by decomposing the map into wavelets or
needlets (Basak & Delabrouille, 2011; Remazeilles et al., 2018; Wagner-Carena
et al., 2019; Irfan et al., 2019), by partitioning the sky into regions
(clusters) according to the similarity of foreground properties and their
location in the sky (Khatri, 2019; Grumitt et al., 2020).

In this work, we do spectral clustering on a proxy of SED properties of the
Galactic foregrounds by considering the spectral parameters and the number of
clouds along the line of sight map. Once the clusters are defined used them
for component-separating multi-frequency observations (Stompor et al., 2008;
Errard et al., 2011; Stompor et al., 2016a; Alonso et al., 2017).

The approach is similar to the one shown in Grumitt et al. (2020) but we
choose a different _clustering method_ to optimally divide the sky. In fact,
Grumitt et al. (2020) used the _mean-shift_ algorithm (Krzanowski & Lai, 1988)
to locate overdensities in the context of image segmentation. They developed a
parametric bayesian component separation algorithm that makes use of a
clustering analysis to forecast the accuracy of component separation for a
LiteBIRD-like CMB space satellite mission (Sugai et al., 2020). They define
clusters in a 5-dimensional parameter space: 3-cartesian coordinates of pixels
of the sphere, and the spectral indices $\beta_{s}$ and $\beta_{d}$ used to
parametrize respectively the emission of synchrotron and dust. However, as
noted by the authors, more work is needed to identify the optimal sky
templates for this kind of study. An obvious area of improvement is to
optimize the definition of sky regions. The characteristics of clusters
obtained in Grumitt et al. (2020) were almost homogeneous across the sky:
showing approximately the same size and similar (convex) shapes. However, the
mean-shift cluster size is highly dependent on the choice of the bandwidth
parameter in the definition of the kernels to identify the clusters. Ideally,
the clustering should reflect the underlying spatial distribution of the
features of interest. More sophisticated algorithms for image segmentation are
available for this. In this work, we focus on improving upon the definition of
clusters on the sky.

We propose in this work a different algorithm to partition the sky via
_spectral clustering_ (Von Luxburg, 2007). In this specific application, the
foreground spectral parameter maps are defined over the whole celestial sphere
and can thus be mathematically treated as real-valued functions defined on the
real manifold $S^{2}$. In the context of unsupervised image segmentation,
spectral clustering has proved to be a powerful technique able to capture
objects with highly non-trivial geometric shapes (Zhang et al., 2018; Zelnik-
Manor & Perona, 2004). Our approach takes into account the degree of
similarity inferred both from geometric positions and from the measurements
performed in different points of the sky. Being based on an eigen-
decomposition, it is less affected by high-dimensionality issues. Furthermore,
we use the signal-to-noise content of the spectral parameter maps to allow the
pixel similarity to be informed by the intrinsic variability of a given
foreground parameter map.

The definition of clusters on a manifold can be divided into two procedures:
(i) building the pixel similarity accounting for the “distances” in a given
metric, (ii) finding the eigenspectrum decomposition for a given set of
features. In Sect. 2 we present the spectral clustering methodology and the
formalism to implement it on $S^{2}$. In Sect. 3, we present our
implementation of spectral clustering that makes use of the HEALPix sky
tesselation. In Sect. 4 we show two applications of how the patches defined
with spectral clustering can be used to improve the performance of parametric
component separation techniques. We finally discuss results and cosmological
implications in Sect. 5 and 6.

##  2 Spectral image segmentation

Image segmentation via spectral clustering is an efficient technique that
combines spatial similarity (or affinity) of pixels with the similarity based
on image-related characteristics of the pixels (Zhang et al., 2018; Zelnik-
Manor & Perona, 2004). There are profound geometric and analytical reasons for
this and are concisely discussed in the dedicated Appendix A.

Any image can be mathematically encoded by means of a real-valued (ideally
smooth) function $f$, which is defined on a portion of a surface and sampled
with a fixed resolution given by the pixel grid. For example a grey-scale
image is encoded by one function, whereas RGB images are encoded by three.

We summarize below the standard steps for a spectral clustering process:

  1. 1.

Translate an image into a graph by defining a suitable _adjacency_ condition.
A typical choice is that any pixel is connected to its nearest-neighbours. For
example, in correspondence to pixel $j$ neighbour of $i$, we can set to a non-
zero value the $ij$ element of the so called _adjacency_ matrix $A_{ij}$. Once
the connections among all pixels are defined via the weights on the adjacency
matrix we can build a graph by choosing a conventional order of the vertices
(pixels) $x_{i}$. The way the adjacency weights are assigned is usually based
on accounting for the geometric distance or the distance between the values
that $f$ takes in each pixel. An example of weights estimated from a distance
is commonly weighting the graph nodes by a Gaussian function of the distance
between two pixels.

  2. 2.

Compute the graph Laplacian matrix from the weighted adjacency matrix. We use
the symmetric random walk Laplacian:

$L_{sym}\equiv I-D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ (1)

where the $I$ is the identity matrix and $D_{i}=\sum_{j}A_{ij}$ is a diagonal
matrix called _degree_ matrix.

  3. 3.

Compute the eigenvectors of the graph Laplacian, $L_{sym}$ and select a subset
of $n$ eigenvectors, in correspondence to the ones that contribute the most to
the eigenspectrum of $L_{sym}$. E.g. for the specific case of 2D images, each
eigenvector can be visualized as images encoding large (small) spatial
variations depending whether the correspondent eigenvalue is small (large).
This process defines an _embedding_ of the graph associated to the image in
$\mathbb{R}^{n}$, i.e. each pixel $i$ in the image is characterized by $n$
features, i.e. the value of each eigenvector in the $i$-th pixel.

  4. 4.

Once the embedding in $\mathbb{R}^{n}$ is defined, the standard Euclidean
distance in $\mathbb{R}^{n}$ is evaluated for all the pair of pixels in order
to run a suitable _Agglomerative_ or _Divisive_ clustering algorithm.

In this work, we combine the spectral clustering and adapted it in the context
of quantities defined in the $S^{2}$ manifold. This is generally referred to
as _manifold learning_ : a discipline that combines statistical methods with
techniques developed in differential geometry. It is based on the assumption
that point clouds of multivariate $n-$dimensional variables are sampled on or
close to smooth compact sub-manifolds of $\mathbb{R}^{n}$ ( the existence of a
boundary is usually ignored). In Appendix A, we outline more in detail how
manifold learning can be combined with spectral clustering. However, we remark
here two relevant and non-standard aspects of the implementation of the method
proposed in this paper that make the geometric nature of spectral clustering
very relevant:

  * •

Maps from very wide astronomical surveys are sampled on the sphere, we
therefore implement the affinity between pixels accounting for the angular
distance on the $S^{2}$ sphere (presented in Subsect. 3.1).

  * •

We choose the partion as the one that minimizes the within- and between-
cluster variances. We present in Subsect. 3.2, how the variance depends on
several partition choices and eigenvector bands.

In order to better understand the details in the following sections, we
briefly summarize few propositions in the context of Graph Theory which will
be used to justify the choice for the adjacency in Sect. 3.

###### Proposition 1

The construction of the graph Laplacian matrix $L$ is known to approximate the
Laplace-Beltrami operator $\Delta_{g}$ on a Riemannian manifold embedded in
$\mathbb{R}^{n}$ and sampled in a grid of points (e.g. our set of pixels).
Therefore, the weighted adjacency matrices can be directly derived from the
integral kernel (or Green function) of the Laplacian operator in the sense of
Fredholm’s theory (Fredholm, 1903).

###### Proposition 2

The Laplace-Beltrami operator on a compact differentiable manifold $M$ has
discrete spectrum
$\lambda_{0}\leq\lambda_{1}\leq\lambda_{2}\leq...\leq\lambda_{n}$ and its
eigenfunctions $\varphi_{i}$ form an orthonormal basis of the space
$L^{2}(M)$. Analogously, every discrete approximation of a function on a
manifold can be expressed with respect to a basis of eigenvectors of the graph
Laplacian.

As a consequence of Proposition 1 and 2, estimating the eigenspectrum of
suitable graph Laplacian matrices sampled in a grid approximates the spectrum
(eigenvalues and eigenfunctions) of the Laplace-Beltrami operator.

Moreover, in defining the optimal image segmentation it is quite common to
define the _Dirichlet energy functional_ , i.e.

$E[f]\equiv\frac{1}{2}\int_{M}\parallel\nabla f(x)\parallel^{2}d\omega_{g}$

that measures the spatial variability (or _spatial frequency_) of a smooth
function on a manifold $M$. Thus, since the Dirichlet energy associated to the
eigenfunctions of the Laplace-Beltrami operator is
$E(\varphi_{i})=\frac{1}{2}\lambda_{i}$ with monotonically increasing
Dirichlet energy111A common example of increasing Dirichlet energy
eigenfunctions can be found among the eigenfunctions of the Laplace-Beltrami
operator in the sphere, i.e. the Spherical Harmonics, $Y_{\ell m}$. We observe
increasingly larger spatial variability of the $Y_{\ell m}$ in correspondence
to larger values of multipole number $\ell$, with fluctuations involving
smaller and smaller angular scales. (Proposition 2), the set of the first $n$
eigenfunctions provides an optimal embedding of $M$ in $\mathbb{R}^{n}$
corresponding to the minimum Dirichlet energy222Similarly, this properties
propagates to the eigenvectors of the Laplacian matrix (eq. (1)).. Generally,
this process can be also addressed by selecting a subset of eigenfunctions
(corresponding to a certain _Dirichlet energy band_) to embed $M$ in a
$n-$dimensional space.

The choice of the eigenfunctions related to a specific Dirichlet energy band
can impact the quality of the image segmentation. This is mainly due to the
fact that the clustering methodology is based on the choice of the
eigenvectors of the Laplacian matrix, to which corresponds a certain degree of
spatial variability, or _granularity_. We devote Sect. 3.2 to describe how
this applies to the case presented in this work and how the choice of the
energy band affects the characteristic size of the clusters.

The key ingredient in this framework is the definition of the adjacency, from
which the Laplacian matrix can be then derived via eq. (1). In the next
section, we describe extensively the adjacency adopted in this work.

##  3 Spectral Clustering on HEALPIX maps

In this section, we present an implementation of spectral clustering applied
on images defined on the full celestial sphere. We adopt maps following the
HEALPix scheme (Górski et al., 2005; Zonca et al., 2019).

###  3.1 Choice of the adjacency

Let us denote by $X=(\mathbf{x}_{1},\dots,\mathbf{x}_{n_{\rm pix}})$, a set of
normed-1 vectors (with $\mathbf{x}_{i}\in\mathbb{R}^{3}$), encoding the
coordinates of each HEALPix pixel. Given $X$, we can construct a matrix as
$\Theta=X^{T}X$, being a $n_{\rm pix}\times n_{\rm pix}$ matrix, with diagonal
equal to $1$, and off-diagonal elements encoding the scalar product estimated
pairwise on the columns of $X$. Notice that $\Theta$ encodes the cosines of
the scalar product as the columns of $X$ are normed-1.

We can therefore define the degree of affinity between pixels in a HEALPix map
from $\Theta$, as:

$A\propto\rm{exp}\left({-(1-\Theta)^{2}/2\sigma_{\rm pix}^{2}}\right),$ (2)

where $\sigma_{\rm pix}$333In this work, we adopted nside=32 which results in
$\sigma_{\rm pix}\sim 46^{\prime}$. is related to the pixel scale of the map.
By expanding eq. (2) in the Harmonic domain in terms of the Legendre
Polynomials $\mathcal{P}_{\ell}$, as

$A\equiv\sum_{\ell=0}^{+\infty}\frac{2\ell+1}{4\pi}e^{-\ell(\ell+1)\sigma_{\rm
pix}^{2}}\mathcal{P}_{\ell}(\Theta),$ (3)

we observe that the scale given by $\sigma_{\rm pix}^{2}$ sets a threshold
$\ell_{\rm max}$ to which we can stop the sum in eq. (3).

In Appendix A.4, we show how the functional form in eq. (3) derives from the
integral kernel (also known as _heat kernel_) of the Laplace-Beltrami
operator, $\Delta_{g}$, in $S^{2}$ Zhao & Song (2018). We show in Fig. 1 a row
of the affinity matrix $A$ as defined in eq. (2).

Figure 1: Adjacency as estimated in eq. (29) relative to the pixel centered at
$(\ell,b)=(0,30)$ deg Galactic coordinates. We recall the reader that $A_{ij}$
is a $n_{\rm pix}\times n_{\rm pix}$, so its rows and/or columns can be
visualized as HEALPix maps. Figure 2: Eigenvalues of Laplacian matrix
estimated as in eq. (1), for several choices of the scaling factor $\alpha$.

Figure 3: The first eigenvectors of $L$ corresponding to $\ell=0,\dots,3$ and
assuming no uncertainty deformation (i.e. $\alpha=0$). The colorbar ranges for
all the maps between $\pm 0.07$.

Figure 4: The first eigenvectors of $L$ with affinity deformation
$\alpha=0.25$. For sake of comparison with Fig. 3 we report also in this plot
the correspondence to $\ell=0,\dots,3$ although the analogy lose specific
meaning for $\alpha\neq 0$. The colorbar ranges for all the maps between $\pm
0.02$.

When clustering methodologies are employed, it is common to ask oneself how
the adjacency in eq. (2) is affected in presence of the uncertainties in the
measurements at a given pixel location. E.g. suppose that we have measurements
of the _modified black-body_ spectral index of dust emission, $\beta_{d}$, and
its uncertainties, $\sigma_{d}\equiv\sigma(\beta_{d})$. We require that two
pixels encoding statistically compatible values of $\beta_{d}$ would be more
easily associated together into a single cluster with respect to pixels with
incompatible values.

We therefore want to identify a way to combine the measurements and the
uncertainties into the heat-kernel adjacency (eq. 2) given realistic
measurements of foreground spectral parameters. For the following application,
we will use the map of the dust spectral index $\beta_{d}$ and the uncertainty
map $\sigma_{d}$ obtained from the _Planck_ data processed with the COMMANDER
separation algorithm Planck Collaboration et al. (2016b).

Intuitively, we can distort $\Theta$ defined in eq. (2) with a certain weight
given by the measurements and the uncertainty of the parameter at a given
pixel location. Given the state of art data on Galactic foregrounds in the
sub-mm regime, we consider the measurements in each pixel $i$ to be normally
distributed around $\beta_{d,i}$ with a width given by its uncertainties,
$\sigma_{d,i}$ (Planck Collaboration et al., 2016b). We therefore generate a
mock sample of 100 Gaussian random numbers given the features
$(\beta_{d,p},\sigma_{d,p})$, and we estimate the two-sample Kolmogorov-
Smirnov (KS) test onto pairs of different pixels $i$ and $j$. This allows us
to test the null hypothesis that the distribution of $\beta_{d}^{i}$ values is
drawn from the same underlying distribution as those of $\beta_{d}^{j}$ within
an assumed confidence level. We repeat the test $100$ times, evaluate the
statistical quantile $Q_{ij}$ of each KS test, and set the median value all
the KS tests, $Q_{ij}^{med}$ as the adjacency weight evaluated for the $ij$
pixel pair. Thus, if $Q_{ij}^{med}$ is small enough to reject the null
hypothesis at a $2\sigma$ confidence level, the connection is weighted with a
very low weight and it is unlikely that the two pixels will be associated into
the same cluster. On the other hand, when $Q_{ij}^{med}\sim 1$ we cannot
reject the null hypothesis and the two pixels can be associated, as they are
_parallel_ in this metric.

In order to combine the adjacency from the KS test and the one defined in eq.
(2), we can treat the KS quantiles as cosine of angles, so that a
straightforward deformation of $\Theta$ can be:

$\Theta^{\prime}=\cos\left(\arccos(\Theta)+\alpha\left(1-Q^{med}\right)\frac{\pi}{2}\right),$
(4)

where $0<\alpha\leq 0.5$ is a _scaling factor_ parameter that weighs the
relative contribution of the KS quantile similarity with respect to the
geometrical one (the proximity of two pixels on the sky). Values of
$\alpha>0.5$ result into too large distortions that break the metric
properties of eq. (4) and make the Laplacian a singular matrix (see next
Subsect. 3.2). The updated cosine matrix $\Theta^{\prime}$ is then finally
inserted into eq. (2) to evaluate the adjacency weights obtained from this
distortion.

Given the definition of adjacency (eq. (2) or (4)), we can therefore estimate
the Laplacian matrix as in eq. (1) and estimate its first $N_{eigen}$
eigenpairs, related to the smallest eigenvalues. We remark that we do not
factorize the whole Laplacian matrix, since we are interested in a subset of
eigenvectors. We thus approximate the eigenpairs by means of the so called
_Ritz approximation_ , a technique which has already been exploited in
previous literature (Szydlarski, M. et al., 2014; Puglisi, G. et al., 2018) to
approximate very well the exact eigenpairs of a matrix.

In Fig. 2, we show the first 100 Ritz eigenvalues $\lambda$ of $L$ for
different choices of $\alpha$. We firstly focus on the $\alpha=0$ case. As
already mentioned above, the functional form of the adjacency in eq. (2)
descends from the integral kernel of the Laplacian operator in $S^{2}$, which
specifically coincides to the _angular momentum_ operator $\hat{L}^{2}$ in
quantum mechanics (see Appendix A.4 for further details). We thus expect the
eigenvectors to be exactly the eigenfunctions of $\hat{L}^{2}$, i.e. the
Spherical Harmonics and the multiplicity of the eigenvalues to respect the
same algebraic multiplicity: for a given multipole number $\ell$, its
multiplicity goes as $2\ell+1$, in correspondence of the $m$ azimuthal
multipole number.

In fact, we notice that for the $\alpha=0$ case (solid blue) line, the
algebraic multiplicity grows exactly with $2\ell+1$. Furthermore, the
eigenvectors shown in Fig. 3 perfectly resemble the spherical harmonics. Both
the eigenvalue degeneracy as well as the morphology of the Laplacian
eigenvectors represent a remarkable validation test for the overall
implementation described above.

On the other hand, the distortion introduced by $\alpha\neq 0$ breaks the
Laplacian eigenfunction algebraic multiplicity, so that the eigenspectrum
becomes gradually strictly monotonic for higher values of $\alpha$. As a
consequence, the eigenvectors as well deviate from the spherical harmonics,
being more and more weighted by the measurement uncertainties. In our specific
case involving Galactic emission, we indeed observe several anisotropies
corresponding to the Galactic plane as shown in Fig. 4 for the case
$\alpha=0.25$.

Algorithm 1 Spectral Clustering Optimization of a field $X$ with uncertainty
$\sigma(X)$

1:procedure Spectral Clustering ($X,\sigma(X)$)

2: Parameter Affinity, $Q$ Initialization with ($X,\sigma(X)$) data

3: for $\alpha$ in $[0,\dots,\,\alpha_{max}]$ do

4: Build Affinity matrix with $\alpha$ distortion (eq. (4) )

5: Build Laplacian matrix, $L$ (eq. (1))

6: Estimate Ritz eigenpairs

7: Select an eigenvector-band, $U$

8: Estimate the euclidean affinity $E$ from the columns of $U$

9: for $\delta$ in $[0,\dots,\delta_{max}]$ do

10: Run Agglomerative clustering

11: Estimate $W(\alpha,\delta)$

12: end for

13: end for

14: Find local minima of $W$, $(\alpha_{*},\delta_{*})$

15:end procedure

###  3.2 Identifying the optimal partition

In this subsection, we aim at discussing how the clusters are estimated and
how to identify the range of partition optimality. The overall spectral
clustering algorithm described in the sections above is summarized in
Algorithm 1. In particular, we employ the implementation of agglomerative
clustering publicly available in the Scikit-Learn python
package444https://scikit-
learn.org/stable/modules/generated/sklearn.cluster.AgglomerativeClustering.html
to estimate the clusters.

Figure 5: Variance $W$ as in eq. (6) in the $\alpha-\delta$ space, evaluated
in clusters constructed with the 3 eigenvector bands as defined in the main
text. For the analysis shown in Sect. 4, we adopted the $\ell_{b3}$ being the
band that accounts for the eigenvectors encoding most of the angular scales to
be clustered.

As shown in Fig. 3, the first eigenvalues correspond to eigenvectors with
longer oscillation in the map, and vice versa. Furthermore, we observe that
the spectrum tends to saturate to a maximum value $\lambda_{max}$ after about
$\sim 100$ values (see Fig. 2). Since we want to identify the smallest
eigenvalues encoding most of the bulk of the Laplacian eigenspectrum, we
consider the first $n_{eig}=256$ Ritz eigenpairs and therefore build the
matrix $U\in\mathbb{R}^{n_{pix}\times n_{eig}}$ containing $n_{eig}$
eigenvectors as columns.

We remark here that the eigenvector matrix $U$ is critical in the context of
spectral clustering because it provides the embedding of $S^{2}$ into
$\mathbb{R}^{n_{eig}}$ space spanned by the Laplacian eigenvectors. Moreover,
it helps in reducing the dimensionality of our problem from $n_{pix}$
features, i.e. the number of features we have in a map, to $n_{eig}$ and in
avoiding the _curse of high dimensionality_ which clustering methodologies
might frequently incur.

We then consider each row of $U$, $(u_{i})_{i=1,\dots,n_{pix}}$ as our data
vector and we perform agglomerative clustering with $n_{pix}$ samples and
$n_{eig}$ features. We measure the pairwise similarity, $E_{ij}$, between two
vectors $u_{i}$ and $u_{j}$ (representing the $i$-th and $j$-th pixels) by
estimating the euclidean distance between them. The two pixels are associated
to a group if $E_{ij}<\delta$, with $\delta$ being a free parameter, commonly
referred as the _linkage distance threshold_. Small values of $\delta$ tend to
typically produce small clusters, and vice versa. This can be intuitively
understood as follows: the more constraining the distance threshold, the
harder it is to associate two pixels together into a cluster.

We notice that for several values of $\alpha$ the distance matrix range
changes in such a way that we get smaller eigenvector similarities for larger
values of $\alpha$. This is a consequence of the fact that the deformation of
adjacency with KS weights (eq. 4) results into a compression of the Laplacian
eigenspectrum to smaller condition numbers. In fact, we notice a trend in the
Laplacian eigenspectrum to i) lose the typical $2\ell+1$ degeneracy for
$\alpha\neq 0$, ii) flatten to smaller condition number, iii) to become
degenerate to a small value numerically close to zero. As a results the
convergence of the Laplacian Ritz eigenpairs estimation time increases for
$\alpha>0.3$. 555This is the reason why we set the maximum value of $\alpha$
to $0.3$.

Furthermore, the eigenvector distance for several values of $\alpha$ needs to
be normalized in such a way that we use use the same range in $\delta$ to
compare the variances estimated from each run of the agglomerative clustering.
We therefore rescale the eigenvector distance $E_{\alpha}$ evaluated for
different choices of $\alpha$ as it follows:

$E^{\prime}_{\alpha}=\frac{E_{\alpha}/M-m}{1-m},$ (5)

with $M=E_{{\rm max},\alpha}$ and $m=E_{{\rm min},\alpha}$. We note that
$E^{\prime}$ ranges from $E^{\prime}_{\rm max}=1$ to $E^{\prime}_{\rm
min}=Cm$, with $C=(1-M)/M\sim{\rm constant}$ for different values of $\alpha$.
The rescaling in eq. (5) can be seen as a sort of _rigid_ scaling of the
distance, making the distance matrices for different choices of $\delta$ and
$\alpha$ to be comparable without loosing too much information. Notice that
eq. (5) scales the maximum of the distance to 1 but we make sure that it does
not affect the shape of distance distribution. We also make sure that the
ratio of minimum distances for different choices of $\alpha$ is the same as
before the rescaling in eq. (5).

Furthermore, we consider 3 bands to identify the eigenvectors with major
contribution to the total Dirichlet energy budget, i.e. the spatial
variability, in the domain partition:

  * •

$\ell_{b1}$, including the first 100 eigenvectors, corresponding to
$1\lesssim\ell\lesssim 10$,

  * •

$\ell_{b2}$, including vectors between the 90-th and 150-th eigenvectors,
corresponding to $9\lesssim\ell\lesssim 12$,

  * •

$\ell_{b3}$, sampling the whole eigenspectrum every other 5 vectors up to the
175-th one, corresponding to $1\lesssim\ell\lesssim 14$.

We create a grid with different points evaluated in the $\alpha-\delta$
parameter space, with $\alpha\in[0,0.5]$, $\delta\in[0.05,0.3]$, and estimate
the cluster variance, $W$, defined as

$W\equiv\sqrt{V_{w}^{2}+V_{b}^{2}},$ (6)

with $V_{w}$ being the within-cluster variance,

$V_{w}\equiv\sum_{k=1}^{K}\sum_{u_{k}\in C_{k}}(u_{k}-\mu_{k})^{2},$

$V_{b}$ the between cluster variance,

$V_{b}\equiv\sum_{k\neq k^{\prime}}(\mu_{k^{\prime}}-\mu_{k})^{2},$

$K$ the total number of clusters and $\mu_{k}$ the centroid of cluster
$C_{k}$.

Fig. 5 shows $W$ defined in $\alpha-\delta$ space evaluated for the 3
different eigenvector bands. The optimal parameters $(\alpha_{*},\delta_{*})$
correspond to local minima of $W$. Notice that the range of optimality changes
for different eigenvector bands as different scales are encoded in each band.
We emphasize that similar optimal ranges are found among the three bands,
indicating that the clustering is stable with respect to the choice of the
band.

Since the $\ell_{b3}$ band encodes scales from most of the Laplacian
eigenspectrum, we adopt this band as a baseline for the applications presented
in the following sections.

Finally, we notice that the optimality, shown in Fig. 5, does not correspond
to a very localized region in the parameter space spanned by $\alpha$ and
$\delta$. Rather it shows as a narrow vertical band at around a certain value
of $\delta$ for which several choices of $\alpha$ can be equally optimal. We
find that intermediate values of $\alpha\sim 0.15$ lead to a balanced trade-
off between defining the clustering given the features of spectral parameters
and the intrinsic adjacency in $S^{2}$. We refer the reader to Appendix C for
further details.

The estimation of the Laplacian adjacency and eigenpairs and the optimization
of spectral clustering has been performed on 300 Cori KNL nodes of the NERSC
Supercomputing facility666https://www.nersc.gov/. It takes $\sim 4000$ cpu-
hours to execute the overall procedure outlined in Algorithm 1 with maps at
nside=32.

##  4 Applications

We present two applications of our spectral clustering algorithm for the
purpose of performing parametric component separation. First, in Sect.4.1, we
apply the spectral clustering algorithm to derive a partition of the sky from
publicly available maps of synchrotron and dust spectral parameters. Once the
patches are defined we perform a parametric component separation on each patch
and assess the post-foreground cleaning residuals in the CMB $B-$mode maps.
Secondly, in Sect. 4.2, we identify the patches by means of an ancillary
dataset: the number of HI clouds along the line of sight, $\mathcal{N}_{c}$,
presented in Panopoulou & Lenz (2020). This dataset is an independent tracer
of the foreground dust emission. We then perform component separation within
these patches and estimate the quality of the reconstruction.

###  4.1 Clustering applied on Galactic foreground spectral parameters

Table 1: Nominal specifics of CMB experiments to run component separation with FGBuster.  | Frequency [GHz] | Sensitivity $\left[\mu\mathrm{K\,arcmin}\right]$ | FWHM $[{\rm arcmin}]$
---|---|---|---
SO-SAT | $27,39,93$ | $35,21,2.6$ | $91,63,30$
$145,220,270$ | $3.3,6.3,16$ | $17,11,9$
LiteBIRD | $40,50,60,68,78$ | $37.5,24,19.9,16.2,13.5,$ | $69,56,48,43,39$
$89,100,119,140,166$ | $11.7,9.2,7.6,5.9,6.5$ | $35,29,25,23,21$
$195,235,280,337,402$ | $5.8,7.7,13.2,19.5,37.5$ | $20,19,24,20,17$
Figure 6: Variance surfaces estimated for several choices of $\alpha$ and
$\delta$, for (left) $\beta_{s}$, (middle) $\beta_{d}$, (right) $T_{d}$. Note
that regions of optimality correspond to local minima (lighter regions) in the
surfaces.

Figure 7: Cluster defined patches obtained from the parameters in the d1s1
model (each color denotes pixels belonging to different clusters). We show
Clusters estimated with $f=0$ (top row) and with $f=1$ (bottom row). Figure 8:
(left column) Foreground parameters from the PySM d1s1 model, (middle column)
parameter maps estimated with FGBuster on each clustering region (see Fig. 7)
using the LiteBIRD frequency channels. (right column) Relative residuals on
the estimated maps with FGBuster and the input ones from PySM.

We perform parametric component separation by fitting for three components:
CMB, dust and synchrotron emission. We consider two different combinations of
data, to be representative of the forthcoming experiments aimed at observing
primordial $B-$modes: the Simons Observatory Small Aperture Telescope (SO-SAT,
Simons Observatory Collaboration et al., 2019) from the ground and the
LiteBIRD space satellite (Sugai et al., 2020). The component separation is
performed respectively on $27\%$ and $60\%$ of the sky (see Fig. 9) with
frequency channels encoding the specifics shown in Table 1. Notice that for
SO-SAT we assume conservatively the _baseline_ configuration as described in
Simons Observatory Collaboration et al. (2019).

We use the Python Sky Model777https://pysm3.readthedocs.io/ (Thorne et al.,
2017, PySM) to simulate full-sky polarized emission of thermal dust, CMB and
synchrotron888The contribution in polarization from AME, CO and free-free is
expected to be very small and it is neglected in the following.. The
synchrotron radiation is commonly parametrized as a power law

$I_{\nu,synch}\propto\nu^{\beta_{s}},$

whereas the thermal dust emission is described by a modified black-body, i.e.

$I_{\nu,dust}\propto\nu^{\beta_{d}}B_{\nu}(T_{d}),$

with $T_{d}$ being the black-body dust temperature.

Among the different models available in PySM, we choose as a reference the
d1s1 model, which accounts for spatial variation in the synchrotron and dust
spectral indices, $\beta_{s},\beta_{d}$ and in $T_{d}$.

The synchrotron spectral index of the s1 model is the Miville-Deschênes et al.
(2008, model 5), estimated by combining the Haslam et al. (1982) map at 408
MHz and the map at 23 GHz from WMAP (Hinshaw et al., 2009). The $\beta_{s}$
map and its uncertainties $\sigma(\beta_{s})$ are available
online999https://lambda.gsfc.nasa.gov/ at $5$ deg resolution.

The d1 model is derived from the _Planck_ data
products101010https://pla.esac.esa.int released in Planck Collaboration et al.
(2016b), with $\beta_{d}$ and $T_{d}$ parameters obtained at about $1$ degree
resolution with the COMMANDER component separation algorithm Planck
Collaboration et al. (2016b).

When multiple frequency channels are simulated with the d1s1 model, algorithms
relying on parametric component separation would tend to reconstruct CMB maps
with large residual foreground bias (commonly referred as _systematic bias_)
if a single constant spectral parameter is fitted for each Galactic component
across the whole observed sky, as the spatial variability of foregrounds is
not taken into account in the fit.

On the other hand, the instrumental noise is responsible for the _statistical_
uncertainties which can be commonly assessed by running component separation
on frequency maps encoding different Monte-Carlo (MC) realizations of
instrumental noise and astrophysical signal.

A partition of the sky where the parameters are fit independently on multiple
regions might mitigate the systematic bias but increase the statistical
uncertainties as the fit is then performed on a smaller number of pixels,
hence encoding a lower signal-to-noise ratio (SNR) than the case with all the
observed pixels (Errard & Stompor, 2019).

Once the set of parameters is defined together with their associated
uncertainties, we can run the spectral clustering optimization procedure as
outlined in Sect. 3. The variance planes for $\beta_{s},\beta_{d},T_{d}$ are
shown in Fig. 6. We show the optimal patches selected in correspondence of
local minimum of the variance in Fig. 7.

We estimate the median value of the number of pixels included in each cluster
and we find about $4$ pixels per cluster for $\beta_{d}$, $8$ for $T_{d}$ and
$16$ for $\beta_{s}$, resulting in a typical angular sizes for the cluster
range from $\sim 3.5,\,7,\,14$ deg respectively for $\beta_{d}$, $T_{d}$, and
$\beta_{s}$. The different morphologies and sizes, as well as the different
number of clusters, result from a trade-off between the intrinsic variability
due to the astrophysical emission, the resolution of each map (see left column
of Fig.8) and the SNR. Specifically, the clustering algorithm finds it harder
to agglomerate pixels related to very high SNR regions (e.g. at low Galactic
latitudes) with respect to regions with larger uncertainties. In fact, this
trend is observed in all the regions: smaller patches where the SNR is high,
larger patches where it is low.

Figure 9:  Masks used for component separation and for the power spectra
estimation with NaMaster (Alonso et al., 2019) for LiteBIRD (top) and SO-SAT
(bottom) encoding respectively $60\%$ and $27\%$ of the sky.

We then implement a new functionality in
FGBuster111111https://fgbuster.github.io/fgbuster(Stompor et al., 2016b), a
package for parametric component separation, aimed at fitting parameters
within the set of multiple partitions of the sky. The emission model proposed
here is slightly different than in a common parametric setting: each spectral
parameter is fitted using data within a given sky region, which can have any
shape. For example, a single $\beta_{s}$ can be fitted within each region
shown in Fig. 7 (left), while at the same time $\beta_{d}$ is fitted on a
different set of regions (see middle panel, Fig. 7). The method is implemented
within the FGBuster framework and further details of the implementation will
be presented in two companion papers Errard et al. (in prep.) and Poletti et
al. (in prep).

For a pixel $p$, the model of emission is given by:

$\displaystyle\small m^{p,X}(\nu)$
$\displaystyle=A^{p,X}_{s}f_{s}(\nu,\beta_{s})+A^{p,X}_{d}f_{d}(\nu,\beta_{d},T_{d})+$
$\displaystyle A^{p,X}_{cmb}f_{cmb}(\nu)+n^{p,X}(\nu),$ (7)

i.e. a linear combination of the amplitudes of each astrophysical component
and the functions
$f_{s}(\nu,\beta_{s}),f_{d}(\nu,\beta_{d},T_{d}),f_{cmb}(\nu)$ encoding the
spectral dependence of synchrotron, dust, and CMB, a at a given frequency
$\nu$ as outlined in Sect. 4.1; $n(\nu)$ accounts for instrumental noise
coadded to the signal. Furthermore, we fit for the same spectral parameter
both for $X=Q,U$ maps.

We generate 20 Monte-Carlo realizations of noisy frequency maps with the
nominal specifications of LiteBIRD and for SO-SAT, see details in Table 1. The
frequency maps are produced with the PySM d1s1 model at nside=64 and include
both the input astrophysical polarized signal and the instrumental noise
(assumed to be only white noise). We show in Fig. 9 the nominal sky patches
where the parametric fit is performed. Notice that we choose a larger
$f_{sky}=0.27$ for SO-SAT with respect to the effective $f_{sky}=0.10$
reported in Simons Observatory Collaboration et al. (2019). The reason for
this choice is mainly motivated by the fact that we want to conservatively
assess the performance of this methodology onto larger regions of the sky
being thus more susceptible to larger Galactic residuals. This choice further
allows us to compare our results with the ones reported in Thorne et al.
(2019).

In Fig. 8 (middle column), we show the parameter maps estimated in each
cluster region with the LiteBIRD frequencies and we evaluate the relative
error of this estimate by considering the difference between the PySM
parameter maps and the parameters estimated with FGBuster (see right column in
Fig. 8).

Figure 10: Histograms evaluated from the relative error maps of
$\beta_{s},\beta_{d},T_{d}$ (right column of Fig. 8) respectively in solid
blue, orange and green. Shown as vertical dashed lines the median values,
namely $0.002,0.006,0.017$ respectively for $\beta_{s},\beta_{d},T_{d}$.

In Fig. 10, we show the distributions of relative errors for the three
parameters together with the median values shown as vertical dashed lines to
be $0.002,0.006,0.017$ respectively for $\beta_{s},\beta_{d},T_{d}$, with the
largest errors, $\leq 20\%$, found in $T_{d}$ estimates. This implies that the
component separation performed on the cluster-defined patches in Fig. 7 (top)
recovers faithfully the spatially varying Galactic components. Furthermore, we
somewhat expect the worst estimate to be the one related to $T_{d}$, since
both LiteBIRD and SO-SAT high frequency channels probe frequency regimes far
from the modified black-body peak ($\nu\sim 600$GHz), making $T_{d}$
essentially degenerate with $\beta_{d}$ at lower frequencies.

Figure 11:  $B-$mode angular power spectra for the (left) LiteBIRD and (right)
SO-SAT cases, estimated with NaMaster onto the recovered CMB maps from
FGBuster performed onto cluster regions. (solid blue thick) Average of
foreground statistical residuals obtained from 20 MC signal+noise simulations
(shaded blue) $1\sigma$ standard deviation of MC simulations, (filled circles)
foreground residual bias estimated with noise-less component separation for
different choices of foreground mis-modelling given by the colorbar. The
shaded grey area in the left panel indicates two extreme cases: the upper
limit is obtained by propagating through FGBuster the $f=1$ mis-modelling to
the spectral parameters, the lower limit spectra (not shown here) is nearly
equal to zero and it is estimated by propagating in the component separation
the $f=0$ case. As a reference we show the primordial B-modes for two
different chosen values of tensor-to-scalar ratio, $r=0.01,0.001$ respectively
in (dashed black) and (solid black). Lensing B-modes are also shown in (solid
thick gray). The spectra in the left and right panels are binned respectively
with $\Delta\ell=5$ and $\Delta\ell=15$.

The case shown in Fig. 8 is an ideal-limit case as the patches have been
derived by clustring onto the same parameter maps that are used to simulate
the frequency channels. In order to introduce a _mis-modelling_ between the
two set of maps, we inject a noise proportional to the per-pixel uncertainties
on $\beta_{s},\beta_{d},T_{d}$, as :

$\tilde{X}=X+f\,\,N_{w}[\sigma(X)],$ (8)

where $X=\beta_{s},\beta_{d},T_{d}$; $N_{w}[\sigma(X)]$ is a random Gaussian
noise map with zero mean and width given by the uncertainty map, and $f$ is a
weighting constant factor, ranging from 0 to 1.

The $f=1$ case represents a pessimistic case where the foregrounds are fully
mis-modelled given the uncertainty budget in the spectral parameters and $f=0$
represents the ideal case without any mis-modelling.

For each mis-modelling case, we estimate the statistical residuals by
performing component separation on maps encoding 20 MC independent
realizations of instrumental noise and astrophysical signal. The systematic
bias is instead estimated by running component separation on noiseless maps.

We therefore run the spectral clustering algorithm on the mis-modelled
parameter maps $\tilde{\beta}_{s},\tilde{\beta}_{d},\tilde{T}_{d}$ for several
values of $f$. E.g. in Fig. 7 (bottom), we show the cluster obtained from
parameter maps with maximum mis-modelling ($f=1$). We note a trend to have
smaller number of clusters with similar sizes and more uniformly distributed
with respect to the $f=0$ case for similar values of $\alpha,\delta$. This is
particularly noticeable for the dust parameters since they are provided with a
higher resolution than the synchrotron one. We conclude that injecting the
mis-modelling noise into the spectral parameter maps tends to homogenize the
typical size of the patches.

We then perform the component separation separately for LiteBIRD and SO-SAT
channels (simulated without mis-modelling) onto the cluster patches derived
with the mis-modelling. This procedure allows us to further set requirements
on the mis-modelling given the residual bias we observe in the recovered CMB
map for different values of $f=[0,0.1,0.5,1]$.

As the output maps of FGBuster are estimated on a partial sky (see Fig. 9), we
estimate the power spectra in the observed regions with NaMaster (Alonso et
al., 2019) to correct the power spectra for the $E-B$ leakage introduced by
masking. We show in Fig.11 the angular power spectra of CMB $B-$mode
polarization anisotropies as estimated from the recovered map outputs of
FGBuster.

We notice that the statistical residuals are not affected by the mis-modelling
since the instrumental noise is the dominant contribution to the post-
component separation residuals. On the other hand, the mis-modelling is
clearly visible when we assess the systematics bias in the residual maps. In
particular, we observe an increase of the systematic bias proportional to $f$,
indicating how the partition is less and less representative with increasing
values of $f$ of the underlying Galactic foreground emission. However, even
with the largest mis-modelling scenario $f=1$, the residuals increase by only
an order of magnitude. As expected, larger $f_{sky}$ (including lower Galactic
latitudes) yields larger residuals in the power spectra, this is the reason
why we observe higher residuals for LiteBIRD with respect to the SO-SAT case).

Furthermore, to show the effective gain of performing parametric component
separation on patches defined with clustering, we consider two limiting cases
where we estimate the systematic residuals by evaluating eq. (7) for LiteBIRD
with the maps of $\beta_{d},\beta_{s},T_{d}$ for $f=0,1$ . Obviously, the case
$f=0$ gives numerically zero residuals as we provide the exact solution to the
parametric fit (i.e. the same spectral parameter maps with which we have
constructed the frequency channels). On the other hand, the case $f=1$
resembles the extreme case where the mis-modelling noise is fully propagated
through the component separation pipeline. In Fig. 11 (left), we show as a
shaded gray area the two extreme cases, with $f=1$ being the largest error
that can be introduced with the mis-modelling defined in eq. (8) (the $f=0$ is
not shown for graphical purposes). The fact that the power spectra obtained
with clustering are well within this range, quantifies the overall benefit of
performing the parametric fitting with the clusters even if clusters are
obtained with the largest mis-modelling case ($f=1$).

Finally, we notice in Fig. 11 that all the $B-$mode systematic residual
spectra tend to converge to a certain residual level at around $\ell\gtrsim
60$. Hereafter, we refer to it as a _partition noise_ being essentially
related to the combination of the typical angular scales of the cluster
patches. This results as a “noise” residual term in the power spectrum, as
that prevents to estimate the variability of the spectral parameters at
multipoles larger than this scale.

###  4.2 Clustering applied on HI clouds maps.

In the previous section, the component separation is performed on partitions
derived from the same maps used in the parametric modeling, i.e.
$\beta_{s},\beta_{d},T_{d}$. However, this relies on the assumption that
spectral parameter maps, from which we derive the cluster patches, are not
contaminated by other Galactic (or extra-galactic) foregrounds. Indeed, this
is commonly the case for both the thermal dust emission, mainly affected by
the Cosmic Infrared Background (CIB) residuals and for synchrotron
contaminated by free-free and AME.

Moreover, given the fact that the use of ancillary datasets for foreground
modeling has been proposed by a number of works (e.g. HI4PI Collaboration et
al. (2016); Clark (2018)), we wish to explore whether such an ancillary
dataset can in principle yield accurate enough results for CMB parameter
estimation purposes. To this end, we use a complementary dataset that traces
foreground dust: HI emission. In particular, we make use of the
parameterization of line-of-sight complexity in the dust distribution
presented by Panopoulou & Lenz (2020). We briefly explain the dataset and our
post-processing of it in the following.

Panopoulou & Lenz (2020) measured the number of clouds along the line of sight
that appear as distinct peaks in HI spectra. They present a peak-finding
algorithm that makes use of the HI4PI survey (HI4PI Collaboration et al.,
2016) over the high Galactic latitude sky and outputs a measure of the number
of clouds in each pixel. They introduced a robust measure of the number of
clouds expected to contribute to the dust emission, which takes into account
the relative contribution of the dust from each cloud along the line of sight,
defined as $\mathcal{N}_{c}$:

$\mathcal{N}_{c}=\sum_{i=1}^{N_{\mathrm{clouds}}}\frac{N_{\mathrm{HI}}^{i}}{N_{\mathrm{HI}}^{\mathrm{max}}}$

where $N_{\mathrm{HI}}^{i}$ refers to the column density of the $i$-th cloud
along the sightline and $N_{\mathrm{HI}}^{\mathrm{max}}$ is the highest value
of column density for a cloud in the same direction. We therefore expect non-
integer values for $\mathcal{N}_{c}$ especially when there is an imbalance in
the column density for clouds along a sightline. The $\mathcal{N}_{c}$ map was
estimated by segmenting the sky into large _superpixels_ at nside=128
(corresponding to $\sim 30$ arcmin), with each superpixel encoding 64 pixels
of the HI4PI map. Recently, Pelgrims et al. (2021) adopted the maps publicly
released by Panopoulou & Lenz (2020) to estimate the line-of-sight frequency
decorrelation of dust polarization in the _Planck_ maps. For the purposes of
our clustering analysis, we require maps at nside=32 and nside=64. We thus
repeat the analysis in Panopoulou & Lenz (2020) to produce maps at these lower
nsides, instead of downgrading the publicly available maps, as recommended by
the authors. The publicly available $\mathcal{N}_{c}$ maps did not include
uncertainties. As our clustering relies on the use of uncertainties, we
calculate the per-pixel uncertainty in the $\mathcal{N}_{c}$ maps that we use
as described in Appendix D.

As discussed in Panopoulou & Lenz (2020, Sect. 5.3.1), the derived maps do not
present any imposed spatial coherence at scales above the superpixel size,
leading to the eventuality of discontinuities between neighbouring pixels.
However, clustering methodologies presented in this work, can implement the
spatial coherence above a given pixel scale.

Figure 12: (top) Maps of $\mathcal{N}_{c}$ (left) and its uncertainties
$\sigma(\mathcal{N}_{c})$ (right) evaluated as in eq. (31). (bottom left)
Partition of $\mathcal{N}_{c}$; each color identifies a different clustering
region. (bottom right) Relative errors $\delta\mathcal{N}_{c}$ estimated with
(eq. (9)).

We employ the $\mathcal{N}_{c}$ and the associated uncertainties as in (31)
with maps at nside=32 (shown in top panel of Fig. 12) and run the spectral
clustering methodology outlined in Sect. 3. We then evaluate the total
variance of the partition in a similar grid to the one shown in Fig.6 and find
the minimum variance to be at $\alpha=0.10$ and $\delta=0.20$, corresponding
to $K=1481$ clusters.

Figure 13: $B-$mode residuals from the recovered CMB maps output of FGBuster
performed within multiple regions defined by clustering of $\mathcal{N}_{c}$
map for (blue) LiteBIRD and (orange) SO-SAT frequency channels, respectively
on $f_{sky}=0.35$ and $0.22$. (filled circles) Foreground residual bias
performed on noiseless simulations, (shaded area) $1\sigma$ standard deviation
of statistical residuals evaluated from 20 MC signal and noise simulations
(the average of the 20 MC realization is shown as solid line). (dashed)
Foreground residual bias estimated from the component separation run performed
by flipping the map of $\mathcal{N}_{c}$ clusters North-to-South. These
spectra show an excess at lower multipoles, demonstrating that the
$\mathcal{N}_{c}$ map provides useful constraints for parameter estimation.

To evaluate how well the clustering morphologies track the input map, we
estimate the median value of pixels of $\mathcal{N}_{c}$ belonging to the same
clusters and produce a _binned_ map, $\mathcal{N}^{bin}_{c}$. We can then
estimate the relative error defined as

$\delta\mathcal{N}_{c}=|\mathcal{N}_{c}-\mathcal{N}^{bin}_{c}|/\mathcal{N}_{c}.$
(9)

The bottom row of Fig.12, shows in the left panel the partition estimated with
spectral clustering and in the right one the relative errors. We observe that
especially around the North and South Galactic Poles
($90^{\circ}<|b|<60^{\circ}$) the clustering optimization tends to produce
smaller clusters, essentially encoding few pixels per cluster. On the
contrary, at intermediate Galactic latitudes, clusters contain a larger number
of pixels. This might be related to the fact that the $\mathcal{N}_{c}$ maps
are produced on too large pixelization (nside=32) resulting in an overall
reduction of the uncertainty budget and a consecutive increase in SNR of the
estimates. As already stated in Subsection 4.1, in the presence of very high
SNR features, the optimal partition leads to _pixel-sized_ clusters. We will
devote a further investigation on this respect in a future work.

Notice that in the regions where the uncertainties are small, the error is
below $\sim 10\%$. Vice versa, in large uncertainties regions the relative
error can be as much as $50\%$ although those regions are localized mainly in
the North Galactic Hemisphere. However, we remark here that these larger
errors should not be compared with the ones shown in Fig. 8 as larger errors
on $\mathcal{N}_{c}$ do not necessarily translate on larger residuals on the
recovered dust parameters.

As already mentioned in Panopoulou & Lenz (2020), one of the interesting
applications of the map of number of clouds per line of sight is to inform
realistic models of polarized dust emission for parametric component
separation methods as the one we adopt in Sect. 4.1. In this specific case, we
utilize the patches derived from the $\mathcal{N}_{c}$ map (bottom left panel
in Fig. 12) to perform the component separation estimation of dust and CMB
(simulated at nside=64). In fact, since the $\mathcal{N}_{c}$ map traces the
dust emission, we do not include the synchrotron as well as all the $\nu<90$
GHz frequency channels in this analysis. This helps in better singling out the
effects of residuals due to the dust component only.

Furthermore, to run NaMaster we regularize the HI4PI observation area
($f_{sky}=0.55\%$) by apodizing it, yielding to $\sim 0.35\%$ of the sky . For
the run with SO-SAT we combine this mask with the SO-SAT one shown in Fig. 9
(bottom). The combined observation patches results in a further reduction,
$f_{sky}\sim 22\%$.

In Fig. 13, we show the $B-$mode power spectra of the residuals in the
recovered CMB map estimated with NaMaster. We observe that the statistical
residuals are reduced by about one order of magnitude compared to the ones
shown in Fig. 11. This is mostly due to the fact that we have excluded the low
frequency channels that have lower sensitivity (see Table 1) yielding larger
post-component separation noise when taken into account. Interestingly, the
residual bias achieves levels comparable to $r\lesssim 10^{-3}$. It therefore
seems that the parameters from the component separation performed on the dust-
emission-agnostic $\mathcal{N}_{c}$-based sky partitioning are accurately
representing the underlying signal.

One might wonder whether this dust-emission-agnostic clustering constitutes a
real benefit in parameter estimation or it is only related to the fact that we
are simply partitioning the northern and southern polar caps with small
patches. In other words: _is there intrinsic value in the information
contained in the number-of-clouds map?_

To answer this question, we disentangle the sky partition from the underlying
physical information by flipping the cluster coordinates North to South.
However, since the North and South Polar caps have different footprints, we
firstly flip the map and then apply the regularized mask to the
$\mathcal{N}_{c}$ map. Encoding the mask on a smaller fraction of sky ensures
that each pixel in the flipped map falls within the $\mathcal{N}_{c}$ map
footprint.

We then perform the parametric fit with the same combination of frequency
channels as before but with the flipped cluster regions. The thick dashed
lines shown in Fig. 13 indicate the power spectra of systematic residuals from
the output of this component separation run and we notice a remarkable
increase of the residuals especially at lower multipoles both for SO-SAT and
LiteBIRD .

A visual inspection of the relative error and standard deviation maps
estimated in both cases for $T_{d}$ and $\beta_{d}$, reveals that both the
error and the standard deviation (particularly in $T_{d}$) increase when we
consider the case with flipped regions. This also indicates that the bias we
observe in Fig. 13 is mostly due to wrong estimation of the $T_{d}$ parameter.

This is not unexpected since $\mathcal{N}_{c}$ has been shown to be correlated
with $T_{d}$ (Panopoulou & Lenz, 2020). The higher values of $\mathcal{N}_{c}$
coincide with regions where the emission from two physically distinct families
of clouds overlap (Low Velocity and Intermediate Velocity clouds). One family
of clouds resides at larger distances from the Galactic plane and is found to
have higher dust temperatures, presumably due to dust shattering as it falls
onto the mid-plane (Planck Collaboration et al., 2011). Therefore, a mis-
modelling of the cloud population (e.g. by flipping the north-south map)
naturally leads to a wrong estimate of $T_{d}$. We thus can conclude that
partitioning the sky observations with patches inferred from the
$\mathcal{N}_{c}$ map can improve the performances of parametric fitting
methodologies.

Finally, we observe also for this case the effects of a partition noise in the
power spectra: with $B-$mode residuals in Fig.13 (dashed lines and filled
circles) converging at around $\ell\sim 60$. A partition with smaller patches,
as the one shown in Fig. 7, results into shifting to higher multipoles the
contribution due to this partition noise. Grumitt et al. (2020) proposed a
hierarchical approach to overcome exactly this limitation and we plan to
integrate it in a future work.

###  4.3 Estimates on $r$

We recall the reader that to specifically assess the bias introduced by
foregrounds, we do not include the CMB emission in the simulated frequency
maps. However, we instead fit for CMB in the parametric component separation.
Given the $B-$mode power spectra shown in Fig. 11 and 13, we estimate the
value of residual bias in terms of $r$ by evaluating the likelihood
$\mathcal{L}(r)$ in the binned multipole domain, $\ell_{b}$, given the fact
that the power spectra have been corrected by the mode coupling with NaMaster.

The likelihood function reads as:

$\ln\mathcal{L}(r)=\sum_{\ell_{b}=\ell_{\rm min}}^{\ell_{\rm
max}}\ln\left(-f_{sky}\Delta\ell\frac{2\ell_{b}+1}{2}\left[\frac{\hat{C}_{\ell_{b}}}{C_{\ell_{b}}}+\ln
C_{\ell_{b}}\right]\right),$ (10)

where $\ell_{\rm max}=200$ , $\ell_{\rm min}=2\,(30)$ and $\Delta\ell=5(15)$
for LiteBIRD (SO-SAT ) and $\hat{C}_{\ell_{b}}$ ($C_{\ell_{b}}$) is the
measured (modeled) B-mode power spectrum. As we are interested in estimating
the bias of systematic residuals we assume the measured $B-$mode spectrum to
encode all the contributions but the primordial one ($r=0$), i.e. :

$\hat{C}_{\ell_{b}}=C_{\ell_{b}}^{\rm lens}+\langle C_{\ell_{b}}^{\rm
fg,tot}\rangle,$ (11)

where $\langle C_{\ell_{b}}^{\rm fg,tot}\rangle$ is the expected statistical
noise power spectrum residuals after component separation obtained by
averaging the spectra of 20 MC signal+noise simulations (Errard & Stompor,
2019) and $C_{\ell_{b}}^{\rm lens}$ is the lensing B-mode power spectrum . The
modeled power spectrum is given as:

$C_{\ell_{b}}=rC_{\ell_{b}}^{\rm tens}+C_{\ell_{b}}^{\rm lens}+\langle
C_{\ell_{b}}^{\rm fg,tot}\rangle-C_{\ell_{b}}^{\rm fg,sys},$ (12)

where $C_{\ell_{b}}^{\rm tens}$ is the tensor mode with $r=1$ and
$C_{\ell_{b}}^{\rm fg,sys}$ is the spectrum encoding the systematics bias due
to mismodelling the foreground emission.

Thus, the bias on $r$, $\Delta r$ is defined as the value that maximizes the
likelihood function:

$\left.\frac{d\mathcal{L}(r)}{dr}\right|_{r=\Delta r}=0,$ (13)

and the error on $r$, $\delta r$ is defined as the value covering the 68% area
of the total likelihood function, i.e.:

$\frac{\int_{0}^{\delta
r}\mathcal{L}(r)dr}{\int_{0}^{\infty}\mathcal{L}(r)dr}=0.68.$ (14)

The values reported in Table 2 are extensively discussed in the following
section.

##  5 Discussion

We devote this section to discuss in more detail the results presented in
Section 4. We estimate how the bias from systematic residuals for different
choices of the mis-modelling parameter $f$ affects the detection of primordial
CMB $B-$modes and compare with the requirements for SO-SAT and LiteBIRD .

As reported in Table 2 and shown in Fig. 14, the effect of foreground mis-
modelling is visible in terms of systematic bias on $r$ estimate but slightly
affects the uncertainties (which are instead driven by the sensitivity of the
experiment, $f_{sky}$, frequency coverage, etc…). Although the bias for
LiteBIRD is always $<10^{-3}$ for most cases, we observe an increase
proportional to the mis-modelling. Without further optimizing the Galactic
mask for polarization the values in the top row of Table 2 show that we can
meet the latest LiteBIRD requirements ($\delta r<10^{-3}$, Sugai et al., 2020)
with a mis-modelling $f\leq 0.5$.

We notice that the bias on $r$ obtained in this work for SO-SAT is slightly
lower than the reported values in Simons Observatory Collaboration et al.
(2019, Table 4), indicating a net improvement thanks to this technique. The
larger uncertainties and the lower systematic bias observed in SO-SAT , with
respect to LiteBIRD ones, are mostly due to the smaller coverage both in
frequency and in the angular scales observed. Moreover, the $r$ estimates for
SO-SAT are assessed by evaluating the likelihood in a limited range of $\ell$
($30<\ell<200$). Vice versa, LiteBIRD benefits of the full-sky angular range
allowing also the large angular scales to be included in the likelihood
($2<\ell<200$).

Clustering on $\beta_{s},\beta_{d},T_{d}$ |  |
---|---|---
FG Mismodelling $f$ | LiteBIRD | SO-SAT
$0$ | $(0.22\pm 0.66)\times 10^{-3}$ | $(0.08\pm 1.21)\times 10^{-3}$
$0.1$ | $(0.28\pm 0.66)\times 10^{-3}$ | $(0.10\pm 1.21)\times 10^{-3}$
$0.5$ | $(0.59\pm 0.67)\times 10^{-3}$ | $(0.14\pm 1.21)\times 10^{-3}$
$1.0$ | $(2.18\pm 0.68)\times 10^{-3}$ | $(0.20\pm 1.19)\times 10^{-3}$
Clustering on $\mathcal{N}_{c}$ | $(0.09\pm 0.26)\times 10^{-3}$ | $(0.04\pm 0.16)\times 10^{-3}$
Table 2: Estimated values for the bias and the uncertainty on $r$ obtained by
evaluating the likelihood in eqs. (10), (13) and (14).  Figure 14: $r$
estimates and the accompanying uncertainties evaluating the likelihood in eqs.
(10), (13) and (14). We adopt the same color scheme as in Fig. 11 to indicate
different amounts of mis-modelling.

Alternatively to foreground parameter maps, we also apply clustering to
partition the map of number of dust clouds along the line of sight,
$\mathcal{N}_{c}$, presented in Panopoulou & Lenz (2020). We perform the same
parametric fitting as before on the patch derived from $\mathcal{N}_{c}$, and
find similar level of residuals post-component separation for both SO-SAT and
LiteBIRD cases when compared with ones obtained from clusters derived with
foreground parameter maps. Not only do we find that the partition obtained is
representative of the dust emission, but also that the patches inferred from
$\mathcal{N}_{c}$ yield lower residuals in the estimates of the thermal dust
temperature parameter, $T_{d}$.

We further estimate the leakage in terms of $r$ also for this case of
application by means of the likelihood defined in (10) and we found
$\hat{r}=(0.09\pm 0.26)\times 10^{-3}$ and $(0.04\pm 0.16)\times 10^{-3}$
respectively for LiteBIRD and SO-SAT. The reason we obtain lower biases with
respect to the one reported in Table 2 for the $f=0$ mismodelling case is
mainly due to the fact that low frequency channels are not accounted for in
the component separation and we assume the synchrotron to be neglible at
$\nu>90$ GHz. Moreover, we get a lower noise bias for SO-SAT with respect to
LiteBIRD mainly because to the fact that we do not account for the $1/f$ noise
in our simulations for both instruments. However, the low bias on $r$ is a
stricking indication that the $\mathcal{N}_{c}$ map can be a very promising
tracer to probe the thermal dust emission in the range of $100-300$ GHz.
Although Pelgrims et al. (2021) showed that the $\mathcal{N}_{c}$ maps can be
used as a tracer for frequency decorrelation in the dust polarization map, the
results shown in this paper offer a further indication of their potential.

As already mentioned, this is not the first time that clustering methodologies
are employed to segment the sky into patches where to perform foreground
cleaning (see Grumitt et al. (2020) where they used the mean-shift algorithm).
However, by comparing the maps of clusters in Fig. 7 with those in Fig. 2 of
(Grumitt et al., 2020), we would like to address several reasons why the two
methods yield very different partitions.

One immediate consideration is the fact that mean-shift and spectral
clustering are based on different algorithms that naturally might yield
different results even when applied on samples with fewer features than the
ones in this context.

The clusters obtained in Grumitt et al. (2020) are larger than the ones shown
in Fig. 7, although we also notice a trend in estimating larger clusters
particularly for the $\beta_{s}$ clusters. This makes the method better suited
for low-resolution and low-frequency channels where synchrotron emission
dominates121212An interesting application in Grumitt et al. (2020) was indeed
the combination of the C-BASS data at 5 GHz with the LiteBIRD frequency
channels., but it was essentially unable to track the finer structure of the
dust parameters.

One of the major differences is that we derive the clusters _separately_ for
each spectral parameter, whereas in Grumitt et al. (2020) the dust and
synchrotron spectral indices $\beta_{s}$ and $\beta_{d}$ are combined together
to derive a single cluster map used to perform the parametric fit (see Fig. 2
of Grumitt et al., 2020).

Furthermore, the features in Grumitt et al. (2020) included the two spectral
indices together with the cartesian coordinates on the sphere $x,y,z$. On the
contrary, our analysis does not rely on the coordinates as features since the
adjacency is defined from the distance between pixels (Section 3.1). We
further account also for $T_{d}$ as an extra-parameter for the thermal dust
emission, leading to a reduction in the foreground bias especially at high
Galactic latitudes.

A novel feature of the technique presented in this work is to perform the
clustering by accounting for the statistical significance of the measurements
reported in the spectral parameter maps and their uncertainties. We weigh the
pixel adjacency in such a way that the significance plays a non-negligible
role in defining the pixel similarities especially in regions of high SNR.
Moreover, we drastically reduce the dimensionality of our problem by
performing the clustering on the eigenspace spanned by the eigenvectors
related to the smallest 256 eigenvalues. This avoids high-dimensionality
problems occurring frequently in clustering algorithms with large number of
features (e.g. the number of pixels in a nside=64 HEALPix map see Section 4 of
Grumitt et al. (2020)).

Finally, in Appendix B, we report also the clustering analysis performed on
the GNILC $\beta_{d}$ and $T_{d}$ maps. The morphologies are slightly
different with respect to the ones observed in the clusters in Fig. 7 as the
spectral parameter features and the uncertainties in the two datasets are
different.

##  6 Conclusions

In this work, we aim at identifying regions in the sky within which the
Galactic foreground emission can be assumed homogeneous. We define a
partitioning of the full celestial sphere by means of a novel technique based
on a spectral clustering algorithm embedded in the $S^{2}$ manifold. In
contrast to previous applications of clustering for CMB studies, this method
works directly on the celestial sphere instead of mapping points onto a
Cartesian grid. It also takes into account the uncertainties in the clustering
procedure. To our knowledge, this is the first time where image segmentation
has been performed on features defined on the $S^{2}$ manifold.

Another key advancement of our algorithm is the use of an objective criterion
for selecting the optimal partition parameters ($\alpha$ and $\delta$). Many
_off-the-shelf_ algorithms rely on subjective measures, such as visual
inspection, to select a ‘reasonable’ segmentation of the domain (e.g. the
maximum pixel size parameter in mean-shift). The algorithm presented here
borrows established metrics for partition measure evaluation in clustering
problems to produce an optimal partition. We rely on the comparison of within-
and between-cluster variances for this, while adding new modifications to
incorporate the presence of uncertainties in the data.

We apply this clustering approach to evaluate Galactic contamination in CMB
polarization studies. We firstly characterize the emission of thermal dust and
synchrotron available from the latest observations of the spectral parameter
maps (together with their uncertainties) to divide the celestial sphere into
regions with spectral clustering. We then perform a parametric component
separation for each region simulating the nominal frequency channels from two
forthcoming CMB experiments observing from the ground and from space. The
optimal clustering partition yields low residuals in the CMB polarization maps
after component separation, even in the case where the Galactic emission is
mis-modelled up to some extent.

We further apply clustering to partition the map of number of Hi clouds along
the line of sight. Although this is one of the first applications using this
kind of tracer, we find that clusters derived from the $\mathcal{N}_{c}$ map
yield a residual bias as low as $r\sim 10^{-3}$ and that the dust spectral
parameters are estimated on these regions up to $<10\%$ relative errors. This
is a remarkable result since the constraints on $T_{d}$ from far infrared data
have been so far complicated by the degeneracy with the spectral index
$\beta_{d}$. On the other hand, the use of ancillary data combined with the
latest data-science techniques shows promising results in order to improve the
modelling of Galactic foregrounds.

In conclusion, this technique could be applied on a wide range of contexts
involving spherical images, e.g. Earth images, identification of Solar
features, wide astronomical surveys and it can be easily extended to higher
dimensional datasets, like 3D surveys or cosmological simulations.

## Acknowledgements

This research used resources of the National Energy Research Scientific
Computing Center (NERSC), a U.S. Department of Energy Office of Science User
Facility located at Lawrence Berkeley National Laboratory, operated under
Contract No. DE-AC02-05CH11231. Giuseppe Puglisi acknowledges support by the
COSMOS & LiteBIRD networks of the Italian Space Agency. The authors thank
Clement Leloup, Carlo Baccigalupi for having read the paper thoroughly.

Giuseppe Puglisi would like to thank: Jonathan Aumont, Mathieu Remazailles,
Jens Chluba, Aditya Rotti, Susanna Azzoni, Leo Vacher, Hans Christian Eriksen,
Nicoletta Krachmalnicoff for useful comments and discussions.

Georgia Panopoulou acknowledges support for this work by NASA through the NASA
Hubble Fellowship grant #HST-HF2-51444.001-A awarded by the Space Telescope
Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

## Data Availability

The whole clustering analysis presented here has been collected into a python
package, ForeGround Clusters (FGClusters
131313https://github.com/giuspugl/fgcluster).

The inputs used throughout this paper together with the outputs obtained with
the Clustering technique have been made publicly available online141414
https://portal.nersc.gov/project/sobs/users/giuspugl/Clustering and in the
Harvard Dataverse :

  * •

$\mathcal{N}_{c}$ map : https://doi.org/10.7910/DVN/XAMJ4X;

  * •

PySM spectral parameters: https://doi.org/10.7910/DVN/WJUEFA;

  * •

GNILC dust parameters: https://doi.org/10.7910/DVN/02SCHB.

## References

  * Alonso et al. (2017) Alonso D., Dunkley J., Thorne B., Næss S., 2017, Phys. Rev. D, 95, 043504
  * Alonso et al. (2019) Alonso D., Sanchez J., Slosar A., 2019, Monthly Notices of the Royal Astronomical Society, 484, 4127–4151
  * Basak & Delabrouille (2011) Basak S., Delabrouille J., 2011, Monthly Notices of the Royal Astronomical Society, 419, 1163
  * Belkin & Niyogi (2003) Belkin M., Niyogi P., 2003, Neural Computation, 15, 1373
  * Bennett et al. (1992) Bennett C. L., et al., 1992, ApJ, 396, L7
  * Berry & Giannakis (2020) Berry T., Giannakis D., 2020, Communications on Pure and Pplied Mathematics, 73, 689
  * Berry & Sauer (2016) Berry T., Sauer T., 2016, Appl.Comp. Harmonic Analysis, 40, 439
  * Bianchini et al. (2020) Bianchini F., et al., 2020, The Astrophysical Journal, 888, 119
  * Bobin et al. (2008) Bobin J., Moudden Y., Starck J. L., Fadili J., Aghanim N., 2008, Statistical Methodology, 5, 307
  * Boyd et al. (2018) Boyd Z., Bae E., Tai X., Bertozzi A. L., 2018, SIAM J. Appl. Math., 78, 2439
  * Chluba et al. (2017) Chluba J., Hill J. C., Abitbol M. H., 2017, Monthly Notices of the Royal Astronomical Society, 472, 1195–1213
  * Choi et al. (2020) Choi S. K., et al., 2020, The Atacama Cosmology Telescope: A Measurement of the Cosmic Microwave Background Power Spectra at 98 and 150 GHz (arXiv:2007.07289)
  * Clark (2018) Clark S. E., 2018, ApJ, 857, L10
  * Coifman & Lafon (2006) Coifman R., Lafon S., 2006, Appl.Comp. Harmonic Analysis, 21, 5
  * Delabrouille et al. (2003) Delabrouille J., Cardoso J. F., Patanchon G., 2003, MNRAS, 346, 1089
  * Dunkley et al. (2009) Dunkley J., et al., 2009, AIP Conference Proceedings, 1141, 222
  * Eriksen et al. (2008) Eriksen H. K., Jewell J. B., Dickinson C., Banday A. J., Górski K. M., Lawrence C. R., 2008, ApJ, 676, 10
  * Errard & Stompor (2019) Errard J., Stompor R., 2019, Physical Review D, 99
  * Errard et al. (2011) Errard J., Stivoli F., Stompor R., 2011, Phys. Rev. D, 84, 063005
  * Fredholm (1903) Fredholm I., 1903, Acta Mathematica, 27, 365
  * Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759
  * Grumitt et al. (2020) Grumitt R. D. P., Jew L. R. P., Dickinson C., 2020, Monthly Notices of the Royal Astronomical Society, 496, 4383–4401
  * Guth (1981) Guth A. H., 1981, Phys. Rev. D, 23, 347
  * H. Hu & Bertozzi (2015) H. Hu J. S., Bertozzi A. L., 2015, Proceedings of 10th International Conference, EMMCVPR, Hong Kong, Energy Minimization Methods in Computer Vision and Pattern Recognition
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration B., et al., 2016, Astronomy & Astrophysics, 594, A116
  * Hansen et al. (2006) Hansen F. K., Banday A. J., Eriksen H. K., Górski K. M., Lilje P. B., 2006, ApJ, 648, 784
  * Haslam et al. (1982) Haslam C. G. T., Salter C. J., Stoffel H., Wilson W. E., 1982, A&AS, 47, 1
  * Hinshaw et al. (2009) Hinshaw G., et al., 2009, The Astrophysical Journal Supplement Series, 180, 225–245
  * Hu & White (1997) Hu W., White M. J., 1997, New Astron., 2, 323
  * Irfan et al. (2019) Irfan M. O., Bobin J., Miville-Deschênes M.-A., Grenier I., 2019, A&A, 623, A21
  * Khatri (2019) Khatri R., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 039
  * Krzanowski & Lai (1988) Krzanowski W. J., Lai Y. T., 1988, Technical Report 1, A Criterion for Determining the Number of Groups in a Data Set Using Sum-of-Squares, https://www.jstor.org/stable/2531893?seq=1&cid=pdf-reference#references_tab_contents. https://www.jstor.org/stable/2531893?seq=1&cid=pdf-reference#references_tab_contents
  * Leach et al. (2008) Leach S. M., et al., 2008, A&A, 491, 597
  * M. Reuter & Peinecke (2005) M. Reuter F.-E. W., Peinecke N., 2005, SPM ’05: Proceedings of the 2005 ACM symposium on Solid and physical modeling, pp 101–106
  * Maino et al. (2007) Maino D., Donzelli S., Banday A. J., Stivoli F., Baccigalupi C., 2007, MNRAS, 374, 1207
  * Mangilli et al. (2021) Mangilli A., Aumont J., Rotti A., Boulanger F., Chluba J., Ghosh T., Montier L., 2021, Astronomy & Astrophysics, 647, A52
  * Meng et al. (2017) Meng Z., Merkurjev E., Koniges A., Bertozzi A. L., 2017, IPOL Journal of Image Processing On Line
  * Mihaylov & Spallanzani (2016) Mihaylov G., Spallanzani G., 2016, International Journal of Prognostics and Health Management, Special Issue Big Data and Analytics
  * Miville-Deschênes et al. (2008) Miville-Deschênes M.-A., Ysard N., Lavabre A., Ponthieu N., Macías-Pérez J. F., Aumont J., Bernard J. P., 2008, Astronomy & Astrophysics, 490, 1093–1102
  * Panopoulou & Lenz (2020) Panopoulou G. V., Lenz D., 2020, Maps of the number of HI clouds along the line of sight at high galactic latitude (arXiv:2004.00647), doi:10.3847/1538-4357/abb6f5
  * Pelgrims et al. (2021) Pelgrims V., Clark S. E., Hensley B. S., Panopoulou G. V., Pavlidou V., Tassis K., Eriksen H. K., Wehus I. K., 2021, Astronomy & Astrophysics, 647, A16
  * Planck Collaboration et al. (2011) Planck Collaboration et al., 2011, A&A, 536, A24
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A12
  * Planck Collaboration et al. (2016a) Planck Collaboration et al., 2016a, A&A, 594, A9
  * Planck Collaboration et al. (2016b) Planck Collaboration et al., 2016b, A&A, 594, A10
  * Planck Collaboration et al. (2016c) Planck Collaboration et al., 2016c, A&A, 594, A10
  * Planck Collaboration et al. (2016d) Planck Collaboration et al., 2016d, A&A, 594, A13
  * Planck Collaboration et al. (2016e) Planck Collaboration et al., 2016e, Astronomy & Astrophysics, 596, A109
  * Planck Collaboration et al. (2020a) Planck Collaboration et al., 2020a, A&A, 641, A4
  * Planck Collaboration et al. (2020b) Planck Collaboration et al., 2020b, Astronomy & Astrophysics, 641, A6
  * Puglisi, G. et al. (2018) Puglisi, G. Poletti, Davide Fabbian, Giulio Baccigalupi, Carlo Heltai, Luca Stompor, Radek 2018, A&A, 618, A62
  * Puglisi et al. (2017) Puglisi G., Fabbian G., Baccigalupi C., 2017, MNRAS, 469, 2982
  * Remazeilles et al. (2018) Remazeilles M., et al., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 023
  * Rustamov (2007) Rustamov R. M., 2007, Proceedings of Eurographics Symposium on Geometry, 15
  * Seljak & Zaldarriaga (1997) Seljak U., Zaldarriaga M., 1997, Physical Review Letters, 78, 2054
  * Simons Observatory Collaboration et al. (2019) Simons Observatory Collaboration et al., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 056–056
  * Singer & Wu (1997) Singer A., Wu H.-T., 1997, Information and Inference: A Journal of the IMA, 6, 58
  * Starobinsky (1982) Starobinsky A. A., 1982, Physics Letters B, 117, 175
  * Stolyarov et al. (2005) Stolyarov V., Hobson M. P., Lasenby A. N., Barreiro R. B., 2005, MNRAS, 357, 145
  * Stompor et al. (2008) Stompor R., Leach S., Stivoli F., Baccigalupi C., 2008, Monthly Notices of the Royal Astronomical Society, 392, 216
  * Stompor et al. (2016a) Stompor R., Errard J., Poletti D., 2016a, Phys. Rev. D, 94, 083526
  * Stompor et al. (2016b) Stompor R., Errard J., Poletti D., 2016b, Phys. Rev. D, 94, 083526
  * Sugai et al. (2020) Sugai H., et al., 2020, Journal of Low Temperature Physics, 199, 1107–1117
  * Szydlarski, M. et al. (2014) Szydlarski, M. Grigori, L. Stompor, R. 2014, A&A, 572, A39
  * Tassis & Pavlidou (2015) Tassis K., Pavlidou V., 2015, Monthly Notices of the Royal Astronomical Society: Letters, 451, L90
  * The Polarbear Collaboration et al. (2017) The Polarbear Collaboration et al., 2017, ApJS
  * The Polarbear Collaboration et al. (2019) The Polarbear Collaboration et al., 2019, A Measurement of the Degree Scale CMB B-mode Angular Power Spectrum with POLARBEAR (arXiv:1910.02608)
  * Thorne et al. (2017) Thorne B., Dunkley J., Alonso D., Næss S., 2017, Monthly Notices of the Royal Astronomical Society, 469, 2821–2833
  * Thorne et al. (2019) Thorne B., et al., 2019, arXiv e-prints, p. arXiv:1905.08888
  * Tristram et al. (2021) Tristram M., et al., 2021, Astronomy & Astrophysics, 647, A128
  * Von Luxburg (2007) Von Luxburg U., 2007, Technical Report 4, A Tutorial on Spectral Clustering, www.springer.com.. (arXiv:0711.0189v1), www.springer.com.
  * Wagner-Carena et al. (2019) Wagner-Carena S., Hopkins M., Rivero A. D., Dvorkin C., 2019, Monthly Notices of the Royal Astronomical Society, 494, 1507
  * Zelnik-Manor & Perona (2004) Zelnik-Manor L., Perona P., 2004, Technical report, Self-Tuning Spectral Clustering, http://www.vision.caltech.edu/lihi/Demos/SelfTuningClustering.html. http://www.vision.caltech.edu/lihi/Demos/SelfTuningClustering.html
  * Zhang et al. (2018) Zhang C., Zhu G., Chen M., Chen H., Wu C., 2018, Technical report, Image Segmentation Based on Multiscale Fast Spectral Clustering. (arXiv:1812.04816v1)
  * Zhao & Song (2018) Zhao C., Song J. S., 2018, Frontiers in Applied Mathematics and Statistics, 4, 1
  * Zonca et al. (2019) Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, Journal of Open Source Software, 4, 1298

##  Appendix A A geometric intuition on spectral image segmentation

The Laplace-Beltrami operator $\Delta_{g}$ is the generalisation of the
Laplace operator $\Delta$ (divergence of the gradient) on a Riemannian
manifold with metric $g$. In Local coordinates it can be written as:

$\Delta_{g}f=\frac{1}{\sqrt{|g|}}\partial_{i}\left(\sqrt{|g|}g^{ij}\partial_{j}f\right)$

The Laplace-Beltrami operator is a second order linear elliptic differential
operator. It is the Hodge (connection) Laplacian acting on functions, it is by
construction self-adjoint and positive-semidefinite.

The spectral properties of $\Delta_{g}$ on a compact differentiable manifold
$M$ are well-known. The spectrum $\lambda_{i}$ is discrete, $\lambda_{1}=0$ is
always an eigenvalue, all eigenvalues $0\leq\lambda_{1}...\leq\lambda_{j}\leq$
are positive. The eigenfunctions $\varphi_{i}$ of the Laplace-Beltrami
operator on $M$ form a complete orthonormal basis of the space $L^{2}(M)$,
i.e. any function in $L^{2}(M)$ can be written as a convergent series in
$L^{2}(M)$ with real coefficients. This is the celebrated Sturm-Liouville
decomposition of $L^{2}$ functions on smooth compact manifolds. The discussion
of the eigenvalue problem on manifolds with boundary (with Dirichlet or
Newmann boundary conditions) is often irrelevant in the discrete applications.

The Dirichlet energy functional of a function on $M$ is defined as

$E[f]=\frac{1}{2}\int_{M}||\nabla f(x)||^{2}d\omega_{g},$

where $\omega_{g}$ denotes the volume form over a Riemannian manifold. The
Dirichlet energy measures the “variability” of $f$ on $M$.

The Dirichlet functional evaluated on the eigenfunctions $\varphi_{i}$ of
$\Delta_{g}$ is monotonously increasing as the corresponding eigenvalues
increase. The Stokes theorem applied on a function $f$ and a vector field $X$
implies that $-div$ and $\nabla$ are formally adjoint operators (see Belkin &
Niyogi (2003)), i.e.:

$\int_{\mathcal{M}}\langle X,\nabla f\rangle=-\int_{\mathcal{M}}(\nabla\cdot
X)f.$

So that:

$\displaystyle\int_{\mathcal{M}}\parallel\nabla
f(x)\parallel^{2}d\omega_{g}=\int_{\mathcal{M}}$ $\displaystyle\langle\nabla
f,\nabla f\rangle d\omega_{g}=$ $\displaystyle-\int_{\mathcal{M}}$
$\displaystyle(\nabla\cdot\nabla
f)fd\omega_{g}=-\int_{\mathcal{M}}f(\Delta_{g}f)d\omega_{g}.$

The former equation implies:

$\displaystyle E[\varphi_{i}]=\frac{1}{2}\int_{M}$
$\displaystyle||\nabla\varphi_{i}(x)||^{2}d\omega_{g}=\frac{1}{2}\int_{M}\langle\nabla\varphi_{i},\nabla\varphi_{i}\rangle
d\omega_{g}=$ $\displaystyle\frac{1}{2}\int_{M}$
$\displaystyle\varphi_{i}(\Delta_{g}\varphi_{i})d\omega_{g}=\frac{1}{2}\lambda_{i}$

For $f=\sum_{i}\alpha_{i}\varphi_{i}$ this expression linearly extends to

$E[f]=\frac{1}{2}\sum_{i}\alpha_{i}^{2}\lambda_{i}.$

The Laplace-Beltrami operator on $L^{2}(M)$ can be alternatively defined
though the quadratic functional $E[f]$:

$E[f]=\langle\nabla f,\nabla
f\rangle_{L^{2}(M)}=\langle\Delta_{g}f,f\rangle_{L^{2}(M)}.$ (15)

Interestingly, $\nabla E[f]=\Delta_{g}f$ i.e. the Laplace-Beltrami operator is
the variational derivative of Dirichlet energy in $L^{2}(M)$.

###### Proposition 3

The function that minimises the Dirichlet energy and is orthogonal in
$L^{2}(M)$ to the space spanned by $\varphi_{0},...,\varphi_{i}$ is
$\varphi_{i+1}$.

###  A.1 Fredholm theory and discrete differential operators

Digitally produced images resemble smooth functions sampled discretely over
regular domains (a grid of pixels), the above construction needs to be
replicated in a discrete form. For some general ideas about the discretisation
of differential operators acting on sections of tensor bundles, we refer the
reader to Mihaylov & Spallanzani (2016)). Manifold learning techniques that
include consistent discretisations of key elements of Riemannian geometry
(vector fields, connections) have been described in Singer & Wu (1997).
Recently this approach has been further developed in Berry & Giannakis (2020)
including differential forms, Laplace-De Rham operators etc.

Let us consider the following differential equation:

$\mathcal{D}f=g,$ (16)

where $\mathcal{D}:L^{2}(M)\longrightarrow L^{2}(M)$ is a linear differential
operator. The Green’s function of a differential operator is a special
integral kernel related to the Dirac’s delta.

$\mathcal{D}G(x,y)=\delta(x-y)$

As an application of Fredholm’s theory, the solution of Equation (16) can be
written in the following equivalent integral form:

$f(x)=\int_{M}G(x,y)g(y)d\omega.$

For operators with discrete spectrum and a complete orthonormal basis of
eigenfunctions $\\{\varphi_{i}\\}$ the Green’s function can be expressed:

$G(x,y)=\sum_{i}\frac{1}{\lambda_{i}}\varphi_{i}(x)\varphi_{i}(y).$ (17)

If we consider a sampling of $n$ points $x_{i}\in M$ and choose a conventional
order for the sampled points, a finite sampling of a function $f$ can be
represented by a real n-dimensional vector. Denote by $G(x,y)$ the Green’s
function of a linear differential operator $\mathcal{D}$ on $M$. In line with
Mihaylov & Spallanzani (2016), we call a _coarse discretisation_ of
$\mathcal{D}$, the linear map $D:\mathbb{R}^{n}\longrightarrow\mathbb{R}^{n}$
defined by:

$DG(x_{i},x_{j})=\delta_{ij},$ (18)

where $\delta_{ij}$ denotes Kronecker’s symbol. With this definition, the
continuous eigenfuction problem $\mathcal{D}f=\lambda f$ is discretised into

$\sum_{i}G(x_{i},x_{j})f(x_{j})=\mu f(x_{i}),$ (19)

The Green’s function, measures the way in which two points are geometrically
related. In the above problem, the contribution of the sampling $f(x_{j})$ to
the computation of $f(x_{i})$ can be further localised by introducing an
adjacency matrix $W_{ij}$ which provides the sampling $x_{i}\in M$ with a
graph structure. The elements of $W_{ij}$ are equal to 1 or 0 based on
different criteria depending on the problem, a distance-based cutoff, a fixed
number of neighbouring points (pixels), etc. Equation (19) becomes:

$Df(x_{i})=\sum_{jk}G(x_{i},x_{j})W_{jk}f(x_{k})=\mu f(x_{i}).$ (20)

Self-adjoint operators give rise to symmetric discretisations and the discrete
version of Equation (17) is:

$G(x_{m},x_{n})=\sum_{i}\frac{1}{\lambda_{i}}\varphi_{i}(x_{m})\varphi_{i}(x_{n}).$

This equation and Equation (18) mean that $D$ is the generalised inverse of
$G$ (or $GW$). By construction, if we increase the density of the sampling as
$n\longrightarrow\infty$, the coarse discretisation converges to the
differential operator.

The coarse discretisation of a known differential operator on a given manifold
is a conceptually straightforward process and it differs from the typical
manifold learning problem of approximating a linear differential operator
$\mathcal{D}$ on a point cloud sampled from a priori unknown sub-manifold of
$\mathbb{R}^{n}$. We call such a linear map the sample discretisation of
$\mathcal{D}$.

The elements used for building a sample discretisation are the same as in
Equation (20) i.e. an affinity graph and a kernel function $K(x_{i},x_{j})$
that satisfy certain symmetry and regularity assumptions (see Singer & Wu
(1997); Berry & Sauer (2016)). In the manifold learning setup the kernel
$K(x_{i},x_{j})$ reproduces the expression of the (unknown) Green’s function
of the manifold in terms of the coordinates of the points $x_{i}$ in the space
of measured variables. A first approximation attempt regarding the expression
of the kernel can be found between the Green’s function of the same operator
on $\mathbb{R}^{n}$ computed with the embedding coordinates $x_{i}$.

This estimate must be corrected by bias terms that include for example the
curvature on $M$, the sample variance etc. Intuitively, the bias of the
geodesic distance on $M$ with respect to the Euclidean distance in the ambient
space is measured by the curvature.

Given a sample discretisation of a differential operator, two types of
convergence problems arise in the limit of infinite sampling:

  * •

point convergence $D\varphi(x_{i})\longrightarrow\mathcal{D}\varphi(x)$

  * •

spectral - eigenvectors of $D$ converge to eigenfunctions $\mathcal{D}$.

Spectral convergence is a stronger condition. Proving the convergence of a
discretised operator is a non-trivial analytical problem that includes the
curvature of the manifold, the density of sampling, normalization choices etc.
Remarkable results in this direction have been achieved in Belkin & Niyogi
(2003); Singer & Wu (1997).

###  A.2 Discrete Laplacian

For a point cloud in $\mathbb{R}^{n}$ we define a graph with adjacency matrix
$A_{ij}$ with weights $K(x_{i},x_{j})$ and compute the (weighted) diagonal
degree matrix $D_{i}=\sum_{j}A_{ij}$. There are several definitions of Laplace
operators on graphs that use these elements. Typically the adopted expressions
for $K(x_{i},x_{j})$ are inspired by the Laplacian Green’s function in the
ambient space of the point cloud. For example the Green’s functions of the
Laplacian in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$ are well known:

$G_{\mathbb{R}^{2}}(x,y)=\frac{1}{2\pi}\ln(|x-y|),\,\,\,\,\,\,\,\,\,\,G_{\mathbb{R}^{3}}(x,y)=\frac{1}{4\pi|x-y|}.$

Berry & Sauer (2016) studied more general kernel constructions that take into
account non-uniform sampling densities.

The most popular construction is the graph Laplacian

$L:=D-A.$ (21)

This definition arises from the Newton’s cooling law on a graph. The heat
transferred between connected nodes $i$ and $j$ is proportional to the
difference $f(x_{i})-f(x_{j})$ and a heat capacity coefficient $k$, i.e.

$\displaystyle\frac{df(x_{i})}{dt}=$
$\displaystyle-k\sum_{j}A_{ij}(f(x_{i})-f(x_{j}))=$
$\displaystyle-k\left(f(x_{i})\sum_{j}A_{ij}-\sum_{j}A_{ij}f(x_{j})\right)$

This is further motivation for the degree matrix $D$ in Equation (21).

The random walk Laplacian is defined as

$L_{rw}:=D^{-1}L=I-D^{-1}A$ (22)

and is related to Markov processes on graphs. The matrix $L_{rw}$ is not
symmetric, but row-stochastic (satisfies the Markov property). This definition
has been further developed in the theory of diffusion maps (see Coifman &
Lafon (2006); Singer & Wu (1997) etc).

The symmetric random walk Laplacian is:

$L_{sym}:=D^{-\frac{1}{2}}LD^{-\frac{1}{2}}=I-D^{-\frac{1}{2}}AD^{-\frac{1}{2}}.$
(23)

The relations between these definitions and the transformations that map their
spectral structures respectively are well-known.

Remarkably, $L_{sym}$ matrix approximates the Laplace-Beltrami operator both
in pointwise and spectral sense (see Belkin & Niyogi (2003) and Singer & Wu
(1997)). More precisely for $n\longrightarrow\infty$,

$Lf(x_{i})\approx\Delta f(x)+O\left(\frac{||\nabla
f(x)||}{\sqrt{n}\epsilon^{1/2+dim(M)/4}}\right),$

where $\epsilon$ is a parameter (see Berry & Sauer (2016)). Moreover, in the
case of uniform sampling, Belkin & Niyogi (2003) have shown that the
eigenvectors of the graph Laplacian converge to the eigenfunctions of the
Laplace-Beltrami operator on the manifold.

Constructions of discrete Laplacian operators from non-uniform distributions
of manifold have been developed. Sample density can be absorbed in the
definition of the kernel itself or a “right normalisation” is introduced Berry
& Sauer (2016). These constructions are less relevant for the purposes of
image segmentation where images are usually sampled on regular girds of
pixels.

The geometric relation between kernels and Riemannian metrics has been
investigated in Berry & Sauer (2016) for uniform and nonuniform sampling. The
expression of the Laplace-Beltrami operator changes if the Riemannian metric
on $M$ varies in a family of metrics, which includes the one induced by the
embedding $M\hookrightarrow\mathbb{R}^{n}$. These changes are captured by
modifications in the form of the Green’s function of the operator and
subsequently in its discretisation. Every symmetric local kernel with
exponential decay corresponds to a Laplacian operator in a Riemannian geometry
and vice versa.

Once a discrete Laplacian is provided, the discrete Dirichlet energy
functional of Equation (15) applied on a sampling $f(x_{i})$ becomes:

$E[f(x_{i})]=\frac{1}{2}\langle f(x_{i}),Lf(x_{i})\rangle.$

The construction of low-energy embedding of a graph in $\mathbb{R}^{m}$ with
$m$ lower or equal to the number of sampled points is directly applied in the
discrete case and is used in spectral image segmentation.

###  A.3 Spectral segmentation by means of the heat propagator

The geometric affinity/similarity between points on a given manifold can be
measured by a diffusion process. Eq.(27) is the heat diffusion differential
equation on a smooth manifold.

The Laplace-Beltrami-based construction described above can be replicated by
by means of a closely related differential operator
$e^{-t\Delta_{g}}:L^{2}(M)\longrightarrow L^{2}(M)$ called (for $t>0$) the
_heat propagator_. Fredholm theory in this case is applied through a Green’s
function called the _heat kernel_ $G(x,y,t)$. Similarly to the Schrödinger
equation, the heat kernel is a special solution of

$\Delta_{g}G(\mathbf{x},\mathbf{y},t)=\partial_{t}G(\mathbf{x},\mathbf{y},t),\,\,\,\text{
such that
}\,\,\,G(\mathbf{x},\mathbf{y},t)=e^{t\Delta_{g}}\delta(\mathbf{x}-\mathbf{y}).$

The heat kernel is used to solve in integral form the diffusion equation with
a source term $f(x)$ and with boundary conditions:

$u(x,t)=e^{-t\Delta_{g}}f(x)=\int_{M}G(x,y,t)f(y)d\omega_{g}$ (24)

The heat propagator inherits from the Laplace-Beltrami operator the properties
of being a self-adjoint, positive-definite and compact. Its eigenfunctions
form a complete orthonormal basis of $L^{2}(M)$. The eigenvalue problem for
the heat propagator becomes:

$e^{-t\Delta_{g}}\varphi_{i}=\beta_{i}^{t}\varphi_{i}$ (25)

Let us set $\lambda_{i}:=-ln(\beta_{i})$. By construction
$e^{-t\Delta_{g}}\varphi_{i}$ satisfies the heat diffusion equation for all
$i$. Substituting Equation 25 we get:

$0=L[e^{-t\Delta_{g}}\varphi_{i}]=e^{-\lambda_{i}t}(\Delta_{g}\varphi_{i}-\lambda_{i}\varphi_{i}).$

So $\varphi_{i}$ are precisely the eigenfunctions of the Laplace-Beltrami
operator with eigenvalues $\lambda_{i}$. As a consequence, the relevant
properties that enable image segmentation through the Laplace-Beltrami
operator (eigenfunctions define a low Dirichlet energy embedding of $M$ in
$\mathbb{R}^{n}$) can be directly extended to the heat propagator. In
particular

$G(x,y,t)=\sum_{i}e^{-\lambda_{i}t}\varphi_{i}(x)\varphi_{i}(y).$ (26)

For instance, the heat kernel in $\mathbb{R}^{3}$ is:

$G(\mathbf{x},\mathbf{y},t)_{\mathbb{R}^{3}}=\left(\frac{1}{4\pi
t}\right)^{3/2}\exp\left(-\frac{\parallel\mathbf{x}-\mathbf{y}\parallel^{2}}{4t}\right)\,\,\,\text{,}\,\,\,G(\mathbf{x},\mathbf{y},0)=\delta(\mathbf{x}-\mathbf{y})$

An argument provided in Belkin & Niyogi (2003) derives Equation (21) by
discretising the heat propagator. By substituting in 24 in 27 we get:

$\displaystyle\Delta_{g}f(x)=\Delta_{g}u(x,0)=$
$\displaystyle\frac{\partial}{\partial
t}\left(\int_{M}G(x,y,t)f(y)d\omega_{g}\right)_{t=0}$ $\displaystyle\approx$
$\displaystyle-\frac{1}{t}\left(f(x)-\int_{M}G(x,y,t)f(y)d\omega_{g}\right)$
$\displaystyle\approx$
$\displaystyle-\frac{1}{t}\left(f(x_{i})-\sum_{j}G(x,y,t)f(x_{j})\right).$

The overall coefficient $1/t$ does not affect the spectrum of the expression
in the brackets. This argument is in line with the above discussion on the
heat propagator and justifies further the diagonal term in the discrete
definitions of the discrete Laplacian based on diffusion arguments.

Pointwise and spectral convergent discretizations of the heat propagator have
been developed in the manifold learning theory of diffusion maps ( Coifman &
Lafon (2006); Singer & Wu (1997); Berry & Sauer (2016)). Given the functional
form form of the heat kernel in $\mathbb{R}^{n}$ this class of discrete
constructions of the Laplace-Beltrami operator predominantly exploit global
and local Gaussian kernels.

###  A.4 Heat-kernel in $S^{2}$

As discussed in Sect.3, we choose the adjacency weights to be the ones derived
from the functional form of the integral kernel of the Laplace-Beltrami
operator in $S^{2}$. In this section, we derive the functional form in eq.
(2).

First, we write the heat equation in spherical coordinates:

$(\partial_{t}-\Delta_{g})u(x,t)=L[u(x,t)]=\rho(x,t),$ (27)

with $\Delta_{g}$ being the Laplace-Beltrami operator in $S^{2}$, $u$ the
unknown function describing the heat diffusion, $\rho$ the source term, $x$
and $t$ respectively the spatial and time variables. We can express
$\Delta_{g}$ in spherical coordinates (assuming a radius $r=1$):

$\Delta\equiv\frac{1}{\sin\theta}\frac{\partial}{\partial{\theta}}\left(\sin\theta\frac{\partial}{\partial{\theta}}\right)+\frac{1}{\sin^{2}\theta}\frac{\partial^{2}}{\partial{\phi^{2}}}.$
(28)

Expressed in this way, we can easily recognize that the Laplace-Beltrami
operator corresponds to $-\hat{L}^{2}$, i.e. the square of the _orbital
angular momentum_ operator in quantum mechanics, whose eigenfunctions are the
_spherical harmonics_ , $Y_{\ell m}$, and whose eigenvalues are
$-\ell(\ell+1)$ ($\ell\in\mathbb{N}$). We further notice that the heat
propagator operator, $\exp{(-\hat{L}^{2}t)}$ shares the same eigenfunctions
and eigenvalues as $-\hat{L}^{2}$. In Appendix A, we show how the integral
kernel of the heat equation is related to the eigenfunctions of
$\exp{(-\hat{L}^{2}t)}$ operator (see 26), we can thus derive the kernel in
spherical coordinates151515Assuming as an initial condition to the source
term, $\rho$ in (27) to be a Dirac’s delta distribution, i.e.
$\rho(x,0)=\delta(x).$:

$\displaystyle G(\mathbf{x},\mathbf{y},t)=$
$\displaystyle\sum_{\ell=0}^{+\infty}e^{-t\ell(\ell+1)}\sum_{m=-\ell}^{+\ell}Y_{\ell
m}(\mathbf{x})Y_{\ell m}^{*}(\mathbf{y})$ $\displaystyle=$
$\displaystyle\sum_{\ell=0}^{+\infty}\frac{2\ell+1}{4\pi}e^{-t\ell(\ell+1)}\mathcal{P}_{\ell}(\Theta_{\mathbf{xy}}),$
(29)

where we exploit in the last equality the properties of spherical harmonics,
see further details in Appendix A. In particular, the last equality is
obtained by expressing the $Y_{\ell m}$ as a function of _Legendre
polynomials_ , $\mathcal{P}_{\ell}$, evaluated across the _cosine matrix_ ,
defined as :

$\Theta_{\mathbf{xy}}\equiv\cos\theta_{\mathbf{xy}}=\frac{\mathbf{x}\cdot\mathbf{y}}{\parallel\mathbf{x}\parallel\parallel\mathbf{y}\parallel}.$
(30)

We devote Sect. 3 to describe how to derive the adjacency weights from eq.
(29), so that we can construct the Laplacian matrix as in eq. (1) and estimate
the eigenpairs to perform the spectral clustering.

###  A.5 Laplacian eigenfunctions and the Ricci curvature

There are two conceptually distinct ways in which eigenfunctions of the
Laplace-Beltrami operator and their discrete analogues are exploited for image
processing.

\- An image is interpreted as a function on a regular domain (typically a
rectangle in $\mathbb{R}^{2}$). This function is decomposed with respect to
the special basis of eigenfunctions. This decomposition provides a useful
characterisation of an image, allows highlighting or filtering global
components with specific Dirichlet energy i.e. spatial “frequency”. As an
elementary example the coefficient of $\varphi_{0}$ captures an overall
intensity. This approach, similar to spectral shape analysis, has been adopted
in M. Reuter & Peinecke (2005) and Rustamov (2007) and is not typically used
for segmentation purposes.

\- An image itself is considered as smooth surface (real 2-dimensional
manifold) sampled in a discrete set of pixels. The geodesic distance between
points on the manifold combines the distance between pixels with the the
“vertical deformation” that encodes the image. Spectral geometry is a field of
differential geometry that exploits the spectral structure of Laplace-Beltrami
operator to describe relevant geometric and topological properties of
manifolds. The spectral structure of the discrete Laplacian approximates the
spectral structure of the Laplace-Beltrami operator. This is the manifold
learning approach to spectral segmentation.

The application of ideas from geometry and theory of partial differential
equations to image segmentation is an active field of research. Variational
calculus based on a variety of energy functionals (Mumford-Shah, Chan-Vese,
elastic energy with additional forcing terms and well potentials) and dynamics
on different time scales have been developed in H. Hu & Bertozzi (2015); Meng
et al. (2017); Boyd et al. (2018). The deep geometric nature of the spectral
clustering methods is yet to be investigated. In fact the formulation of
spectral image segmentation in the smooth case, means that points in plateau-
ing regions of the manifold are close to level sets of the eigenfunctions of
the Laplace-Beltrami operator in a suitable Dirichlet energy band. Relevant
properties of the nodal (and more general level) sets of the eigenfunctions of
$\Delta_{g}$ have been studied in special geometric cases (constant curvature
or curvature bounded from below). Less is know on how the level sets of
eigenfunctions capture the geometry of a manifold in the general case.

Equivalently, we expect that eigenfunctions of the Laplace-Beltrami operator
with given Dirichlet energy are characterised by small variations on plateau-
ing regions of the manifold. The key question is how can we characterise
platteau-ing regions.

Our hypothesis is that these are the regions in which the Ricci curvature of
the manifold is low (either positive or negative). A version of the celebrated
Bochner’s formula establishes a direct relation between the Ricci’s curvature
tensor, the Laplace-Beltrami operator and the norm of the Hessian of a smooth
function $\varphi$ on a Riemannian manifold $M$ with metric $g$:

$\frac{1}{2}\Delta_{g}|\nabla\varphi|^{2}=g(\Delta_{g}\nabla\varphi,\nabla\varphi)+|H(\varphi)|^{2}+Ric(\nabla\varphi,\nabla\varphi).$

For eigenfunctions $\Delta_{g}\varphi=\lambda\varphi$ the Bochner’s formula
becomes:

$Ric(\nabla\varphi,\nabla\varphi)=\frac{1}{2}\Delta_{g}|\nabla\varphi|^{2}-\lambda|\nabla\varphi|^{2}-|H(\varphi)|^{2}$

For smooth surfaces (images) this equation means that the sectional curvature
in the direction of the maximum variation of the eigenfunction depends
directly on the magnitude of the gradient of the eigenfunction and the
eigenvalue.

The above considerations suggest the existence of a deep geometric reason for
the fact that spectral image segmentation is very efficient in detecting
objects with complicated shapes in unsupervised setup. Understanding spectral
image segmentation in the context of spectral geometry is a promising field of
future research.

##  Appendix B Clustering applied on _Planck_ -GNILC dust maps

Figure 15: Variance surfaces estimated for several choices of $\alpha$ and
$\delta$ for (left) $\beta_{d}$, (right) $T_{d}$. Notice that regions of
optimality corresponds to local minima (lighter regions) in the surfaces.
Figure 16: Cluster maps of (left) $\beta_{d}$ and (right) $T_{d}$ GNILC maps.
Figure 17: (first and second column) Dust parameter and uncertainty maps from
_Planck_ -GNILC. (third column) Dust spectral parameter maps evaluated by
taking the median value of the input GNILC maps on each region defined by
spectral clustering procedure. (fourth column) Relative errors estimated by
taking the difference of GNILC maps with the binned ones.

Planck Collaboration et al. (2016e) released a set of maps obtained by
separating the Galactic thermal dust emission from the cosmic infrared
background (CIB) anisotropies. They implemented a tailored component-
separation method, the Generalized Needlet Internal Linear Combination (GNILC)
method, that adopts spatial information (e.g. the angular power spectra) to
separate the Galactic dust emission and CIB anisotropies. We thus considered
the temperature $T_{d}$ and spectral index $\beta_{d}$ maps of thermal dust
and the uncertainties accompanying these maps to partition the sky into
multiple domains accounting both for the geometrical similarities and the
uncertainties.

Given the fact that we are interested to derive patches at the $\sim$ degree
scales, we firstly reduce the resolution of the GNILC maps (released at
nside=2048 and $5^{\prime}$ resolutions), to a coarser beam, $110^{\prime}$
and to a lower pixel resolution, nside=32. This also makes the Laplacian
matrix computation faster with less memory requirements.

We then performed the clustering following the Algorithm 1 separately for
$T_{d}$ and for $\beta_{d}$. In Fig.15, we show the variances estimated from
the GNILC maps in an $\alpha-\delta$ grid. We note the presence of a minimum
variance vertical stripe at around $\delta\sim 0.2$ and for $\alpha\leq 0.2$.
Notice that high variance regions can be identified in the upper right and
lower left corners of the image related respectively to over and under
partition regimes. Moreover, ranges of optimality both for $\beta_{d}$ and
$T_{d}$ present similar variance contours. Since the optimization is performed
independently for the two dust parameters, we interpret this result as an
indication that the methodology might be optimizing for similar Galactic scale
features in both $\beta_{d}$ and $T_{d}$ maps.

We find the optimal number of clusters to be: 1300 and 1415 respectively for
$\beta_{d}$ and $T_{d}$ shown in Fig. 16, corresponding to the parameter
values of $\alpha=0.15,\delta=0.20$.

Finally, we want to assess how well the optimal regions track the observed
features in the GNILC maps. We thus estimate within each cluster patch, the
median value of the pixels identified in a single cluster and build a _binned_
map for $T_{d}$ and $\beta_{d}$, labeled as $\tilde{T}_{d}$ and
$\tilde{\beta}_{d}$. We show in Fig.17 the maps from GNILC and the binned
maps. To quantify the error of this procedure, we estimate the relative
difference of the input and binned maps as in eq. (9).

Relative error maps are shown in the right-most column of Fig.17. We note that
the errors are lower than the $10\%$.

##  Appendix C Defining the optimality range

Figure 18: (top) Cluster patches for $\beta_{d}$ GNILC maps, in regions of
optimality as shown in Fig. 15, for several choices of $\alpha$ and $\delta$.
(bottom) Maps showing the logarithm of the number of pixels within each
cluster shown in the top row.

From the variances shown in Figures 15, 6 and 5, we notice that the range of
optimality corresponds to a narrow range in $\delta$ but seldom involves a
broader range in the $\alpha$ parameter.

In the top row of Fig. 18, we show several cluster configurations chosen in
the range of optimality from the variance in Fig. 15(left). We quantify the
size of the cluster by estimating how many pixels belong to each cluster. The
logarithm of the number of pixels per cluster is shown in the bottom row of
Fig. 18. Values of $\alpha<0.1$, lead to a very homogeneous and isotropic
partition which is intuitively expected as we partition the celestial sphere
given the properties of the symmetric Laplacian adjacency (eq. (2)). On the
other hand, increasingly larger values of $\alpha>0.2$ make the cluster
morphology more and more dictated by the SNR content of the features,
resulting in inhomogenuous sizes and anisotropic shapes in the vicinity of the
Galactic plane as the manifold adjacency gets more and more distorted (see eq.
(4). We find that intermediate values for $\alpha$ in the range around
$0.1\div 0.2$ result in a good trade off between these two extreme cases.

##  Appendix D Estimation of uncertainties for the $\mathcal{N}_{c}$ map

As described in Section 4.2, we apply our clustering algorithm to partition
the sky in regions of near-constant number of clouds per line of sight, as
determined by the analysis of HI data in Panopoulou & Lenz (2020). A critical
component of the clustering analysis is its ability to handle measurement
uncertainties. The public map of $\mathcal{N}_{c}$ does not contain
uncertainties and so we repeat steps of the original analysis to obtain
estimates of the per-pixel uncertainties in the map.

We compute the uncertainty on $\mathcal{N}_{c}$ by taking into account the two
sources of error that enter in the calculation of $\mathcal{N}_{c}$:

  1. 1.

the uncertainty on the number of velocity peaks identified in the HI spectrum,

  2. 2.

the uncertainty on the column density of each cloud.

As explained in Panopoulou & Lenz (2020, Appendix B1), both sources of
uncertainty are driven by the choice of the velocity kernel size (a parameter
that implements smoothing in velocity space, also referred to as the
bandwidth). The choice of bandwidth presents a trade-off between resolving
power (ability to distinguish nearby peaks in velocity) and fidelity (avoiding
spurious peak detections). Furthermore, the kernel size also alters the shape
of the probability distribution, including the position of extrema and saddle
points, used to define the range of each peak.

We therefore estimate the uncertainties by considering maps of
$\mathcal{N}_{c}$ with three different choices of bandwidth: namely 3, 4 and 5
velocity channels, indicated respectively as $\mathcal{N}_{c}^{\mathrm{bw3}}$,
$\mathcal{N}_{c}^{\mathrm{bw4}}$, $\mathcal{N}_{c}^{\mathrm{bw5}}$, with
$\mathcal{N}_{c}\equiv\mathcal{N}_{c}^{\mathrm{bw4}}$ being the default value
in Panopoulou & Lenz (2020). We evaluated the relative differences of
$\mathcal{N}_{c}$ between different bandwidth runs as:

$\Delta_{rel}^{i,4}=(\mathcal{N}_{c}^{\mathrm{bw\,i}}-\mathcal{N}_{c}^{\mathrm{bw4}})/\mathcal{N}_{c}^{\mathrm{bw4}}.$

These differences can be used to quantify the systematic uncertainty of the
cloud identification method, and hence $\mathcal{N}_{c}$.

Figure 19 shows the joint distribution of $\Delta_{rel}^{4,3}$ and
$\Delta_{rel}^{5,4}$ for the map at nside=32. The majority of pixels show
relative differences of $<0.1$ (10%). The 16–percentile and 84–percentiles of
the distribution of $\Delta_{rel}^{4,3}$ are -0.3 and 0.01, respectively. For
$\Delta_{rel}^{5,4}$, these values are -0.1 and 0.01. There are long tails
extending out to -2 for the distribution of $\Delta_{rel}^{4,3}$ and -0.6 for
that of $\Delta_{rel}^{5,4}$. We interpret these distributions as follows: by
using a smaller bandwidth, a larger number of clouds is detected (due to the
increased velocity resolution). The bulk of the observed differences (within
10—30% relative difference) are likely related to changes in the
identification of low-column density clouds. These clouds are prevalent in the
HI decomposition, as discussed in Panopoulou & Lenz (2020), and contribute a
low-level noise to the determination of $\mathcal{N}_{c}$. We will assign this
low-level noise as a floor in the uncertainty on $\mathcal{N}_{c}$. The longer
(negative) tails of the distributions are likely related to edge cases, where
clouds of significant column density were unresolved at low bandwidth but
become resolved with a small increase in the velocity resolution.

Because our clustering algorithm deals with Gaussian uncertainties, we need to
translate these asymmetric differences into an equivalent standard deviation.
We choose to err on the side of overestimating the uncertainties, by adopting
the following estimate of the uncertainty on $\mathcal{N}_{c}$:

$\sigma(\mathcal{N}_{c})=\max\left(0.3,\,\Delta_{rel}^{3,4},\,\Delta_{rel}^{5,4}\right),$
(31)

meaning that we assign the maximum observed difference between different
bandwidth runs as the (symmetric) $\rm 1-\sigma$ uncertainty in each pixel.
The aforementioned low-level noise contributed by faint clouds represents the
uncertainty floor that we assign (and corresponds to the absolute value of the
16-percentile of $\Delta_{rel}^{3,4})$.

Figure 19: Distribution of asymmetric uncertainties on $\mathcal{N}_{c}$ at
nside=32, using the relative difference between $\mathcal{N}_{c}$ measured at
different values of the bandwidth parameter. Vertical axis: difference between
bandwidths of 3 and 4, normalized to the $\mathcal{N}_{c}^{\mathrm{bw4}}$.
Horizontal axis: same for bandwidths 5 and 4. The dashed lines mark the 16-
and 84- percentiles containing 68% of the distribution of $\Delta_{rel}^{3,4}$
and $\Delta_{rel}^{5,4}$.