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

PAPER • OPEN ACCESS

# Morphological analysis of the polarized synchrotron emission with WMAP and Planck

To cite this article: F.A. Martire et al JCAP04(2023)049

View the [article online](#) for updates and enhancements.

You may also like

- The valve filament at constant voltage
  E H W Banner

- Adaptive Optics Sky Coverage
  Calculations for the GeminiNorth
  Telescope
  B. L. Ellerbroek and D. W. Tyler

- Production of Filaments by Surges
  Yu Liu, Hiroki Kurokawa and Kazunari
  Shibata
# Morphological analysis of the polarized synchrotron emission with WMAP and Planck

F.A. Martire,ª A.J. Banday,ᵇ E. Martínez-Gonzálezª and R.B. Barreiroª

Received January 20, 2023

Revised March 17, 2023

Accepted March 25, 2023

Published April 19, 2023

**Abstract.** The bright polarized synchrotron emission, away from the Galactic plane, originates mostly from filamentary structures. We implement a filament finder algorithm which allows the detection of bright elongated structures in polarized intensity maps. We analyse the sky at 23 and 30 GHz as observed respectively by *WMAP* and *Planck*. We identify 19 filaments, 13 of which have been previously observed. For each filament, we study the polarization fraction, finding values typically larger than for the areas outside the filaments, excluding the Galactic plane, and a fraction of about 30% is reached in two filaments. We study the polarization spectral indices of the filaments, and find a spectral index consistent with the values found in previous analysis (about −3.1) for more diffuse regions. Decomposing the polarization signals into the *E* and *B* families, we find that most of the filaments are detected in *P*<sub>*E*</sub>, but not in *P*<sub>*B*</sub>. We then focus on understanding the statistical properties of the diffuse regions of the synchrotron emission at 23 GHz. Using Minkowski functionals and tensors, we analyse the non-Gaussianity and statistical isotropy of the polarized intensity maps. For a sky coverage corresponding to 80% of the fainter emission, and on scales smaller than 6 degrees (ℓ > 30), the deviations from Gaussianity and isotropy are significantly higher than 3σ. The level of deviation decreases for smaller scales, however, it remains significantly high for the lowest analysed scale (~1.5°). When 60% sky coverage is analysed, we find that the deviations never exceed 3σ. Finally, we present a simple data-driven model to generate non-Gaussian and anisotropic simulations of the synchrotron polarized emission. The simulations are fitted in order to match the spectral and statistical properties of the faintest 80% sky coverage of the data maps.

**Keywords:** CMBR experiments, CMBR polarisation

**ArXiv ePrint:** 2301.08041

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

<footnote>ªInstituto de Física de Cantabria, CSIC-Universidad de Cantabria, Avda. de los Castros s/n, E-39005 Santander, Spain</footnote>

<footnote>ᵇIRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France
E-mail: martire@ifca.unican.es, anthony.banday@irap.omp.eu, martinez@ifca.unican.es, barreiro@ifca.unican.es</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Polarized intensity</td><td>2</td></tr><tr><td>2.1</td><td>Data</td><td>2</td></tr><tr><td>2.2</td><td>Smoothing</td><td>3</td></tr><tr><td>2.3</td><td>Debiased estimator</td><td>3</td></tr><tr><td>3</td><td>Filament finder</td><td>4</td></tr><tr><td>3.1</td><td>Algorithm</td><td>4</td></tr><tr><td>3.2</td><td>Band-pass filter</td><td>5</td></tr><tr><td>3.3</td><td>Results</td><td>5</td></tr><tr><td>4</td><td>Filament analysis</td><td>6</td></tr><tr><td>4.1</td><td>Filaments</td><td>6</td></tr><tr><td>4.2</td><td>Polarization fractions</td><td>7</td></tr><tr><td>4.3</td><td>Spectral index</td><td>8</td></tr><tr><td>4.4</td><td>E and B modes</td><td>10</td></tr><tr><td>4.5</td><td>408 MHz Haslam map</td><td>11</td></tr><tr><td>5</td><td>Statistical properties</td><td>13</td></tr><tr><td>5.1</td><td>Minkowski formalism</td><td>13</td></tr><tr><td>5.2</td><td>Masking and filtering</td><td>14</td></tr><tr><td>5.3</td><td>Gaussian simulations</td><td>15</td></tr><tr><td>5.4</td><td>Results</td><td>16</td></tr><tr><td>6</td><td>Non-Gaussian simulations</td><td>16</td></tr><tr><td>6.1</td><td>A model for non-Gaussian emission</td><td>16</td></tr><tr><td>6.2</td><td>Model for anisotropic modulation of the emission</td><td>19</td></tr><tr><td>6.3</td><td>Tuning and results</td><td>19</td></tr><tr><td>7</td><td>Summary and conclusions</td><td>20</td></tr><tr><td>A</td><td>Power spectra</td><td>24</td></tr><tr><td>B</td><td>Finder algorithm</td><td>25</td></tr><tr><td>B.1</td><td>Toy model for filaments</td><td>25</td></tr><tr><td>B.2</td><td>Minimal length criteria</td><td>27</td></tr><tr><td>C</td><td>Planck statistical properties</td><td>28</td></tr><tr><td>C.1</td><td>Gaussian simulations</td><td>28</td></tr><tr><td>C.2</td><td>Non-Gaussian simulations</td><td>28</td></tr></table>
# 1 Introduction

The principle focus of studies of the Cosmic Microwave Background (CMB) in recent years has been, and indeed remains, to detect, and subsequently characterize, its polarized emission. Precision measurements are necessary in order to detect the very weak B modes generated by primordial gravitational waves, as specifically predicted by models of inflation. However, the sky emission at radio/microwave wavelengths also contains various foreground sources of astrophysical emission that completely obscure the cosmological signal. Specifically, at frequencies below a few GHz, the emission is dominated by synchrotron and free-free radiation; above 10 GHz the so-called Anomalous Microwave Emission (AME) contribution becomes significant; and above 70 GHz the thermal dust emission becomes dominant. Synchrotron and dust radiations are highly polarized, the free-free radiation is intrinsically unpolarized, and AME is expected to be only weakly polarized [1].

Synchrotron radiation is due to relativistic cosmic ray (CR) electrons accelerating around the Galactic magnetic field. The emission intensity depends on the density and energy distributions of the electrons, and on the Galactic magnetic field strength. The electron energy distribution can be approximated by a power-law, $N(E) \propto E^{-p}$ [2, 3], with typical values close to $p = 3$ [4]. In a uniform magnetic field, the synchrotron radiation can reach a polarization fractional of $\Pi = (p+1)/(p+7/3) \approx 75\%$. However, we observe much smaller values, because of a geometric depolarization due to tangled magnetic fields and superposition effects along the line-of-sight. In addition, at frequencies typically below 10 GHz, Faraday rotation effects can depolarize the synchrotron emission near the Galactic plane [5].

The synchrotron spectral index is related to $p$ via $\beta = -(p+3)/2$. Typical values around -3 have been reported for spectral index both in intensity [6, 7] and polarization [8–10]. More recent analyses seem to indicate the presence of spatial variations in the spectral index [11–13], with a tendency to steeper values moving from low to high Galactic latitudes [5].

The all-sky map at 408 MHz [14–16], often referred to hereafter as the Haslam map, provides the best full-sky representation of the synchrotron intensity emission. This is mainly due to cosmic rays accelerated by shock fronts in supernova remnants (SNRs) and pulsar wind nebulae (PWN). Outside the Galactic plane, the strong emission originates mostly from filamentary structures. The North Polar Spur (NPS), or Loop I, is the most obvious feature, but others have been observed: the Cetus arc (or Loop II) [17], Loop III [18] and Loop IV [19]. Those filaments are even more visible in the polarized sky. [20] identify and study 11 filaments in the WMAP polarization maps [21]. The true origins of filaments are still poorly understood. The most widely accepted progenitors of these large structures are old and nearby supernova remnants [22].

The presence of complex structures such as loops and filaments makes the statistics of the synchrotron emission strongly non-Gaussian and anisotropic at large scales, even in the diffuse region. However, it is reasonable to suppose that at small scales the emission could approach Gaussianity and isotropy as a manifestation of the central limit theorem. Several models used to simulate synchrotron assume that the small scale fluctuations are statistically isotropic and Gaussian [23–26]. Analyzing the 408 MHz map, [27] showed that the level of the non-Gaussian deviations decreases on smaller scales as expected, but remains significantly high ($> 3\sigma$) on angular scales of $\sim 1.5^\circ$. These results were confirmed by analyzing WMAP and Planck intensity maps [28].

However, little is known in polarization. Knowing the morphological and statistical properties of the polarized foregrounds emission is crucial to face the future challenges detecting
cosmological signals. Some component separation methods used to produce CMB maps require prior knowledge of the foregrounds [8, 29]. Moreover, the increasing sensitivity of on-going and future experiments, such as the Simons Observatory [30] and LiteBIRD [31], requires more realistic foreground models and simulations. From another perspective, characterizing the synchrotron emission can give us useful information for understanding the physical mechanisms behind the Galactic magnetic field [20].

In this work, we characterize some morphological and statistical features of the synchrotron polarization. We analyze the observations of the *WMAP* K-band and *Planck* 30 GHz frequency channels in a region of the sky where the emission is predominantly diffuse. Section 2 describes the data set and simulations used for the analysis. In section 3, we present a filament finder algorithm and demonstrate its performance on data maps. Section 4 contains the description and analysis of the polarization fractions, spectral indices, *E* and *B* nature, and possible intensity counterparts of the detected filaments. In section 5, we test the statistical properties of the *WMAP* polarization maps. Section 6 presents a simple model to generate simulations which better resemble the statistical nature of the polarized synchrotron. We summarize our results and provide discussion about their implications in section 7. Finally, in appendix A we give a brief review of some cosmological quantities, in appendix B we test the accuracy and limits of the finder algorithm, and in appendix C we present the results obtained from the *Planck* statistical analysis.

# 2 Polarized intensity

## 2.1 Data

For our analysis, we will make use of data taken by the *WMAP* and *Planck* satellites. We focus on the lowest frequency data from *WMAP*, specifically the 9-year *WMAP* K-band (centered at 23 GHz) maps, provided in the *HEALPix*<fnref n="1" /> pixelisation scheme with $N_{side} = 512$ and an effective Gaussian beam of 0.88° full-width-at-half-maximum (FWHM). The *WMAP* products have been downloaded from the Legacy Archive for Microwave Background Data Analysis (LAMBDA).<fnref n="2" /> For the *Planck* analysis, we use the 30 GHz frequency maps generated by the *NPIPE* processing pipeline (PR4). The *NPIPE* processing results in improved High Frequency Instrument (HFI) polarization data with reduced systematic artefacts and lower levels of noise. PR4 data from the Low Frequency Instrument (LFI) are also modified with respect to the 2018 *Planck* release. Further details are available in [32]. The frequency maps were downloaded from the *Planck* Legacy Archive<fnref n="3" /> (PLA) at a pixel resolution corresponding to $N_{side} = 1024$ and an effective beam of FWHM = 31.5 arcminutes. Note that the polarization maps, and consequentially the analysis, follow the *HEALPix* convention.

As the synchrotron emission scales with frequency, the foreground signal is higher in the *WMAP* K-band compared to the *Planck* 30 GHz channel, however, the noise level of *Planck* is lower. As result, at a common resolution of 1°, the overall signal-to-noise ratio of the two experiments are similar [33], although one or other map may be better in some sky regions because of the different scanning strategies.

<footnote n="1">¹https://healpix.sourceforge.io.</footnote>

<footnote n="2">²https://lambda.gsfc.nasa.gov/product/map.</footnote>

<footnote n="3">³https://pla.esac.esa.int.</footnote>
## 2.2 Smoothing

We smooth the WMAP and Planck maps to common resolutions of 1° and 3° FWHM. The 1° maps are downgraded to a HEALPix resolution of $N_{side} = 128$ (corresponding to a representative pixel size of ~27 arcmin) that we will use for filament detection and statistical analysis, and the 3° maps to a resolution of $N_{side} = 64$ (a pixel size of ~55 arcmin), that we will use for the filament analysis. The smoothing and downgrading are performed in harmonic space deconvolving the original effective beam and then convolving with a Gaussian beam.<fnref n="4" /> The smoothing process helps to increase the signal-to-noise ratio of the maps and to minimise any effect due to beam-asymmetries in the two experiments.

We estimate the noise level of the data at the 1° and 3° resolutions, including uncertainties due to smoothing and pixel downgrading. For WMAP, we generate 600 Gaussian noise realisations based on the covariance matrices at full resolution. For Planck, we use the 600 noise simulations provided on the PLA [32]. We downgrade and smooth each simulation in the same way as the data. Finally, for each pixel we compute $\sigma_Q^2$, $\sigma_U^2$ and $\sigma_{QU}^2$ from the variance and covariance over all of the simulated Q and U maps.

## 2.3 Debiased estimator

A morphological analysis is applied to the polarized intensity $P = \sqrt{(Q^2 + U^2)}$, which, given its positive nature, is subject to noise bias. In particular, in the low signal-to-noise regime, $P$ will yield a positive estimate even if Q and U are zero. We use the modified asymptotic (MAS) estimator [34] in order to correct the polarized amplitude for the bias. We recall that the debiased polarized amplitude with the MAS estimator is computed as

$$P_{\text{MAS}} = P - \frac{1 - \exp(-P^2/b^2)}{2P} b^2 \quad (2.1)$$

where the noise bias $b$ is function of the pixel variance and the polarization angle $\phi = \arctan(U/Q)$ given by

$$b^2 = \sigma_U^2 \cos^2(\phi - \theta) + \sigma_Q^2 \sin^2(\phi - \theta), \quad \theta = \frac{1}{2} \arctan \left( \frac{\sigma_{QU}}{\sigma_Q^2 - \sigma_U^2} \right). \quad (2.2)$$

The maps are shown in figure 1. [34] demonstrate that in the regime where the signal-to-noise ratio exceeds 2, the estimator is unbiased and essentially Gaussian. An estimate of the variance is then given by

$$\sigma_P^2 = \sigma_Q^2 \cos^2(\phi - \theta) + \sigma_U^2 \sin^2(\phi - \theta). \quad (2.3)$$

Note that the debiased polarized intensity and its variance are defined pixel by pixel. Thus, we do not take into account the correlation between pixels, which is introduced smoothing and downgrading of the maps. Although we use the MAS estimator in our analysis, several tests have been performed using the Wardle & Kronberg estimator [20, 35, 36] instead, finding consistent results.

<footnote n="4">⁴We generate $a_{\ell m}$'s with the `map2alm` healpy routine from the Q, U maps. We convolve the maps with the new Gaussian beam and pixel window function following the method described in appendix A. Finally, we regenerate the Q, U maps with the `alm2map` routine from the convolved $a_{\ell m}$'s.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Debiased Polarized Intensity maps (1°) of the WMAP K-band at 23 GHz (left) and the 30 GHz channel of the Planck PR4 data set (right).</figcaption></figure>

# 3 Filament finder

## 3.1 Algorithm

We develop an algorithm in order to seek elongated structures in polarized intensity maps. The method resembles the two-dimensional version of the Smoothed Hessian Major Axis Filament Finder (SHMAFF) [37, 38]. The SHMAFF algorithm has previously been used to find filaments in the three-dimensional galaxy distribution and in the analysis of polarized dust structures in the Planck 353 GHz maps [39]. The main difference in our implementation is that, while the SHMAFF detection is based on the minimal eigenvalues of the Hessian matrix, our detection is based directly on the polarized intensity. This is because the Hessian matrix fails to find pixels with minimal eigenvalues in detecting the thick and diffuse filaments which we expect to be in the noisier area of the sky, that is outside the Galactic plane.

The algorithm works on a pixel by pixel basis, examining the orientation angle $\psi$ defined in [39] and comparing the polarized intensity value $P$ with respect to a threshold $P_{th}$. The orientation angle is determined from

$$ \psi = \arctan\left(-\frac{H_{\theta\theta} - H_{\phi\phi} + \alpha}{2H_{\theta\phi}}\right) \qquad (3.1) \\ \text{with } \alpha = \sqrt{(H_{\theta\theta} - H_{\phi\phi})^2 + 4H_{\theta\phi}^2} $$

where $H$ is the Hessian matrix computed from the second-order covariant-derivatives with respect to the spherical coordinates $(\theta, \phi)$ [40]. We compute the threshold values $P_{th}$ from the $P$ distribution. We cannot use the mean and the standard deviation because the distribution of $P$ is not Gaussian and exhibits an extended tail, thus, we define the threshold from the median $m_P$ and the median absolute deviation (MAD) $\sigma_m$ [41, 42] as

$$ P_{th} = m_P + \sigma_m = m_P + 1.4826 \cdot \text{median}(|P - m_P|). \qquad (3.2) $$

The algorithm starts by identifying the brightest pixel $P_0$ and denoting its orientation angle $\psi_0$. It then considers its 8 (or 7) neighbouring pixels, identified with the `get_all_neighbours` routine of HEALPix. For each neighbour pixel, two conditions are checked: (i) if its polarized intensity is larger than the fixed threshold of equation (3.2), (ii)
if its orientation angle is coherent with the initial pixel

$$ (i) \quad P_i > P_{\text{th}} \qquad (3.3) \\ (ii) \quad |\psi_i - \psi_0| < \Delta\psi, $$

where we fix $\Delta\psi$ at 10° in order to detect the longest filaments. Nevertheless, we have checked that the results are robust for a reasonable range of $\Delta\psi$. If both conditions in equation (3.3) are satisfied, we accept the pixel as part of the same structure of the initial pixel. Then, we look for its neighbours which also satisfy the two conditions. We continue this friend-of-friend recursive algorithm until the conditions are no longer satisfied. Once a structure is defined, the pixels are masked from the map, the new brightest pixel is identified and the condition-based procedure is repeated.

We finally define a template, which includes all the pixels satisfying the previous conditions, with all the strongly polarized areas in the $P$ map of arbitrary shapes and sizes. We then smooth it with a 3° Gaussian beam, in order to soften the boundaries. In order to allow only elongated structures which can be identified as filaments, we apply one further criteria to the 3° smoothed map, i.e. we reject structures with length smaller than a threshold minimal length $L < L_{\text{th}}$, fixed at 10°. The length is defined as the maximum pixel-pair angular distance. The threshold length has been selected analysing simulations as discussed in appendix B.2.

## 3.2 Band-pass filter

Different filaments have been observed in low frequency (< 1 GHz) radio continuum surveys, and more recently in the WMAP polarization data [20]. The filaments can be divided into two categories: bright and narrow, or weak and diffuse. In order to optimize the detection of these filaments, we filter the maps to focus on specific angular scales of interest before applying the filament finder algorithm. We compute the spherical harmonic coefficients $a_{lm}$ by means of the map2alm routine of HEALPix, multiply by a band-pass filter of the following form

$$ f(l) = \frac{1}{4} \left[ 1 + \tanh \left( \frac{l - l_{\min}}{\Delta l} \right) \right] \left[ 1 - \tanh \left( \frac{l - l_{\max}}{\Delta l} \right) \right] \qquad (3.4) $$

then generate the filtered maps with the alm2map routine. The filter cuts off the amplitudes below a multipole scale $l_{\min}$ and above $l_{\max}$. In order to detect the thinner filaments, we consider multipoles in the range 20-50, whereas for the diffuse filaments the multipoles are restricted to the range 15-20. The cuts are roughly in accordance with the widths of the filaments ($l \sim 180/\theta$). $\Delta l$ is set to 10, but the method is robust for a reasonable range of $\Delta l$.

The filters are shown in figure 13 in appendix 2. By filtering out small-scale modes, we enhance the contrast of larger structures with respect to the diffuse foreground emission, and also reducing the instrumental noise. Moreover, we remove correlations on large scales which can negatively affect the detection. The application of the filter is critical to increase the accuracy of the estimation of the polarization orientations of the filaments, especially in areas where the signal-to-noise is low. A similar filter has been used in [27] for a statistical analysis of the 408 MHz Haslam data.

## 3.3 Results

Starting from the debiased polarized intensity maps obtained as described in section 2.3, we generate two bandpass-filtered maps $P^{20-50}$ and $P^{15-20}$, applying the filter in equation (3.4),
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Filamentary structures detected in the debiased polarized intensity of the *WMAP* K-band at 23 GHz (left) and the *Planck* 30 GHz channel (right).</figcaption></figure>

where the superscripts correspond to the applied multipole ranges. We first mask the bright point sources, both Galactic and extragalactic, which could bias our algorithm. We use the mask derived for the *Planck* PR4 SEVEM component separation pipeline which includes all the point sources that have polarization detection significance levels of 99% or more in the 30 GHz polarized map [32, 43, 44]. The finder algorithm is based on the threshold condition in equation (3.3), where the threshold is computed from the *P* distribution (see equation (3.2)). However, very bright areas such as the Galactic plane could bias the threshold value upwards, preventing filament detection. Therefore, when analysing the *P*<sup>20-50</sup> maps, we apply a Galactic mask excluding pixels at latitudes |*b*| < 3°. A similar argument applies to the *P*<sup>15-20</sup> analysis. However, since we are looking for very faint filaments, the NPS and the Southern Fan regions are also excluded in addition to the Galactic plane. This can be achieved by simply masking the brightest 30% sky fraction. The specific choice of masks was tested on the simulations described in the appendix B.1.

After applying the filament finder algorithm to the two bandpass-filtered maps, we merge the two sets of results into one template, smooth it by 3° and apply the criteria of minimal length as described in section 3.1. The filament templates determined independently from the *WMAP* K-band and *Planck* 30 GHz data are shown in figure 2. The *WMAP* results reveal more compact and elongated structures than for *Planck*, indicating that the finder algorithm performs better when applied to the data with brighter synchrotron emission. However, several similar structures are detected in the same areas of the sky in both maps. The agreement between the independent results corroborates the validity of the algorithm and supports the existence of the filaments as real emission, and not due to noise or systematic effects. In appendix B.1, we explore the accuracy and limits of the filament finder algorithm, testing our method with toy filamentary foreground models.

# 4 Filament analysis

## 4.1 Filaments

Combining the information obtained from *WMAP* and *Planck*, we get a final template of polarized filamentary structures detected in the frequency range 23–30 GHz, shown in figure 3. We use the *WMAP* detection as our benchmark, but only retain those structures which are detected, at least in part, in the *Planck* data. The one exception is filament XI, which is clearly detected in the *Planck* data, but only partially in *WMAP*. Its existence is supported by previous analysis performed on *WMAP* [20]. This method ensures that detected filaments
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Template map showing the filaments of the polarized synchrotron emission detected by combining *WMAP* (23 GHz) and *Planck* (30 GHz) results.</figcaption></figure>

are not due to noise or systematic residuals. Then, we remove those structures which do not show clear elongation, in particular in the Galactic plane, where the emission is more complex, and polar regions, where the signal-to-noise is low.

In order to specify the filaments, we use and expand the nomenclature used in [20]. Filaments I (NPS), IIIn, and IV have been recognised and studied for more than 60 years. These large structures have been observed in X-ray, gamma-ray and other microwave experiments [18, 19, 45]. Filament II (Cetus Arc) was previously detected in the radio sky [17], and found here for the first time in polarization despite its low emission. We detect ten further filaments reported in [20] (Is, GCS, IIIs, VII, IX, X, XI, XII, XIII, XIV) but not filament VIII. We also identify five new filaments (XV, XVI, XVII, XVIII, XIX) that are visible in both *WMAP* and *Planck*. Filament XV is a bright structure at the center of the Galactic Haze [46]. Filament XVI, because of the position and shape, seems to be related to Filament I. Filament XVIII is a bright structure of the Northern Fan region close to the Galactic plane. Finally, Filaments XVII and XIX are new detections in the region below the Galactic plane.

Most of the detected filaments have circular arc-like shapes, supporting the model of supernova remnants expanding into the Galactic magnetic field [22]. Several structures appear to be spatially correlated with each other, although most are stand-alone features. In the NPS, there are several elongated structures which do not resemble loop-like features. They were first identified in radio observations [19], but the NPS complexity is more evident in polarization. Because of their location, there are models which link these structures with the Fermi Bubbles (FB) detected in Fermi data at energies ~10–500 GeV [47, 48].

## 4.2 Polarization fractions

The polarization fraction is defined as the ratio of the polarization amplitude to total intensity

$$ \Pi = \frac{P}{I}. \qquad (4.1) $$
At the low frequencies of interest here, the synchrotron emission largely dominates the polarization signal, thus, we can neglect other physical emission mechanisms and consider directly the frequency data. However, this assumption does not apply in intensity where the CMB, free-free and AME also contribute to the total emission.

Most of the observations suggest that at frequencies above 20 GHz the spectral index of the synchrotron intensity spectrum is $\beta \approx -3$ [33]. However, according to some models [49], it can get much flatter for frequencies below 10 GHz. Moreover, we expect it to exhibit significant spatial variations. In this analysis, we use the diffuse synchrotron intensity map provided in the Planck 2015 [29] release. The template has been generated at a reference frequency of 408 MHz by Commander (a parametric component separation method) applied to the WMAP, Planck and Haslam observations.

We estimate the synchrotron intensity extrapolating the 408 MHz map up to 23 GHz and 30 GHz adopting a fixed spectral index $\beta = -3.0$. The maps are analysed at 3° resolution and $N_{side} = 64$, masking the Galactic plane.<fnref n="5" /> Pixels with a signal-to-noise ratio lower than 2.5 have been excluded.

The polarization fraction maps are showed in figure 4 (top panels). For each pixel, the polarization fraction error, $\sigma_{II}$, is obtained propagating the errors in the polarization and intensity maps. The intensity uncertainty is, in turn, obtained propagating the uncertainty of the intensity spectral index. The weighted average of the polarization fractions over pixels for each filament are listed in table 1. The largest source of uncertainty in $\sigma_{II}$ is due to the uncertainty in the intensity spectral index. We are aware that, locally, $\beta$ can assume values over a very broad range. However, we compute the polarization fraction averaging over extended areas, so it is reasonable to assume that in these areas $\sigma_{\beta} = 0.1$ as found in previous works on partial-sky analysis [10, 13, 50]. We report good agreement between the WMAP and Planck results, the largest discrepancies arising for the more diffuse filaments, e.g., filament XII.

Assuming a uniform spectral index, we find that the polarization fraction of the filaments are typically larger than for external regions outside the filaments. The filaments with the highest polarization fractions, IX and XV, achieve values above 30%, and are both located in the NPS. Loop I has an average value of about 20%, slightly smaller than the value found at the center of the NPS. These results corroborate the results found in the previous analysis [20, 33]. We report a high polarization fraction also for filament XVIII, located in the Fan region. The lowest polarization fractions are found for filaments XI, XIII and XVII.

## 4.3 Spectral index

The synchrotron spectral energy distribution (SED) is generally approximated by a power law<fnref n="6" /> $S_{\nu} \propto \nu^{\beta}$ where $\beta$ is the energy spectral index. Spatial variations of $\beta$ have been reported in the literature [11–13]. In this section, we measure the spectral index of each filament described in 4.1 employing a method based on the Q and U Stokes parameters [5, 51].

Let us define the vector

$$d(\alpha) = Q \cos(2\alpha) + U \sin(2\alpha) \quad (4.2)$$

<footnote n="5">⁵In order to mask the bright pixels along the Galactic plane, we use a Galactic mask obtained combining the 2015 Galactic plane mask which allows the 90 per cent of the sky (provided in the PLA) and a Galactic latitude mask excluding pixels within $\pm 5^{\circ}$ of the Galactic plane.</footnote>

<footnote n="6">⁶Given the sensitivity of Planck and WMAP data, we cannot explore more complex models.</footnote>
which represents, for each pixel, the projection of the Stokes parameters $(Q, U)$ into a reference frame rotated by the angle $\alpha$. We vary $\alpha$ over the range $(0^\circ, 85^\circ)$, in steps of $5^\circ$. For each filament, we compute a linear fit over all the internal pixels to the relation

$$d_{\text{P30}}(\alpha) = m(\alpha) \cdot d_{\text{WK}}(\alpha) + n(\alpha) \qquad (4.3)$$

where $d_{\text{P30}}(\alpha)$ and $d_{\text{WK}}(\alpha)$ are computed respectively from the Planck and WMAP data. Adding the free parameter $n(\alpha)$ gives the advantage of removing any zero level due to possible systematics in the maps. The fit is performed with the orthogonal distance regression code odr<fnref n="7" /> from SciPy, in order to account for the noise variance of both WMAP and Planck. A calibration error of 0.3% has been added in quadrature to both experiments [21, 32]. From the parameter $m(\alpha)$ and its uncertainty, we compute the spectral index for each $\alpha$ as

$$\beta(\alpha) = \frac{\log m(\alpha)}{\log(\nu_{\text{P30}}/\nu_{\text{WK}})} \qquad \sigma_{\beta}(\alpha) = \frac{\sigma_{m}(\alpha)}{m(\alpha) \log(\nu_{\text{P30}}/\nu_{\text{WK}})} \qquad (4.4)$$

where $\nu_{\text{P30}}=28.4$ GHz and $\nu_{\text{WK}}=22.8$ GHz. The final value of the index is recovered from the weighted average

$$\beta = \frac{\sum_{\alpha=0}^{85} \beta(\alpha) \sigma_{\beta}^{-2}(\alpha)}{\sum_{\alpha=0}^{85} \sigma_{\beta}^{-2}(\alpha)}. \qquad (4.5)$$

Since the $\beta(\alpha)$ values are strongly correlated, we take as the uncertainty on the final spectral index the minimum variance among the measurements $\sigma_{\beta} = \min(\sigma_{\beta}(\alpha))$. We checked that this uncertainty is always larger than the intrinsic uncertainty of $\beta(\alpha)$ given by the standard deviation estimated at different rotation angles.

As a consistency check, we also fitted the spectral index from the debiased polarized intensities with the T-T plot approach. We find consistent results with those determined with the method presented above. It has been shown that the synchrotron spectral index is not stable with respect to polarization orientation in the presence of systematics [52]. The $(Q, U)$ method allows marginalization of the result over the polarization angle, making this approach more reliable.

The maps are analysed at $3^\circ$ resolution and $N_{\text{side}} = 64$, masking the Galactic plane. Pixels with a signal-to-noise ratio lower than 2.5 have been excluded. The maps are converted to Rayleigh-Jeans temperature units and corrected by the colour correction using the coefficients given in [53].

In table 1 the spectral indices determined for the different filaments are listed. For each filament, a mean $\chi^2$ value is computed by averaging over all the values given at different rotation angles by the odr routine. Note that we do not take into account the presence of correlated noise between pixels, thus, leading to an underestimation of the uncertainties. This is one of the reasons which would explain some large $\chi^2$ values. We find spectral index results, both inside and outside the filaments, consistent with the values found in previous analysis of about -3.1 [8-10]. As shown in figure 4 (bottom panel), the $\beta$ values span a very broad range, from -3.59 (XVII) to -2.17 (VII). The value for loop I is consistent with the literature [20, 54]. Filaments IX and XV, the most polarized detections located at the NPS, show slightly flatter values (~ -2.5).

<footnote n="7">⁷https://docs.scipy.org/doc/scipy/reference/odr.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 4. Top: polarization fraction of the WMAP K-band at 23 GHz (left) and the 30 GHz channel of the Planck PR4 (right). Bottom: spectral index $\beta$ (left) with error $\sigma_{\beta}$ (right). Outlines of the filaments (black) and the Galactic mask (green) used in the analysis are also shown. Pixels with a signal-to-noise ratio lower than 2.5 have been masked (grey). For illustrative purpose, the maps are smoothed to a resolution of 5°.</figcaption></figure>

## 4.4 E and B modes

As described in appendix A, the polarized emission can be decomposed into E and B modes. The synchrotron polarized angular power spectra has been analysed in the frequency range 2–30 GHz, finding that: both E- and B- modes can be well described by a power law $C_{\ell}^{EE, BB} \propto \ell^{-2.9}$, the B-to-E ratio ranges between 0.2–0.5 and the EB correlation is compatible with zero [9, 10]. In this section we analyze how the filamentary structure relates to the two polarized components.

In [55] a method was proposed for decomposing the Q and U Stokes parameters into the so-called E- and B-mode families. Starting from the (Q, U) maps, we can compute the $a_{lm}^{E,B}$ coefficients using the `map2alm` routine of HEALPix. Setting $a_{lm}^{B} = 0$ and computing the Stokes parameters with the `alm2map` routine determines the contribution to Q and U from the E mode alone. Similarly, setting $a_{lm}^{E} = 0$, we get the contribution from the B mode. Therefore, we can compute the single-mode polarization intensities as

$$P_E = \sqrt{Q_E^2 + U_E^2} \quad P_B = \sqrt{Q_B^2 + U_B^2} \qquad (4.6)$$

For the sake of brevity, here we show an analysis performed on the WMAP K-band, however, the same conclusions can be obtained from the 30 GHz Planck data. We do not use any estimators to correct the polarized amplitude for the noise bias. This is because the estimator in equation (2.1) would require the decomposition of the variance into the E- and
<table><caption>**Table 1.** Polarization fractions and polarization spectral indices of the WMAP K-band and Planck 30 GHz channel. The synchrotron intensity map is extrapolated from the 408 MHz map up to 23 GHz and 30 GHz using a constant spectral index $\beta = -3.0$. The spectral index is computed over the 23–30 GHz frequency range.</caption><thead><tr><th rowspan="3">filament</th><th rowspan="3">f<sub>sky</sub> [%]</th><th colspan="4">polarization fraction</th><th colspan="3">spectral index</th></tr><tr><th colspan="2">WMAP</th><th colspan="2">Planck</th><th rowspan="2">β</th><th rowspan="2">σ<sub>β</sub></th><th rowspan="2">χ<sup>2</sup><sub>r</sub></th></tr><tr><th>Π [%]</th><th>σ<sub>Π</sub> [%]</th><th>Π [%]</th><th>σ<sub>Π</sub> [%]</th></tr></thead><tbody><tr><td>I</td><td>1.5</td><td>22.6</td><td>9.4</td><td>21.2</td><td>9.3</td><td>-3.06</td><td>0.01</td><td>1.92</td></tr><tr><td>Is</td><td>1.2</td><td>18.3</td><td>7.6</td><td>19.4</td><td>8.5</td><td>-2.80</td><td>0.03</td><td>1.87</td></tr><tr><td>II</td><td>0.8</td><td>24.5</td><td>10.4</td><td>20.4</td><td>9.2</td><td>-3.37</td><td>0.06</td><td>3.01</td></tr><tr><td>IIIn</td><td>2.1</td><td>18.9</td><td>8.0</td><td>17.4</td><td>7.7</td><td>-3.07</td><td>0.02</td><td>1.55</td></tr><tr><td>IIIs</td><td>1.8</td><td>19.8</td><td>8.4</td><td>15.5</td><td>7.1</td><td>-3.65</td><td>0.04</td><td>2.22</td></tr><tr><td>IV</td><td>0.5</td><td>18.1</td><td>7.6</td><td>15.1</td><td>6.8</td><td>-3.07</td><td>0.11</td><td>1.95</td></tr><tr><td>GCS</td><td>0.4</td><td>22.6</td><td>9.4</td><td>22.0</td><td>9.6</td><td>-2.80</td><td>0.03</td><td>2.35</td></tr><tr><td>VII</td><td>0.8</td><td>17.3</td><td>7.3</td><td>20.4</td><td>9.1</td><td>-2.17</td><td>0.09</td><td>1.97</td></tr><tr><td>IX</td><td>0.9</td><td>33.1</td><td>13.7</td><td>32.2</td><td>14.1</td><td>-2.64</td><td>0.02</td><td>1.90</td></tr><tr><td>X</td><td>0.6</td><td>17.1</td><td>7.2</td><td>14.2</td><td>6.3</td><td>-3.39</td><td>0.04</td><td>1.48</td></tr><tr><td>XI</td><td>0.8</td><td>13.2</td><td>5.7</td><td>13.3</td><td>6.0</td><td>-2.45</td><td>0.05</td><td>1.60</td></tr><tr><td>XII</td><td>0.7</td><td>20.6</td><td>9.0</td><td>13.9</td><td>6.5</td><td>-3.38</td><td>0.05</td><td>1.21</td></tr><tr><td>XIII</td><td>0.5</td><td>14.0</td><td>5.8</td><td>12.6</td><td>5.6</td><td>-3.28</td><td>0.03</td><td>2.42</td></tr><tr><td>XIV</td><td>0.5</td><td>17.6</td><td>7.3</td><td>14.5</td><td>6.4</td><td>-3.05</td><td>0.04</td><td>1.34</td></tr><tr><td>XV</td><td>0.5</td><td>33.8</td><td>14.0</td><td>34.8</td><td>15.3</td><td>-2.36</td><td>0.04</td><td>2.44</td></tr><tr><td>XVI</td><td>0.4</td><td>22.0</td><td>9.1</td><td>20.7</td><td>9.1</td><td>-3.39</td><td>0.11</td><td>1.51</td></tr><tr><td>XVII</td><td>0.3</td><td>14.6</td><td>6.1</td><td>15.7</td><td>6.9</td><td>-3.59</td><td>0.05</td><td>1.67</td></tr><tr><td>XVIII</td><td>0.5</td><td>30.1</td><td>12.5</td><td>26.4</td><td>11.6</td><td>-3.08</td><td>0.05</td><td>1.28</td></tr><tr><td>XIX</td><td>0.4</td><td>21.1</td><td>9.1</td><td>26.1</td><td>11.6</td><td>-2.28</td><td>0.09</td><td>0.91</td></tr><tr><td>inside filaments</td><td>15.2</td><td>19.0</td><td>8.0</td><td>17.2</td><td>7.7</td><td>-3.08</td><td>0.01</td><td>3.14</td></tr><tr><td>outside filaments</td><td>62.1</td><td>11.1</td><td>5.0</td><td>10.3</td><td>4.9</td><td>-3.15</td><td>0.01</td><td>3.34</td></tr><tr><td>combined</td><td>76.2</td><td>11.8</td><td>5.2</td><td>11.0</td><td>5.2</td><td>-3.10</td><td>0.01</td><td>3.29</td></tr></tbody></table>

B families, which is not a straightforward operation. Nevertheless, we expect qualitatively correct results for the areas where the signal-to-noise is high, in which the bias is negligible.

Applying the filament finder algorithm to the decomposed maps, we find that nearly all of the filamentary structures are detected at least in part in the $P_E$, but not in the $P_B$ map, as shown in figure 5. As expected, the algorithm fails to detect filaments II, IIIs and X, which are either diffuse or strongly affected by noise bias. The complex filamentary structure of the NPS emission is clearly visible in $P_E$, suggesting its E-nature. In the $P_B$ maps, we detect parts of filaments I and IX (even if slightly shifted). Filament XIX is the only structure which is partially detected in B, but not in E.

## 4.5 408 MHz Haslam map

The current best full-sky map of the synchrotron intensity emission at 408 MHz is due to [14, 15]. A more recent version of the data that has been destriped and cleaned of bright point sources is described in [16]. We apply the filament finder algorithm to this map, provided
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Filamentary structures detected in the decomposed polarized intensity *P*<sub>*E*</sub> (left) and *P*<sub>*B*</sub> (right) maps obtained with the *WMAP* K-band at 23 GHz. In red, the filaments detected from the polarized intensity *P*<sub>*MAS*</sub> as described in section 4.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left: full-sky 408 MHz map. Right: Structures detected in the data (black) compared to the filaments found in polarization at 23–30 GHz (red).</figcaption></figure>

as a LAMBDA product,<fnref n="8" /> to study the filamentary detection in intensity at low frequencies.

The map and the detected filaments are shown in figure 6. Filament I is the major structure in intensity as well as in polarization. Part of the detected structure is well matched by what is found at 23–30 GHz, but is more extended at the lower frequency. A similar observation holds for filament II. The majority of the structures in the NPS observed with *WMAP* are also visible, at least in part, in the *Haslam* map. However, in intensity the synchrotron emission is very diffuse, and the algorithm fails to detect the more diffuse filaments, such as III, XI, XII. An interesting result is that filaments IX and XVI, which are the two strong emission structures detected in polarization around Loop I, are not well detected in the 408 MHz map.

Note that the analysis presented in this section provides us with an additional robustness test for our algorithm. However, only a qualitative comparison between *Haslam* and *WMAP/Planck* maps is possible. This is because, in intensity, the synchrotron emission dominates at 408 MHz, but not at 23/30 GHz, where other components become important. On the other hand, the synchrotron emission dominates in polarization, but comprehensive data sets in polarization do not exist at 408 MHz. We hope that the results presented in this section will stimulate future analysis using new data at similar frequencies.

<footnote n="8">⁸https://lambda.gsfc.nasa.gov/product/foreground/fg_2014_haslam_408_info.html.</footnote>
# 5 Statistical properties

## 5.1 Minkowski formalism

According to Hadwinger's Theorem, any morphological property can be expressed as a linear combination of Minkowski Functionals (MFs). These are defined for any field not requiring any prior assumption, making them particularly advantageous for the analysis of fields for which a non-Gaussian nature is known. Several analyses have already been performed on CMB data using MFs to search for evidence of non-Gaussianity [56–58] and residual foreground contamination [59], or to characterize the properties of foregrounds [27, 60]. In this section, we will briefly review the method used for their numerical calculation, following the methodology developed by [61].

Given a map $u$ on the sphere ($S^2$) and a threshold $v$, there are 3 MFs which represent the area ($V_0$), the perimeter ($V_1$) and the integrated geodesic curvature ($V_2$) of an excursion set, that is the region where $u \ge v$, with boundaries defined by $u = v$. For a map in the HEALPix pixelization, we can numerically compute the MFs via a sum over all pixels

$$V_0(v) = \frac{1}{N_{\text{pix}}} \sum_{\text{pixels}} \mathcal{H}(u-v) \qquad (5.1)$$

$$V_1(v) = \frac{1}{4} \frac{1}{N_{\text{pix}}} \sum_{\text{pixels}} \delta(u-v) \sqrt{u_{;\theta}^2 + u_{;\phi}^2} \qquad (5.2)$$

$$V_2(v) = \frac{1}{2\pi} \frac{1}{N_{\text{pix}}} \sum_{\text{pixels}} \delta(u-v) \frac{u_{;\theta} u_{;\phi} u_{;\theta\phi} - u_{;\theta}^2 u_{;\phi\phi} - u_{;\phi}^2 u_{;\theta\theta}}{u_{;\theta}^2 + u_{;\phi}^2} \qquad (5.3)$$

where $u_{;i}$, $u_{;ij}$ ($i, j \in (\theta, \phi)$) are the first and second partial derivatives in spherical harmonic space, $\mathcal{H}$ is the Heaviside step function and $\delta$ is the delta function. We rescale the field $u$ to have zero mean and unit standard deviation. The $\delta$-function is numerically approximated through a discretization

$$\delta(u - v) = \frac{1}{\Delta v} \left[ \mathcal{H} \left( u + \frac{\Delta v}{2} \right) - \mathcal{H} \left( u - \frac{\Delta v}{2} \right) \right] \qquad (5.4)$$

that is $\delta(u - v)$ is equal to $1/\Delta v$ if $u$ is between $v - \Delta v/2$ and $v + \Delta v/2$, and zero elsewhere.

This pixelization method introduces systematic residuals. It has been shown that residuals scale as the square of the bin-size ($\Delta v^2$) [62]. However, if the bin-size is too small the results can be affected by map noise. We pick the value $\Delta v = 0.5$. We find that this bin-size minimizes the residual obtained comparing the numerical equations (5.2) and (5.3) with the MFs analytical equations valid for a perfect Gaussian field.

Minkowski Tensors (MTs) are tensorial quantities that generalize the scalar MFs. MTs have been already used in cosmology to study CMB [63] and foreground [27] anisotropies. There are three rank-two MTs on the sphere, usually denoted as $W_k$. The three scalar MFs, $V_k$, are then given by the traces of $W_k$. We are particularly interested in $W_1$, also called the Contour Minkowski Tensor (CMT), which encodes shape and alignment information for structures. It can be numerically computed as

$$W_1(v) = \frac{1}{4} \frac{1}{N_{\text{pix}}} \sum_{\text{pixels}} \delta(u-v) \frac{1}{|\nabla u|} M \qquad (5.5)$$
where

$$M = \begin{pmatrix} u_{;\phi}^2 & -u_{;\phi}u_{;\theta} \\ -u_{;\phi}u_{;\theta} & u_{;\theta}^2 \end{pmatrix}. \qquad (5.6)$$

$W_1$ is proportional to the identity matrix if the structures have no elongation in any particular direction. $\lambda_+$, $\lambda_-$ are the two eigenvalues of $W_1$ such that $\lambda_+ > \lambda_-$. We define $\alpha$ as

$$\alpha = \frac{\lambda_+}{\lambda_-}. \qquad (5.7)$$

$\alpha = 1$ implies that the field preserves statistical isotropy (SI). In order to quantify the non-Gaussianity and anisotropy of data maps, we need to compare them with a set of suitable simulations. We then define, at each threshold $v$, the quantities

$$\chi_k = \frac{|\Delta V_k|}{\sigma_{V_k}}, \quad \Delta V_k(v) = V_k^{\text{data}}(v) - V_k^{\text{sim}}(v) \qquad (5.8)$$

where $V_k^{\text{data}}$ is the k-th MF computed from the data, $V_k^{\text{sim}}$ and $\sigma_{V_k}$ are the values obtained taking respectively the average and the standard deviation of the functionals determined from simulations. Analogously, we define at each threshold $v$ the quantities

$$\chi_{(W_1)_{ii}} = \frac{|\Delta(W_1)_{ii}|}{\sigma_{(W_1)_{ii}}}, \quad \Delta(W_1)_{ii} = (W_1)_{ii}^{\text{data}} - (W_1)_{ii}^{\text{sim}} \qquad (5.9)$$

where $(W_1)_{ii}$ stands for $(W_1)_{11}$ and $(W_1)_{22}$, which are the diagonal terms of the CMTs. The values $(W_1)_{ii}^{\text{sim}}$ and $\sigma_{(W_1)_{ii}}$ are computed from simulations. The same quantification could not be applied to $\alpha$ because its statistic follows the Beta probability distribution [27]. Thus we define the quantity

$$\chi_{\alpha} = \frac{\Delta\alpha}{\delta_1}\mathcal{H}(-\Delta\alpha) + \frac{\Delta\alpha}{\delta_2}\mathcal{H}(\Delta\alpha), \quad \Delta\alpha(v) = \alpha^{\text{data}}(v) - \alpha^{\text{sim}}(v) \qquad (5.10)$$

where $\delta_1$ and $\delta_2$ denote the 95% confidence interval and $\mathcal{H}$ is the Heaviside step function. Note that a value $|\chi_\alpha| > 1$ implies a deviation from the simulations outside the 95% confidence interval. In this and the following sections we will focus on the analysis and results obtained from the WMAP K-band polarization maps at resolution 1° and $N_{\text{side}} = 128$. The same analysis has been performed with 30 GHz Planck maps and the results presented in appendix C.

## 5.2 Masking and filtering

It is well known that the synchrotron emission is non-Gaussian and anisotropic on the full sky. In this analysis, we are mainly interested in understanding if this behaviour holds in regions of the sky where the emission is more diffuse, that is when the brightest areas (Galactic plane, the Spur and Fan regions) are masked.

We define two different masks for the analysis of the faintest 80% and 60% of the sky. To avoid possible leakage of power which can arise when computing power spectra in the presence of sharp boundaries between the masked and the unmasked regions, we apodize the masks with the `mask_apodization` (“C2”) routine of NaMaster<fnref n="9" /> with apodization scale of 5°. We also take into account a point source mask, as used in the filament finder analysis

<footnote n="9">⁹NaMaster is a public software package providing a general framework to estimate pseudo-$C_\ell$ angular power spectra [64] https://namaster.readthedocs.io.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Debiased WMAP K-band polarized intensity maps at 1° resolution, after band-pass filtering as defined in equation (3.4) with $\ell_{min} = 30$ and $\ell_{max} = 180$. The grey regions correspond to the missing pixels for the 80% (left) and 60% (right) masks.</figcaption></figure>

in section 3.3, apodized at 1°. To minimize the effects of the mask boundary in the sum in equations (5.1)–(5.5), we only include those pixels such that the smoothed mask value is larger than 0.9. For the two masks adopted here, this corresponds to sky fractions of 76.4% and 57.1%. The threshold has been conservatively chosen in order to maximize the statistical significance of our results. Note that we do not construct a specific filament mask, because from prior tests we have noticed that the Minkowski method is more reliable when we use a compact mask rather than a complex mask with many holes and islands.

Besides the region-dependency, in this analysis we want to test the statistical behavior on different scales. The maps are therefore filtered using the band-pass filter presented in section 3.2. We analyze different scale ranges varying $\ell_{min}$ in equation (3.4). We do not study the maps for multipoles smaller than $\ell_{min} = 30$ because the anisotropic nature of the emission at such large scale is well known, so we will only focus on $\ell_{min} \ge 30$. We fix $\ell_{max} = 180$, in accordance with the map resolution at 1°. In figure 7 the WMAP polarized intensity maps, filtered with the ($\ell_{min} = 30, \ell_{max} = 180$) band-pass, are shown. The region in grey corresponds to the pixels masked in the Minkowski analysis. We finally subtract the mean and then rescale with the standard deviation, where mean and standard deviation are computed from the data in the unmasked area of the maps.

## 5.3 Gaussian simulations

In order to quantify the non-Gaussianity and anisotropy of the polarized synchrotron emission, we need to compare the polarization data with a set of suitable simulations. To generate simulations, we first compute with NaMaster the polarization power spectra of the data maps in the unmasked areas. In fact, we compute the cross-spectra between the co-added WMAP 1 to 4 year sky maps and 5 to 9 year data. In this way, we reduce the effects of instrumental noise and systematics. Based on this spectrum, we then generate Gaussian and isotropic simulations of full-sky Q and U Stokes parameter maps at the data resolution using the HEALPix synfast routine. We add noise and compute the debiased polarized intensity of the simulations as described in section 2.3. We obtain 600 total polarization simulations which have a mean power spectrum compatible with the data, but are isotropic and generated from Gaussian-distributed Q and U. Note that the variance of the simulations includes both sample variance and noise. Finally, we compute the MFs and MTs for the simulations using the same band-pass filter and masks as used for the data.
## 5.4 Results

The results for the MFs and the CMT, derived from band-pass filtered data are shown in figure 8. We find that the non-Gaussian and anisotropic deviations are much lower when bright regions are masked. For the larger sky fraction ($f_{sky} = 80\%$), we find that the first MF ($V_0$) shows the largest deviation from Gaussianity. For most thresholds, the deviation is higher than $3\sigma$ and several values exceed $5\sigma$. The deviations of the other two MFs and the CMT diagonal terms exceed $3\sigma$ for several thresholds, few exceed $5\sigma$. $|\chi_\alpha|$ is lower than one except for one threshold, implying a weak deviation from isotropy. For the smaller sky fraction ($f_{sky} = 60\%$), for all the thresholds we find that the deviations of the three MFs and the CMT diagonal terms never exceed $3\sigma$, with average values (over all thresholds) lower than $1.2\sigma$. The value of $\alpha$ never exceeds the 95% interval.

In figure 9 we show results for different values of $\ell_{min}$. For the larger sky fraction ($f_{sky} = 80\%$), we find that the deviation decreases as $\ell_{min}$ increases, however, the MF deviation for some thresholds remains significantly high ($>3\sigma$), even at the smallest scales considered ($\ell_{min} = 130$). For the smaller sky fraction ($f_{sky} = 60\%$), the deviations remain almost constant with $\ell_{min}$, indicating consistency between the data and simulations for all thresholds and multipoles at the $3\sigma$ level.

Referring back to figure 7, we see that the areas where several filaments (IIIIn, IIIs, IV, VII, X, XVI, XVII) have been detected are masked by the 60% mask and not the 80% mask. This suggests that the contribution of the sky regions where these large complex structures are present affects the emission on scales smaller than $\sim 6^\circ$.

The non-Gaussian deviation decreases for smaller scales, corroborating the hypothesis that at small scales the emission tends to be more Gaussian. However, the effect of bright local structures is still not negligible on a scale of $\sim 1.5^\circ$. From this analysis we also learn that the non-Gaussian nature is mainly of the kurtosis type, since $V_0$ shows the largest deviations from simulations, as explained in [27]. These results at 23 GHz are in good agreement with those from the Planck 30 GHz polarization maps. More details can be found in appendix C.1. Moreover, very similar conclusions are obtained for intensity in the analysis of the 408 MHz maps by [27].

# 6 Non-Gaussian simulations

## 6.1 A model for non-Gaussian emission

Using MFs and MTs, we have determined clear statistical differences between the real synchrotron polarized emission and an isotropic and Gaussian-distributed model, even at small scales. In this section, we propose a simple way to generate polarized synchrotron maps which can better resemble the real statistical properties of the emission in the fainter regions of the sky on scales smaller than about $6^\circ$, which corresponds to multipoles $\ell > 30$.

From our statistical analysis, it is impossible to disentangle non-Gaussianity from anisotropy. That is, we are unable to determine whether the measured non-Gaussianity is an intrinsic feature of the synchrotron emission, or is due to an underlying anisotropic pattern of emission on small scales. However, even if we assume the latter case, because of the lack of information about the small scale distribution of the emission, we need a mechanism to simulate the non-Gaussianity.

We generate the simulations as follows. We start with three complex vectors, $\zeta^T$, $\zeta^E$, $\zeta^B$, generated from a Gaussian random distribution with zero mean and unit variance. The
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Upper panel: MFs (top), difference between data and Gaussian simulations (middle) and $\chi_k$ (bottom) as a function of threshold. Lower panel: CMT diagonal terms and $\alpha$ (top), difference between data and Gaussian simulations (middle) and $\chi_k$ (bottom) as a function of threshold. The maps are previously filtered with a band-pass ($\ell = 30-180$), error bars are computed from simulations.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9. The three MFs (top), the CMT diagonal terms and $\alpha$ (bottom) deviations as a function of the lower multipole cut, $l_{min}$, of the applied band-pass filter. The dots and the triangles represent respectively the average and 95% percentile values computed over all threshold values.</figcaption></figure>

$\zeta^E$ and $\zeta^B$ coefficients are then transformed according to the sinh-arcsinh transformation [65]

$$ \tilde{\zeta}^E = \sinh(\delta \arcsinh(\zeta^E) - \epsilon) \qquad (6.1) $$

$$ \tilde{\zeta}^B = \sinh(\delta \arcsinh(\zeta^B) - \epsilon) \qquad (6.2) $$

where $\zeta^E$ and $\zeta^B$ are the initial Gaussian-distributed vectors and $\tilde{\zeta}^E$ and $\tilde{\zeta}^B$ are the transformed ones. The use of this transformation is motivated by the fact that it allows one to control the level of non-Gaussianities with two parameters, $\epsilon$ and $\delta$, whose effects have a simple statistical interpretation. The parameter $\delta$ introduces symmetrically both positive and negative tails to the statistical distribution, which increases the excess kurtosis. The parameter $\epsilon$ controls the level of skewness. We find that the Gaussian case is correctly recovered with $(\epsilon, \delta) \to (0, 1)$, when performing simple consistency tests.

We assume, for the sake of simplicity, that the level of non-Gaussianity is the same for $E$ and $B$ in the part uncorrelated to $T$. This explains why we use the same parameters $(\epsilon, \delta)$ for both transformations in equations (6.1)–(6.2) and we do not transform $\zeta^T$. We are aware that in a more realistic context this could not be the case, but this goes beyond the scope of this analysis. The values used for $(\epsilon, \delta)$ are discussed in section 6.3. We then generate the spherical harmonic coefficients as follows

$$ a_j^T = \sqrt{C_l^{TT}} \zeta^T \qquad (6.3) $$

$$ \tilde{a}_j^E = \frac{C_l^{TE}}{\sqrt{C_l^{TT}}} \zeta^T + \sqrt{C_l^{EE} - \frac{(C_l^{TE})^2}{C_l^{TT}}} \tilde{\zeta}^E \qquad (6.4) $$

$$ \tilde{a}_j^B = \frac{C_l^{TB}}{\sqrt{C_l^{TT}}} \zeta^T + \sqrt{C_l^{BB} - \frac{(C_l^{TB})^2}{C_l^{TT}}} \tilde{\zeta}^B \qquad (6.5) $$
where the index $j$ refers to the $(\ell, m)$ pair. The $C_{\ell}$ values are computed from the WMAP yearly maps as described in section 5.3. Note that after rescaling, the $\tilde{a}_{lm}$-distributions are different from the $\tilde{\zeta}$-distributions, but still non-Gaussian. We generate the $(I, \tilde{Q}, \tilde{U})$ maps, which represent our small-scale template, from the $\tilde{a}_{lm}$ using the `alm2map` routine of HEALPix.

## 6.2 Model for anisotropic modulation of the emission

Besides the non-Gaussian nature, we have shown that at high multipoles the synchrotron emission is not isotropic. We simply assume that this is due only to a modulation caused by the large bright structures which are also visible at large scales. To simulate the effect, we divide the polarization WMAP map into 3 patches ($p_i$) according to the application of thresholds on $P$. We divide the unmasked 80% of the sky into two patches of 20%, $p_1$: 80–60%, $p_2$: 60–40%, and one of 40%, $p_3$: 40–0%, according to the sky brightness. For each patch, we smooth the boundaries with a 5° Gaussian and compute the polarization spectra $C_{\ell}^P(p_i)$. Then we define a spatially varying normalization factor

$$N_i = \sqrt{\left\langle \frac{C_{\ell}^{P}(p_i)}{C_{\ell}^{P}} \right\rangle} \qquad (6.6)$$

where the $C_{\ell}^P$ is computed from 80% of the sky and the average is computed over the multipoles $\ell \in [30, 180]$. $N$ is then smoothed with a 10° beam. As shown in figure 10 (top left), this method naturally introduces the effect at low-scales of most of the filamentary structures presented in section 4. We multiply the non-Gaussian $\tilde{Q}$ and $\tilde{U}$ maps, computed in the previous section, by the normalization factor

$$Q = N\tilde{Q}, \quad U = N\tilde{U}. \qquad (6.7)$$

The resulting $Q$ and $U$ maps form the final set of anisotropic and non-Gaussian simulations.

To our small-scale model, we add a large-scale template generated directly from the WMAP data. In order to match the maps, we smooth the small-scale and the large-scale templates with respectively the functions $W(\ell)$ and $(1 - W(\ell))$, where

$$W(\ell) = \frac{1}{2} \left( 1 - \tanh \left( \frac{\ell - \ell_0}{\Delta\ell_0} \right) \right), \qquad (6.8)$$

with $\ell_0 = 20$ and $\Delta\ell = 5$, chosen such that at the scales of interest in this work, that is smaller than 6°, the simulations are mainly driven by our small-scale model. Finally, in order to compare simulations with data, we add noise and compute the debiased polarized intensity with the MAS estimator.

## 6.3 Tuning and results

We have described above the method we use to generate polarized synchrotron simulations on scales below 6°. The level of non-Gaussianity of the simulations is controlled by the two parameters $(\epsilon, \delta)$ in equation (6.1)–(6.2), which are related respectively to the skewness $S$ and kurtosis $K$ of the maps.

Similarly to equation (5.8), we define the quantities which measure the deviation between data and simulations as

$$\chi_{\mu} = \sum_{k=2}^{4} \frac{\Delta\mu_{k}}{\sigma_{\mu_{k}}}, \quad \Delta\mu_{k} = |\mu_{k}^{\text{data}} - \mu_{k}^{\text{sim}}| \qquad (6.9)$$
where $\mu_2 = \sigma_P = \sqrt{\langle\langle(p - \langle p \rangle)^2\rangle}$ is the standard deviation, $\mu_3 = S = \langle\langle(p - \langle p \rangle)^3\rangle/\sigma_P^3$ is the skewness and $\mu_4 = K = \langle\langle(p - \langle p \rangle)^4\rangle/\sigma_P^4$ is the kurtosis. The values $\mu_k^{\text{sim}}$ and $\sigma_{\mu_k}$ are computed taking respectively the average and standard deviation over the simulations. Moreover, we define a variable which quantifies the spectral deviations as

$$ \chi_{\text{sp}} = \sum_{\ell=30}^{180} \frac{\Delta C_{\ell}}{\sigma_{\ell}}, \quad \Delta C_{\ell} = |C_{\ell}^{\text{data}} - C_{\ell}^{\text{sim}}| \qquad (6.10) $$

where $C_{\ell}^{\text{sim}}$ and $\sigma_{\ell}$ are computed taking respectively the average and standard deviation over the simulations. Note that the $C_{\ell}$ are binned with a range $\Delta\ell = 10$. We vary $\epsilon$ and $\delta$ over the ranges $(-1, 1)$ and $(0.3, 1.5)$, respectively. For each pair of values, we generate 100 total polarization simulations, smoothing them with the band-pass filter ($\ell_{\text{min}} = 30, \ell_{\text{min}} = 180$). We find that the quantity $\chi_{\mu} + \chi_{\text{sp}}$ is minimised for $\epsilon = -0.46$ and $\delta = 0.78$, as shown in figure 10 (top right). These define the values to generate a set of reference simulations. An example of a simulated sky, and the comparison between the simulated and real power spectra are shown in figure 10 (bottom).

The deviation in the MF and CMT measures between data and simulations, when the band-pass ($\ell_{\text{min}} = 30, \ell_{\text{min}} = 180$) filter is applied, are shown in figure 11. In figure 12, we show the deviations for different $\ell_{\text{min}}$. The simulations agree with the data at the $3\sigma$ level for both the 80% and 60% masks. This suggests that the non-Gaussianity introduced with the sinh-arcsinh transformation (equations (6.1)–(6.2)) is more relevant in the bright regions, instead is mitigated in the more faint sky emission. As noted in section 6.1, we do not know if the non-Gaussianity we introduce is related to the intrinsic nature of the synchrotron emission, or due to an underlying anisotropic emission at small scales. However, the mechanism to produce this non-Gaussianity provides a simple method to tune our simulations by means of two simple parameters $(\epsilon, \delta)$. Overall, the non-Gaussian and anisotropic simulations seem to reproduce the statistical properties of the polarized synchrotron emission well.

In appendix C.2, we test how our method performs simulating the synchrotron emission at 30 GHz, as observed by Planck. We keep the parameters $\epsilon, \delta$ and $N$ as fitted from WMAP. We observe that the method performs well at small scales ($>2.5^\circ$), but under-performs for larger scales.

The goal of this section is to show a simple data-driven way to generate more realistic simulations. There are several assumptions and limitations which are negligible at the current data sensitivity, but will not be the case for future experiments. For example, we assume that the E and B modes are equally non-Gaussian, even if we do not have evidence of it. We do not take into account the correlation of the non-Gaussianity between intensity and polarization. We consider $\epsilon$ and $\delta$ as fixed values, but in a more realistic context, they could vary with scale and space. Moreover, a frequency dependence in our model ($\epsilon, \delta, N$), for example due to a possible decorrelation between frequencies as observed in the dust emission [53], could not be excluded. Such effects could explain the poor performance at 30 GHz for the large scales.

# 7 Summary and conclusions

In the analysis we have covered two important aspects of the polarized synchrotron emission: the presence of large filamentary structures outside the Galactic plane, and its statistical properties at small scales. The analysis has been performed on the debiased polarized amplitude maps at 23 and 30 GHz as observed respectively by WMAP and Planck. We
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Top Left: Spatially varying normalization factor used to modulate the simulations in real space. Top Right: $\chi_{sp} + \chi_{\mu}$ as function of the parameters $\epsilon$ and $\delta$, on logarithmic scale. Bottom Left: Simulation of the *WMAP* polarized intensity map. Bottom Right: Total polarization power spectrum of the *WMAP* K-band (black dots) compared to the $1\sigma$ (dark green) and $2\sigma$ (light green) intervals obtained from the variance of the simulations.</figcaption></figure>

developed a filament finder routine based on a friend-of-friend recursive algorithm, which detects elongated coherent emission in the sky. The method has been tested with foreground simulations including a toy model of filamentary structure. We identify 19 filaments which are detected, at least in part, in both *Planck* and *WMAP*. Some of the filaments have been already reported in the literature as observed in radio sky or in previous *WMAP* analysis. Five of them are reported for the first time in this work. We analysed some properties of the detected filaments.

* We compute the polarization fraction as the ratio of the *WMAP* and *Planck* data maps with respect to a **Commander** intensity template rescaled with a spectral index $\beta_I = -3$. Typically, we find that the polarization fractions of the filaments are larger than for the areas outside the filaments, excluding the Galactic plane. For two filaments, both located in the NPS, we find values of about 30%.

* We study the polarization spectral indices of the filaments from the *Q*, *U* maps. We find consistent spectral indices of about -3.1 inside and outside the filaments. However, the $\beta$ values span a very broad range, from -3.59 to -2.17.

* Applying the filament finder algorithm to the $P_E$ and $P_B$ maps, we find that most of the filaments, especially in the NPS, are clearly visible in *E*, but not in *B*.

* Applying the finder method to the Haslam map, we observe that for some bright filaments in polarization we do not detect bright counterparts in intensity.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Upper panel: MFs (top), the difference between data and *non-Gaussian* simulations (middle) and $\chi_k$ (bottom) as a function of threshold. Lower panel: CMT diagonal terms and $\alpha$ (top), difference between data and *non-Gaussian* simulations (middle) and $\chi_k$ (bottom) as a function of threshold. The maps are previously filtered with a band-pass ($\ell = 30-180$), error bars are computed from the standard deviation of the simulations. For the bottom panels, we use the same axis ranges as in figure 8 to allow a direct comparison.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** The three MFs (top), the CMT diagonal terms and $\alpha$ (bottom) deviations from the non-Gaussian simulations as function of the lower multipole cut, $\ell_{min}$, of the applied band-pass filter. The dots and the triangles represent respectively the average and 95% percentile computed over all threshold values. For comparison purposes, we use the same axes as in figure 9.</figcaption></figure>

Using Minkowski functionals and tensors, we have analysed the non-Gaussianity and statistical isotropy of the polarised *WMAP* and *Planck* maps. We focused on the faintest 80% and 60% of the sky. We compared the results obtained from data with a set of Gaussian and isotropic simulations. We summarize our findings.

• Analysing the 80% sky fraction, we see large deviations ($> 3\sigma$) from Gaussianity and isotropy at 6° scale. The deviations decrease towards smaller scales, even if they remain significantly high down to 1.5°.

• Analysing the 60% sky fraction, we obtain consistent results between data and simulations for all the considered thresholds and multipole ranges at the $3\sigma$ level.

• These results suggest that the large filaments are the main source of non-Gaussianity, even at small scales. When those filamentary structures are masked, the Gaussian and isotropic simulations resemble well the diffuse emission at the *WMAP* and *Planck* resolution.

Finally, we present a simple data-driven method used to generate non-Gaussian and anisotropic simulations. We generate non-Gaussian harmonic coefficients by mean of a simple transformation. We account for the anisotropies with a normalization template that resemble the diffuse filamentary structures. The simulations are then fitted in order to match the spectral and statistical properties of the 80% sky coverage of the data maps.

## Acknowledgments

The authors would like to thank the Spanish Agencia Estatal de Investigación (AEI, MICIU) for the financial support provided under the projects with references PID2019-110610RB-C21, ESP2017-83921-C2-1-R and AYA2017-90675-REDC, co-funded with EU FEDER funds,
and acknowledge support from Universidad de Cantabria and Consejería de Universidades, Igualdad, Cultura y Deporte del Gobierno de Cantabria via the “Instrumentación y ciencia de datos para sondear la naturaleza del universo” project as well as from the Unidad de Excelencia María de Maeztu (MDM-2017-0765). FAM is supported by a fellowship funded by the Unidad de Excelencia María de Maeztu. The authors also thank M. López-Caniego for providing point source masks for CMB data. We acknowledge Santander Supercomputacion support group at the University of Cantabria who provided access to the supercomputer Altamira Supercomputer at the Institute of Physics of Cantabria (IFCA-CSIC), member of the Spanish Supercomputing Network, for performing simulations/analyses. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231. Some of the presented results are based on observations obtained with Planck,<fnref n="10" /> an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. We also acknowledge the Legacy Archive for Microwave Background Data Analysis (LAMBDA). Support for LAMBDA is provided by the NASA office of Space Science. Some of the results in this work have been derived using the HEALPix [66, 67], NaMaster [64], PySM [25], Matplotlib [68], NumPy [69] and SciPy [70] Python packages.

## A Power spectra

In this section, we give a very brief review of the statistical quantities we use in this work, motivated by standard cosmological practises. CMB experiments usually produce data in the form of three pixelized maps, $T$ for intensity and $Q$ and $U$ Stokes parameters for polarization. On the sky, these fields are usually expanded in terms of spherical harmonics

$$T = \sum_{\ell m} a_{\ell m} Y_{\ell m} \qquad (A.1)$$

$$(Q \pm iU) = \sum_{\ell m} a_{\pm 2, \ell m} \pm 2 Y_{\ell m} \qquad (A.2)$$

where $Y_{\ell m}$ and $\pm 2 Y_{\ell m}$ are respectively the standard and tensor (spin-2) spherical harmonics on a 2-sphere. The quantities $a_{\ell m}$ are the so-called spherical harmonic coefficients. Details of the mathematical formalism can be found in [71, 72].

If we define the linear combinations

$$a_{E,\ell m} = -\frac{1}{2}(a_{2,\ell m} + a_{-2,\ell m}) \quad a_{B,\ell m} = -\frac{1}{2i}(a_{2,\ell m} - a_{-2,\ell m}) \qquad (A.3)$$

we can decompose the polarization emission into two scalar fields, the gradient-like $E$ mode and the curl-like $B$ mode

$$E = \sum_{\ell m} a_{E, \ell m} Y_{\ell m} \quad B = \sum_{\ell m} a_{B, \ell m} Y_{\ell m}. \qquad (A.4)$$

The harmonic coefficients may be combined into the angular power spectrum

$$C_{\ell}^{XY} = \frac{1}{2\ell + 1} \sum_{m} \langle a_{X,\ell m}^{*} a_{Y,\ell m} \rangle, \quad X, Y = T, E, B \qquad (A.5)$$

<footnote n="10">¹⁰http://www.esa.int/Planck.</footnote>
which represent: the auto-correlations of temperature and polarization modes denoted by $TT$, $EE$, and $BB$, the cross-correlation between temperature and polarization denoted by $TE$ and $TB$, and the cross-correlation between polarization modes denoted by $EB$. For a Gaussian and isotropic field, all the statistical properties are captured by these two-point statistics.

Experimental observations are affected by the instrumental (beam) response and the pixelization process. The observed maps can then be written as the convolution of the actual sky signal with the instrumental beam ($B$) and the pixel window function ($W$). The latter is a function of the resolution at which the maps are produced. In harmonic space, it implies that

$$a_{lm}^{\text{obs}} = a_{lm} B_{l} W_{l} \qquad (A.6)$$

$$C_{l}^{\text{obs}} = C_{l} B_{l}^{2} W_{l}^{2} \qquad (A.7)$$

where $B_l$ and $W_l$ are respectively the harmonic transformations of the instrumental beam and the pixel window function. When analyzing maps produced by different experiments at resolution different from the one in which the original maps are produced, as discussed in this work, it is appropriate to smooth the maps to a common resolution. This can be achieved in harmonic space by

$$a_{lm}^{\text{out}} = a_{lm}^{\text{in}} \frac{B_{l}^{\text{out}} W_{l}^{\text{out}}}{B_{l}^{\text{in}} W_{l}^{\text{in}}} \qquad (A.8)$$

where $B_l^{\text{in}}$ and $B_l^{\text{out}}$ are respectively the instrumental and the required beams, $W_l^{\text{out}}$ and $W_l^{\text{in}}$ are the pixel window functions at the final and initial resolutions.

# B Finder algorithm

## B.1 Toy model for filaments

In this section, we test the performance of the filament finder algorithm using foreground simulations including a toy model of filamentary structure. Each simulation is computed as the sum of different independent components

$$S = S_{\text{Gal}} + S_{\text{dif}} + S_{\text{Loops}} + S_{\text{noise}}, \qquad (B.1)$$

where:

* $S_{\text{Gal}}$ is a Galactic plane simulation. The template is generated from the WMAP K-band P map, smoothed to a resolution of 5° and filtered with a low-pass filter $f(l) = [1 - \tanh((l - 10)/10)]/2$. In this way we preserve the Galactic morphology on large scales ($l < 10$) whilst removing the small scales corresponding to the real filaments.

* $S_{\text{dif}}$ is a diffuse Gaussian template created with the synfast routine using the power spectra model: $C_l \propto (l/80)^{-2.9}$ [10]. The simulated map is filtered with a high-pass filter $f(l) = [1 - \tanh((l + 10)/10)]/2$, which only allows multipoles $l > 10$.

* $S_{\text{Loops}}$ is a template where different loops are projected onto the sphere. The loops are based on filaments observed in WMAP [20]. Each loop has been generated with a width in the range 2–4°. In order to simulate both thin and diffuse filaments, we smooth the loops with a 1.5° or a 3° Gaussian beam. The loops are shown in figure 14 (top right).

* $S_{\text{noise}}$ is a noise simulation with properties estimated from the WMAP noise covariance matrices.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Band-pass filters defined in equation (3.4) used for the detection of the bright (blue) and the weak (orange) filaments.</figcaption></figure>

The simulation *S* is produced at *N*<sub>side</sub> = 128 and a resolution of 1°, as used for the data. The Galactic and diffuse templates have been re-scaled in order to match the data signal-to-noise ratio. We tested different amplitudes, locations and radii for the loops, although in the following we will only refer to the case including Loops I, III, GCS, VII and XI, at a signal-to-noise ratio of 5, as shown in figure 14. We filter the maps with the filters in figure 13 and apply the friends-of-friends algorithm to 100 simulations. Note that each simulation has the same *S*<sub>Gal</sub> and *S*<sub>Loop</sub>, but different realization of *S*<sub>dif</sub> and *S*<sub>noise</sub>. Figure 14 presents an example of a simulation (top left), and the corresponding detected structures before (bottom left) and after (bottom right) the minimal length criteria is applied, as described in section 3.1.

For each simulation, we recover on average 71.0% (±1.4%) of the original filaments. However, we also assign a detection of filamentary structure to around 7.2% (±0.7%) of the sky which is not associated with any input loops. From figure 14, we observe that the filament finder mostly fails to detect parts of filaments close or tangential to the Galactic plane, where the strong Galactic emission dominates. We also point out that the detection can fail in those areas where two or more loops overlap, because the orientation angle in those pixels is the result of the average over different loops. The detections that are not associated with any input loop mainly arise in the regions with the lowest signal-to-noise ratio, suggesting that the noise is the cause. However, we note that it is possible to identify most of these spurious detections by comparing two simulations with different noise realisations. In practice, in our main analysis with real data, this is achieved by comparing the results of two independent maps, from *WMAP* and *Planck*, which allows us to reduce the number of spurious detections.

Note that the quantitative results presented in this appendix are obtained with reference to the *WMAP* data. However, all qualitative considerations also apply to the *Planck* data. Possible differences in the performance of the algorithm are mainly attributable to the fact that the difference of the filament brightness to the diffuse background is greater in the *WMAP* map than in the *Planck* map. In addition, the different distributions of noise for the two experiments could also have an impact on the performance of the algorithm, although we expect it to be subdominant.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Top Left: Toy filamentary foreground simulation. Top Right: Loop template used in the simulation. Bottom Left: All the detections found with the friends-of-friends recursive algorithm. Bottom Right: final result of the filament finder algorithm.</figcaption></figure>

## B.2 Minimal length criteria

The filament finder method presented in section 3.1 is a simple friends-of-friends recursive algorithm based on the properties of single pixels. When a group of coherent and bright pixels is identified, it is not obvious if it is part of a filamentary structure or not. Considering the positive nature of the polarization intensity, regions where the noise is strong can confuse the detection. Moreover, the synchrotron diffuse background can also have a detrimental effect. In order to reduce spurious detections, we reject structures where the maximum pixel-pair angular distance is smaller than a threshold value $L_{\text{th}}$.

We find the best value for $L_{\text{th}}$ by analysing simulations that do not contain loop structures. As in section B.1, we generate 100 diffuse synchrotron simulations from $S = S_{\text{Gal}}+S_{\text{dif}}+S_{\text{noise}}$, where the single components are described in the previous section. Note that in this case the simulations do not include the $S_{\text{Loops}}$ term.

We apply the finder algorithm to each simulation which now can only detect spurious signals due to noise and the diffuse emission. Figure 15 shows an example of a detection (left) and the distribution of the lengths determined from the simulations (right). We find that 68% of detections have a length smaller than roughly 3.1°, 95% smaller than 10.2° and 99% smaller than 17.5°. From this result, we pick the threshold value $L_{\text{th}} = 10°$. Note that this estimate holds for the pessimistic scenario of a loop-less foreground. In a more realistic case, i.e. including bright filaments, the algorithm would rely on a larger $P_{\text{th}}$, so a part of the noise detection would not exceed the threshold, and we would get less spurious detections.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Left: All the detections found with the friends-of-friends recursive algorithm from a loop-less simulation. Right: Distribution of the maximum angular lengths obtained from 100 loop-less simulations.</figcaption></figure>

## C Planck statistical properties

### C.1 Gaussian simulations

In the main text, we analyse the statistical properties of the *WMAP* K-band maps. The choice is motivated by the fact that at 23 GHz the synchrotron emission is much stronger than at the 30 GHz *Planck* frequency channel. However, as a consistency check, in this section we present the results obtained when analysing the *Planck* data. We use the same masks and filters presented in section 5.2.

In order to quantify the non-Gaussianity and anisotropy of the synchrotron emission observed by *Planck*, we need to compare data with a set of suitable simulations. We compute pseudo-spectra in the unmasked regions cross-correlating A/B split maps.<fnref n="11" /> From the spectra, we generate 600 Gaussian and isotropic simulations of *Q* and *U*, add noise, then compute the unbiased polarized intensity.

Results are shown in figure 16. For the larger sky fraction ($f_{sky} = 80\%$), we find that for the cases with $\ell_{min} < 80$, even averaging over the thresholds, the deviation exceeds $3\sigma$. The deviation decreases when $\ell_{min}$ increases, however, for all the quantities (except ($W_1$)$_{22}$), some thresholds remain significantly higher than $3\sigma$. For the smaller sky fraction ($f_{sky} = 60\%$), we generally find consistency between the data and simulations. These results are in substantial agreement with those determined with *WMAP* at 23 GHz, corroborating the discussion in section 5.4.

### C.2 Non-Gaussian simulations

In section 6, we present a data-driven method to simulate the polarized synchrotron emission at 23 GHz. In this section, we show how the model performs in reproducing the emission at 30 GHz. We use the same spatially varying normalization factor (see figure 10) and $(\epsilon, \delta)$ parameters to introduce anisotropies and non-Gaussianities as for *WMAP*. The results are shown in figure 17.

The simulations agree with the data at the $3\sigma$ level for the 80% mask for those cases with $\ell_{min} > 70$, and for all the considered multipole ranges for the 60% mask. The model seems to under-perform for the cases with $\ell_{min} \le 70$ when considering the 80% mask, although the

<footnote n="11">¹¹The PR4 provides A/B splits for data maps and simulations [32]. For the 30 GHz frequency channel, the A and B subsets are obtained respectively combining maps from years 1 and 3, and years 2 and 4.</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Top: The three MFs, and bottom: the CMT diagonal terms and the $\alpha$ deviations from the Gaussian simulations, computed with *Planck* data, as a function of the low multipole cut $l_{min}$ from the applied band-pass filter. The dots and the triangles represent respectively the average and 95% percentile values computed over all threshold values.</figcaption></figure>

deviations from the data are clearly not so pronounced as when using Gaussian simulations. Considering that the largest deviation comes from $(W_1)_{11}$, it is reasonable to think that we are not correctly taking into account the anisotropy of the field.

It has been shown, even in this work, that the polarization spectral index shows spatial variations, and bright structures at 23 GHz are less detectable at 30 GHz. This suggests that the spatially varying normalization factor computed only from the *WMAP* data, could also depend on frequency. In addition, given that the $V_1$ and $V_2$ values computed with simulations deviate from the data when considering multipoles $\ell < 60$, we can not exclude the possibility that the non-Gaussianity level could also depend on frequency, which in our model translates into $\epsilon = \epsilon(\nu)$ and $\delta = \delta(\nu)$.
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Top: The three MFs, and bottom: the CMT diagonal terms and $\alpha$ deviations from the non-Gaussian simulations, computed with Planck data, as a function of the low multipole cut $l_{min}$ from the applied band-pass filter. The dots and the triangles represent respectively the average and 95% percentile computed over all threshold values. For comparison purposes, we use the same ranges as in figure 16.</figcaption></figure>

References

[1] R. Génova-Santos et al., QUIJOTE scientific results — II. Polarisation measurements of the microwave emission in the Galactic molecular complexes W43 and W47 and supernova remnant W44, Mon. Not. Roy. Astron. Soc. 464 (2017) 4107 [arXiv:1605.04741] [INSPIRE].

[2] PAMELA collaboration, The cosmic-ray electron flux measured by the PAMELA experiment between 1 and 625 GeV, Phys. Rev. Lett. 106 (2011) 201101 [arXiv:1103.2880] [INSPIRE].

[3] FERMI-LAT collaboration, Fermi-LAT observations of the diffuse $\gamma$-ray emission: implications for cosmic rays and the interstellar medium, Astrophys. J. 750 (2012) 3 [arXiv:1202.4039] [INSPIRE].

[4] FERMI-LAT collaboration, Fermi LAT observations of cosmic-ray electrons from 7 GeV to 1 TeV, Phys. Rev. D 82 (2010) 092004 [arXiv:1008.3999] [INSPIRE].

[5] U. Fuskeland et al., Constraints on the spectral index of polarized synchrotron emission from WMAP and Faraday-corrected S-PASS data, Astron. Astrophys. 646 (2021) A69 [arXiv:1909.05923] [INSPIRE].

[6] R.D. Davies et al., A determination of the spectra of Galactic components observed by the Wilkinson Microwave Anisotropy Probe, Mon. Not. Roy. Astron. Soc. 370 (2006) 1125 [astro-ph/0511384] [INSPIRE].

[7] BICEP2 and PLANCK collaborations, Joint analysis of BICEP2/Keck array and Planck data, Phys. Rev. Lett. 114 (2015) 101301 [arXiv:1502.00612] [INSPIRE].

[8] PLANCK collaboration, Planck 2018 results. IV. Diffuse component separation, Astron. Astrophys. 641 (2020) A4 [arXiv:1807.06208] [INSPIRE].

[9] N. Krachmalnicoff et al., S-PASS view of polarized Galactic synchrotron at 2.3 GHz as a contaminant to CMB observations, Astron. Astrophys. 618 (2018) A166 [arXiv:1802.01145] [INSPIRE].
[10] F.A. Martire, R.B. Barreiro and E. Martínez-González, Characterization of the polarized synchrotron emission from Planck and WMAP data, JCAP 04 (2022) 003 [arXiv:2110.12803 [INSPIRE]].

[11] T.L. Svalheim et al., BeyondPlanck XIV. Polarized foreground emission between 30 and 70GHz, arXiv:2011.08503 [INSPIRE].

[12] R. de Belsunce, S. Gratton and G. Efstathiou, Testing for spectral index variations in polarized CMB foregrounds, Mon. Not. Roy. Astron. Soc. 517 (2022) 2855 [arXiv:2205.13968] [INSPIRE].

[13] E. de la Hoz et al., QUIJOTE scientific results — VIII. Diffuse polarized foregrounds from component separation with QUIJOTE-MFI, Mon. Not. Roy. Astron. Soc. 519 (2023) 3504 [arXiv:2301.05117] [INSPIRE].

[14] C.G.T. Haslam et al., A 408 MHz all-sky continuum survey. I-Observations at southern declinations and for the North Polar region, Astron. Astrophys. 100 (1981) 209.

[15] C.G.T. Haslam, C.J. Salter, H. Stoffel and W.E. Wilson, A 408 MHz all-sky continuum survey. II. The atlas of contour maps, Astron. Astrophys. Suppl. Ser. 47 (1982) 1 [INSPIRE].

[16] M. Remazeilles et al., An improved source-subtracted and destriped 408 MHz all-sky map, Mon. Not. Roy. Astron. Soc. 451 (2015) 4311 [arXiv:1411.3628] [INSPIRE].

[17] M.I. Large, M.J.S. Quigley and C.G.T. Haslam, A new feature of the radio sky, Mon. Not. Roy. Astron. Soc. 124 (1962) 405.

[18] M.J.S. Quigley and C.G.T. Haslam, Structure of the Radio Continuum Background at High Galactic Latitudes, Nature 208 (1965) 741.

[19] M.I. Large, M.J.S. Quigley and C.G.T. Haslam, A radio study of the north polar spur. II, A survey at low declinations, Mon. Not. Roy. Astron. Soc. 131 (1966) 335.

[20] M. Vidal, C. Dickinson, R.D. Davies and J.P. Leahy, Polarized radio filaments outside the galactic plane, Mon. Not. Roy. Astron. Soc. 452 (2015) 656 [arXiv:1410.4438] [INSPIRE].

[21] WMAP collaboration, Nine-year wilkinson microwave anisotropy probe (WMAP) observations: Final maps and results, Astrophys. J. Suppl. 208 (2013) 20 [arXiv:1212.5225] [INSPIRE].

[22] T.A.T. Spoelstra, Galactic Loops as Supernova Remnants in the Local Galactic Magnetic Field, Astron. Astrophys. 24 (1973) 149.

[23] M. Tegmark, D.J. Eisenstein, W. Hu and A. de Oliveira-Costa, Foregrounds and forecasts for the cosmic microwave background, Astrophys. J. 530 (2000) 133 [astro-ph/9905257] [INSPIRE].

[24] V. Jelić et al., Foreground simulations for the LOFAR-epoch of reionization experiment, Mon. Not. Roy. Astron. Soc. 389 (2008) 1319 [arXiv:0804.1130] [INSPIRE].

[25] B. Thorne, J. Dunkley, D. Alonso and S. Naess, The Python Sky Model: software for simulating the Galactic microwave sky, Mon. Not. Roy. Astron. Soc. 469 (2017) 2821 [arXiv:1608.02841] [INSPIRE].

[26] A. Waelkens et al., Simulating polarized galactic synchrotron emission at all frequencies — The Hammurabi code, Astron. Astrophys. 495 (2009) 697 [arXiv:0807.2262] [INSPIRE].

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

[28] F. Rahman, P. Chingangbam and T. Ghosh, Statistical properties of Galactic synchrotron temperature and polarization maps — a multi-frequency comparison, arXiv:2212.06076 [INSPIRE].

[29] PLANCK collaboration, Planck 2015 results. X. Diffuse component separation: Foreground maps, Astron. Astrophys. 594 (2016) A10 [arXiv:1502.01588] [INSPIRE].

[30] SIMONS OBSERVATORY collaboration, The Simons Observatory: Science goals and forecasts, JCAP 02 (2019) 056 [arXiv:1808.07445] [INSPIRE].
[31] T. Matsumura et al., *Mission design of LiteBIRD*, J. Low Temp. Phys. **176** (2014) 733 [arXiv:1311.2847] [INSPIRE].

[32] PLANCK collaboration, *Planck intermediate results. LVII. Joint Planck LFI and HFI data processing*, Astron. Astrophys. **643** (2020) A42 [arXiv:2007.04997] [INSPIRE].

[33] PLANCK collaboration, *Planck 2015 results. XXV. Diffuse low-frequency Galactic foregrounds*, Astron. Astrophys. **594** (2016) A25 [arXiv:1506.06660] [INSPIRE].

[34] S. Plaszczynski, L. Montier, F. Levrier and M. Tristram, *A novel estimator of the polarization amplitude from normally distributed Stokes parameters*, Mon. Not. Roy. Astron. Soc. **439** (2014) 4048 [arXiv:1312.0437] [INSPIRE].

[35] J.F.C. Wardle and P.P. Kronberg, *The linear polarization of quasi-stellar radio sources at 3.71 and 11.1 centimeters*, Astrophys. J. **194** (1974) 249.

[36] M. Vidal, J.P. Leahy and C. Dickinson, *A new polarization amplitude bias reduction method*, Mon. Not. Roy. Astron. Soc. **461** (2016) 698 [arXiv:1410.4436] [INSPIRE].

[37] N.A. Bond, M.A. Strauss and R. Cen, *Crawling the Cosmic Network: Identifying and Quantifying Filamentary Structure*, Mon. Not. Roy. Astron. Soc. **409** (2010) 156 [arXiv:1003.3237] [INSPIRE].

[38] N. Bond, M. Strauss and R. Cen, *Crawling the Cosmic Network: Exploring the Morphology of Structure in the Galaxy Distribution*, Mon. Not. Roy. Astron. Soc. **406** (2010) 1609 [arXiv:0903.3601] [INSPIRE].

[39] PLANCK collaboration, *Planck intermediate results. XXXVIII. E- and B-modes of dust polarization from the magnetized filamentary structure of the interstellar medium*, Astron. Astrophys. **586** (2016) A141 [arXiv:1505.02779] [INSPIRE].

[40] C. Monteserín, R.B. Barreiro, J.L. Sanz and E. Martínez-González, *Scalar statistics on the sphere: application to the cosmic microwave background*, Mon. Not. Roy. Astron. Soc. **360** (2005) 9 [astro-ph/0511308] [INSPIRE].

[41] F.R. Hampel, *The influence curve and its role in robust estimation*, J. Am. Statist. Assoc. **69** (1974) 383.

[42] R.W. Komm et al., *Multitaper spectral analysis and wavelet denoising applied to helioseismic data*, Astrophys. J. **519** (1999) 407.

[43] F. Argüeso et al., *Detection/estimation of the modulus of a vector. Application to point-source detection in polarization data*, Mon. Not. Roy. Astron. Soc. **395** (2009) 649 [arXiv:0906.0893] [INSPIRE].

[44] M. López-Caniego et al., *Comparison of filters for the detection of point sources in Planck simulations*, Mon. Not. Roy. Astron. Soc. **370** (2006) 2047 [astro-ph/0606199] [INSPIRE].

[45] R. Hanbury Brown, R.D. Davies and C. Hazard, *A curious feature of the radio sky*, Observatory **80** (1960) 191.

[46] C. Dickinson, *Large-Scale Features of the Radio Sky and a Model for Loop I*, Galaxies **6** (2018) 56 [INSPIRE].

[47] G. Dobler et al., *The Fermi haze: a gamma-ray counterpart to the microwave haze*, Astrophys. J. **717** (2010) 825 [arXiv:0910.4583] [INSPIRE].

[48] M. Su, T.R. Slatyer and D.P. Finkbeiner, *Giant gamma-ray bubbles from Fermi-LAT: Agn activity or bipolar galactic wind?*, Astrophys. J. **724** (2010) 1044 [arXiv:1005.5480] [INSPIRE].

[49] E. Orlando and A. Strong, *Galactic synchrotron emission with cosmic ray propagation models*, Mon. Not. Roy. Astron. Soc. **436** (2013) 2127 [arXiv:1309.2947] [INSPIRE].

[50] PLANCK collaboration, *Planck 2015 results. XXVI. The Second Planck Catalogue of Compact Sources*, Astron. Astrophys. **594** (2016) A26 [arXiv:1507.02058] [INSPIRE].
[51] U. Fuskeland, I.K. Wehus, H.K. Eriksen and S.K. Næss, *Spatial variations in the spectral index of polarized synchrotron emission in the 9-yr WMAP sky maps*, *Astrophys. J.* **790** (2014) 104 [arXiv:1404.5323 [INSPIRE]].

[52] I.K. Wehus, U. Fuskeland and H.K. Eriksen, *The effect of systematics on polarized spectral indices*, *Astrophys. J.* **763** (2013) 138 [arXiv:1201.6348 [INSPIRE]].

[53] PLANCK collaboration, *Planck 2018 results. XI. Polarized dust foregrounds*, *Astron. Astrophys.* **641** (2020) A11 [arXiv:1801.04945 [INSPIRE]].

[54] F. Guidi et al., *QUIJOTE scientific results — VI. The Haze as seen by QUIJOTE*, *Mon. Not. Roy. Astron. Soc.* **519** (2023) 3460 [arXiv:2301.05115 [INSPIRE]].

[55] H. Liu, J. Creswell and P. Naselsky, *E and B families of the Stokes parameters in the polarized synchrotron and thermal dust foregrounds*, *JCAP* **05** (2018) 059 [arXiv:1804.10382 [INSPIRE]].

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

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

[58] A. Carones et al., *Minkowski Functionals of CMB polarisation intensity with Pynkowski: theory and application to Planck data*, arXiv:2211.07562 [INSPIRE].

[59] P. Chingangbam and C. Park, *Residual foreground contamination in the WMAP data and bias in non-Gaussianity estimation*, *JCAP* **02** (2013) 031 [arXiv:1210.2250 [INSPIRE]].

[60] S. Rana, T. Ghosh, J.S. Bagla and P. Chingangbam, *Non-Gaussianity of diffuse Galactic synchrotron emission at 408 MHz*, *Mon. Not. Roy. Astron. Soc.* **481** (2018) 970 [arXiv:1806.01565 [INSPIRE]].

[61] P. Chingangbam et al., *Tensor Minkowski Functionals for random fields on the sphere*, *JCAP* **12** (2017) 023 [arXiv:1707.04386 [INSPIRE]].

[62] E.A. Lim and D. Simon, *Can we detect Hot or Cold spots in the CMB with Minkowski Functionals?*, *JCAP* **01** (2012) 048 [arXiv:1103.4300 [INSPIRE]].

[63] P.K. Joby et al., *Search for anomalous alignments of structures in Planck data using Minkowski Tensors*, *JCAP* **01** (2019) 009 [arXiv:1807.01306 [INSPIRE]].

[64] LSST DARK ENERGY SCIENCE collaboration, *A unified pseudo-$C_\ell$ framework*, *Mon. Not. Roy. Astron. Soc.* **484** (2019) 4127 [arXiv:1809.09603 [INSPIRE]].

[65] M.C. Jones and A. Pewsey, *Sinh-arcsinh distributions*, *Biometrika* **96** (2009) 761.

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

[67] A. Zonca et al., *Healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in Python*, *J. Open Source Softw.* **4** (2019) 1298 [INSPIRE].

[68] J.D. Hunter, *Matplotlib: A 2D Graphics Environment*, *Comput. Sci. Eng.* **9** (2007) 90 [INSPIRE].

[69] C.R. Harris et al., *Array programming with NumPy*, *Nature* **585** (2020) 357 [arXiv:2006.10256 [INSPIRE]].

[70] P. Virtanen et al., *SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python*, *Nature Meth.* **17** (2020) 261 [arXiv:1907.10121 [INSPIRE]].

[71] M. Kamionkowski, A. Kosowsky and A. Stebbins, *Statistics of cosmic microwave background polarization*, *Phys. Rev. D* **55** (1997) 7368 [astro-ph/9611125 [INSPIRE]].

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