  [abridged] The interstellar medium is now widely recognized to display
features ascribable to magnetized turbulence. With the public release of Planck
data and the current balloon-borne and ground-based experiments, the growing
amount of data tracing the polarized thermal emission from Galactic dust in the
submillimetre provides choice diagnostics to constrain the properties of this
magnetized turbulence. We aim to constrain these properties in a statistical
way, focusing in particular on the power spectral index of the turbulent
component of the interstellar magnetic field in a diffuse molecular cloud, the
Polaris Flare. We present an analysis framework which is based on simulating
polarized thermal dust emission maps using model dust density (proportional to
gas density) and magnetic field cubes, integrated along the line of sight, and
comparing these statistically to actual data. The model fields are derived from
fBm processes, which allow a precise control of their one- and two-point
statistics. We explore the nine-dimensional parameter space of these models
through a MCMC analysis, which yields best-fitting parameters and associated
uncertainties. We find that the power spectrum of the turbulent component of
the magnetic field in the Polaris Flare molecular cloud scales with wavenumber
as a power law with a spectral index $2.8\pm 0.2$. It complements a uniform
field whose norm in the POS is approximately twice the norm of the fluctuations
of the turbulent component. The density field is well represented by a
log-normally distributed field with a mean gas density $40\,\mathrm{cm}^{-3}$
and a power spectrum with as spectral index $1.7^{+0.4}_{-0.3}$. The agreement
between the Planck data and the simulated maps for these best-fitting
parameters is quantified by a $\chi^2$ value that is only slightly larger than
unity.
