Minkowski Functionals of CMB polarisation intensity with Pynkowski: theory and application to Planck and future data

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 instead 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 MFs formalism to the CMB polarisation 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 which include CMB and instrumental noise residuals. We find no significant deviations from Gaussianity or isotropy in Planck CMB polarisation intensity. However, MFs could play an important role in the analysis of CMB polarisation 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 $\texttt{Pynkowski}$ (https://github.com/javicarron/pynkowski).


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 ΛCDM cosmological model (Planck Collaboration 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;Pitrou et al. 2008;Barrow & Hervik 2006;Chen 2010).Such effects, however, are not associated only to 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 ★ E-mail: alessandro.carones@roma2.infn.itsignals contain relevant information which is not reflected in the angular power spectrum.Non-Gaussian statistics of CMB lensing have already been shown to help constrain cosmological parameters (Zürcher et al. 2021;Chen et al. 2020;Grewal et al. 2022).
Recently, many studies about primordial non-Gaussianity have also involved the analysis of the Large Scale Structure of the Universe (Andrews et al. 2023;Peña & Candlish 2022).In particular, the highorder 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 et al. 1990), and have become a popular tool in CMB analysis after the discussion by Schmalzing & Górski (1998).An extension of the MFs 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).
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 • The polarization modulus traces the local properties of the polarization field, while the E-and B-modes maps are obtained from a non-local decomposition of Q and U, typically performed in harmonic space. • The 2 field can be robustly obtained for any adopted mask without suffering of E-B and B-E leakage contamination (Lewis et al. 2001;Bunn et al. 2003), which may distort the statistics of the reconstructed E-and B-modes.
• The morphology of the polarization modulus is intrinsically different and can be more effective in detecting those deviations from Gaussianity and statistical isotropy which 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 do not exhaust the morphological characterisation of the polarization field.This motivates the usefulness of analysing complementary quantities as  2 .
This paper has 5 key goals: (i) to introduce a more general formalism to compute the theoretical predictions of scalar maps, such as the CMB temperature anisotropies (), (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 ( 2 ) in the case of Gaussianity and statistical isotropy, (iii) to provide the community with a public Python package called Pynkowski1 to easily estimate MFs on HEALPix 2 maps (Górski et al. 2005) and compute their Gaussian isotropic expectations.
(iv) to compute the MFs on Planck  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  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  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 ACT (Aiola et al. 2020), SPT (Sayre et al. 2020), LSPE (Addamo et al. 2021), Simons Observatory (Ade et al. 2019), 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  2 maps to detect deviations from Gaussianity and statistical isotropy for two future and complementary CMB polarization experiments: Simons Observatory (SO) and LiteBIRD.Both of them will target the detection of the so-called primordial CMB polarization  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 (  sky ∼ 30%) and the atmospheric contamination, SO is sensitive to polarization modes with ℓ > 30 and thus will only target the detection of the recombination peak in the primordial -mode power spectrum at ℓ ∼ 100 (Ade et al. 2019).LiteBIRD, instead, is a space-borne L-class mission selected by the Japanaese Aerospace Exploration Agency (JAXA) in 2019.Given its full-sky coverage, it will be sensitive also to modes on the largest angular scales (LiteBIRD Collaboration et al. 2023).For both experiments, the detection of primordial CMB  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 LiteBIRD and SO-SATs polarization simulated data 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:  2 =  2 +  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  2 MFs to a non-Gaussian anisotropic signal for LiteBIRD and SO.Finally, in Section 7 we report our conclusions.

MINKOWSKI FUNCTIONALS ON THE SPHERE
MFs on the sphere are now well-established tools for CMB data analysis (Planck Collaboration 2020d); nonetheless, we will first recall some background notation and definitions in order to make the comparison with our present work more explicit.Given a scalar random field  (.) observed on the unit sphere S 2 (e.g., CMB temperature anisotropies), the excursion set at a threshold  is defined by We will omit the arguments of   to refer to the excursion sets when  and the domain of  are clear by context.As recalled in Schmalzing & Górski (1998), the morphological properties of these excursion sets can be summarised by the three MFs defined as follows: where  denotes a line element along the boundary of the excursion set and (.) denotes the geodesic curvature of said boundary.Thus,  0 represents the excursion area, while  1 is one fourth of the boundary length. 2 is associated with the number of connected regions minus the number of holes, which on the plane corresponds to the Euler-Poincaré characteristic .On a spherical surface the Euler-Poincaré characteristic does no longer correspond exactly to  2 , but includes an extra term which can be computed by means of the generalised Gauss-Bonnet Theorem, which yields: 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 These quantities are defined as follows.Let us define a tube of width  built around the manifold   as all those points at a distance less than  from   .Its volume can be exactly expressed as a finite Taylor expansion on , whose coefficients correspond to the Lipschitz-Killing curvatures of   (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 twodimensional sphere, but it is convenient to recall them from a more general perspective.
Without loss of generality we assume the field normalised to have unit variance.We define  as the derivative of the covariance function at the origin for the field  .This quantity can be computed as follows: where the sequence { ℓ } denotes, as usual, the angular power spectrum of the field  .Furthermore, following Adler & Taylor (2007) we define the flag coefficients as In general,   represents the volume of the -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): where and in general, 2 denotes the sequence of Hermite polynomials.The function Φ 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 normalised by area.In order, we have: We note that this quantity depends only on , not on the field itself as long as it is normalised to have unit variance. Likewise: so that Finally: and hence: 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, (iv) a set of "density" functions   , 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 formthe only change is in the "density" functions   , which nonetheless can be expressed with a fully explicit analytic 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 normalisation of the MFs.

MINKOWSKI FUNCTIONALS ON POLARIZATION FIELDS
A spin  random field   is a complex-valued field that satisfies the following relation at any fixed point : ′  () = exp()  () .where  ′ () represents the value in  ∈ S 2 when the local coordinates are rotated by an angle .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, that can be written as: where  and  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 realisation 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  = 0, we are back to the case of a complexvalued isotropic Gaussian map, which can always be viewed as the complexification of two Gaussian independent and identically distributed real-valued fields  and .
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 ℓ * , the effect of the spin vanishes as ℓ * 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  ≠ 0 field with the spin  = 0 case) are of the form ∝  ℓ (ℓ+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 , and even under other conditions such as the (non-physical) case of a spin parameter growing with ℓ * .We can informally summarise this mathematical result as follows: Up to negligible terms which depend only on the smallest multipoles, the expected values of the Lipschitz-Killing Curvatures for the excursion sets of the norms of spin  random fields are identical, for all values of .In other words, defining for all values of  and .
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 Kasak et al. (2021), Naselsky & Novikov (1998), 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 generalisations.
In particular, assume that we have  and  two independent copies of Gaussian random fields with the same angular power spectrum, expected value 0 and variance normalised to 1 (i.e.
it is immediate to demonstrate that  2 is a chi-squared with two degrees of freedom,  2 ∼  2 2 , and hence it has marginal density We can now consider the excursion sets We can now recall again the GKF, which in the case of chi-square fields reads where as usual we write  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   (.) with  ; 2 (.) defined by Adler & Taylor (2007, Theorem 15.10.1): More explicitly, for  = 1 we obtain whereas for  = 2 we have Plugging these expressions into equation ( 5) we obtain, for the bound-ary length, The excursion area can be immediately computed as: Finally, for the Euler-Poincaré Characteristic Passing from the Lipschitz-Killing curvatures to MFs and normalising by the area of the sphere, we finally obtain the following predictions for a chi-squared map with 2 degrees of freedom, such as  2 : where the threshold is non-negative,  ≥ 0. For negative thresholds, it follows from the definition that  0 is exactly 1, while  1 and  2 are exactly 0. In particular,  2 presents a discontinuity at  = 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  is the derivative of the covariance of the field at the origin, that can be computed as before: where: with   ℓ and   ℓ , respectively, the EE and BB angular power spectra computed from Q and U maps (Kamionkowski et al. 1997;Zaldarriaga & Seljak 1997).

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  (see Section 2), or for  2 fields, like CMB polarized intensity  2 (see Section 3).They can be computed given either an input angular power spectrum or .
(ii) To determine the MFs on any scalar map, either on the full sky or within any mask.
Both  and  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  ().We compute the MFs normalised by the area of the sphere (  =   4  ), so the extension to masked maps is straight-forward.
Without loss of generality, for this work we normalise input maps to have unit variance, as we did in the previous sections when computing the theoretical predictions.In the case of  2 , both the  and  maps are separately normalised to unit variance.The MFs can now be expressed as a function of the a-dimensional threshold .

First MF, v 0
The first MF,  0 , is computed according to the equation: where Θ() is the Heaviside function (1 where  ≥ 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, 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,  0 () is computed as the fraction of the map above the threshold .

Second MF, v 1
The second MF,  1 , is computed according to the equation: where the second equality is the result of a change of coordinates from the line element  of the boundary of the excursion set to a surface element  on the sphere.The gradient of  is denoted by ∇  and  is the Dirac delta.We start by computing |∇  | at every pixel.Then, in order to approximate the integral with a sum over pixels, we provide a bin of Δ.We create a mask   () identifying the pixels where the value of the map is close to : Now,  1 () can be computed as the mean of the pixels of   • |∇  |, divided over Δ, the bin size, and 4, the normalisation factor.
Third MF, v 2 The third MF,  2 , is computed according to the equation: 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   identically and compute the mean of the pixels of This quantity is then divided by the bin size Δ and by 2 to yield the result for  2 ().
In this formula,  is the geodesic curvature of the boundary of the excursion set.Given that this curve is given by the implicit equation  () = , we can use the expression of  for an implicit curve (see, e.g., do Carmo 1976): where the semicolon denotes the covariant derivative.All the quantities involved in this expression can be computed for each pixel.

Computation of the derivatives
The spatial derivatives of the maps (both partial and covariant) are computed with the help of the healpy3 function alm2map_der in harmonic space (Zonca et al. 2019).The covariant derivatives of a field  on the sphere are computed in terms of the partial derivatives: where we consider  ∈ (−/2, /2) as the latitude and  ∈ [0, 2) 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 ℓ  = 3 •  side − 1; i.e., the derivatives can only be computed if the maps are reasonably smooth.This becomes a requirement in order to estimate  1 and  2 , while  0 can be computed for any map as it does not require the computation of the derivatives.
This smoothness condition is satisfied for realistic  and  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 twicedifferentiable almost everywhere (for other cases, such as a fractal behaviour, see Lan et al. 2018).

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 fit angular power spectrum for polarization (E and B modes) as reported by Planck Collaboration (2020b).For these simulations, we fix a resolution of  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  0 , we exclude the modes with ℓ ≤ 4 in the simulated maps in order to avoid spin effects in the computation of the MF (see Section 3).Conversely, when  1 and  2 are estimated, a smoothing of the maps with a Gaussian beam with full width at half maximum (FWHM) equal to 15 ′ is applied in order not to be affected by pixelisation effects in the computation of the spatial derivatives.
MFs have been computed over the full sky maps of  2 =  2 +  2 , where  and  have been previously normalised to unit variance.We analyse the deviation of each MF at different values of the threshold  in these simulated realisations 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 MFs computations on simulated maps.

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: SMICA (Spectral Matching Independent Component Analysis) and SEVEM (Spectral Estimation Via Expectation Maximization), as reported in Planck Collaboration (2020a).
SMICA (Delabrouille et al. 2003) performs a foregrounds and noise subtraction by linearly combining Planck frequency maps in harmonic space with multipole-dependent weights, up to ℓ = 4000, which minimize on each angular scale the variance of the output map.In polarization, it combines the  and  modes independently and then recombines them to obtain the  and  maps.
SEVEM (Martínez-González et al. 2003), instead, is a templatefitting method.Different templates of the foregrounds 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 minimising the variance of the cleaned map outside a mask, which covers the point sources detected in polarization and the 3% 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 ′ and a maximum considered multipole of ℓ = 3000.
The native resolution of these released maps4 is  side = 2048; we degrade them to  side = 1024 to match exactly the pipeline applied for simulations.Similarly, before the computation of  0 we again exclude the lowest multipoles (ℓ ≤ 4) in order to avoid spin effects, while for the computation of  1 and  2 we bring the maps to a resolution of FWHM=15 ′ 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 all the analysis, as recommended by (Planck Collaboration 2020a).This mask was obtained by thresholding the standard deviation between each of the four cleaned Planck CMB maps and it is 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  maps are shown in Fig. 2.

Planck end-to-end simulations
The SMICA and SEVEM maps include residual effects associated to 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-toend simulations 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.
We use both SMICA and SEVEM simulations with  side = 1024 to make all the analysis identical in all maps.Similarly to real Planck polarization maps, we again exclude the lowest multipoles (ℓ ≤ 4) before the computation of  0 , while we bring the maps to a resolution of FWHM=15 ′ for the computation of  1 and  2 .

Needlet maps
The comparison of MFs computed on Planck  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, 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  2 maps with those of realistic end-to-end simulations (introduced, respectively, in Section 5.2 and Section 5.3) in needlet domain.A similar analysis for -mode Planck maps has been reported in Planck Collaboration (2020d).
Needlets are a wavelet system which guarantees simultaneous localisation in harmonic and pixel space; they have been introduced in the statistical literature by Narcowich et al. (2006) and firstly applied to CMB data by Pietrobon et al. (2006).Given any field of spin  defined on the sphere, the corresponding needlet maps    are obtained by filtering its harmonic coefficients   ℓ with a harmonic weighting function   (ℓ) which isolates modes at different angular scales for each needlet scale  (Geller & Marinucci 2010): where  is a direction in the sky and   ℓ 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  whose width is set by the parameter : lower values of  correspond to a tighter localisation 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  and  are components of a spin-2 field on the sphere: Therefore, in this paper, we first filter the  ±2 ℓ coefficients to obtain  and  needlet maps and then compute the corresponding  2 map.As needlet harmonic function , we have adopted the standard needlet construction (Narcowich et al. 2006) with  = 2.The needlet filters have been obtained through the Python package MTNeedlets 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: The adopted configuration of needlet harmonic filters is shown in Fig. 3.We can note that the -th needlet band is non-zero only for ℓ ≤ 2 +4 .Therefore needlet maps are band-limited and it suffices to reconstruct them at  side = 2 +3 .As it was done in Planck Collaboration (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.).The SO is composed of a Large Aperture Telescope (SO-LAT) and several refracting Small Aperture Telescopes (SO-SATs).SO-LAT will observe 40% of the sky and will be devoted to measurements of the CMB small-scale temperature and polarization anisotropies, the CMB lensing, the SZ effects, and extragalactic sources.SO-SATs, instead, aim at a detection of primordial  modes.Given the affinity with the scientific target of LiteBIRD, in this analysis we have considered only SO-SATs.

5.5
For both experiments, we simulate a fiducial Gaussian data set and a test data set which also includes a non-Gaussian anisotropic signal.In the fiducial Gaussian data set,  and  maps are obtained by adding Gaussian and isotropic CMB and noise realisations: 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: where  f is a scaling factor which 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  f to vary in order to test the sensitivity of  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-fit Planck 2018 parameters (Planck Collaboration 2020c).For the instrumental noise component, we first simulate noise maps at the different frequency channels of the considered experiment as Gaussian realisations from the following power spectrum: where  Q,U is the polarization sensitivity in K CMB •arcmin, while ℓ knee and  knee refer to the contribution from 1/f noise.The employed LiteBIRD and SO-SATs specifications are taken from the latest forecasts of the two experiments (LiteBIRD Collaboration et al. 2023;Ade et al. 2019).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 ℓ knee and  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  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 et al. (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  2 maps.For SO-SATs, this statement does not hold, because the hit counts map is very inhomogeneous.Therefore, we restrict the MFs computation to a smaller patch with  sky = 10% (shown in Fig. 4), where the scanning strategy can be confidently assumed uniform.Such patch is obtained by considering only the 10% of the pixels with the highest values in the hit counts map provided by the SO Collaboration.The adopted SO-SATs sensitivities, which are reported in Ade et al. (2019), indeed refer to a homogeneous hit counts map with such a sky fraction.The final noise map, for both experiments, is obtained by linearly combining the central frequency channels noise maps with an inverse noise weighting.We considered the frequency range 93 ≤  ≤ 225 for SO-SATs and 89 ≤  ≤ 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  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 {, } 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  f .This signal is generated with the PySM Python package6 (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 as s1 and d1.The synchrotron Spectral Energy Distribution (SED) follows a power law: where  s is the spectral index,  the position in the sky, and  the considered frequency.{, } s (,  0 ) represents the synchrotron template at a reference frequency  0 and is the WMAP 9-year 23-GHz / 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 black-body (MBB): where   () is the black-body spectrum,  d is the dust spectral index,  d is the dust temperature, and {, }  (,  0 ) represents the dust template at a reference frequency  0 .Both the spectral parameters maps and the template are obtained from the application of the Generalised Needlet Internal Linear Combination (GNILC) to Planck data (Planck Collaboration 2020a) and non-Gaussian smallscale 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, which, 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 MFs detection power for these future CMB experiments.If one wants to expand the assessment of MFs sensitivity at smaller scales, they would have to choose a signal model different than 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 foregrounds maps are the same of LiteBIRD and SO-SATs, while, as noise component, we used the end-to-end SMICA noise residuals maps introduced in Section 5.3.All maps are brought to the resolution of 37.8 ′ for the comparison between Planck and LiteBIRD, of 30 ′ for the comparison between Planck and SO-SATs.

Compatibility analysis
In order to assess the level of deviation between MFs on different maps, we use the  2 statistic.Let  be a vector with the result of any of the MFs with the mean of the simulations subtracted and Σ the covariance matrix of  (representing the correlation of MFs at different thresholds).The  2 statistic is given by: If  is the sample covariance (computed with simulations), an unbiased estimator of the inverse of the covariance is: where  = 300 is the number of simulations and  = 28 is the number of bins in the considered range of thresholds .The prefactor in this equation is known as the Hartlap factor (Hartlap et al. 2007) and it is needed for the estimator to be unbiased.
In this analysis we report the reduced χ2 =  2 / dof , where  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  2 : where  2 data and  2 sims are respectively associated to any MF computed on the data or on a simulation.In equation ( 22), a value close to 50% means that the data closely resemble the expectation, while a value close to 0% or 100% means that the data are highly incompatible with simulations.
Finally, we also report the significance of the deviation Δ  of  2 in the data with respect to the mean χ2 sims among the simulations, reported in units of standard deviation  sims (see last column of Table 1): To assess the  2 MFs detection power of a non-Gaussian anisotropic signal in the simulated data sets introduced in Section 5.5, we compute a  2 for each realisation as in equation ( 20), where  is the MFs difference between the test and the fiducial Gaussian data set for that realisation and the covariance matrix is estimated from MFs of the fiducial data set.We then report for each case the average probability for a  2 distribution with  dof = 28 (the adopted number of thresholds) to exceed the computed  2 s.If this value is close to 50%, it means that no deviations are observed between MFs of the test and Gaussian fiducial data set.Conversely, if such probability is very close to 0%, we are able to detect at high significance the injected non-Gaussian anisotropic signal in the test data set.

RESULTS
In this section, we proceed to validate the code and the theoretical formulae for the expected values of the MFs for  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  2 maps of LiteBIRD and SO.We use 300 simulations for each case.Throughout all the following analysis, we set thresholds between  = 0 and 5.For visualisation purposes, all MFs are plotted in 100 bins, which correspond to a spacing of Δ = 0.05.The statistical analysis described in Section 5.6 is instead performed the number of bins to 28 (Δ = 0.18) so 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  2 =  2 +  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  2 data sets of LiteBIRD and SO and the results are compared with those obtained from analogous Planck simulations.

Validation
We verify the theoretical predictions for the MFs of  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-fit angular power spectra, so  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 and the theory well inside the 1 region, where  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 ( mean =  sims √ 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 (  ())/  () at every threshold  is below 1% and, in the case of  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, ℓ  (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  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.

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 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 MFs 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 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 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 ±2.Only  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.We quantify the compatibility of all curves in terms of the  2 statistic (as well as the probability-to-exceed and the deviation with respect to simulations in terms of Δ  ), 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  0 seem to present a pattern at high values of the threshold for both SMICA and SEVEM solutions.From the  2 analysis, this trend is not significant and can be explained by the fact that the covariance matrix of  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  sky = 40%), thus disfavouring the hypothesis of a residual Galactic contamination origin.
So far, we have focused on Planck  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  ±  SMICA and SEVEM Planck maps and corresponding end-to-end simulations and then construct needlet  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 fullrange 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  modes in Planck Collaboration (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% of the sky along the Galactic plane before applying the needlet transform to avoid leaking of a residual foregrounds component outside the Planck common mask in the first needlet scales.We used the publicly released Planck GAL97 Galactic mask 7 .
In conclusion, it is possible to state that, within the MFs sensitivity in  2 , no unexpected deviation from Gaussianity and statistical isotropy has been detected in the Planck CMB  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.

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 2020b).Conversely, in polarization, only modes at large angular scales have a modest signal-tonoise 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  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  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  2 maps to detect such a signal with varying amplitude in simulated LiteBIRD and SO-SATs 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  2 MFs at different angular scales.For each experiment, the MFs are computed on needlet maps of the fiducial and test data set.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  2 distribution with  dof = 28 (the adopted number of thresholds) to exceed the computed  2 s.

LiteBIRD vs 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  2 simulations of the fiducial and test data set 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  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  modes.The angular power spectra of such rescaled non-Gaussian anisotropic foreground  2 signal is compared with that 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  f = 0.01, the non-Gaussian signal is so low that  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  f = 0.05, LiteBIRD yields a highly significant detection at needlet scale  = 1 with all MFs and a possible detection at the largest angular scales (  = 0) in  1 and  2 .No deviations, instead, are observed either at higher needlet scales or in Planck.When  f = 0.1, we find a clear detection for LiteBIRD at  = 0 in  1 and  2 , at  = 1 in all MFs, and hints for a detection at  = 2 in  1 and  2 , while, for Planck, deviations only at  = 0 in  1 and  2 .Finally, with  f = 0.3, for LiteBIRD we have a detection at  = 0, 1, 2, and 3, while for Planck we can observe a detection at  = 0 and  = 1 and hints for a deviation at  = 2.
The reported results highlight the highly enhanced detection power of MFs if applied to LiteBIRD  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 which induce deviations from Gaussianity and statistical isotropy in the polarization data.We can note that, in general,  1 and  2 can better detect deviations than  0 .Moreover, for all the considered cases, unexpectedly we have more significant detections for  = 1 rather than for  = 0.This is motivated by the fact that the adopted ] . Angular power spectra of  2 maps computed from input  and  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 the 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 ′ .Right panel: Power spectra have been computed within the reduced SO-SATs 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 ′ .The adopted binning scheme is Δℓ = 15.
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  = 0 highly increases in this case.

SO-SATs vs Planck
SO-SATs is one of the most promising near-future experiments for the observation of primordial  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  2 simulations of the fiducial and test data set for both SO-SATs and Planck within the reduced SO-SATs 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 ℓ > 30, we consider only needlet scales with  ≥ 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 (  sky = 78%).Therefore we adopt the following values for  f : {0.1, 0.3, 0.5, 0.7}.The angular power spectra of such rescaled non-Gaussian anisotropic foreground  2 signal is compared with that of simulated CMB maps in the right panel of Fig. 7.All angular power spectra are computed within the reduced SO-SATs 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  f = 0.1,  2 MFs are not sensitive to the injected non-Gaussian anisotropic signal in the SO-SATs patch neither for SO nor for Planck.Instead, by computing  1 and  2 in SO-SATs data, we can detect deviations at  = 2 for  f = 0.3, at  = 2, 3, and 4 for  f = 0.5 and  f = 0.7, and hints for deviations at  = 3 for  f = 0.3.For the two largest considered rescaling amplitudes ( f = 0.5 and 0.7), also  1 and  2 of Planck needlet  2 test data set at  = 2 significantly deviate from the fiducial Gaussian case.Also for SO-SATs the improved sensitivity will highly enhance the detection power of MFs in  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  0 results less sensitive to non-Gaussian anisotropic signals than  1 and  2 .In this analysis we assumed a pessimistic 1/f noise component for SO-SATs.We thus expect even a higher detection power in the first needlet bands in case of improved instrumental properties.

CONCLUSIONS
MFs have been widely applied to the analysis of CMB data (Schmalzing & Górski 1998; Planck Collaboration 2020d), but their application to polarization has been limited to date.In this work, we introduce a generalised framework to compute the MFs on maps of the squared modulus of polarization ( 2 =  2 +  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), 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  map reflects the local morphological properties of the polarization field.Therefore, we expect the introduced formalism to better capture those deviations which 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: 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  2 maps, such an analysis allows us to assess mainly the statistical properties of noise.
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), Simons Observatory (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  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  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;Schwarz et al. 2016;Adhikari et al. 2016).
The introduced formalism may also be useful to study other complex spin fields defined on the sphere, such as the convergence () or the gravitational shear () maps associated to 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 et al. 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 (Jung et al. 2018;Ben-David et al. 2015).
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.

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 scanning strategy (Planck Collaboration 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  2 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.

Figure 1 .
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    = 78%.

Figure 3 .
Figure3.Representation in the harmonic domain of the needlet configuration adopted in this analysis: standard needlets with  = 2.The first four bands have been merged together according to equation (14), while the other bands are left unchanged.

Figure 4 .
Figure 4.In light grey: The SO-SATs full footprint with sky fraction    = 34%.In dark grey: The SO-SATs 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    = 10%.

Figure 5 .
Figure5.Average MFs of  2 =  2 +  2 for  sims = 300 CMB Gaussian isotropic simulations (orange dots) are compared with the theoretical predictions (solid black lines).From top to bottom:  0 ,  1 , and  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, given by the grey contours) and the dispersion of the mean of simulations (computed as / sims ) in the middle and bottom sections of each panel, respectively.See text for details.

Figure 6 .
Figure6.MFs of the Planck SMICA released  2 =  2 +  2 maps (blue dots) are compared with the mean given by realistic end-to-end simulations (pink solid line).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 area represents the ±1 region obtained in the simulations.
port from ASI/INFN grant n. 2021-43-HH.0.Part of this work was also supported by the InDark INFN project.DM acknowledges support from the MIUR Excellence Project awarded to the Department of Mathematics, Università di Roma Tor Vergata, CUP E83C18000100006.DM is also grateful to the Department of Excellence Programme MatModTov for support.
Collaboration et al. 2023on 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 (LiteBIRDCollaboration et al. 2023 LiteBIRD and SO simulations CMB intensity anisotropies have been measured by Planck with very high accuracy (Planck Collaboration 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, 5 https://javicarron.github.io/mtneedletare planned

Table 1 .
Deviations of MFs computed on Planck SMICA and SEVEM CMB  2 maps with respect to those of the corresponding end-to-end simulations, in terms of normalised  2 , probability-to-exceed, and Δ  of Eq. 23.

Table 2 .
Deviations of MFs computed on Planck SMICA and SEVEM CMB needlet  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.

Table 3 .
Deviations of MFs computed on needlet  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(  sky ∼ 78%).The deviations are expressed in terms of the average probability in percentage for a  2 distribution with  dof = 28 (the adopted number of thresholds) to exceed the computed  2 s.Values lower than 10 −6 are approximated with 0.0.Values close to 50% mean no detection of the injected non-Gaussian anisotropic signal, while values very close to 0% mean detection at very high significance.Different magnitudes of the overall rescaling factor  f of the residual foreground contamination are considered.

Table 4 .
Deviations of MFs computed on needlet  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-SATs patch shown in Fig.4and excluding the masked regions of the Planck polarization common mask (  sky ∼ 10%).The deviations are expressed in terms of the average probability in percentage for a  2 distribution with  dof = 28 (the adopted number of thresholds) to exceed the computed  2 s.Values lower than 10 −6 are approximated with 0.0.Values close to 50% mean no detection of the injected non-Gaussian anisotropic signal, while values very close to 0% mean detection at very high significance.Different magnitudes of the overall amplitude  f of the residual foreground contamination are considered.