---
author:
- 'F. Levrier'
- 'J. Neveu'
- 'E. Falgarone'
- 'F. Boulanger'
- 'A. Bracco'
- 'T. Ghosh'
- 'F. Vansyngel'
bibliography:
- 'biblio.bib'
- 'Planck\_bib.bib'
date: 'Received 19 October 2017; accepted 23 February 2018'
title: 'On the statistics of the polarized submillimetre emission maps from thermal dust in the turbulent, magnetized, diffuse ISM'
---

Introduction
============

In recent years, a number of experiments have dramatically increased the amount of data pertaining to the polarized thermal emission from Galactic dust in the submillimetre [e.g. @matthews-et-al-2009; @ward-thompson-et-al-2009; @dotson-et-al-2010; @bierman-et-al-2011; @vaillancourt-matthews-2012; @hull-et-al-2014; @koch-et-al-2014; @fissel-et-al-2016].  Chief among these is [[*Planck*]{}]{}[^1], which provided the first full-sky map of this emission, leading to several breakthrough results.  It was thus found that the polarization fraction $p$ in diffuse regions of the sky can reach values above 20% [@planck2014-XIX], confirming results previously obtained over one fifth of the sky by the [*Archeops*]{} balloon-borne experiment [@benoit-et-al-2004; @ponthieu-et-al-2005].  Furthermore, the polarization fraction is anti-correlated with the local dispersion $\mathcal{S}$ of polarization angles $\psi$ [@planck2014-XIX; @planck2014-XX], and the decrease of the maximum observed $p$ with increasing gas column density ${N_\mathrm{H}}$ may be fully accounted for, at the scales probed by [[*Planck*]{}]{} (5 at 353GHz), by the tangling of the magnetic field on the line of sight (LOS) [@planck2014-XX].  Similar anti-correlations were found by the BLASTPol experiment [@fissel-et-al-2016] at higher angular resolution (a few tens of arcseconds) towards a single Galactic molecular cloud (Vela C).  At 10 scales, and over a larger sample of clouds, [[*Planck*]{}]{} data showed that the relative orientation of the interstellar magnetic field ${\boldsymbol{B}}$ and filamentary structures of dust emission is consistent with simulated observations derived from numerical simulations of sub- or trans-Alfvénic MHD turbulence [@planck2015-XXXV], and starlight polarization data in extinction yield similar diagnostics [@soler-et-al-2016].  In diffuse regions, the preferential alignment of filamentary structures with the magnetic field [@planck2014-XXXII; @planck2015-XXXVIII] is linked to the measured asymmetry between the power spectral amplitudes of the so-called E- and B-modes of polarized emission.  Finally, measurements of the spatial power spectrum of polarized dust emission showed that it must be taken into account in order to obtain reliable estimates of the cosmological polarization signal [@planck2014-XXX]. 

With this wealth of data, we may be able to put constraints on models of magnetized turbulence in the interstellar medium (ISM), provided we can extract the relevant information from polarization maps.  Of particular interest are the statistical properties of the Galactic magnetic field (GMF) ${\boldsymbol{B}}$. Let us write it as a sum ${\boldsymbol{B}}={\boldsymbol{B}}_0+{\boldsymbol{B}}_t$ of a uniform, large-scale component ${\boldsymbol{B}}_0$, and a turbulent component ${\boldsymbol{B}}_t$ with a null spatial average, $\langle{\boldsymbol{B}}_t\rangle=\boldsymbol{0}$. The statistical properties in question are then essentially modelled by two quantities, i) the ratio of the turbulent component to the mean, $y_B=\sigma_B/B_0$, where $\sigma_B^2=\langle {\boldsymbol{B}}_t^2\rangle$ and $B_0=||{\boldsymbol{B}}_0||$, and ii) the spectral index $\beta_B$, which characterizes the distribution of power of ${\boldsymbol{B}}_t$ across spatial scales, through the relationship $P(k)\propto k^{-\beta_B}$, where $k$ is the wavenumber and $P(k)$ is the power spectrum[^2]. 

As already mentioned, [@planck2015-XXXV] studied the relative orientation between the magnetic field, probed by polarized thermal dust emission, and filaments of matter in and around nearby molecular clouds.  They found that this relative orientation changes, from mostly parallel to mostly perpendicular, as the total gas column density ${N_\mathrm{H}}$ increases, which is a trend observed in simulations of trans-Alfvénic or sub-Alfvénic MHD turbulence [@soler-et-al-2013].  Using the Davis-Chandraskehar-Fermi method [@chandrasekhar-fermi-1953] improved by [@falceta-goncalves-et-al-2008] and [@hildebrand-et-al-2009], and from their results, we can estimate the ratio $y_B$ to be in the range 0.3-0.7. [ @planck2014-XXXII] studied that same relative orientation in the diffuse ISM at intermediate and high Galactic latitudes, and their estimate of $y_B$ is in the range 0.6-1.0 with a preferred value at 0.8.  These estimates are confirmed in [@planck2016-XLIV], which presents a fit of the distributions of polarization angles and fractions observed by [[*Planck*]{}]{} towards the southern Galactic cap.  They use a model of the GMF involving a uniform large-scale field ${\boldsymbol{B}}_0$ and a small number ($N_l\simeq 7$) of independent “polarization layers” on the line of sight, each of which accounts for a fraction $1/N_l$ of the total unpolarized emission.  Within each layer, the turbulent component of the magnetic field ${\boldsymbol{B}}_t$, which is used to compute synthetic Stokes $Q$ and $U$ maps, is an independent realization of a Gaussian 2D random field with a prescribed spectral index $\beta_B$. Through these fits, they confirm the near equipartition of large-scale and turbulent components of ${\boldsymbol{B}}$, with $y_B\simeq 0.9$. They also provide a rough estimate of the magnetic field’s spectral index $\beta_B$ in the range 2-3.  This work was complemented in [@vansyngel-et-al-2017], using the same framework, but including observational constraints on the power spectra of polarized thermal dust emission.  These authors were able to constrain $\beta_B\simeq 2.5$, an exponent which is compatible with the rough estimate of [@planck2016-XLIV], close to that measured for the total intensity of the dust emission.  We note that their exploration of the parameter space does not allow for an estimation of the uncertainty on $\beta_B$.

In [@planck2014-XXXII], [@planck2016-XLIV] and [@vansyngel-et-al-2017], the description of structures, in both dust density and magnetic field, along the LOS is reduced to the bare minimum, while statistical properties in the plane of the sky (POS) are modelled through $y_B$ and $\beta_B$. Orthogonal approaches have also been pursued [e.g. @miville-deschenes-et-al-2008; @odea-et-al-2012], in which the turbulent component of the magnetic field is modelled along each LOS independently from the neighbouring ones, as a realization of a one-dimensional Gaussian random field with a power-law power spectrum.  In this type of approach there is no correlation from pixel to pixel on the sky, and such studies seek to exploit the depolarization along the LOS, rather than spatial correlations in the POS, to constrain statistical properties of the interstellar magnetic field. 

We seek to explore another avenue, taking into account statistical correlation properties of ${\boldsymbol{B}}$ in all three dimensions, as well as properties of the dust density field, building on methods developed in [@planck2014-XX] to compare [[*Planck*]{}]{} data with synthetic polarization maps.  In that paper, the synthetic maps were computed from data cubes of dust density ${n_\mathrm{d}}$ and magnetic field ${\boldsymbol{B}}$ produced by numerical simulations of MHD turbulence.  One could think to generalize this approach, taking advantage of the ever-increasing set of such simulations [see, e.g., @hennebelle-et-al-2008; @hennebelle-2013; @hennebelle-iffrig-2014; @inutsuka-et-al-2015; @seifried-walch-2015].  However, this would be impractical for two main reasons : i) these simulations often have a limited inertial range over which the power spectrum has a power-law behaviour, and ii) a systematic study exploring a wide range of physical parameters with sufficient sampling is not possible due to the computational cost of each simulation. 

We therefore propose an alternative approach, which is to build simple, approximate, three-dimensional models for the dust density ${n_\mathrm{d}}$ and the magnetic field ${\boldsymbol{B}}$, allowing us to perfectly control the statistical properties of these 3D fields, and to fully explore the space of parameters characterizing them.  With this approach, we are able to perform a statistically significant number of simulated polarization maps for each set of parameters.  Actual observations may then be compared to these simulated maps, using least-square analysis methods, to extract best-fitting parameters, in particular the spectral index of the magnetic field, $\beta_B$, and the ratio of turbulent to regular field, $y_B$.

The paper is organized as follows : Sec.  \[sec:buildingmaps\] presents the method used to build simulated thermal dust polarized emission maps using prescribed statistical properties for ${n_\mathrm{d}}$ and ${\boldsymbol{B}}$. Observables derived from these maps, serving as statistical diagnostics of the input parameters, are presented in Sec.  \[sec:observables\].  In Sec.  \[sec:method\], we describe the analysis method devised to explore the space of input parameters for a given set of polarization maps.  The validation of the method and its application to actual observations of polarized dust emission from the Polaris Flare molecular cloud observed by [[*Planck*]{}]{} are given in Sec.  \[sec:results\].  Finally, Sec.  \[sec:discussion\] discusses our results and offers conclusions.  Several appendices complement our work.  Appendix \[sec:appendix:nHproperties\] presents further statistical properties of the model dust density fields.  Appendix \[sec:appendix:L2distance\] details the likelihood used in the MCMC analysis.  Finally, appendix \[sec:appendix:chi2\] details the $\chi^2$ parameter used to estimate the goodness-of-fit. 

Building synthetic polarized emission maps {#sec:buildingmaps}
==========================================

In this section, we first describe the synthetic dust density and magnetic field cubes we use in our analysis, then explain how simulated polarized emission maps are built from these cubes. 

Fractional Brownian motions
---------------------------

The basic ingredients to synthetise polarized thermal dust emission maps are three-dimensional cubes of dust density ${n_\mathrm{d}}$ and magnetic field ${\boldsymbol{B}}$, which we build using fractional Brownian motions (fBm) [@falconer-1990].  An $N$-dimensional fBm $X$ is a random field defined on $\mathbb{R}^N$ such that $\langle\left[X\left(\boldsymbol{r}_2\right)-X\left(\boldsymbol{r}_1\right)\right]^2\rangle\propto{||}\boldsymbol{r}_2-\boldsymbol{r}_1{||}^{2H}$, for any pair of points $(\boldsymbol{r}_1,\boldsymbol{r}_2)$. $H$ is called the Hurst exponent.  These fBm fields are usually built in Fourier space[^3], $$\label{eq:Xfbm}
\widetilde{X}(\boldsymbol{k})=A(\boldsymbol{k})\exp{\left[i\phi_X(\boldsymbol{k})\right]},$$ by specifying amplitudes that scale as a power-law of the wavenumber $k=||\boldsymbol{k}||$, $$A(\boldsymbol{k})=A_0k^{-\beta_X/2},$$ with $\beta_X=2H+N$ the spectral index, and phases drawn from a uniform random distribution in $[-\pi,\pi]$, subject to the constraint $\phi_X(-\boldsymbol{k})=-\phi_X(\boldsymbol{k})$ so that $X$ is real-valued.  Their power spectra are therefore power laws, $$P_X(k)=\left<\left|\widetilde{X}(\boldsymbol{k})\right|^2\right>_{||\boldsymbol{k}||=k}\propto k^{-\beta_X},$$ where the average is taken over the locus of constant wavenumber $k$ in Fourier space.  Such fields have been used previously as toy models for the fractal structure of molecular clouds, in both density and velocity space [@stutzki-et-al-1998; @brunt-heyer-2002; @miville-deschenes2003; @correia-et-al-2016]. 

Dust density {#sec:modeldensity}
------------

![ Total gas column density ${N_\mathrm{H}}$, derived from a synthetic density field ${n_\mathrm{H}}$ built by exponentiation of a fBm field with spectral index $\beta_X=2.6$ and size $120 \times 120 \times 120$ pixels.  The volume density fluctuation level is $y_n=1$, and the column density fluctuation level is $y_{{N_\mathrm{H}}}\simeq0.25$.[]{data-label="fig:densitymaps"}]({Density_map_0_120x120x120_b2.6_Y1.0-eps-converted-to}.pdf){width="9cm"}

In our approach, the dust density ${n_\mathrm{d}}$ is taken to be proportional to the total gas density ${n_\mathrm{H}}$, so that the dust optical depth within each cell is also proportional to ${n_\mathrm{H}}$ (see the derivation of polarization maps in Sec.  \[sec:polarmaps\]).  Therefore, we mean to model ${n_\mathrm{H}}$ from numerical realizations of three-dimensional fBm fields built in Fourier space.  These have means defined by the value chosen for the null-wavevector amplitude $A(\boldsymbol{0})$, so if one wished to use such a synthetic random field $X$ directly as a model for the positive-valued ${n_\mathrm{H}}$, one would be required to choose ${n_\mathrm{H}}=X'=X-a$ with $a\geqslant\mathrm{min}(X)$ a constant.  However, since the distributions of these fields in 3D are close to Gaussian, their ratio of standard deviation to mean is typically $\sigma_{X'}/\langle X'\rangle\lesssim 0.3$, which is much too small compared to observational values.  For instance, the total gas column density fluctuation ratios $\sigma_{{N_\mathrm{H}}}/\left<{N_\mathrm{H}}\right>$ in the ten nearby molecular clouds selected for study in [@planck2014-XX] are in the range 0.3-1, and one should keep in mind that these are only lower bounds for fluctuation ratios in the three-dimensional density field ${n_\mathrm{H}}$.

We remedy this shortcoming by taking $X$ to represent the log-density, i.e., ${n_\mathrm{H}}$ is given by $$\label{eq:expfbm}
{n_\mathrm{H}}=n_0\exp{\left(\frac{X}{X_r}\right)},$$ where $X$ is a three-dimensional fBm field with spectral index $\beta_X$, and $X_r$ and $n_0$ are positive parameters.  The ${n_\mathrm{H}}$ fields built in this fashion have simple and well-controlled statistical properties.  First, their probability distribution functions (PDF) are log-normal, which allows, through an adequate choice of $X_r$, to set the fluctuation level of the density field $y_n=\sigma_{{n_\mathrm{H}}}/\langle{n_\mathrm{H}}\rangle$ to any desired value.  Second, their power spectra, as azimuthal averages in Fourier space, retain the power-law behaviour of the original fBm $X$, $$P_{{n_\mathrm{H}}}(k)=\left<\left|\widetilde{{n_\mathrm{H}}}(\boldsymbol{k})\right|^2\right>_{||\boldsymbol{k}||=k} \propto k^{-\beta_n},$$ although the spectral indices $\beta_n$ may deviate significantly from $\beta_X$. An example of such a field is shown in Fig.  \[fig:densitymaps\], which represents the total gas column density ${N_\mathrm{H}}$ derived from a gas volume density ${n_\mathrm{H}}$ built as the exponential of a $120 \times 120 \times 120$ fractional Brownian motion with zero mean, unit variance, and spectral index $\beta_X=2.6$. The parameters of the exponentiation are $X_r=1.2$ and $n_0=20\,\mathrm{cm}^{-3}$, and the grid is chosen so that the extent of the cube is 30pc on each side, corresponding to a pixel size of 0.25pc.  More details on the properties of these fields are given in Appendix \[sec:appendix:nHproperties\]. 

The density fields built in this fashion are of course only a rough statistical approximation for actual interstellar density fields.  For instance, they are unable to reproduce the filamentary structures observed in dust emission maps [@andre-et-al-2010; @miville-deschenes-et-al-2010].  These structures cannot be captured by one- and two-point statistics such as those used here, and require a description involving higher-order moments, or equivalently of the Fourier phases [see, e.g., @levrier-et-al-2004; @burkhart-lazarian-2016]. 

Magnetic field
--------------

![ Synthetic magnetic field ${\boldsymbol{B}}$ built using Eq.  \[eq:Btlambda\].  The spectral index of the vector potential is $\beta_A=5$ and the size of the cubes is $120 \times 120 \times 120$ pixels, corresponding to 30pc on each side.  Shown are 2D slices through the cubes of the components $B_{x}$ ([*top*]{}), $B_{y}$ ([*middle*]{}), and $B_{z}$ ([*bottom*]{}).  The ratio of the fluctuations of the turbulent component ${\boldsymbol{B}}_t$ to the norm of the uniform component ${\boldsymbol{B}}_0$ is $y_B=1$ in this particular case, with angles $\chi_0=\gamma_0=0^\circ$. []{data-label="fig:Bmaps"}]({Magnetic_map_0_120x120x120_b5.0_Y1.0_c0.0_g0.0_BMap_Slices-eps-converted-to}.pdf){width="9cm"}

To obtain a synthetic turbulent component of the magnetic field ${\boldsymbol{B}}_t$ with null divergence and controlled power spectrum, we start from a vector potential $\boldsymbol{A}$ built as a three-dimensional fractional Brownian motion.  To be more precise, each Cartesian component $A_\lambda$ of $\boldsymbol{A}$ is a fBm field, $$\widetilde{A_\lambda}(\boldsymbol{k})=\mathcal{A}_0k^{-\beta_A/2}\exp\left[i\phi_{A_\lambda}(\boldsymbol{k})\right],$$ where the spectral index $\beta_A$ and the overall normalization parameter $\mathcal{A}_0$ are independent of the Cartesian component $\lambda=x,y,z$ considered.  Using the definition of the magnetic field from the vector potential $B_{t,\lambda}=\epsilon_{\lambda\mu\nu}\partial_\mu A_\nu$, where $\epsilon_{\lambda\mu\nu}$ is the Levi-Civita tensor, and the derivation relation in Fourier space $$\widetilde{\partial_\lambda F}=ik_\lambda \widetilde{F}$$ we have the expression of the components of ${\boldsymbol{B}}_t$ in Fourier space $$\label{eq:Btlambda}
\widetilde{B_{t,\lambda}}(\boldsymbol{k})=\mathcal{A}_0\epsilon_{\lambda\mu\nu}ik_\mu k^{-\beta_A/2}\exp\left[i\phi_{A_\nu}(\boldsymbol{k})\right]$$ As it should, this expression corresponds to a divergence-free turbulent magnetic field, $$ik_\lambda\widetilde{B_{t,\lambda}}=0,$$ using the Einstein notation.  Writing $k_\mu=kf_\mu$, with $\boldsymbol{f}=\left(\sin\vartheta\cos\varphi,\sin\vartheta\sin\varphi,\cos\vartheta\right)$, the power spectrum of each component of ${\boldsymbol{B}}_t$ is then $$P_{B_{t,\lambda}}(k)=\mathcal{A}_0^2k^{2-\beta_A}\left<\left|\epsilon_{\lambda\mu\nu}f_\mu\exp\left[i\phi_{A_\nu}(\boldsymbol{k})\right]\right|^2\right>_{||\boldsymbol{k}||=k}.$$ The last factor is essentially independent of the wavenumber $k$, so the spectral index of each component of ${\boldsymbol{B}}_t$ is $\beta_{B_t}=\beta_A-2$. After Fourier-transforming back to real space, ${\boldsymbol{B}}_t$ is shifted and scaled so that it has zero mean and a standard deviation $\sigma_B$ of $5\,\mu\mathrm{G}$, a value typical of the interstellar magnetic field [see, e.g., @haverkorn-et-al-2008 and references therein]. 

The model magnetic field ${\boldsymbol{B}}$ is obtained by adding a uniform[^4] vector field ${\boldsymbol{B}}_0$ to that turbulent magnetic field ${\boldsymbol{B}}_t$. The effect in Fourier space is limited to a modification for $\boldsymbol{k}=\boldsymbol{0}$ only, so the spectral index of each component $B_\lambda$ of the total magnetic field is the same as that of $B_{t,\lambda}$, i.e., $$\label{eq:PB}
P_{B_{\lambda}}(k)\propto k^{2-\beta_A}.$$ Note that this means that the resulting magnetic fields thus only display anisotropy in the $\boldsymbol{k}=\boldsymbol{0}$ mode, and not at the other scales.  This is a limitation of our model, which is thus not fully consistent with observations of the magnetic field structure [@planck2015-XXXV], but it is sufficient for our purposes. 

The uniform field ${\boldsymbol{B}}_0$ which is added to the turbulent field ${\boldsymbol{B}}_t$ is defined by its norm $B_0$ and a pair of angles, $\gamma_0$ and $\chi_0$, which are respectively the angle between the magnetic field and the POS, and the position angle of the projection of ${\boldsymbol{B}}_0$ in the POS, counted positively clockwise from the north-south direction [see figure 14 of @planck2014-XX].  The total magnetic field’s direction in three-dimensional space is characterized by angles $\gamma$ and $\chi$ defined in the same way.  The ratio of the turbulent to mean magnetic field strengths is then defined by $$y_B=\frac{\sigma_B}{B_0}=\frac{\sqrt{\left<{\boldsymbol{B}}_t^2\right>-\left<{\boldsymbol{B}}_t^{\phantom{2}}\right>^2}}{||{\boldsymbol{B}}_{0}||}=\frac{\sqrt{\left<{\boldsymbol{B}}_t^2\right>}}{||{\boldsymbol{B}}_{0}||}.$$

Fig.  \[fig:Bmaps\] shows an example of a synthetic magnetic field ${\boldsymbol{B}}$ generated in this fashion, and defined on the same 120 $\times$ 120 $\times$ 120 pixels grid that was used for the gas density model described in Sec.  \[sec:modeldensity\].  The parameters used for this specific realization are $\beta_A=5$, $y_B=1$, and $\chi_0=\gamma_0=0^\circ$. The PDFs of the components of such model magnetic fields are Gaussian, as shown in Fig.  \[fig:Bpdfs\], and their power spectra are power laws of the wavenumber, as shown in Fig.  \[fig:Bpowerspectra\], with a common spectral index $\beta_B$ that is related to the input $\beta_A$ by $\beta_B = \beta_A - 2$.

Polarization maps {#sec:polarmaps}
-----------------

Once cubes of total gas density ${n_\mathrm{H}}$ and magnetic field ${\boldsymbol{B}}$ are available, maps of Stokes parameters $I$, $Q$, and $U$ at 353GHz (the frequency of the [[*Planck*]{}]{} channel with the best signal-to-noise ratio in polarized thermal dust emission) are built by integrating along the line of sight through the simulation cubes, following the method in [@planck2014-XX] : $$\begin{aligned}
\label{eq:defI}
I_0&=&\int S_\nu\left[1-p_0\left(\cos^2\gamma-\frac{2}{3}\right)\right]\sigma_{353}\,{n_\mathrm{H}}{\,\mathrm{d}}z;\\
\label{eq:defQ}
Q_0&=&\int p_0\,S_\nu\cos\left(2\phi\right)\cos^2\gamma\,\sigma_{353}\,{n_\mathrm{H}}{\,\mathrm{d}}z;\\
\label{eq:defU}
U_0&=&\int p_0\,S_\nu\sin\left(2\phi\right)\cos^2\gamma\,\sigma_{353}\,{n_\mathrm{H}}{\,\mathrm{d}}z.\end{aligned}$$ In these equations, we take the intrinsic polarization fraction parameter $p_0$ to be uniform, and the source function $S_\nu=B_\nu(T_\mathrm{d})$ to be that of a blackbody with an assumed uniform dust temperature $T_\mathrm{d}$. The dust opacity at this frequency $\sigma_{353}$ is taken to vary with ${N_\mathrm{H}}$, following figure 20 from [@planck2013-p06b] for $X_\mathrm{CO}=2\,10^{20}\,\mathrm{H_2\,cm^{-2}\,K^{-1}\,km^{-1}\,s}$, and propagating the associated errors.  The order of magnitude of the dust opacity is around $\sigma_{353}\approx 10^{-26}\,\mathrm{cm}^{2}$. The values of ${N_\mathrm{H}}$ considered in our study are typically at most a few $10^{21}\,\mathrm{cm}^{-2}$, so the optically thin limit applies in the integrals of Eqs.  \[eq:defI\]-\[eq:defU\].  The angle $\phi$ is the local polarization angle, which is related to the position angle[^5] $\chi$ of the magnetic field’s projection on the POS at each position on the LOS by a rotation of $90^\circ$ [see definitions of angles in @planck2014-XX]. 

The ${n_\mathrm{H}}$ and ${\boldsymbol{B}}$ cubes are built on a grid which is $132\times132$ pixels in the POS and $N_z$ pixels in the $z$ direction (that of the LOS).  The cells have a physical size $\delta=0.24\,\mathrm{pc}$ in each direction, so the depth $d=N_z\,\delta$ of the cloud is a free parameter in our analysis, and the Stokes maps built from Eq.  \[eq:defI\]-\[eq:defU\] are $32\,\mathrm{pc}$ in both $x$ and $y$ directions. 

Noise and beam convolution
--------------------------

In order to proceed with the analysis of observational data, one cannot use these model Stokes maps directly : it is necessary to properly take into account noise and beam convolution.  Anticipating somewhat on the description of the [*Planck*]{} data we shall use as an application of the method, the 353GHz noise covariance matrix maps are taken directly from the [*[*Planck*]{} Legacy Archive*]{}[^6] and are part of the 2015 public release of [[*Planck*]{}]{} data [@planck2014-a01], $$\boldsymbol{\Sigma}=\left(\begin{array}{ccc}
\sigma_{II} & \sigma_{IQ} & \sigma_{IU} \\
\sigma_{QI} & \sigma_{QQ} & \sigma_{QU} \\
\sigma_{UI} & \sigma_{UQ} & \sigma_{UU} \\
\end{array}\right).$$ Noise is added to the model Stokes maps pixel by pixel, as $$\label{eq:addnoise}
I_n=I_0+n_I \qquad Q_n=Q_0+n_Q \qquad U_n=U_0+n_U,$$ where $n_I$, $n_Q$, and $n_U$ are random values drawn from a three-dimensional Gaussian distribution with zero mean and characterized by the noise covariance matrix $\boldsymbol{\Sigma}$. To preserve the properties of [[*Planck*]{}]{} noise, the random values are directly drawn from the [Healpix]{} [@gorski_et_al_2005] covariance matrix maps and added to the simulated maps after a gnomonic projection of the region under study, in our case the Polaris Flare molecular cloud. 

The resulting Stokes $I_n$, $Q_n$, and $U_n$ maps are then placed at a distance[^7] $D=140\,\mathrm{pc}$, so that the angular size of each pixel is about $6\arcmin$, and then convolved by a circular 15 full-width at half maximum (FWHM) Gaussian beam $\mathcal{B}$. To avoid edge effects, only the central $120\times120$ pixels of the convolved maps $I_\mathrm{m}=\mathcal{B}\otimes I_n$, $Q_\mathrm{m}=\mathcal{B}\otimes Q_n$, and $U_\mathrm{m}=\mathcal{B}\otimes U_n$ are retained, corresponding to a field-of-view (FoV) of approximately 12$^\circ$. With this approach, we ensure that these model maps $(I_\mathrm{m},Q_\mathrm{m},U_\mathrm{m})$ are fit to be compared to actual [[*Planck*]{}]{} data, which we discuss in Sec.  \[sec:Polaris\]. 

[cccc]{}\
\[-1ex\] Parameter & Prior & Definition\
\[1ex\]\
\[-1ex\] $\beta_B$ & $\left[1,4\right]$ & Spectral index of the 3D turbulent magnetic field\
$\beta_n$ & $\left[1,5\right]$ & Spectral index of the 3D dust density field\
$\log_{10}{y_n}$ & $\left[-1,1\right]$ & Log of the RMS-to-mean ratio of dust density\
$\log_{10}{y_B^{\rm POS}}$ & $\left[-1,1\right]$ & Log of the ratio of the turbulent magnetic field RMS to the mean magnetic field in the POS\
$\chi_0$ & $\left[-90^\circ,90^\circ\right]$ & Position angle of the mean magnetic field in the POS\
$\log_{10}\left(d/1\,\mathrm{pc}\right)  $ & $\left[-0.3,1.5\right]$ & Depth of the simulated cube\
$\log_{10}\left(\langle{n_\mathrm{H}}\rangle/1\,\mathrm{cm}^{-3}\right)  $ & $\left[1,2.7\right]$ & Mean dust density\
$T_\mathrm{d} $ & $\left[5\,\mathrm{K},200\,\mathrm{K}\right]$ & Dust temperature\
$p_0 $ & $\left[0.01,0.5\right]$ & Intrinsic polarisation fraction parameter\
\[1ex\]\
\[-1ex\]

Observables {#sec:observables}
===========

![ image]({Stokes_map_132x132x40_b2.3_A4.6_Y0.6_YB1.0_c-50.0_g0.0_nH30.0_T18.0_a0.2_Beam15_Reso6_NoisePolaris_I-eps-converted-to}.pdf){width="43.00000%"}

![ image]({Stokes_map_132x132x40_b2.3_A4.6_Y0.6_YB1.0_c-50.0_g0.0_nH30.0_T18.0_a0.2_Beam15_Reso6_NoisePolaris_Q-eps-converted-to}.pdf){width="43.00000%"}

![ image]({Stokes_map_132x132x40_b2.3_A4.6_Y0.6_YB1.0_c-50.0_g0.0_nH30.0_T18.0_a0.2_Beam15_Reso6_NoisePolaris_U-eps-converted-to}.pdf){width="43.00000%"}

From the model maps above, we build an ensemble of derived maps, starting with the normalized Stokes maps $$i=\frac{I_\mathrm{m}}{\langle I_\mathrm{m}\rangle} \qquad q=\frac{Q_\mathrm{m}}{I_\mathrm{m}} \qquad u=\frac{U_\mathrm{m}}{I_\mathrm{m}}$$ where $\langle I_\mathrm{m}\rangle$ is the spatial average of the model Stokes $I_\mathrm{m}$ map.  Then we define the polarization fraction, which requires us to note that since our models include noise, we should not use the “naïve” estimator [@montier-et-al-2015a; @montier-et-al-2015b] $$p=\frac{\sqrt{Q_\mathrm{m}^2+U_\mathrm{m}^2}}{I_\mathrm{m}},$$ but rather the modified asymptotic (MAS) estimator proposed by [@plaszczynski_et_al_14] $$p_\mathrm{MAS}=p-b^2\frac{\displaystyle 1-e^{-p^2/b^2}}{2p}$$ where the noise bias parameter $b^2$ derives from the elements of the noise covariance matrix $\boldsymbol{\Sigma}$ [see @montier-et-al-2015b].  Next, we define the polarization angle $$\psi=\frac{1}{2}\mathrm{atan}{\left(U_\mathrm{m},Q_\mathrm{m}\right)}$$ where the two-argument $\mathrm{atan}$ function lifts the $\pi$-degeneracy of the usual $\mathrm{atan}$ function.  Note that this expression means that the polarization angle is defined in the [Healpix]{} convention. 

We also build maps of the polarization angle dispersion function $\mathcal{S}$ [@planck2014-XIX; @planck2014-XX; @alina-et-al-2016], which quantifies the local dispersion of polarization angles at a given lag $\delta$ and is defined by $$\mathcal{S}(\boldsymbol{r},\delta)=\sqrt{\frac{1}{\mathcal{N}}\sum_{i=1}^\mathcal{N}\left[\psi\left(\boldsymbol{r}+\boldsymbol{\delta}_i\right)-\psi\left(\boldsymbol{r}\right)\right]^2}$$ where the sum is performed over the $\mathcal{N}$ pixels $\boldsymbol{r}+\boldsymbol{\delta}_i$ whose distance to the central pixel $\boldsymbol{r}$ lies between $\delta/2$ and $3\delta/2$. For the sake of consistency with the analysis performed on simulated polarization maps in [@planck2014-XX], we take $\delta=16\arcmin$.

Finally we build the column density and optical depth $\tau_{353}$ maps from the dust density cube using $${N_\mathrm{H}}= \int{n_\mathrm{H}}{\,\mathrm{d}}z \quad\text{ and }\quad \tau_{353} = \sigma_{353}\left({N_\mathrm{H}}\right) \times{N_\mathrm{H}}$$ with the $\sigma_{353}\left({N_\mathrm{H}}\right)$ conversion[^8] from [@planck2013-p06b].  A temperature map $T_{\rm obs}$ is also created using the anti-correlation with the column density ${N_\mathrm{H}}$ observed in the data.  This “dust temperature map” does not pretend to model reality but since $T_\mathrm{d}$ is one of the parameters of the model, the fitting algorithm requires a map whose mean value should yield $T_\mathrm{d}$.

Exploring the parameter space {#sec:method}
=============================

The goal of this paper is to constrain the physical parameters of molecular clouds, in particular the spectral indices of the dust density and of the turbulent magnetic field, using [[*Planck*]{}]{} polarization maps and a grid of model maps built as explained in the previous section. 

Parameter space
---------------

The nine physical parameters that are explored in this paper using fBm simulations are summarized in Table \[tab:priors\].  They are sufficient to describe the one-point and two-point statistical properties of the dust density and magnetic field models.  Note that unlike what was done in [@planck2016-XLIV], the field of view of the maps analysed in the following (approximately $12^\circ$) is too small to contain remarkable features that could be used to constrain the angle $\gamma_0$ that the mean magnetic field makes with the POS.  In such small fields of view, there is a degeneracy between $y_B$ and $\gamma_0$ which cannot be lifted.  Consequently, we chose not to try to fit for $\gamma_0$ and $y_B$, but for the ratio of the turbulent magnetic field RMS to the mean magnetic field in the POS, i.e., $$y_B^{\rm POS} = \frac{\sigma_B}{B_0^{\mathrm{POS}}}=\frac{\sqrt{\left<{\boldsymbol{B}}_t^2\right>}}{||{\boldsymbol{B}}_{0}||\cos\gamma_0} = \frac{y_B}{\cos\gamma_0}.$$

This analysis was applied to the Polaris Flare (see section \[sec:Polaris\]), so the priors are chosen to be flat over a reasonably large range, to cover the expected physical values of the molecular cloud under consideration, but they are necessary for the analysis to converge.  The cloud’s average column density is of the order $\langle{N_\mathrm{H}}\rangle\approx 10^{21}\,\mathrm{cm}^{-2}$ [@planck2014-XX].  This value was used to set the range for the prior on the depth $d$ of the cube, with the limits of this range chosen in such a way that the average gas density $\langle{n_\mathrm{H}}\rangle$ lies between 10 and $500\,\mathrm{cm}^{-3}$, a reasonable assumption for the Polaris Flare molecular cloud.  This translates to a total cube depth $d$ between 0.5 and $32.5\,\mathrm{pc}$. The range used for the prior on $\beta_n$ is justified by a number of observational studies [see, e.g., the review by @hennebelle-falgarone-2012], and that on $\beta_B$ is chosen based on the results from [@vansyngel-et-al-2017], but also on numerical studies of MHD turbulence [see, e.g., @perez-et-al-2012; @beresnyak-2014].  The fluctuation ratios $y_n$ and $y_B^{\rm{POS}}$ are explored in a logarithmic scale, as we are mainly interested in order of magnitude estimates for these parameters.  The polarization maps being statistically identical when the angle $\chi_0$ of the POS projection of the mean magnetic field is shifted by $180^\circ$, the prior on this parameter is such that this periodicity is applied when the Metropolis algorithm (see \[sec:MCMC\]) draws values outside the given range.  The priors chosen for $T_\mathrm{d}$ and $p_0$ are very large and do not play a role in the fitting procedure. 

Comparing models with data {#sec:chi2}
--------------------------

To set constraints on the parameters listed in Table \[tab:priors\],

we build a likelihood function, which expresses the probability that a given set of synthetic polarization maps reproduce adequately actual observational data.  From the model Stokes maps $I_\mathrm{m}$, $Q_\mathrm{m}$, and $U_\mathrm{m}$, we derive a set of observables that are used in the likelihood function.  These observables are given in Table \[tab:observables\].  More precisely, we use i) the mean values for the optical depth $\tau_{353}$ and the dust temperature $T_\mathrm{obs}$, ii) the distribution functions (one-point statistics) of the $I_\mathrm{m}$, $Q_\mathrm{m}$, $U_\mathrm{m}$, $p_\mathrm{MAS}$, $\psi$, $\mathcal{S}$, and $\tau_{353}/\left\langle \tau_{353}\right\rangle$ maps, iii) the power spectra (two-point statistics) of the $I_\mathrm{m}$ ,$Q_\mathrm{m}$, and $U_\mathrm{m}$ maps, and iv) the pixel-by-pixel anti-correlation between $\mathcal{S}$ and $p_\mathrm{MAS}$ underlined by [@planck2014-XIX].  Indeed, we have found that the shape of this two-dimensional distribution function also depends on the model parameters.  Many other observables were tested but we have retained only those which bring constraints on the model parameters. 

[cc]{}\
\[-1ex\] Type & From\
\[1ex\]\
\[-1ex\] Mean values & $\tau_{353}$, $T_\mathrm{obs}$\
Distribution function & $I_\mathrm{m}$, $Q_\mathrm{m}$, $U_\mathrm{m}$, $p_\mathrm{MAS}$, $\psi$, $\mathcal{S}$, $\tau_{353}/\left\langle \tau_{353}\right\rangle$\
Power spectrum & $I_\mathrm{m}$, $Q_\mathrm{m}$, $U_\mathrm{m}$\
Correlation & $\left\lbrace\mathcal{S},p_\mathrm{MAS}\right\rbrace$\
\[1ex\]\
\[-1ex\]

On the simulation side, $N_r=60$ model realizations per set of parameter values are generated with their observables, to be compared with data.  The $N_r$ models differ by the random phases $\phi_X$ and $\phi_{A_\lambda}$ used to build the dust density and magnetic field cubes (see Eqs.  \[eq:Xfbm\] and \[eq:Btlambda\]), and by the random realization of the noise applied to the model (Eq.  \[eq:addnoise\]).  We checked that 60 simulations represent a statistically large enough sample to get robust averages and dispersions for the observables.  The statistical properties of the observables derived from the observational polarization data are thus compared with the observables from those 60 models, through the evaluation of a parameter $D^2$ which quantifies the distance between data and one random realisation of the model, averaged over the $N_r$ random realisations, with contributions associated to the various observables listed in Table \[tab:observables\], i.e., $$\label{eq:chi2tot}
D^2 = \frac{1}{N_r}\sum_{i=1}^{N_{\rm r}} \left[ D^2_{ \mu} + \sum_{o}D^2_{{\rm DF}(o)} + \sum_{o} D^2_{P(o)}  + D^2_{\mathcal{S}-p_\mathrm{MAS}}\right].$$ This quantity is subtly different from the usual $\chi^2$ (see Appendices \[sec:appendix:L2distance\] and \[sec:appendix:chi2\]).  The first term in Eq.  \[eq:chi2tot\] covers the observables $\left\langle \tau_{353}\right\rangle$ and $\left\langle T_\mathrm{obs}\right\rangle$, and quantifies the difference between these values in the simulated maps and in the data.  The second sum extends over the observable maps $o$ in the set $\left\{I_\mathrm{m},Q_\mathrm{m},U_\mathrm{m},p_\mathrm{MAS},\psi,\mathcal{S},\tau_{353}/\left\langle \tau_{353}\right\rangle\right\}$ and quantifies the difference between the distribution functions (DF) of these observables in synthetic maps and those of the same observables in the data.  The third sum extends over the observable maps $o$ in the set $\left\{I_\mathrm{m},Q_\mathrm{m},U_\mathrm{m}\right\}$ and quantifies the difference between the power spectra of the simulated maps and those of the same maps in the data.  Finally, the last term quantifies the discrepancy between the two-dimensional joint DFs of $\mathcal{S}$ and $p_\mathrm{MAS}$ in the data and in synthetic maps.  We detail the computation of these various terms in Appendix \[sec:appendix:L2distance\]. 

MCMC chains {#sec:MCMC}
-----------

Given the vast parameter space to explore, we built a Monte Carlo Markov Chain method (MCMC) [see, e.g. @MCMC] that has the advantage to sample specifically the regions of interest in this space.  We used a simple Metropolis-Hastings algorithm to build five Markov chains which sample the posterior probability distribution of the parameters listed in Table \[tab:priors\].  The likelihood $\mathcal{L}$ of a set $s$ of parameters is evaluated thanks to the $D^2$ criterion described in section \[sec:chi2\] and Appendix \[sec:appendix:L2distance\] as $$\mathcal{L}(s) \propto e^{-D^2(s)/2}\pi(s)$$ with $\pi(s)$ the prior associated to the parameters. 

According to the Metropolis-Hastings algorithm, at each step $q$ of the chain, parameters are drawn according to a multivariate probability distribution function whose covariance is set to allow for an efficient exploration of the parameter space, with an average given by the parameter values $s_{q-1}$ at step $q-1$. If the likelihood for the new set of parameters, $s_q$, is larger than for the previous one, then the chain records the new set.  Otherwise, the likelihood ratio $\mathcal{L}\left(s_q\right)/\mathcal{L}\left(s_{q-1}\right)< 1$ is compared to a number $\alpha$ drawn randomly from a uniform distribution over $[0,1]$. If the likelihood ratio is larger than $\alpha$, the $s_q$ set of parameters is kept, otherwise the chain duplicates the $s_{q-1}$ set, i.e., $s_q=s_{q-1}$. The posterior probability distribution function is then given by the occurence frequency of the parameters along the chains, after removal of the initial “burn-in” phase. 

The priors used for each parameter are detailed in Table \[tab:priors\].  For all parameters, flat priors are set covering a reasonable range of physical interest.  If the Metropolis algorithm draws values outside of these priors we set $\pi(s)=0$, except for the position angle of the mean magnetic field, $\chi_0$, for which the $180^\circ$ periodicity is used to bring back the angle inside its definition range when it is perchance drawn outside. 

The convergence of the Markov chains is tested using the Gelman-Rubin statistics $R$ [@gelman1992], which is essentially the ratio of the variance of the chain means to the mean of the chain variances.  We estimate that the chains converged when $R-1<0.03$ for the least-converged parameter.  The convergence is also assessed by checking visually the $D^2$ and parameter evolutions along the chains. 

![ image]({Polaris_120_27_0_132x132_Beam15_Reso6_I-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Polaris_120_27_0_132x132_Beam15_Reso6_Q-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Polaris_120_27_0_132x132_Beam15_Reso6_U-eps-converted-to}.pdf){width="33.00000%"}\
![ image]({Polaris_120_27_0_132x132_Beam15_Reso6_PoverI-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Polaris_120_27_0_132x132_Beam15_Reso6_Psi-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Polaris_120_27_0_132x132_Beam15_Reso6_DeltaPsi_R16-eps-converted-to}.pdf){width="33.00000%"}

The obtained 9D posterior probability distribution is generally not a multivariate Gaussian distribution.  To quote an estimate of the best fit value for any one of the nine parameters and the associated uncertainties, we first marginalize over the other eight parameters to obtain the one-dimensional posterior PDF for the remaining parameter.  In the following, the quoted best fit value for a parameter is the mean over this posterior PDF (which is less sensitive to binning effect than the maximum likelihood).  As the PDFs are usually not Gaussian, we quote asymmetric error bars following the minimum credible interval technique [see, e.g., @hamann2007]. 

Results {#sec:results}
=======

Validation of the method {#sec:validation}
------------------------

To validate the fitting method, we simulated four sets of model cubes and computed the corresponding $I_\mathrm{m}$, $Q_\mathrm{m}$, and $U_\mathrm{m}$ maps, including noise, with different values of the input parameters (these are labelled simulations A, B, C and D hereafter).  The MCMC fitting procedure was run on these mock polarization data sets to check if it was able to recover the statistical properties of the input dust density and magnetic field cubes through the selected observables.  The results are presented in Table \[tab:results\_sim\].  For the four sets of maps, the fitting method recovered the input values within the quoted uncertainties, after the convergence criteria for all the chains were reached[^9].  This shows that this choice of observables is relevant to extract the input values from polarized thermal dust emission data within our model.  To assess the goodness of fit of the model to the data, we use an [*a posteriori*]{} $\chi^2$ test, as explained in Appendix \[sec:appendix:chi2\]. 

In all four cases, we find that the match is very good, since $\left<\chi^2_{\rm{best}}\right>\approx 1$. For illustration, the posterior probability contours for simulation A are presented in Figure \[fig\_triangle\_simA\].  We note that the MCMC procedure reveals correlations between the model parameters, which is not unexpected, e.g., between $y_B^\mathrm{POS}$ and $p_0$, or between $y_n$ and $\langle{n_\mathrm{H}}\rangle$. These trends are best visualized with the correlation matrix, shown in the upper right corner of Fig.  \[fig\_triangle\_simA\]. 

![ image]({Stokes_map_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_nH40.0_T17.5_a0.12_Beam15_Reso6_NoisePolaris_I-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Stokes_map_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_nH40.0_T17.5_a0.12_Beam15_Reso6_NoisePolaris_Q-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Stokes_map_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_nH40.0_T17.5_a0.12_Beam15_Reso6_NoisePolaris_U-eps-converted-to}.pdf){width="33.00000%"}\
![ image]({Stokes_map_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_nH40.0_T17.5_a0.12_Beam15_Reso6_NoisePolaris_PoverI-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Stokes_map_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_nH40.0_T17.5_a0.12_Beam15_Reso6_NoisePolaris_Psi-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Stokes_map_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_nH40.0_T17.5_a0.12_Beam15_Reso6_NoisePolaris_DeltaPsi_R16-eps-converted-to}.pdf){width="33.00000%"}

![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_HistI_PDF-eps-converted-to}.pdf){width="40.00000%"} ![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_HistQ_PDF-eps-converted-to}.pdf){width="40.00000%"}

![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_HistU_PDF-eps-converted-to}.pdf){width="40.00000%"} ![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_HistPoverI_PDF-eps-converted-to}.pdf){width="40.00000%"}

![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_HistPsi_PDF-eps-converted-to}.pdf){width="40.00000%"} ![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_HistDeltaPsi_PDF-eps-converted-to}.pdf){width="40.00000%"}

![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_HistTau353_PDF-eps-converted-to}.pdf){width="40.00000%"}

![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_PSI-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_PSQ-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_PSU-eps-converted-to}.pdf){width="33.00000%"}\

![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_correlation_data-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_correlation_sim-eps-converted-to}.pdf){width="33.00000%"} ![ image]({Stokes_Maps_132x132x50_b2.3_A4.8_Y1.65_YB0.64_c-69.0_g0.0_correlation_residuals-eps-converted-to}.pdf){width="33.00000%"}

Application to the Polaris Flare {#sec:Polaris}
--------------------------------

As an application of our method, we wish to constrain statistical properties of the turbulent magnetic field in the Polaris Flare, a diffuse, highly dynamical, non-starforming molecular cloud.  There are several reasons for choosing this particular field.  First, it has been widely observed : the structures of matter were studied in dust thermal continuum emission by, e.g., [@miville-deschenes-et-al-2010]; the velocity field of the molecular gas was studied down to very small scales through CO rotational lines [@falgarone-et-al-1998; @hily-blant-et-al-2009]; and the magnetic field was probed by optical stellar polarization data in [@panopoulou-et-al-2016].  Second, as this field does not show signs of star formation, the dynamics of the gas and dust are presumably dominated by magnetized turbulence processes, without contamination by feedback from young stellar objects (YSOs).  It therefore seems like an ideal test case for our method. 

To this aim, we use the full-mission [*Planck*]{} maps of Stokes parameters $(I_{353},Q_{353},U_{353})$ at 353GHz and, as already mentioned, the associated covariance matrices from the [[*Planck*]{}]{} Legacy Archive.  We also use the thermal dust model maps $\tau_{353}$ and $T_{\rm obs}$ from the 2013 public release [@planck2013-p06b].  All maps are at a native $4\arcmin.8$ resolution in the [Healpix]{} format with $N_\mathrm{side}=2048$, and the Polaris Flare maps are obtained by projecting these onto a Cartesian grid with 6$\arcmin$ pixels, centered on Galactic coordinates $(l,b)=(120^\circ,27^\circ)$, with a field of view $\Delta l=\Delta b=12^\circ$. The maps of $I_{353}$, $Q_{353}$, $U_{353}$, and $\tau_{353}$ are then smoothed using a circular Gaussian beam, to obtain maps at a 15$\arcmin$ FWHM resolution.  The covariance matrix maps are computed at the same resolution, using a set of Monte-Carlo simulations of pure noise maps, drawn from the original full-resolution covariance maps and smoothed at 15.  The maps of $I_{353}$, $Q_{353}$, $U_{353}$, $p_\mathrm{MAS}$, $\psi$, and $\mathcal{S}$ obtained in this way are shown in Fig.  \[fig\_polaris\_map\]. 

Note that the features of simulated Stokes maps are not located in the same regions as in the Polaris Flare maps.  As the noise covariance matrices are the same for all the simulated maps, this means that the signal-to-noise ratio per pixel in the model $I_{\rm m}$, $Q_{\rm m}$ and $U_{\rm m}$ maps is different for each set of parameters.  However, the MCMC procedure is able to choose the parameter sets that give signal-to-noise ratios similar to those in the [[*Planck*]{}]{} maps. 

The results of the analysis of the [[*Planck*]{}]{} polarized thermal dust emission data towards the Polaris Flare are presented in Table \[tab:results\], and the posterior probability distribution contours are shown in Figure \[fig\_triangle\].  We find in particular that the spectral index of the turbulent component of the magnetic field is $\beta_B=2.8\pm 0.2$, and that the spectral index of the dust density field is around $\beta_n=1.7$ with a rather large uncertainty.  The fluctuation ratio of the density field is about unity, $y_n\approx 1.6$, and the magnitude of the large scale magnetic field in the POS dominates slightly the RMS of the turbulent component, $y_B^{\rm POS}\approx 0.65$, with a position angle $\chi_0 \approx -69^\circ$. The constraint on the depth of the cloud seems to indicate that $d\approx 13\,\mathrm{pc}$, with $\langle{n_\mathrm{H}}\rangle\approx 40\,\mathrm{cm}^{-3}$. The temperature $T_\mathrm{d}$ is 17.5K equal to the average of the $T_{\rm obs}$ [[*Planck*]{}]{} map, and the polarization fraction is $p_0\approx 0.12$. The parameter set for simulation A was chosen a posteriori to give similar best fit parameters and to test our likelihood method in the conditions driven by the Polaris Flare data. 

Using the best fitting parameters from Table \[tab:results\] we performed simulations to visually check the agreement between the model and [[*Planck*]{}]{} data.  Figure \[fig\_best\_fit\_map\] shows the polarization maps from a simulation using these best fitting parameters.  The overall similarity with the data maps from Figure \[fig\_polaris\_map\] is reasonably good, although spatially coherent structures appear in the data maps which cannot be reproduced by the model maps.  The agreement between the best fitting simulation and the data is quantified through plots of the different observables that were used in the fitting procedure (Figs.  \[fig\_best\_fit\_pdf\], \[fig\_best\_fit\_ps\] and \[fig\_best\_fit\_corr\]).  The agreement is excellent for most observables, although substantial deviations are visible in the DFs of the intensity $I_{\rm m}$, of the normalized optical depth $\tau_{353}/\left\langle\tau_{353}\right\rangle$ and of the polarization angle $\psi$. These deviations are due to the simplifying assumptions of our fBm model.  It may be that in the Polaris Flare the large scale magnetic field has two major components that superimposed, with global orientations $\chi_0 \approx -70^\circ$ and $\chi_0 \approx 50^\circ$. Note that the DF in Fig.  \[fig\_best\_fit\_pdf\] is that of the $\psi$ angle, which differs from $\chi_0$ by 90$^\circ$. Also the exponentiation procedure to model the dust field is a good but incomplete approximation of the reality and it is not able to totally reproduce the shapes of the $I_{\rm m}$ and $\tau_{353}/\left\langle\tau_{353}\right\rangle$ DFs together.  These deviations impact the reduced best fit $\left<\chi^2_{\rm{best}}\right> \approx 2.9$, which is somewhat larger than for mock data ($\approx 1$), but still reasonably good. 

Concerning the mean values used as observables, the Polaris Flare has a mean optical depth of $\left\langle\tau_{353}\right\rangle = \left(1.25 \pm 0.05 \right)\times 10^{-5}$ and a mean temperature of $\left\langle T_{\rm obs} \right\rangle = 17.5 \pm 0.4$K. Using the best fitting parameters from Table \[tab:results\] we get optical depth maps with an average of $\left\langle\tau_{353}\right\rangle = \left(1.82 \pm 0.05 \right)\times 10^{-5}$ over 60 realizations which is in tension with the data value, as mentioned above for the $\tau_{353}/\left\langle\tau_{353}\right\rangle$ discrepancy.  However the best fitting parameter for temperature is $T_\mathrm{d} = 17.5\pm 0.5$K which is exactly the same as in data with a width reflecting the data uncertainties. 

The observables we used to extract the statistical properties of the Polaris Flare field are by themselves unable to constrain the $\gamma_0$ angle of the large scale magnetic field on the LOS.  However, the [@planck2016-XLIV] analysis was able to fit the $\chi_0$ and $\gamma_0$ angle is the southern Galactic cap and found an intrinsic polarization fraction of the gas of $p_\mathrm{int}\approx 0.26$. If we believe this latter value is true also in the Polaris Flare, then it is related to our fit as $p_0 \approx p_\mathrm{int} \left\langle\cos^2 \gamma \right\rangle \approx p_\mathrm{int} \cos^2 \gamma_0$. We can thus constrain the $\gamma_0$ angle to be around $45^\circ$.

Discussion and summary {#sec:discussion}
======================

We have presented an analysis framework for maps of polarized thermal dust emission in the diffuse ISM aimed at constraining the statistical properties of the dust density and magnetic field responsible for this emission.  Our framework rests on a set of synthetic models for the dust density and magnetic field, for which we precisely control the one- and two-point statistics, and on a least-squares analysis in which the space of parameters is explored via a MCMC method.  The application of the method to [[*Planck*]{}]{} maps of the Polaris Flare molecular cloud leads to a spectral index of the turbulent component of the magnetic field $\beta_B=2.8\pm 0.2$, which is in very good agreement with the findings of [@planck2016-XLIV] and [@vansyngel-et-al-2017], who used a very different approach over a much larger fraction of the sky.  The dust density field exhibits a much flatter spectrum, $\beta_n=1.7$. This latter exponent is remarkably close to the Kolmogorov index for the velocity field in incompressible hydrodynamical turbulence, but this comparison should be taken with caution, as closer examination of the power spectrum of the model density field shows a spectral break with an exponent closer to 2.2 at the largest scales ($k\lesssim 1\,\mathrm{pc}^{-1}$) while the smaller scales ($k\gtrsim 1\,\mathrm{pc}^{-1}$) have a 1.7 exponent[^10].  What is clear is that the magnetic field power spectrum is much steeper, which underlines the role that the large scale magnetic field plays in the structure of polarized emission maps.  We find that the fluctuation ratio of the dust density field and the ratio of turbulent-to-uniform magnetic field are both around unity.  Finally, our analysis is able to give a constraint on the polarization fraction, $p_0 \approx 0.12$, and on the depth of the Polaris Flare molecular cloud, $d\approx 13\,\mathrm{pc}$, which is about half the transverse extent of the field-of-view, with $\langle{n_\mathrm{H}}\rangle \approx 40 \,\mathrm{cm}^{-3}$. The good visual agreement between the Polaris Flare maps and model maps for the best-fitting parameters (Figs.  \[fig\_polaris\_map\] and \[fig\_best\_fit\_map\]), and the excellent agreement between the two sets of maps for most of the observables used in the analysis (Figs.  \[fig\_best\_fit\_pdf\], \[fig\_best\_fit\_ps\] and \[fig\_best\_fit\_corr\]), all lead us to conclude that our fBm-based model, although limited, provides a reasonable description of the magnetized, turbulent, diffuse ISM. 

In fact, it is quite remarkable to find such a good agreement with the data, considering the limitations of the model.  First, it is statistically isotropic, and therefore cannot reproduce the interstellar filamentary structures observed at many scales and over a large range in column densities [see, e.g. @miville-deschenes-et-al-2010; @arzoumanian-et-al-2011].  Second, our model dust density and magnetic fields are completely uncorrelated, which is clearly not realistic, as it was found that there is a preferential relative orientation between structures of matter and magnetic field, both in molecular clouds [@planck2015-XXXV] and in the diffuse, high-latitude sky [@planck2014-XXXII; @planck2015-XXXVIII].  The change in relative orientation, from mostly parallel to mostly perpendicular, as the total gas column density ${N_\mathrm{H}}$ increases, is also not reproducible with our fully-synthetic models.  Third, it is now commonly acknowledged that two-point statistics such as power spectra are not sufficient to properly describe the structure of interstellar matter.  Improving our synthetic models along these three directions will be the subject of future work. 

For completeness, we have also looked into applying our MCMC approach based on fBm models to synthetic polarization maps built from a numerical simulation of MHD turbulence.  We used simulation cubes from [http://www.mhdturbulence.com]{} [@cho-lazarian-2003; @burkhart-et-al-2009; @burkhart-et-al-2014], basing our choice on the simulation parameters, which seemed more or less consistent with the parameters found for the Polaris Flare data.  We built simulated Stokes $I$, $Q$, and $U$ maps using the same resolution and noise parameters, and launched the MCMC analysis on these simulated Stokes maps.  It turns out that the Markov chains have a much harder time converging than when applying the method to the [[*Planck*]{}]{} data.  It is not yet completely clear why that is so, but we suspect that part of the reason may lie with the limited range of spatial scales over which the fields in the MHD simulation can be accurately described by scale-invariant processes.  Indeed, while the fBm models exhibit power-law power spectra over the full range of accessible scales (basically one decade in our case), the MHD simulations are hampered by effects of numerical dissipation at small scales (possibly over nearly 10 pixels), and the properties at large scales are dependent on the forcing, which is user-defined.  The data, on the other hand, exhibit a much larger « inertial range ».  In that respect, our fBm models, despite all their drawbacks, and despite the fact that they lack the physically realistic content of MHD simulations, provide a better framework for assessing the statistical properties of the [[*Planck*]{}]{} data than current MHD simulations can.  Of course, this conclusion is based on just one simulation, and there would definitely be a point in applying the MCMC approach to assess various instances of MHD simulations with respect to the observational data, based on the same observables, but independently of the grid of fBm models.  This project, however, is clearly beyond the scope of this paper. 

We gratefully acknowledge fruitful discussions with S. Plaszczynski and O. Perdereau. 

Statistical properties of ${n_\mathrm{H}}$ models {#sec:appendix:nHproperties}
=================================================

Probability distribution function
---------------------------------

The PDF $f({n_\mathrm{H}})$ of the density field ${n_\mathrm{H}}$ built using Eq.  \[eq:expfbm\] derives from the Gaussian PDF of $X$, which in all generality has a mean $\langle X\rangle$ and variance $\sigma_X^2$. We thus have a log-normal PDF $$\label{eq:PDFn}
f({n_\mathrm{H}})=\frac{1}{\sqrt{2\pi}\sigma_X}\frac{X_r}{{n_\mathrm{H}}}\exp{\left\{-\frac{1}{2\sigma_X^2}\left[X_r\ln{\left(\frac{{n_\mathrm{H}}}{n_0}\right)}-\langle X\rangle\right]^2\right\}},$$ which is defined for ${n_\mathrm{H}}>0$. Figure \[fig:densitypdfs\] presents the distribution function of the ${n_\mathrm{H}}$ field[^11] used to build Fig.  \[fig:densitymaps\], with the theoretical PDF expected from Eq.  \[eq:PDFn\].  It should be noted that the distribution function for a single realisation over a finite grid such as the ones used here may deviate from the theoretical PDF, especially for large values of $\beta_X$, but the mean distribution function over a sufficiently large sample converges to the lognormal form (Eq.  \[eq:PDFn\]). 

Moments and fluctuation level
-----------------------------

From Eq.  \[eq:PDFn\], we may compute moments of any order $p$ of the PDF of the total gas density ${n_\mathrm{H}}$ $$\left<{n_\mathrm{H}}^p\right>=\int\nolimits_{0}^{\infty} {n_\mathrm{H}}^p f({n_\mathrm{H}}){\,\mathrm{d}}{n_\mathrm{H}}=n_0^p\exp{\left(p\frac{\left<X\right>}{X_r}+p^2\frac{\sigma_X^2}{2X_r^2}\right)},$$ which allows, in particular, to compute its mean value $$\left<{n_\mathrm{H}}\right>=\int\nolimits_{0}^{\infty} {n_\mathrm{H}}f({n_\mathrm{H}}){\,\mathrm{d}}{n_\mathrm{H}}=n_0\exp{\left(\frac{\left<X\right>}{X_r}+\frac{\sigma_X^2}{2X_r^2}\right)},$$ as well as its variance $$\sigma_{{n_\mathrm{H}}}^2=n_0^2\exp{\left(2\frac{\left<X\right>}{X_r}\right)}\left[\exp{\left(\frac{2\sigma_X^2}{X_r^2}\right)}-\exp{\left(\frac{\sigma_X^2}{X_r^2}\right)}\right].$$ The density fluctuation level, which is one of the parameters of our model, is therefore $$y_n=\frac{\sigma_{{n_\mathrm{H}}}}{\langle{n_\mathrm{H}}\rangle}=\sqrt{2}\exp{\left(\frac{\sigma_X^2}{4X_r^2}\right)}\left[\sinh{\left(\frac{\sigma_X^2}{2X_r^2}\right)}\right]^{1/2}.$$ For instance, the ${n_\mathrm{H}}$ field whose distribution function is shown in Fig.  \[fig:densitypdfs\] has values ranging from 0.3$\mathrm{cm}^{-3}$ to $880\,\mathrm{cm}^{-3}$, with a mean and standard deviation of $\langle{n_\mathrm{H}}\rangle=\sigma_{{n_\mathrm{H}}}=28.3\,\mathrm{cm}^{-3}$, resulting in the desired fluctuation level $y_n=1$.

![ Distribution function of the synthetic field ${n_\mathrm{H}}$ used for Fig.  \[fig:densitymaps\] (black histogram).  The red curve shows the theoretical probability distribution function $f({n_\mathrm{H}})$ for the chosen set of parameters.[]{data-label="fig:densitypdfs"}]({nH_PDF_beta2.6_N120_Xr1.2_Yr20.0_single-eps-converted-to}.pdf){width="9.5cm"}

Power spectra
-------------

Figure \[fig\_11\] shows the azimuthally-averaged power spectra of ${n_\mathrm{H}}$ fields obtained through Eq.  \[eq:expfbm\], from a $120 \times 120 \times 120$ pixels fractional Brownian motion $X$ with spectral index $\beta_X=3$, for various fluctuation levels $y_n$. The spectra were normalized differently so as to allow comparison between them. 

![ Evolution of the differences $\beta_n-\beta_X$ between the spectral indices of the original fBm field $X$ and that of the model density field ${n_\mathrm{H}}$ with fluctuation level $y_n$. Each point corresponds to the mean of 20 realisations of the model density field ${n_\mathrm{H}}$, and the error bars represent the standard deviation of the fitted spectral indices $\beta_n$ and fluctuation levels $y_n$. []{data-label="fig:beta_nH_with_yn"}]({beta_nH_with_yn-eps-converted-to}.pdf){width="9.5cm"}

The power-law behaviour is apparent, even at large fluctuation levels, but the spectral index decreases (i.e., the spectra flatten) as the fluctuation level increases.  This is quantified in Fig.  \[fig:beta\_nH\_with\_yn\], which shows the differences $\beta_n-\beta_X$ between the spectral indices of the original fBm field $X$ and that of the model density field ${n_\mathrm{H}}$. The power spectra of the latter are indeed flatter than the original ones ($\beta_n<\beta_X$), with differences that may become large when $\beta_X$ is low, but remain negligible for higher $\beta_X$. We interpret this trend as the exponentiation process amplifying the two-point differences $X(\boldsymbol{r}+\boldsymbol{\delta})-X(\boldsymbol{r})$, for small separations $||\boldsymbol{\delta}||$, that exist in the original field when $\beta_X$ is low, leading to an increase of the small-scale power in the model ${n_\mathrm{H}}$ field, i.e. $\beta_n<\beta_X$. This effect is all the more important than exponentiation stretches these field differences more strongly, i.e. when $y_n$ increases.  On the other hand, for high $\beta_X$, the fBm fields are much smoother, so the exponentiation process has little impact on these two-point statistics at small scales, leading to $\beta_n\simeq \beta_X$. Note that at low fluctuation levels ($y_n\leqslant 0.3$), the differences are smaller than 0.1 for all values of $\beta_X$.

Likelihood terms {#sec:appendix:L2distance}
================

$D^2$ terms for mean values {#sec:chi2mean}
---------------------------

The first term in Eq.  \[eq:chi2tot\] is given by $$D^2_\mu =\frac{\left(\left\langle \tau_{353}^{\rm m} \right\rangle - \left\langle \tau_{353}^{\rm d} \right\rangle \right)^2}{\sigma_{\left\langle \tau_{353}^{\rm m} \right\rangle}^2 + {\sigma_{\left\langle \tau_{353}^{\rm d} \right\rangle}^2}}+\frac{\left( \left\langle T_{\rm obs}\right\rangle - T_\mathrm{d} \right)^2}{\sigma_{\left\langle T_{\rm obs}\right\rangle}^2}.$$

The mean value of the optical depth $\tau_{353}$ is evaluated on the simulated $ \tau_{353}^{\rm m}$ and [[*Planck*]{}]{} data $ \tau_{353}^{\rm d} $ maps.  These are compared through a standard $\chi^2$ test.  The denominator includes the uncertainty $\sigma_{\left\langle \tau_{353}^{\rm d} \right\rangle}$ on the mean $\tau^{\rm d}_{353}$ value propagated from the uncertainty map provided by the [[*Planck*]{}]{} collaboration [@planck2013-p06b], and the uncertainty $\sigma_{\left\langle \tau_{353}^{\rm m} \right\rangle}$ coming from the conversion factor $\sigma_{353}\left({N_\mathrm{H}}\right)$ used to build the simulated map [@planck2013-p06b]. 

In the second term, $\left\langle T_{\rm obs}\right\rangle$ is the mean value of the temperature map $T_{\rm obs}$ from [[*Planck*]{}]{} data, and $\sigma_{\left\langle T_{\rm obs}\right\rangle}$ represents the uncertainty on the averaged value propagated from the uncertainty map provided by the [[*Planck*]{}]{} collaboration.  $T_\mathrm{d}$ is directly the model parameter for the dust temperature. 

$D^2$ terms for distribution functions {#sec:chi2DF}
--------------------------------------

For a given observable map $o$, we compute its DF over an ensemble of $N_{\rm b}$ bins[^12].  When considering the [[*Planck*]{}]{} data, we write this DF as $h_{o,i}^\mathrm{d}$, where $i$ is the bin number, and we estimate the uncertainty on the value of the DF in bin $i$ through the associated Poisson noise $\sigma_{h_{o,i}^\mathrm{d}}$. When considering the model, we write $h_{o,i}^\mathrm{m}$ and $\sigma_{h_{o,i}^\mathrm{m}}$ to be respectively the bin value and the Poisson noise of the DF in bin $i$, independently for each of the $N_r=60$ model realizations. 

The contribution $D^2_{{\rm DF}(o)}$ of the observable’s DF to the total $D^2$ in Eq.  \[eq:chi2tot\] is then built as an average over the $N_{\rm b}$ bins $$\label{eq:D2DF}
D^2_{{\rm DF}(o)} = \frac{1}{N_{\rm b}}  \sum_{i=1}^{N_{\rm b}}  \frac{\left(h_{o,i}^\mathrm{d} - h_{o,i}^\mathrm{m} \right)^2}{(\sigma_{h_{o,i}^\mathrm{d}})^2 + (\sigma_{h_{o,i}^\mathrm{m}})^2}$$ and averaged over the $N_{\rm r}=60$ realizations.  The inner sum is normalized to the number of bins so that $D^2_{{\rm DF}(o)}$ is less sensitive to the binning choice and the map noise.  If the model fits the data correctly as far as the DF of observable $o$ is concerned, then $D^2_{{\rm DF}(o)}$ is minimum. 

This quantity is different from a standard $\chi^2$ test as it compares data with one random[^13] realisation for a given set of model parameters, and this comparison is repeated and averaged $N_r$ times.  A standard $\chi^2$ test would compare data with a model prediction that would be the average of the random realizations (see Eq \[eq:chi2DF\] in Appendix \[sec:appendix:chi2\]).  The latter could not be used in our MCMC analysis due to the mathematical relation between the power spectrum of the map and the variance in each bin of the DF : with steep power spectra (high values of the spectral index $\beta$), only the few large-scale modes effectively contribute to the power, leading to a large variance in each DF bin[^14].  Thus the $\chi^2$ test tends to favor these large values of $\beta$, as they yield large denominators and thus allow for a “good” fit.  This drives the fit towards a region of parameter space yielding mean DFs that fit the data well but with a huge dispersion: one realization of such a model gives a DF with bin values highly scattered even though the data DF is quite smooth.  This means that the data cannot reasonably be interpreted as a random realization using these parameter values, and explains why we had to switch to the $D^2$ function, which directly compares data with one model random realization.  In this fashion, we are able to reach the region of parameter space correctly describing the data. 

$D^2$ terms for power spectra
-----------------------------

For a given observable map $o$, we first compute its two-dimensional power spectrum $$P_o(\boldsymbol{k})=\left|\,\widetilde{o}(\boldsymbol{k})\right|^2$$ then average these within $N'_{\rm b}$ annuli in Fourier space, centred on a set of wavenumbers $\left\{k_i\right\}$. We write $P_{o,i}^\mathrm{d}$ for this azimuthal average in bin $i$ when considering the [[*Planck*]{}]{} data, and $P_{o,i}^\mathrm{m}$ when considering model fields.  The uncertainties affecting these quantities strongly depend on the number $N_i$ of wavevectors $\boldsymbol{k}_n$ in each bin.  The best estimate for the standard deviation $\sigma_{P_{o,i}^\mathrm{d}}$ of the power spectrum of the data is $$\sigma_{P_{o,i}^\mathrm{d}} = t_{N_i}\sqrt{\frac{1}{N_i\left(N_i-1\right)}\sum_{n=1}^{N_i} \left[ P_{o}^{d} (\boldsymbol{k}_n) - P_{o,i}^d \right]^2}$$ where the factor $t_N$ is the Student coefficient.  We thus obtain the best estimate of the true standard deviation in bins with only a few modes (i.e., at large scale).  The standard deviation $\sigma_{P_{o,i}^\mathrm{m}}$ for the 60 model realizations is computed in the same way as for data. 

The contribution $D^2_{P(o)}$ of the observable’s power spectrum to the total $D^2$ in Eq.  \[eq:chi2tot\] is then computed as a sum over the $N'_{\rm b}$ bins[^15] in wavenumber space, $$\label{eq:chi2PS}
D^2_{P(o)} =\frac{1}{N'_{\rm b}} \sum_{i=1}^{N'_{\rm b}} \frac{\left(P_{o,i}^\mathrm{d} - P_{o,i}^\mathrm{m} \right)^2}{(\sigma_{P_{o,i}^\mathrm{d}})^2 + (\sigma_{P_{o,i}^\mathrm{m}})^2}.$$

$D^2$ term for the $\left\lbrace\mathcal{S},p_\mathrm{MAS}\right\rbrace$ anti-correlation
-----------------------------------------------------------------------------------------

To use the $\left\lbrace\mathcal{S},p_\mathrm{MAS}\right\rbrace$ anti-correlation [@planck2014-XIX; @planck2014-XX], we compute the joint distribution function of the $\mathcal{S}$ and $p_\mathrm{MAS}$ maps, which we write $h_{ij}^\mathrm{d}$ and $h_{ij}^\mathrm{m}$ for the [[*Planck*]{}]{} data and model maps respectively, with $1\leqslant i\leqslant N_{{\rm b},1}$ and $1\leqslant j\leqslant N_{{\rm b},2}$ the binning scheme used for the two maps.  The standard deviations $\sigma_{ij}^\mathrm{d}$ and $\sigma_{ij}^\mathrm{m}$ are defined in the same way as for the one-dimensional DFs in \[sec:chi2DF\], and the contribution $D^2_{\mathcal{S}-p_\mathrm{MAS}}$ to the total $D^2$ is then $$D^2_{\mathcal{S}-p_\mathrm{MAS}} = \frac{1}{N_{\rm b,tot}} \sum_{i=1}^{N_{\rm b,1}} \sum_{j=1}^{N_{\rm b,2}} \frac{\left(h_{ij}^\mathrm{d} - h_{ij}^\mathrm{m} \right)^2}{(\sigma_{ij}^\mathrm{d})^2 + (\sigma_{ij}^\mathrm{m})^2}.$$ In this expression, it should be noted that the total number $N_{\rm b,tot}$ of two-dimensional bins considered is less than the product $N_{\rm b,1}N_{\rm b,2}$ of the number of bins in each dimension, which we set to $N_{\rm b,1}=N_{\rm b,2}=50$. The reason for this is that we discard the empty bins and those with a signal-to-noise ratio below three[^16].  We thus keep only the significantly populated bins that can drive the fit and contribute to the total $D^2$.

Goodness-of-fit {#sec:appendix:chi2}
===============

To assess the goodness of the fit, we use an [*a posteriori*]{} $\chi^2$ test, which we define as $$\label{eq:chi2tot_gof}
\chi^2 = \frac{1}{N_o}\left[ \chi^2_{ \mu} + \sum_{o}\chi^2_{{\rm DF}(o)} + \sum_{o} \chi^2_{P(o)}  + \chi^2_{\mathcal{S}-p_\mathrm{MAS}}\right]$$ with $N_o=13$ the total number of observables.  Each term is a $\chi^2$ test comparing data with the mean of the $N_r=60$ realisations.  The first term from Eq \[eq:chi2tot\_gof\] is $$\chi^2_\mu = \frac{\left(\overline{\left\langle \tau_{353}^{\rm m} \right\rangle} - \left\langle \tau_{353}^{\rm d} \right\rangle \right)^2}{\Sigma_{\left\langle \tau_{353}^{\rm m} \right\rangle}^2 + {\sigma_{\left\langle \tau_{353}^{\rm d} \right\rangle}^2}} +\frac{\left( \left\langle T_{\rm obs}\right\rangle - T_\mathrm{d} \right)^2}{\sigma_{\left\langle T_{\rm obs}\right\rangle}^2}$$ where $\overline{\left\langle \tau_{353}^{\rm m} \right\rangle} $ is the ensemble average over the 60 model realisations of the (spatial) mean of the optical depth.  In the following, the brackets stand for an average on the pixels while the upper bar represents the average over the $N_r$ realisations.  The $\Sigma_{\left\langle \tau_{353}^{\rm m} \right\rangle}^2$ quantity is the variance of $\left\langle \tau_{353}^{\rm m} \right\rangle$ over the $N_r$ realisations (capital $\Sigma$ denotes the variance over the random realisations). 

The second term is $$\label{eq:chi2DF}
\chi^2_{\rm{DF}(o)} = \frac{1}{N_{\rm b}} \sum_{i=1}^{N_{\rm b}} \frac{\left(h_{o,i}^\mathrm{d} - \overline{ h_{o,i}^\mathrm{m}} \right)^2}{(\sigma_{h_{o,i}}^\mathrm{d})^2 + (\Sigma_{h_{o,i}}^\mathrm{m})^2}$$ where $$\overline{ h_{o,i}^\mathrm{m}}=\frac{1}{N_r}\sum_{k=1}^{N_r}h_{o,i}^\mathrm{m}$$ is the ensemble average of the $i^\mathrm{th}$ bin of the DF for the observable $o$, over the $N_r=60$ model realizations, and $(\Sigma_{ h_{o,i}}^\mathrm{m})^2$ is the associated variance.  Note that while $D^2_{\rm DF(o)}$ in Eq.  \[eq:D2DF\] is the average of the observable $\chi^2$, $\chi^2_{\rm{DF}(o)}$ is the $\chi^2$ of the averaged observable. 

The third term is $$\chi^2_{P(o)} = \frac{1}{N'_{\rm b}} \sum_{i=1}^{N'_{\rm b}} \frac{\left(P_{o,i}^\mathrm{d} - \overline{P_{o,i}^\mathrm{m}} \right)^2}{(\sigma_{P_{o,i}^\mathrm{d}})^2 + (\Sigma_{{P_{o,i}^\mathrm{m}}})^2}$$ where $\overline{P_{o,i}^\mathrm{m}}$ is the averaged power spectrum in bin $i$ and $(\Sigma_{{P_{o,i}^\mathrm{m}}})^2$ its variance. 

Finally, the fourth term is $$\chi^2_{\mathcal{S}-p_\mathrm{MAS}} = \frac{1}{N_{\rm b,tot}} \sum_{i=1}^{N_{\rm b,1}} \sum_{j=1}^{N_{\rm b,2}} \frac{\left(h_{ij}^\mathrm{d} - \overline{h_{ij}^\mathrm{m}} \right)^2}{(\Sigma_{ij}^\mathrm{d})^2 + (\sigma_{ij}^\mathrm{m})^2}$$ with the same notation conventions as above. 

To quantify the goodness of fit, once the MCMC procedure has converged, we perform 100 fits for the set of best fitting parameters, each of these fits comprising 60 model realizations and providing a value of the $\chi^2$ quantity defined in equation \[eq:chi2tot\_gof\].  The average of these 100 $\chi^2$ values[^17] is listed as $\left<\chi^2_{\rm{best}}\right>$ in Tables \[tab:results\_sim\] and \[tab:results\]. 

[^1]: [[*Planck*]{}]{} (<http://www.esa.int/Planck>) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA). 

[^2]: In all generality, several spectral indices may be defined, as one may consider the power spectrum of any one of the three cartesian components of ${\boldsymbol{B}}_t$, or of the modulus $|{\boldsymbol{B}}_t|$. Assuming that ${\boldsymbol{B}}_t$ is isotropic, which we will, all of these spectral indices are identical. 

[^3]: In all of this paper, for any field $F$ the notation $\widetilde{F}$ represents its Fourier transform. 

[^4]: Note that we do not consider an ordered random or striated random component of the field [@jaffe-et-al-2010; @jansson-farrar-2012], which we justify by the smallness of the field-of-view considered. 

[^5]: Not to be confused with the corresponding position angle $\chi_0$ of the uniform component of the magnetic field ${\boldsymbol{B}}_0$.

[^6]: http://pla.esac.esa.int/pla/

[^7]: A more recent determination of the distance to the Polaris Flare places it at 350-400pc [@schlafly2014].  For the demonstration of the method presented here, this is not a critical issue, as the power-law power spectra underline self-similar behaviours, so that a change of the distance can be compensated by a change in pixel size. 

[^8]: The conversion factor is given for a map of ${N_\mathrm{H}}$ at a resolution of $30\arcmin$. Thus, before applying it pixel by pixel, we smooth the simulated ${N_\mathrm{H}}$ map to $30\arcmin$ resolution, apply the conversion, then resample the resulting $\tau_{353}$ map at the original resolution. 

[^9]: And after removing the burn-in phase, which is quite short in our case ($\lesssim 30\%$ of the chain lengths in general). 

[^10]: Incidentally, from the [*Planck*]{} maps, we can measure the spectral index of the total intensity for the Polaris Flare to be $\beta_I=2.84\pm0.10$, in excellent agreement with the measurement by [@stutzki-et-al-1998] on CO integrated emission at a similar angular resolution. 

[^11]: In that case, $\langle X\rangle=0$, $\sigma_X^2=1$, $X_r=1.2$, and $n_0=20\,\mathrm{cm}^{-3}$.

[^12]: We set $N_{\rm b}=50$ for the $120 \times 120$ pixel maps. 

[^13]: We recall that the randomness comes from the fBm itself and the noise addition to the Stokes maps. 

[^14]: This is akin to the cosmic variance problem in cosmology. 

[^15]: To have reliable and smooth power spectra with $120\times 120$ pixel maps, we set initially $N'_{\rm b}=100$ but later cut off wavenumbers larger than $k_{\rm{max}}=2\pi/(3\times 15\arcmin)$, which corresponds to scales smaller than 3 beam sizes.  Indeed, for bins $k_i > k_{\rm{max}}$ the power spectrum is completely washed out by the beam convolution (see Fig.  \[fig\_maps\_and\_ps\_simA\]) and contains no information about the underlying interesting parameters.  This uninformative part is thus removed from $D^2_{P(o)}$, and thus $N'_{\rm b}<100$.

[^16]: i.e., bins where $h_{ij}^\mathrm{d}=0$, $h_{ij}^\mathrm{m}=0$, $h_{ij}^\mathrm{d}/\sigma_{ij}^\mathrm{d}<3$, or $h_{ij}^\mathrm{m}/\sigma_{ij}^\mathrm{m}<3$.

[^17]: Each simulation is a random realization of the fBm and the noise, then the $\chi^2$ values have an intrinsic dispersion even when compiling 60 simulations.  To check the agreement with data, we thus average 100 computations of the best $\chi^2$ value.