EDP Sciences
Free Access
Issue
A&A
Volume 614, June 2018
Article Number A124
Number of page(s) 19
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201732128
Published online 25 June 2018

© 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; 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 Planck1, 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% (Planck Collaboration Int. XIX 2015), confirming results previously obtained over one fifth of the sky by the Archeops balloon-borne experiment (Benoît et al. 2004; Ponthieu et al. 2005). Furthermore, the polarization fraction is anti-correlated 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 NH 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 anti-correlations 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 trans-Alfvé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 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 (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 = B0 + Bt of a uniform, large-scale component B0, and a turbulent component Bt that has a null spatial average, ⟨Bt⟩ = 0. The statistical properties in question are then essentially modelled by two quantities, i) the ratio of the turbulent component to the mean, yB = σBB0, where and B0 = ||B0||, and ii) the spectral index βB, which characterizes the distribution of power of Bt 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 NH 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-Gonçalves et al. (2008) and Hildebrand et al. (2009), and from their results, we estimated the ratio yB 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 yB 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 large-scale field B0 and a small number (Nl ≃ 7) of independent “polarization layers” on the LOS, each of which accounts for a fraction 1∕Nl of the total unpolarized emission. Within each layer, the turbulent component Bt of the magnetic field, which is used to compute synthetic Stokes Q and U maps, is an independent realization of a Gaussian two-dimensional random field with a prescribed spectral index βB. Through these fits, they confirm the near equipartition of large-scale and turbulent components of B, with yB ≃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 yB and βB. Orthogonal approaches have also been pursued (e.g. Miville-Deschê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 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, 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 nd and magnetic field B produced bynumerical simulations of MHD turbulence. One could imagine generalizing 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 nd and the magnetic field B, allowing us to perfectly control the statistical properties of these three-dimensional 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, βB, and the ratio of turbulent to regular field, yB.

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 nd 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 three-dimensional cubes of dust density nd and magnetic field B, which we build using fractional Brownian motions (fBm) (Falconer 1990). An N-dimensional fBm X is a random field defined on such that , for any pair of points (r1, r2). H is called the Hurst exponent. These fBm fields are usually built in Fourier space3, (1)

by specifying amplitudes that scale as a power-law 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 real-valued. 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; Miville-Deschênes et al. 2003; Correia et al. 2016).

2.2 Dust density

In our approach, the dust density nd is taken to be proportional to the total gas density nH, so that the dust optical depth within each cell is also proportional to nH (see the derivation of polarization maps in Sect. 2.4). Therefore, we intended to model nH from numerical realizations of three-dimensional fBm fields built in Fourier space. These have means that are defined by the value chosen for the null-wavevector amplitude A(0), so if one wished to use such a synthetic random field X directly as a model for the positive-valued nH, one would be required to choose nH = X′ = Xa 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 three-dimensional density field nH.

We remedied this shortcoming by taking X to represent the log-density, i.e., nH is given by (4)

where X is a three-dimensional fBm field with spectral index βX, and Xr and n0 are positive parameters. The nH 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 Xr, 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 power-law 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 NH derived from a gas volume density nH 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 Xr = 1.2 and n0 = 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; Miville-Deschênes 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 the Fourier phases (see, e.g., Levrier et al. 2006; Burkhart & Lazarian 2016).

thumbnail Fig. 1

Total gas column density NH, derived froma synthetic density field nH 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 yn = 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 Bt with null divergence and controlled power spectrum, we started from a vector potential A built as a three-dimensional 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 Bt,λ = ϵλμνμAν, with the Einstein notation, where ϵλμν is the Levi-Civita tensor, and the derivation relation in Fourier space (7)

we have the expression of the components of Bt in Fourier space (8)

As it should, this expression corresponds to a divergence-free turbulent magnetic field, (9)

Writing kμ = kfμ, with , the power spectrum of each component of Bt is then (10)

The last factor is essentially independent of the wavenumber k, so the spectral index of each component of Bt is . After Fourier-transforming back to real space, Bt 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 uniform4 vector field B0 to that turbulent magnetic field Bt. 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 Bt,λ, 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 B0 which is added to the turbulent field Bt is defined byits norm B0 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 B0 in the POS, counted positively clockwise from the north-south direction (see Fig. 14 of Planck Collaboration Int. XX 2015). The total magnetic field’s direction in three-dimensional 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, yB = 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.

thumbnail 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 two-dimensional slices through the cubes of the components Bx (top), By (middle), and Bz (bottom). The ratio of the fluctuations of the turbulent component Bt to the norm of the uniform component B0 is yB =1 in this particular case, with angles χ0 = γ0 = 0°.

Open with DEXTER
thumbnail Fig. 3

Distribution functions of the components Bx, By, and Bz of a model magnetic field B = B0 + Bt built on a grid 120 × 120 × 120 pixels using Eq. (8) with βA = 5, and a mean, large-scale magnetic field B0 defined by the angles χ0 = 0° and γ0 = 60°, and a norm B0 = 50 μG such that the fluctuation level is yB = 0.1. The vertical lines represent the projected values of the large scale magnetic field B0x = B0 sinχ0 cosγ0, B0y = −B0 cosχ0 cosγ0 and B0z = B0 sinγ0. See Fig. 14 in Planck Collaboration Int. XX (2015) for the definition of angles.

Open with DEXTER
thumbnail Fig. 4

Powerspectra of the components Bx, By and Bz of a model magnetic field B = B0 + Bt 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 power-laws 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
thumbnail Fig. 5

Left column row: example maps (from top to bottom: total intensity Im on a logarithmic scale, Stokes Qm, and Stokes Um) from simulation A (see Sect. 5 and Table 3). Right column: corresponding power spectra. The grey points represent the two-dimensional 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 power-law fit to the black points.

Open with DEXTER

2.4 Polarization maps

Once cubes of total gas density nH 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 signal-to-noise 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 p0 to be uniform, and the source function Sν = Bν(Td) to be that of a black body with an assumed uniform dust temperature Td. The dust opacity at this frequency σ353 is taken to vary with NH, following Fig. 20 in Planck Collaboration XI (2014) for XCO = 2 1020 H2 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 cm2. The values of NH considered in our study are typically at most a few 1021 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 angle5 χ 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 nH and B cubes are built on a grid which is 132 × 132 pixels in the POS and Nz 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 = Nz δ 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 Archive6 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 nI, nQ, and nU are random values drawn from a three-dimensional 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 In, Qn, and Un maps are then placed at a distance7 D = 140 pc, so that the angular size of each pixel is about 6′, and then convolved by a circular 15′ full-width 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 (Im, Qm, Um) 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 ⟨Im⟩ is the spatial average of the model Stokes Im 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 b2 derives from the elements of the noise covariance matrix Σ (see Montier et al. 2015b). Next, we define the polarization angle (21)

where the two-argument 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 conversion8 from Planck Collaboration XI (2014). A temperature map Tobs is also created using the anti-correlation with the column density NH observed in the data. This “dust temperature map” does not pretend to model reality but since Td is one of the parameters of the model, the fitting algorithm requires a map whose mean value should yield Td.

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 one-point and two-point 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 yB and γ0 which cannot be lifted. Consequently, we chose not to try to fit for γ0 and yB, 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 ⟨NH⟩≈ 1021 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 ⟨nH ⟩ 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 yn 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 Td and p0 are very large and do not play a role in the fitting procedure.

Table 1

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 Im, Qm, and Um, 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 Tobs, ii) the distribution functions (one-point statistics) of the Im, Qm, Um, pMAS, ψ, , and maps, iii) the power spectra (two-point statistics) of the Im,Qm, and Um maps, and iv) the pixel-by-pixel anti-correlation between and pMAS underlined by Planck Collaboration Int. XIX (2015). Indeed, we find 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.

On the simulation side, Nr = 60 model realizations per set of parameter values are generated with their observables, to be compared with data. The Nr 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 D2 which quantifies the distance between data and one random realisation of the model, averaged over the Nr 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 two-dimensional joint DFs of and pMAS in the data and in synthetic maps. We detail the computation of these various terms in Appendix B.

Table 2

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 Metropolis-Hastings 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 D2 criterion described in Sect. 4.2 and Appendix B as (26)

with π(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 sq−1 at step q − 1. If the likelihood for the new set of parameters, sq, 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 sq set of parameters is kept, otherwise the chain duplicates the sq−1 set, i.e., sq = sq−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 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 Gelman-Rubin 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 least-converged parameter. The convergence is also assessed by visually checking the D2 and parameter evolutions along the chains.

The obtained nine-dimensional 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 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 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).

Table 3

Best fit values from four fBm simulations using the observables from Table 2.

5 Results

5.1 Validation of the method

To validate the fitting method, we simulated four sets of model cubes and computed the corresponding Im, Qm, and Um 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 reached9. 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 p0, or between yn and ⟨nH⟩. These trends are best visualized with the correlation matrix, shown in the upper right corner of Fig. 6.

thumbnail 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 two-dimensional 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
thumbnail Fig. 7

Planck 353 GHz maps of the Polaris Flare molecular cloud. The top row shows, from left to right, the total intensity I353 on a logarithmic scale, the Stokes Q353 map and the Stokes U353 map, while the bottom row shows the polarization fraction pMAS, the polarization angle ψ, and the polarization angle dispersion function . The τ353 and Tobs maps have the same aspects as the I353 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, 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, for example, Miville-Deschê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; Hily-Blant & 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 full-mission Planck maps of Stokes parameters (I353, Q353, U353) 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 Tobs from the 2013 public release (Planck Collaboration XI 2014). All maps are at a native 4′.8 resolution in the Healpix format with Nside = 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 ofI353, Q353, U353, 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 Monte-Carlo simulations of pure noise maps, drawn from the original full-resolution covariance maps and smoothed at 15′. The maps of I353, Q353, U353, pMAS, ψ, 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 signal-to-noise ratio per pixel in the model Im, Qm and Um 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 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, yn ≈ 1.6, and the magnitude of the large-scale 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 ⟨nH⟩≈ 40 cm−3. The temperature Td is 17.5 K equal to the average of the Tobs Planck map, and the polarization fraction is p0 ≈ 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. 1012). The agreement is excellent for most observables, although substantial deviations are visible in the DFs of the intensity Im, 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 large-scale 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 Im 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 Td = 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 large-scale 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 pint ≈ 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°.

thumbnail 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 two-dimensional 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
thumbnail 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
Table 4

Best fit values for the Planck Polaris Flare maps, using the observables from Table 2.

thumbnail 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
thumbnail Fig. 11

Comparison of the Im, Qm, and Um power spectra extracted from the Planck Polaris Flare maps (grey points representing the two-dimensional 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
thumbnail Fig. 12

Two-dimensional distribution function of and polarization fraction pMAS 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 pMAS and the dashed black line is a linear fit of that curve, restricted to bins in pMAS 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 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 β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 exponent10. 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 turbulent-to-uniform magnetic field are both around unity. Finally, our analysis is able to give a constraint on the polarization fraction, p0 ≈ 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 ⟨nH⟩≈ 40 cm−3. The good visual agreement between the Polaris Flare maps and model maps for the best-fitting parameters (Figs. 7and 9), and the excellent agreement between the two sets of maps for most of the observables used in the analysis (Figs. 1012), 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-Deschê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, high-latitude 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 NH 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 polarizationmaps built from a numerical simulation of MHD turbulence. We used simulation cubes11 (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 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 ten pixels), and the properties at large scales are dependent on the forcing, which is user-defined. 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 co-funded by CEA and CNES.

Appendix A Statistical properties of nH models

A.1 Probability distribution function

The PDF f(nH) of the density field nH 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 log-normal PDF (A.1)

which is defined for nH > 0. Figure A.1 presents the distribution function of the nH field12 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 nH (A.2)

which allows, in particular, to compute its mean value (A.3)

as well as its variance (A.4)

The density fluctuation level, which is one of the parameters of our model, is therefore (A.5)

For instance, the nH 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 yn = 1.

thumbnail Fig. A.1

Distribution function of the synthetic field nH used for Fig. 1 (black histogram). The red curve shows the theoretical probability distribution function f(nH) for the chosen set of parameters.

Open with DEXTER
thumbnail Fig. A.2

Power spectra of synthetic nH fields obtained by exponentiation of a 120 × 120 × 120 pixels fractional Brownian motion with spectral index βX = 3. The fluctuation levels yn 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 azimuthally-averaged power spectra of nH fields obtained through Eq. (4), from a 120 × 120 × 120 pixels fractional Brownian motion X with spectral index βX = 3, for various fluctuation levels yn. The spectra were normalized differently so as to allow comparison between them.

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. 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 nH. 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 two-point differences X(r+ δ) − X(r), for small separations ||δ ||, that exist in the original field when βX is low, leading to an increase of the small-scale power in the model nH field, i.e. βn < βX. This effect is all the more important when exponentiation stretches these field differences more strongly, i.e. when yn increases. On the other hand, for high βX, the fBm fields are much smoother, so the exponentiation process has little impact on these two-point statistics at small scales, leading to βnβX. We note that at low fluctuation levels (yn ≤ 0.3), the differences are smaller than 0.1 for all values of βX.

thumbnail 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 nH with fluctuation level yn. Each point corresponds to the mean of 20 realisations of the model density field nH, and the error bars represent the standard deviation of the fitted spectral indices βn and fluctuation levels yn.

Open with DEXTER

Appendix B: Likelihood terms

B.1 D2 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 Tobs from Planckdata, and represents the uncertainty on the averaged value propagated from the uncertainty map provided by the Planck collaboration. Td is directly the model parameter for the dust temperature.

B.2 D2 terms for distribution functions

For a given observable map o, we compute its DF over an ensemble of Nb bins13. 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 Nr = 60 model realizations.

The contribution of the observable’s DF to the total D2 in Eq. (25) is then built as an average over the Nb 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 random14 realisation for a given set of model parameters, and this comparison is repeated and averaged Nr 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 large-scale modes effectively contribute to the power, leading to a large variance in each DF bin15. 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 D2 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 D2 terms for power spectra

For a given observable map o, we first compute its two-dimensional 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 Ni of wavevectors kn in each bin. The best estimate for the standard deviation of the power spectrum of the data is (B.4)

where the factor tN 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 D2 in Eq. (25) is then computed as an average over the bins16 in wavenumberspace, (B.5)

B.4: D2 term for the anti-correlation

To use the anti-correlation (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XX 2015), we compute the joint distribution function of the and pMAS maps, which we write and for the Planck data and model maps respectively, with 1 ≤ iNb,1 and 1 ≤ jNb,2 the binning scheme used for the two maps. The standard deviations and are defined in the same way as for the one-dimensional DFs in Sect. B.2, and the contribution to the total D2 is then (B.6)

In this expression, we note that the total number Nb,tot of two-dimensional bins considered is less than the product Nb,1Nb,2 of the number of bins in each dimension, which we set to Nb,1 = Nb,2 = 50. The reason for this is that we discard the empty bins and those with a signal-to-noise ratio below three17. We thus keep only the significantly populated bins that can drive the fit and contribute to the total D2.

Appendix C: Goodness-of-fit

To assess the goodness of the fit, we use an a posteriori χ2 test, which we define as (C.1)

with No = 13 the total number of observables. Each term is a χ2 test comparing data with the mean of the Nr = 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 Nr realisations. The quantity is the variance of over the Nr realisations (capital Σ denotes the variance over the random realisations).

The second term is (C.3)

where (C.4)

is the ensemble average of the ith bin of the DF for the observable o, over the Nr = 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.

The third term is (C.5)

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 values18 is listed as in Tables 3 and 4.

References


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 Bt, or of the modulus |Bt|. Assuming that Bt is isotropic, which we will, all of these spectral indices are identical.

3

Throughout this paper, for any field F the notation represents its Fourier transform.

4

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 field-of-view considered.

5

Not to be confused with the corresponding position angle χ0 of the uniform component of the magnetic field B0.

7

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 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 NH at a resolution of 30′. Thus, before applying it pixel by pixel, we smooth the simulated NH map to 30′ resolution, apply the conversion, then resample the resulting τ353 map at the original resolution.

9

And after removing the burn-in phase, which is quite short in our case (≲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 β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.

12

In that case, ⟨X⟩ = 0, , Xr = 1.2, and n0 = 20 cm−3.

13

We set Nb = 50 for the 120 × 120 pixel maps.

14

We remind the reader that the randomness comes from the fBm itself and the noise addition to the Stokes maps.

15

This is akin to the cosmic variance problem in cosmology.

16

To have reliable and smooth power spectra with 120 × 120 pixel maps, we set initially but later cut off wavenumbers larger than kmax = 2π∕(3 × 15′), which corresponds to scales smaller than three beam sizes. Indeed, for bins ki > kmax 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 .

17

i.e., bins where , , , or .

18

Each simulation is a random realization of the fBm and the noise, then the χ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 χ2 value.

All Tables

Table 1

Parameter space explored in the grid of model polarization maps.

Table 2

Observables from polarization maps used to fit data.

Table 3

Best fit values from four fBm simulations using the observables from Table 2.

Table 4

Best fit values for the Planck Polaris Flare maps, using the observables from Table 2.

All Figures

thumbnail Fig. 1

Total gas column density NH, derived froma synthetic density field nH 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 yn = 1, and the column density fluctuation level is .

Open with DEXTER
In the text
thumbnail 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 two-dimensional slices through the cubes of the components Bx (top), By (middle), and Bz (bottom). The ratio of the fluctuations of the turbulent component Bt to the norm of the uniform component B0 is yB =1 in this particular case, with angles χ0 = γ0 = 0°.

Open with DEXTER
In the text
thumbnail Fig. 3

Distribution functions of the components Bx, By, and Bz of a model magnetic field B = B0 + Bt built on a grid 120 × 120 × 120 pixels using Eq. (8) with βA = 5, and a mean, large-scale magnetic field B0 defined by the angles χ0 = 0° and γ0 = 60°, and a norm B0 = 50 μG such that the fluctuation level is yB = 0.1. The vertical lines represent the projected values of the large scale magnetic field B0x = B0 sinχ0 cosγ0, B0y = −B0 cosχ0 cosγ0 and B0z = B0 sinγ0. See Fig. 14 in Planck Collaboration Int. XX (2015) for the definition of angles.

Open with DEXTER
In the text
thumbnail Fig. 4

Powerspectra of the components Bx, By and Bz of a model magnetic field B = B0 + Bt 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 power-laws 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
thumbnail Fig. 5

Left column row: example maps (from top to bottom: total intensity Im on a logarithmic scale, Stokes Qm, and Stokes Um) from simulation A (see Sect. 5 and Table 3). Right column: corresponding power spectra. The grey points represent the two-dimensional 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 power-law fit to the black points.

Open with DEXTER
In the text
thumbnail 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 two-dimensional 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
thumbnail Fig. 7

Planck 353 GHz maps of the Polaris Flare molecular cloud. The top row shows, from left to right, the total intensity I353 on a logarithmic scale, the Stokes Q353 map and the Stokes U353 map, while the bottom row shows the polarization fraction pMAS, the polarization angle ψ, and the polarization angle dispersion function . The τ353 and Tobs maps have the same aspects as the I353 map but with their own scales.

Open with DEXTER
In the text
thumbnail 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 two-dimensional 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
thumbnail 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
thumbnail 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
thumbnail Fig. 11

Comparison of the Im, Qm, and Um power spectra extracted from the Planck Polaris Flare maps (grey points representing the two-dimensional 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
thumbnail Fig. 12

Two-dimensional distribution function of and polarization fraction pMAS 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 pMAS and the dashed black line is a linear fit of that curve, restricted to bins in pMAS which contain at least 1% of the total number of points (120 × 120).

Open with DEXTER
In the text
thumbnail Fig. A.1

Distribution function of the synthetic field nH used for Fig. 1 (black histogram). The red curve shows the theoretical probability distribution function f(nH) for the chosen set of parameters.

Open with DEXTER
In the text
thumbnail Fig. A.2

Power spectra of synthetic nH fields obtained by exponentiation of a 120 × 120 × 120 pixels fractional Brownian motion with spectral index βX = 3. The fluctuation levels yn 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
thumbnail 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 nH with fluctuation level yn. Each point corresponds to the mean of 20 realisations of the model density field nH, and the error bars represent the standard deviation of the fitted spectral indices βn and fluctuation levels yn.

Open with DEXTER
In the text

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

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

Initial download of the metrics may take a while.