Open Access
Issue
A&A
Volume 642, October 2020
Article Number A217
Number of page(s) 20
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202038044
Published online 22 October 2020

© B. Regaldo-Saint Blancard et al. 2020

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 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 ≃ 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 maps1. 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 PyWST2.

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 ν 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ν(Td) with Td the dust temperature, τν is the optical depth, p0 is an intrinsic polarization fraction parameter, γ is the anglethat 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 HEALPix3 (Górski et al. 2005) convention (see Fig. 1).

In the following, we assume uniform values for the intrinsic polarization p0 and the dust temperature Td, an optically thin medium, and an infinitesimaloptical depth equal to dτν = σνnHdz where σν is the dust opacity (assumed uniform) and nH is the total gas density. In this case, Eqs. (1)–(3) become:

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 nH, 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 follows4 :

These relations can be written in a more compact form introducing complex quantities: (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.

thumbnail Fig. 1

Definition of the angles γ and ϕ related tothe 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ϕ is counted positively clockwise from the x direction in the HEALPix convention.

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 , polarization fraction p and complex polarization angle exp(2).

We first define the normalized Stokes complex variable as follows: (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): (11)

which is independent of the density field nH.

The polarization fraction p and the polarization angle ψ follow the usual definitions: (12)

and (13)

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 nonlinearity 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 (2). This variable also avoids dealing with discontinuities of ψ maps that appear where the polarization angle approaches ± π∕2. Moreover, the statistical properties of exp(2) are easier to compare to those of since when PI (which is typically the case).

2.3 WST statistical description

The WST isa 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 maps discussed in the previous section: complex maps and exp (2), and real map p.

Let us briefly recall the definition of the WST of a given 2D real field I(x) that is statistically homogeneous5. 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 ≤ jJ − 1 and θ ∈{∕Θ, 0 ≤ k ≤ Θ − 1}, with J and Θ two integers. These wavelets ψj,θ are the result of the dilation by a factor 2j and the rotation by an angle θ of a mother Morlet wavelet Ψ: (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 Ψ: (15)

with α and β two constants that are adjusted to ensure a zero mean and a unit L1 norm, k0 = k0ex the wave vector of the plane wave factor, and σ the standard deviation of the Gaussian envelope6,7. 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 . 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 ≤ j1J − 1 and 0 ≤ θ1 < π, (18)

with 0 ≤ j1, j2J − 1 such that j2 > j1 and 0 ≤ θ1, θ2 < π, where the averaging operator ⟨⋅⟩ here is simply an average over a map, meaning with A the area of integration, and where the ⋆ symbol stands for the convolution operator. The number of WST coefficients depends on the J and Θ parameters: we have a single coefficient, J × Θ coefficients, and Θ2 × J × (J − 1)∕2 coefficients.

These coefficients can be interpreted in the following manner: is simply the mean of the field, is a measure of the amplitude of oscillation of the normalized I∕⟨I⟩ field for modes centered on , and finally characterizes how the normalized amplitude of oscillation of the field for a mode is modulated by a mode of oscillation . Accordingly, the coefficient characterizes the amplitude at a single oriented scale (j1, θ1), while the coefficient measures the coupling between two oriented scales (j1, θ1) and (j2, θ2).

We note that even if the 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 L2 norms of the wavelet transform, the coefficients are computed with L1 norms. One can however recover the power spectrum at a given scale from a quadratic sum of and coefficients8. Notably, we expect a non-Gaussian field to have higher coefficients compared to those of a Gaussian field with identical power spectrum, and this should be counterbalanced by smaller 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 or exp (2) as follows:

where we still have 0 ≤ j1, j2J − 1 with j2 > j1 for 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 the variable. Fora real signal I(x), one can easily show that , 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). We note that the range of angles for θ2 is unchanged because is a real signal. Thus, for a complex field, we end up with twice as many and coefficients as for the WST of a real field.

Let us note that this definition of the WST on complex polarization maps does not depend on the global phase of maps, meaning that the WST coefficients associated to 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 and exp(2) 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 , exp (2), 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.

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 box divided into 5123 cells with periodic boundary conditions. At t = 0 the gas has uniform properties with a density nH = 1.5 cm−3 and a temperature T = 8000 K, and the magnetic field is also uniform, with B0 = B0ex and B0 ~ 3.8 μG. In steady state, the turbulent forcing leads to an approximate velocity dispersion σv ~ 2.6 km s-1. 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 G0 = 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 δτ = 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 l2∕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 δτ 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)–(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 p0 = 0.2 and arbitrary values for σν and Td 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, , Ũ, 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 is . Because the direction of integration chosen to compute the previous Stokes maps is orthogonal to the direction of the mean magnetic field 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, , Ũ, p, cos (2ψ), and sin (2ψ). We use the ∥ 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: , , and . In addition, we also analyze phase randomized data sets and built respectively from the and data sets.

We define R[⋅] to be the operatorthat 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π). 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[⋅] operator is to severely damage the structure of the image. We use phase randomization as an approximation to Gaussianization9. 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 (respectively ) refers to the set of 14 maps produced byrandomizing separately and Ũ maps from the data set (respectively ). Technical details about phase randomization can be found in Appendix B, as well as an example of a phase randomized map for in Fig. B.1.

Figure 2 shows examples of maps for I, Q, U, P, , Ũ, p, , and 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 p0), hence we have approximately so that p roughly behaves as the modulus of , 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 map that the anisotropy of the magnetic field in the simulation leads to values that are concentrated close to 1.

thumbnail Fig. 2

Examples of I, Q, U, P, , Ũ, p, , and maps (from top to bottom and left to right) that are built from a given snapshot of the MHD simulation.

3.3 Regularity in the WST coefficients

We now focus on the WST coefficients associated with the data set, but the following reasoning remains valid for the other data sets, including those related to exp (2) and p.

For the 14 maps of the data set, we compute the WST coefficients for J = 7 and Θ = 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).

These J and Θ values lead for each map to 112 coefficients and 2688 coefficients. Adding to this the coefficient, we end up with 2801 coefficients per map. We define mean , and 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 coefficients and a representative subsample of coefficients (for j1 = 1) on a logarithmic scale, for . In both figures, we represent the multivariate functions and in a lexicographical order for the (j1, θ1) and (j1, θ1, j2, θ2) variables, respectively. Vertical gray lines help to mark increments of these variables.

In Fig. 4, we see that for a fixed scale j1 a smooth pattern emerges with respect to the angular variable θ1. Because of the definition of coefficients, at fixed j1, 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 nevertheless explain it by noticing that: (22)

because forMorlet wavelets we have 10. Hence, saying that functions are π-periodic amounts to saying that and have the same WST first order statistics.

We can identify the same kind of regularity properties for coefficients in Fig. 5. Since 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 plotting coefficients as a function of angular variables θ1 and θ2 when j1 = 1 and j2 = 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 constant coefficients. On this example we thus expect most of the angular modulation to depend on the θ2θ1 variable.

Table 1

Correspondence between scale index j and related wavelengths on simulated maps, both in pixel units and dimensional units.

thumbnail Fig. 3

Same as Fig. 2 but for the I, Q, U, P, , Ũ, p, , and data sets.

thumbnail Fig. 4

WST coefficients on a logarithmic scale for the data set, presented in a lexicographical order on (j1, θ1). Vertical dashed lines delimit distinct j1 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.

thumbnail Fig. 5

j1 = 1 selection of WST coefficients on a logarithmic scale for the data set presented in a lexicographical order on (j1, θ1, j2, θ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 θ1 and j2 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.

thumbnail Fig. 6

Surface representation of WST coefficients 1, θ1, j2 = 3, θ2) for the data set as a function of θ1 and θ2 variables only.

3.4 Definition of a RWST statistical description

The smooth periodic patterns identified in the 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 to 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 coefficients are written as: (23)

where , , and θref,1(j1) are the parameters of this angular model for scale j1. This model thus depends on 3 × J parameters. We also enforce in order to lift a degeneracy between the and θref,1 parameters. quantifies the isotropic component of the data fluctuations at a scale j1, while defines a measure of the angular modulation of the coefficients at scale j1, introducing a reference angle θref,1(j1) that defines a privileged direction in the maps11.

Similarly, the RWST model for coefficients reads: (24)

where , , , , and θref,2(j1, j2) are the parameters of this angular model for each pair of scales (j1, j2). As we have j2 > j1, we end up with 5 × J × (J − 1)∕2 parameters for this model. Here again we make sure that to avoid any parameter degeneracy. measures the overall amplitude ofcoupling between the scales j1 and j2. represents the amplitude of the modulation due to the relative orientation of the wavelets and , and we interpret it as a signature of filamentary structures at a given scale. Indeed for an oriented filamentary structure, we expect the coefficients to peak when θ2 = θ1 and to reach a minimum when θ2 = θ1 + π∕2. Finally, and are measures of the anisotropic properties of the data in second order WST coefficients, here decoupling θ1 and θ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 j1, and to the second order WST coefficients for every pair of scales (j1, j2) (with j2 > j1). The accuracy of these multiple fits is quantified with 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 j1 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 uncertainties12. This mitigates most of the correlation between WST coefficients at the same scale j1.

Figure 7 shows the values for both and fits for the data set. The values are globally close to 1 except at low j1 for and at low j2j1 for . 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 (2) and p data sets, and we get for all of them similar 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 j1 and j2j1 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 (including coefficient in the description) compared to the original 2801 WST coefficients in the case of complex variables and exp (2)13, thus providing a very large compression of the statistical information contained in the WST coefficients.

thumbnail Fig. 7

Reduced chi square (top) and (bottom) associated with the RWST fit of the WST coefficients (see Eqs. (23) and (24)) for the data set. Bottom panel: each curve corresponds to a fixed j1 value with j2 varying fromj1 + 1 to J − 1 = 6. For j1 = J − 2, the curve is reduced to a single dot on the figure. We use logarithmic scales for better visibility.

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

Figure 8 presents coefficients as a function of scale for , and data sets. We note that the term cancels the normalization of the WST coefficients defined in Eq. (20) in order to analyze the statistics of (and not the statistics of ). We see that levels whichmeasure an amplitude of fluctuations of the signal, per scale, are higher for and 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 within maps when the mean magnetic field is in the plane of the sky.

The differences between the data set and its corresponding randomized counterpart (in both ⊥ and ∥ cases) illustrate the difference between coefficients and the power spectrum. Indeed, and maps share the same power spectrum but have different coefficients. The fact that the R[⋅] operator increases the 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 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 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 coefficients to those of p and exp (2) for respectively ∥ and ⊥ data sets. Since we have we would like to compare the relative contributions of coefficients of maps, but this seems more complicated than expected as we found out that a proper comparison through coefficients highly depends on the choice of normalization of the WST coefficients for p.

coefficients roughly range from − 4.5 to − 3.0 for while they range from − 5.8 to − 4.2 for . 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 alongthe line of sight ( for p compared to for p) this feature is consistent with the anti-correlation observed with polarization data between the angle dispersion function and the polarization fraction p (Planck Collaboration Int. XX 2015; Planck Collaboration XII 2020; Fissel et al. 2016).

thumbnail Fig. 8

RWST coefficients for , and 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 = 2j1 +1πk0).

thumbnail Fig. 9

RWST coefficients for and p data sets in the ⊥ (top, solid lines) and ∥ (bottom, dotted lines) cases.

4.2 Anisotropic fluctuations in first order coefficients

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 the data set corresponding to snapshots that are sufficiently distant in time to rightfully assume the independence (we choose 6 × δτ instead of δτ). The increase of 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 j1. 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 ∥ 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 coefficients are relatively small while θref,1 coefficients clearly deviate from zero.

thumbnail Fig. 10

(top) and θref,1(j1) (bottom) RWST coefficients for , p, and data sets in the ⊥ (solid lines) and ∥ (dotted lines) cases.

thumbnail Fig. 11

(top) and (bottom) RWST coefficients for and data sets. Each curve corresponds to a fixed j1 value with j2 varying fromj1 + 1 to J − 1 = 6. For j1 = J − 2, the curve is reduced to a single dot on the figure.

4.3 Second order coefficients and non-Gaussianity of the data

maps are Gaussian approximations of 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 and data sets. The two dominant second orderRWST coefficients and 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, and coefficients are globally smaller for the data sets, which is in line with what we had foreseen in Sect. 4.1. In addition, coefficients for show a scale invariance property: coefficients only depend on the difference j2j1. 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).

coefficients also show two distinct trends: the coefficients quickly tend to zero when j2j1 increases for while coefficients tend to strictly positive values for the largest j2j1 values for non randomized data sets . 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. We see the same trend in Fig. 12 for the coefficients of the other non randomized data sets which all present filamentary structures. For all data sets, and fixed j1 values, coefficients decrease as a function of j2. This property seems to be general. We also expect a signal due to the imprint of the wavelets in the coefficients that would increase their values when j2 is close to j1 + 1.

We see in Fig. 11 that the differences between the data set and its randomized counterpart decrease for the highest j1 values. Here again, this shows that the non-Gaussianity of maps decreases at large scales. This is consistent with what we have already pointed out in Sect. 4.1 for coefficients and we interpret this trend similarly.

Second order anisotropic coefficients , 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 and . We do not show these coefficients here as we do not discuss them any further in this work.

5 Random syntheses of 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 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 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.

thumbnail Fig. 12

(left column) and (right column) RWST coefficients for , exp (2) and p data sets in the ⊥ case. Each curve corresponds to a fixed j1 value with j2 varying fromj1 + 1 to J − 1 = 6. For j1 = J − 2, the curve is reduced to a single dot on the figure.

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 . 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 in pixel space. Hence this optimization problem is defined in practice in a space of dimension 5122. Starting from a realization of a complex Gaussian white noise map , successive maps {Xk} 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 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 X0 with those of Xr to address this gap in a deterministic way14. Formally, noting the Fourier transform of X0 initially drawn as a realization of a Gaussian white noise we set: (25)

where kmin is the wave number corresponding to the largest scale probed by the WST. For J = 7, using Table 1 we have kmin ≈ 2π∕171 ≈ 0.037 pixel−1. In practice, 25 Fourier modes out of 5122 are set by this procedure. Even if in our case this value of kmin 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 Xr one-point statistics as a reference.

We define the loss function of a complex image as follows: (26)

with the loss function constrainingthe WST coefficients of X, (respectively ) that constraining a few one-point moments of (respectively Ũ), and μ a weighting coefficient balancing the importance of one-point moments constraints with that of WST coefficients constraints. More specifically we set: (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 Θ = 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 maps (of which Xr 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: (28)

where is an estimation of the variance of the distribution of pixels values in the map, is an estimation of its skewness, and is an estimation of its kurtosis.

We generate synthetic maps using the RWST description as a target but we also provide equivalent results for the 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 . 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 (kN = π 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 that does not properly constrain these modes. Figure 13 shows a synthetic 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 maps15, but we also see at intermediate and small scales, which are the truly synthetic scales, consistent filamentary patterns and dynamic ranges.

thumbnail Fig. 13

Synthesis of a map (right column) built from its corresponding RWST description with additional constraints on large-scale components and on a few one-point moments of and Ũ maps, to be consistent with those of a reference map (left column). The reference maps shown here are the same as in Fig. 3.

thumbnail Fig. 14

Power spectra of the synthetic map shown in Fig. 13, compared to those of the reference map, for , Ũ, and (left, middle and right respectively). The vertical dashed lines mark the wavelet central wave numbers corresponding to the scale indices j = 0, …, J − 1.

thumbnail Fig. 15

One-point distribution functions of the synthetic map shown in Fig. 13, compared to those of the reference map, for , Ũ, and (left, middle and right respectively).

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 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 , Ũ, and up to the one dimensional Nyquist wave number kN 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 , Ũ, and 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 kN. We interpret this as the result of poorly constrained modes close to kN in the optimization process. It is not surprising to reproduce the power spectrum of maps as the and coefficients constrain the power spectrum of (see discussion in Sect. 2.3). However we point out that we correctly reproduce the power spectra of and Ũ taken separately. Still we note that we did not investigate cases where and Ũ have very different power spectra.

In Fig. 15, we compare the one-point distribution functions of the reference maps and the synthetic maps for , Ũ, and . 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 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 and Ũ maps, we end up with a total of 158 coefficients to generate 512 × 512 complex 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.

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 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 (2). 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 , p, and exp (2) 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 maps, coefficients are larger when the mean magnetic field is in the plane of the sky, while for exp(2) maps is larger when the mean magnetic field is along the line of sight.

  • coefficients quantify the statistical anisotropy of the maps. When the mean magnetic field is parallel to the line of sight coefficients are negligible, while when the mean magnetic field is in the plane of the sky they allow us to identify the directionof 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 and maps.

  • Second order RWST coefficients clearly exhibit the non-Gaussianity of maps (although this is also visible to a lesser extent in first order coefficients). While the randomized data setspresent characteristic properties of scale invariant Gaussian fields (invariance of as a function of the scale difference j2j1 and a quick decrease to zero of as this scale difference increases), the and coefficients for the corresponding data setsshow clearly different patterns. In particular, the strictly positive values of at large j2j1 are interpreted as signatures of the filamentary structure of the maps.

We have used the RWST approach to synthesize 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 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).

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 coefficients at a given scale j1 and write the corresponding angular model16. If one assumes that there are no more than one privileged direction θref,1(j1) in the maps we are dealing with, we can write generally as a Fourier series expansion using the natural 2π-periodicity of this function: (A.1)

Assuming a mirror symmetry with respect to the potential reference direction, we expect {bk } coefficients to vanish (which is the case in practice for our data). Adding to that the π-periodicity identified in Sect. 3.3 the angular model reduces to: (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 and 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 with Ω = 0, N − 12. Its discrete Fourier transform reads for a given mode k = (kx, ky) ∈ Ω: (B.1)

where rk = xkx + yky. We also recall the inverse discrete Fourier transform relation for the chosen convention: (B.2)

The phase randomization of f consists in defining a new field g in Fourier space such that: (B.3)

for every mode k ∈ Ω and where ϕ is a realization of a uniform random phase such as defined in Galerne et al. (2011), that is, defining Ω0 = {(0, N∕2), (N∕2, 0), (N∕2, N∕2)} if N is even, and Ω0 = ∅ otherwise, ϕ verifies:

  • 1.

    k ∈ Ω \ Ω0, ϕ(−k) = − ϕ(k) (we extend the domain of ϕ to using periodic boundary conditions).

  • 2.

    k ∈ Ω, ϕ(k) is drawn uniformly and independently in [0, 2π) (the independence holds to the extent of the first relation).

  • 3.

    k ∈ Ω0, ϕ(k) is drawn uniformly and independently in the set {0, π}.

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 ϕ(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 map in Fig. B.1.

thumbnail Fig. B.1

Example of a phase randomized map (right) next to its corresponding original map (left).

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 coefficients: (C.1)

where the additional lattice terms and quantify angular modulations that are respectively π∕2 and π∕4-periodic and aligned with the main directions of the lattice. Similarly the enhanced RWST model of coefficients is: (C.2)

where the additional term measures the amplitude ofan additional harmonic of the θ1θ2 angular modulation that is π∕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 and 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 improved normalized residuals compared to the previous ones shown in Figs. 4 and 5. In particular at j1 = 0 we no longer observe the strong angular pattern seen in Fig. 4.

thumbnail Fig. C.1

Reduced chi square (top) and (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 data set. Each curve in the plot corresponds to a fixed j1 value while j2 ranges from j1 + 1 to J − 1 = 6. For j1 = J − 2, the curve is reduced to a single dot on the figure. We use logarithmic scales for better visibility.

thumbnail Fig. C.2

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

For completeness, we show in Fig. C.4 the additional RWST terms given for the example of the data set. C.2:

C.2 Syntheses for Q̃ +

We show in Figs. C.5C.7 the maps, power spectra, and distribution functions of synthetic maps in the ⊥ case that are analogous to those of Figs. 1315.

thumbnail Fig. C.3

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

thumbnail Fig. C.4

(left column), (middle column) and (right column) RWST coefficients for the data set. In the plot, each curve corresponds to a fixed j1 value while j2 ranges from j1 + 1 to J − 1 = 6. For j1 = J − 2, the curve is reduced to a single dot on the figure.

thumbnail Fig. C.5

Same as Fig. 13 but for the ⊥ case. The reference maps shown here are the same maps as in Fig. 2.

thumbnail Fig. C.6

Same as Fig. 14 but for the ⊥ case.

thumbnail Fig. C.7

Same as Fig. 15 but for the ⊥ case.

References

  1. Abazajian, K., Addison, G., Adshead, P., et al. 2019, ArXiv e-prints, [arXiv:1907.04473] [Google Scholar]
  2. Adak, D., Ghosh, T., Boulanger, F., et al. 2020, A&A 640, A100 [CrossRef] [EDP Sciences] [Google Scholar]
  3. Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, J. Cosmol. Astropart. Phys., 2019, 056 [CrossRef] [Google Scholar]
  4. Allys, E., Levrier, F., Zhang, S., et al. 2019, A&A, 629, A115 [CrossRef] [EDP Sciences] [Google Scholar]
  5. Andersson, B. G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501 [NASA ADS] [CrossRef] [Google Scholar]
  6. Aylor, K., Haq, M., Knox, L., Hezaveh, Y., & Perreault-Levasseur, L. 2019, ArXiv e-prints, [arXiv:1909.06467] [Google Scholar]
  7. BICEP2 Collaboration 2018, Phys. Rev. Lett., 121, 221301 [CrossRef] [PubMed] [Google Scholar]
  8. BICEP2/Keck Array and Planck Collaborations 2015, Phys. Rev. Lett., 114, 101301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  9. Bruna, J., & Mallat, S. 2013, IEEE Trans. Pattern Anal. Mach. Intell., 35, 1872 [CrossRef] [Google Scholar]
  10. Bruna, J., & Mallat, S. 2019, Math. Stat. Learn., 1, 257 [CrossRef] [Google Scholar]
  11. Bruna, J., Mallat, S., Bacry, E., & Muzy, J.-F. 2015, Ann. Stat., 43, 323 [CrossRef] [Google Scholar]
  12. Byrd, R., Lu, P., Nocedal, J., & Zhu, C. 1995, SIAM J. Sci. Comput., 16, 1190 [CrossRef] [MathSciNet] [Google Scholar]
  13. Chuss, D. T., Andersson, B. G., Bally, J., et al. 2019, ApJ, 872, 187 [NASA ADS] [CrossRef] [Google Scholar]
  14. Clark, S. E., & Hensley, B. S. 2019, ApJ, 887, 136 [CrossRef] [Google Scholar]
  15. Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton: Princeton University Press) [Google Scholar]
  16. Draine, B. T., & Fraisse, A. A. 2009, ApJ, 696, 1 [NASA ADS] [CrossRef] [Google Scholar]
  17. Farge, M. 1992, Ann. Rev. Fluid Mech., 24, 395 [NASA ADS] [CrossRef] [Google Scholar]
  18. Federrath, C., Roman-Duval, J., Klessen, R., Schmidt, W., & Mac Low, M. M. 2009, A&A, 512, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Fissel, L. M., Ade, P. A. R., Angilè, F. E., et al. 2016, ApJ, 824, 134 [NASA ADS] [CrossRef] [Google Scholar]
  20. Fletcher, R. 1987, Practical Methods of Optimization, 2nd edn. (Chichester: John Wiley & Sons) [Google Scholar]
  21. Fraisse, A. A., Ade, P. A. R., Amiri, M., et al. 2013, J. Cosmol. Astropart. Phys., 2013, 047 [NASA ADS] [CrossRef] [Google Scholar]
  22. Frisch, U. 1995, Turbulence, the Legacy of A.N. Kolmogorov (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
  23. Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Galerne, B., Gousseau, Y., & Morel, J.-M. 2011, IEEE Trans. Image Process., 20, 257 [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
  25. Ghosh, T., Boulanger, F., Martin, P. G., et al. 2017, A&A, 601, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
  27. Guillet, V., Fanciullo, L., Verstraete, L., et al. 2018, A&A, 610, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Guth, A. H. 1981, Phys. Rev. D, 23, 347 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  29. Haverkorn, M., Katgert, P., & de Bruyn, A. G. 2004, A&A, 427, 549 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Hennebelle, P., & Falgarone, E. 2012, A&ARv, 20, 55 [NASA ADS] [CrossRef] [Google Scholar]
  31. Hennebelle, P., & Inutsuka, S. I. 2019, Front. Astron. Space Sci., 6 [Google Scholar]
  32. Hull, C. L. H., Mocz, P., Burkhart, B., et al. 2017, ApJ, 842, L9 [NASA ADS] [CrossRef] [Google Scholar]
  33. Iffrig, O., & Hennebelle, P. 2017, A&A, 604, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Ishino, H., Akiba, Y., Arnold, K., et al. 2016, SPIE Conf. Ser., 9904, 99040X [Google Scholar]
  35. Julesz, B. 1981, Nature, 290, 91 [CrossRef] [Google Scholar]
  36. Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. Lett., 78, 2058 [NASA ADS] [CrossRef] [Google Scholar]
  37. Kim, C.-G., Choi, S. K., & Flauger, R. 2019, ApJ, 880, 106 [NASA ADS] [CrossRef] [Google Scholar]
  38. Krizhevsky, A., Sutskever, I., & Hinton, G. E. 2012, Adv. Neural Inf. Process. Syst. 25, 1097 [Google Scholar]
  39. Kroese, D. P.,& Botev, Z. I. 2015, Spatial Process Simulation, ed. V. Schmidt (Cham: Springer International Publishing), 369 [Google Scholar]
  40. Landi Degl’innocenti, E., & Landolfi, M. 2004, Polarization in Spectral Lines (Netherlands: Springer Netherlands) [CrossRef] [Google Scholar]
  41. LeCun, Y., Kavukcuoglu, K., & Farabet, C. 2010, in Proceedings of 2010 IEEE International Symposium on Circuits and Systems, 253 [CrossRef] [Google Scholar]
  42. Levrier, F., Neveu, J., Falgarone, E., et al. 2018, A&A, 614, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Linde, A. D. 1982, Phys. Lett. B, 108, 389 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  44. Mallat, S. 2012, Commun. Pure Appl. Math., 65, 1331 [CrossRef] [Google Scholar]
  45. Mangilli, A., Aumont, J., Bernard, J. P., et al. 2019, A&A, 630, A74 [CrossRef] [EDP Sciences] [Google Scholar]
  46. Miville-Deschênes, M. A., Martin, P. G., Abergel, A., et al. 2010, A&A, 518, L104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Naess, S., Hasselfield, M., McMahon, J., et al. 2014, J. Cosmol. Astropart. Phys., 2014, 007 [NASA ADS] [CrossRef] [Google Scholar]
  48. Oppenheim, A. V., & Lim, J. S. 1981, Proc. IEEE, 69, 529 [CrossRef] [Google Scholar]
  49. Planck Collaboration XI. 2020, A&A, 641, A11 [CrossRef] [EDP Sciences] [Google Scholar]
  50. Planck Collaboration XII. 2020, A&A, 641, A12 [CrossRef] [EDP Sciences] [Google Scholar]
  51. Planck Collaboration Int. XIX. 2015, A&A, 576, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Planck Collaboration Int. XX. 2015, A&A, 576, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Planck Collaboration Int. XXX. 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Planck Collaboration Int. XLII. 2016, A&A, 596, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Reissl, S., Guillet, V., Brauer, R., et al. 2020, A&A 640, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Ritacco, A., Adam, R., Ade, P., et al. 2020, Eur. Phys. J. Web Conf., 228, 00022 [CrossRef] [Google Scholar]
  57. Schmidt, W., Hillebrandt, W., & Niemeyer, J. C. 2006, Comput. Fluids, 35, 353 [CrossRef] [Google Scholar]
  58. Siebenmorgen, R., Voshchinnikov, N. V., & Bagnulo, S. 2014, A&A, 561, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Sutin, B. M., Alvarez, M., Battaglia, N., et al. 2018, SPIE Conf. Ser., 10698, 106984F [Google Scholar]
  60. Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Vansyngel, F., Boulanger, F., Ghosh, T., et al. 2017, A&A, 603, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Wandelt, B. D. 2013, Astrostatistical Challenges for the New Astronomy, ed. J. M. Hilbe (New York: Springer New York), 87 [CrossRef] [Google Scholar]

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).

4

These relations hold for an angle φ 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.

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).

6

In practice we choose k0 = 3π∕4 pixel−1 and σ = 0.8 pixel as in Bruna & Mallat (2013).

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.

8

See Eq. (9) in Allys et al. (2019).

9

Stationary Gaussian random fields (GRFs) do have uniformly distributed phases on [0, 2π) but this property alone does not define them (see Wandelt 2013).

10

The * symbol stands for the complex conjugate.

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 to the isotropic amplitude of fluctuations . We use a base 2 logarithm to be consistent with the base 2 definition of scales j1 and j2.

12

This decreases the effective number of degrees of freedom by one when computing values.

13

For p maps, the WST descriptions consist of 1401 coefficients.

14

We note that such an approach cannot supply the statistical information concerning the couplings between the largest scales with smaller ones.

15

Since the WST constraints leave the largest scales of X0 unchanged, this shows that the impact of the one-point moments constraints on these scales is minor.

16

The same kind of reasoning holds for coefficients.

All Tables

Table 1

Correspondence between scale index j and related wavelengths on simulated maps, both in pixel units and dimensional units.

All Figures

thumbnail Fig. 1

Definition of the angles γ and ϕ related tothe 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ϕ is counted positively clockwise from the x direction in the HEALPix convention.

In the text
thumbnail Fig. 2

Examples of I, Q, U, P, , Ũ, p, , and maps (from top to bottom and left to right) that are built from a given snapshot of the MHD simulation.

In the text
thumbnail Fig. 3

Same as Fig. 2 but for the I, Q, U, P, , Ũ, p, , and data sets.

In the text
thumbnail Fig. 4

WST coefficients on a logarithmic scale for the data set, presented in a lexicographical order on (j1, θ1). Vertical dashed lines delimit distinct j1 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.

In the text
thumbnail Fig. 5

j1 = 1 selection of WST coefficients on a logarithmic scale for the data set presented in a lexicographical order on (j1, θ1, j2, θ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 θ1 and j2 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.

In the text
thumbnail Fig. 6

Surface representation of WST coefficients 1, θ1, j2 = 3, θ2) for the data set as a function of θ1 and θ2 variables only.

In the text
thumbnail Fig. 7

Reduced chi square (top) and (bottom) associated with the RWST fit of the WST coefficients (see Eqs. (23) and (24)) for the data set. Bottom panel: each curve corresponds to a fixed j1 value with j2 varying fromj1 + 1 to J − 1 = 6. For j1 = J − 2, the curve is reduced to a single dot on the figure. We use logarithmic scales for better visibility.

In the text
thumbnail Fig. 8

RWST coefficients for , and 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 = 2j1 +1πk0).

In the text
thumbnail Fig. 9

RWST coefficients for and p data sets in the ⊥ (top, solid lines) and ∥ (bottom, dotted lines) cases.

In the text
thumbnail Fig. 10

(top) and θref,1(j1) (bottom) RWST coefficients for , p, and data sets in the ⊥ (solid lines) and ∥ (dotted lines) cases.

In the text
thumbnail Fig. 11

(top) and (bottom) RWST coefficients for and data sets. Each curve corresponds to a fixed j1 value with j2 varying fromj1 + 1 to J − 1 = 6. For j1 = J − 2, the curve is reduced to a single dot on the figure.

In the text
thumbnail Fig. 12

(left column) and (right column) RWST coefficients for , exp (2) and p data sets in the ⊥ case. Each curve corresponds to a fixed j1 value with j2 varying fromj1 + 1 to J − 1 = 6. For j1 = J − 2, the curve is reduced to a single dot on the figure.

In the text
thumbnail Fig. 13

Synthesis of a map (right column) built from its corresponding RWST description with additional constraints on large-scale components and on a few one-point moments of and Ũ maps, to be consistent with those of a reference map (left column). The reference maps shown here are the same as in Fig. 3.

In the text
thumbnail Fig. 14

Power spectra of the synthetic map shown in Fig. 13, compared to those of the reference map, for , Ũ, and (left, middle and right respectively). The vertical dashed lines mark the wavelet central wave numbers corresponding to the scale indices j = 0, …, J − 1.

In the text
thumbnail Fig. 15

One-point distribution functions of the synthetic map shown in Fig. 13, compared to those of the reference map, for , Ũ, and (left, middle and right respectively).

In the text
thumbnail Fig. B.1

Example of a phase randomized map (right) next to its corresponding original map (left).

In the text
thumbnail Fig. C.1

Reduced chi square (top) and (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 data set. Each curve in the plot corresponds to a fixed j1 value while j2 ranges from j1 + 1 to J − 1 = 6. For j1 = J − 2, the curve is reduced to a single dot on the figure. We use logarithmic scales for better visibility.

In the text
thumbnail Fig. C.2

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

In the text
thumbnail Fig. C.3

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

In the text
thumbnail Fig. C.4

(left column), (middle column) and (right column) RWST coefficients for the data set. In the plot, each curve corresponds to a fixed j1 value while j2 ranges from j1 + 1 to J − 1 = 6. For j1 = J − 2, the curve is reduced to a single dot on the figure.

In the text
thumbnail Fig. C.5

Same as Fig. 13 but for the ⊥ case. The reference maps shown here are the same maps as in Fig. 2.

In the text
thumbnail Fig. C.6

Same as Fig. 14 but for the ⊥ case.

In the text
thumbnail Fig. C.7

Same as Fig. 15 but for the ⊥ case.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.