# Statistical description of dust polarized emission from the diffuse interstellar medium – A RWST approach

Bruno Regaldo-Saint Blancard, François Levrier, Erwan Allys, Elena Bellomi, François Boulanger

► To cite this version:

Bruno Regaldo-Saint Blancard, François Levrier, Erwan Allys, Elena Bellomi, François Boulanger. Statistical description of dust polarized emission from the diffuse interstellar medium – A RWST approach. Astronomy & Astrophysics - A&A, 2020, 642, pp.A217. 10.1051/0004-6361/202038044. hal-02914518

HAL Id: hal-02914518

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

Submitted on 23 Aug 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Statistical description of dust polarized emission from the diffuse interstellar medium

## A RWST approach*

B. Regaldo-Saint Blancard<fnref n="1" />,<fnref n="2" />, F. Levrier<fnref n="1" />, E. Allys<fnref n="1" />, E. Bellomi<fnref n="1" />,<fnref n="2" />, and F. Boulanger<fnref n="1" />

Received 27 March 2020 / Accepted 27 August 2020

## ABSTRACT

The statistical characterization of the diffuse magnetized interstellar medium (ISM) and Galactic foregrounds to the cosmic microwave background (CMB) poses a major challenge. To account for their non-Gaussian statistics, we need a data analysis approach capable of efficiently quantifying statistical couplings across scales. This information is encoded in the data, but most of it is lost when using conventional tools, such as one-point statistics and power spectra. The wavelet scattering transform (WST), a low-variance statistical descriptor of non-Gaussian processes introduced in data science, opens a path towards this goal. To establish the methodology, we applied the WST to noise-free maps of dust polarized thermal emission computed from a numerical simulation of magnetohydrodynamical turbulence in the diffuse ISM. We analyzed normalized complex Stokes maps and maps of the polarization fraction and polarization angle. The WST yields a few thousand coefficients; some of them measure the amplitude of the signal at a given scale, and the others characterize the couplings between scales and orientations. The dependence on orientation can be fitted with the reduced wavelet scattering transform (RWST), an angular model introduced in previous works for total intensity maps. The RWST provides a statistical description of the polarization maps, quantifying their multiscale properties in terms of isotropic and anisotropic contributions. It allowed us to exhibit the dependence of the map structure on the orientation of the mean magnetic field and to quantify the non-Gaussianity of the data. We also used RWST coefficients, complemented by additional constraints, to generate random synthetic maps with similar statistics. Their agreement with the original maps demonstrates the comprehensiveness of the statistical description provided by the RWST. This work is a step forward in the analysis of observational data and the modeling of CMB foregrounds. We also release PyWST, a public Python package to perform WST and RWST analyses of two-dimensional data.

**Key words.** dust, extinction – ISM: magnetic fields – turbulence – methods: statistical – polarization

## 1. Introduction

The interstellar medium (ISM) is a beautifully complex physical system, in which gas particles and dust grains, coupled to a pervasive magnetic field, experience turbulent motions across a vast range of scales (Draine 2011; Hennebelle & Falgarone 2012). Nonspherical grains tend to align locally with their longest axis perpendicular to the Galactic magnetic field (Andersson et al. 2015; Reissl et al. 2020), leading to polarization in extinction in the visible and in emission in the submillimeter (Draine & Fraisse 2009; Guillet et al. 2018). The multiscale filamentary structure of diffuse interstellar matter is spectacularly illustrated by Herschel observations of dust emission (Miville-Deschênes et al. 2010). Observations of dust polarization provide an additional perspective because they probe the orientation of magnetic fields. Planck has provided us with the first all-sky survey of polarized dust emission, opening the path to statistical studies (Planck Collaboration XII 2020). This broad view is being complemented by observations at higher angular resolution, which are carried out by the balloon-borne experiments

BLASTPOL (Fissel et al. 2016) and PILOT (Mangilli et al. 2019), the far-IR HAWC+ camera onboard SOFIA (Chuss et al. 2019) and imaging at sub-mm/mm wavelengths from large single-dish telescopes (Ritacco et al. 2020) and ALMA (Hull et al. 2017). These observations all contribute to a common scientific goal: understanding the role turbulence and magnetic fields play along the star formation process, from the diffuse interstellar medium to molecular clouds and protostellar cores. In this ambitious endeavor, all projects face the difficulty of inferring information on magnetic fields from observations of Stokes parameters.

Over the last few years, studies of the magnetized and turbulent ISM have become closely entwined with a major goal of observational cosmology that is the detection of the primordial gravitational wave (GW) signal from the early Universe's inflationary era (Guth 1981; Linde 1982). The reason for this entanglement lies in the superposition between the expected signal from these GWs on the one hand, imprinted in the curl-like “B-mode” part of the polarization of the cosmic microwave background (CMB, Kamionkowski et al. 1997), and the polarized emission from the Galaxy on the other hand (BICEP2/Keck Array and Planck Collaborations 2015; Planck Collaboration Int. XXX 2016; Planck Collaboration XI 2020). The current

<footnote n="1">¹ Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France
e-mail: bruno.regaldo@phys.ens.fr</footnote>

<footnote n="2">² Observatoire de Paris, PSL University, Sorbonne Université, LERMA, 75014 Paris, France</footnote>

<footnote>* The public Python package PyWST is available at https://github.com/bregaldo/pywst</footnote>
constraints on the cosmological B-mode tensor-to-scalar ratio $r < 0.06$ (95% confidence level, BICEP2 Collaboration 2018) are expected to be significantly improved by the next generation of CMB experiments, such as ACT (Naess et al. 2014), SPIDER (Fraisse et al. 2013), LiteBIRD (Ishino et al. 2016), PICO (Sutin et al. 2018), the Simons Observatory (Ade et al. 2019), and CMB-S4 (Abazajian et al. 2019), with a detection limit goal of $r \approx 10^{-3}$. However, any claim to the detection of the cosmological B-mode signal will have to be critically assessed against alternative explanations involving Galactic foregrounds.

An overarching challenge for studying both interstellar turbulence and CMB foregrounds is the need for a new approach to data analysis, which is able to efficiently capture statistical couplings across scales. This information is encoded in the data, but most of it is lost when one resorts to such classical tools as one-point statistics and power spectra.

To characterize magnetized interstellar turbulence, as well as to produce realistic synthetic dust polarization sky maps for CMB data analysis, we aim at a statistical description of the polarized Galactic emission maps that encompasses non-Gaussian characteristics. Recent advances in data science open up a new path towards this goal. Allys et al. (2019) provided the first astrophysical application of the wavelet scattering transform (WST), a low-variance statistical description of non-Gaussian processes (Mallat 2012) inspired by convolutional neural networks but that does not require any training stage (Bruna & Mallat 2013). They applied the WST to column density maps inferred from magnetohydrodynamical (MHD) simulations and to an *Herschel* observation of the thermal emission from Galactic dust. The physical regularity of the maps led them to introduce the reduced wavelet scattering transform (RWST), a statistical description of reduced dimensionality, obtained through a fit of the angular dependencies of the WST coefficients.

In this paper, we extend the WST and RWST analyses to maps of linearly polarized emission, which is represented by Stokes Q and U maps<fnref n="1" />. We begin by presenting the properties of these Stokes maps and define three polarization variables of interest built from the full (I, Q, U) maps. We then recall the definition of the WST for real maps and extend it to complex maps. In order to focus on the properties of this new statistical description of polarization data, we work with data sets built from a noise-free MHD simulation, from which we compute simulated Stokes maps. We identify regular patterns in the angular dependencies of the WST coefficients for these maps, leading us to define a RWST model for polarization maps. We give interpretations of the RWST coefficients that highlight the impact of the orientation of the mean magnetic field on the statistical properties of polarization maps. We also show how these coefficients quantitatively exhibit the non-Gaussianity of these maps. We finally assess the exhaustiveness of RWST descriptions of polarization maps by generating random maps from these statistical coefficients.

The paper is organized as follows: in Sect. 2, we recall the basic properties of Stokes maps of polarized thermal emission from dust and introduce the set of polarization variables on which we define the WST. In Sect. 3, we present the MHD simulation and the polarization maps derived from it, exhibit regularities in the WST coefficients for these maps, and define

the RWST model. In Sect. 4, we give interpretations for a subset of the RWST coefficients. We present synthetic realizations of random polarization maps based on the RWST coefficients in Sect. 5, and summarize our conclusions in Sect. 6. We also provide a public Python package to perform WST and RWST analyses of two-dimensional data called PyWST<fnref n="2" />.

## 2. Statistical description of polarization maps with the WST

We first discuss the properties of Stokes I, Q, U maps of the polarized emission of dust before defining convenient transformations of these observables. We then explain how to apply the WST to these transformed maps.

### 2.1. Properties of I, Q, U maps

In the diffuse interstellar medium and at the frequencies of a few hundred GHz, we observe the thermal emission of large grains of dust (typically grains of radii greater than 50 nm, see e.g., Draine 2011). These large grains are in thermal equilibrium, absorbing the light of surrounding stars and emitting a black body radiation at typical temperatures of a few tens of K. This emission is polarized due to the statistical alignment between the orientation of aspherical dust grains and the local magnetic field (Andersson et al. 2015), leading to preferential directions of absorption and emission of light, and thus to a polarization of the signal (Planck Collaboration Int. XIX 2015).

Taking these processes into account, a physical model of radiative transfer has been built (see derivation and references in Planck Collaboration Int. XX 2015). This model provides analytical expressions at a given submillimeter frequency $v$ for the Stokes parameters I, Q, and U:

$$I = \int S_{\nu} e^{-\tau_{\nu}} \left[ 1 - p_{0} \left( \cos^{2} \gamma - \frac{2}{3} \right) \right] d\tau_{\nu}, \quad (1)$$

$$Q = \int p_{0} S_{\nu} e^{-\tau_{\nu}} \cos(2\phi) \cos^{2} \gamma d\tau_{\nu}, \quad (2)$$

$$U = \int p_{0} S_{\nu} e^{-\tau_{\nu}} \sin(2\phi) \cos^{2} \gamma d\tau_{\nu}, \quad (3)$$

where $S_{\nu}$ is the source function which is assumed here to be that of a black body so that $S_{\nu} = B_{\nu}(T_{d})$ with $T_{d}$ the dust temperature, $\tau_{\nu}$ is the optical depth, $p_{0}$ is an intrinsic polarization fraction parameter, $\gamma$ is the angle that the local magnetic field makes with the plane of the sky, and $\phi$ is the angle that the projection of the local magnetic field on the plane of the sky makes with the $x$ direction when $z$ is the line of sight in the HEALPix<fnref n="3" /> (Górski et al. 2005) convention (see Fig. 1).

In the following, we assume uniform values for the intrinsic polarization $p_{0}$ and the dust temperature $T_{d}$, an optically thin medium, and an infinitesimal optical depth equal to $d\tau_{\nu} = \sigma_{\nu}n_{H}dz$ where $\sigma_{\nu}$ is the dust opacity (assumed uniform) and $n_{H}$ is the total gas density. In this case, Eqs. (1)–(3) become:

$$I = S_{\nu} \sigma_{\nu} \left[ N_{H} - p_{0} \int n_{H} \left( \cos^{2} \gamma - \frac{2}{3} \right) dz \right], \quad (4)$$

$$Q = p_{0} \sigma_{\nu} S_{\nu} \int n_{H} \cos(2\phi) \cos^{2} \gamma dz, \quad (5)$$

<footnote n="1">¹ The last Stokes parameter $V$ measuring the intensity of circularly polarized light is ignored in this paper as it is expected to be negligible in the diffuse interstellar medium at frequencies of a few hundred GHz (Siebenmorgen et al. 2014).</footnote>

<footnote n="2">² https://github.com/bregaldo/pywst</footnote>

<footnote n="3">³ https://healpix.jpl.nasa.gov</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Definition of the angles $\gamma$ and $\phi$ related to the orientation of the local magnetic field with respect to the plane of the sky when the z-axis corresponds to the line of sight (adapted from Planck Collaboration Int. XX 2015). The line of sight points away from the observer and $\phi$ is counted positively clockwise from the x direction in the HEALPix convention.</figcaption></figure>

$$U = p_0 \sigma_v S_v \int n_H \sin(2\phi) \cos^2 \gamma dz. \qquad (6)$$

From these formulæ, Stokes parameters $I$, $Q$, and $U$ appear to be proportional to an integration along the line of sight of quantities depending on the same three variables: the total gas density $n_H$, and the angles $\gamma$ and $\phi$. Therefore we expect a priori some statistical dependencies between $I$, $Q$, and $U$ maps, and a statistical description of these Stokes maps should take into account these dependencies.

These expressions also underline how closely related Stokes $Q$ and $U$ are. Indeed they are defined with respect to a given reference direction (see for example Landi Degl'innocenti & Landolfi 2004), and any rotation of this reference direction by a given angle $\varphi$ transforms $Q$ and $U$ into rotated quantities $Q'$ and $U'$ as follows<fnref n="4" />:

$$Q' = \cos(2\varphi)Q + \sin(2\varphi)U, \qquad (7)$$

$$U' = -\sin(2\varphi)Q + \cos(2\varphi)U. \qquad (8)$$

These relations can be written in a more compact form introducing complex quantities:

$$Q' + iU' = (Q + iU)e^{-2i\varphi}. \qquad (9)$$

As a consequence, the global phase of $Q + iU$ complex maps is directly related to the reference frame in which we define the angular variables measuring the orientation of the local magnetic field within the plane of the sky. Because of the previous transformation properties, the complex variable $Q + iU$ is more apt to represent linear polarization than $Q$ and $U$ separately (see e.g., Haverkorn et al. 2004). Moreover it will be useful to define a statistical description of polarization maps that is independent from this choice of reference frame as we discuss in Sect. 2.3.

Finally, let us discuss the non-Gaussianity of Stokes maps $I$, $Q$, and $U$. We know that the gas and the magnetic field in the diffuse interstellar medium are highly turbulent (Hennebelle & Falgarone 2012), and the nonlinearity involved in the MHD equations that describe their dynamics couples scales. Therefore we expect the statistical properties of the sky maps of $I$, $Q$, and $U$ to be highly non-Gaussian. Unless one can find a mapping of the data that leads to Gaussian statistics, a mere statistical description in terms of power spectra measurements cannot be sufficient to statistically characterize these maps.

## 2.2. Transformation of Stokes maps

Stokes maps $I$ and $Q + iU$ are raw observables of the polarized emission of dust that we first want to transform. Our goal is twofold: to simplify the statistical properties of these observables and to ease the physical interpretation of their properties. In this work, we consider three dimensionless variables based on Stokes maps $I$ and $Q + iU$ to build and interpret a statistical description of polarized light: normalized Stokes variable $\tilde{Q} + i\tilde{U}$, polarization fraction $p$ and complex polarization angle $\exp(2i\psi)$.

We first define the normalized Stokes complex variable $\tilde{Q} + i\tilde{U}$ as follows:

$$\tilde{Q} + i\tilde{U} = \frac{Q + iU}{I + P}, \qquad (10)$$

with $P = |Q + iU|$ the polarized intensity. Such a definition roughly disentangles the structure of the magnetic field and the structure of dust density at the lowest order. Indeed, assuming a constant orientation of the magnetic field along the line of sight, one gets using Eqs. (4)–(6):

$$\tilde{Q} + i\tilde{U} = \frac{3p_0}{3 + 2p_0} \cos^2 \gamma \exp(2i\phi), \qquad (11)$$

which is independent of the density field $n_H$.

The polarization fraction $p$ and the polarization angle $\psi$ follow the usual definitions:

$$p = \frac{|Q + iU|}{I}, \qquad (12)$$

and

$$\psi = \frac{1}{2} \operatorname{atan2}(U, Q), \qquad (13)$$

where atan2(b, a) is the function that returns the angle restricted to $(-\pi, \pi]$ between the positive x-axis and the ray to the point (a, b) in the Euclidean plane. Because of the nonlinearity of the atan2 function, we expect the statistics of $\psi$ maps to be unnecessarily complicated compared to those of the original Stokes maps. We thus choose to analyze maps of the complex variable $\exp(2i\psi)$. This variable also avoids dealing with discontinuities of $\psi$ maps that appear where the polarization angle approaches $\pm\pi/2$. Moreover, the statistical properties of $\exp(2i\psi)$ are easier to compare to those of $\tilde{Q} + i\tilde{U}$ since $\tilde{Q} + i\tilde{U} \approx p \exp(2i\psi)$ when $P \ll I$ (which is typically the case).

## 2.3. WST statistical description

The WST is a tool from data science that computes statistical descriptions of images that can efficiently discriminate between images which have identical power spectra but distinct higher order moments (Bruna & Mallat 2013). An advantage of this description compared to the direct estimation of higher order moments is that the estimators of the WST coefficients have a low variance and thus do not need a large amount of samples for an accurate estimation. This is because the WST relies on nonexpansive operators (see Bruna & Mallat 2013). Originally designed to understand the properties of deep convolutional neural networks which have proved successful for classification problems (LeCun et al. 2010; Krizhevsky et al. 2012), this tool also provides a convenient statistical description of astronomical observations of the interstellar medium for total intensity

<footnote n="4">⁴ These relations hold for an angle $\varphi$ that is counted positively clockwise from the x direction when the z-axis is the line of sight in the HEALPix convention. See Fig. 1.</footnote>
maps (Allys et al. 2019). In this subsection we use the WST to describe the statistical properties of the polarization related maps discussed in the previous section: complex maps $\tilde{Q} + i\tilde{U}$ and $\exp(2i\psi)$, and real map $p$.

Let us briefly recall the definition of the WST of a given 2D real field $I(\mathbf{x})$ that is statistically homogeneous<fnref n="5" />. For a more complete presentation we refer to Allys et al. (2019) and Bruna & Mallat (2013). The WST relies on convolutions of the target field $I(\mathbf{x})$ with a set of $J \times \Theta$ Morlet wavelets $\psi_{j,\theta}$ with $0 \le j \le J-1$ and $\theta \in \{k\pi/\Theta, 0 \le k \le \Theta-1\}$, with $J$ and $\Theta$ two integers. These wavelets $\psi_{j,\theta}$ are the result of the dilation by a factor $2^j$ and the rotation by an angle $\theta$ of a mother Morlet wavelet $\Psi$:

$$ \psi_{j,\theta}(\mathbf{x}) = 2^{-2j} \Psi(2^{-j} r_{\theta}^{-1} \mathbf{x}). \quad (14) $$

A Morlet wavelet is essentially a plane wave modulated by a Gaussian envelope. It is a quite general kind of wavelet to study physical fields and it provides a good angular selectivity (Farge 1992). We recall the definition of the mother Morlet wavelet $\Psi$:

$$ \Psi(\mathbf{x}) = \alpha (e^{i k_0 \cdot \mathbf{x}} - \beta) e^{-|\mathbf{x}|^2 / 2\sigma^2}, \quad (15) $$

with $\alpha$ and $\beta$ two constants that are adjusted to ensure a zero mean and a unit $L^1$ norm, $k_0 = k_0 e_x$ the wave vector of the plane wave factor, and $\sigma$ the standard deviation of the Gaussian envelope<fnref n="6" />,<fnref n="7" />. Convolving a field $I(\mathbf{x})$ with a wavelet $\psi_{j,\theta}$ corresponds to a local band-pass filtering of the field at frequencies centered on a mode $k_{j,\theta} = k_0 2^{-j} (\cos(\theta) \hat{e}_x + \sin(\theta) \hat{e}_y)$. $J$ and $\Theta$ constants thus correspond respectively to the number of scales and angles that discretize the 2D Fourier space. For a statistically homogeneous real field $I(\mathbf{x})$, normalized WST coefficients are defined as follows:

$$ \bar{S}_0 = \langle I \rangle, \quad (16) $$

$$ \bar{S}_1(j_1, \theta_1) = \frac{\langle |I \star \psi_{j_1, \theta_1}| \rangle}{\langle I \rangle}, \quad (17) $$

with $0 \le j_1 \le J-1$ and $0 \le \theta_1 < \pi$,

$$ \bar{S}_2(j_1, \theta_1, j_2, \theta_2) = \frac{\langle |I \star \psi_{j_1, \theta_1} | \star \psi_{j_2, \theta_2} | \rangle}{\langle |I \star \psi_{j_1, \theta_1} | \rangle}, \quad (18) $$

with $0 \le j_1, j_2 \le J-1$ such that $j_2 > j_1$ and $0 \le \theta_1, \theta_2 < \pi$, where the averaging operator $\langle \cdot \rangle$ here is simply an average over a map, meaning $\langle I \rangle = \frac{1}{A} \int I(\mathbf{x}) d^2\mathbf{x}$ with $A$ the area of integration, and where the $\star$ symbol stands for the convolution operator. The number of WST coefficients depends on the $J$ and $\Theta$ parameters: we have a single $\bar{S}_0$ coefficient, $J \times \Theta \bar{S}_1$ coefficients, and $\Theta^2 \times J \times (J-1)/2 \bar{S}_2$ coefficients.

These coefficients can be interpreted in the following manner: $\bar{S}_0$ is simply the mean of the field, $\bar{S}_1(j_1, \theta_1)$ is a measure of the amplitude of oscillation of the normalized $I/\langle I \rangle$ field for modes centered on $k_{j_1, \theta_1}$, and finally $\bar{S}_2(j_1, \theta_1, j_2, \theta_2)$ characterizes how the normalized amplitude of oscillation of the field for a mode $k_{j_1, \theta_1}$ is modulated by a mode of oscillation $k_{j_2, \theta_2}$. Accordingly, the $\bar{S}_1(j_1, \theta_1)$ coefficient characterizes the amplitude at a

single oriented scale ($j_1, \theta_1$), while the $\bar{S}_2(j_1, \theta_1, j_2, \theta_2)$ coefficient measures the coupling between two oriented scales ($j_1, \theta_1$) and ($j_2, \theta_2$).

We note that even if the $\bar{S}_1$ coefficients characterize the Fourier amplitude of the field under study in the spectral band of the wavelets, they differ in practice from the power spectrum. While the power spectrum can be computed from $L^2$ norms of the wavelet transform, the $\bar{S}_1$ coefficients are computed with $L^1$ norms. One can however recover the power spectrum at a given scale from a quadratic sum of $\bar{S}_1$ and $\bar{S}_2$ coefficients<fnref n="8" />. Notably, we expect a non-Gaussian field to have higher $\bar{S}_2$ coefficients compared to those of a Gaussian field with identical power spectrum, and this should be counterbalanced by smaller $\bar{S}_1$ coefficients. This is related to the sparsity of the wavelet representation of the data (see Allys et al. 2019 and Bruna & Mallat 2013 for further discussions).

We can extend the previous definition of normalized WST coefficients to statistically homogeneous complex fields such as $\tilde{Q} + i\tilde{U}$ or $\exp(2i\psi)$ as follows:

$$ \bar{S}_0 = \langle |\tilde{Q} + i\tilde{U}| \rangle, \quad (19) $$

$$ \bar{S}_1(j_1, \theta_1) = \frac{\langle |(\tilde{Q} + i\tilde{U}) \star \psi_{j_1, \theta_1}| \rangle}{\langle |\tilde{Q} + i\tilde{U}| \rangle}, \quad (20) $$

$$ \bar{S}_2(j_1, \theta_1, j_2, \theta_2) = \frac{\langle |(\tilde{Q} + i\tilde{U}) \star \psi_{j_1, \theta_1} | \star \psi_{j_2, \theta_2}| \rangle}{\langle |(\tilde{Q} + i\tilde{U}) \star \psi_{j_1, \theta_1}| \rangle}, \quad (21) $$

where we still have $0 \le j_1, j_2 \le J-1$ with $j_2 > j_1$ for $\bar{S}_2$ coefficients, but here the main difference is that $0 \le \theta_1 < 2\pi$ and $0 \le \theta_2 < \pi$. This new range of angles for $\theta_1$ is directly related to the complex nature of the $\tilde{Q} + i\tilde{U}$ variable. For a real signal $I(\mathbf{x})$, one can easily show that $|I \star \psi_{j_1, \theta_1+\pi}| = |I \star \psi_{j_1, \theta_1}|$, which explains the range of $\theta_1$ in Eq. (17). For a complex signal the previous relation does not hold anymore, and we have to consider rotations of the mother wavelet $\Psi$ for angles between 0 and $2\pi$ (see Mallat 2012, for more details). We note that the range of angles for $\theta_2$ is unchanged because $|(\tilde{Q} + i\tilde{U}) \star \psi_{j_1, \theta_1}|$ is a real signal. Thus, for a complex field, we end up with twice as many $\bar{S}_1$ and $\bar{S}_2$ coefficients as for the WST of a real field.

Let us note that this definition of the WST on complex polarization maps $\tilde{Q} + i\tilde{U}$ does not depend on the global phase of $\tilde{Q} + i\tilde{U}$ maps, meaning that the WST coefficients associated to $\tilde{Q} + i\tilde{U}$ maps do not depend on the reference frame mentioned in Sect. 2.1.

In the following, the statistical properties of $p$ maps will be described with the WST coefficients defined in Eqs. (16)–(18), while the statistical properties of $\tilde{Q} + i\tilde{U}$ and $\exp(2i\psi)$ maps will be described with the WST coefficients defined in Eqs. (19)–(21).

### 3. Simplification of the WST statistical descriptions for simulated Stokes maps with the RWST

The WST coefficients of polarization variables $\tilde{Q} + i\tilde{U}$, $\exp(2i\psi)$, and $p$ define statistical descriptions of polarization maps. We expect some form of regularity in these coefficients that should lead to possible simplifications. This was the case for total intensity maps, which led to defining a reduced set of coefficients (Allys et al. 2019). In this section, we similarly define the reduced wavelet scattering transform (hereafter RWST) to characterize the statistics of polarization data.

<footnote n="5">⁵ In this work, we restrict our analysis to statistically homogeneous data for simplicity, but this is not a limitation of the methodology as the WST can also be computed locally (see e.g., Allys et al. 2019).</footnote>

<footnote n="6">⁶ In practice we choose $k_0 = 3\pi/4$ pixel$^{-1}$ and $\sigma = 0.8$ pixel as in Bruna & Mallat (2013).</footnote>

<footnote n="7">⁷ Actually, the Gaussian envelope has an elliptical shape that is not included in Eq. (15) for simplicity. This elliptical shape enhances the angular selectivity of the wavelets.</footnote>

<footnote n="8">⁸ See Eq. (9) in Allys et al. (2019).</footnote>
We first present the simulated sets of polarization maps we built from a numerical simulation of the diffuse interstellar medium. Then, we show what kind of statistical regularity arises from the WST descriptions of these data sets. We finally define RWST descriptions of these maps which are based on the modeling of the angular dependencies of the WST coefficients.

## 3.1. Building simulated Stokes maps

Throughout this paper we work with Stokes maps which are directly computed as integrations along a given direction of 3D simulated cubes of data extracted from a numerical simulation of magnetized interstellar turbulence (this given direction of integration corresponds to the line of sight for observations). As a first step, working with MHD simulations instead of observational data is a stepping-stone to relate our statistical descriptions to physics where we avoid the difficulty of dealing with data noise.

We use a MHD simulation designed to study the biphasic nature of the diffuse interstellar medium (Bellomi et al., in prep.). The simulation uses the adaptive mesh refinement code RAMSES (Teyssier 2002; Fromang et al. 2006) to solve the equations of ideal MHD as described in Iffrig & Hennebelle (2017), neglecting self-gravity and taking into account heating and cooling processes of the gas. Turbulent forcing is applied to inject kinetic energy and balance numerical dissipation so that the simulation reaches a statistical steady state. In practice this turbulent forcing consists in a large-scale stochastic force field (a three component Ornstein-Uhlenbeck process defined in spectral space). Details on the turbulent forcing model can be found in Schmidt et al. (2006).

The simulation volume is a (50 pc)³ box divided into 512³ cells with periodic boundary conditions. At $t=0$ the gas has uniform properties with a density $n_H = 1.5 \text{ cm}^{-3}$ and a temperature $T = 8000$ K, and the magnetic field is also uniform, with $B_0 = B_0 e_x$ and $B_0 \sim 3.8 \mu\text{G}$. In steady state, the turbulent forcing leads to an approximate velocity dispersion $\sigma_v \sim 2.6 \text{ km s}^{-1}$. This gives a turnover time at large scale $\tau_L \approx 18.8$ Myr. Finally, an isotropic Habing radiation field is applied at the boundaries of the box. Its intensity is scaled by a factor $G_0 = 1$ (for a definition, see Draine 2011).

Once the simulation has reached a statistical steady state, we extract 14 snapshots that are approximately statistically independent using an approach similar to Federrath et al. (2009). In practice we make sure that two consecutive snapshots correspond to a minimal time evolution of $\delta\tau = 1.25$ Myr which is roughly ten percent of $\tau_L$. The phenomenology of turbulence in the sense of Kolmogorov shows that the turnover time $\tau_L$ at a given scale $l$ scales as $l^{2/3}$ (Frisch 1995). This scaling holds for incompressible hydrodynamical turbulence only but we assume for the sake of simplicity that it extends reasonably well to ideal compressible MHD. At the range of scales considered in the following, we find that $\delta\tau$ is about fives times smaller than the corresponding largest value of $\tau_L$. Although this value is not completely satisfactory, we assume that the snapshots are statistically independent.

We can compute a set of Stokes maps $I$, $Q$, and $U$ for each of these snapshots using Eqs. (4)–(6) and choosing the z-axis of the cubes as the direction of integration. We note that this integration procedure is identical to the one used in Planck Collaboration Int. XX (2015). We choose a typical value for the polarization fraction parameter $p_0 = 0.2$ and arbitrary values for $\sigma_v$ and $T_d$ as these only determine the global amplitude of $I$, $Q$ and $U$ maps but do not impact the analysis. These maps are relevant for any

frequency $\nu$ provided that the dust emission remains optically thin. Then from these Stokes maps we compute the associated maps $P$, $\tilde{Q}$, $\tilde{U}$, $p$, $\cos(2\psi)$, and $\sin(2\psi)$ which are the transformations of the Stokes observables defined in Sect. 2.2. Finally, we end up with a set of 14 maps for each of these variables.

We point out that the initial conditions of the MHD simulations are anisotropic because of the initial direction of the uniform magnetic field. This anisotropy remains once the simulation has reached steady state, due to magnetic flux conservation, and the value of the mean magnetic field is $\tilde{B} \approx B_0$. Because the direction of integration chosen to compute the previous Stokes maps is orthogonal to the direction of the mean magnetic field $\tilde{B}$ in the simulation, we refer to the previous sets of maps using the $\perp$ symbol. For instance, the $p_\perp$ data set refers to the set of 14 maps of polarization fraction $p$ computed from the 14 maps $I$, $Q$, and $U$ for which the z-axis was the axis of integration. We similarly compute Stokes maps integrating along the x-axis which is the direction of the mean magnetic field. This results in a set of maps that are statistically isotropic and from which we also compute the associated maps $P$, $\tilde{Q}$, $\tilde{U}$, $p$, $\cos(2\psi)$, and $\sin(2\psi)$. We use the $\parallel$ symbol to refer to these sets of maps.

## 3.2. Presentation of the data sets

We use for this work eight data sets, each one comprising 14 statistically independent maps. The various maps presented in the last subsection define 6 data sets: $\exp(2i\psi)_\perp$, $p_\perp$, $\tilde{Q}_\perp + i\tilde{U}_\perp$, $\exp(2i\psi)_\parallel$, $p_\parallel$, and $\tilde{Q}_\parallel + i\tilde{U}_\parallel$. In addition, we also analyze phase randomized data sets $R[\tilde{Q}_\perp + i\tilde{U}_\perp]$ and $R[\tilde{Q}_\parallel + i\tilde{U}_\parallel]$ built respectively from the $\tilde{Q}_\perp + i\tilde{U}_\perp$ and $\tilde{Q}_\parallel + i\tilde{U}_\parallel$ data sets.

We define $R[\cdot]$ to be the operator that acts on a map by randomizing the phases of the map in Fourier space, meaning that the new phases are drawn from a uniform distribution on $[0, 2\pi)$. We note that the modulus of the Fourier coefficients are retained so that the power spectrum is also unchanged. The phase information of an image is tightly bound to its structure (Oppenheim & Lim 1981) and the main effect of the $R[\cdot]$ operator is to severely damage the structure of the image. We use phase randomization as an approximation to Gaussianization<fnref n="9" />. We could have used a standard Gaussian random field generation (see e.g., Kroese & Botev 2015), however we found out that the two approaches give similar results. In practice $R[\tilde{Q}_\perp + i\tilde{U}_\perp]$ (respectively $R[\tilde{Q}_\parallel + i\tilde{U}_\parallel]$) refers to the set of 14 maps produced by randomizing separately $\tilde{Q}$ and $\tilde{U}$ maps from the $\tilde{Q}_\perp + i\tilde{U}_\perp$ data set (respectively $\tilde{Q}_\parallel + i\tilde{U}_\parallel$). Technical details about phase randomization can be found in Appendix B, as well as an example of a phase randomized map for $R[\tilde{Q}_\perp + i\tilde{U}_\perp]$ in Fig. B.1.

Figure 2 shows examples of maps for $I_\perp$, $Q_\perp$, $U_\perp$, $P_\perp$, $\tilde{Q}_\perp$, $\tilde{U}_\perp$, $p_\perp$, $\cos(2\psi)_\perp$, and $\sin(2\psi)_\perp$ data sets, while Fig. 3 shows examples of maps for the corresponding $\parallel$ data sets. These figures show maps that all rely on the same snapshot of a MHD simulation. We draw the attention of the reader to a few points. First, we clearly see filamentary patterns on these maps that will demand a statistical description involving higher orders statistics compared to simple power spectra. We also note that $P$ is an order of magnitude lower than the intensity $I$ (this is due to the value of the polarization fraction parameter $p_0$), hence we have approximately $\tilde{Q} + i\tilde{U} \approx (Q + iU)/I$ so that $p$ roughly behaves as the modulus of $\tilde{Q} + i\tilde{U}$, and $2\psi$ roughly behaves as its complex argument. Next, we note that the magnitude of $p_\parallel$ is

<footnote n="9">⁹ Stationary Gaussian random fields (GRFs) do have uniformly distributed phases on $[0, 2\pi)$ but this property alone does not define them (see Wandelt 2013).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Examples of $I_{\perp}$, $Q_{\perp}$, $U_{\perp}$, $P_{\perp}$, $\tilde{Q}_{\perp}$, $\tilde{U}_{\perp}$, $p_{\perp}$, $\cos(2\psi)_{\perp}$, and $\sin(2\psi)_{\perp}$ maps (from top to bottom and left to right) that are built from a given snapshot of the MHD simulation.</figcaption></figure>

unsurprisingly much lower than that of $p_{\perp}$ as the direction of the projection of the local magnetic field on the plane of the sky is much less coherent along the line of sight when the line of sight corresponds to the direction of the mean magnetic field. Finally, we see on the $\cos(2\psi)_{\perp}$ map that the anisotropy of the magnetic field in the simulation leads to values that are concentrated close to 1.

### 3.3. Regularity in the WST coefficients

We now focus on the WST coefficients associated with the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ data set, but the following reasoning remains valid for the other data sets, including those related to $\exp(2i\psi)$ and $p$.

For the 14 maps of the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ data set, we compute the WST coefficients for $J=7$ and $\Theta=8$ (see Eqs. (19)–(21)). $J$ could have been fixed to a higher value for 512×512 maps, but because of our limited number of maps we chose to restrict our statistical description to scales for which we have a sufficient number of modes for reliable estimations. We give in Table 1 the correspondence between the scale index $j$ and the central wavelength of the related dilated Morlet wavelet both in pixel units and in dimensional units (related to the MHD simulation).

<table><caption>**Table 1.** Correspondence between scale index $j$ and related wavelengths on simulated maps, both in pixel units and dimensional units.</caption><thead><tr><th>j</th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>λ̃ [pixels]</td><td>2.67</td><td>5.33</td><td>10.7</td><td>21.3</td><td>42.7</td><td>85.3</td><td>171</td></tr><tr><td>λ [pc]</td><td>0.26</td><td>0.52</td><td>1.04</td><td>2.08</td><td>4.17</td><td>8.33</td><td>16.7</td></tr></tbody></table>

**Notes.** Those wavelengths $\tilde{\lambda}$ come from the definition of dilated Morlet wavelets $\psi_{j,\theta}$ (for more details, see Bruna & Mallat 2013).

These $J$ and $\Theta$ values lead for each map to 112 $\tilde{S}_1$ coefficients and 2688 $\tilde{S}_2$ coefficients. Adding to this the $\tilde{S}_0$ coefficient, we end up with 2801 coefficients per map. We define mean $\tilde{S}_0$, $\tilde{S}_1$ and $\tilde{S}_2$ coefficients as means of the WST coefficients over the 14 maps and we also compute the standard deviation of the mean for each of these coefficients assuming that the maps are statistically independent.

Figures 4 and 5 represent (in blue) respectively $\tilde{S}_1(j_1, \theta_1)$ coefficients and a representative subsample of $\tilde{S}_2(j_1, \theta_1, j_2, \theta_2)$ coefficients (for $j_1=1$) on a logarithmic scale, for $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Same as Fig. 2 but for the $I_\parallel$, $Q_\parallel$, $U_\parallel$, $P_\parallel$, $\tilde{Q}_\parallel$, $\tilde{U}_\parallel$, $p_\parallel$, $\cos(2\psi)_\parallel$, and $\sin(2\psi)_\parallel$ data sets.</figcaption></figure>

In both figures, we represent the multivariate functions $\tilde{S}_1$ and $\tilde{S}_2$ in a lexicographical order for the ($j_1, \theta_1$) and ($j_1, \theta_1, j_2, \theta_2$) variables, respectively. Vertical gray lines help to mark increments of these variables.

In Fig. 4, we see that for a fixed scale $j_1$ a smooth pattern emerges with respect to the angular variable $\theta_1$. Because of the definition of $\tilde{S}_1$ coefficients, at fixed $j_1$, $\log_2(\tilde{S}_1)(j_1, \cdot)$ functions must be $2\pi$-periodic, but here the smooth pattern appears to be $\pi$-periodic. While it is not surprising to get smooth angular patterns that reflect the regularity of physical processes, this $\pi$-periodicity with respect to the angular variable $\theta_1$ was unexpected.

We can nevertheless explain it by noticing that:

$$|(\tilde{Q} + i\tilde{U}) \star \psi_{j,\theta}| = |[(\tilde{Q} + i\tilde{U}) \star \psi_{j,\theta}]^*| = |(\tilde{Q} - i\tilde{U}) \star \psi_{j,\theta+\pi}|, \quad (22)$$

because for Morlet wavelets we have $\psi_{j,\theta}^* = \psi_{j,\theta+\pi}^{10}$. Hence, saying that $\log_2(\tilde{S}_1)(j_1, \cdot)$ functions are $\pi$-periodic amounts to saying that $\tilde{Q} + i\tilde{U}$ and $\tilde{Q} - i\tilde{U}$ have the same WST first order statistics.

We can identify the same kind of regularity properties for $\tilde{S}_2$ coefficients in Fig. 5. Since $\tilde{S}_2$ coefficients depend on four

<figure><img src="image_5.png" /><figcaption>Fig. 4. WST coefficients $\tilde{S}_1(j_1, \theta_1)$ on a logarithmic scale for the $\tilde{Q}_\perp + i\tilde{U}_\perp$ data set, presented in a lexicographical order on ($j_1, \theta_1$). Vertical dashed lines delimit distinct $j_1$ values. Top panel: original data (solid lines) and the RWST fit (dotted lines) corresponding to the model of Eq. (23), while bottom panel: normalized residuals of the fit.</figcaption></figure>

<footnote n="10">¹⁰ The * symbol stands for the complex conjugate.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5. $j_1 = 1$ selection of WST coefficients $\bar{S}_2(j_1, \theta_1, j_2, \theta_2)$ on a logarithmic scale for the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ data set presented in a lexicographical order on $(j_1, \theta_1, j_2, \theta_2)$. This specific selection of coefficients is arbitrary, and we find similar results for other scales and the other data sets. Vertical dashed and dotted lines delimit distinct $\theta_1$ and $j_2$ values, respectively. Top panel: original data (solid lines) and the RWST fit (dotted lines) corresponding to the model of Eq. (24), while bottom panel: normalized residuals of the fit.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Surface representation of WST coefficients $\bar{S}_2(j_1 = 1, \theta_1, j_2 = 3, \theta_2)$ for the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ data set as a function of $\theta_1$ and $\theta_2$ variables only.</figcaption></figure>

variables, two of which are angular variables $\theta_1$ and $\theta_2$, it is more complicated to see this angular regularity, but Fig. 6 helps us to exhibit this regularity by plotting $\bar{S}_2$ coefficients as a function of angular variables $\theta_1$ and $\theta_2$ when $j_1 = 1$ and $j_2 = 3$. Smooth angular patterns appear on this surface which is parameterized by $\theta_1$ and $\theta_2$. In particular we see that $\theta_1 - \theta_2 = c$ cuts of this surface for arbitrary constants $c$ give roughly constant $\bar{S}_2$ coefficients. On this example we thus expect most of the angular modulation to depend on the $\theta_2 - \theta_1$ variable.

### 3.4. Definition of a RWST statistical description

The smooth periodic patterns identified in the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ WST coefficients suggest that a simplification of the WST statistical description is possible, through an adequate modeling of its

angular dependencies. The purpose of such a modeling of the WST coefficients is twofold: (1) lowering the dimensionality of the statistical description of our data, and (2) providing an interpretable representation of these angular dependencies in terms of considerations of isotropy and anisotropy of the data.

We model the regularity of WST coefficients with respect to angular variables $\theta_1$ and $\theta_2$ using the RWST model introduced in Allys et al. (2019). It is remarkable that a model developed for total intensity maps may be applied to $\tilde{Q} + i\tilde{U}$ complex Stokes maps of polarized thermal emission of dust, without any modification. This highlights the generality of such an angular modeling of the WST coefficients for maps of dust emission, and this model may surely be extended to other kinds of complex-valued fields in physics. We present in Appendix A the RWST model in terms of Fourier series expansions and rephrase the characteristics of this generality.

We now briefly recall the RWST model, and refer the reader to Allys et al. (2019) for more details.

In the RWST model, the $\bar{S}_1$ coefficients are written as:

$$ \log_2[\bar{S}_1(j_1, \theta_1)] = \hat{S}_1^{\text{iso}}(j_1) + \hat{S}_1^{\text{aniso}}(j_1) \cos(2[\theta_1 - \theta^{\text{ref},1}(j_1)]), \quad (23) $$

where $\hat{S}_1^{\text{iso}}(j_1)$, $\hat{S}_1^{\text{aniso}}(j_1)$, and $\theta^{\text{ref},1}(j_1)$ are the parameters of this angular model for scale $j_1$. This model thus depends on $3 \times J$ parameters. We also enforce $\hat{S}_1^{\text{aniso}}(j_1) \ge 0$ in order to lift a degeneracy between the $\hat{S}_1^{\text{aniso}}(j_1)$ and $\theta^{\text{ref},1}$ parameters. $\hat{S}_1^{\text{iso}}(j_1)$ quantifies the isotropic component of the data fluctuations at a scale $j_1$, while $\hat{S}_1^{\text{aniso}}(j_1)$ defines a measure of the angular modulation of the coefficients at scale $j_1$, introducing a reference angle $\theta^{\text{ref},1}(j_1)$ that defines a privileged direction in the maps$^{11}$.

<footnote n="11">11 We note that Eq. (23) defines an angular model of the logarithm of the WST coefficients that turns angular modulations of the WST coefficients into additive corrections $\hat{S}_1^{\text{aniso}}(j_1) \cos(2[\theta_1 - \theta^{\text{ref},1}(j_1)])$ to the</footnote>
Similarly, the RWST model for $\tilde{S}_2$ coefficients reads:

$$ \log_2 [\tilde{S}_2(j_1, \theta_1, j_2, \theta_2)] = \hat{S}_2^{\text{iso},1}(j_1, j_2) \\ + \hat{S}_2^{\text{iso},2}(j_1, j_2) \cos(2[\theta_1 - \theta_2]) \\ + \hat{S}_2^{\text{aniso},1}(j_1, j_2) \cos(2[\theta_1 - \theta^{\text{ref},2}(j_1, j_2)]) \\ + \hat{S}_2^{\text{aniso},2}(j_1, j_2) \cos(2[\theta_2 - \theta^{\text{ref},2}(j_1, j_2)]), \quad (24) $$

where $\hat{S}_2^{\text{iso},1}(j_1, j_2)$, $\hat{S}_2^{\text{iso},2}(j_1, j_2)$, $\hat{S}_2^{\text{aniso},1}(j_1, j_2)$, $\hat{S}_2^{\text{aniso},2}(j_1, j_2)$, and $\theta^{\text{ref},2}(j_1, j_2)$ are the parameters of this angular model for each pair of scales $(j_1, j_2)$. As we have $j_2 > j_1$, we end up with $5 \times J \times (J-1)/2$ parameters for this model. Here again we make sure that $\hat{S}_2^{\text{aniso},1}(j_1, j_2) \ge 0$ to avoid any parameter degeneracy. $\hat{S}_2^{\text{iso},1}$ measures the overall amplitude of coupling between the scales $j_1$ and $j_2$. $\hat{S}_2^{\text{iso},2}$ represents the amplitude of the modulation due to the relative orientation of the wavelets $\psi_{j_1,\theta_1}$ and $\psi_{j_2,\theta_2}$, and we interpret it as a signature of filamentary structures at a given scale. Indeed for an oriented filamentary structure, we expect the $\tilde{S}_2$ coefficients to peak when $\theta_2 = \theta_1$ and to reach a minimum when $\theta_2 = \theta_1 + \pi/2$. Finally, $\hat{S}_2^{\text{aniso},1}$ and $\hat{S}_2^{\text{aniso},2}$ are measures of the anisotropic properties of the data in second order WST coefficients, here decoupling $\theta_1$ and $\theta_2$ contributions.

In practice, for a given data set, this RWST model of the angular dependencies is fitted to the first order WST coefficients for every scale $j_1$, and to the second order WST coefficients for every pair of scales $(j_1, j_2)$ (with $j_2 > j_1$). The accuracy of these multiple fits is quantified with $\chi_r^2$ statistics as described in Allys et al. (2019). Since it is not possible to properly estimate the full covariance matrix with only 14 samples per coefficient, the uncertainties affecting the WST coefficients for a given data set are simply estimated from the sample variance across the various simulation snapshots.

However an important correlation of the first order WST coefficients across angles for each scale $j_1$ needs to be taken into account to properly estimate statistical uncertainties. For each sample we compute a mean coefficient across angles for a given scale and subtract this mean before computing the statistical uncertainties<fnref n="12" />. This mitigates most of the correlation between WST coefficients at the same scale $j_1$.

Figure 7 shows the $\chi_r^2$ values for both $\log_2(\tilde{S}_1)$ and $\log_2(\tilde{S}_2)$ fits for the $\tilde{Q}_\perp + i\tilde{U}_\perp$ data set. The $\chi_r^2$ values are globally close to 1 except at low $j_1$ for $\chi_r^{2,S_1}$ and at low $j_2 - j_1$ for $\chi_r^{2,S_2}$. This deterioration of the goodnesses of the fits is due to a pixellization effect at small scales, and may be fixed by adding adequate lattice terms in the RWST model, as described in the Appendix C of Allys et al. (2019). This is shown for our data in Appendix C.1. We note that these additional terms do not affect the values of the coefficients corresponding to the original RWST model. Finally, the same RWST model applies to the other data sets used in this work, including $\exp(2i\psi)$ and $p$ data sets, and we get for all of them similar $\chi_r^2$ values.

In Figs. 4 and 5 we show the RWST fit overplotted on a selection of first and second order WST coefficients, and also show the corresponding normalized residuals. The curves are globally in good agreement and the flaws of the RWST model due to numerical effects at low $j_1$ and $j_2 - j_1$ appear as clear patterns in the residuals. Additional figures in Appendix C.1 show how these artifacts may be taken into account.

isotropic amplitude of fluctuations $\hat{S}_1^{\text{iso}}(j_1)$. We use a base 2 logarithm to be consistent with the base 2 definition of scales $j_1$ and $j_2$.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Reduced chi square $\chi_r^{2,S_1}(j_1)$ (top) and $\chi_r^{2,S_2}(j_1, j_2)$ (bottom) associated with the RWST fit of the WST coefficients (see Eqs. (23) and (24)) for the $\tilde{Q}_\perp + i\tilde{U}_\perp$ data set. Bottom panel: each curve corresponds to a fixed $j_1$ value with $j_2$ varying from $j_1 + 1$ to $J - 1 = 6$. For $j_1 = J - 2$, the curve is reduced to a single dot on the figure. We use logarithmic scales for better visibility.</figcaption></figure>

In the end, the RWST coefficients define statistical descriptions of the data sets in terms of simple considerations of isotropy and anisotropy with respect to a reference direction. Furthermore, these statistical descriptions exhibit much lower dimensionality, with a total of 127 coefficients (including $\tilde{S}_0$ coefficient in the description) compared to the original 2801 WST coefficients in the case of complex variables $\tilde{Q} + i\tilde{U}$ and $\exp(2i\psi)^{13}$, thus providing a very large compression of the statistical information contained in the WST coefficients.

# 4. RWST data analysis and interpretation

In this section we analyze and interpret the RWST statistical descriptions for the 8 data sets built from the same MHD simulation. We relate the coefficients of these descriptions to the physical properties of the simulation.

## 4.1. Isotropic fluctuations in first order coefficients $\hat{S}_1^{\text{iso}}$

Figure 8 presents $\hat{S}_1^{\text{iso}} + \log_2(\tilde{S}_0)$ coefficients as a function of scale for $\tilde{Q}_\perp + i\tilde{U}_\perp$, $\tilde{Q}_\parallel + i\tilde{U}_\parallel$, $R[\tilde{Q}_\perp + i\tilde{U}_\perp]$ and $R[\tilde{Q}_\parallel + i\tilde{U}_\parallel]$ data sets. We note that the $\log_2(\tilde{S}_0)$ term cancels the normalization of the $\tilde{S}_1$ WST coefficients defined in Eq. (20) in order to analyze the statistics of $\tilde{Q} + i\tilde{U}$ (and not the statistics of $(\tilde{Q} + i\tilde{U})/\langle|\tilde{Q} + i\tilde{U}|\rangle$). We see that $\hat{S}_1^{\text{iso}} + \log_2(\tilde{S}_0)$ levels which measure an amplitude of fluctuations of the signal, per scale, are higher for $\tilde{Q}_\perp + i\tilde{U}_\perp$ and $R[\tilde{Q}_\parallel + i\tilde{U}_\parallel]$ than for the corresponding $\parallel$ data sets. This shows that the amplitudes differ depending on the orientation of the mean magnetic field with respect to the line

<footnote n="12">¹² This decreases the effective number of degrees of freedom by one when computing $\chi_r^2$ values.</footnote>

<footnote n="13">¹³ For $p$ maps, the WST descriptions consist of 1401 coefficients.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** $\hat{S}_1^{\text{iso}}(j_1) + \log_2(\hat{S}_0)$ RWST coefficients for $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}, \tilde{Q}_{\parallel} + i\tilde{U}_{\parallel}$, $R[\tilde{Q}_{\perp} + i\tilde{U}_{\perp}]$ and $R[\tilde{Q}_{\parallel} + i\tilde{U}_{\parallel}]$ data sets. We use solid lines for $\perp$ data sets and dotted lines for $\parallel$ data sets. For reference, we show at the top of this figure the correspondence between $j$ scale indices and the related wavelengths on the maps (in practice, we have $\lambda_1 = 2^{j_1+1}\pi/k_0$).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** $\hat{S}_1^{\text{iso}}(j_1)$ RWST coefficients for $\tilde{Q} + i\tilde{U}, \exp(2i\psi)$ and $p$ data sets in the $\perp$ (top, solid lines) and $\parallel$ (bottom, dotted lines) cases.</figcaption></figure>

of sight. We have more fluctuations within $\tilde{Q} + i\tilde{U}$ maps when the mean magnetic field is in the plane of the sky.

The differences between the $\tilde{Q} + i\tilde{U}$ data set and its corresponding randomized counterpart $R[\tilde{Q} + i\tilde{U}]$ (in both $\perp$ and $\parallel$ cases) illustrate the difference between $\tilde{S}_1$ coefficients and the power spectrum. Indeed, $\tilde{Q} + i\tilde{U}$ and $R[\tilde{Q} + i\tilde{U}]$ maps share the same power spectrum but have different $\tilde{S}_1$ coefficients. The fact that the $R[\cdot]$ operator increases the $\tilde{S}_1$ values leaving the power spectrum unchanged shows that it reduces the sparsity of these maps (see discussion in Sect. 2.3). This feature underlines the non-Gaussianity of the $\tilde{Q} + i\tilde{U}$ data set and we therefore expect its second order coefficients to be higher compared to those of the corresponding randomized data set. We note that these differences wear off at large scales. We interpret this as statistical evidence that the non-Gaussianity of $\tilde{Q} + i\tilde{U}$ maps decreases at large scales. This result may reflect a characteristic of interstellar turbulence but could also follow from the fact that we start to probe the Gaussian distribution of the turbulent forcing of the simulation.

In Fig. 9, we compare $\tilde{Q} + i\tilde{U}$ $\hat{S}_1^{\text{iso}}$ coefficients to those of $p$ and $\exp(2i\psi)$ for respectively $\parallel$ and $\perp$ data sets. Since we have $\tilde{Q} + i\tilde{U} \approx p \exp(2i\psi)$ we would like to compare the relative contributions of $\hat{S}_1^{\text{iso}}$ coefficients of $\tilde{Q} + i\tilde{U}$ maps, but this seems more complicated than expected as we found out that a proper comparison through $\hat{S}_1^{\text{iso}}$ coefficients highly depends on the choice of normalization of the WST coefficients for $p$.

$\hat{S}_1^{\text{iso}}$ coefficients roughly range from -4.5 to -3.0 for $\exp(2i\psi)_{\parallel}$ while they range from -5.8 to -4.2 for $\exp(2i\psi)_{\perp}$. These differences indicate larger fluctuations of the polarization angle when the mean magnetic field is along the line of sight compared to when the mean magnetic field is in the plane of the sky. Since the average polarization fraction $p$ is lower when the mean magnetic field is along the line of sight ($\tilde{S}_0 \approx 0.03$ for $p_{\parallel}$ compared to $\tilde{S}_0 \approx 0.1$ for $p_{\perp}$) this feature is consistent with the anti-correlation observed with polarization data between the angle dispersion function $S$ and the polarization fraction $p$ (Planck Collaboration Int. XX 2015; Planck Collaboration XII 2020; Fissel et al. 2016).

## 4.2. Anisotropic fluctuations in first order coefficients $\hat{S}_1^{\text{aniso}}$

$\hat{S}_1^{\text{aniso}}$ coefficients measure the angular modulation of the first order WST coefficients. They are presented in Fig. 10 for all data sets. We see that this anisotropy is much larger for $\perp$ data sets than for $\parallel$ ones. The $p_{\perp}$ data set constitutes an exception among the $\perp$ data sets as it has a rather low anisotropy level. These differences are not surprising as we expect a stronger anisotropy when the mean magnetic field is in the plane of the sky, while statistical isotropy is expected when integrating along the mean magnetic field. For larger scales, we see an increase of these coefficients for $\perp$ data sets. This trend has already been observed on observational data of the Polaris flare in total intensity (Allys et al. 2019) and deserves a closer examination.

As explained in Sect. 3.1 the large scales of consecutive snapshots could be correlated to some extent. To assess the potential impact of these correlations on our analysis, we have computed separate RWST statistics for three maps of the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ data set corresponding to snapshots that are sufficiently distant in time to rightfully assume the independence (we choose $6 \times \delta\tau$ instead of $\delta\tau$). The increase of $\hat{S}_1^{\text{aniso}}$ coefficients remains significant for each map, which demonstrates that this trend is not a consequence of correlations among snapshots.

Reference angles $\theta^{\text{ref},1}$ also presented in Fig. 10 show that the preferential direction identified for anisotropic $\perp$ data sets is the direction corresponding to $\theta^{\text{ref},1} = 0$ for all scales $j_1$. Such a value of the reference angle indicates a statistical tendency for structures, including filaments, to be elongated vertically rather than horizontally, that is, along the $y$ axis in Fig. 2. This corresponds to an elongation which is orthogonal to the mean magnetic field. This result is to be compared in further works to the abundant literature on the relative orientation between magnetic fields and structures traced by interstellar dust (for a review, see Hennebelle & Inutsuka 2019).

The reference angle found for $\parallel$ data set is well defined and approximately equal to $\pi/4$ for all scales while the anisotropy
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** $\hat{S}_1^{\text{aniso}}(j_1)$ (top) and $\theta^{\text{ref},1}(j_1)$ (bottom) RWST coefficients for $\tilde{Q} + i\tilde{U}$, $\exp(2i\psi)$, $p$, and $R[\tilde{Q} + i\tilde{U}]$ data sets in the $\perp$ (solid lines) and $\parallel$ (dotted lines) cases.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** $\hat{S}_2^{\text{iso},1}(j_1, j_2)$ (top) and $\hat{S}_2^{\text{iso},2}(j_1, j_2)$ (bottom) RWST coefficients for $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ and $R[\tilde{Q}_{\perp} + i\tilde{U}_{\perp}]$ data sets. Each curve corresponds to a fixed $j_1$ value with $j_2$ varying from $j_1 + 1$ to $J - 1 = 6$. For $j_1 = J - 2$, the curve is reduced to a single dot on the figure.</figcaption></figure>

levels are close to zero. These values of $\theta^{\text{ref},1}$ are surprising because we were not expecting any anisotropy for these data sets. By examining the RWST statistics separately for each map, for each of the corresponding data sets, we found out that these surprising values correspond to an intermittent rise of the anisotropy level that appears in a few consecutive snapshots of the simulation. In this case where the level and direction of anisotropy are not coherent over snapshots, the mean coefficient gives an incomplete view of the anisotropic properties of the simulation. Notably, even with significant levels of anisotropy per map, if the reference angles are incoherent, we expect the mean level of anisotropy to be small. We found out that this is what happens for the $p_{\perp}$ data set, where $\hat{S}_1^{\text{aniso}}$ coefficients are relatively small while $\theta^{\text{ref},1}$ coefficients clearly deviate from zero.

### 4.3. Second order coefficients and non-Gaussianity of the data

$R[\tilde{Q} + i\tilde{U}]$ maps are Gaussian approximations of $\tilde{Q} + i\tilde{U}$ maps, and we have already exhibited differences between these data sets in their first order RWST coefficients in Sect. 4.1. Similarly, second order RWST coefficients show remarkable differences that highlight the non-Gaussianity of the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ and $\tilde{Q}_{\parallel} + i\tilde{U}_{\parallel}$ data sets. The two dominant second order RWST coefficients $\hat{S}_2^{\text{iso},1}$ and $\hat{S}_2^{\text{iso},2}$ presented in Fig. 11 display clearly distinct patterns between the original data sets and the randomized ones on the example of the $\perp$ data sets. First, $\hat{S}_2^{\text{iso},1}$ and $\hat{S}_2^{\text{iso},2}$ coefficients are globally smaller for the $R[\tilde{Q} + i\tilde{U}]$ data sets, which is in line with what we had foreseen in Sect. 4.1. In addition, $\hat{S}_2^{\text{iso},1}$ coefficients for $R[\tilde{Q} + i\tilde{U}]$ show a scale invariance property: $\hat{S}_2^{\text{iso},1}$ coefficients only depend on the difference $j_2 - j_1$. We point out that these scale invariant patterns are signatures of scale invariant Gaussian processes (Bruna et al. 2015) and have already been observed for fractional Brownian motions processes in Allys et al. (2019).

$\hat{S}_2^{\text{iso},2}$ coefficients also show two distinct trends: the coefficients quickly tend to zero when $j_2 - j_1$ increases for $R[\tilde{Q} + i\tilde{U}]$ while coefficients tend to strictly positive values for the largest $j_2 - j_1$ values for non randomized data sets $\tilde{Q} + i\tilde{U}$. This result is related to the filamentary structure of the non randomized maps, because a filamentary structure involves a modulation of the RWST coefficients as a function of the angle difference $\theta_2 - \theta_1$. We see the same trend in Fig. 12 for the $\hat{S}_2^{\text{iso},2}$ coefficients of the other non randomized data sets which all present filamentary structures. For all data sets, and fixed $j_1$ values, $\hat{S}_2^{\text{iso},2}$ coefficients decrease as a function of $j_2$. This property seems to be general. We also expect a signal due to the imprint of the wavelets in the $\hat{S}_2^{\text{iso},2}$ coefficients that would increase their values when $j_2$ is close to $j_1 + 1$.

We see in Fig. 11 that the differences between the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ data set and its randomized counterpart decrease for the highest $j_1$ values. Here again, this shows that the non-Gaussianity of $\tilde{Q} + i\tilde{U}$ maps decreases at large scales. This is consistent with what we have already pointed out in Sect. 4.1 for $\hat{S}_1^{\text{iso}}$ coefficients and we interpret this trend similarly.

Second order anisotropic coefficients $\hat{S}_2^{\text{aniso},1}$, $\hat{S}_2^{\text{aniso},2}$, and $\theta^{\text{ref},2}$ essentially show consistent results with first order anisotropic coefficients both in terms of amplitude and direction of anisotropy. They have generally smaller values compared to those of isotropic coefficients $\hat{S}_2^{\text{iso},1}$ and $\hat{S}_2^{\text{iso},2}$. We do not show these coefficients here as we do not discuss them any further in this work.

## 5. Random syntheses of $\tilde{Q} + i\tilde{U}$ polarization maps

The RWST coefficients are not only statistical descriptors of the data but can also be used as a set of constraints for generating multiple random realizations of polarization maps. By
<figure><img src="image_13.png" /><figcaption>Fig. 12. $\hat{S}_{2}^{\text{iso},1}(j_1, j_2)$ (left column) and $\hat{S}_{2}^{\text{iso},2}(j_1, j_2)$ (right column) RWST coefficients for $\tilde{Q} + i\tilde{U}$, $\exp(2i\psi)$ and $p$ data sets in the $\perp$ case. Each curve corresponds to a fixed $j_1$ value with $j_2$ varying from $j_1 + 1$ to $J - 1 = 6$. For $j_1 = J - 2$, the curve is reduced to a single dot on the figure.</figcaption></figure>

comparing the original and synthetic maps one can assess the exhaustiveness of the statistical description. A visual agreement would give more confidence in the relevance of the statistical information captured by a RWST description. However, this qualitative assessment has its limits as it is unclear what kind of statistics human eyes are sensitive to (Julesz 1981). A quantitative comparison using summary statistics other than the RWST is thus also needed.

The synthesis is also a means of simulating noise-free maps artificially and realistically augmenting demanding MHD simulations, or producing realistic foregrounds for component separation methods for CMB data analysis. This approach differs from previous work where dust polarization maps are computed from a phenomenological model and ancillary observations (Planck Collaboration Int. XLII 2016; Ghosh et al. 2017; Vansyngel et al. 2017; Levrier et al. 2018; Adak et al. 2020; Clark & Hensley 2019) or directly from numerical simulations (Kim et al. 2019). It follows what has been done for dust total intensity by Allys et al. (2019) using the RWST, and Aylor et al. (2019) using a Generative Adversarial Network.

In this section, we present random synthetic $\tilde{Q} + i\tilde{U}$ maps generated from RWST coefficients derived from the analysis of our MHD simulation and complemented by additional constraints on the large-scale components of the maps and on the one-point probability distribution function. We then compare synthetic and original maps using one-point and two-point statistics.

## 5.1. Generation of synthetic polarization maps from a RWST description

If synthetic total intensity maps generated from a RWST description have already been produced in Allys et al. (2019), in this paper we extend and improve this procedure for complex polarization maps $\tilde{Q} + i\tilde{U}$. We give in this subsection some technical

details on the implementation but details on the mathematical formalism can be found in Bruna & Mallat (2019).

The generation of synthetic maps is an iterative process that consists in the minimization of a loss function $\mathcal{L}$ in pixel space. Hence this optimization problem is defined in practice in a space of dimension $512^2$. Starting from a realization of a complex Gaussian white noise map $X_0 = \tilde{Q}_0 + i\tilde{U}_0$, successive maps $\{X_k\}$ are built with a quasi-Newton method. We use a L-BFGS-B algorithm, although we do not impose any boundary to the values of the pixels. Technical details on quasi-Newton methods and L-BFGS-B algorithm can be found in Fletcher (1987) and Byrd et al. (1995).

In the following, one specific map $X^t = \tilde{Q}^t + i\tilde{U}^t$ of the data set on which the RWST description is built serves as a reference for comparison with the synthetic maps. Because our WST and RWST analyses does not extend to the largest scales, we have no statistical information on these. Thus, we replace the largest scales of $X_0$ with those of $X^t$ to address this gap in a deterministic way<fnref n="14" />. Formally, noting $\mathcal{F}[X_0]$ the Fourier transform of $X_0$ initially drawn as a realization of a Gaussian white noise we set:

$$ \mathcal{F}[X_0](k) = \mathcal{F}[X^t](k) \text{ for } |k| < k_{\min}, \qquad (25) $$

where $k_{\min}$ is the wave number corresponding to the largest scale probed by the WST. For $J=7$, using Table 1 we have $k_{\min} \approx 2\pi/171 \approx 0.037 \text{ pixel}^{-1}$. In practice, 25 Fourier modes out of $512^2$ are set by this procedure. Even if in our case this value of $k_{\min}$ is related to an ad hoc modeling decision (J value), we point out that statistical approaches are not always relevant at all scales for the analysis of the diffuse ISM. For example, in a turbulent medium we know that injection scales generally correspond to specific events (e.g., a supernova) for which a statistical

<footnote n="14">¹⁴ We note that such an approach cannot supply the statistical information concerning the couplings between the largest scales with smaller ones.</footnote>
<figure><img src="image_14.png" /><figcaption>Fig. 13. Synthesis of a $\tilde{Q}_{\parallel} + i\tilde{U}_{\parallel}$ map (right column) built from its corresponding RWST description with additional constraints on large-scale components and on a few one-point moments of $\tilde{Q}$ and $\tilde{U}$ maps, to be consistent with those of a reference map (left column). The reference maps shown here are the same as in Fig. 3.</figcaption></figure>

description has little meaning, while a statistical description is relevant to describe the inertial range. Moreover, the WST of an image does not adequately characterize its one-point statistics, so in the following we also constrain these for the synthetic maps using $X^r$ one-point statistics as a reference.

We define the loss function $\mathcal{L}$ of a complex image $X = \tilde{Q} + i\tilde{U}$ as follows:

$$ \mathcal{L}[X] = \mathcal{L}_{\text{WST}}[X] + \mu (\mathcal{L}_{\text{one-point}}[\tilde{Q}] + \mathcal{L}_{\text{one-point}}[\tilde{U}]), \quad (26) $$

with $\mathcal{L}_{\text{WST}}[X]$ the loss function constraining the WST coefficients of $X$, $\mathcal{L}_{\text{one-point}}[\tilde{Q}]$ (respectively $\mathcal{L}_{\text{one-point}}[\tilde{U}]$) that constraining a few one-point moments of $\tilde{Q}$ (respectively $\tilde{U}$), and $\mu$ a weighting coefficient balancing the importance of one-point moments constraints with that of WST coefficients constraints. More specifically we set:

$$ \mathcal{L}_{\text{WST}}[X] = \frac{1}{N} \left( \begin{aligned} & (S_0[X] - S_0^t)^2 + \sum_{j_1, \theta_1} (S_1(j_1, \theta_1)[X] - S_1^t(j_1, \theta_1))^2 \\ & + \sum_{j_1, j_2, \theta_1, \theta_2} (S_2(j_1, \theta_2, j_2, \theta_2)[X] - S_2^t(j_1, \theta_1, j_2, \theta_2))^2 \end{aligned} \right), \quad (27) $$

where the “t” superscript refers to the target WST coefficients that the synthetic map should have, and $N$ is the total number of WST coefficients. In our case ($J=7$ and $\Theta=8$) we recall that we have $N=2801$. These target WST coefficients are computed from the RWST coefficients derived from a given data set of $\tilde{Q} + i\tilde{U}$ maps (of which $X^r$ is part) using the RWST model defined in Eqs. (23) and (24). We see on this loss function that none of the WST coefficient is privileged and that we do not weigh the differences of WST coefficients by any uncertainty on

the target coefficients. This is certainly something that can be improved in future work on WST syntheses. We finally define:

$$ \mathcal{L}_{\text{one-point}}[\tilde{Q}] = \frac{1}{3} \left( \left( \frac{M_2[\tilde{Q}]}{M_2[\tilde{Q}^r]} - 1 \right)^2 + \left( M_3[\tilde{Q}] - M_3[\tilde{Q}^r] \right)^2 + \left( M_4[\tilde{Q}] - M_4[\tilde{Q}^r] \right)^2 \right), \quad (28) $$

where $M_2[\tilde{Q}] = \langle (\tilde{Q} - \langle \tilde{Q} \rangle)^2 \rangle$ is an estimation of the variance of the distribution of pixels values in the $\tilde{Q}$ map, $M_3[\tilde{Q}] = \langle (\tilde{Q} - \langle \tilde{Q} \rangle)^3 \rangle / M_2[\tilde{Q}]^{3/2}$ is an estimation of its skewness, and $M_4[\tilde{Q}] = \langle (\tilde{Q} - \langle \tilde{Q} \rangle)^4 \rangle / M_2[\tilde{Q}]^2$ is an estimation of its kurtosis.

We generate synthetic maps using the $\tilde{Q}_{\parallel} + i\tilde{U}_{\parallel}$ RWST description as a target but we also provide equivalent results for the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ RWST description in Appendix C. The optimization algorithm stops when the values of the loss function stagnates between two consecutive iterations. In practice we choose $\mu = 5 \times 10^{-8}$ and the algorithm stops when $|\mathcal{L}[X_k] - \mathcal{L}[X_{k+1}]| \lesssim 10^{-12}$. These numerical values correspond to a reasonable trade-off between the quality of the syntheses and the execution time of the algorithm. They lead to constrained statistical coefficients with an accuracy better than a few percent on average. Once the optimization is done, we finally filter the modes of the resulting maps at higher wave numbers than the one dimensional Nyquist wave number ($k_N = \pi \text{ pixels}^{-1}$) to avoid unwanted numerical artifacts. Indeed, the WST relies on a bank of Morlet wavelets that does not properly cover this range of frequencies, which results in a total loss function $\mathcal{L}$ that does not properly constrain these modes. Figure 13 shows a synthetic $\tilde{Q}_{\parallel} + i\tilde{U}_{\parallel}$ map (right column) next to its reference map (left column). The overall appearance of the synthetic maps is satisfactory. The largest scales are roughly consistent with
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Power spectra of the synthetic $\tilde{Q}_{\parallel} + i\tilde{U}_{\parallel}$ map shown in Fig. 13, compared to those of the reference map, for $\tilde{Q}, \tilde{U}$, and $|\tilde{Q} + i\tilde{U}|$ (left, middle and right respectively). The vertical dashed lines mark the wavelet central wave numbers corresponding to the scale indices $j=0, \dots, J-1$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** One-point distribution functions of the synthetic $\tilde{Q}_{\parallel} + i\tilde{U}_{\parallel}$ map shown in Fig. 13, compared to those of the reference map, for $\tilde{Q}, \tilde{U}$, and $|\tilde{Q} + i\tilde{U}|$ (left, middle and right respectively).</figcaption></figure>

those of the reference maps<fnref n="15" />, but we also see at intermediate and small scales, which are the truly synthetic scales, consistent filamentary patterns and dynamic ranges.

## 5.2. One-point and two-point statistics of synthetic maps

By construction, the synthetic maps we built have the same WST statistics as the one prescribed by the RWST description of the $\tilde{Q}_{\parallel} + i\tilde{U}_{\parallel}$ data set, similar large scales and some identical one-point moments compared to the reference maps. We may wonder if elementary one-point and two-point statistics are fully consistent with the ones of the reference maps.

Figure 14 shows the power spectra of $\tilde{Q}, \tilde{U}$, and $\tilde{Q} + i\tilde{U}$ up to the one dimensional Nyquist wave number $k_N$ for both the reference maps and the syntheses shown in Fig. 13. The power spectra are computed by binning the squared amplitudes of Fourier modes with respect to the modulus of the corresponding wave number $k$. We use a regular binning in $k$ and the estimations of the power spectra are computed as the means for each bin. We also represent standard deviations of the mean per bin. We see that the power spectra as well as their standard deviations are in good agreement for the three variables $\tilde{Q}, \tilde{U}$, and $\tilde{Q} + i\tilde{U}$ for all scales except the smallest ones. These discrepancies at small scales take the form of a lack followed by an excess of power in the syntheses for wave numbers approaching $k_N$. We interpret this as the result of poorly constrained modes close to

$k_N$ in the optimization process. It is not surprising to reproduce the power spectrum of $\tilde{Q} + i\tilde{U}$ maps as the $\tilde{S}_1$ and $\tilde{S}_2$ coefficients constrain the power spectrum of $\tilde{Q} + i\tilde{U}$ (see discussion in Sect. 2.3). However we point out that we correctly reproduce the power spectra of $\tilde{Q}$ and $\tilde{U}$ taken separately. Still we note that we did not investigate cases where $\tilde{Q}$ and $\tilde{U}$ have very different power spectra.

In Fig. 15, we compare the one-point distribution functions of the reference maps and the synthetic maps for $\tilde{Q}, \tilde{U}$, and $|\tilde{Q} + i\tilde{U}|$. These are also in fairly good agreement for all the variables. We notably successfully reproduce the tails of these distributions, this must be due to the combined constraints on WST coefficients and one-point moments. One could surely enhance the agreement between the reference and synthetic maps by taking into account higher order moments in the $L_{\text{one-point}}$ loss functions.

We point out that such syntheses were generated using RWST descriptions that comprise 127 coefficients (see Sect. 3.4). Adding to that the largest scales that were set in a deterministic way as well as the constraints on the one-point moments of $\tilde{Q}$ and $\tilde{U}$ maps, we end up with a total of 158 coefficients to generate 512×512 complex $\tilde{Q} + i\tilde{U}$ maps that are in good visual agreement with the maps of the original data set and successfully reproduce their one-point and two-point statistics. Excluding the largest scales, these sets of statistical coefficients not only describe the statistical properties of these maps, but they also define a statistical model of the mechanism that generated these data. In a mathematical wording, these coefficients model the probability distribution of the random field of which these maps are realizations. Although these statistical syntheses will

<footnote n="15">¹⁵ Since the WST constraints leave the largest scales of $X_0$ unchanged, this shows that the impact of the one-point moments constraints on these scales is minor.</footnote>
definitely benefit from more dedicated work, they already are a
promising avenue to address the construction of the statistical
model of polarization maps that we seek.

**6. Conclusions and perspectives**

In this paper, we extended the WST analysis to maps of polarized thermal emission from interstellar dust, using 512 × 512 pixels Stokes *I*, *Q*, and *U* maps built from a numerical simulation of MHD turbulence designed to reproduce typical properties of the diffuse ISM. To alleviate the fact that Stokes *Q* and *U* rely on the definition of an arbitrary reference frame, and to remove the zeroth-order impact of the matter distribution on their properties and thus focus on the statistics of the magnetic field, the WST was applied to complex Stokes maps $\tilde{Q} + i\tilde{U}$ that are normalized by *I* + *P*. To study the contributions of the polarization fraction *p* and of the polarization angle $\psi$ to the statistical properties of these complex Stokes maps, we also applied the WST to the corresponding maps of *p* and $\exp(2i\psi)$. We finally analyzed “Gaussianized” complex Stokes maps obtained after phase randomization.

The WST gives a low-variance statistical description of these complex and real maps through typically a few thousand coefficients indexed in terms of orientations and scales. These coefficients capture the power spectra of the maps and characterize couplings between oriented scales. WST coefficients for maps of $\tilde{Q} + i\tilde{U}$, *p*, and $\exp(2i\psi)$ present a striking regularity when taken as functions of the sole angular variables. This is very much in line with what we observed in Allys et al. (2019) for column density and total intensity maps, and in fact the same functional form introduced in that paper can be used to fit the angular dependencies of the WST coefficients of polarization maps studied here, thus extending the RWST model introduced in Allys et al. (2019). The RWST yields a statistical description of polarization maps that quantifies their multiscale properties in terms of isotropic and anisotropic contributions, all the while requiring more than one order of magnitude fewer coefficients than the WST.

In the rest of this section, we summarize the main results of
our work, then highlight a few perspectives. The RWST analysis
allowed us to identify statistical characteristics that exhibit the
dependence of the map structure on the orientation of the mean
magnetic field and quantify the non-Gaussianity of data.

- The overall level of first order coefficients depends on the orientation of the mean magnetic field with respect to the line of sight. For $\tilde{Q} + i\tilde{U}$ maps, $\hat{S}_1^{\text{iso}} + \log_2(\hat{S}_0)$ coefficients are larger when the mean magnetic field is in the plane of the sky, while for $\exp(2i\psi)$ maps $\hat{S}_1^{\text{iso}}$ is larger when the mean magnetic field is along the line of sight.

- $\hat{S}_1^{\text{aniso}}$ coefficients quantify the statistical anisotropy of the maps. When the mean magnetic field is parallel to the line of sight $\hat{S}_1^{\text{aniso}}$ coefficients are negligible, while when the mean magnetic field is in the plane of the sky they allow us to identify the direction of anisotropy at each scale. For the MHD simulation we analyzed, this direction is orthogonal to the direction of the mean magnetic field for both the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ and $\exp(2i\psi)_{\perp}$ maps.

- Second order RWST coefficients clearly exhibit the non-Gaussianity of $\tilde{Q} + i\tilde{U}$ maps (although this is also visible to a lesser extent in first order coefficients). While the randomized $R[\tilde{Q} + i\tilde{U}]$ data sets present characteristic properties of scale invariant Gaussian fields (invariance of $\hat{S}_2^{\text{iso},1}$ as a function of the scale difference $j_2 - j_1$ and a quick decrease

to zero of $\hat{S}_2^{\text{iso},2}$ as this scale difference increases), the $\hat{S}_2^{\text{iso},1}$ and $\hat{S}_2^{\text{iso},2}$ coefficients for the corresponding $\tilde{Q} + i\tilde{U}$ data sets show clearly different patterns. In particular, the strictly positive values of $\hat{S}_2^{\text{iso},2}$ at large $j_2 - j_1$ are interpreted as signatures of the filamentary structure of the maps.

We have used the RWST approach to synthesize $\tilde{Q} + i\tilde{U}$ maps. Combining the RWST coefficients with additional constraints, we obtained synthetic maps that statistically match the original maps. The agreement demonstrates the comprehensive-ness of the statistical description provided by the RWST. The additional constraints include large scale modes that cannot be described statistically with the limited amount of samples we have worked with, as well as statistical constraints on the one-point distribution function.

In this paper, to establish the methodology, we have worked with noise-free maps computed from numerical simulations of MHD turbulence. A future extension would be to apply it to observations of polarized thermal emission from dust. To do this, we need to learn how to handle data noise. Indeed, while signal-to-noise ratios in total intensity for both *Herschel* and *Planck* maps are quite high, this is not the case for currently available polarization data. Studying how data noise affects the WST coefficients would demand to repeat the analysis of MHD simulations with noise added to the dust signal. Once this difficulty is overcome, we may use the RWST to define a metric to compare observations with simulations and phenomenological models. This will be a stepping stone towards a more refined physical interpretation of the RWST coefficients. A main motivation would be to use the RWST to characterize magnetized interstellar turbulence.

Throughout this work we chose to work with Stokes *I*, *Q*, and *U* maps to analyze the polarization of dust thermal emission as astronomers do for Galactic astrophysics. In the framework of CMB data analysis, polarization is more usually characterized through E and B modes. Thus, it would be interesting to apply the RWST to E and B maps. This will lead to a physically motivated statistical model that would include observational constraints such as the E/B power asymmetry and the correlation between E-modes of the polarization and total intensity (Planck Collaboration XI 2020). Once this is achieved, the RWST may be used to synthesize dust polarization maps that would help to assess and optimize component separation methods for CMB data analysis.

*Acknowledgements.* We thank the anonymous referee for a very helpful report. We gratefully acknowledge V. Jeli? for pointing out to us the complex-valued description of polarization maps. We also acknowledge fruitful discussions with T. Marchand, J.-F. Cardoso, A. Frolov, S. Zhang, and S. Mallat that enriched this work in various ways. This research was supported by the Agence Nationale de la Recherche (project BxB: ANR-17-CE31-0022). F.L. and F.B. acknowledge support from the European Research Council under the European Union's Horizon 2020 Research & Innovation Framework Programme / ERC grant agreement ERC-2016-ADG-742719. B.R.-S.B. acknowledges support from the Centre National d'Etudes Spatiales (CNES).

**References**

Abazajian, K., Addison, G., Adshead, P., et al. 2019, ArXiv e-prints,
[arXiv:1907.04473]

Adak, D., Ghosh, T., Boulanger, F., et al. 2020, A&A 640, A100

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

Allys, E., Levrier, F., Zhang, S., et al. 2019, A&A, 629, A115

Andersson, B. G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501

Aylor, K., Haq, M., Knox, L., Hezaveh, Y., & Perreault-Levasseur, L. 2019, ArXiv e-prints, [arXiv:1909.06467]
BICEP2 Collaboration 2018, Phys. Rev. Lett., 121, 221301

BICEP2/Keck Array and Planck Collaborations 2015, Phys. Rev. Lett., 114, 101301

Bruna, J., & Mallat, S. 2013, IEEE Trans. Pattern Anal. Mach. Intell., 35, 1872

Bruna, J., & Mallat, S. 2019, Math. Stat. Learn., 1, 257

Bruna, J., Mallat, S., Bacry, E., & Muzy, J.-F. 2015, Ann. Stat., 43, 323

Byrd, R., Lu, P., Nocedal, J., & Zhu, C. 1995, SIAM J. Sci. Comput., 16, 1190

Chuss, D. T., Andersson, B. G., Bally, J., et al. 2019, ApJ, 872, 187

Clark, S. E., & Hensley, B. S. 2019, ApJ, 887, 136

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton: Princeton University Press)

Draine, B. T., & Fraisse, A. A. 2009, ApJ, 696, 1

Farge, M. 1992, Ann. Rev. Fluid Mech., 24, 395

Federrath, C., Roman-Duval, J., Klessen, R., Schmidt, W., & Mac Low, M. M. 2009, A&A, 512, A81

Fissel, L. M., Ade, P. A. R., Angilè, F. E., et al. 2016, ApJ, 824, 134

Fletcher, R. 1987, Practical Methods of Optimization, 2nd edn. (Chichester: John Wiley & Sons)

Fraisse, A. A., Ade, P. A. R., Amiri, M., et al. 2013, J. Cosmol. Astropart. Phys., 2013, 047

Frisch, U. 1995, Turbulence, the Legacy of A.N. Kolmogorov (Cambridge: Cambridge University Press)

Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371

Galerne, B., Gousseau, Y., & Morel, J.-M. 2011, IEEE Trans. Image Process., 20, 257

Ghosh, T., Boulanger, F., Martin, P. G., et al. 2017, A&A, 601, A71

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Guillet, V., Fanciullo, L., Verstraete, L., et al. 2018, A&A, 610, A16

Guth, A. H. 1981, Phys. Rev. D, 23, 347

Haverkorn, M., Katgert, P., & de Bruyn, A. G. 2004, A&A, 427, 549

Hennebelle, P., & Falgarone, E. 2012, A&ARv, 20, 55

Hennebelle, P., & Inutsuka, S. I. 2019, Front. Astron. Space Sci., 6

Hull, C. L. H., Mocz, P., Burkhart, B., et al. 2017, ApJ, 842, L9

Iffrig, O., & Hennebelle, P. 2017, A&A, 604, A70

Ishino, H., Akiba, Y., Arnold, K., et al. 2016, SPIE Conf. Ser., 9904, 99040X

Julesz, B. 1981, Nature, 290, 91

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

Kim, C.-G., Choi, S. K., & Flauger, R. 2019, ApJ, 880, 106

Krizhevsky, A., Sutskever, I., & Hinton, G. E. 2012, Adv. Neural Inf. Process. Syst. 25, 1097

Kroese, D. P., & Botev, Z. I. 2015, Spatial Process Simulation, ed. V. Schmidt (Cham: Springer International Publishing), 369

Landi Degl'innocenti, E., & Landolfi, M. 2004, Polarization in Spectral Lines (Netherlands: Springer Netherlands)

LeCun, Y., Kavukcuoglu, K., & Farabet, C. 2010, in Proceedings of 2010 IEEE International Symposium on Circuits and Systems, 253

Levrier, F., Neveu, J., Falgarone, E., et al. 2018, A&A, 614, A124

Linde, A. D. 1982, Phys. Lett. B, 108, 389

Mallat, S. 2012, Commun. Pure Appl. Math., 65, 1331

Mangilli, A., Aumont, J., Bernard, J. P., et al. 2019, A&A, 630, A74

Miville-Deschênes, M. A., Martin, P. G., Abergel, A., et al. 2010, A&A, 518, L104

Naess, S., Hasselfield, M., McMahon, J., et al. 2014, J. Cosmol. Astropart. Phys., 2014, 007

Oppenheim, A. V., & Lim, J. S. 1981, Proc. IEEE, 69, 529

Planck Collaboration XI. 2020, A&A, 641, A11

Planck Collaboration XII. 2020, A&A, 641, A12

Planck Collaboration Int. XIX. 2015, A&A, 576, A104

Planck Collaboration Int. XX. 2015, A&A, 576, A105

Planck Collaboration Int. XXX. 2016, A&A, 586, A133

Planck Collaboration Int. XLII. 2016, A&A, 596, A103

Reissl, S., Guillet, V., Brauer, R., et al. 2020, A&A 640, A118

Ritacco, A., Adam, R., Ade, P., et al. 2020, Eur. Phys. J. Web Conf., 228, 00022

Schmidt, W., Hillebrandt, W., & Niemeyer, J. C. 2006, Comput. Fluids, 35, 353

Siebenmorgen, R., Voshchinnikov, N. V., & Bagnulo, S. 2014, A&A, 561, A1

Sutin, B. M., Alvarez, M., Battaglia, N., et al. 2018, SPIE Conf. Ser., 10698, 106984F

Teyssier, R. 2002, A&A, 385, 337

Vansyngel, F., Boulanger, F., Ghosh, T., et al. 2017, A&A, 603, A62

Wandelt, B. D. 2013, Astrostatistical Challenges for the New Astronomy, ed. J. M. Hilbe (New York: Springer New York), 87
# Appendix A: The RWST model in terms of Fourier series expansions

The RWST model can be rephrased in terms of truncated Fourier series. Let us consider the $log_2(\tilde{S}_1)$ coefficients at a given scale $j_1$ and write $f_{j_1}(\theta_1)$ the corresponding angular model<fnref n="16" />. If one assumes that there are no more than one privileged direction $\theta^{\text{ref},1}(j_1)$ in the maps we are dealing with, we can write generally $f_{j_1}(\theta_1)$ as a Fourier series expansion using the natural $2\pi$-periodicity of this function:

$$f_{j_1}(\theta_1) = a_0(j_1) + \sum_{k=1}^{+\infty} [a_k(j_1) \cos(k[\theta_1 - \theta^{\text{ref},1}(j_1)]) \\ + b_k(j_1) \sin(k[\theta_1 - \theta^{\text{ref},1}(j_1)])]. \quad (A.1)$$

Assuming a mirror symmetry with respect to the potential reference direction, we expect {$b_k$} coefficients to vanish (which is the case in practice for our data). Adding to that the $\pi$-periodicity identified in Sect. 3.3 the angular model reduces to:

$$f_{j_1}(\theta_1) = a_0(j_1) + \sum_{k=2}^{+\infty} a_k(j_1) \cos(k[\theta_1 - \theta^{\text{ref},1}(j_1)]). \quad (A.2)$$

Finally, the smoothness of the patterns presented in Sect. 3.3 implies a fast decrease of the amplitudes of the harmonics. Truncating the expansion after the second term and writing $a_0(j_1) = \hat{S}_1^{\text{iso}}(j_1)$ and $a_2(j_1) = \hat{S}_1^{\text{aniso}}(j_1)$ we end up with the RWST model of Eq. (23).

# Appendix B: Phase randomization

We give some technical details concerning the phase randomization process of a discretized scalar field.

Let us consider a 2D discretized real scalar field $f: \Omega \to \mathbb{R}$ with $\Omega = [0, N-1]^2$. Its discrete Fourier transform reads for a given mode $\mathbf{k} = (k_x, k_y) \in \Omega$:

$$\hat{f}(\mathbf{k}) = \sum_{\mathbf{r} \in \Omega} f(\mathbf{r}) \exp \left( -2\pi i \frac{\mathbf{r} \cdot \mathbf{k}}{N} \right), \quad (B.1)$$

where $\mathbf{r} \cdot \mathbf{k} = xk_x + yk_y$. We also recall the inverse discrete Fourier transform relation for the chosen convention:

$$f(\mathbf{r}) = \frac{1}{N^2} \sum_{\mathbf{k} \in \Omega} \hat{f}(\mathbf{k}) \exp \left( 2\pi i \frac{\mathbf{r} \cdot \mathbf{k}}{N} \right). \quad (B.2)$$

The phase randomization of $f$ consists in defining a new field $g$ in Fourier space such that:

$$\hat{g}(\mathbf{k}) = |\hat{f}(\mathbf{k})| e^{i\phi(\mathbf{k})}, \quad (B.3)$$

for every mode $\mathbf{k} \in \Omega$ and where $\phi$ is a realization of a uniform random phase such as defined in Galerne et al. (2011), that is, defining $\Omega_0 = \{(0, N/2), (N/2, 0), (N/2, N/2)\}$ if $N$ is even, and $\Omega_0 = \emptyset$ otherwise, $\phi$ verifies:

1. $\forall \mathbf{k} \in \Omega \setminus \Omega_0, \phi(-\mathbf{k}) = -\phi(\mathbf{k})$ (we extend the domain of $\phi$ to $\mathbb{Z}^2$ using periodic boundary conditions).

2. $\forall \mathbf{k} \in \Omega, \phi(\mathbf{k})$ is drawn uniformly and independently in $[0, 2\pi)$ (the independence holds to the extent of the first relation).

<figure><img src="image_17.png" /><figcaption>Fig. B.1. Example of a phase randomized $R[\tilde{Q}_{\perp}]$ map (right) next to its corresponding original $\tilde{Q}_{\perp}$ map (left).</figcaption></figure>

3. $\forall \mathbf{k} \in \Omega_0, \phi(\mathbf{k})$ is drawn uniformly and independently in the set $\{0, \pi\}$.

These relations ensure that $g$ will be a real field, and that the mean of $g$ will be equal to the absolute value of the mean of $f$ (since condition 1. gives $\phi(0) = 0$). Such a process obviously defines a field $g$ which preserves the power spectrum of $f$. For reference, we show an example of a phase randomized $R[\tilde{Q}_{\perp}]$ map in Fig. B.1.

# Appendix C: Additional results

## C.1. Additional terms in the RWST model

Following the appendix C of Allys et al. (2019), we enhance the RWST model defined in Eqs. (23) and (24) by adding so-called lattice terms related to pixellization effects at small scales for first order coefficients and a supplementary harmonic of the angular modulation of the second order WST coefficients. This enhanced RWST model of the angular dependency of the WST coefficients becomes, for $\tilde{S}_1$ coefficients:

$$\log_2[\tilde{S}_1(j_1, \theta_1)] = \hat{S}_1^{\text{iso}}(j_1) \\ + \hat{S}_1^{\text{aniso}}(j_1) \cos(2[\theta_1 - \theta^{\text{ref},1}(j_1)]) \\ + \hat{S}_1^{\text{lat},1}(j_1) \cos(4\theta_1) + \hat{S}_1^{\text{lat},2}(j_1) \cos(8\theta_1), \quad (C.1)$$

where the additional lattice terms $\hat{S}_1^{\text{lat},1}$ and $\hat{S}_1^{\text{lat},2}$ quantify angular modulations that are respectively $\pi/2$ and $\pi/4$-periodic and aligned with the main directions of the lattice. Similarly the enhanced RWST model of $\tilde{S}_2$ coefficients is:

$$\log_2[\tilde{S}_2(j_1, \theta_1, j_2, \theta_2)] = \hat{S}_2^{\text{iso},1}(j_1, j_2) \\ + \hat{S}_2^{\text{iso},2}(j_1, j_2) \cos(2[\theta_1 - \theta_2]) \\ + \hat{S}_2^{\text{iso},3}(j_1, j_2) \cos(4[\theta_1 - \theta_2]) \\ + \hat{S}_2^{\text{aniso},1}(j_1, j_2) \cos(2[\theta_1 - \theta^{\text{ref},2}(j_1, j_2)]) \\ + \hat{S}_2^{\text{aniso},2}(j_1, j_2) \cos(2[\theta_2 - \theta^{\text{ref},2}(j_1, j_2)]), \quad (C.2)$$

where the additional term $\hat{S}_2^{\text{iso},3}$ measures the amplitude of an additional harmonic of the $\theta_1 - \theta_2$ angular modulation that is $\pi/2$-periodic.

These additional terms do not affect the values of the RWST coefficients discussed in the main body of the paper and significantly improve reduced chi square $\chi_r^{2,S_1}$ and $\chi_r^{2,S_2}$ at small scales as shown in Fig. C.1. These chi square functions are to be compared to those of Fig. 7. Figures C.2 and C.3 show greatly

<footnote n="16">¹⁶ The same kind of reasoning holds for $log_2(\tilde{S}_2)$ coefficients.</footnote>
<figure><img src="image_18.png" /></figure>

**Fig. C.1.** Reduced chi square $\chi_r^{2,S_1}(j_1)$ (top) and $\chi_r^{2,S_2}(j_1, j_2)$ (bottom) associated with the RWST fits of the WST coefficients that take into account lattice terms (see Eqs. (C.1) and (C.2)) for the $\tilde{Q}_\perp + i\tilde{U}_\perp$ data set. Each curve in the $\chi_r^{2,S_2}(j_1, j_2)$ plot corresponds to a fixed $j_1$ value while $j_2$ ranges from $j_1 + 1$ to $J - 1 = 6$. For $j_1 = J - 2$, the curve is reduced to a single dot on the figure. We use logarithmic scales for better visibility.

improved normalized residuals compared to the previous ones shown in Figs. 4 and 5. In particular at $j_1=0$ we no longer observe the strong angular pattern seen in Fig. 4.

For completeness, we show in Fig. C.4 the additional RWST terms given for the example of the $\tilde{Q}_\perp + i\tilde{U}_\perp$ data set.

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

**Fig. C.2.** Same as Fig. 4 but for the RWST fit corresponding to the model of Eq. (C.1) including lattice terms.

### C.2. Syntheses for $\tilde{Q}_\perp + i\tilde{U}_\perp$

We show in Figs. C.5–C.7 the maps, power spectra, and distribution functions of synthetic maps in the $\perp$ case that are analogous to those of Figs. 13–15.
<figure><img src="image_20.png" /><figcaption>Fig. C.3. Same as Fig. 5 but for the RWST fit corresponding to the model of Eq. (C.1) including lattice terms.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. C.4. $\hat{S}_1^{\text{lat},1}(j_1)$ (left column), $\hat{S}_1^{\text{lat},2}(j_1)$ (middle column) and $\hat{S}_2^{\text{iso},3}(j_1, j_2)$ (right column) RWST coefficients for the $\tilde{Q}_{\perp} + i\tilde{U}_{\perp}$ data set. In the $\hat{S}_2^{\text{iso},3}$ plot, each curve corresponds to a fixed $j_1$ value while $j_2$ ranges from $j_1 + 1$ to $J - 1 = 6$. For $j_1 = J - 2$, the curve is reduced to a single dot on the figure.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. C.5.** Same as Fig. 13 but for the ⊥ case. The reference maps shown here are the same maps as in Fig. 2.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. C.6.** Same as Fig. 14 but for the ⊥ case.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. C.7.** Same as Fig. 15 but for the ⊥ case.</figcaption></figure>