  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,
such as one-point statistics and power spectra. The wavelet scattering
transform (WST), a low-variance statistical descriptor of non-Gaussian
processes introduced in data science, opens a path towards this goal. We
applied the WST to noise-free maps of dust polarized thermal emission computed
from a numerical simulation of MHD turbulence. We analyzed normalized complex
Stokes maps and maps of the polarization fraction and polarization angle. The
WST yields a few thousand coefficients; some of them measure the amplitude of
the signal at a given scale, and the others characterize the couplings between
scales and orientations. The dependence on orientation can be fitted with the
reduced WST (RWST), an angular model introduced in previous works. The RWST
provides a statistical description of the polarization maps, quantifying their
multiscale properties in terms of isotropic and anisotropic contributions. It
allowed us to exhibit the dependence of the map structure on the orientation of
the mean magnetic field and to quantify the non-Gaussianity of the data. We
also used RWST coefficients, complemented by additional constraints, to
generate random synthetic maps with similar statistics. Their agreement with
the original maps demonstrates the comprehensiveness of the statistical
description provided by the RWST. This work is a step forward in the analysis
of observational data and the modeling of CMB foregrounds. We also release
PyWST, a Python package to perform WST/RWST analyses at:
https://github.com/bregaldo/pywst.
