# Minkowski Functionals in SO(3) for the spin-2 CMB polarisation field

J. Carrón Duque, A. Carones, D. Marinucci, M. Migliaccio, N. Vittorio

► To cite this version:

J. Carrón Duque, A. Carones, D. Marinucci, M. Migliaccio, N. Vittorio. Minkowski Functionals in SO(3) for the spin-2 CMB polarisation field. Journal of Cosmology and Astroparticle Physics, 2024, 2024 (01), pp.039. 10.1088/1475-7516/2024/01/039. [hal-04863136]

HAL Id: hal-04863136

https://hal.science/hal-04863136v1

Submitted on 3 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Minkowski Functionals in (3) for the spin-2 CMB polarisation field

To cite this article: J. Carrón Duque et al JCAP01(2024)039

View the article online for updates and enhancements.

## You may also like

- The VLA/ALMA Nascent Disk And Multiplicity (VANDAM) Survey of Orion Protostars. V. A Characterization of Protostellar Multiplicity
John J. Tobin, Stella S. R. Offner, Kaitlin M. Kratter et al.

- Enhancing Morphological Measurements of the Cosmic Web with Delaunay Tessellation Field Estimation
Yu Liu, Yu Yu, Pengjie Zhang et al.

- STATISTICAL AND SYSTEMATIC ERRORS IN THE MEASUREMENT OF WEAK-LENSING MINKOWSKI FUNCTIONALS: APPLICATION TO THE CANADA-FRANCE-HAWAII LENSING SURVEY
Masato Shirasaki and Naoki Yoshida
# Minkowski Functionals in SO(3) for the spin-2 CMB polarisation field

J. Carrón Duque,$^{a,b,*}$ A. Carones,$^{a,b}$ D. Marinucci,$^c$ M. Migliaccio$^{a,b}$ and N. Vittorio$^{a,b}$

$^a$Dipartimento di Fisica, Università di Roma “Tor Vergata”, via della Ricerca Scientifica 1, Roma I-00133, Italy

$^b$INFN, Sezione di Roma 2, via della Ricerca Scientifica 1, Roma I-00133, Italy

$^c$Dipartimento di Matematica, Università di Roma Tor Vergata, via della Ricerca Scientifica 1, Roma I-00133, Italy

E-mail: javier.carron@roma2.infn.it, alessandro.carones@roma2.infn.it, marinucc@mat.uniroma2.it, marina.migliaccio@roma2.infn.it, nicola.vittorio@uniroma2.it

Received February 13, 2023

Revised September 26, 2023

Accepted November 30, 2023

Published January 17, 2024

**Abstract.** The study of the angular power spectrum of Cosmic Microwave Background (CMB) anisotropies, both in intensity and in polarisation, has led to the tightest constraints on cosmological parameters. However, this statistical quantity is not sensitive to any deviation from Gaussianity and statistical isotropy in the CMB data. Minkowski Functionals (MFs) have been adopted as one of the most powerful statistical tools to study such deviations, since they characterise the topology and geometry of the field of interest. In this paper, we extend the application of MFs to CMB polarisation data by introducing a new formalism, where we lift the spin 2 polarisation field to a scalar function in a higher-dimensional manifold: the group of rotations of the sphere, SO(3). Such a function is defined as $f = Q \cos(2\zeta) - U \sin(2\zeta)$. We analytically obtain the expected values for the MFs of $f$ in the case of Gaussian isotropic polarisation maps. Furthermore, we present a new pipeline which estimates these MFs from input HEALPix polarisation maps. We apply it to CMB simulations in order to validate the theoretical results and the methodology. The pipeline is to be included in the publicly available Python package *Pynkowski*.

**Keywords:** CMBR polarisation, non-gaussianity

**ArXiv ePrint:** 2301.13191

<footnote>*Corresponding author.</footnote>
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Spin field as a scalar field in SO(3)</td><td>2</td></tr><tr><td>3 Minkowski Functionals in SO(3)</td><td>4</td></tr><tr><td>3.1 Definition</td><td>5</td></tr><tr><td>3.2 Gaussian Kinematic Formula</td><td>6</td></tr><tr><td>3.3 Theoretical predictions</td><td>8</td></tr><tr><td>4 Implementation: Pynkowski</td><td>11</td></tr><tr><td>4.1 General considerations</td><td>14</td></tr><tr><td>5 Simulations</td><td>14</td></tr><tr><td>5.1 Monochromatic maps</td><td>14</td></tr><tr><td>5.2 Realistic angular power spectrum</td><td>14</td></tr><tr><td>6 Results</td><td>15</td></tr><tr><td>6.1 Normalisation constants</td><td>15</td></tr><tr><td>6.2 Scaling relations</td><td>15</td></tr><tr><td>6.3 Gaussian CMB simulations</td><td>16</td></tr><tr><td>7 Conclusions</td><td>17</td></tr><tr><td>A Polarization as a random field on SO(3)</td><td>19</td></tr><tr><td>B Metric and derivatives in SO(3)</td><td>20</td></tr></table>

## 1 Introduction

The Cosmic Microwave Background (CMB) encodes information from the Early Universe, both in the intensity and polarisation of the light. The CMB polarisation field is usually decomposed into two distinct rotationally-invariant fields: the *E* and *B* modes [see 1, 2, for details]. These fields are commonly studied through their angular power spectra (equivalently, 2-point correlation functions). However, this tool is not sensitive to the possible presence of non-Gaussianities or departures from statistical isotropy of the CMB anisotropy fields.

Non-Gaussianity is predicted by many inflationary models [3–5] and could shed new light on our knowledge of the primordial Universe. There is also a growing amount of literature on a possible large-scale anisotropy of the Universe, with dipoles being measured in several observables. Furthermore, the CMB maps contain foregrounds contamination because of Galactic emission and the lensing of CMB photons due to their interaction with the Large Scale Structure. These effects significantly deviate from the hypothesis of Gaussianity and isotropy, and thus have to be carefully considered when analysing the data. These effects are especially important in CMB polarisation.

Minkowski Functionals (MFs) are one of the tools adopted by the Cosmology community to study possible deviations from Gaussianity or statistical isotropy. These functionals
encode geometrical and topological information of the field, not reflected in the power spec-
tra. Other tools include the bispectrum and trispectrum, or, equivalently, the 3- and 4-points
correlation functions [3, 6, 7], the distribution of maxima and minima [8–10], or the distribu-
tion of nonpolarised points in polarisation fields [11]. MFs present several advantages with
respect to the bispectrum and the trispectrum, such as the computational cost, the ease of
masking or weighting data, and the possibility of studying deviations at different thresholds.
The last one makes MFs naturally suited to study non-Gaussianities that are not optimally
expressed in terms of momenta expansion ($f_{NL}, g_{NL}, ...$); this is the case, for example, in
inflationary models that can produce primordial black holes, such as Stochastic Inflation, as
this introduces non-Gaussianity mostly at high values of the field [12, 13].

The application of MFs has been mostly limited to scalar maps so far, such as CMB
temperature [14, 15], and weak lensing [16–18]. They have also been used to study the
morphological properties of Galactic emission, like thermal dust [19], and synchrotron [20, 21].
However, the CMB polarisation field is a complex spin 2 quantity and MFs have not been
defined for this kind of maps. In polarisation studies, MFs are usually applied to the *E* and
*B* scalar maps independently [15, 22, 23], or directly to the *Q* and *U* maps, ignoring spin
effects [24–26].

In a previous work [27], we focused on the application of MFs to the squared polarised
intensity of the CMB, $P^2 = Q^2 + U^2$. We introduced the formalism and computed the
theoretical expectations in the Gaussian isotropic case by making use of the Gaussian Kine-
matic Formula. We also developed a Python package to estimate the MFs on arbitrary
HEALPix scalar maps and compare them with the theoretical predictions; this software,
called *Pynkowski*, is now publicly available.<fnref n="1" />

In this work, we introduce a new theoretical framework to analyse the full information
of the polarisation data, i.e., without limiting the analysis to scalar quantities defined on the
sphere such as *P*, or the *E* or *B* modes. Such an approach provides more complete information
of the statistics of CMB polarisation data, while also avoiding leakage contamination in the
decomposition of masked *Q* and *U* maps into *E* and *B* modes. This is explicitly done by
lifting the polarisation field to a three-dimensional space. We will therefore use MFs on three
dimensions, for which there is a large body of literature [28–31]. However, previous works
are limited to fields defined on $\mathbb{R}^3$, while for our purposes we shall generalise this formalism
for arbitrary manifolds.

The paper has the following structure. In section 2 we introduce the field on which the MFs will be computed and we explain some of the technical aspects needed to perform such computations. In section 3 we obtain the theoretical expectations of MFs for Gaussian and statistically isotropic spin maps. In section 4 we present the pipeline to estimate the MFs on arbitrary HEALPix spin maps from the *Q* and *U* data. In section 5 we introduce the simulations we use to validate the formalism and the pipeline, while in section 6 we present the results of applying this framework and software to the aforementioned CMB polarisation simulated maps. Finally, in section 7 we summarise our conclusions.

## 2 Spin field as a scalar field in SO(3)

The CMB temperature anisotropies map can be seen as a real scalar field defined on the
sphere. Thus, its statistical properties can be analysed with plenty of tools, such as the MFs
formalism, first introduced in the Cosmological literature in [14]. This tool is used to describe

<footnote n="1">¹https://github.com/javicarron/pynkowski.</footnote>
several characteristics of the excursion sets of the fields at different thresholds, which define their geometry and topology. However, the CMB polarisation has a different geometrical structure: it constitutes a complex spin 2 field on the sphere [2], for which excursion sets cannot be directly defined. To overcome this issue, we lift the field to a higher-dimensional space where it can be seen as a scalar field, following the framework introduced in [32]. See also [33] for further mathematical discussion on spin random fields.

Let $Q(\phi, \theta)$ and $U(\phi, \theta)$ be the maps of Stokes parameters for linear polarisation in the usual base. We define $f(\phi, \theta, \zeta) : SO(3) \to \mathbb{R}$ as:

$$f(\phi, \theta, \zeta) = Q(\phi, \theta) \cos(2\zeta) - U(\phi, \theta) \sin(2\zeta), \qquad (2.1)$$

which can be interpreted as the linear polarisation that one would observe at the point on the sky $(\phi, \theta)$ when observing along the polarisation direction $\zeta$, which is just an additional coordinate of the ambient space. Specifically, $f$ corresponds to the value of $Q$ if the local reference frame was rotated by an angle $\zeta$. We note that the value of $\zeta$ where $f$ is maximum for fixed $(\phi, \theta)$ is the physical polarisation angle $\psi$. The field $f$ is a three-dimensional scalar field, for which excursion sets and MFs are properly defined, as we will see in the next section.

The domain of this function, $SO(3)$, requires further technical discussion, to which we dedicate the remainder of this section. Keeping in mind the interpretation of the function variables $(\phi, \theta, \zeta)$ as the position and the polarisation direction, it can be seen that the domain of $f$ must cover all points of the sphere and all possible polarisation directions. The domain must then be contained in a three-dimensional hypersphere, $S^3$. Since the CMB polarisation is a spin 2 field, we can perform the identification $(\phi, \theta, \zeta) \cong (\phi, \theta, \zeta + \pi)$. Therefore, the domain can actually be seen as half a 3-sphere; this space is diffeomorphic to $SO(3)$, the set of rotations of the sphere. The coordinates are the longitude $\phi \in [0, 2\pi]$, the latitude $\theta \in [-\frac{\pi}{2}, \frac{\pi}{2}]$, and the polarisation direction $\zeta \in [0, \pi]$. We note that any parametrisation of $SO(3)$ must present singularities. In this case, it specifically fails at $\theta = \pm \frac{\pi}{2}$ (corresponding to the poles of $S^2$). However, these points constitute a zero-measure set, and therefore they do not affect the computation of the MFs, since, as discussed in the next section, they are integrated quantities.

Following this framework, we can lift the complex spin field on the sphere to a complex scalar field on $SO(3)$. The real and imaginary parts of this field are Gaussian fields and they are just translations of each other, so it is enough to study only the real part in order to characterise the geometry and topology of the polarisation field. This real part is what we have called $f$ in equation (2.1). A more detailed discussion of such lifted field can be found in appendix A, including the equivalence between its real and imaginary parts together with a discussion about their Gaussianity.

An important consequence of this construction is that if we consider an isotropic spin 2 field on the sphere, it does not constitute an isotropic field on $S^3$ nor $SO(3)$. Physically, this can be seen as a consequence of the different behaviour of the polarisation direction coordinate and the sky coordinates. In particular, the function $f$ has a deterministic (sinusoidal) dependence on the $\zeta$ variable; this will have important consequences in the predictions of the MFs of this field, as we will see in section 3.2. The treatment of this anisotropy in the MF formalism represents one of the main novelties introduced in this work. Mathematically, an isotropic random field on $S^3$ has to be invariant in law to the action of any isometry of $S^3$; it can be proven that this is only satisfied if every multipole component of the field is a sum of fields with spin $s = -\ell, -\ell+1, \dots, \ell-1, \ell$, each with equal power. This cannot be the case for fields produced by lifting a spin field, such as CMB polarisation, since all multipoles are
constituted only by $s = \pm 2$. The details of this construction and the consequence on random fields can be found in [32]; some statistical properties of fields where the spin increases with the multipole can be found in [33].

Although $f$ is not isotropic in SO(3), we still have enough information to produce accurate predictions for its MFs. We need simply to assume isotropy on the sphere ($S^2$) and the fact that we have a spin 2 field (i.e., knowledge of the behaviour of the polarisation direction coordinate). We will exploit these aspects in the next section.

To correctly compute the derivatives of the field needed to estimate the MFs, we have to take into account the geometry of the ambient manifold. The metric of SO(3) in this framework is given by:

$$g_{\mu\nu} = \begin{pmatrix} 2 & 0 & 2\sin(\theta) \\ 0 & 2 & 0 \\ 2\sin(\theta) & 0 & 2 \end{pmatrix}, \qquad (2.2)$$

where the order is $(\phi, \theta, \zeta)$; see appendix B for details on the computation. With this metric, the volume element and the total volume of SO(3) are:

$$dV = \cos(\theta) d\phi d\theta d\zeta, \qquad (2.3)$$

$$\int_{0}^{2\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \int_{0}^{\pi} \cos(\theta) d\phi d\theta d\zeta = 4\pi^2. \qquad (2.4)$$

A consequence of the geometry of this manifold is that the base of the tangential space given by the derivatives $\{\frac{\partial}{\partial\phi}, \frac{\partial}{\partial\theta}, \frac{\partial}{\partial\zeta}\}$ is not orthonormal. Therefore, it will be useful to work in the orthonormal basis $\{e_1, e_2, e_3\}$, built considering the square root of the inverse of the metric tensor:

$$e_1 = \frac{\sqrt{1 - \sin\theta} + \sqrt{\sin\theta + 1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\phi} + \frac{\sqrt{1 - \sin\theta} - \sqrt{\sin\theta + 1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta}, \qquad (2.5a)$$

$$e_2 = \frac{1}{\sqrt{2}} \frac{\partial}{\partial\theta}, \qquad (2.5b)$$

$$e_3 = \frac{\sqrt{1 - \sin\theta} - \sqrt{\sin\theta + 1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\phi} + \frac{\sqrt{1 - \sin\theta} + \sqrt{\sin\theta + 1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta}. \qquad (2.5c)$$

In this basis, the gradient and the Hessian of a function $f$ can be expressed as

$$\nabla f = \begin{pmatrix} e_1 \\ e_2 \\ e_3 \end{pmatrix} f, \qquad (2.6)$$

$$\mathcal{H}_{ij}(f) = \nabla^2 f(e_i, e_j) = e_i e_j f - (\nabla_{e_i} e_j) f. \qquad (2.7)$$

The expressions of these operators in terms of the usual spatial derivatives $\{\frac{\partial}{\partial\phi}, \frac{\partial}{\partial\theta}, \frac{\partial}{\partial\zeta}\}$ can be found in appendix B, together with the computations to obtain them. These results are needed to derive the MFs of arbitrary polarisation maps, as we will see in section 4.

### 3 Minkowski Functionals in SO(3)

MFs are statistical tools that quantify the morphology produced by a scalar function. They are higher-order statistics, meaning that their value cannot be fully predicted from any n-point correlation function. This renders MFs a useful complementary tool to the angular
power spectrum. Although typically used on the sphere, they can be employed on any manifold, as we show in this section for the case of SO(3).

## 3.1 Definition

Let $f: SO(3) \to \mathbb{R}$ be a $\mathcal{C}^2$ function, and $u \in \mathbb{R}$. The excursion set of $f$ at threshold $u$, called $A_u(f, SO(3))$, is defined as the collection of all points in the domain where the function is larger than (or equal to) $u$:

$$A_u(f, SO(3)) = \{x \in SO(3) : f(x) \ge u\}. \qquad (3.1)$$

In the following, we will refer to it simply as $A_u$ to simplify the notation. We will use the standard notation $\partial A_u$ to refer to the boundary of $A_u$, corresponding to the points with a value of exactly $u$.

In the case of the sphere or any 2-dimensional manifold, we can define three independent MFs to fully characterise the morphology of the excursion set, as shown in [14]. In a 3-dimensional manifold, like SO(3), we can define four independent MFs:

$$V_0(A_u) = \int_{A_u} dx, \qquad (3.2a)$$

$$V_1(A_u) = \frac{1}{6} \int_{\partial A_u} da, \qquad (3.2b)$$

$$V_2(A_u) = \frac{1}{6\pi} \int_{\partial A_u} H(a) da, \qquad (3.2c)$$

$$V_3(A_u) = \frac{1}{4\pi} \int_{\partial A_u} K(a) da, \qquad (3.2d)$$

where $da$ denotes an area element along the boundary surface $\partial A_u$, $H(a)$ and $K(a)$ are, respectively, the mean curvature (average of the two principal curvatures) and the Gaussian curvature (product of the two principal curvatures) at a point $a$ on the boundary surface. In general, in a $n$-dimensional space, we can define $V_0$ as the hypervolume of the excursion set, and $n$ independent MFs as integrals on the boundary of a combination of the principal curvatures; see [14] for more details. We follow the normalisation convention in this reference, noting that other normalisations have been used before, more notably in [34, 35].

MFs can be interpreted as geometrical descriptors of the excursion set as a function of threshold, and therefore they represent a statistical characterisation of the original function. $V_0$ is the total volume of the excursion set. $V_1$ is the area of its boundary. To our knowledge, $V_2$ does not have a further interpretation beyond the average mean curvature of the boundary of the excursion set. Finally, $V_3$ is connected to the Euler-Poincaré characteristic $\chi$ due to an extension of the Gauss-Bonnet theorem known as the Chern-Gauss-Bonnet theorem [originally proven in 36], which states:

$$\chi(A_u) = V_3(A_u) + \frac{1}{4} V_1(A_u). \qquad (3.3)$$

Similarly to [27], we introduce the Lipschitz-Killing curvatures, as the theoretical predictions are obtained for these quantities (see section 3.3), which are defined as follows. Let us consider a tube of width $\rho$ built around the manifold $A_u$ as all the points at a distance less than $\rho$ from $A_u$. Its volume can be exactly expressed as a finite Taylor expansion on $\rho$, whose
coefficients correspond to the Lipschitz-Killing curvatures of $A_u$ (see [34, 37] for a detailed description). These quantities are connected to the MFs in the following way:

$$ \mathcal{L}_3(A_u) = V_0(A_u), \qquad (3.4a) $$

$$ \mathcal{L}_2(A_u) = 3V_1(A_u), \qquad (3.4b) $$

$$ \mathcal{L}_1(A_u) = 6V_2(A_u), \qquad (3.4c) $$

$$ \mathcal{L}_0(A_u) = V_3(A_u). \qquad (3.4d) $$

In order to compute the theoretical predictions for a Gaussian field, we have to invoke the Gaussian Kinematic Formula, as we shall see in the following section.

## 3.2 Gaussian Kinematic Formula

The theoretical expected values of the Lipschitz-Killing curvatures (and therefore of MFs) for a Gaussian field can be computed with the Gaussian Kinematic Formula, as done in [27] for the polarised intensity of the CMB.

In our case, we want to study the function $f = Q(\phi, \theta) \cos(2\zeta) - U(\phi, \theta) \sin(2\zeta)$. This function is Gaussian, as discussed in appendix A, but it is not fully isotropic, as the behaviour in the $\theta$ and $\phi$ directions is different from the behaviour in the $\zeta$ direction. This is connected to the fact that spin $s$ functions on the sphere correspond to functions on SO(3) with a single spin component $s$, but isotropic functions of SO(3) must have components at all spins, and with the same angular power spectrum, as explained in section 2. We shall see that this anisotropy introduces some non-negligible changes in the computations of the expected values for the MFs.

Without loss of generality, we normalise the function $f$ to have unit variance; if $Q$ and $U$ have the same variance (as required by isotropy on the sphere), this is equivalent to normalising both of them to unit variance. As a consequence of spherical isotropy, the variances of all first derivatives are equal, and the covariance between orthogonal directions is zero:

$$ \mathbb{E}\left[\left(\frac{\partial Q}{\partial \theta}\right)^2\right] = \mathbb{E}\left[\left(\frac{\partial U}{\partial \theta}\right)^2\right] = \mathbb{E}\left[\left(\frac{\partial Q}{\partial \phi}\right)^2\right] = \mathbb{E}\left[\left(\frac{\partial U}{\partial \phi}\right)^2\right] = \mu, \qquad (3.5a) $$

$$ \mathbb{E}\left[\frac{\partial Q}{\partial \theta} \cdot \frac{\partial Q}{\partial \phi}\right] = \mathbb{E}\left[\frac{\partial U}{\partial \theta} \cdot \frac{\partial U}{\partial \phi}\right] = 0. \qquad (3.5b) $$

The value of $\mu$ can be directly computed as follows:

$$ \mu = \sum_{\ell} \frac{2\ell + 1}{4\pi} \frac{(\ell - s)(\ell + s + 1)}{2} C_{\ell}, \qquad (3.6) $$

where $s = 2$ and:

$$ C_{\ell} = \frac{1}{2} (C_{\ell}^{\text{EE}} + C_{\ell}^{\text{BB}}), $$

with $C_{\ell}^{\text{EE}}$ and $C_{\ell}^{\text{BB}}$ the EE and BB angular power spectra, respectively, computed from $Q$ and $U$ maps. We note that $\mu \gg 1$, as a direct consequence of normalising the field to unit variance.

The last quantity we need to introduce is the covariance matrix of the derivatives of $f$, which can be defined as $\Sigma_{ij} = \mathbb{E}[\frac{\partial f}{\partial e_i} \cdot \frac{\partial f}{\partial e_j}]$, using the orthonormal basis $\{e_1, e_2, e_3\}$,
introduced in equation (2.5). Specifically, its determinant is:

$$|\Sigma|^{\frac{1}{2}} = \frac{5}{2\sqrt{2}} \mu. \qquad (3.7)$$

We note that, in a fully isotropic case, $|\Sigma|^{\frac{1}{2}}$ would scale as $\mu^{3/2}$. Thus, the anisotropy of $f$ will introduce different $\mu$-scaling factors in the prediction of MFs. This anisotropy is associated to the different behaviour of the field $f$ in the $(\theta, \phi)$ coordinates (random Gaussian) and in the $\zeta$ coordinate (deterministic sinusoidal).

At this point, the Gaussian Kinematic Formula can be formulated for these non-isotropic Gaussian functions with a caveat: the Lipschitz-Killing curvatures must be computed not with the usual metric of the ambient manifold (see [14, 27]), but with a metric given by the covariance matrix $\Sigma$ introduced above. We denote $\mathcal{L}_i^f(A_u)$ the Lipschitz-Killing curvatures computed with this metric. For the Gaussian function $f: M \to \mathbb{R}$, the Gaussian Kinematic Formula takes the following form:

$$\mathbb{E}[\mathcal{L}_i^f(A_u)] = \sum_{k=0}^{\dim(M)-1} \begin{bmatrix} k+i \\ k \end{bmatrix} \rho_k(u) \mathcal{L}_{k+i}^f(M), \qquad (3.8)$$

where the flag coefficients are [34]:

$$\begin{bmatrix} k+j \\ k \end{bmatrix} = \frac{\omega_{k+j}}{\omega_k \omega_j} \begin{pmatrix} k+j \\ k \end{pmatrix}, \qquad \omega_j = \frac{\pi^{j/2}}{\Gamma\left(\frac{j}{2}+1\right)},$$

with $\omega_j$ representing the volume of the $j$-dimensional unit ball, while the functions $\rho_k(u)$ are defined as:

$$\rho_k(u) = \frac{1}{(2\pi)^{k/2}} \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{u^2}{2}\right) H_{k-1}(u),$$

$$H_{-1}(u) = \sqrt{2\pi} \exp\left(\frac{u^2}{2}\right) (1 - \Phi(u)),$$

$$H_0(u) = 1,$$

$$H_1(u) = u,$$

$$H_2(u) = u^2 - 1,$$

$$H_k(u) = (-1)^k \exp\left(\frac{u^2}{2}\right) \frac{d^k}{du^k} \exp\left(-\frac{u^2}{2}\right).$$

The function $\Phi$ represents the cumulative normal distribution and $H_k$ are the Hermite polynomials.

In the case of an isotropic function with the normalisation described above, it can be proven that $\mathcal{L}_j^f(M) = \mathcal{L}_j(M)\mu^{j/2}$, with $0 \le j \le \dim(M)$. This identity simplifies equation (3.8) to the Gaussian Kinematic Formula described in [27].

In the non-isotropic case, the relation between the Lipschitz-Killing curvatures computed with both metrics is not trivial, but there are two important observations. First, $\mathcal{L}_0$ corresponds to the Euler-Poincaré characteristic, which is a topological invariant. This means that this quantity depends only on the topological structure and not on the metric. Therefore, $\mathcal{L}_0^f(M) = \mathcal{L}_0(M)$. Second, there is a useful identity for $\mathcal{L}_{\dim(M)}(M)$: this quantity corresponds to the total volume of the manifold $M$ or, more technically, to its Hausdorff measure. Computing it with a different metric is analogous to a change of variable in the
integral of the measure. Therefore, $\mathcal{L}_{\text{dim}(M)}^f(M) = |\Sigma|^{1/2}\mathcal{L}_{\text{dim}(M)}(M)$, which, thus, in our case can be rephrased by: $\mathcal{L}_3^f(M) = \left(\frac{5}{2\sqrt{2}}\mu\right) \mathcal{L}_3(M)$.

Note that if we ignore the anisotropy of $f$, we would obtain that the factor is $\propto \mu^{3/2}$ instead of $\propto \mu$, yielding an incorrect behaviour with $\mu$ and therefore with the angular power spectrum of the polarisation maps. This is one of the most remarkable advantages of adopting this general formalism.

Finally, the cases of $\mathcal{L}_1^f(M)$ and $\mathcal{L}_2^f(M)$ are more complicated, as they are related to the change of integrating submanifolds with lower dimensionality when the global (non-isotropic) metric changes. Our conjecture is that the scaling relations should be $\mathcal{L}_1^f(M) \propto \mathcal{L}_1(M)$ and $\mathcal{L}_2^f(M) \propto \mathcal{L}_2(M)\mu^{1/2}$. In section 6.2 we will see that these guesses are highly compatible with simulations, and we will compute the proportionality factors.

## 3.3 Theoretical predictions

We are now in the position to explicitly compute the theoretical predictions for the MFs of the excursion set. In our case, the global manifold is $M = \text{SO}(3)$. We note that the right hand side of the Gaussian Kinematic Formula, equation (3.8), always contains the term proportional to $\mathcal{L}_3(\text{SO}(3))$ as the leading term. As mentioned above, the last Lipschitz-Killing curvature is equal to the Hausdorff measure (or volume) of the manifold, which in this case is $\text{Vol}(\text{SO}(3)) = 4\pi^2$, as shown in equation (2.4).

Only $\mathcal{L}_3$ is needed to compute all MFs at leading order but we report the rest of the Lipschitz-Killing curvatures of $\text{SO}(3)$ for completeness:

$$ \mathcal{L}_3(\text{SO}(3)) = 4\pi^2, \qquad (3.9a) $$

$$ \mathcal{L}_2(\text{SO}(3)) = 0, \qquad (3.9b) $$

$$ \mathcal{L}_1(\text{SO}(3)) = 6\pi, \qquad (3.9c) $$

$$ \mathcal{L}_0(\text{SO}(3)) = 0. \qquad (3.9d) $$

For the topology-inclined readers, the last line implies that the Euler-Poincaré characteristic vanishes: $\chi(\text{SO}(3)) = 0$; this is a consequence of the 3-sphere being a (2-fold) cover of $\text{SO}(3)$ and $\chi(\mathbb{S}^3) = 0$.

We can now compute the Lipschitz-Killing curvatures of the excursion set of $f$ using the Gaussian Kinematic Formula, equation (3.8). We assume no mask in these computations, but they can be readily introduced by modifying the Lipschitz-Killing curvatures of the global manifold, most notably by multiplying the total volume by the sky fraction.

**Volume of the excursion set, $V_0$.** We use the Gaussian Kinematic Formula with $j = 3$ to compute the expected value of the volume of the excursion set, along with the expression $\mathcal{L}_3^f(M) = \left(\frac{5}{2\sqrt{2}}\mu\right) \mathcal{L}_3(M)$ derived above:

$$ \begin{aligned} \mathbb{E}[\mathcal{L}_3(A_u)] &= \frac{2\sqrt{2}}{5\mu}\mathbb{E}[\mathcal{L}_3^f(A_u)] \\ &= \frac{2\sqrt{2}}{5\mu}\sum_{k=0}^{3-3} \begin{bmatrix} k+3 \\ k \end{bmatrix} \rho_k(u) \mathcal{L}_{k+3}(\text{SO}(3)) \qquad (3.10) \\ &= \frac{2\sqrt{2}}{5\mu} \begin{bmatrix} 3 \\ 0 \end{bmatrix} \rho_0(u) \left[ \frac{5}{2\sqrt{2}} \mu \mathcal{L}_3(\text{SO}(3)) \right] = 4\pi^2 [1 - \Phi(u)], \end{aligned} $$
where $\Phi(u)$ is the cumulative distribution function of the standard normal distribution. It can be seen that $V_0$ is equal to the total volume at low thresholds (since the field at all points is greater than $u$), and it is 0 at high thresholds (since the field in all points is lower than $u$). Interestingly, the definition of $f$ implies that at $u = 0$, $\mathcal{L}_3(A_0) = 2\pi^2$ with no scattering whatsoever. We note that the factor $\mu$ cancels out, and the theoretical expectation depends only on the threshold $u$, not on the field itself (as long as it is normalised to have unit variance). The summation in equation (3.10) has exactly one term, corresponding to the volume of SO(3), so no approximation is needed for this prediction.

**Area of the boundary of the excursion set, $V_1$.** We can compute the theoretical expectation by setting $j = 2$ in the Gaussian Kinematic Formula and assuming the scaling relation explained above:

$$ \mathcal{L}_2^f(M) = K_1^{-1} \mu^{1/2} \mathcal{L}_2(M), \qquad (3.11) $$

where we have introduced the unknown constant $K_1$ (as it is related to $V_1$). Thus:

$$ \begin{aligned} \mathbb{E}[\mathcal{L}_2(A_u)] &= K_1 \mu^{-1/2} \mathbb{E}[\mathcal{L}_2^f(A_u)] \\ &= K_1 \mu^{-1/2} \sum_{k=0}^{3-2} \begin{bmatrix} k+2 \\ k \end{bmatrix} \rho_k(u) \mathcal{L}_{k+2}^f(\mathrm{SO}(3)) \\ &= K_1 \mu^{-1/2} \begin{bmatrix} 3 \\ 1 \end{bmatrix} \rho_1(u) \begin{bmatrix} \frac{5}{2\sqrt{2}} \mu \mathcal{L}_3(\mathrm{SO}(3)) \end{bmatrix} \\ &= K_1 \frac{5}{\sqrt{2}} \sqrt{\mu} 4\pi^2 \left[ \frac{1}{2\pi} \exp\left(\frac{-u^2}{2}\right) \right], \end{aligned} \qquad (3.12) $$

where the square bracket in the fourth line corresponds to $\rho_1(u)$ and the term $k = 0$ in the summation is zero because of equation (3.9).

**Mean curvature of the boundary of the excursion set, $V_2$.** In a similar way, we can compute the theoretical expectation of $V_2$ by setting $j = 1$ in the Gaussian Kinematic Formula and assuming the scaling relation explained above:

$$ \mathcal{L}_1^f(M) = K_2^{-1} \mathcal{L}_1(M), \qquad (3.13) $$

where we again introduce an unknown constant $K_2$. Thus:

$$ \begin{aligned} \mathbb{E}[\mathcal{L}_1(A_u)] &= K_2 \mathbb{E}[\mathcal{L}_1^f(A_u)] \\ &= K_2 \sum_{k=0}^{3-1} \begin{bmatrix} k+1 \\ k \end{bmatrix} \rho_k(u) \mathcal{L}_{k+1}^f(\mathrm{SO}(3)) \\ &= K_2 \begin{bmatrix} 3 \\ 2 \end{bmatrix} \rho_2(u) \begin{bmatrix} \frac{5}{2\sqrt{2}} \mu \mathcal{L}_3(\mathrm{SO}(3)) \end{bmatrix} + \mathcal{O}(\mathcal{L}_1^f(\mathrm{SO}(3))) \\ &= K_2 \frac{5}{\sqrt{2}} \mu 4\pi^2 \left[ \frac{u}{(2\pi)^{3/2}} \exp\left(\frac{-u^2}{2}\right) \right] + \mathcal{O}(\mu^0), \end{aligned} \qquad (3.14) $$

where the square bracket in the fourth line corresponds to $\rho_2(u)$. Only two terms in the summation are not zero, one of order $\mu$ and another of order 1. Given that $\mu \gg 1$ (in the case of cosmological fields, we typically have $\mu \sim 10^5$), we neglect the second order terms in the computations.
**Euler-Poincaré characteristic of the excursion set, V3.** Finally, we compute the theoretical prediction for the Euler-Poincaré characteristic of the excursion sets of $f$. As explained before, $\mathcal{L}_0^f(M) = \mathcal{L}_0(M)$, since it is a topological invariant. Thus:

$$ 
\begin{align}
\mathbb{E}[\mathcal{L}_0(A_u)] &= \mathbb{E}[\mathcal{L}_0^f(A_u)] \nonumber \\
&= \sum_{k=0}^{3-0} \begin{bmatrix} k+0 \\ k \end{bmatrix} \rho_k(u) \mathcal{L}_{k+0}^f(\mathrm{SO}(3)) \nonumber \\
&= \begin{bmatrix} 3 \\ 3 \end{bmatrix} \rho_3(u) \left[ \frac{5}{2\sqrt{2}} \mu \mathcal{L}_3(\mathrm{SO}(3)) \right] + \begin{bmatrix} 1 \\ 1 \end{bmatrix} \rho_1(u) \left[ \mathcal{L}_1^f(\mathrm{SO}(3)) \right] \tag{3.15} \\
&= \frac{5}{2\sqrt{2}} \mu 4\pi^2 \left[ \frac{(u^2-1)}{(2\pi)^2} \exp\left(-\frac{u^2}{2}\right) \right] + \mathcal{O}(\mu^0), \nonumber
\end{align}
$$

where the square bracket in the fourth line corresponds to $\rho_3(u)$. Again, only two terms in the summation are not zero, one of order $\mu$ and another of order 1. The latter, the term corresponding to $\mathcal{L}_1^f(\mathrm{SO}(3))$, can be safely ignored.

We note that neglecting the anisotropy of $f$ would yield an incorrect exponent for $\mu$ in the theoretical prediction for $\mathcal{L}_0$ (i.e., the Euler-Poincaré characteristic).

**Predictions for the MFs.** We can now convert the predictions for Lipschitz-Killing curvatures into MFs. Additionally, in order to ease the comparison with data and the interpretation, we work with normalised MFs, where the quantities are divided over the volume (in our case this is $4\pi^2$, possibly multiplied by the sky fraction $f_{\text{sky}}$ if we impose a mask):

$$ v_i = \frac{V_i}{4\pi^2}. \tag{3.16} $$

The theoretical predictions of the normalised MFs are as follows:

$$ \mathbb{E}[v_0] = 1 - \Phi(u), \tag{3.17a} $$

$$ \mathbb{E}[v_1] = K_1 \frac{5}{6\pi\sqrt{2}} \mu^{1/2} \exp\left(\frac{-u^2}{2}\right), \tag{3.17b} $$

$$ \mathbb{E}[v_2] = K_2 \frac{5}{24\pi^{3/2}} \mu u \exp\left(\frac{-u^2}{2}\right), \tag{3.17c} $$

$$ \mathbb{E}[v_3] = \frac{5}{8\sqrt{2}\pi^2} \mu(u^2 - 1) \exp\left(\frac{-u^2}{2}\right). \tag{3.17d} $$

We recall that the first two formulae are exact, while the last two are correct at the leading order in $\mu$. One of the main advantages of this formalism is that the expected values of the normalised MFs do not depend on the use of masks. These theoretical predictions show the power of this approach. We are able to predict a topological feature of the field knowing only the parameter $\mu$, which can be computed from the angular power spectra of the E and B modes. Alternatively, this parameter can also be obtained directly on the maps, which may be recommended when they are masked to avoid leakage effects, which distort the estimation of the polarisation angular power spectra.
Similar results to equation (3.17) are known in the literature for both Gaussian and weakly non-Gaussian isotropic fields defined on $\mathbb{R}^3$ [30, 31]. The main differences in our approach are the following:

* The use of the Gaussian Kinematic Formula, which unifies the theoretical computations for all MFs.

* We consider a space with a non-trivial metric given in equation (2.2); this approach can be easily generalised to arbitrary manifolds. We note that this only affects the normalisation factors in equation (3.17).

* The studied field $f$ given in equation (2.1) is not isotropic, leading to additional complications. We showed how to deal with them when the anisotropic behaviour is known. In our case, this implied a distinct trend with $\mu$ for the Euler-Poincaré characteristic.

In this work, we present only the predictions for Gaussian fields. The extension to weakly non-Gaussian fields could be obtained by adopting an approach similar to that presented in [30], by considering higher-order cumulants of the probability density function.

# 4 Implementation: Pynkowski

We implement both the theoretical predictions (presented in section 3) and the computation of the MFs of $f$ on data (to be explained in this section). Both aspects are to be included in the publicly available Python package called Pynkowski, which can be found in https://github.com/javicarron/pynkowski and was first introduced in [27] in the context of MFs for scalar maps, such as $T$ and $P^2 = Q^2 + U^2$.

We develop a software to compute all the necessary quantities related to $f$, including its derivatives in all directions, the covariant gradient, and the Hessian. We use this code on Gaussian CMB simulations in order to compare the values of their MFs with the theoretical predictions. After this validation, it can be applied to arbitrary maps to assess any deviation from statistical isotropy and Gaussianity. The actual implementation of the computation of these MFs on maps is described in this section.

For a map $f$, we store the values of $Q$ and $U$ in the HEALPix pixelisation scheme [38]. The variable $\zeta$ is not pixelised, it is calculated exactly for each pixel $\xi_i$ as $f(\xi_i, \zeta) = Q(\xi_i) \cos(2\zeta) - U(\xi_i) \sin(2\zeta)$.

The first and second spatial derivatives of $Q$ and $U$ with respect to $\theta$ and $\phi$ are also obtained when needed, with the help of the healpy<fnref n="2" /> function alm2map_der, which performs such computation in harmonic space. The first and second spatial derivatives of $f$ with respect to its three variables are simplified analytically and expressed as exact functions of the polarisation maps and their derivatives. This means that the only pixelisation is due to the original one of $Q$ and $U$ and their derivatives, while the treatment in $\zeta$ is always exact in order to avoid numerical artefacts.

The first step in our procedure after loading the map is to normalise $f$ to have unit variance. We do not assume that $Q$ and $U$ have similar statistical properties, but if they do, this step is equivalent to normalising both of them to have unit variance. With such a normalisation of $f$, we can express the MFs as a function of an adimensional threshold $u$.

<footnote n="2">²https://github.com/healpy/healpy.</footnote>
For convenience, we compute the MFs divided by the total volume of SO(3) to have the normalised MFs ($v_i = \frac{V_i}{4\pi^2}$).

Before explaining the specific computations for each MF, we shall introduce some common notation. Let $P(\xi_i) = \sqrt{Q(\xi_i)^2 + U(\xi_i)^2}$, i.e., the maximum value of $f$ in each pixel $\xi_i$. For a fixed $\xi_i$, we note that $f$ has a simple sinusoidal behaviour in its variable $\zeta$: it is always below $u$ if $u > P(\xi_i)$, always above $u$ if $u < -P(\xi_i)$, and above $u$ in a single segment of $\zeta$ otherwise. We use this observation to split all the integrals involved in the computation of the MFs. In the case where $f(\xi_i) > u$ only in a range of $\zeta$, let $\zeta_M(\xi_i)$ be the angle at which $f(\xi_i, \zeta)$ is maximum, while $\zeta_1(\xi_i)$ and $\zeta_2(\xi_i)$ the angles for which $f(\xi_i) = u$; they can be calculated as:

$$ \zeta_M(\xi_i) = \frac{1}{2} \arctan \left( -\frac{U(\xi_i)}{Q(\xi_i)} \right), \qquad (4.1a) $$

$$ \zeta_1(\xi_i) = \zeta_M(\xi_i) - \frac{1}{2} \arccos \left( \frac{u}{P(\xi_i)} \right), \qquad (4.1b) $$

$$ \zeta_2(\xi_i) = \zeta_M(\xi_i) + \frac{1}{2} \arccos \left( \frac{u}{P(\xi_i)} \right), \qquad (4.1c) $$

where the arctan in the first line is defined to be in the quadrant corresponding to the coordinate $(Q(\xi_i), U(\xi_i))$. Polarisation directions $\zeta$ and their differences are always defined between $0$ and $\pi$, due to the geometry of SO(3).

**First MF, $v_0$.** It can be seen as the volume fraction of the manifold for which $f > u$:

$$ v_0(u) = \frac{1}{4\pi^2} \int_{SO(3)} \Theta(f(x) - u) \, dx, \qquad (4.2) $$

where $\Theta(r)$ is the Heaviside function ($1$ where $r \ge 0$; $0$ otherwise). We can split the integrand by pixel according to the cases explained above and integrate with respect to $\zeta$ first. The integrand will be $0$ for pixels where $f$ is always smaller than the threshold ($u > P(\xi_i)$) and it will be $1$ if $u < -P(\xi_i)$. For all the other pixels, it will be $1$ as long as $\zeta$ is between $\zeta_1$ and $\zeta_2$ and $0$ otherwise, so the integral in $\zeta$ will yield the length of the interval between these two angles, i.e., $\arccos\left(\frac{u}{P(\xi_i)}\right) \frac{1}{\pi}$.

Therefore, to compute equation (4.2), we just need the length of this interval for every pixel and average it over all pixels.

**Second MF, $v_1$.** It can be computed as the area of the boundary of the excursion sets, i.e., the manifold defined by $f = u$:

$$ \begin{aligned} v_1(u) &= \frac{1}{4\pi^2} \frac{1}{6} \int_{\partial A_u} da \\ &= \frac{1}{4\pi^2} \frac{1}{6} \int_{SO(3)} \delta(f(x) - u) \cdot |\nabla f(x)| \, dx, \end{aligned} \qquad (4.3) $$

where the second equality comes from a change of coordinates from the surface element of the boundary $da$ to a volume element on SO(3), $dx$. The gradient of $f$ is denoted by $\nabla f$ and can be computed with equation (2.5); $\delta$ is the Dirac delta.

In a fully pixelised field (such as $T$ or $P$), one typically has to bin the threshold and approximate the delta function (see [27] for a discussion in that case). However, in the
present framework we can exactly determine the points where $f(\xi_i, \zeta) = u$ for each pixel (if they exist) as $(\xi_i, \zeta_1)$ and $(\xi_i, \zeta_2)$.

To compute the integral, it suffices to evaluate $|\nabla f(x)|$ at these points (adding the values corresponding to both $\zeta_1$ and $\zeta_2$ for each $\xi_i$), fill with 0 the rest of the pixels, and compute the average value of this quantity across the entire map. This result is then divided by 6, the normalisation factor.

**Third MF, v2.** It is the average mean curvature ($H$) of the boundary of the excursion sets:

$$ 
\begin{aligned}
v_2(u) &= \frac{1}{4\pi^2} \frac{1}{6\pi} \int_{\partial A_u} H(a) da \\
&= \frac{1}{4\pi^2} \frac{1}{6\pi} \int_{\text{SO}(3)} \delta(f(x) - u) \cdot H(x) \cdot |\nabla f(x)| dx,
\end{aligned} 
\quad (4.4) $$

where we use the same change of coordinates as in the previous case, which introduces the $\delta$ and the factor $|\nabla f(x)|$.

In general, evaluating the mean curvature of a surface is far from trivial. In this case, however, the surface is defined implicitly by the function $F(x) := f(x) - u = 0$, so we can use the expression for the mean curvature in an implicitly defined surface [39]:

$$ H = \frac{\nabla F \cdot \mathcal{H}(F) \nabla F^T - |\nabla F|^2 \operatorname{Tr}(\mathcal{H}(F))}{2|\nabla F|^3}, \quad (4.5) $$

where all quantities are functions of $x$.

As before, we compute the integral by evaluating $H(x) \cdot |\nabla f(x)|$ in the points where $f = u$ (adding for each pixel the values corresponding to both $\zeta_1$ and $\zeta_2$), we fill with 0 the rest of the pixels, and calculate the average value of this quantity across the entire map. This result is then divided by $6\pi$, the normalisation factor.

**Fourth MF, v3.** It can be computed as the average Gaussian curvature ($K$) of the boundary of the excursion sets:

$$ 
\begin{aligned}
v_3(u) &= \frac{1}{4\pi^2} \frac{1}{4\pi} \int_{\partial A_u} K(a) da \\
&= \frac{1}{4\pi^2} \frac{1}{4\pi} \int_{\text{SO}(3)} \delta(f(x) - u) \cdot K(x) \cdot |\nabla f(x)| dx,
\end{aligned} 
\quad (4.6) $$

where we use the same strategy as in the previous case, now computing the Gaussian curvature $K(x)$ instead of the mean curvature $H(x)$. In the case of an implicitly defined surface (again $F(x) := f(x) - u = 0$), the Gaussian curvature can be computed as [39]:

$$ K = - \frac{\begin{vmatrix} H(F) & \nabla F^T \\ \nabla F & 0 \end{vmatrix}}{|\nabla F|^4} = - \frac{\begin{vmatrix} F_{xx} & F_{xy} & F_{xz} & F_x \\ F_{xy} & F_{yy} & F_{yz} & F_y \\ F_{xz} & F_{yz} & F_{zz} & F_z \\ F_x & F_y & F_z & 0 \end{vmatrix}}{|\nabla F|^4}. \quad (4.7) $$

We then compute the integral by evaluating $K(x) \cdot |\nabla f(x)|$ at the points where $f = u$ (adding for each pixel the values corresponding to both $\zeta_1$ and $\zeta_2$), we fill with 0 the rest of the pixels, and calculate the average value of this quantity across the entire map. This result is then divided by $4\pi$, the normalisation factor.
## 4.1 General considerations

All the above computations are performed without pixelising the variable $\zeta$, only $\theta$ and $\phi$. However, numerical issues can arise, especially in $v_3$ and, to a lesser extent, $v_2$, due to the computation of the second spatial derivatives in pixelised maps. In order to avoid artefacts, we must work with smooth maps, where the value of the field does not change abruptly between neighbouring pixels; in other words, the derivatives of the map must exist and be reasonably smooth. This can be seen as a requirement on the maps $Q$ and $U$ to have negligible power at high multipoles or on the needed pixel resolution (the $N_{side}$ parameter in HEALPix maps). This is not only a computational limitation but, in fact, a mathematical requirement. Indeed, the Gaussian Kinematic Formula holds if the boundaries of the excursion sets are twice differentiable (see [34]). In practice, this means that the angular power spectrum must decay fast enough. For a study where this condition does not hold, see, e.g., [40].

Nevertheless, we note that the theoretical predictions for $v_2$ and $v_3$ are correct up to leading order in $\mu$. Given that $\mu$ increases with the considered range of multipoles, these approximations are increasingly accurate at higher multipoles. Therefore, it may be convenient to apply this formalism on needlet components. Indeed, needlets naturally filter out the low and high multipoles, thus reducing both the pixelisation effects and the error in the approximation in $\mu$. MFs in needlet domain were studied for scalar fields on the sphere in, e.g., [15, 35].

# 5 Simulations

In this section we introduce the different polarisation fields that we analyse with MFs.

## 5.1 Monochromatic maps

The exact value of $v_1$ and $v_2$, introduced in section 3.3, depends on the unknown normalisation factors $K_1$ and $K_2$, respectively.

We generate simulated maps, using the `healpy` Python package [41], in order to verify the predicted dependence with $\mu$ of the four MFs, and compute the two unknown normalisation factors. We create 1600 monochromatic maps, i.e., maps with power in a single multipole $\ell_m$:

$$C_{\ell}^{\mathrm{EE}} = C_{\ell}^{\mathrm{BB}} = \begin{cases} 1 & \ell = \ell_m, \\ 0 & \text{otherwise.} \end{cases} \qquad (5.1)$$

We simulate 8 sets of 200 maps, each set with a different $\ell_m$. We have considered $\ell_m = \{100, 485, 675, 825, 950, 1060, 1165, 1255\}$. These values of $\ell_m$ are chosen in order to have a set of approximately equally-spaced values of $\mu_m = (\ell_m - s)(\ell_m + s + 1)/2$, where $s = 2$.

## 5.2 Realistic angular power spectrum

In order to validate the theoretical predictions and the implemented computation of MFs on realistic maps, we generate 300 Gaussian isotropic CMB simulations with the best fit angular power spectrum for CMB polarisation (E and B modes) reported by [42]. We use a map resolution of $N_{side} = 1024$. We smooth the maps with a Gaussian beam of FWHM = 15' to avoid pixelisation effects in the computation of the spatial derivatives.

In this way, we are able to simultaneously achieve three goals: assess the accuracy of the theoretical formulae obtained in section 3.3, verify the implementation of MFs computation on simulated maps, and check that the normalisation constants extend correctly to non-monochromatic maps.
# 6 Results

## 6.1 Normalisation constants

In this section we determine the normalisation constants for $v_1$ and $v_2$, i.e., $K_1$ and $K_2$. As mentioned in section 3.3, the exact values for the predicted MFs can only be computed for $v_0$ and $v_3$, while the prediction for $v_1$ and $v_2$ is correct except for a global normalisation factor related to the anisotropic $\zeta$ direction, see equation (3.17). We note that both $K_1$ and $K_2$ are truly constants and therefore do not depend on the threshold, mask, or value of $\mu$.

We estimate the values of these constants by computing the MFs on the monochromatic maps introduced in section 5.1 with the software presented in section 4 and soon to be included in Pynkowski. Then, $K_1$ and $K_2$ are obtained as the ratio between the computation on the maps and theoretical predictions, averaging over all simulations and all thresholds. We find that the constants are (average and standard deviation):

$$K_1 = 0.31912 \pm 0.00001, \qquad (6.1a)$$

$$K_2 = 0.7088 \pm 0.0002. \qquad (6.1b)$$

We have verified that these values do not significantly vary with the threshold or value of $\mu$ (i.e., multipole). We have also checked that the constants are valid for non-monochromatic maps, as we will see in section 6.3.

Lastly, we note that the theoretical predictions yield the correct normalisation factor for $v_0$ and $v_3$, as expected. If we repeat the procedure with additional multiplicative factors, $K_{v_0}$ and $K_{v_3}$, for these MFs, we see that both are compatible with 1:

$$K_{v_0} = 1.0000 \pm 0.0001, \qquad (6.2a)$$

$$K_{v_3} = 1.001 \pm 0.003. \qquad (6.2b)$$

## 6.2 Scaling relations

The theoretical expectations of the MFs in section 3 assume a particular dependence of $\mu$ for the Lipschitz-Killing computed with different metrics (see the last paragraph in section 3.2). In this section, we show that this scaling is strongly supported by simulations.

The assumed $\mu$ dependence has a direct impact on the predictions for $v_1$ and $v_2$. If the assumption was wrong, these MFs would present a different scaling with the parameter $\mu$. We test for deviations by using the MFs computed on monochromatic maps with a wide range of values of $\mu$.

In order to assess the dependence of $v_1$ and $v_2$ on $\mu$, we consider $\bar{v}_i = \langle \frac{v_i}{v_i(\ell_m=950)} \rangle$, where $\langle \cdot \rangle$ represents the average over the thresholds between $u = -3$ and $u = 3$. By considering this quantity, we remove the dependence on threshold and constant factors, leaving only the dependence on $\mu$.

In figure 1 we show the trend of $\bar{v}_1$ and $\bar{v}_2$ for different values of $\mu$. We observe that they are perfectly compatible with the scaling assumed in the theory: $\mu^{1/2}$ and $\mu$, respectively.

Additionally, we have verified with the same procedure that $v_0$ and $v_3$ scale as $\mu^0$ and $\mu$, respectively, as predicted by the theory. All scaling relations remain unchanged by considering different thresholds and non-monochromatic maps, as expected.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Dependence on $\mu$ of the normalised MFs $\bar{v}_1$ (blue dots) and $\bar{v}_2$ (red dots), defined in the text. We can see in dashed lines the curves for $\mu^{1/2}$ (blue) and $\mu$ (red), corresponding to the expected scaling relations.</figcaption></figure>

## 6.3 Gaussian CMB simulations

In this section we verify that the theoretical predictions for the MFs (see section 3.3) agree with the results computed on CMB Gaussian isotropic maps with a realistic angular power spectrum. We use the 300 simulations introduced in section 5.2, which are generated with the Planck best-fit polarisation angular power spectra; we use the corresponding value of $\mu$ for the theoretical predictions. We compute the MFs (see section 4) at thresholds between $u = -4$ and $u = 4$, with a spacing of $\Delta u = 0.2$. We use the values of the normalisation factors $K_1$ and $K_2$ obtained in monochromatic maps, i.e., equation (6.1). This means that there are no free parameters in the analysis of this section.

Figure 2 presents the comparison between the MFs computed on these maps and the corresponding theoretical predictions. We find full compatibility between theory and simulations for all MFs. In all cases, the residuals are well within the $1\sigma$ region, where $\sigma$ is the dispersion of the MFs among the 300 simulations.

Additionally, to explore any possible systematic deviation, we compare the average residual between simulations and theoretical expectations. This is done considering the standard deviation of the mean ($\sigma_m = \frac{\sigma_{\text{sims}}}{\sqrt{300}}$). This comparison is shown in the bottom sections of each panel in figure 2, in terms of $\Delta v_i/\sigma_m$, i.e., deviations over $\sigma_m$. It can be seen that no point deviates more than $2\sigma_m$, and there are no significant systematic residuals.

The four MFs are all perfectly compatible with the theory, in the case of both individual simulations and of average behaviour. We verify that the results hold for arbitrary angular power spectra, as long as there is no significant power in the scales corresponding to the pixel size.

We also note the low statistical variation of these curves: the relative uncertainty of every point is below 1 part in 1000 for $v_0$, below 1% for $v_1$ and $v_2$ and below 10% for $v_3$. These values of the statistical standard deviation depend on the angular power spectrum of the studied map: we observe that the variance decreases when increasing the maximum multipole considered, $\ell_{\text{max}}$ (i.e., when increasing the resolution of the maps). This is also known to be the case for temperature maps, as quantitatively studied by [35].
Gaussianities or deviations from statistical isotropy. The detection of such effects would reveal physical information about the Early Universe.

In the past, MFs have been applied mainly to CMB temperature anisotropies and, more generally, to scalar fields in order to search for such signatures in Cosmological data [14, 15, 27].

In this work, we have extended the MFs formalism to the polarisation field by lifting it to the higher-dimensional manifold SO(3), thus considering the full information embedded in this observable. In this framework, the complex part of the field is just a translation of the real one. Therefore, without loss of generality, we take only the real part, $f(\phi, \theta, \zeta) : SO(3) \to \mathbb{R}$ as:

$$f(\phi, \theta, \zeta) = Q(\phi, \theta) \cos(2\zeta) - U(\phi, \theta) \sin(2\zeta), \quad (7.1)$$

which is a scalar field that can be interpreted as the linear polarisation observed at the point $(\phi, \theta)$ on the sky in the polarisation direction $\zeta$.

To summarise:

* We introduce a formalism to predict the values of MFs for $f$ in the case of Gaussian Stokes parameters through the use of the Gaussian Kinematic Formula for anisotropic fields (see sections 3.2 and 3.3).

* We implement a code to compute MFs of $f$ from input Q and U maps pixelised with the HEALPix convention and to be released in the already publicly available Python package Pynkowski.

* We find that theoretical predictions are compatible with the computations on Gaussian CMB maps generated with the Planck best-fit angular power spectra.

The application on CMB polarisation data of statistical tools beyond the angular power spectrum (such as the MFs formalism introduced in this paper) will grow in importance with upcoming much more sensitive experiments from both ground and space, e.g., ACT [44], SPT [45], Simons Observatory [46], LiteBIRD [47], PICO [48].

We are currently applying this formalism to the CMB polarisation data from Planck and to simulated future experiments. The corresponding results will be presented in an upcoming paper.

Furthermore, this kind of analysis can help to blindly detect the presence of non-Gaussian contaminants in the maps like residual foregrounds contamination. It can also be used to characterise the morphology of Galactic emission, as recently done in [20, 21] for synchrotron radiation.

The application proposed in this paper is versatile and can be performed on other spin-2 quantities. Moreover, in this work we have developed the predictions for Gaussian fields, but this framework can be extended to the weakly non-Gaussian case by following the approach adopted in [30] for fields defined on Euclidean spaces.

## Acknowledgments

The authors thank Michele Stecconi and Maurizia Rossi for insightful discussions. JCD and MM acknowledge support from ASI/INFN grant n. 2021-43-HH.0. MM and NV acknowledge support by ASI/COSMOS grant n. 2016-24-H.0 and ASI/LiteBIRD grant n. 2020-9-HH.0. Part of this work was also supported by the InDark INFN project. DM acknowledges support
from the MIUR Excellence Project awarded to the Department of Mathematics, Università di Roma Tor Vergata, CUP E83C18000100006. DM is also grateful to the Department of Excellence Programme MatModTov for support. This paper is supported by the Fondazione ICSC, Spoke 3 Astrophysics and Cosmos Observations, National Recovery and Resilience Plan (Pi-ano Nazionale di Ripresa e Resilienza, PNRR) Project ID CN_00000013 “Italian Research Center on High-Performance Computing, Big Data and Quantum Computing” funded by “MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di campioni nazionali di R&S (M4C2-19) — Next Generation EU (NGEU)”.

## A Polarization as a random field on SO(3)

In section 2, we argued that the CMB polarisation spin field can be lifted to a complex function of SO(3) and that its real and the imaginary parts are just translations of each other. Therefore, it suffices to consider only the real part, which corresponds to the scalar function $f$ introduced in equation (2.1). In this appendix we will prove this statement. We note that this argument generalises to spherical fields of arbitrary spin.

We start by recalling the usual spherical harmonic decomposition of a spin-$s$ field:

$$Q(\theta, \phi) + iU(\theta, \phi) = \sum_{\ell m} a_{\ell m; s} Y_{\ell m}(\theta, \phi),$$

with $sY_{\ell m}(\theta, \phi)$ the spin-weighted spherical harmonics. We can rewrite this decomposition by using the Wigner $D$-matrix $D_{\ell,m;s}$:

$$Q(\theta, \phi) + iU(\theta, \phi) = \sum_{\ell m} a_{\ell m; s} D_{\ell, m; s}(\phi, \theta, 0).$$

The generalisation of this expression allows for the introduction of a three-dimensional complex field $\mathcal{P}$:

$$\mathcal{P}(\phi, \theta, \zeta) = Q(\phi, \theta, \zeta) + iU(\phi, \theta, \zeta) = \sum_{\ell m} a_{\ell m; s} D_{\ell, m; s}(\phi, \theta, \zeta),$$

where $Q$ and $U$ are the Stokes parameters $Q$ and $U$ when the local reference frame is rotated by an angle $\zeta$:

$$Q(\phi, \theta, \zeta) = Q(\theta, \phi) \cos(s\zeta) - U(\theta, \phi) \sin(s\zeta),$$

$$U(\phi, \theta, \zeta) = U(\theta, \phi) \cos(s\zeta) + Q(\theta, \phi) \sin(s\zeta).$$

$\mathcal{P}$ is a Gaussian complex field since it is given by a linear combination of Gaussian complex-valued coefficients $a_{\ell m; s}$ times deterministic functions. Therefore its real and imaginary parts, $Q$ and $U$, are also Gaussian fields. Such deterministic functions are the Wigner $D$-matrices, whose dependence on their variables can be explicitly expressed by using the Wigner's $d$ function:

$$D_{\ell, m; s}(\phi, \theta, \zeta) = \exp(-im\phi) d_{\ell, m; s}(\theta) \exp(is\zeta).$$
We can see now that the following symmetry holds for the function $\mathcal{P}$:

$$ 
\begin{aligned}
\mathcal{P}\left(\phi, \theta, \zeta + \frac{\pi}{2s}\right) &= \sum_{\ell m} a_{\ell m; s} D_{\ell, m; s}(\phi, \theta, \zeta) \\
&= \sum_{\ell m} a_{\ell m; s} \exp(-im\phi) d_{\ell, m; s}(\theta) \exp\left(is\left(\zeta + \frac{\pi}{2s}\right)\right) \\
&= i \sum_{\ell m} a_{\ell m; s} \exp(-im\phi) d_{\ell, m; s}(\theta) \exp(is\zeta) \\
&= i \mathcal{P}(\phi, \theta, \zeta).
\end{aligned}
$$

Likewise:

$$ \mathcal{P}\left(\phi, \theta, \zeta - \frac{\pi}{2s}\right) = -i \mathcal{P}(\phi, \theta, \zeta). $$

In words, this means that a rotation of $\frac{\pi}{2s}$ of the local reference frame is equal to multiplying the function by a factor $i$ or, equivalently, swapping the real and imaginary components of $\mathcal{P}: \mathcal{Q}$ and $\mathcal{U}$. These components are therefore linked by the following expression:

$$ \mathcal{Q}\left(\phi, \theta, \zeta + \frac{\pi}{2s}\right) = \mathrm{Re}\left[\mathcal{P}\left(\phi, \theta, \zeta + \frac{\pi}{2s}\right)\right] = -\mathrm{Im}\left[\mathcal{P}(\phi, \theta, \zeta)\right] = -\mathcal{U}(\phi, \theta, \zeta). $$

Consequently, both the real and imaginary part of the lifted polarisation field $\mathcal{P}$ convey the same information, so it suffices to study one of them. We choose to study the real part $\mathcal{Q}$, which we call $f$ in the main text. As noted before, this function is a real-valued Gaussian field on SO(3).

## B Metric and derivatives in SO(3)

In this paper, the introduced MFs formalism is based on the definition of the field $f(\phi, \theta, \zeta) = \mathcal{Q}(\phi, \theta) \cos(2\zeta) - \mathcal{U}(\phi, \theta) \sin(2\zeta)$, with domain in SO(3). The geometry of this manifold is not trivial. In particular, we were unable to find in the literature the expressions for an orthonormal basis and for the Hessian of a function defined on this manifold. We performed these computations and we report them in this appendix in case they are useful to the community.

SO(3) can be seen as the group of rotations of the sphere, which can be parameterised by the Euler angles. These accept many different definitions, depending on the order of the axis selected to perform the rotations. To make the angles compatible with the usual coordinates on the sphere, we must select the definition in the zyx convention. Then, a generic rotation of the sphere can be expressed as [see, e.g., 49]:

$$ R(\phi, \theta, \zeta) = R_x(\zeta) R_y(\theta) R_z(\phi), $$

where:

$$ R_x(\zeta) = \begin{pmatrix} 1 & 0 & 0 \\ 0 \cos(\zeta) & -\sin(\zeta) \\ 0 \sin(\zeta) & \cos(\zeta) \end{pmatrix}, $$

$$ R_y(\theta) = \begin{pmatrix} \cos(\theta) & 0 & \sin(\theta) \\ 0 & 1 & 0 \\ -\sin(\theta) & 0 & \cos(\theta) \end{pmatrix}, $$

$$ R_z(\phi) = \begin{pmatrix} \cos(\phi) & -\sin(\phi) & 0 \\ \sin(\phi) & \cos(\phi) & 0 \\ 0 & 0 & 1 \end{pmatrix}. $$
It is then possible to compute the derivative of a random rotation with respect to the three parameters. Now, taking the Frobenius inner product between matrices (the sum of the element-wise multiplication of matrices), the metric of SO(3) reads:

$$g_{ij} = \langle \partial_i R, \partial_j R \rangle = 2 \begin{pmatrix} 1 & 0 & \sin(\theta) \\ 0 & 1 & 0 \\ \sin(\theta) & 0 & 1 \end{pmatrix},$$

in the order $(\phi, \theta, \zeta)$. We can now determine the Christoffel symbols in this manifold:

$$\Gamma_{\theta\zeta}^{\phi} = \Gamma_{\zeta\theta}^{\phi} = \Gamma_{\phi\theta}^{\zeta} = \Gamma_{\theta\phi}^{\zeta} = \frac{1}{2\cos(\theta)},$$

$$\Gamma_{\theta\phi}^{\phi} = \Gamma_{\phi\theta}^{\phi} = \Gamma_{\zeta\theta}^{\zeta} = \Gamma_{\theta\zeta}^{\zeta} = \frac{-\tan(\theta)}{2},$$

$$\Gamma_{\phi\zeta}^{\theta} = \Gamma_{\zeta\phi}^{\theta} = \frac{-\cos(\theta)}{2},$$

with the remaining ones equal to 0. Using the Christoffel symbols, the covariant derivatives are as follows:

$$\nabla_{\partial_a} \partial_b = \Gamma^c_{ab} \partial_c,$$

which, in our case, implies:

$$\nabla_{\partial_{\phi}} \partial_{\phi} = 0,$$

$$\nabla_{\partial_{\theta}} \partial_{\theta} = 0,$$

$$\nabla_{\partial_{\zeta}} \partial_{\zeta} = 0,$$

$$\nabla_{\partial_{\phi}} \partial_{\theta} = \nabla_{\partial_{\theta}} \partial_{\phi} = -\frac{\tan(\theta)}{2} \partial_{\phi} + \frac{1}{2\cos(\theta)} \partial_{\zeta},$$

$$\nabla_{\partial_{\zeta}} \partial_{\theta} = \nabla_{\partial_{\theta}} \partial_{\zeta} = \frac{1}{2\cos(\theta)} \partial_{\zeta} - \frac{\tan(\theta)}{2} \partial_{\phi},$$

$$\nabla_{\partial_{\phi}} \partial_{\zeta} = \nabla_{\partial_{\zeta}} \partial_{\phi} = -\frac{\cos(\theta)}{2} \partial_{\theta}.$$

These covariant derivatives can then be employed to compute the components of the Hessian according to:

$$\mathcal{H}_{ij}(f) = \nabla^2 f(e_i, e_j) = e_i e_j f - \nabla_{e_i} e_j f.$$

For the first part of each term, we obtain:

$$\begin{aligned} e_1 e_1 &= \left( \frac{\sqrt{1 - \sin\theta} + \sqrt{\sin\theta + 1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\phi} + \frac{\sqrt{1 - \sin\theta} - \sqrt{\sin\theta + 1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta} \right) \\ &\quad \times \left( \frac{\sqrt{1 - \sin\theta} + \sqrt{\sin\theta + 1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\phi} + \frac{\sqrt{1 - \sin\theta} - \sqrt{\sin\theta + 1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta} \right) \\ &= \frac{1}{4\cos^2(\theta)} \left[ (1 + \cos(\theta)) \frac{\partial^2}{\partial\phi^2} + (1 - \cos(\theta)) \frac{\partial^2}{\partial\zeta^2} - 2\sin(\theta) \frac{\partial^2}{\partial\zeta\partial\phi} \right], \end{aligned}$$
$$e_2 e_2 = \left( \frac{1}{\sqrt{2}} \frac{\partial}{\partial \theta} \right) \left( \frac{1}{\sqrt{2}} \frac{\partial}{\partial \theta} \right) = \frac{\partial^2}{2 \partial \theta^2},$$

$$e_3 e_3 = \left( \frac{\sqrt{1 - \sin \theta} - \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \phi} + \frac{\sqrt{1 - \sin \theta} + \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \zeta} \right) \\ \times \left( \frac{\sqrt{1 - \sin \theta} - \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \phi} + \frac{\sqrt{1 - \sin \theta} + \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \zeta} \right) \\ = \frac{1}{4 \cos^2(\theta)} \left[ (1 - \cos(\theta)) \frac{\partial^2}{\partial \phi^2} + (1 + \cos(\theta)) \frac{\partial^2}{\partial \zeta^2} - 2 \sin(\theta) \frac{\partial^2}{\partial \zeta \partial \phi} \right],$$

$$e_1 e_2 = e_2 e_1 = \left( \frac{\sqrt{1 - \sin \theta} + \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \phi} + \frac{\sqrt{1 - \sin \theta} - \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \zeta} \right) \left( \frac{1}{\sqrt{2}} \frac{\partial}{\partial \theta} \right) \\ = \frac{1}{4 \cos(\theta)} \left[ \left( \sqrt{1 - \sin(\theta)} - \sqrt{\sin(\theta) + 1} \right) \frac{\partial^2}{\partial \theta \partial \zeta} + \left( \sqrt{1 - \sin(\theta)} + \sqrt{\sin(\theta) + 1} \right) \frac{\partial^2}{\partial \theta \partial \phi} \right],$$

$$e_1 e_3 = e_3 e_1 = \left( \frac{\sqrt{1 - \sin \theta} + \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \phi} + \frac{\sqrt{1 - \sin \theta} - \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \zeta} \right) \\ \times \left( \frac{\sqrt{1 - \sin \theta} - \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \phi} + \frac{\sqrt{1 - \sin \theta} + \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \zeta} \right) \\ = \frac{1}{4 \cos^2(\theta)} \left[ -\sin(\theta) \frac{\partial^2}{\partial \phi^2} - \sin(\theta) \frac{\partial^2}{\partial \theta^2} + 2 \frac{\partial^2}{\partial \phi \partial \zeta} \right],$$

$$e_2 e_3 = e_2 e_3 = \left( \frac{1}{\sqrt{2}} \frac{\partial}{\partial \theta} \right) \left( \frac{\sqrt{1 - \sin \theta} - \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \phi} + \frac{\sqrt{1 - \sin \theta} + \sqrt{\sin \theta + 1}}{2\sqrt{2} \cos \theta} \frac{\partial}{\partial \zeta} \right) \\ = \frac{1}{4 \cos(\theta)} \left[ \left( \sqrt{1 - \sin(\theta)} - \sqrt{\sin(\theta) + 1} \right) \frac{\partial^2}{\partial \theta \partial \zeta} + \left( \sqrt{1 - \sin(\theta)} + \sqrt{\sin(\theta) + 1} \right) \frac{\partial^2}{\partial \theta \partial \phi} \right],$$

while the second part of each term is given by:

$$\nabla_{e_1} e_1 = \nabla \frac{\sqrt{1-\sin\theta}+\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\phi} + \frac{\sqrt{1-\sin\theta}-\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta} \left( \frac{\sqrt{1-\sin\theta}+\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\phi} + \frac{\sqrt{1-\sin\theta}-\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta} \right) = \frac{\sin(2\theta)}{8\cos^2\theta} \frac{\partial}{\partial\theta},$$

$$\nabla_{e_2} e_2 = \nabla \frac{1}{\sqrt{2}} \frac{\partial}{\partial \theta} \left( \frac{1}{\sqrt{2}} \frac{\partial}{\partial \theta} \right) = 0,$$

$$\nabla_{e_3} e_3 = \nabla \frac{\sqrt{1-\sin\theta}-\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\phi} + \frac{\sqrt{1-\sin\theta}+\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta} \left( \frac{\sqrt{1-\sin\theta}-\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\phi} + \frac{\sqrt{1-\sin\theta}+\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta} \right) = \frac{\sin(2\theta)}{8\cos^2\theta} \frac{\partial}{\partial\theta},$$

$$\nabla_{e_1} e_2 = \nabla_{e_2} e_1 = \nabla \frac{\sqrt{1-\sin\theta}+\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\phi} + \frac{\sqrt{1-\sin\theta}-\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta} \left( \frac{1}{\sqrt{2}} \frac{\partial}{\partial\theta} \right) \\ = \frac{1}{8 \cos^2(\theta)} \left[ (P(\theta) \sin(\theta) - M(\theta)) \frac{\partial}{\partial\phi} + (M(\theta) \sin(\theta) - P(\theta)) \frac{\partial}{\partial\zeta} \right],$$
$$ \nabla_{e_1} e_3 = \nabla_{e_3} e_1 = \nabla_{\frac{\sqrt{1-\sin\theta}+\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta}} \left( \frac{\sqrt{1-\sin\theta}-\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta} + \frac{\sqrt{1-\sin\theta}+\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta} \right) = -\frac{1}{4\cos(\theta)} \frac{\partial}{\partial\theta}, $$

$$ \begin{aligned} \nabla_{e_2} e_3 = \nabla_{e_3} e_2 &= \nabla \left( \frac{\sqrt{1-\sin\theta}-\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\phi} + \frac{\sqrt{1-\sin\theta}+\sqrt{\sin\theta+1}}{2\sqrt{2}\cos\theta} \frac{\partial}{\partial\zeta} \right) \left( \frac{1}{\sqrt{2}} \frac{\partial}{\partial\theta} \right) \\ &= \frac{1}{8\cos^2(\theta)} \left[ (M(\theta)\sin(\theta) - P(\theta)) \frac{\partial}{\partial\phi} + (P(\theta)\sin(\theta) - M(\theta)) \frac{\partial}{\partial\zeta} \right], \end{aligned} $$

where:

$$ \begin{aligned} P(\theta) &:= \sqrt{1 - \sin(\theta)} + \sqrt{\sin(\theta) + 1}, \\ M(\theta) &:= \sqrt{1 - \sin(\theta)} - \sqrt{\sin(\theta) + 1}. \end{aligned} $$

Therefore, the components of the Hessian are the following:

$$ \mathcal{H}_{11} = \frac{1}{4 \cos^2(\theta)} \left[ (1 + \cos(\theta)) \frac{\partial^2}{\partial \phi^2} + (1 - \cos(\theta)) \frac{\partial^2}{\partial \zeta^2} - 2 \sin(\theta) \frac{\partial^2}{\partial \zeta \partial \phi} - \frac{\sin(2\theta)}{2} \frac{\partial}{\partial \theta} \right], $$

$$ \mathcal{H}_{22} = \frac{\partial^2}{2\partial\theta^2}, $$

$$ \mathcal{H}_{33} = \frac{1}{4 \cos^2(\theta)} \left[ (1 - \cos(\theta)) \frac{\partial^2}{\partial \phi^2} + (1 + \cos(\theta)) \frac{\partial^2}{\partial \zeta^2} - 2 \sin(\theta) \frac{\partial^2}{\partial \zeta \partial \phi} - \frac{\sin(2\theta)}{2} \frac{\partial}{\partial \theta} \right], $$

$$ \mathcal{H}_{12} = \mathcal{H}_{21} = \frac{1}{8 \cos^2(\theta)} \left[ 2 \cos(\theta) \left( M(\theta) \frac{\partial^2}{\partial\theta\partial\zeta} + P(\theta) \frac{\partial^2}{\partial\theta\partial\phi} \right) + (P(\theta)\sin(\theta) - M(\theta)) \frac{\partial}{\partial\phi} + (M(\theta)\sin(\theta) - P(\theta)) \frac{\partial}{\partial\zeta} \right], $$

$$ \mathcal{H}_{13} = \mathcal{H}_{31} = \frac{1}{4 \cos^2(\theta)} \left[ -\sin(\theta) \frac{\partial^2}{\partial \phi^2} - \sin(\theta) \frac{\partial^2}{\partial \zeta^2} + 2 \frac{\partial^2}{\partial \zeta \partial \phi} + \cos(\theta) \frac{\partial}{\partial \theta} \right], $$

$$ \mathcal{H}_{23} = \mathcal{H}_{32} = \frac{1}{8 \cos^2(\theta)} \left[ 2 \cos(\theta) \left( M(\theta) \frac{\partial^2}{\partial\theta\partial\phi} + P(\theta) \frac{\partial^2}{\partial\theta\partial\zeta} \right) + (M(\theta)\sin(\theta) - P(\theta)) \frac{\partial}{\partial\phi} + (P(\theta)\sin(\theta) - M(\theta)) \frac{\partial}{\partial\zeta} \right]. $$

## References

[1] M. Kamionkowski, A. Kosowsky and A. Stebbins, A Probe of primordial gravity waves and vorticity, Phys. Rev. Lett. 78 (1997) 2058 [astro-ph/9609132] [INSPIRE].

[2] M. Zaldarriaga and U. Seljak, An all sky analysis of polarization in the microwave background, Phys. Rev. D 55 (1997) 1830 [astro-ph/9609170] [INSPIRE].

[3] N. Bartolo, E. Komatsu, S. Matarrese and A. Riotto, Non-Gaussianity from inflation: Theory and observations, Phys. Rep. 402 (2004) 103 [astro-ph/0406398] [INSPIRE].
[4] C. Pitrou, T.S. Pereira and J.-P. Uzan, *Predictions from an anisotropic inflationary era*, JCAP **04** (2008) 004 [arXiv:0801.3596] [INSPIRE].

[5] J.D. Barrow and S. Hervik, *Anisotropically inflating universes*, Phys. Rev. D **73** (2006) 023007 [gr-qc/0511127] [INSPIRE].

[6] D.N. Spergel and D.M. Goldberg, *Microwave background bispectrum. Part 1. Basic formalism*, Phys. Rev. D **59** (1999) 103001 [astro-ph/9811252] [INSPIRE].

[7] C.T. Byrnes, M. Gerstenlauer, S. Nurmi, G. Tasinato and D. Wands, *Scale-dependent non-Gaussianity probes inflationary physics*, JCAP **10** (2010) 004 [arXiv:1007.4277] [INSPIRE].

[8] D. Pogosyan, C. Pichon and C. Gay, *Non Gaussian extrema counts for CMB maps*, Phys. Rev. D **84** (2011) 083510 [arXiv:1107.1863] [INSPIRE].

[9] J. Carrón Duque, A. Buzzelli, Y. Fantaye, D. Marinucci, A. Schwartzman and N. Vittorio, *Point Source Detection and False Discovery Rate Control on CMB Maps*, Astron. Comput. **28** (2019) 100310 [arXiv:1902.06636] [INSPIRE].

[10] A. Vafaee Sadr and S.M.S. Movahed, *Clustering of local extrema in Planck CMB maps*, Mon. Not. Roy. Astron. Soc. **503** (2021) 815 [arXiv:2003.07364] [INSPIRE].

[11] J. Kasak, J. Creswell, H. Liu and P. Naselsky, *Statistics of nonpolarized points in the CMB polarization maps*, Phys. Rev. D **104** (2021) 023502 [arXiv:2012.15811] [INSPIRE].

[12] D.G. Figueroa, S. Raatikainen, S. Räsänen and E. Tomberg, *Non-Gaussian Tail of the Curvature Perturbation in Stochastic Ultraslow-Roll Inflation: Implications for Primordial Black Hole Production*, Phys. Rev. Lett. **127** (2021) 101302 [arXiv:2012.06551] [INSPIRE].

[13] D. Cruces, *Review on Stochastic Approach to Inflation*, Universe **8** (2022) 334 [arXiv:2203.13852] [INSPIRE].

[14] J. Schmalzing and K.M. Górski, *Minkowski functionals used in the morphological analysis of cosmic microwave background anisotropy maps*, Mon. Not. Roy. Astron. Soc. **297** (1998) 355 [astro-ph/9710185] [INSPIRE].

[15] PLANCK collaboration, *Planck 2018 results. Part VII. Isotropy and Statistics of the CMB*, Astron. Astrophys. **641** (2020) A7 [arXiv:1906.02552] [INSPIRE].

[16] M. Vicinanza, V.F. Cardone, R. Maoli, R. Scaramella, X. Er and I. Tereno, *Minkowski functionals of convergence maps and the lensing figure of merit*, Phys. Rev. D **99** (2019) 043534 [arXiv:1905.00410] [INSPIRE].

[17] J.M. Zorrilla Matilla, M. Sharma, D. Hsu and Z. Haiman, *Interpreting deep learning models for weak lensing*, Phys. Rev. D **102** (2020) 123506 [arXiv:2007.06529] [INSPIRE].

[18] D. Zürcher, J. Fluri, R. Sgier, T. Kacprzak and A. Refregier, *Cosmological Forecast for non-Gaussian Statistics in large-scale weak Lensing Surveys*, JCAP **01** (2021) 028 [arXiv:2006.12506] [INSPIRE].

[19] C. Hervías-Caimapo and K.M. Huffenberger, *Full-sky, Arcminute-scale, 3D Models of Galactic Microwave Foreground Dust Emission Based on Filaments*, Astrophys. J. **928** (2022) 65 [arXiv:2107.08317] [INSPIRE].

[20] F. Rahman, P. Chingangbam and T. Ghosh, *The nature of non-Gaussianity and statistical isotropy of the 408 MHz Haslam synchrotron map*, JCAP **07** (2021) 026 [arXiv:2104.00419] [INSPIRE].

[21] F.A. Martire, A.J. Banday, E. Martínez-González and R.B. Barreiro, *Morphological analysis of the polarized synchrotron emission with WMAP and Planck*, JCAP **04** (2023) 049 [arXiv:2301.08041] [INSPIRE].
[22] V. Ganesan, P. Chingangbam, K.P. Yogendran and C. Park, *Primordial non-Gaussian signatures in CMB polarization*, JCAP 02 (2015) 028 [arXiv:1411.5256] [INSPIRE].

[23] L. Santos, K. Wang and W. Zhao, *Probing the statistical properties of CMB B-mode polarization through Minkowski Functionals*, JCAP 07 (2016) 029 [arXiv:1510.07779] [INSPIRE].

[24] P. Chingangbam, V. Ganesan, K.P. Yogendran and C. Park, *On Minkowski Functionals of CMB polarization*, Phys. Lett. B 771 (2017) 67 [arXiv:1705.04454] [INSPIRE].

[25] G. Puglisi and X. Bai, *Inpainting Galactic Foreground Intensity and Polarization Maps Using Convolutional Neural Networks*, Astrophys. J. 905 (2020) 143 [arXiv:2003.13691] [INSPIRE].

[26] N. Krachmalnicoff and G. Puglisi, *ForSE: A GAN-based Algorithm for Extending CMB Foreground Models to Subdegree Angular Scales*, Astrophys. J. 911 (2021) 42 [arXiv:2011.02221] [INSPIRE].

[27] A. Carones, J. Carrón Duque, D. Marinucci, M. Migliaccio and N. Vittorio, *Minkowski Functionals of CMB polarisation intensity with Pynkowski: theory and application to Planck data*, Mon. Not. Roy. Astron. Soc. 527 (2024) 756 [arXiv:2211.07562] [INSPIRE].

[28] J.R. Gott III, A.L. Melott and M. Dickinson, *The Sponge-like topology of large-scale structure in the Universe*, Astrophys. J. 306 (1986) 341 [INSPIRE].

[29] H. Tomita, *Curvature Invariants of Random Interface Generated by Gaussian Fields*, Prog. Theor. Phys. 76 (1986) 952.

[30] T. Matsubara, *Statistics of Smoothed Cosmic Fields in Perturbation Theory. Part 1. Formulation and Useful Formulae in Second Order Perturbation Theory*, Astrophys. J. 584 (2003) 1 [INSPIRE].

[31] C. Hikage, E. Komatsu and T. Matsubara, *Primordial Non-Gaussianity and Analytical Formula for Minkowski Functionals of the Cosmic Microwave Background and Large-scale Structure*, Astrophys. J. 653 (2006) 11 [astro-ph/0607284] [INSPIRE].

[32] M. Stecconi, *Isotropic random spin weighted functions on $S^2$ vs. isotropic random fields on $S^3$*, arXiv:2108.00736.

[33] A. Lerario, D. Marinucci, M. Rossi and M. Stecconi, *Geometry and topology of spin random fields*, arXiv:2207.08413.

[34] R.J. Adler and J.E. Taylor, *Random fields and geometry*, in Springer Monographs in Mathematics, Springer (2007) [DOI:10.1007/978-0-387-48116-6].

[35] Y. Fantaye, D. Marinucci, F. Hansen and D. Maino, *Applications of the Gaussian kinematic formula to CMB data analysis*, Phys. Rev. D 91 (2015) 063501.

[36] S.-S. Chern, *A Simple Intrinsic Proof of the Gauss-Bonnet Formula for Closed Riemannian Manifolds*, Ann. Math. 45 (1944) 747.

[37] Y.T. Fantaye and D. Marinucci, *Cosmological Applications of the Gaussian Kinematic Formula*, Proc. Int. Astron. Union 10 (2014) 135 [arXiv:1406.5420].

[38] K.M. Górski et al., *HEALPix — A Framework for high resolution discretization, and fast analysis of data distributed on the sphere*, Astrophys. J. 622 (2005) 759 [astro-ph/0409513] [INSPIRE].

[39] M.P. do Carmo, *Differential geometry of curves and surfaces*, Prentice Hall (1976).

[40] X. Lan, D. Marinucci and Y. Xiao, *Strong local nondeterminism and exact modulus of continuity for spherical Gaussian fields*, Stochastic Processes Appl. 128 (2018) 1294 [arXiv:1603.03642].

[41] A. Zonca et al., *healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in Python*, J. Open Source Softw. 4 (2019) 1298 [INSPIRE].
[42] PLANCK collaboration, Planck 2018 results. Part V. CMB power spectra and likelihoods, Astron. Astrophys. 641 (2020) A5 [arXiv:1907.12875] [INSPIRE].

[43] PLANCK collaboration, Planck 2018 results. Part VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6 [Erratum ibid. 652 (2021) C4] [arXiv:1807.06209] [INSPIRE].

[44] S. Aiola et al., The Atacama Cosmology Telescope: DR4 Maps and Cosmological Parameters, JCAP 12 (2020) 047 [arXiv:2007.07288] [INSPIRE].

[45] SPT collaboration, Measurements of B-mode Polarization of the Cosmic Microwave Background from 500 Square Degrees of SPTpol Data, Phys. Rev. D 101 (2020) 122003 [arXiv:1910.05748] [INSPIRE].

[46] P. Ade et al., The Simons Observatory: Science goals and forecasts, JCAP 02 (2019) 056 [arXiv:1808.07445] [INSPIRE].

[47] T. Matsumura et al., Mission design of LiteBIRD, J. Low Temp. Phys. 176 (2014) 733 [arXiv:1311.2847] [INSPIRE].

[48] S. Hanany et al., PICO: Probe of Inflation and Cosmic Origins, arXiv:1902.10541 [INSPIRE].

[49] I.M. Gelfand, R.A. Minlos and Z.Y. Shapiro, Representations of the Rotation and Lorentz Groups and their Applications, Courier Dover Publications (1963).