Statistics of the polarized submillimetre emission maps from thermal dust in the turbulent, magnetized, diffuse ISM
^{1}
Sorbonne Université, Observatoire de Paris, Université PSL,
École normale supérieure, CNRS, LERMA,
75005
Paris, France
email: francois.levrier@ens.fr
^{2}
Université ParisSud, LAL, UMR 8607, 91898 Orsay Cedex, France & CNRS/IN2P3,
91405
Orsay, France
^{3}
Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay,
Bât. 121,
91405
Orsay cedex, France
^{4}
Laboratoire AIM, IRFU/Service d’Astrophysique, CEA/DSMCNRS  Université Paris Diderot,
Bât. 709, CEASaclay,
91191
GifsurYvette Cedex, France
^{5}
Nordita, KTH Royal Institute of Technology and Stockholm University,
Roslagstullsbacken 23,
10691
Stockholm, Sweden
^{6}
School of Physical Sciences, National Institute of Science Education and Research, HBNI,
Jatni
752050,
Odisha, India
^{7}
Instituto de Astrofísica de Canarias,
38200
La Laguna,
Santa Cruz de Tenerife, Spain
Received:
19
October
2017
Accepted:
23
February
2018
Context. The interstellar medium (ISM) is now widely acknowledged to display features ascribable to magnetized turbulence. With the public release of Planck data and the current balloonborne and groundbased 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. Aims. We aim to constrain these properties in a statistical way, focussing in particular on the power spectral index β_{B} of the turbulent component of the interstellar magnetic field in a diffuse molecular cloud, the Polaris Flare.
Methods. We present an analysis framework based on simulating polarized thermal dust emission maps using model dust density (proportional to gas density n_{H}) and magnetic field cubes, integrated along the line of sight (LOS), and comparing these statistically to actual data. The model fields are derived from fractional Brownian motion (fBm) processes, which allows a precise control of their one and twopoint statistics. The parameters controlling the model are (1)–(2) the spectral indices of the density and magnetic field cubes, (3)–(4) the RMStomean ratios for both fields, (5) the mean gas density, (6) the orientation of the mean magnetic field in the plane of the sky (POS), (7) the dust temperature, (8) the dust polarization fraction, and (9) the depth of the simulated cubes. We explore the ninedimensional parameter space through a Markov chain Monte Carlo analysis, which yields bestfitting parameters and associated uncertainties.
Results. We find that the power spectrum of the turbulent component of the magnetic field in the Polaris Flare molecular cloud scales with wavenumber as k^{−βB} with a spectral index β_{B} = 2.8 ± 0.2. It complements a uniform field whose norm in the POS is approximately twice the norm of the fluctuations of the turbulent component, and whose position angle with respect to the northsouth direction is χ_{0} ≈−69°. The density field n_{H} is well represented by a lognormally distributed field with a mean gas density 〈n_{H}〉≈40 cm^{−3}, a fluctuation ratio σ_{nH}/〈_{nH}〉≈1.6, and a power spectrum with an index β_{n}=1.7_{−0.3}^{+0.4}. We also constrain the depth of the cloud to be d ≈ 13 pc, and the polarization fraction p_{0} ≈ 0.12. The agreement between the Planck data and the simulated maps for these bestfitting parameters is quantified by a χ^{2} value that is only slightly larger than unity.
Conclusions. We conclude that our fBmbased model is a reasonable description of the diffuse, turbulent, magnetized ISM in the Polaris Flare molecular cloud, and that our analysis framework is able to yield quantitative estimates of the statistical properties of the dust density and magnetic field in this cloud.
Key words: ISM: magnetic fields / ISM: structure / ISM: individual objects: Polaris Flare / polarization / turbulence
© ESO 2018
1 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; WardThompson 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 fullsky 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% (Planck Collaboration Int. XIX 2015), confirming results previously obtained over one fifth of the sky by the Archeops balloonborne experiment (Benoît et al. 2004; Ponthieu et al. 2005). Furthermore, the polarization fraction is anticorrelated with the local dispersion of polarization angles ψ (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XX 2015), and the decrease of the maximum observed p with increasing gas column density N_{H} may be fully accounted for, at the scales probed by Planck (5′ at 353 GHz), by the tangling of the magnetic field on the line of sight (LOS; Planck Collaboration Int. XX 2015). Similar anticorrelations were found with 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 B and filamentary structures of dust emission is consistent with simulated observations derived from numerical simulations of sub or transAlfvénic magnetohydrodynamical (MHD) turbulence (Planck Collaboration Int. XXXV 2016), 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 (Planck Collaboration Int. XXXII 2016; Planck Collaboration Int. XXXVIII 2016) is linked to the measured asymmetry between the power spectral amplitudes of the socalled E and Bmodes 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 (Planck Collaboration Int. XXX 2016).
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) B. Let us write this as a sum B = B_{0} + B_{t} of a uniform, largescale component B_{0}, and a turbulent component B_{t} that has a null spatial average, ⟨B_{t}⟩ = 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} = σ_{B}∕B_{0}, where and B_{0} = B_{0}, and ii) the spectral index β_{B}, which characterizes the distribution of power of B_{t} across spatial scales, through the relationship P(k) ∝ k^{−βB}, where k is the wavenumber and P(k) is the power spectrum ^{2}.
As mentioned above, Planck Collaboration Int. XXXV (2016) 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 find that this relative orientation changes, from mostly parallel to mostly perpendicular, as the total gas column density N_{H} increases, which is a trend observed in simulations of transAlfvénic or subAlfvénic MHD turbulence (Soler et al. 2013). Using the DavisChandraskeharFermi method (Chandrasekhar & Fermi 1953) improved by FalcetaGonçalves et al. (2008) and Hildebrand et al. (2009), and from their results, we estimated the ratio y_{B} to be in the range 0.3–0.7. Planck Collaboration Int. XXXII (2016) 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 Planck Collaboration Int. XLIV (2016), 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 largescale field B_{0} and a small number (N_{l} ≃ 7) of independent “polarization layers” on the LOS, each of which accounts for a fraction 1∕N_{l} of the total unpolarized emission. Within each layer, the turbulent component B_{t} of the magnetic field, which is used to compute synthetic Stokes Q and U maps, is an independent realization of a Gaussian twodimensional random field with a prescribed spectral index β_{B}. Through these fits, they confirm the near equipartition of largescale and turbulent components of B, with y_{B} ≃0.9. They also provide a rough estimate of the magnetic field’s spectral index β_{B} in the range two to three. 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 are able to constrain β_{B} ≃ 2.5, an exponent which is compatible with the rough estimate of Planck Collaboration Int. XLIV (2016), and 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 β_{B}.
In Planck Collaboration Int. XXXII (2016), Planck Collaboration Int. XLIV (2016) 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 β_{B}. Orthogonal approaches have also been pursued (e.g. MivilleDeschênes et al. 2008; O’Dea 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 onedimensional Gaussian random field with a powerlaw 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, ratherthan spatial correlations in the POS, to constrain statistical properties of the interstellar magnetic field.
We sought to explore another avenue, taking into account statistical correlation properties of B in all three dimensions, as well as properties of the dust density field, building on methods developed in Planck Collaboration Int. XX (2015) to compare Planck data with synthetic polarization maps. In that paper, the synthetic maps are computed from data cubes of dust density n_{d} and magnetic field B produced bynumerical simulations of MHD turbulence. One could imagine generalizing this approach, taking advantage of the everincreasing 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 powerlaw 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, threedimensional models for the dust density n_{d} and the magnetic field B, allowing us to perfectly control the statistical properties of these threedimensional 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 leastsquare analysis methods, to extract bestfitting parameters, in particular the spectral index of the magnetic field, β_{B}, and the ratio of turbulent to regular field, y_{B}.
The paper is organized as follows: Sect. 2 presents the method used to build simulated thermal dust polarized emission maps using prescribed statistical properties for n_{d} and B. Observables derived from these maps, serving as statistical diagnostics of the input parameters, are presented in Sect. 3. In Sect. 4, we describe the Markov chain Monte Carlo (MCMC) 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 Sect. 5. Finally, Sect. 6 discusses our results and offers conclusions. Several appendices complement our work. Appendix A presents further statistical properties of the model dust density fields. Appendix B details the likelihood used in the MCMC analysis. Finally, Appendix C details the χ^{2} parameter used to estimate the goodness of fit.
2 Building synthetic polarized emission maps
In this section, we first describe the synthetic dust density and magnetic field cubes we have used in our analysis. We then explain how simulated polarized emission maps are built from these cubes.
2.1 Fractional Brownian motions
The basic ingredients to synthesize polarized thermal dust emission maps are threedimensional cubes of dust density n_{d} and magnetic field B, which we build using fractional Brownian motions (fBm) (Falconer 1990). An Ndimensional fBm X is a random field defined on such that , for any pair of points (r_{1}, r_{2}). H is called the Hurst exponent. These fBm fields are usually built in Fourier space^{3}, (1)
by specifying amplitudes that scale as a powerlaw of the wavenumber k = k, (2)
with β_{X} = 2H + N the spectral index, and phases drawn from a uniform random distribution in [−π, π], subject to the constraint ϕ_{X}(−k) = −ϕ_{X}(k) so that X is realvalued. Their power spectra are therefore power laws, (3)
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; MivilleDeschênes et al. 2003; Correia et al. 2016).
2.2 Dust density
In our approach, the dust density n_{d} is taken to be proportional to the total gas density n_{H}, so that the dust optical depth within each cell is also proportional to n_{H} (see the derivation of polarization maps in Sect. 2.4). Therefore, we intended to model n_{H} from numerical realizations of threedimensional fBm fields built in Fourier space. These have means that are defined by the value chosen for the nullwavevector amplitude A(0), so if one wished to use such a synthetic random field X directly as a model for the positivevalued n_{H}, one would be required to choose n_{H} = X′ = X − a with a ≥min(X) a constant. However, since the distributions of these fields in three dimensions are close to Gaussian, their ratio of standard deviation to mean is typically , which is much too small compared to observational values. For instance, the total gas column density fluctuation ratios in the ten nearby molecular clouds selected for study in Planck Collaboration Int. XX (2015) are in the range 0.3–1, and one should keep in mind that these are only lower bounds for fluctuation ratios in the threedimensional density field n_{H}.
We remedied this shortcoming by taking X to represent the logdensity, i.e., n_{H} is given by (4)
where X is a threedimensional fBm field with spectral index β_{X}, and X_{r} and n_{0} are positive parameters. The n_{H} fields built in this fashion have simple and wellcontrolled statistical properties. First, their probability distribution functions (PDF) are lognormal, which allows, through an adequate choice of X_{r}, to set the fluctuation level of the density field to any desired value. Second, their power spectra, as azimuthal averages in Fourier space, retain the powerlaw behaviour of the original fBm X, (5)
although the spectral indices β_{n} may deviate significantly from β_{X}. An example of such a field is shown in Fig. 1, which represents the total gas column density N_{H} derived from a gas volume density n_{H} built as the exponential of a 120 × 120 × 120 fractional Brownian motion with zero mean, unit variance, and spectral index β_{X} = 2.6. The parameters of the exponentiation are X_{r} = 1.2 and n_{0} = 20 cm^{−3}, and the gridis chosen so that the extent of the cube is 30 pc on each side, corresponding to a pixel size of 0.25 pc. More detail on the properties of these fields are given in Appendix A.
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 (André et al. 2010; MivilleDeschênes et al. 2010). These structures cannot be captured by one and twopoint statistics such as those used here, and require a description involving higherorder moments or the Fourier phases (see, e.g., Levrier et al. 2006; Burkhart & Lazarian 2016).
Fig. 1 Total gas column density N_{H}, derived froma synthetic density field n_{H} built by exponentiation of a fBm field with spectral index β_{X} = 2.6 and size 120 × 120 × 120 pixels. The volume density fluctuation level is y_{n} = 1, and the column density fluctuation level is . 

Open with DEXTER 
2.3 Magnetic field
To obtain a synthetic turbulent component of the magnetic field B_{t} with null divergence and controlled power spectrum, we started from a vector potential A built as a threedimensional fractional Brownian motion. To be more precise, each Cartesian component A_{λ} of A is a fBm field, (6)
where the spectral index β_{A} and the overall normalization parameter are independent of the Cartesian component λ = x, y, z considered. Using the definition of the magnetic field from the vector potential B_{t,λ} = ϵ_{λμν}∂_{μ}A_{ν}, with the Einstein notation, where ϵ_{λμν} is the LeviCivita tensor, and the derivation relation in Fourier space (7)
we have the expression of the components of B_{t} in Fourier space (8)
As it should, this expression corresponds to a divergencefree turbulent magnetic field, (9)
Writing k_{μ} = kf_{μ}, with , the power spectrum of each component of B_{t} is then (10)
The last factor is essentially independent of the wavenumber k, so the spectral index of each component of B_{t} is . After Fouriertransforming back to real space, B_{t} is shifted and scaled so that it has zero mean and a standard deviation σ_{B} of 5 μG, a value typical of the interstellar magnetic field (see, e.g., Haverkorn et al. 2008, and references therein).
The model magnetic field B is obtained by adding a uniform^{4} vector field B_{0} to that turbulent magnetic field B_{t}. The effect in Fourier space is limited to a modification for k = 0 only, so the spectral index of each component B_{λ} of the total magnetic field is the same as that of B_{t,λ}, i.e., (11)
This means that the resulting magnetic fields thus only display anisotropy in the k =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 (Planck Collaboration Int. XXXV 2016), but it is sufficient for our purposes.
The uniform field B_{0} which is added to the turbulent field B_{t} is defined byits norm B_{0} and a pair of angles, γ_{0} and χ_{0}, which are respectively the angle between the magnetic field and the POS, and the position angle of the projection of B_{0} in the POS, counted positively clockwise from the northsouth direction (see Fig. 14 of Planck Collaboration Int. XX 2015). The total magnetic field’s direction in threedimensional space is characterized by angles γ and χ defined in the same way. The ratio of the turbulent to mean magnetic field strengths is then defined by (12)
Figure 2 shows an example of a synthetic magnetic field B generated in this fashion, and defined on the same 120 × 120 × 120 pixels grid that was used for the gas density model described in Sect. 2.2. The parameters used for this specific realization are β_{A} = 5, y_{B} = 1, and χ_{0} = γ_{0} = 0°. The distribution functions (DF) of the components of such model magnetic fields are Gaussian, as shown in Fig. 3, and their power spectra are power laws of the wavenumber, as shown in Fig. 4, with a common spectral index β_{B} that is related to the input β_{A} by β_{B} = β_{A} − 2.
Fig. 2 Synthetic magnetic field B built using Eq. (8). The spectral index of the vector potential is β_{A} = 5 and the size of the cubes is 120 × 120 × 120 pixels, corresponding to 30 pc on each side. Shown are twodimensional 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 B_{t} to the norm of the uniform component B_{0} is y_{B} =1 in this particular case, with angles χ_{0} = γ_{0} = 0°. 

Open with DEXTER 
Fig. 3 Distribution functions of the components B_{x}, B_{y}, and B_{z} of a model magnetic field B = B_{0} + B_{t} built on a grid 120 × 120 × 120 pixels using Eq. (8) with β_{A} = 5, and a mean, largescale magnetic field B_{0} defined by the angles χ_{0} = 0° and γ_{0} = 60°, and a norm B_{0} = 50 μG such that the fluctuation level is y_{B} = 0.1. The vertical lines represent the projected values of the large scale magnetic field B_{0x} = B_{0} sinχ_{0} cosγ_{0}, B_{0y} = −B_{0} cosχ_{0} cosγ_{0} and B_{0z} = B_{0} sinγ_{0}. See Fig. 14 in Planck Collaboration Int. XX (2015) for the definition of angles. 

Open with DEXTER 
Fig. 4 Powerspectra of the components B_{x}, B_{y} and B_{z} of a model magnetic field B = B_{0} + B_{t} built on a grid 120 × 120 × 120 pixels using Eq. (8) with β_{A} = 3 (different shades of blue for the three components) and β_{A} = 5 (different shades of red for the three components). The power spectra are normalized differently so as to allow comparison between them. The fitted powerlaws shown as solid lines yield spectral indices β_{B} = 1 and β_{B} = 3, in agreement with Eq. (11). These are the power spectra of the same particular realizations shown in Fig. 2. 

Open with DEXTER 
Fig. 5 Left column row: example maps (from top to bottom: total intensity I_{m} on a logarithmic scale, Stokes Q_{m}, and Stokes U_{m}) from simulation A (see Sect. 5 and Table 3). Right column: corresponding power spectra. The grey points represent the twodimensional power spectra, while the black dots represent the azimuthal averages in Fourier space in a set of wavenumber bins, and the blue line is a powerlaw fit to the black points. 

Open with DEXTER 
2.4 Polarization maps
Once cubes of total gas density n_{H} and magnetic field B are available, maps of Stokes parameters I, Q, and U at 353 GHz (the frequency of the Planck channel with the best signaltonoise ratio in polarized thermal dust emission) are built by integrating along the LOS through the simulation cubes, following the method in Planck Collaboration Int. XX (2015):
In these equations, we take the intrinsic polarization fraction parameter p_{0} to be uniform, and the source function S_{ν} = B_{ν}(T_{d}) to be that of a black body with an assumed uniform dust temperature T_{d}. The dust opacity at this frequency σ_{353} is taken to vary with N_{H}, following Fig. 20 in Planck Collaboration XI (2014) for X_{CO} = 2 10^{20} H_{2} cm^{2} K^{1} km^{1} s, and propagating the associated errors. The order of magnitude of the dust opacity is around σ_{353} ≈ 10^{−26} cm^{2}. The values of N_{H} considered in our study are typically at most a few 10^{21} cm^{−2}, so the optically thin limit applies in the integrals of Eqs. (13)–(15). The angle ϕ is the local polarization angle, which is related to the position angle^{5} χ of the magnetic field’s projection on the POS at each position on the LOS by a rotation of 90° (see definitions of angles in Planck Collaboration Int. XX 2015).
The n_{H} and B cubes are built on a grid which is 132 × 132 pixels in the POS and N_{z} pixels in the z direction (that of the LOS). The cells have a physical size δ = 0.24 pc in each direction, so the depth d = N_{z} δ of the cloud is a free parameter in our analysis, and the Stokes maps built from Eqs. (13)–(15) are 32 pc in both x and y directions.
2.5 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 the description of the Planck data we shall use as an application of the method, the 353 GHz noise covariance matrix maps are taken directly from the Planck Legacy Archive^{6} and are part of the 2015 public release of Planck data (Planck Collaboration I 2016), (16)
Noise is added to the model Stokes maps pixel by pixel, as (17)
where n_{I}, n_{Q}, and n_{U} are random values drawn from a threedimensional Gaussian distribution with zero mean and characterized by the noise covariance matrix Σ. To preserve the properties of Planck noise, the random values are directly drawn from the Healpix (Górski 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 pc, so that the angular size of each pixel is about 6′, and then convolved by a circular 15′ fullwidth at half maximum (FWHM) Gaussian beam . To avoid edge effects, only the central 120 × 120 pixels of the convolved maps , , and are retained, corresponding to a field of view (FoV) of approximately 12°. With this approach, the model maps (I_{m}, Q_{m}, U_{m}) are fit to be compared to actual Planck data, which we discuss in Sect. 5.2.
3 Observables
From the model maps above, we build an ensemble of derived maps, starting with the normalized Stokes maps (18)
where ⟨I_{m}⟩ is the spatial average of the model Stokes I_{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,b) (19)
but rather the modified asymptotic (MAS) estimator proposed by Plaszczynski et al. (2014) (20)
where the noise bias parameter b^{2} derives from the elements of the noise covariance matrix Σ (see Montier et al. 2015b). Next, we define the polarization angle (21)
where the twoargument atan function lifts the πdegeneracy of the usual atan function. We 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 (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XX 2015; Alina et al. 2016), which quantifies the local dispersion of polarization angles at a given lag δ and is defined by (22)
where the sum is performed over the pixels r + δ_{i} whose distance to the central pixel r lies between δ∕2 and 3δ∕2. For the sake of consistency with the analysis performed on simulated polarization maps in Planck Collaboration Int. XX (2015), we take δ = 16′.
Finally we build the column density and optical depth τ_{353} maps from the dust density cube using (23)
with the conversion^{8} from Planck Collaboration XI (2014). A temperature map T_{obs} is also created using the anticorrelation with the column density N_{H} observed in the data. This “dust temperature map” does not pretend to model reality but since T_{d} is one of the parameters of the model, the fitting algorithm requires a map whose mean value should yield T_{d}.
4 Exploring the parameter space
The goal of this paper is to constrain the physical parameters of molecular clouds. In particular we aim to constrain 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 inthe previous section.
4.1 Parameter space
The nine physical parameters that are explored in this paper using fBm simulations are summarized in Table 1. They are sufficient to describe the onepoint and twopoint statistical properties of the dust density and magnetic field models. We note that contrary to the method in Planck Collaboration Int. XLIV (2016), the FoV of the maps analysed in the following (approximately 12°) is too small to contain remarkable features that could be used to constrain the angle γ_{0} that the mean magnetic field makes with the POS. In such small FoVs, there is a degeneracy between y_{B} and γ_{0} which cannot be lifted. Consequently, we chose not to try to fit for γ_{0} and y_{B}, but for the ratio of the turbulent magnetic field RMS to the mean magnetic field in the POS, i.e., (24)
This analysis was applied to the Polaris Flare (see Sect 5.2), so the priors were 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 ⟨N_{H}⟩≈ 10^{21} cm^{−2} (Planck Collaboration Int. XX 2015). 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 ⟨n_{H} ⟩ lies between 10 and 500 cm^{−3}, a reasonable assumption for the Polaris Flare molecular cloud. This translates to a total cube depth d of between 0.5 and 32.5 pc. The range used for the prior on β_{n} is justified by a number of observational studies (see, e.g., the review by Hennebelle & Falgarone 2012), and that on β_{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 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 χ_{0} of the POS projection of the mean magnetic field is shifted by 180°, the prior on this parameter is such that this periodicity is applied when the Metropolis algorithm (see Sect. 4.3) draws values outside the given range. The priors chosen for T_{d} and p_{0} are very large and do not play a role in the fitting procedure.
Parameter space explored in the grid of model polarization maps.
4.2 Comparing models with data
To set constraints on the parameters listed in Table 1, we built 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_{m}, Q_{m}, and U_{m}, we derived a set of observables that are used in the likelihood function. These observables are given in Table 2. More precisely, we used i) the mean values for the optical depth τ_{353} and the dust temperature T_{obs}, ii) the distribution functions (onepoint statistics) of the I_{m}, Q_{m}, U_{m}, p_{MAS}, ψ, , and maps, iii) the power spectra (twopoint statistics) of the I_{m},Q_{m}, and U_{m} maps, and iv) the pixelbypixel anticorrelation between and p_{MAS} underlined by Planck Collaboration Int. XIX (2015). Indeed, we find that the shape of this twodimensional 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.
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 ϕ_{X} and used to build the dust density and magnetic field cubes (see Eqs. (1) and (8)), and by the random realization of the noise applied to the model (Eq. (17)). 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 2, i.e., (25)
This quantity is subtly different from the usual χ^{2} (see Appendices B and C). The first term in Eq. (25) covers the observables and , 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 and quantifies the difference between the distribution functions 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 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 twodimensional joint DFs of and p_{MAS} in the data and in synthetic maps. We detail the computation of these various terms in Appendix B.
Observables from polarization maps used to fit data.
4.3 MCMC analysis
Given the vast parameter space to explore, we built a MCMC analysis (see, e.g. Brooks et al. 2011) that has the advantage to sample specifically the regions of interest in this space. We used a simple MetropolisHastings algorithm to build five Markov chains which sample the posterior probability distribution function (PDF) of the parameters listed in Table 1. The likelihood of a set s of parameters is evaluated thanks to the D^{2} criterion described in Sect. 4.2 and Appendix B as (26)
with π(s) the prior associated to the parameters.
According to the MetropolisHastings 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 is compared to a number α drawn randomly from a uniform distribution over [0, 1]. If the likelihood ratio is larger than α, 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 “burnin” phase.
The priors used for each parameter are detailed in Table 1. 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 π(s) = 0, except for the position angle of the mean magnetic field, χ_{0}, for which the 180° periodicity is used to bring back the angle inside its definition range when it is by chance drawn outside.
The convergence of the Markov chains is tested using the GelmanRubin statistics R (Gelman & Rubin 1992), which is essentially the ratio of the variance of the chain means to the mean of the chain variances. We estimated that the chains converged when R − 1 < 0.03 for the leastconverged parameter. The convergence is also assessed by visually checking the D^{2} and parameter evolutions along the chains.
The obtained ninedimensional 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 marginalize over the other eight parameters to obtain the onedimensional 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 effects 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., Hamann et al. 2007).
5 Results
5.1 Validation of the method
To validate the fitting method, we simulated four sets of model cubes and computed the corresponding I_{m}, Q_{m}, and U_{m} maps, including noise, with different values of the input parameters (hereafter simulations A, B, C, and D). 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 3. 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 χ^{2} test, as explained in Appendix C. In all four cases, we find that the match is very good, since . For illustration, the posterior probability contours for simulation A are presented in Fig. 6. We note that the MCMC procedure reveals correlations between the model parameters, which is not unexpected, for example, between and p_{0}, or between y_{n} and ⟨n_{H}⟩. These trends are best visualized with the correlation matrix, shown in the upper right corner of Fig. 6.
Fig. 6 Constraints (posterior probability contours and marginalized PDFs) on the statistical properties of the dust density and magnetic field for the simulation A maps. On the posterior probability contours, the filled dark and light blue regions respectively enclose 68.3% and 95.4% of the probability, the black stars indicate the averages over the twodimensional posterior PDFs, and the red circles indicate the input values for the simulation. In the plots showing the marginalized posterior PDFs, the light blue regions enclose 68.3% of the probability, the dashed blue lines indicate the averages over the posterior PDFs, and the solid red lines indicate the input values for the simulation. The upper right plot displays the correlation matrix between the fitted parameters. 

Open with DEXTER 
Fig. 7 Planck 353 GHz maps of the Polaris Flare molecular cloud. The top row shows, from left to right, the total intensity I_{353} on a logarithmic scale, the Stokes Q_{353} map and the Stokes U_{353} map, while the bottom row shows the polarization fraction p_{MAS}, the polarization angle ψ, and the polarization angle dispersion function . The τ_{353} and T_{obs} maps have the same aspects as the I_{353} map but with their own scales. 

Open with DEXTER 
5.2 Application to the Polaris Flare
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, nonstarforming 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, for example, MivilleDeschênes 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; HilyBlant & Falgarone 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 thisaim, we used the fullmission Planck maps of Stokes parameters (I_{353}, Q_{353}, U_{353}) at 353 GHz and, as already mentioned, the associated covariance matrices from the Planck Legacy Archive. We also used the thermal dust model maps τ_{353} and T_{obs} from the 2013 public release (Planck Collaboration XI 2014). All maps are at a native 4′.8 resolution in the Healpix format with N_{side} = 2048, and the Polaris Flare maps were obtained by projecting these onto a Cartesian grid with 6′ pixels, centred on Galactic coordinates (l, b) = (120°, 27°), with an FoVΔl = Δb = 12°. The maps ofI_{353}, Q_{353}, U_{353}, and τ_{353} were then smoothed using a circular Gaussian beam, to obtain maps at a 15′ FWHM resolution. The covariance matrix maps were computed at the same resolution, using a set of MonteCarlo simulations of pure noise maps, drawn from the original fullresolution covariance maps and smoothed at 15′. The maps of I_{353}, Q_{353}, U_{353}, p_{MAS}, ψ, and obtained in this way are shown in Fig. 7.
We notethat 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 signaltonoise ratio per pixel in the model I_{m}, Q_{m} and U_{m} maps is different for each set of parameters. However, the MCMC procedure is able to choose the parameter sets that give signaltonoise 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 4, and the posterior probability distribution contours are shown in Fig. 8. We find in particular that the spectral index of the turbulent component of the magnetic field is β_{B} = 2.8 ± 0.2, and that the spectral index of the dust density field is around β_{n} = 1.7 with a rather large uncertainty. The fluctuation ratio of the density field is about unity, y_{n} ≈ 1.6, and the magnitude of the largescale magnetic field in the POS dominates slightly the RMS of the turbulent component, , with a position angle χ_{0} ≈−69°. The constraint on the depth of the cloud seems to indicate that d ≈ 13 pc, with ⟨n_{H}⟩≈ 40 cm^{−3}. The temperature T_{d} is 17.5 K equal to the average of the T_{obs} Planck map, and the polarization fraction is p_{0} ≈ 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 4 we performed simulations to visually check the agreement between the model and Planck data. Figure 9 shows the polarization maps from a simulation using these best fitting parameters. The overall similarity with the data maps from Fig. 7 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. 10–12). The agreement is excellent for most observables, although substantial deviations are visible in the DFs of the intensity I_{m}, of the normalized optical depth and of the polarization angle ψ. These deviations are due to the simplifying assumptions of our fBm model. It may be that in the Polaris Flare the largescale magnetic field has two major components, with global orientations χ_{0} ≈−70° and χ_{0} ≈ 50°. We note that the DF in Fig. 10 is that of the ψ angle, whichdiffers from χ_{0} by 90°. 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_{m} and DFs together. These deviations impact the reduced bestfit , which is somewhat larger than for mock data (≈1), but still reasonably good.
Concerning the mean values used as observables, the Polaris Flare has a mean optical depth of and a mean temperatureof K. Using the best fitting parameters from Table 4 we obtain optical depth maps with an average of over 60 realizations which is not fully consistent with the data value, as mentioned above for the discrepancy. However the best fitting parameter for temperature is T_{d} = 17.5 ± 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 γ_{0} angle of the largescale magnetic field on the LOS. However, the Planck Collaboration Int. XLIV (2016) analysis was able to fit the χ_{0} and γ_{0} angle is the southern Galactic cap and found an intrinsic polarization fraction of the gas of p_{int} ≈ 0.26. If we believe this latter value is true also in the Polaris Flare, then it is related to our fit as . We can thus constrain the γ_{0} angle to be around 45°.
Fig. 8 Constraints (posterior probability contours and marginalized PDFs) on the statistical properties of the dust density and magnetic field for the Planck maps of the Polaris Flare. On the posterior probability contours, the filled dark and light blue regions respectively enclose 68.3% and 95.4% of the probability, and the black stars indicate the averages over the twodimensional posterior PDFs. In the plots showing the marginalized posterior PDFs, the light blue regions enclose 68.3% of the probability, and the dashed blue lines indicate the averages over the posterior PDFs. The upper right plot displays the correlation matrix between the fitted parameters. 

Open with DEXTER 
Fig. 9 Same as Fig. 7 with the same colour scales, but for model maps using the best fitting parameters to the Polaris Flare data. 

Open with DEXTER 
Fig. 10 Comparison of the DFs extracted from the Planck Polaris Flare maps (black points) with the observables computed from simulations using the best fitting parameters (blue curves). The latter curves are averaged over 60 realizations, as describedin Sect. 4.2: the average is given by the central blue curve and the shaded bands give the 1σ and 2σ standard deviation in each bin. 

Open with DEXTER 
Fig. 11 Comparison of the I_{m}, Q_{m}, and U_{m} power spectra extracted from the Planck Polaris Flare maps (grey points representing the twodimensional power spectra, and black dots representing the azimuthal averages in Fourier space) with the observables computed from simulations using the best fitting parameters (blue curves). The latter curves are averaged over 60 realizations as described in Sect 4.2: the average is given by the central blue curve and the narrow shaded bands give the 1σ and 2σ standard deviation in each bin. 

Open with DEXTER 
Fig. 12 Twodimensional distribution function of and polarization fraction p_{MAS} for the Polaris Flare maps (left), for the model maps using the best fitting parameters averaged over 60 realizations (middle), and residuals (right). The polarization angle dispersion function is computed at a lag δ = 16′. The solid black line shows the mean for each bin in p_{MAS} and the dashed black line is a linear fit of that curve, restricted to bins in p_{MAS} which contain at least 1% of the total number of points (120 × 120). 

Open with DEXTER 
6 Discussion and summary
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 twopoint statistics, and on a leastsquares 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 β_{B} = 2.8 ± 0.2, which is in very good agreement with the findings of Planck Collaboration Int. XLIV (2016) 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, β_{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 ≲ 1 pc^{−1}) while the smaller scales (k ≳ 1 pc^{−1}) have a 1.7 exponent^{10}. It is, however, clear 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 turbulenttouniform magnetic field are both around unity. Finally, our analysis is able to give a constraint on the polarization fraction, p_{0} ≈ 0.12, and on the depth of the Polaris Flare molecular cloud, d ≈ 13 pc, which is about half the transverse extent of the FoV, with ⟨n_{H}⟩≈ 40 cm^{−3}. The good visual agreement between the Polaris Flare maps and model maps for the bestfitting parameters (Figs. 7and 9), and the excellent agreement between the two sets of maps for most of the observables used in the analysis (Figs. 10–12), all lead us to conclude that our fBmbased 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. MivilleDeschênes 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 (Planck Collaboration Int. XXXV 2016) and in the diffuse, highlatitude sky (Planck Collaboration Int. XXXII 2016; Planck Collaboration Int. XXXVIII 2016). The change in relative orientation, from mostly parallel to mostly perpendicular, as the total gas column density N_{H} increases, is also not reproducible with our fullysynthetic models. Third, it is now commonly acknowledged that twopoint 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 polarizationmaps built from a numerical simulation of MHD turbulence. We used simulation cubes^{11} (Cho & Lazarian 2003; Burkhart et al. 2009, 2014), basing our choice on the simulation parameters, which seemed more or lessconsistent 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 it is more difficult for the Markov chains to converge in this case 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 scaleinvariant processes. Indeed, while the fBm models exhibit powerlaw 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 ten pixels), and the properties at large scales are dependent on the forcing, which is userdefined. The data, on theother 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.
Acknowledgements
We gratefully acknowledge fruitful discussions with S. Plaszczynski and O. Perdereau. This work was supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP cofunded by CEA and CNES.
Appendix A Statistical properties of n_{H} models
A.1 Probability distribution function
The PDF f(n_{H}) of the density field n_{H} built using Eq. (4) derives from the Gaussian PDF of X, which in all generality has a mean ⟨X⟩ and variance . We thus have a lognormal PDF (A.1)
which is defined for n_{H} > 0. Figure A.1 presents the distribution function of the n_{H} field^{12} used to build Fig. 1, with the theoretical PDF expected from Eq. (A.1). We note 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 β_{X}, but the mean distribution function over a sufficiently large sample converges to the lognormal form (Eq. (A.1)).
A.2 Moments and fluctuation level
From Eq. (A.1), we may compute moments of any order p of the PDF of the total gas density n_{H} (A.2)
which allows, in particular, to compute its mean value (A.3)
The density fluctuation level, which is one of the parameters of our model, is therefore (A.5)
For instance, the n_{H} field whose distribution function is shown in Fig. A.1 has values ranging from 0.3 cm^{−3} to 880 cm^{−3}, with a mean and standard deviation of , resulting in the desired fluctuation level y_{n} = 1.
Fig. A.1 Distribution function of the synthetic field n_{H} used for Fig. 1 (black histogram). The red curve shows the theoretical probability distribution function f(n_{H}) for the chosen set of parameters. 

Open with DEXTER 
Fig. A.2 Power spectra of synthetic n_{H} fields obtained by exponentiation of a 120 × 120 × 120 pixels fractional Brownian motion with spectral index β_{X} = 3. The fluctuation levels y_{n} are specified next to each curve, and the original field’s power spectrum is represented as a dashed line. 

Open with DEXTER 
Power spectra
Figure A.2 shows the azimuthallyaveraged power spectra of n_{H} fields obtained through Eq. (4), from a 120 × 120 × 120 pixels fractional Brownian motion X with spectral index β_{X} = 3, for various fluctuation levels y_{n}. The spectra were normalized differently so as to allow comparison between them.
The powerlaw 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. A.3, which shows the differences β_{n} − β_{X} between the spectral indices of the original fBm field X and that of the model density field n_{H}. The power spectra of the latter are indeed flatter than the original ones (β_{n} < β_{X}), with differences that may become large when β_{X} is low, but remain negligible for higher β_{X}. We interpret this trend as the exponentiation process amplifying the twopoint differences X(r+ δ) − X(r), for small separations δ , that exist in the original field when β_{X} is low, leading to an increase of the smallscale power in the model n_{H} field, i.e. β_{n} < β_{X}. This effect is all the more important when exponentiation stretches these field differences more strongly, i.e. when y_{n} increases. On the other hand, for high β_{X}, the fBm fields are much smoother, so the exponentiation process has little impact on these twopoint statistics at small scales, leading to β_{n} ≃ β_{X}. We note that at low fluctuation levels (y_{n} ≤ 0.3), the differences are smaller than 0.1 for all values of β_{X}.
Fig. A.3 Evolution of the differences β_{n} − β_{X} between the spectral indices of the original fBm field X and that of the model density field n_{H} with fluctuation level y_{n}. Each point corresponds to the mean of 20 realisations of the model density field n_{H}, and the error bars represent the standard deviation of the fitted spectral indices β_{n} and fluctuation levels y_{n}. 

Open with DEXTER 
Appendix B: Likelihood terms
B.1 D^{2} terms for mean values
The first term in Eq. (25) is given by (B.1)
The mean value of the optical depth τ_{353} is evaluated on the simulated and Planck data maps. These are compared through a standard χ^{2} test. The denominator includes the uncertainty on the mean value propagated from the uncertainty map provided by the Planck collaboration (Planck Collaboration XI 2014), and the uncertainty coming from the conversion factor used to build the simulated map (Planck Collaboration XI 2014).
In the second term, is the mean value of the temperature map T_{obs} from Planckdata, and represents the uncertainty on the averaged value propagated from the uncertainty map provided by the Planck collaboration. T_{d} is directly the model parameter for the dust temperature.
B.2 D^{2} terms for distribution functions
For a given observable map o, we compute its DF over an ensemble of N_{b} bins^{13}. When considering the Planck data, we write this DF as , 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 . When considering the model, we write and 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 of the observable’s DF to the total D^{2} in Eq. (25) is then built as an average over the N_{b} bins (B.2)
The sum is normalized to the number of bins so that 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 is minimum.
This quantity is different from a standard χ^{2} test as it compares data with one random^{14} realisation for a given set of model parameters, and this comparison is repeated and averaged N_{r} times. A standard χ^{2} test would compare data with a model prediction that would be the average of the random realizations (see Eq (C.3) in Appendix C). 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 β), only the few largescale modes effectively contribute to the power, leading to a large variance in each DF bin^{15}. Thus the χ^{2} test tends to favour these large values of β, 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.
B.3 D^{2} terms for power spectra
For a given observable map o, we first compute its twodimensional power spectrum (B.3)
then average these within annuli in Fourier space, centred on a set of wavenumbers . We write for this azimuthal average in bin i when considering the Planck data, and when considering model fields. The uncertainties affecting these quantities strongly depend on the number N_{i} of wavevectors k_{n} in each bin. The best estimate for the standard deviation of the power spectrum of the data is (B.4)
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 for the 60 model realizations is computed in the same way as for data.
The contribution of the observable’s power spectrum to the total D^{2} in Eq. (25) is then computed as an average over the bins^{16} in wavenumberspace, (B.5)
B.4: D^{2} term for the anticorrelation
To use the anticorrelation (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XX 2015), we compute the joint distribution function of the and p_{MAS} maps, which we write and for the Planck data and model maps respectively, with 1 ≤ i ≤ N_{b,1} and 1 ≤ j ≤ N_{b,2} the binning scheme used for the two maps. The standard deviations and are defined in the same way as for the onedimensional DFs in Sect. B.2, and the contribution to the total D^{2} is then (B.6)
In this expression, we note that the total number N_{b,tot} of twodimensional bins considered is less than the product N_{b,1}N_{b,2} of the number of bins in each dimension, which we set to N_{b,1} = N_{b,2} = 50. The reason for this is that we discard the empty bins and those with a signaltonoise ratio below three^{17}. We thus keep only the significantly populated bins that can drive the fit and contribute to the total D^{2}.
Appendix C: Goodnessoffit
To assess the goodness of the fit, we use an a posteriori χ^{2} test, which we define as (C.1)
with N_{o} = 13 the total number of observables. Each term is a χ^{2} test comparing data with the mean of the N_{r} = 60 realisations. The first term from Eq (C.1) is (C.2)
where 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 quantity is the variance of over the N_{r} realisations (capital Σ denotes the variance over the random realisations).
is the ensemble average of the i^{th} bin of the DF for the observable o, over the N_{r} = 60 model realizations, and is the associated variance. We note that while in Eq. (B.2) is the average of the observable χ^{2}, is the χ^{2} of the averaged observable.
where is the averaged power spectrum in bin i and its variance.
Finally, the fourth term is (C.6)
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 χ^{2} quantity defined in Eq. (C.1). The average of these 100 χ^{2} values^{18} is listed as in Tables 3 and 4.
References
 Alina, D., Montier, L., Ristorcelli, I., et al. 2016, A&A, 595, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arzoumanian, D., André, P., Didelon, P., et al. 2011, A&A, 529, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benoît, A., Ade, P., Amblard, A., et al. 2004, A&A, 424, 571 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beresnyak, A. 2014, ApJ, 784, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Bierman, E. M., Matsumura, T., Dowell, C. D., et al. 2011, ApJ, 741, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Brooks, S., Gelman, A., Jones, G. L., & Meng, X. L. 2011, Handbook of Markov Chain Monte Carlo (Boca Raton: Chapman & Hall/CRC) [CrossRef] [Google Scholar]
 Brunt, C. M., & Heyer, M. H. 2002, ApJ, 566, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Burkhart, B., & Lazarian, A. 2016, ApJ, 827, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Burkhart, B., FalcetaGonçalves, D., Kowal, G., & Lazarian, A. 2009, ApJ, 693, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Burkhart, B., Lazarian, A., Leão, I. C., de Medeiros, J. R., & Esquivel, A. 2014, ApJ, 790, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Cho, J., & Lazarian, A. 2003, MNRAS, 345, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, C., Lazarian, A., Burkhart, B., Pogosyan, D., & De Medeiros J. R. 2016, ApJ, 818, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Dotson, J. L., Vaillancourt, J. E., Kirby, L., et al. 2010, ApJS, 186, 406 [NASA ADS] [CrossRef] [Google Scholar]
 FalcetaGonçalves, D., Lazarian, A., & Kowal, G. 2008, ApJ, 679, 537 [NASA ADS] [CrossRef] [Google Scholar]
 Falconer, K. 1990, Fractal Geometry: Mathematical Foundations and Applications (Chichester: John Wiley & Sons) [Google Scholar]
 Falgarone, E., Panis, J.F., Heithausen, A., et al. 1998, A&A, 331, 669 [NASA ADS] [Google Scholar]
 Fissel, L. M., Ade, P. A. R., Angilè, F. E., et al. 2016, ApJ, 824, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Gelman, A., & Rubin, D. B. 1992, Statist. Sci., 7, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Hamann, J., Hannestad, S., Raffelt, G. G., & Wong, Y. Y. Y. 2007, JCAP, 8, 021 [NASA ADS] [CrossRef] [Google Scholar]
 Haverkorn, M., Brown, J. C., Gaensler, B. M., & McClureGriffiths, N. M. 2008, ApJ, 680, 362 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P. 2013, A&A, 556, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennebelle, P., & Falgarone, E. 2012, A&ARv, 20, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P., & Iffrig, O. 2014, A&A, 570, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennebelle, P., Banerjee, R., VázquezSemadeni, E., Klessen, R. S., & Audit, E. 2008, A&A, 486, L43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hildebrand, R. H., Kirby, L., Dotson, J. L., Houde, M., & Vaillancourt, J. E. 2009, ApJ, 696, 567 [NASA ADS] [CrossRef] [Google Scholar]
 HilyBlant, P.,& Falgarone, E. 2009, A&A, 500, L29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, ApJS, 213, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Inutsuka, S.i., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, A&A, 580, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jaffe, T. R., Leahy, J. P., Banday, A. J., et al. 2010, MNRAS, 401, 1013 [NASA ADS] [CrossRef] [Google Scholar]
 Jansson, R., & Farrar, G. R. 2012, ApJ, 761, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Koch, P. M., Tang, Y.W., Ho, P. T. P., et al. 2014, ApJ, 797, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Levrier, F., Falgarone, E., & Viallefond, F. 2006, A&A, 456, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matthews, B. C., McPhee, C. A., Fissel, L. M., & Curran, R. L. 2009, ApJS, 182, 143 [NASA ADS] [CrossRef] [Google Scholar]
 MivilleDeschênes, M.A., Levrier, F., & Falgarone, E. 2003, ApJ, 593, 831 [NASA ADS] [CrossRef] [Google Scholar]
 MivilleDeschênes, M.A., Ysard, N., Lavabre, A., et al. 2008, A&A, 490, 1093 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MivilleDeschênes, M.A., Martin, P. G., Abergel, A., et al. 2010, A&A, 518, L104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Montier, L., Plaszczynski, S., Levrier, F., et al. 2015a, A&A, 574, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Montier, L., Plaszczynski, S., Levrier, F., et al. 2015b, A&A, 574, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 O’Dea, D. T., Clark, C. N., Contaldi, C. R., & MacTavish, C. J. 2012, MNRAS, 419, 1795 [NASA ADS] [CrossRef] [Google Scholar]
 Panopoulou, G. V., Psaradaki, I., & Tassis, K. 2016, MNRAS, 462, 1517 [NASA ADS] [CrossRef] [Google Scholar]
 Perez, J. C., Mason, J., Boldyrev, S., & Cattaneo, F. 2012, Phys. Rev. X, 2, 041005 [Google Scholar]
 Planck Collaboration I. 2016, A&A, 594, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XIX. 2015, A&A, 576, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XX. 2015, A&A, 576, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXX. 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXII. 2016, A&A, 586, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXV. 2016, A&A, 586, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXVIII. 2016, A&A, 586, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLIV. 2016, A&A, 596, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Plaszczynski, S., Montier, L., Levrier, F., & Tristram, M. 2014, MNRAS, 439, 4048 [NASA ADS] [CrossRef] [Google Scholar]
 Ponthieu, N., MacíasPérez, J. F., Tristram, M., et al. 2005, A&A, 444, 327 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 786, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Seifried, D., & Walch, S. 2015, MNRAS, 452, 2410 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, J. D., Hennebelle, P., Martin, P. G., et al. 2013, ApJ, 774, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, J. D., Alves, F., Boulanger, F., et al. 2016, A&A, 596, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., & Zielinsky, M. 1998, A&A, 336, 697 [NASA ADS] [Google Scholar]
 Vaillancourt, J. E., & Matthews, B. C. 2012, ApJS, 201, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Vansyngel, F., Boulanger, F., Ghosh, T., et al. 2017, A&A, 603, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 WardThompson, D., Sen, A. K., Kirk, J. M., & Nutter, D. 2009, MNRAS, 398, 394 [NASA ADS] [CrossRef] [Google Scholar]
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).
We do not consider an ordered random or striated random component of the field (Jaffe et al. 2010; Jansson & Farrar 2012), which we justify with the smallness of the fieldofview considered.
A more recent determination of the distance to the Polaris Flare places it at 350–400 pc (Schlafly et al. 2014). For the demonstration of the method presented here, this is not a critical issue, as the powerlaw power spectra underline selfsimilar behaviours, so that a change of the distance can be compensated by a change in pixel size.
Incidentally, from the Planck maps, we can measure the spectral index of the total intensity for the Polaris Flare to be β_{I} = 2.84 ± 0.10, in excellent agreement with the measurement by Stutzki et al. (1998) on CO integrated emission at a similar angular resolution.
To have reliable and smooth power spectra with 120 × 120 pixel maps, we set initially but later cut off wavenumbers larger than k_{max} = 2π∕(3 × 15′), which corresponds to scales smaller than three beam sizes. Indeed, for bins k_{i} > k_{max} the power spectrum is completely washed out by the beam convolution (see Fig. 5) and contains no information about the underlying interesting parameters. This uninformative part is thus removed from , and thus .
All Tables
All Figures
Fig. 1 Total gas column density N_{H}, derived froma synthetic density field n_{H} built by exponentiation of a fBm field with spectral index β_{X} = 2.6 and size 120 × 120 × 120 pixels. The volume density fluctuation level is y_{n} = 1, and the column density fluctuation level is . 

Open with DEXTER  
In the text 
Fig. 2 Synthetic magnetic field B built using Eq. (8). The spectral index of the vector potential is β_{A} = 5 and the size of the cubes is 120 × 120 × 120 pixels, corresponding to 30 pc on each side. Shown are twodimensional 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 B_{t} to the norm of the uniform component B_{0} is y_{B} =1 in this particular case, with angles χ_{0} = γ_{0} = 0°. 

Open with DEXTER  
In the text 
Fig. 3 Distribution functions of the components B_{x}, B_{y}, and B_{z} of a model magnetic field B = B_{0} + B_{t} built on a grid 120 × 120 × 120 pixels using Eq. (8) with β_{A} = 5, and a mean, largescale magnetic field B_{0} defined by the angles χ_{0} = 0° and γ_{0} = 60°, and a norm B_{0} = 50 μG such that the fluctuation level is y_{B} = 0.1. The vertical lines represent the projected values of the large scale magnetic field B_{0x} = B_{0} sinχ_{0} cosγ_{0}, B_{0y} = −B_{0} cosχ_{0} cosγ_{0} and B_{0z} = B_{0} sinγ_{0}. See Fig. 14 in Planck Collaboration Int. XX (2015) for the definition of angles. 

Open with DEXTER  
In the text 
Fig. 4 Powerspectra of the components B_{x}, B_{y} and B_{z} of a model magnetic field B = B_{0} + B_{t} built on a grid 120 × 120 × 120 pixels using Eq. (8) with β_{A} = 3 (different shades of blue for the three components) and β_{A} = 5 (different shades of red for the three components). The power spectra are normalized differently so as to allow comparison between them. The fitted powerlaws shown as solid lines yield spectral indices β_{B} = 1 and β_{B} = 3, in agreement with Eq. (11). These are the power spectra of the same particular realizations shown in Fig. 2. 

Open with DEXTER  
In the text 
Fig. 5 Left column row: example maps (from top to bottom: total intensity I_{m} on a logarithmic scale, Stokes Q_{m}, and Stokes U_{m}) from simulation A (see Sect. 5 and Table 3). Right column: corresponding power spectra. The grey points represent the twodimensional power spectra, while the black dots represent the azimuthal averages in Fourier space in a set of wavenumber bins, and the blue line is a powerlaw fit to the black points. 

Open with DEXTER  
In the text 
Fig. 6 Constraints (posterior probability contours and marginalized PDFs) on the statistical properties of the dust density and magnetic field for the simulation A maps. On the posterior probability contours, the filled dark and light blue regions respectively enclose 68.3% and 95.4% of the probability, the black stars indicate the averages over the twodimensional posterior PDFs, and the red circles indicate the input values for the simulation. In the plots showing the marginalized posterior PDFs, the light blue regions enclose 68.3% of the probability, the dashed blue lines indicate the averages over the posterior PDFs, and the solid red lines indicate the input values for the simulation. The upper right plot displays the correlation matrix between the fitted parameters. 

Open with DEXTER  
In the text 
Fig. 7 Planck 353 GHz maps of the Polaris Flare molecular cloud. The top row shows, from left to right, the total intensity I_{353} on a logarithmic scale, the Stokes Q_{353} map and the Stokes U_{353} map, while the bottom row shows the polarization fraction p_{MAS}, the polarization angle ψ, and the polarization angle dispersion function . The τ_{353} and T_{obs} maps have the same aspects as the I_{353} map but with their own scales. 

Open with DEXTER  
In the text 
Fig. 8 Constraints (posterior probability contours and marginalized PDFs) on the statistical properties of the dust density and magnetic field for the Planck maps of the Polaris Flare. On the posterior probability contours, the filled dark and light blue regions respectively enclose 68.3% and 95.4% of the probability, and the black stars indicate the averages over the twodimensional posterior PDFs. In the plots showing the marginalized posterior PDFs, the light blue regions enclose 68.3% of the probability, and the dashed blue lines indicate the averages over the posterior PDFs. The upper right plot displays the correlation matrix between the fitted parameters. 

Open with DEXTER  
In the text 
Fig. 9 Same as Fig. 7 with the same colour scales, but for model maps using the best fitting parameters to the Polaris Flare data. 

Open with DEXTER  
In the text 
Fig. 10 Comparison of the DFs extracted from the Planck Polaris Flare maps (black points) with the observables computed from simulations using the best fitting parameters (blue curves). The latter curves are averaged over 60 realizations, as describedin Sect. 4.2: the average is given by the central blue curve and the shaded bands give the 1σ and 2σ standard deviation in each bin. 

Open with DEXTER  
In the text 
Fig. 11 Comparison of the I_{m}, Q_{m}, and U_{m} power spectra extracted from the Planck Polaris Flare maps (grey points representing the twodimensional power spectra, and black dots representing the azimuthal averages in Fourier space) with the observables computed from simulations using the best fitting parameters (blue curves). The latter curves are averaged over 60 realizations as described in Sect 4.2: the average is given by the central blue curve and the narrow shaded bands give the 1σ and 2σ standard deviation in each bin. 

Open with DEXTER  
In the text 
Fig. 12 Twodimensional distribution function of and polarization fraction p_{MAS} for the Polaris Flare maps (left), for the model maps using the best fitting parameters averaged over 60 realizations (middle), and residuals (right). The polarization angle dispersion function is computed at a lag δ = 16′. The solid black line shows the mean for each bin in p_{MAS} and the dashed black line is a linear fit of that curve, restricted to bins in p_{MAS} which contain at least 1% of the total number of points (120 × 120). 

Open with DEXTER  
In the text 
Fig. A.1 Distribution function of the synthetic field n_{H} used for Fig. 1 (black histogram). The red curve shows the theoretical probability distribution function f(n_{H}) for the chosen set of parameters. 

Open with DEXTER  
In the text 
Fig. A.2 Power spectra of synthetic n_{H} fields obtained by exponentiation of a 120 × 120 × 120 pixels fractional Brownian motion with spectral index β_{X} = 3. The fluctuation levels y_{n} are specified next to each curve, and the original field’s power spectrum is represented as a dashed line. 

Open with DEXTER  
In the text 
Fig. A.3 Evolution of the differences β_{n} − β_{X} between the spectral indices of the original fBm field X and that of the model density field n_{H} with fluctuation level y_{n}. Each point corresponds to the mean of 20 realisations of the model density field n_{H}, and the error bars represent the standard deviation of the fitted spectral indices β_{n} and fluctuation levels y_{n}. 

Open with DEXTER  
In the text 