Statistical description of dust polarized emission from the diffuse interstellar medium -- A WST/RWST approach

The statistical characterization of the diffuse magnetized ISM and Galactic foregrounds to the 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 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. We apply the WST to noise-free maps of dust polarized thermal emission computed from a numerical simulation of MHD turbulence. We analyze normalized complex Stokes maps, and maps of 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 same Reduced WST (RWST) angular model introduced by Allys+2019 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 allows 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 use 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 Python package to perform WST/RWST analyses at: https://github.com/bregaldo/pywst.


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). Non-spherical 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 provided us with the first all-sky survey of polarized dust emission, opening the path to statistical studies (Planck Collaboration XII 2019). This broad view is being complemented by observations at higher angular resolution, carried out by the balloon-borne experiments BLAST-POL (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 endeavour, 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 GW on the one hand, imprinted in the curl-like "Bmode" 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 2019). The current constraints on the cosmological B-mode tensor-to-scalar ratio r < 0.06 (95% confidence level, BICEP2 Collaboration 2018) Article number, page 1 of 20 arXiv:2007.08242v1 [astro-ph.CO] 16 Jul 2020 A&A proofs: manuscript no. Main 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 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, 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/RWST analysis to maps of linearly polarized emission, which is represented by Stokes Q and U maps 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 noisefree 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/RWST analyses of two-dimensional data called PyWST at: https://github.com/bregaldo/pywst.

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.

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 ν for the Stokes parameters I, Q, and U: where S ν is the source function which is assumed here to be that of a black body so that S ν = B ν (T d ) with T d the dust temperature, τ ν is the optical depth, p 0 is an intrinsic polarization fraction parameter, γ is the angle that the local magnetic field makes with the plane of the sky, and φ 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 2 (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τ ν = σ ν n H dz where σ ν is the dust opacity (assumed uniform) and n H is the total gas density. In this case, Eqs. (1)-(3) become: U = p 0 σ ν S ν n H sin (2φ) cos 2 γdz. From these formulae, 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 γ and φ. 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 ϕ transforms Q and U into rotated quantities Q and U as follows 3 : These relations can be written in a more compact form introducing complex quantities: 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 will 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 non-linearity 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.

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 variableQ + iŨ, polarization fraction p and complex polarization angle exp(2iψ).
We first define the normalized Stokes complex variablẽ Q + iŨ as follows: 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), (5), and (6): which is independent of the density field n H . The polarization fraction p and the polarization angle ψ follow the usual definitions: and where atan2(b, a) is the function that returns the angle restricted to (−π, π] between the positive x-axis and the ray to the point (a, b) in the Euclidean plane. Because of the non-linearity of the atan2 function, we expect the statistics of ψ 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ψ). This variable also avoids dealing with discontinuities of ψ maps that appear where the polarization angle approaches ±π/2. Moreover, the statistical properties of exp(2iψ) are easier to compare to those ofQ + iŨ sinceQ + iŨ ≈ p exp(2iψ) when P I (which is typically the case).

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 maps (Allys et al. 2019). In this subsection we use the WST to describe the statistical properties of the polarization related Article number, page 3 of 20 A&A proofs: manuscript no. Main maps discussed in the previous section: complex mapsQ + iŨ and exp(2iψ), and real map p.
Let us briefly recall the definition of the WST of a given 2D real field I(x) that is statistically homogeneous 4 . 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(x) with a set of J × Θ Morlet wavelets ψ j,θ with 0 ≤ j ≤ J − 1 and θ ∈ {kπ/Θ, 0 ≤ k ≤ Θ − 1}, with J and Θ two integers. These wavelets ψ j,θ are the result of the dilation by a factor 2 j and the rotation by an angle θ of a mother Morlet wavelet Ψ: 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 Ψ: with α and β 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 σ the standard deviation of the Gaussian envelope 5,6 . Convolving a field I(x) with a wavelet ψ j,θ corresponds to a local band-pass filtering of the field at frequencies centered on a mode k j,θ = k 0 2 − j cos (θ)ê x + sin (θ)ê y . J and Θ constants thus correspond respectively to the number of scales and angles that discretize the 2D Fourier space. For a statistically homogeneous real field I(x), normalized WST coefficients are defined as follows: with 0 ≤ j 1 ≤ J − 1 and 0 ≤ θ 1 < π, with 0 ≤ j 1 , j 2 ≤ J − 1 such that j 2 > j 1 and 0 ≤ θ 1 , θ 2 < π, where the averaging operator · here is simply an average over a map, meaning I = 1 A I(x)d 2 x with A the area of integration, and where the symbol stands for the convolution operator. Note that the number of WST coefficients depends on the J and Θ parameters: we have a singleS 0 coefficient, J × ΘS 1 coefficients, and Θ 2 × J × (J − 1)/2S 2 coefficients.
These coefficients can be interpreted in the following manner:S 0 is simply the mean of the field,S 1 ( j 1 , θ 1 ) is a measure of the amplitude of oscillation of the normalized I/ I field for modes centered on k j 1 ,θ 1 , and finallyS 2 ( j 1 , θ 1 , j 2 , θ 2 ) characterizes how the normalized amplitude of oscillation of the field for a mode k j 1 ,θ 1 is modulated by a mode of oscillation k j 2 ,θ 2 . Accordingly, theS 1 ( j 1 , θ 1 ) coefficient characterizes the amplitude at a single oriented scale ( j 1 , θ 1 ), while theS 2 ( j 1 , θ 1 , j 2 , θ 2 ) coefficient measures the coupling between two oriented scales ( j 1 , θ 1 ) and ( j 2 , θ 2 ).
Note that even if theS 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, theS 1 coefficients are computed with L 1 norms. One can however recover the power spectrum at a given scale from a quadratic sum ofS 1 andS 2 coefficients 7 . Notably, we expect a non-Gaussian field to have higherS 2 coefficients compared to those of a Gaussian field with identical power spectrum, and this should be counterbalanced by smallerS 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 Q + iŨ or exp(2iψ) as follows: where we still have 0 ≤ j 1 , j 2 ≤ J − 1 with j 2 > j 1 forS 2 coefficients, but here the main difference is that 0 ≤ θ 1 < 2π and 0 ≤ θ 2 < π. This new range of angles for θ 1 is directly related to the complex nature of theQ + iŨ variable. For a real signal I(x), one can easily show that |I ψ j 1 ,θ 1 +π | = |I ψ j 1 ,θ 1 |, which explains the range of θ 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 Ψ for angles between 0 and 2π (see Mallat 2012, for more details). Note that the range of angles for θ 2 is unchanged because |(Q + iŨ) ψ j 1 ,θ 1 | is a real signal. Thus, for a complex field, we end up with twice as manyS 1 andS 2 coefficients as for the WST of a real field. Let us note that this definition of the WST on complex polarization mapsQ + iŨ does not depend on the global phase of Q + iŨ maps, meaning that the WST coefficients associated tõ Q + iŨ 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), (17), and (18), while the statistical properties ofQ + iŨ and exp(2iψ) maps will be described with the WST coefficients defined in Eqs. (19), (20), and (21).

Simplification of the WST statistical descriptions for simulated Stokes maps with the RWST
The WST coefficients of polarization variablesQ + iŨ, exp(2iψ), 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. 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.

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) 3 box divided into 512 3 cells with periodic boundary conditions. At t = 0 the gas has uniform properties with a density n H = 1.5 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 ∼ 3.8 µG. In steady state, the turbulent forcing leads to an approximate velocity dispersion σ v ∼ 2.6 km/s. This gives a turnover time at large scale τ L ≈ 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 in-A&A proofs: manuscript no. Main dependent 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 δτ = 1.25 Myr which is roughly ten percent of τ L . The phenomenology of turbulence in the sense of Kolmogorov shows that the turnover time τ 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 we will be considering in the following, we find that δτ is about fives times smaller than the corresponding largest value of τ 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), (5), and (6) and choosing the z-axis of the cubes as the direction of integration. 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 σ ν 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 ν provided that the dust emission remains optically thin. Then from these Stokes maps we compute the associated maps P,Q,Ũ, p, cos(2ψ), and sin(2ψ) 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 isB ≈ B 0 . Because the direction of integration chosen to compute the previous Stokes maps is orthogonal to the direction of the mean magnetic field B in the simulation, we refer to the previous sets of maps using the ⊥ symbol. For instance, the p ⊥ 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,Q,Ũ, p, cos(2ψ), and sin(2ψ). We use the symbol to refer to these sets of maps.

Presentation of the data sets
We use for this work 8 data sets, each one comprising 14 statistically independent maps. The various maps presented in the last subsection define 6 data sets: In addition, we also analyze phase randomized data sets is an 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π). 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[·] operator is to severely damage the structure of the image. We use phase randomization as an approximation to Gaussianization 8 . 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 refers to the set of 14 maps produced by randomizing separatelỹ Q andŨ maps from theQ ⊥ +iŨ ⊥ data set (respectivelyQ +iŨ ). Technical details about phase randomization can be found in Appendix B, as well as an example of a phase randomized map for Figure 2 shows examples of maps for I ⊥ , Q ⊥ , U ⊥ , P ⊥ ,Q ⊥ , U ⊥ , p ⊥ , cos(2ψ) ⊥ , and sin(2ψ) ⊥ data sets, while Fig. 3 shows examples of maps for the corresponding 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 approximatelyQ + iŨ ≈ (Q + iU)/I so that p roughly behaves as the modulus ofQ + iŨ, and 2ψ roughly behaves as its complex argument. Next, we note that the magnitude of p is unsurprisingly much lower than that of p ⊥ 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ψ) ⊥ map that the anisotropy of the magnetic field in the simulation leads to values that are concentrated close to 1.

Regularity in the WST coefficients
We now focus on the WST coefficients associated with thẽ Q ⊥ + iŨ ⊥ data set, but the following reasoning remains valid for the other data sets, including those related to exp(2iψ) and p. For the 14 maps of theQ ⊥ + iŨ ⊥ data set, we compute the WST coefficients for J = 7 and Θ = 8 (see Eqs. (19), (20), and (21)). Note that 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).
These J and Θ values lead for each map to 112S 1 coefficients and 2688S 2 coefficients. Adding to this theS 0 coefficient, we end up with 2801 coefficients per map. We define meanS 0 , S 1 andS 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.
In Fig. 4, we see that for a fixed scale j 1 a smooth pattern emerges with respect to the angular variable θ 1 . Because of the definition ofS 1 coefficients, at fixed j 1 , log 2 (S 1 )( j 1 , ·) functions must be 2π-periodic, but here the smooth pattern appears to be πperiodic. While it is not surprising to get smooth angular patterns that reflect the regularity of physical processes, this π-periodicity with respect to the angular variable θ 1 was unexpected.
We can identify the same kind of regularity properties forS 2 coefficients in Fig. 5. SinceS 2 coefficients depend on four variables, two of which are angular variables θ 1 and θ 2 , it is more complicated to see this angular regularity, but Fig. 6 helps us to exhibit this regularity by plottingS 2 coefficients as a function of angular variables θ 1 and θ 2 when j 1 = 1 and j 2 = 3. Smooth angular patterns appear on this surface which is parameterized by θ 1 and θ 2 . In particular we see that θ 1 − θ 2 = c cuts of this surface for arbitrary constants c give roughly constantS 2 coefficients.
On this example we thus expect most of the angular modulation to depend on the θ 2 − θ 1 variable.

Definition of a RWST statistical description
The smooth periodic patterns identified in theQ ⊥ + iŨ ⊥ 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 θ 1 and θ 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 toQ + iŨ 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, theS 1 coefficients are written as: parameters. We also enforceŜ aniso 1 ( j 1 ) ≥ 0 in order to lift a degeneracy between theŜ aniso 1 ( j 1 ) and θ ref,1 parameters.Ŝ iso 1 ( j 1 ) quantifies the isotropic component of the data fluctuations at a scale j 1 , whileŜ aniso 1 ( j 1 ) defines a measure of the angular modulation of the coefficients at scale j 1 , introducing a reference angle θ ref,1 ( j 1 ) that defines a privileged direction in the maps 10 .
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 χ 2 r 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 11 . This mitigates most of the correlation between WST coefficients at the same scale j 1 . Figure 7 shows the χ 2 r values for both log 2 (S 1 ) and log 2 (S 2 ) fits for theQ ⊥ + iŨ ⊥ data set. The χ 2 r values are globally close to 1 except at low j 1 for χ 2,S 1 r and at low j 2 − j 1 for χ 2,S 2 r . 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. Note that these additional terms do not affect the values of the coefficients corresponding to the original RWST model. Finally, 10 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Ŝ aniso 1 ( j 1 ) cos 2 θ 1 − θ ref,1 ( j 1 ) to the isotropic amplitude of fluctuationsŜ iso 1 ( j 1 ). We use a base 2 logarithm to be consistent with the base 2 definition of scales j 1 and j 2 . 11 This decreases the effective number of degrees of freedom by one when computing χ 2 r values.  (23) and (24)) for theQ ⊥ + iŨ ⊥ data set. In the 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. the same RWST model applies to the other data sets used in this work, including exp(2iψ) and p data sets, and we get for all of them similar χ 2 r 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.
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 (includingS 0 coefficient in the description) compared to the original 2801 WST coefficients in the case of complex variablesQ + iŨ and exp(2iψ) 12 , thus providing a very large compression of the statistical information contained in the WST coefficients.

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.  Figure 8 presentsŜ iso 1 + log 2 (S 0 ) coefficients as a function of scale forQ ⊥ + iŨ ⊥ ,Q + iŨ , R[Q ⊥ + iŨ ⊥ ] and R[Q + iŨ ] data Fig. 8.Ŝ iso 1 ( j 1 ) + log 2 (S 0 ) RWST coefficients forQ ⊥ + iŨ ⊥ ,Q + iŨ , R[Q ⊥ + iŨ ⊥ ] and R[Q + iŨ ] data sets. We use solid lines for ⊥ data sets and dotted lines for 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 λ 1 = 2 j 1 +1 π/k 0 ). sets. Note that the log 2 (S 0 ) term cancels the normalization of thē S 1 WST coefficients defined in Eq. (20) in order to analyze the statistics ofQ + iŨ (and not the statistics of (Q + iŨ)/ |Q + iŨ| ). We see thatŜ iso 1 + log 2 (S 0 ) levels which measure an amplitude of fluctuations of the signal, per scale, are higher forQ ⊥ + iŨ ⊥ and R[Q ⊥ + iŨ ⊥ ] than for the corresponding data sets. This shows that the amplitudes differ depending on the orientation of the mean magnetic field with respect to the line of sight. We have more fluctuations withinQ + iŨ maps when the mean magnetic field is in the plane of the sky.
The differences between theQ + iŨ data set and its corresponding randomized counterpart R[Q + iŨ] (in both ⊥ and cases) illustrate the difference betweenS 1 coefficients and the power spectrum. Indeed,Q + iŨ and R[Q + iŨ] maps share the same power spectrum but have differentS 1 coefficients. The fact that the R[·] operator increases theS 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 theQ + iŨ data set and we therefore expect its second order coefficients to be higher compared to those of the corresponding randomized data set. Note that these differences wear off at large scales. We interpret this as statistical evidence that the non-Gaussianity ofQ + iŨ 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 compareQ + iŨŜ iso 1 coefficients to those of p and exp(2iψ) for respectively and ⊥ data sets. Since we havẽ Q + iŨ ≈ p exp(2iψ) we would like to compare the relative contributions ofŜ iso 1 coefficients ofQ+iŨ maps, but this seems more complicated than expected as we found out that a proper comparison throughŜ iso 1 coefficients highly depends on the choice of normalization of the WST coefficients for p.
S iso 1 coefficients roughly range from −4.5 to −3.0 for exp(2iψ) while they range from −5.8 to −4.2 for exp(2iψ) ⊥ . 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 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 ⊥ data sets than for ones. The p ⊥ data set constitutes an exception among the ⊥ 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 ⊥ 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 thẽ Q ⊥ + iŨ ⊥ data set corresponding to snapshots that are sufficiently distant in time to rightfully assume the independence (we choose 6 × δτ instead of δτ). The increase ofŜ aniso 1 coefficients remains significant for each map, which demonstrates that this trend is not a consequence of correlations among snapshots.
Reference angles θ ref,1 also presented in Fig. 10 show that the preferential direction identified for anisotropic ⊥ data sets is the direction corresponding to θ 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, i.e., along the y axis in Fig. 2. This corre- sponds 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 data set is well defined and approximately equal to π/4 for all scales while the anisotropy levels are close to zero. These values of θ 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 ⊥ data set, whereŜ aniso 1 coefficients are relatively small while θ ref,1 coefficients clearly deviate from zero.

R[Q+iŨ]
maps are Gaussian approximations ofQ+iŨ 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 theQ ⊥ + iŨ ⊥ andQ + iŨ data sets. The two dominant second order RWST coefficientsŜ iso,1 2 andŜ iso,2 2 presented in Fig. 11 display clearly distinct patterns between the original data sets and the randomized ones on the example of the ⊥ data sets. First,Ŝ iso,1 2 andŜ iso,2 2 coefficients are globally smaller for the R[Q + iŨ] data sets, which is in line with Fig. 11.Ŝ iso,1 2 ( j 1 , j 2 ) (top) andŜ iso,2 2 ( j 1 , j 2 ) (bottom) RWST coefficients forQ ⊥ + iŨ ⊥ and R[Q ⊥ + iŨ ⊥ ] 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.
what we had foreseen in Sect. 4.1. In addition,Ŝ iso,1 2 coefficients for R[Q+iŨ] show a scale invariance property:Ŝ iso,1 2 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).
S iso,2 2 coefficients also show two distinct trends: the coefficients quickly tend to zero when j 2 − j 1 increases for R[Q + iŨ] while coefficients tend to strictly positive values for the largest j 2 − j 1 values for non randomized data setsQ + iŨ. This result is related to the filamentary structure of the non randomized maps, because a filamentary structure involves a modulation of the WST coefficients as a function of the angle difference θ 2 −θ 1 . Note that we see the same trend in Fig. 12 for theŜ iso,2 2 coefficients of the other non randomized data sets which all present filamentary structures. For all data sets, and fixed j 1 values,Ŝ iso,2 2 coefficients decrease as a function of j 2 . This property seems to be general. Although note that we expect a signal due to the imprint of the wavelets in theŜ iso,2 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 theQ ⊥ + iŨ ⊥ data set and its randomized counterpart decrease for the highest j 1 values. Here again, this shows that the non-Gaussianity ofQ + iŨ maps decreases at large scales. This is consistent with what we have already pointed out in Sect. 4.1 forŜ iso 1 coefficients and we interpret this trend similarly.
Second order anisotropic coefficientsŜ aniso,1 2 ,Ŝ aniso,2 2 , and θ 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Ŝ iso,1 2 andŜ iso,2 2 . We do not show these coefficients here as we do not discuss them any further in this work.
A&A proofs: manuscript no. Main ( j 1 , j 2 ) (right column) RWST coefficients forQ + iŨ, exp(2iψ) and p data sets in the ⊥ 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.

Random syntheses ofQ + iŨ 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 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 WST/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. 2019;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.
In this section, we present random syntheticQ + iŨ 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.

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 mapsQ + iŨ. 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 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 =Q 0 + iŨ 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 r =Q r + iŨ r of the data set on which the RWST description is built serves as a reference for comparison with the synthetic maps. Because our WST/RWST analysis 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 r to address this gap in a deterministic way 13 . Formally, noting F [X 0 ] the Fourier transform of X 0 initially drawn as a realisation of a Gaussian white noise we set: where k min is the wave number corresponding to the largest scale probed by the WST. For J = 7, using Table 1   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 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 L of a complex image X =Q+iŨ as follows: with L WST [X] the loss function constraining the WST coefficients of X, L one−point [Q] (respectively L one−point [Ũ]) that constraining a few one-point moments ofQ (respectivelyŨ), and µ a weighting coefficient balancing the importance of one-point moments constraints with that of WST coefficients constraints. More specifically we set: 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 Θ = 8) we recall that we have N = 2801. These target WST coefficients are computed from the RWST coefficients derived from a given data set ofQ + iŨ 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: where M  We generate synthetic maps using theQ + iŨ RWST description as a target but we also provide equivalent results for theQ ⊥ + iŨ ⊥ 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 µ = 5 × 10 −8 and the algorithm stops when |L[X k ] − L[X k+1 ]| 10 −12 . These numerical values correspond Fig. 14. Power spectra of the syntheticQ + iŨ map shown in Fig. 13, compared to those of the reference map, forQ,Ũ, and |Q + iŨ| (left, middle and right respectively). The vertical dashed lines mark the wavelet central wave numbers corresponding to the scale indices j = 0, . . . , J − 1. 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 = π 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 L that does not properly constrain these modes. Figure 13 shows a syn-theticQ + iŨ 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 those of the reference maps 14 , but we also see at intermediate and small scales, which are the truly synthetic scales, consistent filamentary patterns and dynamic ranges.

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 thẽ Q + iŨ data set, similar large scales and some identical onepoint 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 ofQ,Ũ, andQ + iŨ 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 variablesQ,Ũ, andQ + iŨ 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 ofQ + iŨ maps as theS 1 andS 2 coefficients constrain the power spectrum ofQ+iŨ (see discussion in Sect. 2.3). However we point out that we correctly reproduce the power spectra ofQ andŨ taken separately. Still note that we did not investigate cases whereQ andŨ have very different power spectra.
In Fig. 15, we compare the one-point distribution functions of the reference maps and the synthetic maps forQ,Ũ, and |Q + iŨ|. 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 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 ofQ andŨ maps, we end up with a total of 158 coefficients to generate 512 × 512 complexQ + iŨ 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 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.

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 mapsQ + iŨ that are normalized by I + P. To study the contributions of the polarization fraction p and of the polarization angle ψ to the statistical properties of these complex Stokes maps, we also applied the WST to the corresponding maps of p and exp(2iψ). 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 Q + iŨ, p, and exp(2iψ) 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. ForQ + iŨ maps,Ŝ iso 1 + log 2 (S 0 ) coefficients are larger when the mean magnetic field is in the plane of the sky, while for exp(2iψ) mapsŜ iso 1 is larger when the mean magnetic field is along the line of sight.
-Ŝ aniso 1 coefficients quantify the statistical anisotropy of the maps. When the mean magnetic field is parallel to the line of sightŜ aniso 1 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 thẽ Q ⊥ + iŨ ⊥ and exp(2iψ) ⊥ maps.
-Second order RWST coefficients clearly exhibit the non-Gaussianity ofQ + iŨ maps (although this is also visible to a lesser extent in first order coefficients). While the randomized R[Q + iŨ] data sets present characteristic properties of scale invariant Gaussian fields (invariance ofŜ iso,1 2 as a function of the scale difference j 2 − j 1 and a quick decrease to zero ofŜ iso,2 2 as this scale difference increases), theŜ iso,1 2 and S iso,2 2 coefficients for the correspondingQ+iŨ data sets show clearly different patterns. In particular, the strictly positive values ofŜ iso,2 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 synthesizeQ+iŨ maps. Combining the RWST coefficients with additional constraints, we obtained synthetic maps that statistically match the original maps. The agreement demonstrates the comprehensiveness 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 will 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 2019). 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.   1 ( j 1 ) (left column),Ŝ lat,2 1 ( j 1 ) (middle column) andŜ iso,3 2 ( j 1 , j 2 ) (right column) RWST coefficients for theQ ⊥ + iŨ ⊥ data set. In thê S iso,3 2 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.