© ESO, 2015

1. Introduction

Planck¹ (Tauber et al. 2010; Planck Collaboration I 2011) is the third generation space-mission aimed at mapping the anisotropies of the cosmic microwave background (CMB). With its unprecedented sensitivity and large spectral coverage (nine channels from 30 GHz to 857 GHz) it has provided exquisite maps of that relic radiation (Planck Collaboration I 2014). With its polarimetric capabilities up to 353 GHz, Planck will also provide clues on the physics of the early Universe, by measuring the CMB polarization. However, dominant foreground emission is also partially polarized, masking the primordial signal. In the range of the High Frequency Instrument (HFI, Lamarre et al. 2010), from 100 GHz to 857 GHz, the main contribution to the observed radiation, besides point sources, is thermal emission from dust grains.

The angular momenta of aspherical and spinning grains tend to align with the local magnetic field, although the details of how this alignment proceeds are still the subject of study: see for instance Andersson (2012) for a review on observational constraints regarding grain alignment with respect to current dust models. Submillimetre thermal dust emission is therefore polarized and represents a powerful tool to study interstellar magnetic fields and dust properties. Ideally, we would like to know where in interstellar clouds, and with what efficiency the dust emission and extinction is polarized. This would allow us to use polarization data to infer the spatial structure of the magnetic field. There is an extensive literature on this topic based on observations of starlight polarization, which have been interpreted from two different viewpoints, i.e., grain alignement and magnetic field structure, without achieving a clear understanding of the respective roles of these processes in accounting for variations of polarization across the sky. A number of papers (e.g., Pereyra & Magalhães 2007; Alves et al. 2008; Marchwinski et al. 2012) use the data to infer the magnetic field strength using the Chandreskar-Fermi method (Chandrasekhar & Fermi 1953). Other papers focus on the observed decrease of polarization fraction p with N_H to interpret the data as a decrease of the dust alignment efficiency in dense clouds (Lazarian et al. 1997; Whittet et al. 2008; Chapman et al. 2011).

Magnetohydrodynamical (MHD) simulations provide a theoretical framework to consider both aspects in the interpretation of polarization datasets. Ostriker et al. (2001) were among the first to present simulated polarization maps from MHD simulations, for comparison with data and to study the field structure beyond the simple Chandrasekhar-Fermi method. Falceta-Gonçalves et al. (2008) used a similar technique to study the effect of the Alfvénic Mach number, while Pelkonen et al. (2009) added to this approach the modelling of the alignment process by radiative torques (Hoang & Lazarian 2008).

Planck has mapped the polarized dust emission with great sensitivity and resolution (Planck Collaboration Int. XIX 2014), allowing us to characterize spatial variations of dust polarization and compare data with MHD simulations with unprecedented statistics. This paper is the second in a series of four dealing with a first presentation of the Planck polarized thermal emission from Galactic dust. The other three are the following: Planck Collaboration Int. XIX (2014) describes the polarized dust emission at 353 GHz as seen by Planck over the whole sky and shows in particular that the maximum polarization fraction p_max at a given total gas column density N_H decreases as N_H increases, and that there is an anti-correlation between polarization fractions p and angle dispersion functions , an effect which has also been seen with starlight polarization data (Hatano et al. 2013). Planck Collaboration Int. XXI (2014) compares polarized thermal emission from dust at 353 GHz to polarization in extinction in the visible towards a sample of stars. Finally, Planck Collaboration Int. XXII (2014) discusses the variation of polarized thermal emission from dust with frequency, from 70 to 353 GHz. Both Planck Collaboration Int. XXI (2014) and Planck Collaboration Int. XXII (2014) aim at providing constraints for models of interstellar dust.

In this paper, we use Planck polarization data at 353 GHz to present statistics of polarization fractions and angles in nearby interstellar clouds seen outside the Galactic plane. We then compare the Planck results with simulated observations of polarized thermal dust emission at 353 GHz built from a three-dimensional MHD simulation of the formation of a molecular cloud within colliding flows (Hennebelle et al. 2008).

In these simulated observations, we work under the assumption that the optical properties and the intrinsic polarization fraction of dust grains are constant. At this stage we do not aim at testing models of grain alignment. In this picture, it is expected that the polarization fraction should be maximal when the magnetic field is in the plane of the sky and should, in this case, yield valuable information on the intrinsic polarization fraction. That is why we first focus on the decrease of the maximum value of p, rather than its mean or median values, with increasing column density. We then consider the correlation between polarization fractions and local measures of the dispersion in polarization angles, as it is expected that larger angular dispersions should lower the observed polarization fraction.

The paper is organized as follows. Section 2 describes the Planck data used and the statistics drawn from them in the selected regions. Section 3 presents simulated polarized emission observations based on an MHD simulation of interstellar turbulence and compares their statistical properties with those found towards similar fields in the Planck data. Conclusions are given in Sect. 4. Appendix A presents supplementary figures, and Appendix B details the derivation of the equations yielding the Stokes parameters for dust emission.

2. Planck observations of polarized dust emission

2.1. Planck all-sky data post-processing

The data processing of Planck HFI is presented in Planck Collaboration VI (2014), Planck Collaboration VII (2014), Planck Collaboration VIII (2014), Planck Collaboration IX (2014), and Planck Collaboration X (2014). The specifics of the data processing in terms of polarization are given in Planck Collaboration Int. XIX (2014). We use the same Planck data set as that presented in Planck Collaboration Int. XIX (2014), i.e., full 5-survey HFI mission data for Stokes I, Q, and U at 353 GHz (which is the Planck channel offering the best signal-to-noise ratio for dust polarization) from the “DR3” internal data release. Bandpass mismatch between individual elements of a pair of polarization sensitive bolometers (PSBs) is corrected using in-flight measurements for the dust emission but not for the negligible CO J = 3→2 emission (Planck Collaboration IX 2014). From the total intensity map we subtract the offset I_offset = 0.0887 MJy sr^-1 to set the Galactic zero level at 353 GHz (Planck Collaboration XI 2014). Note that this value includes the cosmic infrared background (CIB) monopole and is slightly different from the one given in Planck Collaboration XI (2014), as the maps are not the same (full mission vs. nominal mission). We do not correct for zodiacal light emission, nor for the residual dipole identified by Planck Collaboration XI (2014) at 353 GHz. CMB and CIB fluctuations are ignored, since the regions selected in this study are outside the CMB-CIB mask described in Planck Collaboration Int. XIX (2014), so the polarized emission there is dominated by the dust.

The Planck polarization and intensity data that we use in this analysis have been generated in exactly the same manner as the data publicly released in March 2013 and described in Planck Collaboration I (2014) and associated papers. Note, however, that the publicly available data include only temperature maps based on the first two surveys. Planck Collaboration XVI (2014) shows the very good consistency of cosmological models derived from intensity only with polarization data at small scales (high CMB multipoles). However, as detailed in Planck Collaboration VI (2014; see their Fig. 27), the 2013 polarization data are known to be affected by systematic effects at low multipoles which were not yet fully corrected, and thus these data were not used for cosmology². We have been careful to check that the Galactic science results in this paper are robust with respect to these systematics³.

We focus in this paper on the polarization fractions p and the polarization angles ψ derived from the Stokes I, Q, and U maps obtained by Planck at 353 GHz and at an angular resolution of 15′. In the absence of noise, p and ψ are defined by (1)and (2)Note that ψ is here defined in the HEALPix⁴ convention (Górski et al. 2005), which means that angles are counted positively clockwise from the north-south direction. Working in that convention instead of the IAU one, which is anti-clockwise (Planck Collaboration Int. XXI 2014), has however no impact on the results presented here. Additionally, since we work on ratios of Stokes parameters, no colour correction is necessary.

When (possibly correlated) noise affects the Stokes parameters, the polarization fraction computed directly using Eq. (1) is biased. We call this one the “naïve” estimator of p, but various methods have been devised to correct for the bias (Montier et al. 2014a), and their respective efficiencies are compared in Montier et al. (2014b). Among them is the modified asymptotic (MAS) estimator introduced by Plaszczynski et al. (2014), which is computed from the naive estimator and the noise covariance matrix pertaining to Q and U. Another estimator of the polarization fraction and angle is the Bayesian estimator described in Montier et al. (2014a) and Planck Collaboration Int. XIX (2014), which has the advantage of taking into account the full noise covariance matrix in I, Q and U, and also taking into account the uncertainty on the zero-level offset for I. In the rest of this paper, except where noted, the maps of polarization fraction p and polarization angle ψ at 353 GHz refer to these Bayesian estimators. The Bayesian method also provides maps of the polarization fraction and angle uncertainties, σ_p and σ_ψ.

For the total hydrogen column density map N_H, we use a conversion from the optical depth at 353 GHz, τ₃₅₃, derived from Planck Collaboration XI (2014): for N_H ≳ 2 × 10²¹ cm^-2, the dust opacity is approximately constant, with σ₃₅₃ = τ₃₅₃/N_H ≃ 1.2 × 10^-26 cm². We are aware that this conversion is crude, with possible variations in dust opacity of the order of 20% to 25%, but our findings do not critically depend on that calibration.

All of the maps used in this study have a HEALPix resolution N_side = 1024.

2.2. Overview of the statistics of polarized emission in various fields

We have selected ten regions, each 12° × 12° in size, that are highlighted in Fig. 1 and whose locations are given in Table 1. These are the same as some of the individual regions mentioned in Planck Collaboration Int. XIX (2014). All of these fields are outside the Galactic plane and probe nearby interstellar material, but they exhibit very different physical conditions, from the diffuse, turbulent ISM with little to no star-forming activity (Polaris Flare), to self-gravitating, star-forming clouds (Orion). They also differ in terms of polarized emission. Some diffuse regions have high polarization fractions (e.g., Pavo), while some have low polarization fractions (e.g., Polaris Flare). This variety of conditions in terms of polarization fraction and gas content is emphasized in Fig. 2, which shows the distribution of p and N_H in these regions, compared with the large-scale distribution shown in Planck Collaboration Int. XIX (2014). The latter is represented by its upper and lower envelopes, computed from the 0.01% and 99.99% percentiles of the p distribution within each bin in column density. All the envelopes of two-dimensional distribution functions shown in this paper are computed in this fashion. Note that to facilitate the comparison with Planck Collaboration Int. XIX (2014), Fig. 2 uses maps at 1° resolution. In the rest of the paper, as already stated, we use 15′ resolution maps.

It appears that for column densities between a few times 10²⁰ cm^-2 and a few times 10²² cm^-2, the selected fields probe most of the range of polarization fractions observed over the whole sky in this range of column densities. The diffuse Polaris Flare field shows low polarization, while high polarization fractions are reached at similar column densities in the Chamaeleon-Musca complex, which, being closer to the Galactic plane, is threaded by the large-scale Galactic magnetic field. Another notable feature of Fig. 2 is the fact that in regions with the largest column densities (Taurus, Orion, and Ophiuchus) the maximum polarization fraction decreases with increasing N_H, and that the slopes are comparable to the large-scale trend.

In the following, we perform statistical analyses of the polarization data in these nearby fields by simply selecting HEALPix pixels whose centres fall within the region of interest, directly from the large-scale maps. Only pixels for which p/σ_p> 3 are retained. This threshold is a reasonable value above which the polarization signal-to-noise ratio is properly estimated (Montier et al. 2014b). Note that some of the fields in Table 1 are quite diffuse (e.g., Pavo), so that the dynamic range in column densities is too small to exhibit a significant relationship between p_max and N_H. These diffuse fields are therefore discarded in the later analysis.

We also build local maps of polarized emission using gnomic projections of the HEALPix maps. These are shown in the middle row panels of Fig. 3 for the Ophiuchus and Chamaeleon-Musca fields. Similar figures for all other fields are given in Appendix A. On all these maps, which share the same scale, we show the polarization fractions p at 353 GHz (colour scale) overlaid with contours of the total gas column density and bars of constant length giving the orientation of the apparent projection of the magnetic field on the plane of the sky. These are built by rotating the 353 GHz polarization bars by 90° so as to recover the average magnetic field orientation in the plane of the sky. In the rest of the paper, we will refer to the rotated polarization bars as the magnetic orientation bars. Note that although they are plotted once every few pixels only, to improve visibility, each of these bars represents the orientation at the given pixel. In other words, beyond the 15′ smoothing performed on the Stokes maps, no further averaging is done to plot the orientation bars on Fig. 3 and similar plots.

The large-scale structure of the Galactic magnetic field appears clearly (see e.g., the top part of the Chamaeleon-Musca field, Fig. 3). There is also a strong correlation between the coherence of the polarization orientation and the level of polarization fraction, in the sense that more ordered regions have higher polarization fractions. This feature, which is already seen at 1° resolution in Planck Collaboration Int. XIX (2014), is discussed later on in Sect. 2.5.

A final qualitative aspect of these maps is that regions with higher column densities tend to be less polarized than their surroundings. An example of this effect can be seen in the Chamaeleon-Musca field (Fig. 3, center right panel) near (l,b) = (301°, −9°), where p ≃ 10%, while it is surrounded by more diffuse material with p ≃ 15%. A future paper (Planck Collaboration Int. XIII 2014) will discuss in more detail the structure of the polarized thermal emission with respect to the morphology of the clouds themselves.

2.3. Maximum polarization fraction

We give in Table 2 the maximum polarization fractions p_max in all the selected fields. Note that for the most diffuse fields Microscopium, Pisces, Perseus, Ara, and Pavo, the quoted values should be taken with caution, since most pixels in these regions have , which corresponds roughly to , and therefore the effect of the (uncertain) zero-level offset on the polarization fraction p may not be negligible.

In the less diffuse fields, the values of p_max are noticeably larger than those found in the same fields at 1° resolution in Planck Collaboration Int. XIX (2014)⁵, which shows the strong effect of spatial resolution on polarization measurements. The uncertainties σ_pmax on the maximum polarization fractions, listed in Table 2, are derived from the various sources of uncertainty involved.

First, the noise properties on the Stokes parameters I, Q, and U in each pixel are described in the data by the noise covariance matrices, which are input in the Bayesian method of Montier et al. (2014a) and Planck Collaboration Int. XIX (2014), and lead to a map of the uncertainty σ_p on the polarization fraction. This includes the 0.0068 MJy sr^-1 uncertainty on the zero-level offset. We then compute the difference σ_pmax,p between the maximum polarization fractions found in the maps of p − σ_p and p + σ_p.

Second, there is a part of the uncertainty related to the method used to debias the data (Montier et al. 2014a). We have computed the standard deviation σ_pmax,d of the maximum polarization fractions obtained in each field when using the “naïve” , modified asymptotic (MAS, Plaszczynski et al. 2014) and Bayesian estimators of p.

Third, we have computed the standard deviation σ_pmax,s of the maximum polarization fractions obtained in each field when considering subsets of the data, namely half-ring maps (one half of each stable pointing period) and detector set maps (one half of the detectors).

The final uncertainty quoted in Table 2 is then given by the quadratic sum (3)It should be noted that the last contribution is usually the dominant one in the selected fields, and that the uncertainty related to the debiasing method is much smaller than the other two.

2.4. Polarization fraction vs. column density

We show the distributions of p and N_H for the Ophiuchus and Chamaeleon-Musca fields in Fig. 4 and for all other fields in Appendix A. The decrease in maximum polarization fraction p_max at higher column densities is apparent for all fields, above a given threshold in N_H that depends on the field and is of the order of 10²¹ to 3 × 10²¹ cm^-2, corresponding to visual extinctions A_V ≃ 0.6 to 1.7, for the fields that are not too diffuse (Polaris Flare, Taurus, Orion, Chamaeleon-Musca, and Ophiuchus). Below this threshold, the polarization fraction may be related to the background more than to the clouds themselves. To quantify the decrease in maximum polarization fraction p_max with increasing N_H, we consider the upper envelope of the distribution of p and N_H, computed as described in Sect. 2.2, and fit this curve with a function , restricted to a range of column densities that depends on the field considered (see Table 2). Note that we perform this fit for the above five fields only, for which there is a large enough dynamic range in column density.

Results of these fits are shown as solid black lines on each panel of Fig. 4, and values of the slopes m and intercepts c are listed in Table 2. Uncertainties on these parameters are derived in the same way as for the maximum polarization fractions p_max in the previous section. The slopes m range between −0.068 for Orion and −0.140 for Taurus, and regions exhibiting stronger column density peaks (e.g., Orion) tend to have shallower slopes than more diffuse molecular clouds (e.g., Polaris Flare).

As mentioned before, the pixels selected for plotting Fig. 4 and performing the fits are those for which the polarization signal-to-noise ratio is p/σ_p> 3. We have checked that modifying this threshold does not change our results, as can be seen in Fig. 5, which shows the same as the top panel of Fig. 4 but with a signal-to-noise ratio threshold p/σ_p> 10. The effect of that stricter selection is to remove points below the original lower envelope, but leaves the upper envelope unchanged. Consequently, both the absolute maximum polarization fraction p_max and the slope of the decrease of p_max at the high end of column densities are quite robust.

2.5. Polarization angle coherence vs. polarization fraction

We show in the bottom row panels of Fig. 3 the maps of the angle dispersion functions for the Ophiuchus and Chamaeleon-Musca fields. Similar maps for all other fields are shown in Appendix A. We recall that this function, defined in Planck Collaboration Int. XIX (2014), is (4)where the sum extends over pixels whose distances from the central pixel r are between δ/ 2 and 3δ/ 2. Here they are computed at a lag δ = 16′, comparable to the size of the beam’s FWHM. One can readily see filamentary structures that correspond to regions where the polarization angle is less ordered or where it changes abruptly. These filaments are already noted at 1° resolution in Planck Collaboration Int. XIX (2014) over several degrees. These regions of large angular dispersions correspond to regions of low polarization fraction, as can be seen for instance by comparing the middle and bottom row panels of Fig. 3.

When increasing the value of the lag δ, we obtain maps of such as that presented in Fig. 6 for the Ophiuchus field at δ = 34′ (approximately twice the FWHM). It appears that the overall value of increases with lag, as already noted in Hildebrand et al. (2009) and Planck Collaboration Int. XIX (2014). However, since has an upper limit of 90°, this means that the anti-correlation with p (see below) will flatten out at large lags. Note however that a completely random sample yields (Planck Collaboration Int. XIX 2014). Values larger than this are few, but they do exist, as can be seen on the maps of in Figs. 3 and 6. They may be linked to sharp boundaries between two well-ordered regions: for instance, the angle dispersion function at the interface between two half-planes with orthogonal magnetic orientations is .

To confirm the visual impression that the spatial coherence of the polarization angle is anti-correlated with the polarization fraction, we show the distribution function of these two quantities for the Ophiuchus and Chamaeleon-Musca fields in Figs. 7 and 8, respectively, and for all other fields in Appendix A.

The large-scale anti-correlation seen in Planck Collaboration Int. XIX (2014) at 1° resolution and δ = 30′ is also present when using a lag close to the beam size. With , we find it to be , where is measured in degrees. Since in this case the ratio δ/ FWHM is the same as for our higher resolution maps (FWHM = 15′ and δ = 16′), we compare the anti-correlations found in the selected fields to this law. Note that the slope −0.75 is similar to the value −0.834 quoted in Planck Collaboration Int. XIX (2014), but the intercept is larger (−0.06 vs. −0.504). This points to a global increase of at larger δ/ FWHM values, which we interpret as a decorrelation of polarization angles at larger lags.

The distributions of p and in the various fields considered show an anti-correlation very similar to the large-scale trend, with slopes and intercepts of the fits through the data points that are very close to the large-scale fit values. When increasing the lag at the same resolution, however, increases and the anti-correlation with p flattens out, as can be seen in Fig. 9. The linear fits for the individual fields are listed in Table 2. The uncertainties on the parameters m′ and c′ are the quadratic sums of uncertainties obtained in three ways: (i) by performing the linear regression using the three estimators of p, i.e., the “naïve”, MAS and Bayesian ones; (ii) by using half-ring maps and detector set maps; (iii) via a Monte-Carlo simulation using the maps of polarization fraction uncertainty σ_p and angle dispersion function uncertainty (Planck Collaboration Int. XIX 2014).

3. Simulations of polarized emission

3.1. Simulations of MHD turbulence

We aim to compare the observed polarization statistics in the selected fields to predictions built on the results of a numerical simulation of MHD turbulence. This simulation is described in detail in Hennebelle et al. (2008)⁶. It follows the formation of clumps of dense and cold gas (cold neutral medium, CNM) out of magnetized warm neutral atomic gas (warm neutral medium, WNM) in an open box of 50 pc on each side, without reaching the stage when cold cores of column density larger than 2 × 10²² cm^-2 form.

The simulation cube initially contains a uniform distribution of WNM with density n_H = 1 cm^-3 and temperature T = 8000 K, and two converging flows of that same gas are injected from opposing faces along the x axis with a velocity ΔV_x ≃ 40 km s^-1 relative to each other. Spatial modulations of the velocity are imposed on the incoming flows, with amplitudes relative to the mean flow of about unity and a periodicity of about 10 pc. Periodic boundary conditions are applied on the remaining four faces. The total mass contained in the cube continuously increases with time. The magnetic field’s initial direction is along that of the incoming flows, and its intensity is about 5 μG, consistent with observational values at these densities (Crutcher et al. 2010). There is therefore a large-scale anisotropic component of the magnetic field throughout the simulation, as well as a turbulent component linked to the velocity perturbations imposed on the converging flows.

These flows collide near the midplane, where the combined effects of cooling and self-gravity eventually lead to the formation of dense (n_H> 100 cm^-3) clumps of cold gas (T of the order of 10–50 K; Hennebelle & Audit 2007). To follow that condensation, the grid is adaptively refined, with an effective (maximum) resolution of 0.05 pc.

In this paper, we select a cubic subset (18 pc × 18 pc × 18 pc) of the density and magnetic field in the simulation snapshot timed at t = 10.9 Myr, which corresponds to an evolved state of the simulation, given the crossing time t_c ≃ 2.4 Myr. The structures present in the simulation are due to the collision of the incoming flows and not to a pure gravitational collapse, since the initial free-fall time is t_ff ≃ 44 Myr. However, some of the densest structures (n_H> 10⁴ cm^-3) may have had time to collapse.

The chosen subset is located near the midplane, so that the influence of boundary conditions is minimal. It contains approximately 3200 M_⊙ of gas; its physical properties are listed in Table 3, and the distribution functions of total gas density n_H and magnetic field components B_x, B_y, B_z are shown in Figs. 10 and 11, respectively. The standard deviations are very similar for all three magnetic field components, but only the x component has a significant mean value, which shows that the mean magnetic field within the cube is approximatively aligned with the x axis, that is with the incoming flows.

We would like to stress here that the MHD simulations we use for comparison with the Planck polarization data do not faithfully reproduce the whole range of densities and column densities spanned by the cloud sample of Table 1, i.e., from diffuse molecular clouds (Polaris Flare) to massive star-forming clouds (Orion). However, as shown in Table 1, only a few percent of the pixels (at most 5.4% in Orion) have column densities larger than 10²² cm^-2 in these fields, the regions of star formation filling only a small fraction of the area in each field. The MHD simulations with their broad range of densities (Fig. 10) and column densities reaching⁷N_H = 1.6 × 10²² cm^-2 are therefore representative of the dynamics of the bulk of the gas. Together with their anisotropy, due to the large-scale magnetic field pervading the cube, these simulations are particularly well suited to analyse the polarization properties of nearby molecular clouds immersed in their low density and large-scale environment.

To compute simulated polarization fractions p, the local gas density n_H and magnetic field components B_x, B_y, B_z are extracted from the simulation and interpolated on a regular grid at the next-to-highest spatial resolution available, so that pixel sizes are approximately Δx = 0.1 pc. These cubes are used in the following section to build simulated polarized emission maps. However, they are first rotated around the y axis, as sketched out in Fig. 13, to explore the full range of possible angles between the mean magnetic field and the line of sight, and therefore to test the effects of the large-scale magnetic field’s anisotropy. The viewing angle α introduced in Fig. 13 is such that the mean magnetic field is approximately in the plane of the sky for α = 0°, and along the line of sight for α = 90°.

3.2. Simulated Planck observations

We build simulated Stokes I, Q, and U maps by integrating along the line of sight (z′ in Fig. 13) through the rotated simulation cube, following the method in Wardle & Königl (1990), Fiege & Pudritz (2000), Pelkonen et al. (2009), and Padovani et al. (2012). Because of a number of inconsistencies in the literature, we give the correct derivation in Appendix B, drawing on the works of Lee & Draine (1985) and Wardle & Königl (1990). This results in: Here p₀ is a polarization fraction parameter related to the intrinsic polarization fraction (see Eq. (8) and Appendix B), γ is the angle that the local magnetic field makes with the plane of the sky, and φ is the local polarization angle in the HEALPix convention. This angle differs by 90° from the angle χ of the plane of the sky projection of the magnetic field, as defined in Fig. 14, and should not be confused with the actual polarization angle ψ. These angles are equal (φ = ψ) only for a uniform magnetic field along the line of sight.

Note that the corrective term in Eq. (5) is incorrectly written in Fiege & Pudritz (2000), Gonçalves et al. (2005), Pelkonen et al. (2009), and Padovani et al. (2012), with p₀/ 2 instead of p₀.

The hypotheses made here, besides the absence of background radiation, are that p₀ = 0.2 is uniform, that the source function S_ν = B_ν(T_d) is that of a blackbody with an assumed uniform dust temperature T_d = 18 K, and that since we are working at 353 GHz the optical depth is simply given by dτ_ν = σ₃₅₃n_H dz′. We use the value σ₃₅₃ = 1.2 × 10^-26 cm² (see Sect. 2.1), and n_H is the total gas density in the simulation. Given the maximum gas column density in the simulation subset computed over all possible viewing angles α, N_H,max = 1.6 × 10²² cm^-2, the maximum optical depth at 353 GHz using this conversion factor is τ_max = 1.9 × 10^-4, so we may safely neglect optical depth effects and take e^−τ_ν = 1 in the I, Q, and U integrals. We are aware (Planck Collaboration XI 2014) that the opacity actually varies with N_H, but the variation is at most a factor of 3 from the value assumed here, so the optical depth is in any case much lower than unity. Moreover, the choice of the conversion factor has no impact on the simulated maps of polarization fractions and angles, provided that a constant value is assumed along each line of sight.

We note that the dense cores that exist in our simulated cube are only weakly shielded from the ambient UV radiation field. Indeed, the mean column density through the cube is about 10²¹ cm^-2 (corresponding to A_V ≃ 0.6), which is comparable to the values in the simulation of Pelkonen et al. (2009), but over a much larger volume (18 pc box compared to less than 1 pc); the bulk of the gas is therefore more fragmented and radiation penetrates more easily (Levrier et al. 2012). That is why we take a uniform parameter p₀.

The maps of Stokes parameters are placed at a distance of D = 100 pc and convolved with a circular 15′ FWHM Gaussian beam (corresponding to a physical size 0.44 pc). The resulting field of view is a little less than 10° across, which is comparable to the selected Planck fields, and small enough that separate smoothing of Stokes I, Q, and U is not an issue (see Appendix A of Planck Collaboration Int. XIX 2014). Maps of polarization fractions and angles are then built from these convolved Stokes parameter maps using Eqs. (1), (2) for consistency with the Planck data. Let us stress that ψ is defined in the HEALPix convention, which means that it is counted positively clockwise from the north-south direction, and not in the IAU convention (anti-clockwise).

Figure 15 (middle row) shows the maps of polarization fraction p and magnetic orientation in these simulated observations, when integrating along the mean magnetic field (α = 90°), and perpendicular to it (α = 0°). The large-scale component of the magnetic field is clearly visible in several regions, for instance in the lower right corner of the α = 0° case: it leads to long-range coherence in the polarization angle, which correlates with the highest polarization fractions and lowest column densities. Conversely, when integrating along the direction of the large-scale field (α = 90°, right column), p is on average much lower, and no such long-range ordering of χ is visible, although some local correlations are present. These effects are expected from the vectorial nature of the polarization: with the magnetic field more or less aligned with the line of sight, only its transverse fluctuations lead to a signal in polarization, and these fluctuations are isotropic in the plane of the sky, so they cancel out in the integration (along the line of sight and also through beam dilution). This correlation between p and spatial coherence of the polarization angle is discussed later on (Sect. 3.4).

Statistics of simulated maps of the polarization fraction (maximum, mean and standard deviation) are shown as a function of the viewing angle α in Fig. 16. We find the maximum polarization fraction to be p_max ≃ 0.14–0.21 (depending on the viewing angle α). On some lines of sight, in the most tenuous parts of the map integrated perpendicularly to the large-scale B (e.g., in the lower right corner of the map in the α = 0° case), p_max almost reaches the theoretical maximum value possible, which is the intrinsic polarization fraction, (8)obtained when the medium is homogeneous and the magnetic field is uniform and parallel to the plane of the sky (γ = 0°). Figure 16 emphasizes the importance of the magnetic field geometry on the measured p_max, as that value varies by about 40% over the range of viewing angles.

3.3. Polarization fraction vs. column density

We show in Fig. 17 the joint distribution function of polarization fractions p and total gas column densities N_H in the simulated observations when integrating along both directions used in Fig. 15, and in the intermediate case α = 45°. The most striking feature of the plots in Fig. 17 is the different behaviour at low column densities N_H < 10²⁰ cm^-2. Along these lines of sight, the density is essentially uniform, with n_H of about 2 cm^-3, so the computed polarization is entirely due to magnetic field geometry; when we integrate with α = 0° the mean magnetic field is almost in the plane of the sky, γ ≃ 0°, and polarized emission is at its highest, while when we integrate with α = 90°, then the ordered field is almost along the line of sight, so γ ≃ 90° and no polarized emission appears. In fact, for each value of α, polarization fractions observed towards the most diffuse lines of sight are well reproduced by the formula for a homogeneous medium, easily derived from Eqs. (5)–(7), (9)as can be seen in Fig. 16. We may therefore only derive the polarization fraction parameter p₀ from the maximum observed value p_max if the angle between the magnetic field and the plane of the sky is known, which is a strong assumption.

The second striking feature of Fig. 17 is the decrease of the maximum polarization fraction with increasing column density, as observed in the data. The same linear fit yields slopes that span values from −0.025 (for α = 80°) to −0.15 (for α = −15°), the latter being comparable to those found in the data for the selected fields.

For a global comparison between simulations and observations, we show in Fig. 18 the distribution of p and N_H for all the simulated fields, with their upper and lower envelopes, together with the envelope for the selected sky fields. Linear fits to the distributions’ upper envelopes are performed, restricted to a common range of column densities 2 × 10²¹ cm^-2<N_H < 2 × 10²² cm^-2. They yield similar values in terms of both slopes (m = −0.109 for simulations, compared to m = −0.113 for the selected fields) and intercepts (c = 2.52 for simulations, compared to c = 2.59 for the selected fields). Note that the “ripple” pattern in the density plot at low N_H is due to the sampling in viewing angles α, and is a signature of the decrease of p with viewing angle for the most diffuse lines of sight, as already noted in Fig. 16.

3.4. Polarization angle coherence vs. polarization fraction

The angle dispersion function is computed from the simulated ψ maps, using a lag δ = 16′, as we did for the data. We first note that the mean angle dispersion function is larger when the large-scale magnetic field is oriented along the line of sight, with for α = 0° and for α = 90°, a result that is consistent with the findings of Falceta-Gonçalves et al. (2008). Maps of (for the α = 0° and α = 90° cases) can be seen in the lower row panels of Fig. 15, exhibiting filamentary patterns similar to those found in observations. These filaments of high also correspond to regions where the polarization angle rotates on small scales, and are correlated with regions of low polarization fraction p (compare with the middle row panels of Fig. 15). This anti-correlation is clearly seen in distribution functions of log (p) and , as shown in Fig. 19 for the α = 0° case. A linear fit to the mean per bin of log (p) is performed, restricted to bins which contain at least 1% of the total number of points and limited to p<p₀ to avoid the most diffuse lines of sight. The slope and intercept of the anti-correlation observed in the data are fairly well reproduced (m′ = −1.0 ± 0.3 and c′ = 0.02 ± 0.34 over the range of α, compared to m′ = −0.75 and c′ = −0.06 in observations) with steeper slopes for viewing angles α ≃ 0° and shallower slopes for viewing angles α ≃ 90° (see Fig. 20). However, since the slopes in simulations are generally steeper than what is observed, but with very similar intercepts at p = 1, the angle dispersion function in simulations is globally higher than in observations for a given polarization fraction.

This result suggests that, in the simulations, the angle dispersion function is too large for a given polarization fraction, i.e., that the magnetic field is too tangled. Since the physical processes one can think of to reduce the field’s tangling (e.g., larger field intensity with respect to turbulence or partial ion-neutral decoupling) would also affect p, we propose that this difference comes from the lack of power in the low frequency modes of the simulated turbulence, as illustrated by the fact that the power spectra of the velocity and magnetic field components flatten out at small wavenumber k. In reality, molecular clouds are organized in a self-similar structure over a broad range of scales and that is therefore not properly reproduced in the simulations we used. In short, the large-scale fluctuations of the magnetic field are closer to random in simulations than in reality.

3.5. Statistics on the magnetic field fluctuations in the simulations

We investigate here the possible causes of the variations in the polarization fraction p and the dispersion of the polarization angle in the simulations, i.e., what are the respective roles of the field tangling and the orientation of the large-scale field in the variations of p and . To quantify these roles, we compute the average and dispersion along the line of sight of both cos²γ and sinχ (see Fig. 14 for the definition of angles). These quantities are computed for different viewing angles. In the following, we write the magnetic field as B = B₀ + ΔB, where B₀ is the large-scale ordered field and ΔB is the fluctuating part of B.

The role of the average values of the angles γ and χ along the line of sight is illustrated in Figs. 21 and 22. First, the role of the large-scale field B₀ is clear: the largest values of p are obtained when B₀ is viewed in the plane of the sky (α ≃ 0°). The largest p values are obtained when the average of cos²γ along the line of sight stays close to unity. In that case, the field perturbations are such that they keep the field close to the plane of the sky, on average, hence the large p. The same effect is visible in the top panel of Fig. 22 where the largest polarization fractions are obtained for average values of χ close to 90°.

However, even in this configuration (α = 0°), small values of p are obtained. The fraction of low p values is clearly larger when the large-scale field is viewed along the line of sight (α = 90°). The remarkable feature visible in Fig. 21 (bottom panel) is the proportionality of p_max with the average of cos²γ: the smaller this average, the closer γ is to 90°, therefore the closer the field is aligned with the line of sight, and the smaller the resulting value of p_max. One also sees in Fig. 21 that reaches much smaller values when B₀ is along the line of sight (bottom panel), producing lower values of p than in the case where B₀ is in the plane of the sky (top panel).

We note, interestingly, that the same effect is not visible in Fig. 22, which displays the line of sight average of sinχ versus p: there is no such upper value of p that would scale with the average of sinχ because this fluctuation of the field direction is measured in the plane of the sky and does not affect the maximal polarization fraction that can be obtained. Instead, when B₀ is along the line of sight for instance, the scatter of ⟨ sinχ ⟩ along the line of sight is the largest and the resulting values of p are low.

Figure 23 also illustrates the effect of the field tangling: the larger the dispersion of cos²γ along the line of sight (and the larger the scatter of this dispersion), the smaller p is. Obviously, when the line of sight is dominated by the large-scale field, the scatter is the lowest.

Figure 24 shows the joint distribution of the average of cos²γ and , where one recognizes the role of the large-scale field when α = 0°: the lowest values of are obtained when γ stays close to 0°, meaning that the field is more or less in the plane of the sky. Clearly, the largest values of are obtained when the influence of the large-scale field is minimized (α = 90°, bottom panel).

4. Conclusions

To summarize, the maximum polarization fraction p_max observed towards the sample of nearby fields selected in this study is reached in the most diffuse fields. The large-scale decrease of p_max with increasing N_H is seen in the individual fields considered here, as soon as N_H> 10²¹ cm^-2. This trend is fairly well reproduced by numerical simulations of anisotropic MHD turbulence, even assuming uniform dust temperatures and grain alignment efficiencies in the gas weakly shielded from the UV radiation. The polarization of thermal dust emission observed by Planck towards these regions is essentially related to the geometry of the magnetic field and in particular to its orientation at large scales with respect to the line of sight. We do not discuss the evolution of polarization fractions at large column densities N_H> 3 × 10²² cm^-2, for which the MHD simulation considered is not suitable. It is clear, however, that additional processes must be at work to achieve the change of slope in the p_max vs. relation observed towards the most opaque lines of sight. This change is probably related to variations in the properties of dust alignment, as pointed out by Soler et al. (2013). We also find that polarization fractions observed by Planck towards these nearby regions correlate well with the local coherence of the polarization angle, which is measured using the angle dispersion function . This correlation is also found in simulations, with slopes that are very close to observational values. In simulations, however, values of for a given polarization fraction are globally too high compared to observations, which points to a possible limitation of the specific MHD simulation used.

Online material

Appendix A: Additional figures

In the main body of the paper, we showed maps and plots for the Chamaeleon-Musca and Ophiuchus fields. In this appendix we show similar figures for the remaining eight fields, in the same order as in Tables 1 and 2. We first show maps similar to Fig. 3 (Figs. A.1 to A.8), then distribution functions of p and N_H similar to Fig. 4 (Figs. A.9 to A.16), and finally distribution functions of and p similar to Fig. 7 (Figs. A.17 to A.24).

Appendix B: Derivation of the Stokes parameters for emission

The derivation of the Stokes equations Eqs. (5)–(7), as presented by Wardle & Königl (1990) based upon Lee & Draine (1985), considers the extinction cross sections C_∥ and C_⊥ for light that is polarized parallel or perpendicular to the grain symmetry axis, and distinguishes oblate and prolate grains. Say that at each point M on the line of sight we define a reference frame (Mx₀y₀z₀) such that z₀ points to the observer, and the local magnetic field B is in the (My₀z₀) plane. With β the angle between B and the angular momentum J of a rotating grain at M, and γ the angle between B and the plane of the sky, as defined in Fig. 14, Lee & Draine (1985) give, for oblate grains and for prolate grains For spherical grains, all these cross-sections are of course equal, C_x0 = C_y0 = C_⊥ = C_∥. The expressions for the Stokes parameters in terms of the cross-sections are where the average ⟨ ... ⟩ is performed on the possible angles β. The equivalent expressions given by Wardle & Königl (1990) are incorrect in omitting the factor 1 / 2 (it is easily checked that our expressions match the expected form of I in the case of spherical grains, and of P/I in the case of fully polarizing grains: 100% polarization when C_y0 = 0.

Computation of the sums and differences of C_x0 and C_y0 for both grain geometries leads to the same expressions for the Stokes parameters where we have introduced the average cross-section (B.11)and the polarization cross section These expressions match those in Martin (1972), Martin (1974), Martin (1975), and Draine & Fraisse (2009); those adopted by Lee & Draine (1985) are a factor 2 larger. The parameter p₀ is then given by (B.14)with R a Rayleigh reduction factor accounting for the chosen form of imperfect alignment (Lee & Draine 1985). Writing the equations for I, Q and U using the optical depth τ_ν (which is small in the submillimetre) in place of the physical position s on the line of sight, one is led to Eqs. (5)–(7).

The intrinsic polarization fraction is easily computed for both grain geometries:

---

Acknowledgments

The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN, JA and RES (Spain); Tekes, AoF and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and PRACE (EU). A description of the Planck Collaboration and a list of its members, including the technical or scientific activities in which they have been involved, can be found at http://www.sciops.esa.int/index.php?project=planck&page=Planck_Collaboration. Some of the results in this paper have been derived using the HEALPix package. The authors would like to thank Charles Beichman for his careful reading of the manuscript and useful comments. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement No. 267934.

References

All Tables

All Figures