# Minkowski Functionals of CMB polarization intensity with Pynkowski: Theory and application to Planck and future data

Alessandro Carones, Javier Carrón Duque, Domenico Marinucci, Marina Migliaccio, Nicola Vittorio

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

## ► To cite this version:

Alessandro Carones, Javier Carrón Duque, Domenico Marinucci, Marina Migliaccio, Nicola Vittorio. Minkowski Functionals of CMB polarization intensity with Pynkowski: Theory and application to Planck and future data. Monthly Notices of the Royal Astronomical Society, 2023, 527, pp.756 - 773. 10.1093/mnras/stad3002. hal-04863108

HAL Id: hal-04863108

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

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Minkowski functionals of CMB polarization intensity with PYNKOWSKI: theory and application to Planck and future data

Alessandro Carones,$^{1,2*}$ Javier CarrónDuque,$^{1,2}$ Domenico Marinucci,$^3$ Marina Migliaccio$^{1,2}$ and Nicola Vittorio$^{1,2}$

Accepted 2023 September 26. Received 2023 September 7; in original form 2023 February 9

## ABSTRACT

The angular power spectrum of the cosmic microwave background (CMB) anisotropies is a key tool to study the Universe. However, it is blind to the presence of non-Gaussianities and deviations from statistical isotropy, which can be detected with other statistics such as Minkowski functionals (MFs). These tools have been applied to CMB temperature and E-mode anisotropies with no detection of deviations from Gaussianity and isotropy. In this work, we extend the MF formalism to the CMB polarization intensity, $P^2 = Q^2 + U^2$. We use the Gaussian kinematic formula to derive the theoretical predictions of MFs for Gaussian isotropic fields. We develop a software that computes MFs on $P^2$ HEALPIX maps and apply it to simulations to verify the robustness of both theory and methodology. We then estimate MFs of $P^2$ maps from Planck, both in pixel space and needlet domain, comparing them with realistic simulations that include CMB and instrumental noise residuals. We find no significant deviations from Gaussianity or isotropy in Planck CMB polarization intensity. However, MFs could play an important role in the analysis of CMB polarization measurements from upcoming experiments with improved sensitivity. Therefore, we forecast the ability of MFs applied to $P^2$ maps to detect much fainter non-Gaussian anisotropic signals than with Planck data for two future complementary experiments: the LiteBIRD satellite and the ground-based Simons Observatory. We publicly release the software to compute MFs in arbitrary scalar HEALPIX maps as a fully documented PYTHON package called PYNKOWSKI.

**Key words:** methods: statistical – software: public release – cosmic background radiation – cosmology: observations.

## 1 INTRODUCTION

A wide variety of powerful statistical and geometrical tools have been adopted for the analysis of cosmological data. Among them, the angular power spectrum of the anisotropies of the cosmic microwave background (CMB), both in temperature and polarization, has played a key role in constraining the values of the parameters of the Lambda cold dark matter cosmological model (Planck Collaboration VI 2020c). The angular power spectrum contains all the statistical information of a Gaussian isotropic field, but it is blind to the presence of non-Gaussianities or physical anisotropies in the data; in other words, two observed maps can have the same angular power spectrum but still have widely different statistical properties.

Searching for the presence of small non-Gaussianities or deviations from statistical isotropy in the primordial Universe is an active field of research, as they could encode information about cosmic inflation and physics of the very early Universe, if detected (Bartolo et al. 2004; Barrow & Hervik 2006; Pitrou, Pereira & Uzan 2008; Chen 2010). Such effects, however, are not associated only with the physics of the primordial Universe: CMB maps

are contaminated by residual foregrounds, complex instrumental systematics, and non-primordial cosmological effects such as CMB lensing, which can induce significant deviations from the Gaussianity and isotropy hypotheses. These signals contain relevant information that is not reflected in the angular power spectrum. Non-Gaussian statistics of CMB lensing have already been shown to help constrain cosmological parameters (Chen et al. 2020; Zürcher et al. 2021; Grewal et al. 2022).

Recently, many studies about primordial non-Gaussianity have also involved the analysis of the large-scale structure of the Universe (Peña & Candlish 2022; Andrews et al. 2023). In particular, the high-order statistics of galaxy clustering with future galaxy surveys could improve the constraints on primordial non-Gaussianity (Desjacques & Seljak 2010).

Minkowski functionals (MFs) are one of the most popular tools used to analyse information complementary to that contained in the angular power spectrum. MFs have been widely explored in the mathematical literature (Tomita 1986; Coles & Barrow 1987; Gott III et al. 1990), and have become a popular tool in CMB analysis after the discussion by Schmalzing & Górski (1998). An extension of the MF formalism to study the full information embedded in the spin-2 CMB polarization field has recently been proposed in Carrón Duque et al. (2023).

<footnote>$^1$Dipartimento di Fisica, Università di Roma 'Tor Vergata', via della Ricerca Scientifica 1, I-00133 Roma, Italy</footnote>

<footnote>$^2$Sezione INFN Roma 2, via della Ricerca Scientifica 1, I-00133 Roma, Italy</footnote>

<footnote>$^3$Dipartimento di Matematica, Universita' di Roma Tor Vergata, via della Ricerca Scientifica 1, I-00133 Roma, Italy</footnote>

<footnote>* E-mail: alessandro.carones@roma2.infn.it</footnote>
In the case of scalar Gaussian isotropic fields such as the CMB temperature anisotropies, the theoretical predictions of these functionals can be accurately computed. Given their low variance, any deviation from Gaussianity or statistical isotropy can be detected at high significance in a model-independent manner. MFs are now widely used to measure possible deviations from primordial Gaussianity (Planck Collaboration VII 2020d) and, beyond that, e.g. to constrain cosmological parameters (Zürcher et al. 2021). Moreover, these tools can be applied to detect contamination by foreground and instrumental systematics in the reconstructed CMB map in a fully blind fashion, given their highly non-Gaussian and anisotropic nature. The impact of such systematic and astrophysical effects on the constraining power for primordial non-Gaussianity in temperature data has also been explored by Ducout et al. (2012). Up to now, MFs have been largely applied to the scalar fields associated with the CMB polarization, E modes and B modes (Ganesan et al. 2015; Planck Collaboration XVI 2016; Santos, Wang & Zhao 2016; Planck Collaboration VII 2020d), while they remain mostly unexploited in the analysis of spin fields. Therefore, in this work, we aim to extend the formalism of MFs to the squared modulus of the CMB polarization, $P^2 = Q^2 + U^2$, which hereafter we will refer to simply as polarization modulus or polarized intensity.

This new formalism conveys complementary information on the CMB polarization field with respect to the study of E and B modes in several aspects:

(i) The polarization modulus traces the local properties of the polarization field, while the E- and B-mode maps are obtained from a non-local decomposition of $Q$ and $U$, typically performed in harmonic space.

(ii) The $P^2$ field can be robustly obtained for any adopted mask without suffering from E-B and B-E leakage contamination (Lewis, Challinor & Turok 2001; Bunn et al. 2003), which may distort the statistics of the reconstructed E and B modes.

(iii) The morphology of the polarization modulus is intrinsically different and can be more effective in detecting those deviations from Gaussianity and statistical isotropy that primarily affect the intensity rather than the direction of the polarized emission.

We note that, although E- and B-mode maps embed the full information on the polarization field (in the full-sky case), analysing MFs of these fields independently does not exhaust the morphological characterization of the polarization field. This motivates the usefulness of analysing complementary quantities as $P^2$.

This paper has five key goals:

(i) To introduce a more general formalism to compute the theoretical predictions of scalar maps, such as the CMB temperature anisotropies (T).

(ii) To apply this extended formalism to obtain the expected values of the MFs of the modulus of a spin field, such as the CMB polarization ($P^2$) in the case of Gaussianity and statistical isotropy.

(iii) To provide the community with a public PYTHON package called PYNKOWSKI<fnref n="1" /> to easily estimate MFs on HEALPIX<fnref n="2" /> maps (Górski et al. 2005) and compute their Gaussian isotropic expectations.

(iv) To compute the MFs on Planck $P^2$ maps and compare them with those from realistic simulations in order to assess any possible unexpected deviation from Gaussianity or isotropy.

(v) To forecast the sensitivity of MFs applied to $P^2$ maps to deviations from Gaussianity and statistical isotropy for future CMB experiments.

We note that, up to second-order terms, theoretical predictions of MFs for the sum of two squared Gaussian fields have already been presented in the literature (see Naselsky & Novikov 1998, and the references therein). However, in this work, these results are further justified and are derived by means of a more general approach. Furthermore, for the first time, these predictions are compared with realistic CMB simulations, and the MFs of the $P^2$ field are estimated on Planck data. These CMB polarization maps are dominated by the anisotropic noise; however, we expect these tools to grow in importance for the analysis of data with much greater sensitivity from ongoing or future suborbital and satellite experiments such as Atacama Cosmology Telescope (Aiola et al. 2020), South Pole Telescope (SPT, Sayre et al. 2020), Large Scale Polarization Explorer (LSPE, Addamo et al. 2021), Simons Observatory (SO; Ade et al. 2019), CMB Stage-4 (CMB-S4, Abazajian et al. 2022), and LiteBIRD (Matsumura et al. 2014). To further investigate this point, we forecast the capability of MFs computed on simulated realistic $P^2$ maps to detect deviations from Gaussianity and statistical isotropy for two future and complementary CMB polarization experiments: SO and LiteBIRD. Both of them will target the detection of the so-called primordial CMB polarization B modes, supposed to be generated by a stochastic background of gravitational waves from cosmic inflation. SO is ground-based and located in the Atacama Desert in Chile. Due to its limited observed sky fraction ($f_{sky} \sim 0.3$) and the atmospheric contamination, SO is sensitive to polarization modes with $\ell > 30$ and thus will only target the detection of the recombination peak in the primordial B-mode power spectrum at $\ell \sim 100$ (Ade et al. 2019). LiteBIRD, instead, is a space-borne L-class mission selected by the Japanese Aerospace Exploration Agency in 2019. Given its full-sky coverage, it will be sensitive also to modes on the largest angular scales (LiteBIRD Collaboration 2023). For both experiments, the detection of primordial CMB B modes will be hampered by residual contamination of the Galactic emission, which is highly non-Gaussian and anisotropic. Therefore, we have included such a signal in the polarization simulated data of LiteBIRD and SO Small Aperture Telescopes (SO-SATs) as non-Gaussian anisotropic component.

In this paper, we also present, to our knowledge, the first public software needed to compute these MFs on any HEALPIX map. The structure of this article is as follows. In Section 2, we present the definition of MFs and their application to stochastic fields defined on the sphere. In Section 3, we derive the theoretical expectations of the MFs for CMB polarized intensity: $P^2 = Q^2 + U^2$. In Section 4, we present the PYNKOWSKI PYTHON package and explain the procedures implemented to estimate these statistical quantities on the maps. In Section 5, we introduce the data and simulations we use in this paper. In Section 6, we present the results of applying this framework and software to CMB polarization simulations and Planck observations, as well as the forecasted sensitivity of $P^2$ MFs to a non-Gaussian anisotropic signal for LiteBIRD and SO. Finally, in Section 7, we report our conclusions.

## 2 MINKOWSKI FUNCTIONALS ON THE SPHERE

MFs on the sphere are now well-established tools for CMB data analysis (Planck Collaboration VII 2020d); none the less, we will first recall some background notation and definitions in order to make the comparison with this work more explicit. Given a scalar random

<footnote n="1">¹https://github.com/javicarron/pynkowski</footnote>

<footnote n="2">²https://github.com/healpy/healpy</footnote>
field $f(.)$ observed on the unit sphere $\mathbb{S}^2$ (e.g. CMB temperature anisotropies), the excursion set at a threshold $u$ is defined by

$$A_u(f, .), \mathbb{S}^2) = \{x \in \mathbb{S}^2 : f(x) \geq u\}.$$

We will omit the arguments of $A_u$ to refer to the excursion sets when $f$ and the domain of $f$ are clear by context.

As recalled in Schmalzing & Górski (1998), the morphological properties of these excursion sets can be summarized by the three MFs defined as follows:

$$ 
\begin{aligned}
V_0(A_u) &= \int_{A_u} dx, \\
V_1(A_u) &= \frac{1}{4} \int_{\partial A_u} dr, \\
V_2(A_u) &= \frac{1}{2\pi} \int_{\partial A_u} \kappa(r) dr,
\end{aligned}
$$

where $dr$ denotes a line element along the boundary of the excursion set and $\kappa(.)$ denotes the geodesic curvature of said boundary. Thus, $V_0$ represents the excursion area, while $V_1$ is one-fourth of the boundary length. $V_2$ is associated with the number of connected regions minus the number of holes, which on the plane corresponds to the Euler-Poincaré characteristic $\chi$. On a spherical surface, the Euler-Poincaré characteristic does no longer correspond exactly to $V_2$, but includes an extra term that can be computed by means of the generalized Gauss-Bonnet theorem, which yields

$$\chi(A_u) = V_2(A_u) + \frac{1}{2\pi} V_0(A_u).$$

From the mathematical point of view, it turns out to be notationally more convenient to replace MFs with the equivalent notion of Lipschitz-Killing curvatures, which are defined by

$$\mathcal{L}_2(A_u) = V_0(A_u), \qquad (1)$$

$$\mathcal{L}_1(A_u) = 2V_1(A_u), \qquad (2)$$

$$\mathcal{L}_0(A_u) = \chi(A_u). \qquad (3)$$

These quantities are defined as follows. Let us define a tube of width $\rho$ built around the manifold $A_u$ as all those 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 Adler & Taylor 2007; Fantaye & Marinucci 2014, for a detailed discussion).

Let us now recall how to compute the expected values of these statistical quantities for Gaussian isotropic maps. These results are well known in the cosmological literature in the case of the two-dimensional sphere, but it is convenient to recall them from a more general perspective.

Without loss of generality, we assume the field normalized to have unit variance. We define $\mu$ as the derivative of the covariance function at the origin for the field $f$. This quantity can be computed as follows:

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

where the sequence $\{C_{\ell}\}$ denotes, as usual, the angular power spectrum of the field $f$. Furthermore, following Adler & Taylor (2007), we define the 'flag' coefficients as

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

with $\omega_1 = 2$, $\omega_2 = \pi$, and $\omega_3 = \frac{4\pi}{3}$. In general, $\omega_j$ represents the volume of the $j$-dimensional unit ball. For an isotropic Gaussian field, the expected value of the Lipschitz-Killing curvatures on the sphere is then given by the Gaussian kinematic formula (GKF):

$$\mathbb{E}[\mathcal{L}_j(A_u)] = \sum_{k=0}^{2-j} \begin{bmatrix} k+j \\ k \end{bmatrix} \rho_k(u) \mathcal{L}_{k+j}(\mathbb{S}^2) \mu^{k/2}, \qquad (4)$$

where

$$ 
\begin{aligned}
\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,
\end{aligned}
$$

and in general, $H_k(u) = (-1)^k \exp\left(\frac{u^2}{2}\right) \frac{d^k}{du^k} \exp\left(-\frac{u^2}{2}\right)$ denotes the sequence of Hermite polynomials. The function $\Phi$ is the cumulative normal distribution.

In particular, by equation (4), we get the well-known results for the expected values of the Lipschitz-Killing curvatures and, using equation (3), of MFs, which we report normalized by area. In order, we have

$$ 
\begin{aligned}
\mathbb{E}[\mathcal{L}_2(A_u)] &= \rho_0(u)\mathcal{L}_2(\mathbb{S}^2) = 4\pi\{1 - \Phi(u)\}, \\
\frac{\mathbb{E}[V_0]}{4\pi} &= 1 - \Phi(u).
\end{aligned}
$$

We note that this quantity depends only on $u$, not on the field itself as long as it is normalized to have unit variance.

Likewise,

$$ 
\begin{aligned}
\mathbb{E}[\mathcal{L}_1(A_u)] &= 2 \frac{\omega_2}{\omega_1^2} \rho_1(u) \mathcal{L}_2(\mathbb{S}^2) \mu^{1/2} \\
&= \frac{2\pi}{4} \frac{1}{(2\pi)^{1/2}} \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{u^2}{2}\right) 4\pi \mu^{1/2} \\
&= \exp\left(-\frac{u^2}{2}\right) \pi \mu^{1/2}
\end{aligned}
$$

so that

$$\frac{\mathbb{E}[V_1]}{4\pi} = \frac{1}{8} \exp\left(-\frac{u^2}{2}\right) \mu^{1/2}.$$

Finally,

$$ 
\begin{aligned}
\mathbb{E}[\mathcal{L}_0(A_u)] &= \rho_0(u)\mathcal{L}_0(\mathbb{S}^2) + \rho_2(u)\mathcal{L}_2(\mathbb{S}^2)\mu \\
&= \{1 - \Phi(u)\} 2 + \frac{2}{\sqrt{2\pi}} \mu \exp\left(-\frac{u^2}{2}\right) u,
\end{aligned}
$$

and hence

$$\frac{\mathbb{E}[V_2]}{4\pi} = \frac{1}{\sqrt{(2\pi)^3}} \mu \exp\left(-\frac{u^2}{2}\right) u.$$

The interpretation of the previous formulae is worth discussing – the expected values of the MFs (and Lipschitz-Killing curvatures) are expressed by means of a product of four different independent ingredients:

(i) a set of universal coefficients (called 'flag' coefficients),

(ii) the Lipschitz-Killing curvatures evaluated on the original manifold (in our case, the unitary sphere),

(iii) a power of the derivative of the covariance function of the field at the origin, and
(iv) a set of 'density' functions $\rho_k$, dependent only on the threshold at which the MFs (or Lipschitz-Killing curvatures) are evaluated.

The remarkable fact is that in the case of polarization, the expected values for Lipschitz-Killing curvatures take exactly the same form – the only change is in the 'density' functions $\rho_k$, which none the less can be expressed with a fully explicit analytical form. We give more details in the section to follow.

We also note that the presence of a mask only affects the manifold where the map is defined, and therefore, at leading order, it has an impact only on the normalization of the MFs.

### 3 MINKOWSKI FUNCTIONALS ON POLARIZATION FIELDS

A spin $s$ random field $P_s$ is a complex-valued field that satisfies the following relation at any fixed point $x$:

$$P'_s(x) = \exp(is\theta)P_s(x),$$

where $P'(x)$ represents the value in $x \in \mathbb{S}^2$ when the local coordinates are rotated by an angle $\theta$. This notion can be made mathematically rigorous viewing spin fields as a section of a so-called spin fibre bundle, first introduced in mathematical physics by Newman & Penrose (1966).

Under this definition, the CMB polarization is a Gaussian complex-valued spin-2 random field, which can be written as

$$P(x) = Q(x) + iU(x), \quad P'(x) = \exp(i2\theta)P(x),$$

where $Q$ and $U$ are real maps. It is important to note that, for any fixed choice of a coordinate system, polarization fields can be viewed as complex-valued Gaussian random variables at a single point, but they can no longer be considered as the realization of an isotropic map. Indeed, the correlation function at any two points is not invariant to rotations, and hence it does not only depend on the geodesic distance between the two points (as required by isotropy). Only in the special case where $s = 0$, we are back to the case of a complex-valued isotropic Gaussian map, which can always be viewed as the complexification of two Gaussian, independent, and identically distributed real-valued fields $Q$ and $U$.

The behaviour of the expected values for Lipschitz-Killing curvatures and other topological functionals of spin fields has been very recently investigated in the mathematical literature by Lerario et al. (2022). The main finding of these authors is that for 'band-limited' spin-valued random fields with power spectrum concentrated around multipoles equal or larger than $\ell_*$, the effect of the spin vanishes as $\ell_*$ grows. In that case, these geometric functionals have the same behaviour for spin Gaussian fields and for complex-valued scalar Gaussian fields. More precisely, the remainder terms (neglected when approximating the geometry of a spin $s \neq 0$ field with the spin $s = 0$ case) are of the form $\propto \frac{s}{\ell(t+1)}$. Thus, spin effects can be safely ignored for the random fields of cosmological interest, as there is very little power at small multipoles. Incidentally, from the mathematical point of view, these results continue to hold for any fixed value of the spin parameter $s$, and even under other conditions such as the (non-physical) case of a spin parameter growing with $\ell_*$. We can informally summarize this mathematical result as follows:

Up to negligible terms that depend only on the smallest multipoles, the expected values of the Lipschitz-Killing curvatures for the excursion sets of the norms of spin $s$ random fields are identical, for all values of $s$. In other words, defining

$$|P_s(x)|^2 := Q_s^2(x) + U_s^2(x),$$

we have

$$\mathbb{E}[\mathcal{L}_j(A_u(|P_s(.)|^2, \mathbb{S}^2)] = \mathbb{E}[\mathcal{L}_j(A_u(|P_0(.)|^2, \mathbb{S}^2)],$$

for all values of $j$ and $s$.

We are now going to derive the expected values of Lipschitz-Killing curvatures using the GKF, in an analogous way to the previous section. It should be noted that, up to second-order terms, results equivalent to ours have been given before in the literature (see Naselsky & Novikov 1998; Kasak et al. 2021, and the references therein). Apart from putting these results into more solid grounds by means of the approximation that we just explained, we believe that our derivation by means of the GKF is clearer and more amenable to generalizations.

In particular, assume that we have $Q$ and $U$ as two independent copies of Gaussian random fields with the same angular power spectrum, expected value 0, and variance normalized to 1 (i.e. $\mathbb{E}[Q^2(x)] = \mathbb{E}[U^2(x)] = 1$). Defining

$$P^2(x) = Q^2(x) + U^2(x),$$

it is immediate to demonstrate that $P^2$ is a chi-squared with 2 degrees of freedom, $P^2 \sim \chi_2^2$, and hence it has marginal density

$$f_{P^2}(u) = \frac{1}{2} \exp \left( -\frac{u}{2} \right).$$

We can now consider the excursion sets

$$A_u(P^2, \mathbb{S}^2) = \{x \in \mathbb{S}^2 : P^2(x) \ge u\}.$$

We can now recall again the GKF, which in the case of chi-square fields reads

$$\mathbb{E}[\mathcal{L}_j(A_u(P^2, \mathbb{S}^2))] = \sum_{k=0}^{d-j} \begin{bmatrix} k+j \\ k \end{bmatrix} \rho_{k;P^2}(u) \mathcal{L}_{k+j}(\mathbb{S}^2) \mu^{k/2}, \quad (5)$$

where as usual we write $\mu$ for the derivative of the covariance function evaluated at the origin. It should be noted that the expression for the expected values takes exactly the same form as in equation (4), simply replacing the functions $\rho_{k,\cdot}$ with $\rho_{k;P^2}(\cdot)$ defined by Adler & Taylor (2007, theorem 15.10.1):

$$\rho_{j,P^2}(u) = \frac{u^{(k-j)/2} \exp^{-u/2}}{(2\pi)^{j/2} \Gamma(k/2) 2^{(k-2)/2}} \\ \times \sum_{l=0}^{\lfloor \frac{j-1}{2} \rfloor} \sum_{m=0}^{j-1-2l} \mathbb{1}_{\{k \ge j-m-2l\}} \begin{pmatrix} k-1 \\ j-1-m-2l \end{pmatrix} \\ \times \frac{(-1)^{j-1+m+l}}{m!l!2^l} u^{m+l}.$$

More explicitly, for $j=1$ we obtain

$$\rho_{1,P^2}(u) = \frac{u^{1/2} \exp^{-u/2}}{(2\pi)^{1/2}},$$

whereas for $j=2$ we have

$$\begin{align*} \rho_{2,P^2}(u) &= \frac{\exp^{-u/2}}{(2\pi)} \sum_{m=0}^{1} \binom{1}{1-m} (-1)^{1+m} u^m \\ &= -\frac{\exp^{-u/2}}{(2\pi)} + \frac{\exp^{-u/2} u}{(2\pi)} \\ &= \frac{(u-1) \exp(-u/2)}{2\pi}. \end{align*}$$
Plugging these expressions into equation (5), we obtain, for the boundary length

$$
\begin{align*}
\mathbb{E}[2\mathcal{L}_1(A_u(P^2(\cdot), \mathcal{S}^2))] &= 2 \sum_{k=0}^{d-j} \begin{bmatrix} k+j \\ k \end{bmatrix} \rho_{P^2}(u) \mathcal{L}_{k+j}(\mathcal{S}^2) \sqrt{\mu} \\
&= 2 \begin{bmatrix} 2 \\ 1 \end{bmatrix} \frac{u^{1/2} \exp^{-u/2}}{(2\pi)^{1/2}} \mathcal{L}_2(\mathcal{S}^2) \sqrt{\mu} \\
&= 4 \frac{\omega_2}{\omega_1^2} \frac{u^{1/2} \exp^{-u/2}}{(2\pi)^{1/2}} 4\pi \sqrt{\mu} \\
&= 4 \frac{\pi}{4} \frac{u^{1/2} \exp^{-u/2}}{(2\pi)^{1/2}} 4\pi \sqrt{\mu} \\
&= \frac{\sqrt{2\pi}}{2} u^{1/2} \exp^{-u/2} 4\pi \sqrt{\mu}.
\end{align*}
$$

The excursion area can be immediately computed as

$$
\begin{align*}
\mathbb{E}[\mathcal{L}_2(A_u(P^2(\cdot), \mathcal{S}^2))] &= \mathbb{E} \left[ \int_{\mathcal{S}^2} \mathbb{P}(P^2(x) \ge u) dx \right] \\
&= \int_{\mathcal{S}^2} \mathbb{E} \left[ \mathbb{P}(P^2(x) \ge u) \right] dx \\
&= 4\pi \exp(-u/2).
\end{align*}
$$

Finally, for the Euler-Poincaré characteristic

$$
\mathbb{E}[\mathcal{L}_0(A_u(P^2, \mathcal{S}^2))] = 2\mu(u-1)\exp(-u/2) \\
\qquad + 2\exp\left(-\frac{u}{2}\right).
$$

Passing from the Lipschitz-Killing curvatures to MFs and normalizing by the area of the sphere, we finally obtain the following predictions for a chi-squared map with 2 degrees of freedom, such as $P^2$:

$$
\frac{\mathbb{E}[V_0(A_u)]}{4\pi} = \exp\left(-\frac{u}{2}\right), \qquad (6a)
$$

$$
\frac{\mathbb{E}[V_1(A_u)]}{4\pi} = \frac{\sqrt{2\pi}}{8} \sqrt{\mu u} \exp\left(-\frac{u}{2}\right), \qquad (6b)
$$

$$
\frac{\mathbb{E}[V_2(A_u)]}{4\pi} = \mu \frac{(u-1)\exp(-u/2)}{2\pi}, \qquad (6c)
$$

where the threshold is non-negative, $u \ge 0$. For negative thresholds, it follows from the definition that $V_0$ is exactly 1, while $V_1$ and $V_2$ are exactly 0. In particular, $V_2$ presents a discontinuity at $u=0$ connected with the number of non-polarized points (see Kasak et al. 2021).

In all these expressions, the variance of the field is unity and $\mu$ is the derivative of the covariance of the field at the origin, which can be computed as before:

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

where

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

with $C_{\ell}^{\mathrm{EE}}$ and $C_{\ell}^{\mathrm{BB}}$, respectively, the EE and BB angular power spectra computed from $Q$ and $U$ maps (Kamionkowski, Kosowsky & Stebbins 1997; Zaldarriaga & Seljak 1997).

4 PYNKOWSKI

We develop a PYTHON package called PYNKOWSKI to compute MFs on spherical maps formatted in the HEALPIX convention (Górski et al. 2005). We make this package publicly available to the research community in https://github.com/javicarron/pynkowski.

PYNKOWSKI can be used in two different ways:

(i) To compute the expected values of the MFs for Gaussian scalar fields, like CMB temperature anisotropies $T$ (see Section 2), or for $\chi^2$ fields, like CMB polarized intensity $P^2$ (see Section 3). They can be computed given either an input angular power spectrum or $\mu$.

(ii) To determine the MFs on any scalar map, either on the full sky or within any mask.

Both $T$ and $P^2$ are scalar fields on the sphere, and therefore the MFs can be computed with the same code for both quantities, although their expected results are different. We describe the implementation in this section, where we refer generally to an arbitrary (smooth) field on the sphere $f(x)$. We compute the MFs normalized by the area of the sphere ($v_i = V_i / 4\pi$), so the extension to masked maps is straightforward.

Without loss of generality, for this work, we normalize input maps to have unit variance, as we did in the previous sections when computing the theoretical predictions. In the case of $P^2$, both the $U$ and $Q$ maps are separately normalized to unit variance. The MFs can now be expressed as a function of the $a$-dimensional threshold $u$.

4.1 Computation of the MFs

4.1.1 First MF, $v_0$

The first MF, $v_0$, is computed according to the equation

$$
v_0(u) = \frac{1}{4\pi} \int_{\mathcal{S}^2} \Theta(f(x) - u) \, dx, \qquad (7)
$$

where $\Theta(r)$ is the Heaviside function (1 where $r \ge 0$; 0 otherwise). In the computation, we approximate the integral with a sum over all valid pixels. The integral divided by the area of the sphere, $4\pi$, can be seen as the average value of the argument inside the integral. Therefore, it can be approximated with the average over all pixels. Thus, $v_0(u)$ is computed as the fraction of the map above the threshold $u$.

4.1.2 Second MF, $v_1$

The second MF, $v_1$, is computed according to the equation

$$
\begin{align}
v_1(u) &= \frac{1}{4\pi} \frac{1}{4} \int_{\partial A_u} dr \nonumber \\
&= \frac{1}{4\pi} \frac{1}{4} \int_{\mathcal{S}^2} \delta(f(x) - u) |\nabla f| \, dx, \tag{8}
\end{align}
$$

where the second equality is the result of a change of coordinates from the line element $dr$ of the boundary of the excursion set to a surface element $dx$ on the sphere. The gradient of $f$ is denoted by $\nabla f$ and $\delta$ is the Dirac delta.

We start by computing $|\nabla f|$ at every pixel. Then, in order to approximate the integral with a sum over pixels, we provide a bin of $\Delta u$. We create a mask $M_u(x)$ identifying the pixels where the value of the map is close to $u$:

$$
M_u(x) = \begin{cases} 1 & \text{if } u - \frac{\Delta u}{2} \le f(x) \le u + \frac{\Delta u}{2}, \\ 0 & \text{otherwise.} \end{cases} \qquad (9)
$$

Now, $v_1(u)$ can be computed as the mean of the pixels of $M_u \cdot |\nabla f|$, divided over $\Delta u$, the bin size, and 4, the normalization factor.
### 4.1.3 Third MF, $v_2$

The third MF, $v_2$, is computed according to the equation

$$v_2(u) = \frac{1}{4\pi} \frac{1}{2\pi} \int_{\partial A_u} \kappa(r) dr \\ = \frac{1}{4\pi} \frac{1}{2\pi} \int_{S^2} \delta(f(x) - u) |\nabla f| \kappa \, dx. \qquad (10)$$

The change of coordinates to surface element and the approximation to a finite sum of pixels is analogous to the previous case. We define the mask $M_u$ identically and compute the mean of the pixels of $M_u$. $|\nabla f| \cdot \kappa$. This quantity is then divided by the bin size $\Delta u$ and by $2\pi$ to yield the result for $v_2(u)$.

In this formula, $\kappa$ is the geodesic curvature of the boundary of the excursion set. Given that this curve is given by the implicit equation $f(x) = u$, we can use the expression of $\kappa$ for an implicit curve (see e.g. do Carmo 1976):

$$\kappa = \frac{2 f_{;\phi} f_{;\theta} f_{;\phi\theta} - f_{;\theta}^2 f_{;\phi\phi} - f_{;\phi}^2 f_{;\theta\theta}}{(f_{;\phi}^2 + f_{;\theta}^2)^{\frac{3}{2}}}, \qquad (11)$$

where the semicolon denotes the covariant derivative. All the quantities involved in this expression can be computed for each pixel.

## 4.2 Computation of the derivatives

The spatial derivatives of the maps (both partial and covariant) are computed with the help of the HEALPY<fnref n="3" /> function alm2map_der in harmonic space (Zonca et al. 2019). The covariant derivatives of a field $f$ on the sphere are computed in terms of the partial derivatives:

$$\begin{align*} f_{;\theta} &= f_{,\theta}, \\ f_{;\phi} &= \frac{f_{,\phi}}{\cos \theta}, \\ f_{;\theta\theta} &= f_{,\theta\theta}, \\ f_{;\phi\phi} &= \frac{f_{,\phi\phi}}{\cos \theta} - \frac{\sin \theta f_{,\phi}}{\cos^2 \theta}, \\ f_{;\phi\phi} &= -\frac{f_{,\phi\phi}}{\cos^2 \theta} + f_{,\theta} \tan \theta, \end{align*}$$

where we consider $\theta \in (-\pi/2, \pi/2)$ as the latitude and $\phi \in [0, 2\pi)$ as the longitude.

The need of estimating the spatial derivatives of the map means that the computations are reliable only for maps whose angular power spectrum is negligible at multipoles close to $\ell_{max} = 3 \cdot N_{side} - 1$; i.e. the derivatives can only be computed if the maps are reasonably smooth. This becomes a requirement in order to estimate $v_1$ and $v_2$, while $v_0$ can be computed for any map as it does not require the computation of the derivatives.

This smoothness condition is satisfied for realistic $T$ and $P^2$ maps, as they are reasonably smooth, but it can break if one tries to analyse maps of pure noise or heavily dominated by the modes at the highest multipoles (corresponding to pixel-size scales). We note that our theoretical results only hold for functions that are at least twice differentiable almost everywhere (for other cases, such as a fractal behaviour, see Lan, Marinucci & Xiao 2018).

# 5 SIMULATIONS AND DATA

## 5.1 Gaussian simulations

In order to validate the implementation of the computation on HEALPIX maps and the theoretical predictions, we generate 300 Gaussian isotropic CMB maps of the Stokes parameters ($Q, U$). We use the HEALPY function synfast with the best-fitting angular power spectrum for polarization (E and B modes) as reported by Planck Collaboration V (2020b). For these simulations, we fix a resolution of $N_{side} = 1024$; this is half of the Planck value in order to reduce the computational time without a significant loss of resolution; we have verified that the results remain unchanged by this choice.

For the estimation of $v_0$, we exclude the modes with $\ell \le 4$ in the simulated maps in order to avoid spin effects in the computation of the MF (see Section 3). Conversely, when $v_1$ and $v_2$ are estimated, a smoothing of the maps with a Gaussian beam with full width at half maximum (FWHM) equal to 15 arcmin is applied in order not to be affected by pixelization effects in the computation of the spatial derivatives.

MFs have been computed over the full-sky maps of $P^2 = Q^2 + U^2$, where $Q$ and $U$ have been previously normalized to unit variance. We analyse the deviation of each MF at different values of the threshold $u$ in these simulated realizations of Gaussian isotropic CMB with respect to the theoretical expectations (see Section 6.1). In this way, we are able to simultaneously assess the accuracy of the formulae obtained in Section 3 and verify the implementation of MF computations on simulated maps.

## 5.2 Planck polarization maps

After validating the machinery with Gaussian isotropic CMB simulations, we apply the MFs to Planck CMB polarization data. In order to verify that the results are consistent and do not depend on the component separation algorithm used to reconstruct the CMB, we use the polarization maps released by the Planck Collaboration obtained with two different pipelines: Spectral Matching Independent Component Analysis (SMICA) and Spectral Estimation Via Expectation Maximization (SEVEM), as reported in Planck Collaboration IV (2020a).

SMICA (Delabrouille, Cardoso & Patanchon 2003) performs a foreground and noise subtraction by linearly combining Planck frequency maps in harmonic space with multipole-dependent weights, up to $\ell = 4000$, which minimize on each angular scale the variance of the output map. In polarization, it combines the E and B modes independently and then recombines them to obtain the $Q$ and $U$ maps.

SEVEM (Martínez-González et al. 2003), instead, is a template-fitting method. Different templates of the foreground emission are estimated as the difference of maps from neighbouring frequencies. This procedure removes the CMB signal while preserving the Galactic one. A linear combination of these templates is then subtracted from the input data to produce a cleaned CMB map at a specific Planck frequency. The coefficients of the linear combination are determined by minimizing the variance of the cleaned map outside a mask, which covers the point sources detected in polarization and the 3 per cent brightest Galactic emission. The 100, 143, and 217 GHz cleaned maps are then combined in harmonic space, using E and B decomposition, to produce the final CMB maps for the $Q$ and $U$ components at a resolution of FWHM = 5 arcmin and a maximum considered multipole of $\ell = 3000$.

<footnote n="3">³https://github.com/healpy/healpy</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** 2018 Component Separation Common Planck mask for polarization. It excludes the most contaminated regions near the Galactic plane and the polarized point sources. It keeps a total sky fraction of $f_{sky} = 0.78$.</figcaption></figure>

The native resolution of these released maps<fnref n="4" /> is $N_{side} = 2048$; we degrade them to $N_{side} = 1024$ to match exactly the pipeline applied for simulations. Similarly, before the computation of $V_0$ we again exclude the lowest multipoles ($\ell \le 4$) in order to avoid spin effects, while for the computation of $V_1$ and $V_2$ we bring the maps to a resolution of FWHM = 15 arcmin to not be affected by pixelization effects.

Planck CMB polarization maps still present residual contamination due to Galactic foregrounds and point sources. The most contaminated regions of the sky must be excluded from any cosmological analysis. Therefore, we adopt the 2018 Component Separation Common polarization mask (shown in Fig. 1) throughout the analysis, as recommended by Planck Collaboration IV (2020a). This mask was obtained by thresholding the standard deviation between each of the four cleaned Planck CMB maps and it was then augmented with the Commander and SEVEM confidence masks, as well as with the SEVEM and SMICA in-painting masks. The masked SMICA and SEVEM P maps are shown in Fig. 2.

## 5.3 Planck end-to-end simulations

The SMICA and SEVEM maps include residual effects associated with Planck instrumental noise, which are especially relevant in polarization. Due to the scanning strategy of the satellite, this component presents an anisotropic pattern. Such contamination has to be properly taken into account for any cosmological analysis of the Planck CMB maps. To this end, the Planck Collaboration released end-to-end simulations,<fnref n="4" /> which include the expected CMB signal and noise residuals.

These simulations have been obtained by generating realistic CMB (Gaussian and isotropic) and noise maps at the different Planck frequency channels, which are then processed by each component separation algorithm. The resulting maps are thus able to reflect the expected statistical properties across the sky of noise residuals and CMB given by each method.

We compare the MFs of the SMICA- and SEVEM-released CMB polarization maps with those computed on these realistic simulations. Compatible results would suggest that any non-Gaussianity or deviation from statistical isotropy present in CMB polarization maps is compatible with the ones produced by noise residuals. A significant deviation might signal the presence of residual foregrounds, unmodelled systematic effects, or a hint of non-Gaussianity or departure from statistical isotropy of the CMB with primordial origin.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Planck CMB maps for $P = \sqrt{Q^2 + U^2}$ obtained with the SMICA (top) and SEVEM (bottom) pipelines. They have been smoothed with a Gaussian beam of FWHM = 30 arcmin for visualization purposes.</figcaption></figure>

We use both SMICA and SEVEM simulations with $N_{side} = 1024$ to make all the analyses identical in all maps. Similarly to real Planck polarization maps, we again exclude the lowest multipoles ($\ell \le 4$) before the computation of $V_0$, while we bring the maps to a resolution of FWHM = 15 arcmin for the computation of $V_1$ and $V_2$.

## 5.4 Needlet maps

The comparison of MFs computed on Planck $P^2$ maps with those of end-to-end simulations allows us to perform a statistical analysis over the full range of angular scales observed by the Planck satellite. However, any non-Gaussian or statistically anisotropic signal may be relevant only in a specific range of multipoles (e.g. Galactic residuals on large scales or point sources at high multipoles). Therefore, assessing the compatibility between maps and end-to-end simulations separately at different angular scales may provide a more detailed statistical inspection of the data. We thus compare MFs of Planck SMICA and SEVEM $P^2$ maps with those of realistic end-to-end simulations (introduced, respectively, in Sections 5.2 and 5.3) in needlet domain. A similar analysis for E-mode Planck maps has been reported in Planck Collaboration VII (2020d).

Needlets are a wavelet system that guarantees simultaneous localization in harmonic and pixel space; they have been introduced in the statistical literature by Narcowich, Petrushev & Ward (2006) and first applied to CMB data by Pietrobon, Balbi & Marinucci (2006). Given any field of spin $s$ defined on the sphere, the corresponding needlet maps $\beta_j^s$ are obtained by filtering its harmonic coefficients $a_{\ell m}^s$ with a harmonic weighting function $b_j(\ell)$, which isolates modes at different angular scales for each needlet scale $j$ (Geller & Marinucci 2010):

$$ \beta_j^s(x) = \sum_{\ell,m} (a_{\ell m}^s \cdot b_j(\ell)) s Y_{\ell m}(x), \qquad (12) $$

<footnote n="4">⁴https://pla.esac.esa.int/pla/</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3. Representation in the harmonic domain of the needlet configuration adopted in this analysis: standard needlets with $B = 2$. The first four bands have been merged together according to equation (14), while the other bands are left unchanged.</figcaption></figure>

where $x$ is a direction in the sky and $s_{Y_{\ell m}}$ are the spin-weighted spherical harmonics. Such procedure in harmonic space is equivalent to performing a convolution of the map in real domain. The shape of the needlet bands is defined by the choice of the harmonic function $b$ whose width is set by the parameter $B$: lower values of $B$ correspond to a tighter localization in harmonic space (less multipoles entering into any needlet coefficient), whereas larger values ensure wider harmonic bands.

As already mentioned in Section 3, the Stokes parameters $Q$ and $U$ are components of a spin-2 field on the sphere:

$$Q \pm iU(x) = \sum_{\ell,m} a_{\ell m}^{\pm 2} \cdot \pm 2 Y_{\ell m}(x). \qquad (13)$$

Therefore, in this paper, we first filter the $a_{\ell m}^{\pm 2}$ coefficients to obtain $Q$ and $U$ needlet maps and then compute the corresponding $P^2$ map. As needlet harmonic function $b$, we have adopted the standard needlet construction (Narcowich et al. 2006) with $B = 2$. The needlet filters have been obtained through the PYTHON package MTNEEDLETS<fnref n="5" /> (Carrón Duque et al. 2019).

Since the first needlet bands would select only very few large-scale modes, we merge together the first four bands as follows:

$$b_0^{\text{new}}(\ell) = \sqrt{\sum_{j=0}^{3} b_j^2(\ell)}. \qquad (14)$$

The adopted configuration of needlet harmonic filters is shown in Fig. 3. We can note that the $j$th needlet band is non-zero only for $\ell \le 2^{j+4}$. Therefore, needlet maps are band-limited and it suffices to reconstruct them at $N_{\text{side}} = 2^{j+3}$. As it was done in Planck Collaboration VII (2020d), at each needlet scale, corresponding masks are obtained by first downgrading the mask at the targeted resolution as an intensity map. The resulting smoothed downgraded mask is then thresholded by setting pixels with values lower than 0.9 to zero and all others to unity so that we again have a binary mask.

## 5.5 LiteBIRD and SO simulations

CMB intensity anisotropies have been measured by Planck with very high accuracy (Planck Collaboration V 2020b). Conversely,

polarization data are mostly noise-dominated and contaminated by several instrumental systematic effects. As already discussed in Section 1, several experiments, either from the ground or from space, are planned to measure CMB polarization with unprecedented sensitivity. Therefore, in this paper, we forecast the capability of MFs of the polarization modulus to detect a deviation from Gaussianity and statistical isotropy in realistic simulated maps for two future and complementary experiments: SO (Ade et al. 2019) and the LiteBIRD satellite (LiteBIRD Collaboration 2023). The SO is composed of a Large Aperture Telescope (SO-LAT) and several refracting SO-SATs. SO-LAT will observe 40 per cent of the sky and will be devoted to measurements of the CMB small-scale temperature and polarization anisotropies, the CMB lensing, the Sunyaev-Zel'dovich effects, and extragalactic sources. SO-SATs, instead, aim at a detection of primordial B modes. Given the affinity with the scientific target of LiteBIRD, in this analysis, we have considered only SO-SATs.

For both experiments, we simulate a fiducial Gaussian data set and a test data set that also includes a non-Gaussian anisotropic signal. In the fiducial Gaussian data set, $Q$ and $U$ maps are obtained by adding Gaussian and isotropic CMB and noise realizations:

$$\{Q, U\}_{\text{fid}} = \{Q, U\}_{\text{CMB}} + \{Q, U\}_{\text{noise}}. \qquad (15)$$

In the test data set, instead, we also add a map of simulated Galactic foreground emission as an example of a generic non-Gaussian and anisotropic component:

$$\{Q, U\}_{\text{test}} = \{Q, U\}_{\text{CMB}} + \{Q, U\}_{\text{noise}} + \alpha_f \{Q, U\}_{\text{fgd}}, \qquad (16)$$

where $\alpha_f$ is a scaling factor that modulates the overall amplitude of this non-Gaussian contamination. Such rescaling roughly parametrizes the level of residual Galactic foregrounds after a component separation analysis has been performed on CMB multifrequency data. We let $\alpha_f$ to vary in order to test the sensitivity of $P^2$ MFs to different magnitudes of the simulated non-Gaussian and anisotropic signal across the sky.

The CMB maps are generated with the HEALPIX PYTHON package (Zonca et al. 2019) from the angular power spectra of the best-fitting Planck 2018 parameters (Planck Collaboration VI 2020c). For the instrumental noise component, we first simulate noise maps at the different frequency channels of the considered experiment as Gaussian realizations from the following power spectrum:

$$N_{\ell} = \left( \sigma_{Q,U} \frac{\pi}{180 \ 60} \right)^2 \left[ 1 + \left( \frac{\ell}{\ell_{\text{knee}}} \right)^{\alpha_{\text{knee}}} \right], \qquad (17)$$

where $\sigma_{Q,U}$ is the polarization sensitivity in $\mu K_{\text{CMB}} \cdot \text{arcmin}$, while $\ell_{\text{knee}}$ and $\alpha_{\text{knee}}$ refer to the contribution from $1/f$ noise. The employed LiteBIRD and SO-SAT specifications are taken from the latest forecasts of the two experiments (Ade et al. 2019; LiteBIRD Collaboration 2023). For LiteBIRD, we do not consider the $1/f$ term whose contribution is expected to be fully suppressed by the modulation of the half wave plate. For SO-SATs, we have adopted the values for $\ell_{\text{knee}}$ and $\alpha_{\text{knee}}$ reported in Ade et al. (2019) for the pessimistic $1/f$ case. In this way, we will obtain the most conservative constraints on the detection power of $P^2$ MFs for SO-SATs, which is likely to improve for more realistic instrumental configurations. For both experiments, we neglect for simplicity the scanning strategy of the telescopes and, thus, the statistical properties of the noise are assumed constant across the sky. This is a fair assumption for LiteBIRD, since, as also discussed in LiteBIRD Collaboration (2023), the scanning strategy parameters are chosen so as to have the distribution of hit counts as uniform as possible and are expected to have a completely negligible effect in $P^2$ maps. For SO-SATs, this statement does not hold, because the map of hit counts is very

<footnote n="5">⁵https://javicarron.github.io/mtneedlet</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** In light grey: The SO-SAT full footprint with sky fraction $f_{sky} = 0.34$. In dark grey: The SO-SAT patch considered in this analysis, which corresponds to the sky area with the highest signal-to-noise ratio in the full footprint. This was obtained considering pixels with the largest values in the map of hit counts for a final sky fraction $f_{sky} = 0.1$.</figcaption></figure>

inhomogeneous. Therefore, we restrict the MF computation to a smaller patch with $f_{sky} = 0.1$ (shown in Fig. 4), where the scanning strategy can be confidently assumed uniform. Such patch is obtained by considering only the 10 per cent of the pixels with the highest values in the map of hit counts provided by the SO Collaboration. The adopted SO-SAT sensitivities, which are reported in Ade et al. (2019), indeed refer to a homogeneous map of hit counts with such a sky fraction. The final noise map, for both experiments, is obtained by linearly combining noise maps of the central frequency channels with an inverse noise weighting. We considered the frequency range $93 \le \nu \le 225$ for SO-SATs and $89 \le \nu \le 235$ for LiteBIRD.

In this paper, we aim at providing a first simplified forecast of the capability of the MFs to detect non-Gaussian and anisotropic signals in $P^2$ maps for experiments with much higher sensitivity with respect to Planck. A more detailed analysis would require considering noise residuals given by component separation pipelines. However, we do not want to add the additional degree of freedom of the choice of specific component separation methods to be applied to the considered experiments and we leave such analysis for future work.

The non-Gaussian and anisotropic signal ${Q, U}_{fgd}$ injected in the test data set of equation (16) is a simulated polarized Galactic foreground emission at 100 GHz multiplied by a varying rescaling factor $\alpha_f$. This signal is generated with the PYSM PYTHON package<fnref n="6" /> (Thorne et al. 2017) and includes the two main polarized Galactic components: synchrotron and thermal dust emission. They are, respectively, modelled with the PYSM s5 and d10. The choice of these models is driven by the fact that non-Gaussianity is injected in them up to much smaller scales than in other analogous commonly used PYSM models, such as s1 and d1. The synchrotron spectral energy distribution (SED) follows a power law:

$$ \{Q, U\}_{s}(x, \nu) = \{Q, U\}_{s}(x, \nu_0) \left(\frac{\nu}{\nu_0}\right)^{\beta_s(x)}, \qquad (18) $$

where $\beta_s$ is the spectral index, $x$ the position in the sky, and $\nu$ the considered frequency. ${Q, U}_{s}(x, \nu_0)$ represents the synchrotron template at a reference frequency $\nu_0$ and is the Wilkinson Microwave Anisotropy Probe 9-yr 23-GHz Q/U map (Bennett et al. 2013) smoothed with a Gaussian kernel of FWHM = 5° with non-Gaussian small scales added through the Logarithm of the Polarization Fraction

Tensor (logpoltens) formalism. The thermal dust SED is assumed to be a modified blackbody:

$$ \{Q, U\}_{d}(x, \nu) = \{Q, U\}_{d}(x, \nu_0) \left(\frac{\nu}{\nu_0}\right)^{\beta_d(x)} \frac{B_{\nu}(T_d(x))}{B_{\nu_0}(T_d(x))}, \qquad (19) $$

where $B_{\nu}(T)$ is the blackbody spectrum, $\beta_d$ is the dust spectral index, $T_d$ is the dust temperature, and ${Q, U}_{d}(x, \nu_0)$ represents the dust template at a reference frequency $\nu_0$. Both the maps of the spectral parameters and the template are obtained from the application of the Generalized Needlet Internal Linear Combination to Planck data (Planck Collaboration IV 2020a) and non-Gaussian small-scale fluctuations have been added in the logpoltens formalism.

We note that the injected non-Gaussian anisotropic signal is particularly relevant on large and intermediate angular scales. We thus expect to be able to constrain effectively the detection power of the proposed tools only in that range of multipoles that, however, is exactly that targeted by LiteBIRD and SO-SATs for accurate measurements of the polarization anisotropies. Moreover, small scales are usually noise-dominated. Therefore, Galactic emission represents a good probe to forecast the MF detection power for these future CMB experiments. If one wants to expand the assessment of MF sensitivity at smaller scales, one would have to choose a signal model different from Galactic emission with non-Gaussianities and anisotropies at those scales.

The capability of MFs to detect deviations from Gaussianity and statistical isotropy in the test data set for SO-SATs and LiteBIRD is compared with that of Planck. We thus generated also for Planck a fiducial and a test data set, where CMB and Galactic foreground maps are the same of LiteBIRD and SO-SATs, while, as noise component, we used the end-to-end SMICA noise residual maps introduced in Section 5.3. All maps are brought to the resolution of 37.8 arcmin for the comparison between Planck and LiteBIRD, of 30 arcmin for the comparison between Planck and SO-SATs.

## 5.6 Compatibility analysis

In order to assess the level of deviation between MFs on different maps, we use the $\chi^2$ statistic. Let $X$ be a vector with the result of any of the MFs with the mean of the simulations subtracted and $\Sigma$ the covariance matrix of $X$ (representing the correlation of MFs at different thresholds). The $\chi^2$ statistic is given by

$$ \chi^2 = X^T \Sigma^{-1} X. \qquad (20) $$

If $S$ is the sample covariance (computed with simulations), an unbiased estimator of the inverse of the covariance is

$$ \Sigma^{-1} = \frac{N - d - 2}{N - 1} S^{-1}, \qquad (21) $$

where $N = 300$ is the number of simulations and $d = 28$ is the number of bins in the considered range of thresholds $u$. The pre-factor in this equation is known as the Hartlap factor (Hartlap, Simon & Schneider 2007) and it is needed for the estimator to be unbiased.

In this analysis, we report the reduced $\hat{\chi}^2 = \chi^2/N_{dof}$, where $N_{dof} = 28$ is the number of degrees of freedom (i.e. the number of thresholds in our computation). We also compute the probability-to-exceed a given $\chi^2$:

$$ p_{\text{exc}}(\chi_{\text{data}}^2) = P(\chi_{\text{sims}}^2 > \chi_{\text{data}}^2), \qquad (22) $$

where $\chi_{\text{data}}^2$ and $\chi_{\text{sims}}^2$ are, respectively, associated with any MF computed on the data or on a simulation. In equation (22), a value close to 50 per cent means that the data closely resemble the

<footnote n="6">⁶https://github.com/galsci/pysm</footnote>
expectation, while a value close to 0 or 100 per cent means that the data are highly incompatible with simulations.

Finally, we also report the significance of the deviation $\Delta_{\chi}$ of $\chi^2$ in the data with respect to the mean $\bar{\chi}_{\text{sims}}^2$ among the simulations, reported in units of standard deviation $\sigma_{\text{sims}}$:

$$ \Delta_{\chi} = \frac{\chi_{\text{data}}^2 - \bar{\chi}_{\text{sims}}^2}{\sigma_{\text{sims}}}. \qquad (23) $$

To assess the $P^2$ MF detection power of a non-Gaussian anisotropic signal in the simulated data sets introduced in Section 5.5, we compute a $\chi^2$ for each realization as in equation (20), where $X$ is the MF difference between the test and the fiducial Gaussian data set for that realization and the covariance matrix is estimated from MFs of the fiducial data set. We then report for each case the average probability for a $\chi^2$ distribution with $N_{\text{dof}} = 28$ (the adopted number of thresholds) to exceed the computed $\chi^2$-s. If this value is close to 50 per cent, it means that no deviations are observed between MFs of the test and Gaussian fiducial data sets. Conversely, if such probability is very close to 0 per cent, we are able to detect at high significance the injected non-Gaussian anisotropic signal in the test data set.

# 6 RESULTS

In this section, we proceed to validate the code and the theoretical formulae for the expected values of the MFs for $P^2$ maps. In order to assess the capability of these tools to detect deviations from Gaussianity and statistical isotropy, we then apply the pipeline to Planck CMB polarization data and simulated $P^2$ maps of LiteBIRD and SO. We use 300 simulations for each case.

Throughout all the following analysis, we set thresholds between $u = 0$ and 5. For visualization purposes, all MFs are plotted in 100 bins, which correspond to a spacing of $\Delta u = 0.05$. The statistical analysis described in Section 5.6 is instead performed lowering the number of bins to 28 ($\Delta u = 0.18$) so as to get a more robust estimate of the covariance matrix given in equation (21).

In Section 6.1, we compare the theoretical expectation of MFs for a $P^2 = Q^2 + U^2$ field under the assumption of Gaussianity and isotropy with those computed on CMB simulations generated with the same statistical properties. In Section 6.2, MFs are estimated on the Planck SMICA and SEVEM CMB polarization modulus map and compared with those of the corresponding end-to-end simulations. In Section 6.3, we assess the detection power of a non-Gaussian and anisotropic signal at different angular scales for MFs computed on realistic simulated $P^2$ data sets of LiteBIRD and SO and the results are compared with those obtained from analogous Planck simulations.

## 6.1 Validation

We verify the theoretical predictions for the MFs of $P^2$ introduced in Section 3 by comparing them with the results of simulated CMB Gaussian isotropic maps. We simulate 300 such maps as described in Section 5.1 and compute the MFs on them as explained in Section 4. Maps are generated with Planck best-fitting angular power spectra, so $\mu$ can be exactly computed. We use this value to obtain the theoretical expectations of the MFs according to equation (6). 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.

The results of the MFs for these simulated maps can be seen in Fig. 5, together with the theoretical predictions. We find that the simulations are fully compatible with the theoretical expectations for all three MFs, with differences between the average of the simulations

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Average MFs of $P^2 = Q^2 + U^2$ for $N_{\text{sims}} = 300$ CMB Gaussian isotropic simulations (dots) are compared with the theoretical predictions (solid lines). From top to bottom: $V_0$, $V_1$, and $V_2$. The values of the MFs are shown in the top of each panel. The difference between average MFs from simulations and theory is compared with the dispersion among simulations ($1\sigma$, given by the grey contours) and the dispersion of the mean of simulations (computed as $\sigma/N_{\text{sims}}$) in the middle and bottom sections of each panel, respectively. See text for details.</figcaption></figure>
and the theory well inside the $1\sigma$ region, where $\sigma$ is the dispersion of the MFs among the 300 simulations.

Furthermore, in order to explore any possible systematic deviation, we compare the average residuals in the simulations with respect to the theoretical expectations considering the standard deviation of the mean ($\sigma_{\text{mean}} = \sigma_{\text{sims}}/\sqrt{300}$). This comparison is shown in the bottom section of each panel in Fig. 5. The three MFs are all perfectly compatible with theory, both in the case of individual simulations and of the average behaviour.

We also note the low statistical variation of these curves: the relative uncertainty $\sigma(V_i(u))/V_i(u)$ at every threshold $u$ is below 1 per cent and, in the case of $V_0$, it is below 1 part in 1000. 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, $l_{\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 in Fantaye et al. (2015). Such low levels of fluctuations in the values of MFs make them suitable tools to detect deviations from Gaussianity and isotropy in real data.

These results provide a double validation: on the one hand, they verify the mathematical framework used to predict the expected values of the MFs of $P^2$ maps; on the other hand, they validate our implementation in the PYNKOWSKI package. We have verified that masking the maps has a negligible effect on these results beyond slightly increasing the noise due to the smaller sky fraction. This robustness is expected because, unlike the angular power spectrum, MFs are purely local quantities. Thus, the effect of masking is completely negligible.

## 6.2 Analysis of Planck 2018 polarization maps

Once both the theoretical predictions and our implementation to compute MFs have been validated, we can apply this formalism to the observed CMB polarization maps in order to assess the presence of any non-Gaussianity or deviation from statistical isotropy. In this work, as a proof of concept, we analyse the CMB maps reported by Planck (Planck Collaboration IV 2020a), produced with two different methods: SMICA (shown in this section) and SEVEM (shown in Appendix A) (see Section 5 for details). We verify that our results are consistent on both complementary pipelines.

Planck polarization maps are affected by two known contaminating signals that would bias our MF estimates if they are unaccounted for. The first is the residual emission from the Galactic foregrounds and point sources. In order to minimize their impact, we adopt the 2018 Component Separation Common polarization mask (Planck Collaboration IV 2020a), shown in Fig. 1. The second contaminant is the residual noise whose statistical properties are not homogeneous across the sky due to the scanning strategy of the Planck satellite (Planck Collaboration I 2011), as can be seen in Fig. 2. This effect is consistently included in the end-to-end simulations provided by the collaboration, as explained in Section 5.3. The impact of the anisotropic distribution of the statistical properties of noise is discussed in Appendix B.

In Fig. 6, we report the comparison between the MFs computed on the Planck polarization modulus map and on end-to-end simulations, both produced by SMICA. We see for all MFs a remarkable agreement between them with all residuals within $\pm 2\sigma$. Only $V_0$ appears to be consistently below the average values in simulations at high thresholds, but this deviation is small and consistent with a statistical fluctuation.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** MFs of the Planck SMICA-released $P^2 = Q^2 + U^2$ maps (dots) are compared with the mean given by realistic end-to-end simulations (solid lines). The first, second, and third MFs are shown (top to bottom), with the values (top section of each panel) and the residual with respect to the mean of simulations (bottom section of each panel). In all plots, the pink shaded area represents the $\pm 1\sigma$ region obtained in the simulations.</figcaption></figure>

We quantify the compatibility of all curves in terms of the $\chi^2$ statistic (as well as the probability-to-exceed and the deviation with respect to simulations in terms of $\Delta\chi$), as explained in Section 5.6. The results can be found in Table 1 and show that no deviation is statistically significant or consistently detected across different component separation algorithms.

However, we can highlight some hints of deviation. The residuals of $V_0$ seem to present a pattern at high values of the threshold for both
<table><caption>**Table 1.** Deviations of MFs computed on *Planck* SMICA and SEVEM CMB $P^2$ maps with respect to those of the corresponding end-to-end simulations, in terms of normalized $\chi^2$, probability-to-exceed, and $\Delta_\chi$ of equation (23).</caption><thead><tr><th></th><th></th><th>χ<sup>2</sup></th><th>P<sub>exc</sub><br>(per cent)</th><th>Δ<sub>χ</sub></th></tr></thead><tbody><tr><td rowspan="2">V<sub>0</sub></td><td>SMICA</td><td>1.02</td><td>43.7</td><td>0.04</td></tr><tr><td>SEVEM</td><td>0.81</td><td>75.3</td><td>-0.71</td></tr><tr><td rowspan="2">V<sub>1</sub></td><td>SMICA</td><td>0.91</td><td>60.7</td><td>-0.31</td></tr><tr><td>SEVEM</td><td>1.36</td><td>8.3</td><td>1.48</td></tr><tr><td rowspan="2">V<sub>2</sub></td><td>SMICA</td><td>1.22</td><td>21.0</td><td>0.72</td></tr><tr><td>SEVEM</td><td>0.82</td><td>68.0</td><td>-0.60</td></tr></tbody></table>

SMICA and SEVEM solutions. From the $\chi^2$ analysis, this trend is not significant and can be explained by the fact that the covariance matrix of $V_0$ presents significant off-diagonal terms, unlike the other two functionals. We note that this pattern is still present when adopting a more aggressive Galactic mask (with $f_{sky} = 0.4$), thus disfavouring the hypothesis of a residual Galactic contamination origin.

So far, we have focused on *Planck* $P^2$ maps, which include the full range of multipoles. However, as already discussed in Section 5.4, applying an equivalent procedure in the needlet domain may provide a more detailed assessment separately at different angular scales. We apply a needlet transform on $Q \pm iU$ SMICA and SEVEM *Planck* maps and corresponding end-to-end simulations and then construct needlet $P^2$ maps. This is done by adopting the configuration of needlet bands shown in Fig. 3. We compute MFs of such obtained needlet maps and perform the same compatibility analysis of the full-range pixel-space case. The results for all needlet bands are shown in Table 2, where we report the probability-to-exceed to ease the comparison with the analogous analysis performed for E modes in Planck Collaboration VII (2020d).

As in the full-range pixel-space analysis, we do not observe statistically significant deviations from the MFs of the simulated CMB and noise components at any scale. For SEVEM, we had to mask a 3 per cent of the sky along the Galactic plane before applying the needlet transform to avoid leaking of a residual foreground component outside the *Planck* common mask in the first needlet scales. We used the publicly released *Planck* GAL97 Galactic mask.<fnref n="7" />

In conclusion, it is possible to state that, within the MF sensitivity in $P^2$, no unexpected deviation from Gaussianity and statistical isotropy has been detected in the *Planck* CMB $P^2$ maps. However, given the noise level in the *Planck* polarization observations, such an analysis allows us to assess mainly the statistical properties of noise rather than those of the CMB signal. This is expected to change in the near future with upcoming high-sensitivity CMB experiments, as we shall show in the following section.

## 6.3 Forecasts for LiteBIRD and SO

The *Planck* Collaboration provided cosmic variance limited measurements of the CMB temperature anisotropies up to very high multipoles (Planck Collaboration V 2020b). Conversely, in polarization, only modes at large angular scales have a modest signal-to-noise ratio. This is not the case for future surveys, which will provide much more accurate polarization data. In this section, we present a forecast of the sensitivity of $P^2$ MFs to non-Gaussian anisotropic signals for two CMB experiments: *LiteBIRD* and SO-SATs. One of the main challenges these experiments will need to face in their hunt

for primordial B modes will be an exquisite control of residual contamination by Galactic foregrounds, which is highly non-Gaussian and anisotropic. Therefore, we test the ability of MFs applied to $P^2$ maps to detect such a signal with varying amplitude in simulated *LiteBIRD* and SO-SAT data. The obtained results are compared with those of analogous *Planck* simulated polarization maps.

The *LiteBIRD* and SO simulated data sets have been described in Section 5.5. The analysis is performed in the needlet domain to assess the detection power of $P^2$ MFs at different angular scales. For each experiment, the MFs are computed on needlet maps of the fiducial and test data sets. The fiducial data set includes only Gaussian CMB and noise components, while in the test data set we injected a non-Gaussian and anisotropic diffused Galactic emission with a rescaled overall amplitude (see Section 5.5). Since such signal is particularly relevant on large and intermediate angular scales, we report results only for the first six needlet scales. As explained in Section 5.6, we report for each case the average probability for a $\chi^2$ distribution with $N_{dof} = 28$ (the adopted number of thresholds) to exceed the computed $\chi^2$s.

### 6.3.1 LiteBIRD versus Planck

*LiteBIRD* will be the successor of the *Planck* satellite with much improved sensitivity to CMB polarization. In this section, we report the improvement in detecting any deviation from Gaussianity and statistical isotropy, which is expected to be reached with *LiteBIRD* when MFs are applied on maps of the polarization modulus. We compute MFs on 300 $P^2$ simulations of the fiducial and test data sets for both *LiteBIRD* and *Planck* within the *Planck* 2018 common polarization mask shown in Fig. 1. This analysis is performed in needlet domain. To match the HEALPIX resolution of the needlet maps, the mask is downgraded accordingly, if needed, as described in Section 5.4.

In this case, we consider the following values for $\alpha_f$: {0.01, 0.05, 0.1, 0.3}. Some of these amplitudes are in the range of the required subtraction of the Galactic emission to possibly detect CMB primordial B modes. The angular power spectra of such rescaled non-Gaussian anisotropic foreground $P^2$ signal are compared with those of simulated CMB maps in the left panel of Fig. 7. All angular power spectra are computed adopting the *Planck* common mask of Fig. 1.

The obtained results for all needlet scales are shown in Table 3. For $\alpha_f = 0.01$, the non-Gaussian signal is so low that $P^2$ MFs of the test data set are fully compatible with the Gaussian and isotropic hypotheses in both *LiteBIRD* and *Planck*. For a signal with $\alpha_f = 0.05$, *LiteBIRD* yields a highly significant detection at needlet scale $j = 1$ with all MFs and a possible detection at the largest angular scales ($j = 0$) in $V_1$ and $V_2$. No deviations, instead, are observed either at higher needlet scales or in *Planck*. When $\alpha_f = 0.1$, we find a clear detection for *LiteBIRD* at $j = 0$ in $V_1$ and $V_2$, at $j = 1$ in all MFs, and hints for a detection at $j = 2$ in $V_1$ and $V_2$, while, for *Planck*, deviations only at $j = 0$ in $V_1$ and $V_2$. Finally, with $\alpha_f = 0.3$, for *LiteBIRD* we have a detection at $j = 0, 1, 2$, and 3, while for *Planck* we can observe a detection at $j = 0$ and $j = 1$ and hints for a deviation at $j = 2$.

The reported results highlight the highly enhanced detection power of MFs if applied to *LiteBIRD* $P^2$ data in comparison with *Planck* at all the angular scales of interest. Therefore, with *LiteBIRD*, we will have much better sensitivity to any possible residual contamination from foregrounds or other systematics and to cosmological effects that induce deviations from Gaussianity and statistical isotropy in the polarization data. We can note that, in general, $V_1$ and $V_2$ can better detect deviations than $V_0$. Moreover, for all the considered cases,

<footnote n="7">⁷https://pla.esac.esa.int</footnote>
<table><caption>**Table 2.** Deviations of MFs computed on *Planck* SMICA and SEVEM CMB needlet $P^2$ maps with respect to those of the corresponding end-to-end simulations. They are reported in terms of the probability-to-exceed defined in Section 5.6.</caption><thead><tr><th rowspan="3">Needlet scale</th><th rowspan="3">ℓ range</th><th colspan="6">p<sub>exc</sub> (per cent)</th></tr><tr><th colspan="2">V<sub>0</sub></th><th colspan="2">V<sub>1</sub></th><th colspan="2">V<sub>2</sub></th></tr><tr><th>SMICA</th><th>SEVEM</th><th>SMICA</th><th>SEVEM</th><th>SMICA</th><th>SEVEM</th></tr></thead><tbody><tr><td>0</td><td>2–16</td><td>46.7</td><td>43.0</td><td>15.3</td><td>12.3</td><td>11.3</td><td>4.3</td></tr><tr><td>1</td><td>8–32</td><td>60.7</td><td>35.7</td><td>37.0</td><td>43.3</td><td>67.7</td><td>31.7</td></tr><tr><td>2</td><td>16–64</td><td>47.0</td><td>75.3</td><td>83.0</td><td>87.3</td><td>79.7</td><td>86.3</td></tr><tr><td>3</td><td>32–128</td><td>88.0</td><td>36.3</td><td>83.0</td><td>93.3</td><td>75.7</td><td>2.7</td></tr><tr><td>4</td><td>64–256</td><td>64.3</td><td>90.7</td><td>82.3</td><td>53.7</td><td>31.3</td><td>38.0</td></tr><tr><td>5</td><td>128–512</td><td>7.3</td><td>33.7</td><td>59.3</td><td>39.3</td><td>31.0</td><td>74.3</td></tr><tr><td>6</td><td>256–1024</td><td>82.3</td><td>30.0</td><td>20.0</td><td>38.3</td><td>45.7</td><td>64.3</td></tr><tr><td>7</td><td>512–1024</td><td>52.7</td><td>10.0</td><td>29.7</td><td>6.3</td><td>40.3</td><td>16.7</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Angular power spectra of $P^2$ maps computed from input *Q* and *U* maps employed to simulate the fiducial and test data sets for *LiteBIRD*, SO-SATs, and *Planck*. Blue solid lines represent input CMB, while other solid lines show power spectra of injected Galactic foreground emission with different values of the overall rescaling amplitude (readers can refer to Section 5.5 for details). Left panel: Power spectra have been computed within the *Planck* common polarization mask (see Fig. 1) and input maps have been smoothed with a Gaussian beam with FWHM = 37.8 arcmin. Right panel: Power spectra have been computed within the reduced SO-SAT patch (see Fig. 4) and excluding the masked regions in the *Planck* common polarization mask. In this case, input maps have been smoothed with a Gaussian beam with FWHM = 30.0 arcmin. The adopted binning scheme is $\Delta\ell = 15$.</figcaption></figure>

unexpectedly we have more significant detections for $j = 1$ rather than for $j = 0$. This is motivated by the fact that the adopted *Planck* common polarization mask better removes foreground signal on the largest angular scales. We have repeated the same analysis computing the MFs over the full sky and, indeed, the detection power at $j = 0$ highly increases in this case.

### 6.3.2 SO-SATs versus Planck

SO-SATs are one of the most promising near-future experiments for the observation of primordial B modes from the ground. In this section, we compare the capability of detecting non-Gaussian anisotropic signals in the polarization modulus maps with MFs for SO-SATs and *Planck*. We compute MFs on 300 $P^2$ simulations of the fiducial and test data sets for both SO-SATs and *Planck* within the reduced SO-SAT patch shown in Fig. 4 excluding only the masked sources in the *Planck* common polarization mask. This analysis is performed in needlet domain. Since SO-SATs will be sensitive to modes at $\ell > 30$, we consider only needlet scales with $j \ge 2$. On these smaller scales, the residual foreground emission, to be detectable, has to be more intense. Moreover, SO-SATs will observe only regions at high Galactic latitudes (see Fig. 4), where the injected non-Gaussian anisotropic signal is expected to be fainter than that in the mask used

for *LiteBIRD* ($f_{sky} = 0.78$). Therefore, we adopt the following values for $\alpha_f$: {0.1, 0.3, 0.5, 0.7}. The angular power spectra of such rescaled non-Gaussian anisotropic foreground $P^2$ signal are compared with those of simulated CMB maps in the right panel of Fig. 7. All angular power spectra are computed within the reduced SO-SAT sky patch of Fig. 4 and excluding the masked regions in the *Planck* common polarization mask.

The obtained results are shown in Table 4. In the case of an overall rescaling factor $\alpha_f = 0.1$, $P^2$ MFs are not sensitive to the injected non-Gaussian anisotropic signal in the SO-SAT patch neither for SO nor for *Planck*. Instead, by computing $V_1$ and $V_2$ in SO-SAT data, we can detect deviations at $j = 2$ for $\alpha_f = 0.3$, at $j = 2, 3$, and 4 for $\alpha_f = 0.5$ and $\alpha_f = 0.7$, and hints for deviations at $j = 3$ for $\alpha_f = 0.3$. For the two largest considered rescaling amplitudes ($\alpha_f = 0.5$ and 0.7), also $V_1$ and $V_2$ of *Planck* needlet $P^2$ test data set at $j = 2$ significantly deviate from the fiducial Gaussian case.

Therefore, for SO-SATs the improved sensitivity will highly enhance the detection power of MFs in $P^2$ data with respect to the current *Planck* maps within the SO-observed patch. As in the comparison between *LiteBIRD* and *Planck*, the first MF $V_0$ results less sensitive to non-Gaussian anisotropic signals than $V_1$ and $V_2$. In this analysis, we assumed a pessimistic 1/f noise component for SO-SATs. We thus expect even a higher detection
<table><caption>**Table 3.** Deviations of MFs computed on needlet $P^2$ maps of the test data set with respect to those of the Gaussian fiducial data set for *LiteBIRD* and *Planck*. MFs are computed within the *Planck* common mask for polarization shown in Fig. 1 ($f_{sky} \sim 0.78$). The deviations are expressed in terms of the average probability in percentage for a $\chi^2$ distribution with $N_{dof} = 28$ (the adopted number of thresholds) to exceed the computed $\chi^2$s. Values lower than $10^{-6}$ are approximated with 0.0. Values close to 50 per cent mean no detection of the injected non-Gaussian anisotropic signal, while values very close to 0 per cent mean detection at very high significance. Different magnitudes of the overall rescaling factor $\alpha_f$ of the residual foreground contamination are considered.</caption><thead><tr><th rowspan="2">Needlet scale</th><th rowspan="2">ℓ range</th><th colspan="2">V<sub>0</sub></th><th colspan="2">V<sub>1</sub></th><th colspan="2">V<sub>2</sub></th></tr><tr><th>Planck</th><th>LiteBIRD</th><th>Planck</th><th>LiteBIRD</th><th>Planck</th><th>LiteBIRD</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">α<sub>f</sub> = 0.01</td></tr><tr><td>0</td><td>2–16</td><td>50.2</td><td>47.5</td><td>51.7</td><td>49.4</td><td>50.2</td><td>49.6</td></tr><tr><td>1</td><td>8–32</td><td>52.2</td><td>51.4</td><td>54.5</td><td>52.3</td><td>50.5</td><td>50.0</td></tr><tr><td>2</td><td>16–64</td><td>49.1</td><td>49.0</td><td>51.3</td><td>49.0</td><td>50.6</td><td>50.8</td></tr><tr><td>3</td><td>32–128</td><td>49.9</td><td>50.2</td><td>51.9</td><td>50.0</td><td>50.1</td><td>49.4</td></tr><tr><td>4</td><td>64–256</td><td>49.6</td><td>50.8</td><td>50.8</td><td>50.0</td><td>49.9</td><td>48.8</td></tr><tr><td>5</td><td>128–512</td><td>49.1</td><td>50.4</td><td>49.1</td><td>50.2</td><td>49.5</td><td>51.9</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">α<sub>f</sub> = 0.05</td></tr><tr><td>0</td><td>2–16</td><td>49.3</td><td>51.2</td><td>45.5</td><td>14.2</td><td>43.2</td><td>10.8</td></tr><tr><td>1</td><td>8–32</td><td>47.7</td><td>0.0</td><td>54.2</td><td>0.01</td><td>53.7</td><td>0.4</td></tr><tr><td>2</td><td>16–64</td><td>49.4</td><td>45.7</td><td>52.9</td><td>48.1</td><td>54.8</td><td>50.3</td></tr><tr><td>3</td><td>32–128</td><td>52.1</td><td>49.4</td><td>52.3</td><td>51.6</td><td>51.1</td><td>49.1</td></tr><tr><td>4</td><td>64–256</td><td>48.9</td><td>50.1</td><td>50.7</td><td>49.1</td><td>50.6</td><td>48.6</td></tr><tr><td>5</td><td>128–512</td><td>47.4</td><td>51.3</td><td>49.8</td><td>51.0</td><td>49.8</td><td>50.2</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">α<sub>f</sub> = 0.1</td></tr><tr><td>0</td><td>2–16</td><td>49.2</td><td>53.5</td><td>20.8</td><td>0.07</td><td>20.4</td><td>0.6</td></tr><tr><td>1</td><td>8–32</td><td>42.7</td><td>0.0</td><td>43.9</td><td>0.0</td><td>49.9</td><td>0.0</td></tr><tr><td>2</td><td>16–64</td><td>45.3</td><td>37.9</td><td>53.9</td><td>23.7</td><td>52.3</td><td>23.5</td></tr><tr><td>3</td><td>32–128</td><td>51.2</td><td>49.6</td><td>51.6</td><td>49.6</td><td>50.6</td><td>48.2</td></tr><tr><td>4</td><td>64–256</td><td>45.2</td><td>50.8</td><td>51.7</td><td>49.7</td><td>51.2</td><td>48.8</td></tr><tr><td>5</td><td>128–512</td><td>45.9</td><td>52.4</td><td>48.6</td><td>50.8</td><td>49.2</td><td>50.7</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">α<sub>f</sub> = 0.3</td></tr><tr><td>0</td><td>2–16</td><td>52.1</td><td>29.3</td><td>0.002</td><td>0.0</td><td>0.7</td><td>0.001</td></tr><tr><td>1</td><td>8–32</td><td>0.0002</td><td>0.0</td><td>0.00004</td><td>0.0</td><td>1.6</td><td>0.0</td></tr><tr><td>2</td><td>16–64</td><td>26.0</td><td>0.0</td><td>10.0</td><td>0.0</td><td>11.9</td><td>0.0</td></tr><tr><td>3</td><td>32–128</td><td>52.5</td><td>43.2</td><td>53.9</td><td>0.6</td><td>53.4</td><td>3.5</td></tr><tr><td>4</td><td>64–256</td><td>45.3</td><td>53.2</td><td>52.6</td><td>36.7</td><td>51.3</td><td>33.2</td></tr><tr><td>5</td><td>128–512</td><td>47.0</td><td>51.5</td><td>48.7</td><td>48.2</td><td>50.1</td><td>48.3</td></tr></tbody></table>

power in the first needlet bands in case of improved instrumental properties.

## 7 CONCLUSIONS

MFs have been widely applied to the analysis of CMB data (Schmalzing & Górski 1998; Planck Collaboration VII 2020d), but their application to polarization has been limited to date. In this work, we introduce a generalized framework to compute the MFs on maps of the squared modulus of polarization ($P^2 = Q^2 + U^2$). MFs are able to provide information complementary to that obtained from angular power spectra and, in particular, it is sensitive to deviations from the hypothesis of a Gaussian isotropic field. Therefore, they are expected to be able to blindly detect (i) early- and late-time cosmological effects (e.g. primordial non-Gaussianity or CMB lensing) and (ii) significant residual contamination by foregrounds and instrumental systematics.

The statistical analysis of the polarized intensity is complementary to that performed on E or B modes, as the CMB P map reflects the local morphological properties of the polarization field. Therefore, we expect the introduced formalism to better capture those deviations that primarily affect the intensity rather than the direction of the polarization field. Furthermore, it allows us to easily study polarization

data in the presence of any mask without suffering from any E-B and B-E leakage contamination (Lewis et al. 2001; Bunn et al. 2003). An alternative formalism to study the spin-2 CMB polarization field has recently been introduced in Carrón Duque et al. (2023).

Our main contributions are the following:

(i) We introduce a general formalism to predict the value of MFs for $P^2$ in the case of a Gaussian and isotropic field, justifying how the effects of the spin are negligible in the cosmological case.

(ii) We implement and release PYNKOWSKI, a public PYTHON package to compute MFs on HEALPIX maps of any scalar field, such as CMB T and $P^2$.

(iii) We validate theoretical predictions for MFs from $P^2$ on Gaussian isotropic CMB simulations.

(iv) We apply the formalism to Planck polarization modulus maps (produced by SMICA and SEVEM) both in pixel space and in needlet domain. We find compatibility (with less than $1.5\sigma$ deviation in all cases) between their MFs and those computed on realistic end-to-end simulations, which take into account the anisotropic noise residuals due to the scanning strategy of the satellite. We can conclude that we do not observe any hint for deviations from Gaussianity or statistical isotropy in Planck CMB polarization maps at any angular scale. However, given the noise level in the Planck $P^2$ maps, such an analysis allows us to assess mainly the statistical properties of noise.
<table><caption>**Table 4.** Deviations of MFs computed on needlet $P^2$ maps of the test data set with respect to those of the Gaussian fiducial data set for SO-SATs and *Planck*. MFs are computed within the reduced SO-SAT patch shown in Fig. 4 and excluding the masked regions of the *Planck* polarization common mask ($f_{sky} \sim 0.1$). The deviations are expressed in terms of the average probability in percentage for a $\chi^2$ distribution with $N_{dof} = 28$ (the adopted number of thresholds) to exceed the computed $\chi^2$s. Values lower than $10^{-6}$ are approximated with 0.0. Values close to 50 per cent mean no detection of the injected non-Gaussian anisotropic signal, while values very close to 0 per cent mean detection at very high significance. Different magnitudes of the overall amplitude $\alpha_f$ of the residual foreground contamination are considered.</caption><thead><tr><th rowspan="2">Needlet scale</th><th rowspan="2">ℓ range</th><th colspan="2">V<sub>0</sub></th><th colspan="2">V<sub>1</sub></th><th colspan="2">V<sub>2</sub></th></tr><tr><th>Planck</th><th>SO-SATs</th><th>Planck</th><th>SO-SATs</th><th>Planck</th><th>SO-SATs</th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">α<sub>f</sub> = 0.1</td></tr><tr><td>2</td><td>30–64</td><td>47.9</td><td>53.1</td><td>54.8</td><td>33.2</td><td>56.0</td><td>39.0</td></tr><tr><td>3</td><td>32–128</td><td>50.6</td><td>51.4</td><td>52.6</td><td>51.3</td><td>49.7</td><td>50.3</td></tr><tr><td>4</td><td>64–256</td><td>48.9</td><td>46.4</td><td>50.1</td><td>49.6</td><td>50.5</td><td>49.5</td></tr><tr><td>5</td><td>128–512</td><td>49.4</td><td>48.6</td><td>49.2</td><td>50.4</td><td>49.2</td><td>48.4</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">α<sub>f</sub> = 0.3</td></tr><tr><td>2</td><td>30–64</td><td>50.3</td><td>54.2</td><td>35.3</td><td>0.0</td><td>26.2</td><td>0.003</td></tr><tr><td>3</td><td>32–128</td><td>52.8</td><td>49.7</td><td>56.1</td><td>18.0</td><td>54.0</td><td>31.9</td></tr><tr><td>4</td><td>64–256</td><td>52.6</td><td>49.6</td><td>49.7</td><td>48.6</td><td>49.6</td><td>49.2</td></tr><tr><td>5</td><td>128–512</td><td>48.5</td><td>50.0</td><td>49.8</td><td>48.6</td><td>49.0</td><td>49.3</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">α<sub>f</sub> = 0.5</td></tr><tr><td>2</td><td>30–64</td><td>51.6</td><td>47.5</td><td>0.04</td><td>0.0</td><td>0.4</td><td>0.0</td></tr><tr><td>3</td><td>32–128</td><td>50.3</td><td>51.4</td><td>54.0</td><td>0.0</td><td>50.7</td><td>0.3</td></tr><tr><td>4</td><td>64–256</td><td>50.7</td><td>49.7</td><td>51.5</td><td>0.0</td><td>52.3</td><td>4.2</td></tr><tr><td>5</td><td>128–512</td><td>47.8</td><td>48.8</td><td>49.0</td><td>47.7</td><td>48.9</td><td>48.0</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">α<sub>f</sub> = 0.7</td></tr><tr><td>2</td><td>30–64</td><td>49.3</td><td>26.8</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr><tr><td>3</td><td>32–128</td><td>50.1</td><td>47.1</td><td>26.0</td><td>0.0</td><td>26.0</td><td>0.0</td></tr><tr><td>4</td><td>64–256</td><td>50.6</td><td>51.0</td><td>52.1</td><td>0.0</td><td>52.1</td><td>0.0</td></tr><tr><td>5</td><td>128–512</td><td>46.7</td><td>47.6</td><td>48.8</td><td>46.3</td><td>48.8</td><td>46.8</td></tr></tbody></table>

In the coming years, we expect new high-quality CMB polarization observations from ongoing or planned experiments such as SPT (Sayre et al. 2020), LSPE (Addamo et al. 2021), SO (Ade et al. 2019), CMB-S4 (Abazajian et al. 2022), and *LiteBIRD* (Hazumi et al. 2019). MFs can play a key role in detecting unexpected non-Gaussian components or departures from statistical isotropy in these polarization maps. To further assess this point, we forecast the detection power of $P^2$ MFs for two future and complementary CMB experiments, the *LiteBIRD* satellite and the ground-based SO-SATs, which will target the observation of primordial CMB polarization B modes. We show that, when applied to simulated maps of the polarization modulus of *LiteBIRD* and SO, MFs are much more sensitive to the presence of a non-Gaussian anisotropic signal than *Planck*. We have shown this by using a residual Galactic foreground contamination model with varying amplitudes as an example of non-Gaussian anisotropic component.

Furthermore, MFs might be especially relevant in light of the large-scale anomalies that are being detected in multiple observables (Gruppuso et al. 2013; Adhikari, Shandera & Erickcek 2016; Schwarz et al. 2016). The introduced formalism may also be useful to study other complex spin fields defined on the sphere, such as the convergence ($\kappa$) or the gravitational shear ($\sigma$) maps associated with weak lensing. Lensing maps, indeed, are found to be highly non-Gaussian, as they are produced by the non-linear distribution of matter in the Universe (Zürcher et al. 2021).

Finally, MFs can be employed to study polarized Galactic and extragalactic foregrounds, as they strongly deviate from a Gaussian distribution and are (in the Galactic case) highly anisotropic. MFs can complement other higher order statistics that are already used to analyse foregrounds, such as peak statistics (Carrón Duque et al.

2019) or neural networks (Farsian, Krachmalnicoff & Baccigalupi 2020). As an example of this synergy, in Krachmalnicoff & Puglisi (2021) they train a neural network to generate a high-resolution dust foreground map with a given shape of their MFs, as dust emission is highly non-Gaussian (Ben-David, von Hausegger & Jackson 2015; Jung, Racine & van Tent 2018).

The PYNKOWSKI PYTHON package we have developed to study MFs is fully documented and can be found in https://github.com/javicarron/pynkowski, along with examples of its use.

## ACKNOWLEDGEMENTS

MM and NV acknowledge support by Agenzia Spaziale Italiana (ASI)/COSMOS grant no. 2016-24-H.0 and grant no. 2020-9-HH.0. JCD and MM acknowledge support from ASI/Instituto Nazionale di Fisica Nucleare (INFN) grant no. 2021-43-HH.0. Part of this work was also supported by the InDark INFN project. DM acknowledges support from the Ministero dell'Università e della Ricerca 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.

## DATA AVAILABILITY

The *Planck* data used in this article are available in the *Planck* Legacy Archive, at https://pla.esac.esa.int/pla/. The results produced in this article will be shared on reasonable request to the corresponding author. The software is available at https://github.com/javicarron/py-nkowski.
REFERENCES

Abazajian K. et al., 2022, preprint (arXiv:2203.08024)

Addamo G. et al., 2021, J. Cosmol. Astropart. Phys., 2021, 008

Ade P. et al., 2019, J. Cosmol. Astropart. Phys., 2019, 056

Adhikari S., Shandera S., Erickcek A. L., 2016, Phys. Rev. D, 93, 023524

Adler R. J., Taylor J. E., 2007, Random Fields and Geometry. Springer, New York

Aiola S. et al., 2020, J. Cosmol. Astropart. Phys., 2020, 047

Andrews A., Jasche J., Lavaux G., Schmidt F., 2023, MNRAS, 520, 5746

Barrow J. D., Hervik S., 2006, Phys. Rev. D, 73, 023007

Bartolo N., Komatsu E., Matarrese S., Riotto A., 2004, Phys. Rep., 402, 103

Ben-David A., von Hausegger S., Jackson A. D., 2015, J. Cosmol. Astropart. Phys., 2015, 019

Bennett C. L. et al., 2013, ApJS, 208, 20

Bunn E. F., Zaldarriaga M., Tegmark M., de Oliveira-Costa A., 2003, Phys. Rev. D, 67, 023501

Carrón Duque J., Buzzelli A., Fantaye Y., Marinucci D., Schwartzman A., Vittorio N., 2019, Astron. Comput., 28, 100310

Carrón Duque J., Carones A., Marinucci D., Migliaccio M., Vittorio N., 2023, preprint (arXiv:2301.13191)

Chen X., 2010, Adv. Astron., 2010, 638979

Chen Z., Yu Y., Liu X., Fan Z., 2020, ApJ, 897, 14

Coles P., Barrow J. D., 1987, MNRAS, 228, 407

Delabrouille J., Cardoso J. F., Patanchon G., 2003, MNRAS, 346, 1089

Desjacques V., Seljak U., 2010, Class. Quantum Gravity, 27, 124011

do Carmo M. P., 1976, Differential Geometry of Curves and Surfaces. Prentice Hall, Englewood Cliffs, NJ, USA

Ducout A., Bouchet F. R., Colombi S., Pogosyan D., Prunet S., 2012, MNRAS, 429, 2104

Fantaye Y. T., Marinucci D., 2014, in Heavens A., Starck J.L., Krone-Martins A., eds, Proc. IAU Symp. 306, Statistical Challenges in 21st Century Cosmology. Cambridge Univ. Press, Cambridge, p. 135

Fantaye Y., Marinucci D., Hansen F., Maino D., 2015, Phys. Rev. D, 91, 063501

Farsian F., Krachmalnicoff N., Baccigalupi C., 2020, J. Cosmol. Astropart. Phys., 2020, 017

Ganesan V., Chingangbam P., Yogendran K. P., Park C., 2015, J. Cosmol. Astropart. Phys., 2015, 028

Geller D., Marinucci D., 2010, J. Fourier Anal. Appl., 16, 840

Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

Gott J. R., III, Park C., Juszkiewicz R., Bies W. E., Bennett D. P., Bouchet F. R., Stebbins A., 1990, ApJ, 352, 1

Grewal N., Zuntz J., Tröster T., Amon A., 2022, Open J. Astrophys., 5, 13

Gruppuso A., Natoli P., Paci F., Finelli F., Molinari D., De Rosa A., Mandolesi N., 2013, J. Cosmol. Astropart. Phys., 2013, 047

Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399

Hazumi M. et al., 2019, J. Low Temp. Phys., 194, 443

Jung G., Racine B., van Tent B., 2018, J. Cosmol. Astropart. Phys., 2018, 047

Kamionkowski M., Kosowsky A., Stebbins A., 1997, Phys. Rev. Lett., 78, 2058

Kasak J., Creswell J., Naselsky P., Liu H., 2021, Phys. Rev. D, 104, 023502

Krachmalnicoff N., Puglisi G., 2021, ApJ, 911, 42

Lan X., Marinucci D., Xiao Y., 2018, Stoch. Process. Appl., 128, 1294

Lerario A., Marinucci D., Rossi M., Stecconi M., 2022, preprint (arXiv:2207.08413)

Lewis A., Challinor A., Turok N., 2001, Phys. Rev. D, 65, 023505

LiteBIRD Collaboration, 2023, Prog. Theor. Exp. Phys., 2023, 042F01

Martínez-González E., Diego J. M., Vielva P., Silk J., 2003, MNRAS, 345, 1101

Matsumura T. et al., 2014, J. Low Temp. Phys., 176, 733

Narcowich F. J., Petrushev P., Ward J. D., 2006, SIAM J. Math. Anal., 38, 574

Naselsky P. D., Novikov D. I., 1998, ApJ, 507, 31

Newman E. T., Penrose R., 1966, J. Math. Phys., 7, 863

Peña G. A., Candlish G. N., 2022, MNRAS, 511, 2259

Pietrobon D., Balbi A., Marinucci D., 2006, Phys. Rev. D, 74, 043524

Pitrou C., Pereira T. S., Uzan J.-P., 2008, J. Cosmol. Astropart. Phys., 2008, 004

Planck Collaboration I, 2011, A&A, 536, A1

Planck Collaboration XVI, 2016, A&A, 594, A16

Planck Collaboration IV, 2020a, A&A, 641, A4

Planck Collaboration V, 2020b, A&A, 641, A5

Planck Collaboration VI, 2020c, A&A, 641, A6

Planck Collaboration VII, 2020d, A&A, 641, A7

Santos L., Wang K., Zhao W., 2016, J. Cosmol. Astropart. Phys., 2016, 029

Sayre J. T. et al., 2020, Phys. Rev. D, 101, 122003

Schmalzing J., Górski K. M., 1998, MNRAS, 297, 355

Schwarz D. J., Copi C. J., Huterer D., Starkman G. D., 2016, Class. Quantum Gravity, 33, 184001

Thorne B., Dunkley J., Alonso D., Næss S., 2017, MNRAS, 469, 2821

Tomita H., 1986, Prog. Theor. Phys., 76, 952

Zaldarriaga M., Seljak U., 1997, Phys. Rev. D, 55, 1830

Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, J. Open Source Softw., 4, 1298

Zürcher D., Fluri J., Sgier R., Kacprzak T., Refregier A., 2021, J. Cosmol. Astropart. Phys., 2021, 028

APPENDIX A: ROBUSTNESS TO CMB COMPONENT SEPARATION ALGORITHM

In Section 6.2, we showed the results when MFs are computed on CMB Planck polarization maps and realistic end-to-end simulations, both produced by the SMICA cleaning algorithm, assuming primordial Gaussianity and statistical isotropy, obtaining a complete compatibility between them. We perform the same analysis with the SEVEM reconstructed maps in order to verify that the results are robust to the specific method used to extract the CMB signal.

In Fig. A1, we show the comparison between MFs of the Planck polarization modulus map and the ones of realistic simulations obtained with the SEVEM pipeline. As with SMICA, we can observe full compatibility between them with deviations always within ±2σ, where σ represents the dispersion among the simulations. Therefore, we find no indication of primordial non-Gaussianity or deviations from statistical isotropy of the observed CMB polarization modulus.
<figure><img src="image_10.png" /><figcaption>**Figure A1.** MFs of the Planck SEVEM-released $P^2 = Q^2 + U^2$ maps (dots) are compared with the mean given by realistic end-to-end simulations (solid lines). The first, second, and third MFs are shown (top to bottom), with the values (top section of each panel) and the residual with respect to the mean of simulations (bottom section of each panel). In all plots, the pink shaded area represents the ±1σ region obtained in the simulations.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure A2.** Deviations of theoretical MFs (dots) in the Gaussianity and statistical isotropy hypotheses with respect to those computed on the Planck SMICA $P^2 = Q^2 + U^2$ end-to-end simulations (solid lines). The first, second, and third MFs are shown (top to bottom). In all plots, the pink shaded area represents the ±1σ region obtained in the simulations.</figcaption></figure>
APPENDIX B: IMPACT OF ANISOTROPIC
NOISE

Planck CMB polarization maps are known to be affected by noise
residuals that are significantly anisotropic due to the mission scan-
ning strategy (Planck Collaboration I 2011), as it can be seen in
Fig. 2. Some portions of the sky are better observed; this reduces
the noise level and, consequently, the variance of the map in these
regions. The anisotropic noise pattern is clearly visible in the SMICA
and SEVEM P² maps. Hence, it needs to be properly modelled and
accounted for when analysing the data looking for non-Gaussianity
and departure from isotropy of primordial origin. The impact of
anisotropic noise has been studied for Planck CMB temperature by
Ducout et al. (2012). To stress this point, let us study the impact of
anisotropic noise contamination by comparing the MFs of SMICA
and SEVEM realistic simulations with the theoretical expected
values under the Gaussianity and isotropy hypotheses computed
considering the average value of μ among all simulations and using
equation (6).

We follow an analogous procedure to the one in Section 6.2. The
deviations of theoretical expectations with respect to MFs computed
on SMICA end-to-end simulations are reported in Fig. A2. The
deviations are strongly significant at a level of 20σ–45σ, detectable in
all MFs and both component separation algorithms. It is remarkable
that this strong deviation can be reproduced exactly by realistic
noise simulations (compare Figs 6 and A1, where noise + CMB
simulations are perfectly compatible with observations). This result
highlights the potential of the application of MFs to future data,
which will not be dominated by anisotropic noise as strongly as
Planck maps. It also underlines the importance of calibration to
realistic simulations, especially in polarization data. Any possible
future claim of detection of primordial non-Gaussianity or deviations
from statistical isotropy of polarized CMB signal will have to be
calibrated not only to realistic noise, but also to possible realistic
residuals of Galactic and extragalactic foregrounds.

This paper has been typeset from a TeX/LA TeX file prepared by the author.